Springe direkt zu Inhalt

Comparison and Analysis of the Results

To close this chapter of monitoring NDVI and NBR, we will look at a summarized example on how to efficiently monitor, analyze and compare the NDVI and NBR of any research area. This will give us a chance to test and strenghten the skills that we acquired by now.


To do so, we will have a look at wildfires that occured during the period of 09/01/2020 to 10/15/2020 in Syria.
As Lattakia was severely affected by the wildfires, we will have a look at this region to see if we can see the impacts of it in our data.


This time, we will be using Sentinel 2 data, as it has a higher spatial resolution and forces us to not just simply copy the previous code, but adapt our knowledge to another sensors properties.

Monitoring and Analysing the NDVI and NBR - Wildfires - Part 1

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Let us start by collecting all the parameters, functions and variable that we will need often later on.

//Visualization parameters:
var visParams = {'min': 400,'max': [4000,3000,3000],   'bands':'B8,B4,B3'};
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
var nbrParams = {min: -1, max: 1, palette: ['red', 'white', 'green']};
var dndviParams = {min: -1, max: 1, palette: ['DarkGreen', 'green', 'LimeGreen', 'white', 'burlywood', 'brown', 'maroon']};
var dnbrParams = {min: -1, max: 1, palette: ['DarkGreen', 'green', 'LimeGreen', 'white', 'burlywood', 'brown', 'maroon']};

//Function to add NDVI to an ImageCollection
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

//Function to add NBR to an ImageCollection
var addNBR = function(image) {
  var nbr = image.normalizedDifference(['B8', 'B12']).rename('NBR');
  return image.addBands(nbr);
};

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//Access Geometry data of the Region of Lattakia
var extent_aoi = ee.FeatureCollection("FAO/GAUL/2015/level1")
                  .filter(ee.Filter.eq('ADM1_NAME', 'Lattakia')); //filter for entry that equals the UN country name 'Syrian Arab Rebuplic'
print(extent_aoi, 'Extent AOI');
Map.addLayer(extent_aoi, {},'Extent AOI');
Map.centerObject(extent_aoi, 10)


//Step 2: Access Sentinel-2 data from before and after the wildfires. Filter for less than 10% cloud cover
var s2_pre = ee.ImageCollection('COPERNICUS/S2_SR')
        .filterBounds(extent_aoi)
        .filterDate('2020-08-01', '2020-08-25')
        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10));
print(s2_pre, 'Image Collection Before Wildfire 2020 Lattakia <10% clouds'); //use the console to look at the amount of elements accumulated
Map.addLayer(s2_pre, visParams ,'Image Collection Before Wildfire 2020 Lattakia <10% clouds', false);   // Visualize the ImageCollection in the Map panel and assign a name to it; Hide the layer by default, as it takes long to render

var s2_post = ee.ImageCollection('COPERNICUS/S2_SR')
        .filterBounds(extent_aoi)
        .filterDate('2020-10-15', '2020-10-31')
        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10));
print(s2_post, 'Image Collection After Wildfire 2020 Lattakia <10% clouds'); //use the console to look at the amount of elements accumulated
Map.addLayer(s2_post, visParams ,'Image Collection After Wildfire 2020 Lattakia <10% clouds', false);   // Visualize the ImageCollection in the Map panel and assign a name to it; Hide the layer by default, as it takes long to render


//As Wildfires tend to correlate with obstructed view due to smoke and clouds, we will use the Quality Assessment bands of our S2-dataset to mask out clouds.
//To do so, we will need to write a short function which is adapted from one of the example scripts included in the GEE; browse Examples -> Cloud Masking -> Sentinel2
function maskS2clouds(image) {
  var qa = image.select('QA60')

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0))

  // Return the masked and scaled data, without the QA bands.
  return image.updateMask(mask)
      .select("B.*")
      .copyProperties(image, ["system:time_start"])
}

//Apply the cloud mask to the ImageCollection.
var s2_pre_masked = s2_pre.map(maskS2clouds);
Map.addLayer(s2_pre_masked, visParams, 'Cloudmask Image Collection Before Wildfire 2020 Lattakia', false)

var s2_post_masked = s2_post.map(maskS2clouds);
Map.addLayer(s2_post_masked, visParams, 'Cloudmask Image Collection After Wildfire 2020 Lattakia', false)

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


//Calculate NDVI and NBR and add it to the masked ImageCollections
var s2_pre_indices = s2_pre_masked.map(addNDVI).map(addNBR);
print(s2_pre_indices, 'Image Collection Before Wildfire with NDVI and NBR');

