Smoothing Land Cover with Morphological Processing in GEE

Apply morphological filtering to MODIS land cover in GEE to remove noise and enhance classification consistency using erosion and dilation.

Land cover data, derived from satellite imagery, often contains noise and small, isolated regions that can hinder analysis. Morphological image processing techniques, such as erosion and dilation, can help to smooth and generalize land cover classifications, removing these artifacts. This article will guide you through a Google Earth Engine (GEE) script that demonstrates how to apply morphological processing to a MODIS land cover image to reduce noise and improve the spatial coherence of the classification.

What This Script Does

This Google Earth Engine script processes a MODIS land cover image to smooth the classification and remove small, isolated regions. It employs a series of spatial filtering operations, including focal mode filtering, erosion, and dilation. The script also demonstrates how to force these operations to be performed at the native resolution of the input data using reprojection. This ensures that the processing is applied consistently and avoids artifacts that can arise from processing at varying pixel sizes.

Smothing Land Cover GEE Code Sample

// Morphological processing of land cover. 
// Force projection of 500 meters/pixel, which is the native MODIS resolution.
var SCALE = 500;
// Load a 2001 MODIS land cover image.
var image1 = ee.Image('MODIS/051/MCD12Q1/2001_01_01');
// Select the classification band of interest.
var image2 = image1.select(['Land_Cover_Type_1']);
// Reproject to WGS84 to force the image to be reprojected on load.
// This is just for display purposes, to visualize the input to
// the following operations.  The next reproject is sufficient
// to force the computation to occur at native scale.
var image3 = image2.reproject('EPSG:4326', null, SCALE);
// Smooth with a mode filter.
var image4 = image3.focal_mode();
// Use erosion and dilation to get rid of small islands.
var image5 = image4.focal_max(3).focal_min(5).focal_max(3);
// Reproject to force the operations to be performed at SCALE.
var image6 = image5.reproject('EPSG:4326', null, SCALE);
// Define display parameters with appropriate colors for the MODIS
// land cover classification image.
var PALETTE = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40', // shrub, grass, savannah
    '111149', // wetlands
    'cdb33b', // croplands
    'cc0013', // urban
    '33280d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
].join(',');
var vis_params = {min: 0, max: 17, palette: PALETTE};
// Display each step of the computation.
Map.setCenter(-113.41842, 40.055489, 6);
Map.addLayer(image2, vis_params, 'IGBP classification');
Map.addLayer(image3, vis_params, 'Reprojected');
Map.addLayer(image4, vis_params, 'Mode');
Map.addLayer(image5, vis_params, 'Smooth');
Map.addLayer(image6, vis_params, 'Smooth');

Step-by-Step Explanation

Here's a breakdown of the code:

1. Define the output scale:

   var SCALE = 500;

   This line defines a constant SCALE with a value of 500 meters. This represents the desired pixel size for the output images, which corresponds to the native resolution of the MODIS data.

2. Load the MODIS land cover image:

   var image1 = ee.Image('MODIS/051/MCD12Q1/2001_01_01');

   This line loads a MODIS land cover image from the MCD12Q1 dataset for the year 2001. The ee.Image object represents the image data.

3. Select the land cover classification band:

   var image2 = image1.select(['Land_Cover_Type_1']);

   The select() method is used to extract the 'Land_Cover_Type_1' band, which contains the primary land cover classification information.

4. Reproject the image (for display):

   var image3 = image2.reproject('EPSG:4326', null, SCALE);

   This line reprojects the image to the WGS84 coordinate system ('EPSG:4326') and sets the pixel scale to 500 meters. This reprojection is primarily for display purposes, allowing you to visualize the original image in a standard geographic coordinate system. The null parameter indicates that no specific resampling method is specified, so the default method will be used.

5. Apply a mode filter:

   var image4 = image3.focal_mode();

   The focal_mode() method applies a mode filter to the image. For each pixel, it determines the most frequent land cover value within a neighborhood (the default is a 3x3 neighborhood). This helps to smooth the image and reduce speckle or isolated pixels.

6. Apply erosion and dilation:

   var image5 = image4.focal_max(3).focal_min(5).focal_max(3);

   This line applies a sequence of morphological operations:

   • focal_max(3): Dilation with a radius of 3 pixels. This operation expands the boundaries of land cover classes.
   • focal_min(5): Erosion with a radius of 5 pixels. This operation shrinks the boundaries of land cover classes.
   • focal_max(3): Dilation with a radius of 3 pixels. Another dilation to restore the size of the features after erosion.

   The erosion and dilation operations together remove small, isolated regions of land cover, further smoothing the classification. The order and radii of these operations are important for achieving the desired smoothing effect.

7. Reproject the smoothed image:

   var image6 = image5.reproject('EPSG:4326', null, SCALE);

   This line reprojects the smoothed image (image5) to the WGS84 coordinate system with a pixel scale of 500 meters. This reprojection is crucial to ensure that all previous spatial operations are performed at the native resolution of the MODIS data. This prevents artifacts.

8. Define visualization parameters:

   This section defines a color palette (PALETTE) that maps land cover class values to specific colors, and visualization parameters (vis_params) that specify the minimum and maximum pixel values and the palette to use when displaying the land cover images.

9. Display the results:

   The Map.setCenter() line sets the center and zoom level of the map.

   The Map.addLayer() lines add the original and processed images to the map.

Datasets Used

MODIS Land Cover Type Yearly Global 500m (MCD12Q1), Collection 051: This dataset provides global land cover classifications derived from MODIS data.

Applications

This script demonstrates a valuable technique for improving the quality of land cover data, which has applications in various fields, including:

• Land Use Change Analysis: Smoothing land cover maps can improve the accuracy of detecting changes over time.
• Ecosystem Modeling: More accurate land cover data can lead to better estimates of ecosystem parameters.
• Agricultural Monitoring: Removing noise from land cover classifications can help in identifying agricultural areas.
• Urban Planning: Smoothing urban land cover classifications can improve the delineation of urban boundaries.
• Climate Change Research: Accurate land cover data is essential for modeling.

Notes

• The order and radii of the focal_max() and focal_min() operations are important.
• The reproject() operation is crucial.
• The PALETTE variable defines a color palette.

Related Articles :  

Calculating the Enhanced Vegetation Index (EVI) with GEE

Visualizing Hillshades in Google Earth Engine

Calculating NDVI Using MODIS in Google Earth Engine

---

Prepared by: Jamal Chaaouan | GEE Academy @ GeoJamal.com