Google Earth Engine NDVI – Calculate and Extract Layer of NDVI Change

google-earth-enginendviremote sensingvegetation-index

My experience and expertise in Google Earth Engine are super low. My study area is Sundarban, Bangladesh. I have 2 NDVI from 2 years (2016 & 2017). I want to find the vegetation change by subtracting these two layers. But my study area also contains unwanted water bodies. So, I need to calculate vegetational changes by excluding water bodies. After calculating the difference I need a layer of this changed vegetative area.

Any kind soul to help me out solving this problem?

/// Masking clouds
// // Create the cloud masking function.
var maskClouds = function(image){
  var cloudScore = ee.Algorithms.Landsat.simpleCloudScore(image);
  var cloudLikelihood = cloudScore.select('cloud');
  var cloudPixels = cloudLikelihood.lt(cloud_thresh);
   return image.updateMask(cloudPixels);
};
////For 2016 image
var landsat8_collection_2016 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                                  .filterBounds(Sundarban)
                                  .filterDate('2016-01-01', '2016-01-31');
var landsat8_collection_2016_NoClouds = landsat8_collection_2016.map(maskClouds);
// Reduce the collection to the median value per pixel.
var median_L8_2016 = landsat8_collection_2016_NoClouds.median()
                                                      .clip(Sundarban);
// Create an NDVI image using bands the NIR and red bands (5 and4) for 2016.
var NDVI = median_L8_2016.normalizedDifference(['B5', 'B4']);

////For 2016 image
var landsat8_collection_2017 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                                  .filterBounds(Sundarban)
                                  .filterDate('2017-01-01', '2017-01-30');

var landsat8_collection_2017_NoClouds = landsat8_collection_2017.map(maskClouds);
                                  .filterDate('2017-01-01', '2017-01-30');

// // Reduce the collection to the median value per pixel.
var median_L8_2017 = landsat8_collection_2017_NoClouds.median()
                                                      .clip(Sundarban);

  // Create an NDVI image using bands the NIR and red bands (5 and4) for 2017.
var NDVI17 = median_L8_2017.normalizedDifference(['B5', 'B4']);

Best Answer

I used the water masking approach described here, which uses the Hansen forest dataset to mask out water. That is a global dataset so it may be suitable for your study area. Since you didn't provide a study area, I used a county in Maine as a stand-in (and since you didn't provide a cloud threshold, I used 50).

This approach will use a constant mask for water across the years of your study, rather than adjusting the size of the mask with respect to changing lake conditions. After mapping the cloud mask to your ImageCollections and finding the median NDVI, we can then mask the median images using the water mask. After creating the masked Images, use mathematical operations to find the band difference. Set some visualization parameters, make the subtraction, and check your work with Map.addLayer().

var maineCounties = ee.FeatureCollection('TIGER/2016/Counties')
  .filter(ee.Filter.eq('NAME', 'Waldo'));
print(maineCounties);
var Sundarban = maineCounties;

/// Water mask
// Load or import the Hansen et al. forest change dataset.
var hansenImage = ee.Image('UMD/hansen/global_forest_change_2015');
// Select the land/water mask.
var datamask = hansenImage.select('datamask');
// Create a binary mask.
var watermask = datamask.eq(1);
Map.addLayer(watermask, {}, 'water');

/// Masking clouds
// // Create the cloud masking function.
var maskClouds = function(image){
  var cloudScore = ee.Algorithms.Landsat.simpleCloudScore(image);
  var cloudLikelihood = cloudScore.select('cloud');
  var cloudPixels = cloudLikelihood.lt(50);
   return image.updateMask(cloudPixels);
};
////For 2016 image
var landsat8_collection_2016 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                                  .filterBounds(Sundarban)
                                  .filterDate('2016-01-01', '2016-01-31');
var landsat8_collection_2016_NoClouds = landsat8_collection_2016
  .map(maskClouds);
// Reduce the collection to the median value per pixel.
var median_L8_2016 = landsat8_collection_2016_NoClouds
  .median()
  .clip(Sundarban)
  .updateMask(watermask);

// Create an NDVI image using bands the NIR and red bands (5 and4) for 2016.
var NDVI16 = median_L8_2016.normalizedDifference(['B5', 'B4']);

////For 2017 image
var landsat8_collection_2017 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                                  .filterBounds(Sundarban)
                                  .filterDate('2017-01-01', '2017-01-30');

var landsat8_collection_2017_NoClouds = landsat8_collection_2017
  .map(maskClouds)
  .filterDate('2017-01-01', '2017-01-30');

// // Reduce the collection to the median value per pixel.
var median_L8_2017 = landsat8_collection_2017_NoClouds
  .median()
  .clip(Sundarban)
  .updateMask(watermask);
// Create an NDVI image using bands the NIR and red bands (5 and4) for 2017.
var NDVI17 = median_L8_2017.normalizedDifference(['B5', 'B4']);

