Description code for article “Using the Google Earth Engine (GEE) for Detection of Burned Areas”

I want to continue my article “Using the Google Earth Engine (GEE) for Detection of Burned Areas” (link) and describe in detail script for detection burned areas. I decided to post here code of this script with comments, shchems and illustrations of wotk this scrip

Link to script on Google Earth Engine

This script shows two variants for detection burned area:

1. Calculating spectral index NBR for before and after forest fire images, download on your computer and compare these scenes in software using function change detection.
2. Calculating spectral index NBR for after forest fire images and select burned areas using threshold. My empirical threshold for NBR index is 0.3

/////////////////////      Block setting  parameters of script   /////////////////////     

/* In this part of the script We can choose setting work of script:    

1. Year for the analysis
2. Choose region for analysis (Russian region and forestry of Arkhangelsk region)
3. Period for creation of cloudless composite. Cloudless composite is used for creation of water mask and cloud mask for removing these elements in final product
4. Period for creation composites after and before forest fires
5. Set threshold of NBR value for determination of burned areas*/

Setting main parameters

//Select year for analysis of burned areas
var year = 2011;

//Select varian 1 or 2
//Variant 1 is trigger for selection of Russian regions
//Variant 2 is trigger for selection of forestries of Arkhangelsk region
var variant = 1;

//reg - number of region (1-89). 29 is number of Arkhangelsk region
var reg = 29;

//NAME_LES - forestry name
var NAME_LES = 'Мезенское';

//set threshold for NBR image. If pixel has less value than 0.3 than this pixel defined as burned area
var BurnedValue = 0.3;

The algorithm for choosing the satellite data and data set

//Processing of selected variant
if (variant==1){
//Upload boundaries of Russian regions
var Russia = ee.FeatureCollection('ft:1z7BRawAQRGo2kznrwVBCD74UuPJIICrvZPEbGIfA');
//Using filter for selection of region
var region = Russia.filterMetadata('name', 'equals', reg);
} 

