Image visualization and Land Use Land Cover (LULC) classification with GEE

Landsat 8 median composite for Punjab state was developed with GEE in this project, followed by several band combinations (true and false color), NDVI, and NDWI for the clipped image for visualization and land use land cover classification.

  1. Before importing the necessary Landsat tiles, a function named maskL8sr() was constructed to mask the overcast pixel on each tile while simultaneously scaling all the Landsat bands (optical and thermal) to their standard values.

  2. // firstly import the shapefile of your desired region which in our case is Punjab
    // you can upload your own shapefile as an asset in GEE and can be imported from assets taskbar on the left side of GEE
    var pt = ee.Geometry.Point(75.46, 31.7008);
    var maskL8sr = function(image) {
        // Bit 0 - Fill
        // Bit 1 - Dilated Cloud
        // Bit 2 - Cirrus
        // Bit 3 - Cloud
        // Bit 4 - Cloud Shadow
        var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
        var saturationMask = image.select('QA_RADSAT').eq(0);
    
        // Apply the scaling factors to the appropriate bands.
        var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
        var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
    
        // Replace the original bands with the scaled ones and apply the masks.
        return image.addBands(opticalBands, null, true)
        .addBands(thermalBands, null, true)
        .updateMask(qaMask)
        .updateMask(saturationMask);
    };
  3. To collect required Landsat tiles, the GEE built-in ImageCollection() function was called, in which filterBounds() selects the tiles related to the input geometry, which in our case was a Punjab shapefile, and filter to the specific date was provided by means of the filterDate() function, and finally, the maskL8sr() function created in step 1 is applied to image collection. Following that, the geometry was reduced to median values, and the composite Landsat visual for the Punjab shapefile is clipped.

  4. var landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                            .filterBounds(pb.geometry()) // here pb is the shapefile needed to import from assets menu bar on the left side of GEE
                            .filterDate('2022-01-01', '2022-12-31')
                            .map(maskL8sr);
    landsat = landsat.median().clip(pb.geometry());
  5. Step 3 estimates a normalized difference for vegetation and water by utilizing the inbuilt normalized difference function. For NDVI & NDWI, bands 5 & 4 are 5 & 6 were used respectively.

  6. var ndvi = landsat.normalizedDifference(['SR_B5','SR_B4']).rename('NDVI');
    var ndwi = landsat.normalizedDifference(['SR_B5','SR_B6']).rename('NDWI');
    landsat = landsat.addBands(ndvi)
                     .addBands(ndwi);
  7. All of these layers, including Landsat true color and false color composites, NDVI, and NDWI, have been put on the map for visualization.

  8. Map.addLayer(landsat.select('NDVI'), {min:-1,max:1,palette:['red', 'white', 'green']}, 'NDVI');
    Map.addLayer(landsat.select('NDWI'), {min:-0.5,max:1,palette:['white', 'blue']}, 'NDWI');
    Map.centerObject(pt, 7);
    Map.addLayer(landsat, 
        {bands: ['SR_B4','SR_B3','SR_B2'],min: 0,max:0.3},
        'True color');
    Map.addLayer(landsat, 
        {bands: ['SR_B5','SR_B4','SR_B3'],min: 0,max:0.3},
        'False color');
    Landsat Image Color Combinations
    Normalized Difference Vegetation and Water Index
  9. MODIS land cover vector data was added to the map to make it easier to identify and confirm sites for various classes of land cover classification. (If you are confident in your ability to select multiple classes appropriately on your own, you can skip this stage.)

  10. var LC = (ee.ImageCollection("MODIS/061/MCD12Q1")
                            .filterDate('2021-01-01','2022-12-31')
                            .sort('system:index',false)
                            .select("LC_Type1")
                            .first()).clip(pb.geometry());
                            
    var LCVis = {
    min: 1.0,
    max: 17.0,
    palette: [
        '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
        'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
        '69fff8', 'f9ffa4', '1c0dff'
    ]};
    //Map.addLayer(LC,LCVis,'Landcover');
  11. Following the selection of 150 points for each class (forest, developed, water, and cropland), all of these points were compiled into a single feature collection for training purposes. Classifications were based on prediction bands 1–7, 10, NDVI, and NDWI. Then variable classifierTraining and param were defined which were used as input for various supervised learning method.

  12. var trainingFeatures = ee.FeatureCollection([Forest, Developed, Water, Cropland]).flatten();
    
    var predictionBands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7',
        'ST_B10', 'NDVI', 'NDWI'];
    
    var classifierTraining = landsat.select(predictionBands)
        .sampleRegions({collection: trainingFeatures,
                    properties: ['class'],
                    scale: 30
        });
    print(classifierTraining);
    
    var param = {features: classifierTraining,
                classProperty: 'class',
                inputProperties: predictionBands
    };
  13. Landsat composite images were classified using Classification and Regression Tress (SmileCART) and Random Forest (RF) algorithms based on defined training points and prediction bands. The generated LULC maps were represented on the map.

  14. var classifierCART = ee.Classifier.smileCart().train(param);
    
    var classified = landsat.select(predictionBands).classify(classifierCART);
    
    var classificationVis = {min:0,
                max:4,
                palette: ['green', 'red', 'blue', 'yellow','white']};
    
    Map.addLayer(classified, classificationVis, 'CART classified');
    
    var classifierRF = ee.Classifier.smileRandomForest(100).train(param);
    
    var RFclassified = landsat.select(predictionBands).classify(
        classifierRF);
    
    // Add classified image to the map.
    Map.addLayer(RFclassified, classificationVis, 'RF classified');
    
    
    // Import the reference dataset.
    var data = ee.FeatureCollection(classifierTraining);
    
    // Define the prediction bands.
    var predictionBands = [
        'SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7',
        'ST_B10', 'NDVI', 'NDWI'
    ];
    
    // Split the dataset into training and testing sets.
    var trainingTesting = data.randomColumn();
    var trainingSet = trainingTesting
        .filter(ee.Filter.lessThan('random', 0.8));
    var testingSet = trainingTesting
        .filter(ee.Filter.greaterThanOrEquals('random', 0.8));
        
    // Train the Random Forest Classifier with the trainingSet.
    var RFclassifier = ee.Classifier.smileRandomForest(100).train({
        features: trainingSet,
        classProperty: 'class',
        inputProperties: predictionBands
    });
    
    // Now, to test the classification (verify model's accuracy), 
    // we classify the testingSet and get a confusion matrix.
    var confusionMatrix = testingSet.classify(RFclassifier)
        .errorMatrix({
            actual: 'class',
            predicted: 'classification'
        });
        
    // Print the results.
    print('Confusion matrix:', confusionMatrix);
    print('Overall Accuracy:', confusionMatrix.accuracy());
    print('Producers Accuracy:', confusionMatrix.producersAccuracy());
    print('Consumers Accuracy:', confusionMatrix.consumersAccuracy());
    print('Kappa:', confusionMatrix.kappa());
    
    LULC Classification with CART and RF Algorithm

    Click to view code