In March 2023, the serene forests of Goa experienced forest fires at an alarming scale. While there are many theories regarding the cause of the fires, mapping is helpful to know the scale of the fires and understand the extent of the damage. The sites identified at this stage are also helpful to monitor recovery and possible spread of invasive species in the aftermath of fires. This can aid in efforts to restore and maintain forest health.
In this blogpost, we look at the use of Google Earth Engine (GEE) to locate burn scars and hence identify areas affected by fires. I’ll be outlining the steps taken and sharing the relevant code for each stage of the analysis.
Importing Satellite Imagery
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 image.updateMask(mask).divide(10000);
}
var dataset = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterDate('2023-02-01', '2023-03-31')
// Pre-filter to get less cloudy granules.
.filterBounds(goa)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',1))
.map(maskS2clouds);
print(dataset)
To begin, I implemented the ‘maskS2clouds’ function to filter out cirrus and cloudy pixels. This function uses the QA60 bitmask band of Sentinel imagery to create a mask that identifies cloudless pixels suitable for further processing. This process is applied across the entire image collection dataset.
Creating a burn scar map, necessitates pre-fire and post-fire images, for an estimation of the scarred area. For this purpose, I selected appropriate images of post and pre-fire scenarios from the image collection. I wrote the below code to be able to display all the images date wise in the image collection to the map view.
Display All Images
//Add date attribute to all images to display the images and select the required //
dataset = dataset.map(function(image){
var str = ee.String(image.get('system:index'))
var sdate = str.slice(0,8)
//var sdate = ee.Date(image.get('system:time_start'));
image = ee.Image(image).set('date',sdate)
return image;
})
print(dataset, 'with date')
var listOfImages = dataset.toList(dataset.size())
print(listOfImages)
//Add all the images in the dataset to the Map view //////
for(var i = 0; i < 28 ; i++){
var image = ee.Image(listOfImages.get(i));
var date = image.get('date')
Map.addLayer(image, band_viz, ''+date.getInfo())
}
By attributing a date to each image and displaying them by date, I could visually assess the quality of the images in the collection. The first part of the code adds date as an attribute to all the images in the image collections. The second part displays the images using their date attribute. The output of the code is displayed below. This format of display helped me to gauge the quality of all the images in the collection. This further guided my selection of the suitable images for the pre and post-fires scenarios.
Example of all images displayed in the map view.
With all the images on display, I selected the 13th of February 2023 as the pre-fire image and the 10th of March 2023 as the post-fire image. This approach as in the above display also highlighted if the selected images need to be mosaicked. Mosaicking is joining two or more contiguous satellite images for analysis in order to cover the entire area of interest. Accordingly, I proceeded to convert the images to a list, then selected the required images and mosaicked them.
Selecting Pre-fire image
dataset = dataset.toList(dataset.size()) print(dataset, 'imageCollection') //mosaic of 13/02/2023 - pre fire image. var preFire = ee.ImageCollection([ee.Image(dataset.get(3)),ee.Image(dataset.get(6)),ee.Image(dataset.get(8)),ee.Image(dataset.get(11))]) preFire = (preFire.mosaic()).clip(goa) print(preFire,'preFire') Map.addLayer(preFire,band_viz,'preFire')
Selecting Post-fire image
//mosaic of 10/03/2023 - post fire image. var postFire = ee.ImageCollection([ee.Image(dataset.get(26)),ee.Image(dataset.get(25)),ee.Image(dataset.get(24)),ee.Image(dataset.get(23))]) postFire = (postFire.mosaic()).clip(goa) print(postFire,'postFire') Map.addLayer(postFire,band_viz,'postFire')
The next step involved calculating Normalized Burn Ratio (NBR) for both pre and post-fire images. NBR relies on Near Infrared (NIR) and Short Wave Infrared (SWIR2) reflectance data. Burnt areas exhibit high reflectance in the SWIR2 region and very low reflectance in the NIR region of the electromagnetic spectrum. The NBR is designed to exploit this difference in reflectance along the spectrum. NBR is formulated as follows:
NBR = (NIR - SWIR2) / (NIR + SWIR2)
The NBR was calculated for both pre and post-fire images, using the below code.
Normalized Burn Ratio (NBR) Calculation
//Finding Pre and post fire NBR = normalized burn index ratio = NIR-SWR2/NIR+SWIR2
var PreNumer = preFire.select('B8').subtract(preFire.select('B12'))
var PreDenomin = preFire.select('B8').add(preFire.select('B12'))
var PreNBR = PreNumer.divide(PreDenomin)
//Map.addLayer(PreNBR,{},'Pre-NBR')
var PostNumer = postFire.select('B8').subtract(postFire.select('B12'))
var PostDenomin = postFire.select('B8').add(postFire.select('B12'))
var PostNBR = PostNumer.divide(PostDenomin)
The ensuing step involved determining the difference between pre and post-fire NBR images. The resultant output shows us the extent of burnt area in the region.
Difference between Pre and Post Fire NBR
//Subtracting postNBR fron PreNBR for scar map var scar =PreNBR.subtract(PostNBR) var palette = ['1a9850','91cf60','d9ef8b','ffffbf','fee08b','fc8d59', 'd73027'] Map.addLayer(scar, ,'ScarMap')
In subsequent stages, I incorporated the application of a threshold to isolate the scarred regions, which were then vectorized using a mask.
Define thresholds
// Define thresholds var zones = scar.lt(0.09).add(scar.gt(0.09)).add(scar.gte(0.2)) zones = zones.updateMask(zones.neq(0));
I chose the thresholds by closely examining pixel values in the scar map using the Inspector Tool in GEE Code Editor interface. The ‘Inspector’ displays the pixel value at the point chosen on the map, for every displayed layer.
Inspector tool highlighted in the screengrab.
The output of the ‘define threshold’ step here was a raster with two values, 1 and 2. By this means, all values less than 0.2 were classified as 1 and while those greater than 0.2 were assigned the value 2. Using the code below, I converted the classified raster to a vector, yielding a downloadable KML file to create a comprehensive map of the affected area.
Reduce to Vector
// Convert the zones of the thresholded scars to vectors.
var vectors = zones.addBands(scar).reduceToVectors({
geometry: goa,
//crs: nl2012.projection(),
scale: 10,
geometryType: 'polygon',
eightConnected: false,
labelProperty: 'zone',
reducer: ee.Reducer.mean(),
bestEffort: true,
maxPixels: 3784216672400
});
// Display the thresholds.
Map.setCenter(74.13215, 15.46816, 10);
Map.addLayer(zones, {min: 1, max: 3, palette: ['0000FF', '00FF00', 'FF0000']}, 'raster');
// Make a display image for the vectors, add it to the map.
var display = ee.Image(0).updateMask(0).paint(vectors, '000000', 3);
Map.addLayer(display, {palette: '000000'}, 'vectors');
In summary, employing this method was effective to identify forest fire scars in Goa. The inclusion of ground points of scarred areas were valuable for ground truthing. Thresholding plays a very important role in delineating scarred areas, and hence ground truth points are helpful. If you’re involved in fire scar mapping, this demonstration could be a beneficial reference. Kindly get in touch via our contact form, if you have any further questions.
Link to the code: https://code.earthengine.google.com/c1010c05dd2762d18626b76a90112223
(Note: This code does not include the required shapefiles for the code to work, this is for reference only.)