else if(variant==2){
//Upload boundaries of forestries of Arkhangelsk region
var Arh_les = ee.FeatureCollection('ft:1Z55isSh-K22pFEL0i6O0NeyLJnj6B5ckrDQ6CnNi');
//Using filter for selection of forestry by name
var region = Arh_les.filterMetadata('NAME_LES', 'equals', NAME_LES);
//Upload Hot Spots Arkhangelsk region from 2006 to 2015 years (Product of MODIS)
var HotSpotsKML = ee.FeatureCollection('ft:1i1_y8xKqQqNCZj4jXtw9hNNRbaGymmiU2x1NejNw');
//Using parameter ‘year’ for filter
var HotSpots = HotSpotsKML.filterMetadata('YEAR', 'equals', year).filterBounds(region);

//If you select year less than 2012 script uses images of Landsat 5 TM
var B2 = 'B2';
var B3 = 'B3';
var B4 = 'B4';
var B5 = 'B5';
var B7 = 'B7';
//If you select 2012 year script uses images of Landsat 7 ETM+
if(year==2012){
var satellite = 'Landsat 7 ETM+';
var sat = 'LE7_L1T_TOA';

//If you select 2013 year or next years script uses images of Landsat 8 OLI
} else if(year>=2013){
var satellite = 'Landsat 8 OLI';
var sat = 'LC8_L1T_TOA';
//Bands order of Landsat 8 OLI distinguishes from Landsat 5 TM and Landsat 7 ETM +
var B2 = 'B3';
var B3 = 'B4';
var B4 = 'B5';
var B5 = 'B6';
var B7 = 'B7';
//If you select year 2001 or not more than 2011 year script uses images of Landsat 5 TM
} else if (year>=2001){
var satellite = 'Landsat 5 TM';
var sat = 'LT5_L1T_TOA';}

//If you select year less than 2001, because in this case GEE doesn’t have FIRMS burned areas data
else {
var satellite = 'FIRMS does not have data';
var sat = 'LT5_L1T_TOA';}  

To print name of Satellite and determination of period for satellite images

//Show selected satellite name
print(satellite);
//Show selected collection name
print(sat);
//Define date for FIRMS. This is time range for definition of burned areas.
var StarFire = year + '-05-01';
var EndFire = year + '-10-01';
//Period for creation of cloudless composite
var Start = '2000-7-1';
var End = '2009-9-1';
//Period for creation composite before forest fires
var StartUntilFire = (year-1) + '-08-01';
var EndUntilFire = (year-1) + '-10-01';
//Period for creation composite after forest fires
var StartAfterFire = year + '-08-01';
var EndAfterFire = year + '-10-01';

///////////////////////////////          FIRMS          ///////////////////////////////

/*  Select Product T21 from FIRMS dataset and creation mosaic on region

    – T21 is product of MODIS includes burned areas and has resolution of 500 x 500 m */

Creation buffer using FIRMS data set for

//dataset of FIRMS
var FIRMS = ee.ImageCollection('FIRMS')
 .filterBounds(region)
 .filterDate(StarFire, EndFire);
var T21 = FIRMS.mosaic()
 .clip(region)
 .select('T21');

//Creation mask for vectorization
var zones = T21.gt(0);
zones = zones.updateMask(zones.neq(0));

//Creation of vector
var vectors = zones.addBands(T21).reduceToVectors({
 geometry: region,
 crs: T21.projection(),
 scale: 500,
 geometryType: 'polygon',
 eightConnected: false,
 labelProperty: 'zone',
 reducer: ee.Reducer.mean()});

//Function for creating buffer
var buffer = function(feature) {
 return feature.buffer(3000);};
//Creation buffer
var bounds = vectors.map(buffer);

////////////////////////         Processing Landsat  scenes       ////////////////////////

/*It is main block of script, because here landsat scenes is procceded.

1. Creation freeloud composit, before and after forest fire composites
2. Creation cloud and water masks
3. Removing masks from composites
4. Calculating index NBR for before and after forest fires composites

We can define burned areas using two ways:

1. Calculating NBR for before and after forest fires composites, download NBR layers to computer and  process in program using function change detection.
2. Calculating NBR for after forest fires composte, define bured area using thremhold (var BurnedValue). */

Creation cloudless composite

var L5 = ee.ImageCollection('LANDSAT/LT5_L1T');
var free_clouds = ee.Algorithms.Landsat.simpleComposite({
collection: L5.filterDate(Start, End),
asFloat: true});

Creation composition before forest fires

var collection_1 = ee.ImageCollection('LANDSAT/'+sat)
 .filterBounds(region)
 .filterDate(StartUntilFire, EndUntilFire)
 .sort('CLOUD_COVER',false);

Creation composition after forest fires

var collection_2 = ee.ImageCollection('LANDSAT/'+sat)
.filterBounds(region)
.filterDate(StartAfterFire, EndAfterFire)
.sort('CLOUD_COVER',false);

Removing clouds from copmosits

//Function for removing of clouds
var fDeleteClouds = function(image) {
 var cloud =  ee.Algorithms.Landsat.simpleCloudScore(image);
 var score =  cloud.select(['cloud']).lte(20);
 return image.updateMask(score);
};

//Removing clouds in composite

var composite_1 = collection_1.map(fDeleteClouds)
 .mosaic()
 .clip(region);

var composite_2 = collection_2.map(fDeleteClouds)
 .mosaic()
 .clip(region);

Creation cloud mask and water mask

var clouds_1 = composite_1.select(B3).subtract(free_clouds.select(B3));
var cloud_mask_1 = clouds_1.expression('float(b("'+B3+'"))>0.03 ? 1:0');
var clouds_2 = composite_2.select(B3).subtract(free_clouds.select(B3));
var cloud_mask_2 = clouds_2.expression('float(b("'+B3+'"))>0.03 ? 1:0').clip(bounds);
var water_mask = free_clouds.expression('float(b("'+B4+'"))<0.05&&float(b("'+B5+'"))<0.05?1:0');

Calculating NBR layers and removing clouds

//Calculation of index NBR
var NBR_1 = composite_1.expression('float(b("'+B4+'") - b("'+B7+'")) / (b("'+B4+'") + b("'+B7+'"))');
var NBR_2 = composite_2.expression('float(b("'+B4+'") - b("'+B7+'")) / (b("'+B4+'") + b("'+B7+'"))');

//Removig from compostions NBR_1 and NBR_2 clouds
var other_mask_1 = cloud_mask_1.expression('b("constant")>0?0:1').clip(geometry);
var NBR_1 = NBR_1.updateMask(other_mask_1);
var other_mask_2 = cloud_mask_2.expression('b("constant")>0?0:1').clip(geometry);
var NBR_2 = NBR_2.updateMask(other_mask_2);

Creation mask of burned areas

var burned_mask = NBR_2.lt(BurnedValue);
burned_mask = burned_mask.updateMask(burned_mask.eq(1));

Creation vector from burned mask

var burned_vectors = burned_mask.addBands(NBR_2).reduceToVectors({
 geometry: geometry,    //geometry-region
 crs: NBR_2.projection(),
 scale: 30,
 geometryType: 'polygon',
 eightConnected: false,
 labelProperty: 'zone',
 reducer: ee.Reducer.mean(),
 maxPixels: 5000000000
});

Show received data

var region = ee.Image(0).updateMask(0).paint(region, '000000',2);
Map.addLayer(region,{palette:'#000000'},'Boundary region',0);
var display = ee.Image(0).updateMask(0).paint(vectors, '000000',2);
Map.addLayer(display,{palette: '#4B0082'},'FIRMS'+(year), 0);
Map.addLayer(HotSpots,{color:'ff0000'},'hot spots '+year,0);
Map.addLayer(composite_1, {bands: [B7, B4, B2],min:0, max: 0.3}, 'collection '+(year-1),0);
Map.addLayer(composite_2, {bands: [B7, B4, B2],min:0, max: 0.3}, 'collection '+year,1);
Map.addLayer(NBR_1,{},'NBR '+(year-1),0);
Map.addLayer(NBR_2,{},'NBR '+(year),0);

Map.addLayer(cloud_mask_2,{},'cloud_mask '+(year),0);
var burned = ee.Image(0).updateMask(0).paint(burned_vectors, '000000',2);
Map.addLayer(burned,{palette: 'ff0000'},'burned '+(year),0);

Export to Google Drive

//Export vector of burned areas to Google Drive
Export.table.toDrive({
 collection: burned_vectors,
 description:'burned',
 fileFormat: 'KML'
});

//Export layer NBR to Google Drive
Export.image.toDrive({
image: NBR_1,
description: 'NBR_'+(year-1),
scale: 30,
region: geometry,
maxPixels: 2697651327
});

Export.image.toDrive({
image: NBR_2,
description: 'NBR_'+year,
scale: 30,
region: geometry,
maxPixels: 2697651327
});

//Export vector of burned areas to Google Drive
Export.table.toDrive({
 collection: bounds,
 description:'bounds '+year,
 fileFormat: 'KML'