var s2_post_indices = s2_post_masked.map(addNDVI).map(addNBR);
print(s2_post_indices, 'Image Collection After Wildfire with NDVI and NBR');

//Mosaic the Imagecollection and clip it to the region of Lattakia
var s2_pre_mosaic = s2_pre_indices.mosaic().clip(extent_aoi);
print(s2_pre_mosaic, 'Mosaic Image Before Wildfire with NDVI and NBR');
Map.addLayer(s2_pre_mosaic, visParams, 'Mosaic Image Before Wildfire with NDVI and NBR', false);

var s2_post_mosaic = s2_post_indices.mosaic().clip(extent_aoi);
print(s2_post_mosaic, 'Mosaic Image After Wildfire with NDVI and NBR');
Map.addLayer(s2_post_mosaic, visParams, 'Mosaic Image After Wildfire with NDVI and NBR', false);

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//Create a Difference Image from both Images and create a separate Variable for NDVI and NBR
var diff_img = s2_pre_mosaic.subtract(s2_post_mosaic);
print(diff_img, 'Difference Image Wildfire 2020 Lattakia');
Map.addLayer(diff_img, visParams, 'Difference Image Wildfire 2020 Lattakia', false);
//Here, we can already see where differences between both Images occur, but we have no basis to interprete the differences

var dndvi = diff_img.select('NDVI');
print(dndvi, 'NDVI Difference Image 2020 Lattakia');
Map.addLayer(dndvi, dndviParams, 'NDVI Difference Image 2020 Lattakia', false);
//NDVI Difference Image: Brown areas are less vegetated after the wildfire, green areas are more vegetated. Watch out, as most of these areas are actually clouds we weren't able to mask!

var dnbr = diff_img.select('NBR');
print(dnbr, 'NBR Difference Image 2020 Lattakia');
Map.addLayer(dnbr, dndviParams, 'NBR Difference Image 2020 Lattakia', false);
//NBR Difference Image: Brown areas are either bare soil or recently burnt, green areas are healthy vegetation. We can see that the patterns are the same, which is good!
//Burnt areas stand out a lot clearer than for the NDVI Difference Image, as the NBR is specifically designed to detect them. The NDVI only detects healthy vegetation, which is logically similar, but not the same.
//The NBR is also less prone to cloud disturbances!
//We can also clearly see that water bodies are marked as green areas in both Difference Images, which is due to water having really low pixel values at both times. This does not necessarily mean that there is a lot of vegetation present!

//Use the Inspector to look at NDVI- and NBR-values for specific pixels of the pre-event, post-event and difference image-NDVI!

Monitoring and Analysing the NDVI and NBR - Wildfires - Part 2

//Charting the NDVI and NBR
//To keep processing times low and reduce the chance of producing errors, we are going to downsample our data to a spatial resolution of 200m
var chart_pre_event = ui.Chart.image.series({
  imageCollection: s2_pre_indices.select('NDVI','NBR'),
  region: extent_aoi,
  reducer: ee.Reducer.mean(),
  scale: 200
})
.setOptions({
          title: 'Mean NDVI and NBR Values by Months Pre Event',
          hAxis: {title: 'Month', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'NDVI',titleTextStyle: {italic: false, bold: true}},
  });
print(chart_pre_event)

var chart_post_event = ui.Chart.image.series({
  imageCollection: s2_post_indices.select('NDVI','NBR'),
  region: extent_aoi,
  reducer: ee.Reducer.mean(),
  scale: 200
})
.setOptions({
          title: 'Mean NDVI and NBR Values by Months Post Event',
          hAxis: {title: 'Month', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'NDVI',titleTextStyle: {italic: false, bold: true}},
  });
print(chart_post_event)


//Add a day-of-year-chart of the mean NDVI values
var doychart_pre_event = ui.Chart.image
  .doySeries({
    imageCollection: s2_pre_indices.select('NDVI','NBR'),
    region: extent_aoi,
    regionReducer: ee.Reducer.mean(),
    scale: 200,
    yearReducer: ee.Reducer.mean(),
    startDay: 1,
    endDay: 365
  })
  .setOptions({
          title: 'Mean NDVI and NBR Values by Date Pre Event',
          hAxis: {title: 'Day of year', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'NDVI',titleTextStyle: {italic: false, bold: true}},
  });

print(doychart_pre_event);


var doychart_post_event = ui.Chart.image
  .doySeries({
    imageCollection: s2_post_indices.select('NDVI','NBR'),
    region: extent_aoi,
    regionReducer: ee.Reducer.mean(),
    scale: 200,
    yearReducer: ee.Reducer.mean(),
    startDay: 1,
    endDay: 365
  })
  .setOptions({
          title: 'Mean NDVI and NBR Values by Date Post Event',
          hAxis: {title: 'Day of year', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'NDVI',titleTextStyle: {italic: false, bold: true}},
  });
