flood mapping using time series

using landsat 8

In this exercise we will use landsat 8 to map surface water in Cambodia. We will apply cloud removal and calculate the NDWI.

Step 1: import the administrative boundaries of Cambodia

//Get the study area country boundary
var StudyArea = ee.FeatureCollection("projects/servir-mekong/Cambodia-Dashboard-tool/boundaries/cambodia_country");

//Add Shape file layer to the Map
Map.addLayer(StudyArea,{},"Cambodia");

//Center the Map with the Study Area
Map.centerObject(StudyArea, 8);

Step 2: import landsat 8 and filter for Cambodia

// Get a Landsat-8 TOA collection. 
var L8_collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA');

//Filter L8 image collection using shape file (or) boundary
var L8_StudyArea = L8_collection.filterBounds(StudyArea);

Step 3: we filter for the year 2020

// Define time range
var startyear = 2019;
var endyear = 2020;
 
// Set date in ee date format
var startdate = ee.Date.fromYMD(startyear,1,1);
var enddate = ee.Date.fromYMD(endyear,12,31);

//Filter L8 Image Collection with Date
var L8_SA_Date = L8_StudyArea.filterDate(startdate,enddate);
print("Total Images:",L8_SA_Date.size());

Step 4: we use a cloud function to remove clouds from the image

//Create Cloud Masking Function to remove clouds
var cloudfunction = function(image){
  // set cloud threshold
  var cloud_thresh = 40;
  //use add the cloud likelihood band to the image
  var CloudScore = ee.Algorithms.Landsat.simpleCloudScore(image);
  //isolate the cloud likelihood band
  var quality = CloudScore.select('cloud');
  //get pixels above the threshold
  var cloud01 = quality.gt(cloud_thresh);
  //create a mask from high likelihood pixels
  var cloudmask = image.mask().and(cloud01.not());
  //mask those pixels from the image
  return image.updateMask(cloudmask);
};

// mask all clouds in the image collection
var L8_NoCloud = L8_SA_Date.map(cloudfunction);

Map.addLayer(L8_NoCloud.mosaic().clip(StudyArea), { min:0.05, max: 0.8, bands: ['B6', 'B5', 'B4']},'Landsat 8 in study  region'); 

Step 5: we calculate the NDWI for each image in the collection

// Create function to calculate NDWI from landsat 8
//NDWI is Normalized Difference Water Index of NIR (B5) and Green (B3) bands
function Calculate_ndwi(img) {
  var ndwi = img.normalizedDifference(['B3', 'B5']).rename('NDWI');
  return img.addBands(ndwi);
}

var L8_ndwi = L8_NoCloud.map(Calculate_ndwi);

Step 6: we calculate 2 percentiles and use that information to derive the inundated area

//Set the threshold to difference between Dry and Wet Seasons
var DIFF_THRESHOLD = 0.4;

// select dry and wet conditions 
var dry = L8_ndwi.select("NDWI").reduce(ee.Reducer.percentile([10]));
var wet = L8_ndwi.select("NDWI").reduce(ee.Reducer.percentile([90]));

var diff = wet.subtract(dry);

var indundatedArea = diff.updateMask(diff.gt(DIFF_THRESHOLD)).clip(StudyArea);

// add layers 
Map.addLayer(dry.clip(StudyArea),{min:-0.3, max:0.4, palette:"white,blue,darkblue"},"dry");
Map.addLayer(wet.clip(StudyArea),{min:-0.3, max:0.8, palette:"white,blue,darkblue"}, "wet");
Map.addLayer(indundatedArea,{palette:"purple"}, "Indundated Area");

step 7 : we overlay the inundate area with the road network

//import road feature
var roads = ee.FeatureCollection('projects/servir-mekong/osm/cambodia/gis_osm_roads');

// create lists to categorize the road networks
var primary = ["primary","primary_link"];
var secondary = ["secondary","secondary_link"];
var tertiary = ["tertiary","tertiary_link"];
var other = ee.List(roads.aggregate_histogram("fclass").keys()).removeAll(primary).removeAll(secondary).removeAll(tertiary);
 
// filter for primary, secondary and tertiary roads
var primaryRoads = roads.filter(ee.Filter.inList("fclass",primary));
var secondaryRoads = roads.filter(ee.Filter.inList("fclass",secondary));
var tertiaryRoads = roads.filter(ee.Filter.inList("fclass",tertiary));
var otherRoads = roads.filter(ee.Filter.inList("fclass",other));

Map.addLayer(primaryRoads.draw("red"),{},"primary roads")
Map.addLayer(secondaryRoads.draw("blue"),{},"secondary roads")
Map.addLayer(tertiaryRoads.draw("green"),{},"tertiary roads")
Map.addLayer(otherRoads.draw("gray"),{},"other roads")

Step 8: add the dataset on buildings

var buildings = ee.FeatureCollection("projects/servir-mekong/undp/KHM_buildings/buildingsKhv1")
Map.addLayer(buildings,{}, "buildings"

One comment

Leave a Reply