// Create some visualization parameters
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
var diffParams = {min: -1, max: 1, palette:['red','white','green']};
// Center the map layer
Map.centerObject(Sundarban, 7);
Map.addLayer(Sundarban, {}, "ROI");
Map.addLayer(NDVI16, ndviParams, "2016 NDVI");
Map.addLayer(NDVI17, ndviParams, "2017 NDVI");

var NDVIdiff = NDVI17.subtract(NDVI16);
Map.addLayer(NDVIdiff, diffParams, "2017-2016 change");
// Adding a mask to retain just the negative values
var lowNDVImask = NDVIdiff.lt(0);
var lowNDVI = NDVIdiff.updateMask(lowNDVImask);
Map.addLayer(lowNDVI,diffParams,"low NDVI only");
// Export a geoTIFF to drive
Export.image.toDrive({
  image: lowNDVI,
  description: 'lowNDVI',
  folder: 'stackExchange',
  region: Sundarban,
  fileFormat: 'GeoTIFF'
});

Here is the product of the above, for my example Maine county. Red areas have a lower NDVI in 2017 than 2016, while green areas have a higher NDVI in 2017 than 2016:

Related Solutions

Google Earth Engine – Fixing NDVI Chart Projection Error

The issue is in used region. Turkey is a Feature Collection with a Multi Polygon including 25561 vertices in 4 individual Polygons. You need to find out an approach that simplifies that geometry. I'm going to use an easy way to test that it works by using Turkey bounds instead complete geometry. Code looks as follows:

var countries=ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');
var Turkey = countries.filter(ee.Filter.eq('country_na','Turkey'));
var vegIndices = ee.ImageCollection('MODIS/006/MOD13A1')
                     .filterBounds(Turkey)
                     .filter(ee.Filter.date('2012-01-01', '2022-03-30'))
                     .select(['NDVI']);

Map.centerObject(Turkey);
Map.addLayer(Turkey);
print("Turkey", Turkey);
print("Turkey bounds", Turkey.geometry().bounds());

print('vegIndices',vegIndices);

var ndviScaled = vegIndices.map(function(image) {
  return image.multiply(0.0001)
   .copyProperties(image, ['system:time_start','system:time_end','system:index']);
});

print('scaled_ndvi',ndviScaled);
var chart =
    ui.Chart.image
        .series({
          imageCollection: ndviScaled,
          region: Turkey.geometry().bounds(),
          reducer: ee.Reducer.mean(),
          scale: 500,
          xProperty: 'system:time_start'
});
print(chart);

After running above code in GEE code editor, I got following result with produced chart for Turkey bounds (only 5 vertices):

Google Earth Engine – Calculate Transition of Different NDVI Classes

There are two steps embedded in your question. Step one is to identify all pixels of a given pre-event NDVI class that changed to a different NDVI class after the event. Step two is to calculate the total area of those pixels.

Step one: You can do this by starting with making a binary map of changed pixels (look for all pixels where you saw a change in NDVI), then looking for pixels that have a value of 1 and that also overlap with a given NDVI class.

var changes = ee.Image(0).where(ndvi1.eq(ndvi2).not(),1)
var just1 = ee.Image(0).where(ndvi1.eq(1).and(changes.eq(1)),1)
var just2 = ee.Image(0).where(ndvi1.eq(2).and(changes.eq(1)),1)
var just3 = ee.Image(0).where(ndvi1.eq(3).and(changes.eq(1)),1)
var just4 = ee.Image(0).where(ndvi1.eq(4).and(changes.eq(1)),1)

Step two: You can do this using ee.Image.pixelArea(), which returns pixel areas as an image object. Since the map of changed pixels is binary, you can multiply it by the pixelArea map to get areas of changed pixels, then use a sum Reducer to add up all that area within your ROI:

var areas1 = just1.multiply(ee.Image.pixelArea())
var areas2 = just2.multiply(ee.Image.pixelArea())
var areas3 = just3.multiply(ee.Image.pixelArea())
var areas4 = just4.multiply(ee.Image.pixelArea())



Map.addLayer(areas1, {}, 'areas1')
var changed1 = areas1.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi.geometry(),
  scale: 5,
  maxPixels: 1e9
});
print(changed1)

The output (changed1) will give you the total area (in square meters) of pixels that used to be class 1 NDVI, but are now not class 1 NDVI. You could do that for all your classes. (Or, if you want, you could use a grouped reducer. This link has some good info about doing that.)

(Also, I would not say that your GEE competency is low, for what it's worth. I would not say this is a beginner-level script, and if you wrote it, or even if someone else wrote it and you're able to manipulate it, that's a pretty high level of competency.)

Best Answer

Related Solutions

Google Earth Engine – Fixing NDVI Chart Projection Error

Google Earth Engine – Calculate Transition of Different NDVI Classes

Related Question