print(doychart_post_event);


//Histogram
//Plotting histograms of multiple variables can be buggy in the GEE, so it is advised to plot seperate histograms for NDVI and NBR
//The added option 'explorer: {} allows to pan and zoom the histograms in the Console.

//NDVI Histograms
var histogram_pre_event_ndvi =
    ui.Chart.image
        .histogram({
          image: s2_pre_mosaic.select('NDVI'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of NDVI Before Wildfires',
          hAxis: {title: 'NDVI'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        });
print(histogram_pre_event_ndvi);

var histogram_post_event_ndvi =
    ui.Chart.image
        .histogram({
          image: s2_post_mosaic.select('NDVI'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of NDVI After Wildfires',
          hAxis: {title: 'NDVI'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        });
print(histogram_post_event_ndvi);

var histogram_diff_img_ndvi =
    ui.Chart.image
        .histogram({
          image: diff_img.select('NDVI'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of Wildfires NDVI-Difference Image',
          hAxis: {title: 'NDVI'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        })
;
print(histogram_diff_img_ndvi);


//NBR Histograms
var histogram_pre_event_nbr =
    ui.Chart.image
        .histogram({
          image: s2_pre_mosaic.select('NBR'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of NBR Before Wildfires',
          hAxis: {title: 'NBR'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        });
print(histogram_pre_event_nbr);

var histogram_post_event_nbr =
    ui.Chart.image
        .histogram({
          image: s2_post_mosaic.select('NBR'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of NBR After Wildfires',
          hAxis: {title: 'NBR'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        });
print(histogram_post_event_nbr);

var histogram_diff_img_nbr =
    ui.Chart.image
        .histogram({
          image: diff_img.select('NBR'),
          region: extent_aoi,
          scale: 200,
          minBucketWidth: 0.1,
        })
        .setOptions({
          title: 'Histogram of Wildfires NBR-Difference Image',
          hAxis: {title: 'NBR'},
          vAxis: {title: 'Frequency'},
          explorer: {}
        })
;
print(histogram_diff_img_nbr);

Monitoring and Analysing the NDVI and NBR - Wildfires - Part 3

//Classification
//This time, we will classify the Difference Images for the NDVI and the NBR. This way, we can instantly see which regions changed in which way.

//Remap NDVI-Difference-Image into 7 classes
var thresholds_ndvi = ee.Image([-0.4, -0.2,-0.1, 0.1, 0.3, 0.5, 1, 1.5]);
var classified_ndvi = diff_img.select('NDVI').gt(thresholds_ndvi).reduce('sum').toInt();
print(classified_ndvi, 'Classified NDVI Difference Image Wildfires Lattakia');
var classifiedParams_ndvi = {min: 0, max: 7, palette: ['blue', 'green','white', 'white', 'burlywood', 'Maroon', 'brown']};
Map.addLayer(classified_ndvi.clip(extent_aoi), classifiedParams_ndvi, 'Classified NDVI Difference Image Wildfires Lattakia');

// Now that we have classified our NDVI-Image, let us have a look at area statistics for the single classes.
// The following approach is adapted from the UN Spider Burn Severity Mapping Recommended Practices. Make sure to have a look at it, as this is highly informative and really well written!
// https://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/burn-severity-earth-engine


// First, we want to count the number of pixels in the entire layer for future reference
var allpix_ndvi =  classified_ndvi.updateMask(classified_ndvi);  // mask the entire layer
var pixstats_ndvi = allpix_ndvi.reduceRegion({
  reducer: ee.Reducer.count(),               // count pixels in a single class
  geometry: extent_aoi,
  scale: 10,
  maxPixels: 1e10
  });
var allpixels_ndvi = ee.Number(pixstats_ndvi.get('sum')); // extract pixel count as a number


// Then, we want to create an empty list to store the area values we will calculate in
var arealist_ndvi = [];

// Now, we can create a function to derive the extent of one NDVI class
// The arguments are class number (cnr) and class name (name)
var areacount_ndvi = function(cnr, name) {
 var singleMask_ndvi =  classified_ndvi.updateMask(classified_ndvi.eq(cnr));  // mask a single class
 var stats_ndvi = singleMask_ndvi.reduceRegion({
  reducer: ee.Reducer.count(),               // count pixels in a single class
  geometry: extent_aoi,
  scale: 10,
  maxPixels: 1e10
  });
var pix_ndvi =  ee.Number(stats_ndvi.get('sum'));
var hect_ndvi = pix_ndvi.multiply(100).divide(10000);                // Landsat pixel = 30m x 30m --> 900 sqm
var perc_ndvi = pix_ndvi.divide(allpixels_ndvi).multiply(10000).round().divide(100);   // get area percent by class and round to 2 decimals
arealist_ndvi.push({Class: name, Pixels: pix_ndvi, Hectares: hect_ndvi, Percentage: perc_ndvi});
};

// Create a list that contains the NDVI class names (7 classes, ranging from [-0.4, -0.2,-0.1, 0.1, 0.3, 0.5, 1, 1.5])
var names2_ndvi = ['Water / Clouds', 'Vegetation Gain', 'Very Low Vegetation Gain',
'Very Low Vegetation Loss', 'Moderate Vegetation Loss','Moderate-high Vegetation Loss', 'High Vegetation Loss'];

// execute function for each class
for (var i = 0; i < 7; i++) {
  areacount_ndvi(i, names2_ndvi[i]);
  }

print('Vegetated Area by NDVI Class', arealist_ndvi, '--> click list objects for individual classes');

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//For the NBR, we will work according to the official UN-SPIDER recommended practices.
//This includes scaling the data by the factor 1000.
//Remap NBR-Difference-Image into 8 classes
var thresholds_nbr = ee.Image([-1000, -251, -101, 99, 269, 439, 659, 2000]);
var classified_nbr = diff_img.select('NBR').multiply(1000).gt(thresholds_nbr).reduce('sum').toInt();
print(classified_nbr, 'Classified NBR Difference Image Wildfires Lattakia');
var classifiedParams_nbr = {min: 0, max: 8, palette: ['#ffffff','#7a8737', '#acbe4d','#0ae042', '#fff70b', '#ffaf38', '#ff641b', '#a41fd6']};
Map.addLayer(classified_nbr.clip(extent_aoi), classifiedParams_nbr, 'Classified NBR Difference Image Wildfires Lattakia');

//Now that we have classified our NBR-Image, let us have a look at area statistics for the single classes.
//The following approach is adapted from the UN Spider Burn Severity Mapping Recommended Practices. Make sure to have a look at it, as this is highly informative and really well written!
//https://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/burn-severity-earth-engine


// First, we want to count the number of pixels in the entire layer for future reference
var allpix_nbr =  classified_nbr.updateMask(classified_nbr);  // mask the entire layer
var pixstats_nbr = allpix_nbr.reduceRegion({
  reducer: ee.Reducer.count(),               // count pixels in a single class
  geometry: extent_aoi,
  scale: 20,
  maxPixels: 1e10
  });
var allpixels_nbr = ee.Number(pixstats_nbr.get('sum')); // extract pixel count as a number


// Then, we want to create an empty list to store the area values we will calculate in
var arealist_nbr = [];

// Now, we can create a function to derive the extent of one NBR class
// The arguments are class number (cnr) and class name (name)
var areacount_nbr = function(cnr, name) {
 var singleMask_nbr =  classified_nbr.updateMask(classified_nbr.eq(cnr));  // mask a single class
 var stats_nbr = singleMask_nbr.reduceRegion({
  reducer: ee.Reducer.count(),               // count pixels in a single class
  geometry: extent_aoi,
  scale: 20,
  maxPixels: 1e10
  });
var pix_nbr =  ee.Number(stats_nbr.get('sum'));
var hect_nbr = pix_nbr.multiply(400).divide(10000);                // Landsat pixel = 30m x 30m --> 900 sqm
var perc_nbr = pix_nbr.divide(allpixels_nbr).multiply(10000).round().divide(100);   // get area percent by class and round to 2 decimals
arealist_nbr.push({Class: name, Pixels: pix_nbr, Hectares: hect_nbr, Percentage: perc_nbr});
};

// Create a list that contains the NDVI class names (7 classes, ranging from [-1000, -251, -101, 99, 269, 439, 659, 2000])
var names2_nbr = ['Enhanced Regrowth, High', 'Enhanced Regrowth, Low', 'Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity', 'NA'];

// execute function for each class
for (var i = 0; i < 8; i++) {
  areacount_nbr(i, names2_nbr[i]);
  }

print('Burned Area by NBR Class', arealist_nbr, '--> click list objects for individual classes');

We can now visually and numerically compare the results of our monitoring analysis.

We have calculated the NDVI and NBR as well as their respective Difference Images of our research area, displayed them in maps and charts, classified them and calculated the percentage and area of the single classes.