web123456

GEE Case—Using the Sentinel-1 Sentinel-1 Satellite Data in GEE to Draw Flood (Body of Water) Flooding Diagram

  • //Collect information from Copernicus and filter it according to different aspects
  • var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
  • .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  • .filter(ee.Filter.eq('instrumentMode', 'IW'))
  • .filter(ee.Filter.or(
  • ee.Filter.eq('orbitProperties_pass', 'DESCENDING'),
  • ee.Filter.eq('orbitProperties_pass', 'ASCENDING')
  • ));
  • //Filter according to different dates before and after the flood period. We used 2021 data
  • var before = collection
  • .filter(ee.Filter.date('2021-03-01', '2021-03-20'))
  • .filterBounds(square);
  • var after = collection
  • .filter(ee.Filter.date('2021-05-01', '2021-05-20'))
  • .filterBounds(square);
  • //print the collection
  • print(before);
  • print(after);
  • //Mosaic the collected content and cut it to the area of ​​interest. In this case, "square" is the area of ​​interest
  • var before_image = before.select('VH').mosaic().clip(square);
  • var after_image = after.select('VH').mosaic().clip(square);
  • //function
  • function toNatural(img){
  • //"""Function to convert from dB"""
  • return ee.Image(10.0).pow(img.select(0).divide(10.0))
  • }
  • function toDB(img){
  • //"""Function to convert to dB"""
  • return ee.Image(img).log10().multiply(10.0)
  • }
  • function RefinedLee(img) {
  • // """The RL speckle filter
  • // img must be in natural units, . not in dB!
  • // Set up 3x3 kernels"""
  • var bandNames = img.bandNames();
  • var img = toNatural(img);
  • var weights3 = ee.List.repeat(ee.List.repeat(1,3),3);
  • var kernel3 = ee.Kernel.fixed(3,3, weights3, 1, 1, false);
  • var mean3 = img.reduceNeighborhood(ee.Reducer.mean(), kernel3);
  • var variance3 = img.reduceNeighborhood(ee.Reducer.variance(), kernel3);
  • //# Use a sample of the 3x3 windows inside a 7x7 windows to determine gradients and directions
  • var sample_weights = ee.List([[0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0], [0,1,0,1,0,1,0], [0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0]]);
  • var sample_kernel = ee.Kernel.fixed(7,7, sample_weights, 3,3, false);
  • //# Calculate mean and variance for the sampled windows and store as 9 bands
  • var sample_mean = mean3.neighborhoodToBands(sample_kernel);
  • var sample_var = variance3.neighborhoodToBands(sample_kernel);
  • //# Determine the 4 gradients for the sampled windows
  • var gradients = sample_mean.select(1).subtract(sample_mean.select(7)).abs();
  • var gradients = gradients.addBands(sample_mean.select(6).subtract(sample_mean.select(2)).abs());
  • var gradients = gradients.addBands(sample_mean.select(3).subtract(sample_mean.select(5)).abs());
  • var gradients = gradients.addBands(sample_mean.select(0).subtract(sample_mean.select(8)).abs());
  • //# And find the maximum gradient amongst gradient bands
  • var max_gradient = gradients.reduce(ee.Reducer.max());
  • //# Create a mask for band pixels that are the maximum gradient
  • var gradmask = gradients.eq(max_gradient);
  • //# duplicate gradmask bands: each gradient represents 2 directions
  • var gradmask = gradmask.addBands(gradmask);
  • //# Determine the 8 directions
  • var directions = sample_mean.select(1).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(7))).multiply(1);
  • directions = directions.addBands(sample_mean.select(6).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(2))).multiply(2));
  • directions = directions.addBands(sample_mean.select(3).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(5))).multiply(3));
  • directions = directions.addBands(sample_mean.select(0).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(8))).multiply(4));
  • // The next 4 are the not() of the previous 4
  • directions = directions.addBands(directions.select(0).not().multiply(5));
  • directions = directions.addBands(directions.select(1).not().multiply(6));
  • directions = directions.addBands(directions.select(2).not().multiply(7));
  • directions = directions.addBands(directions.select(3).not().multiply(8));
  • //# Mask all values that are not 1-8
  • var directions = directions.updateMask(gradmask);
  • //# "collapse" the stack into a singe band image (due to masking, each pixel has just one value (1-8) in it's directional band, and is otherwise masked)
  • var directions = directions.reduce(ee.Reducer.sum());
  • //#pal = ['ffffff','ff0000','ffff00', '00ff00', '00ffff', '0000ff', 'ff00ff', '000000']
  • //#((()), {min:1, max:8, palette: pal}, 'Directions', False)
  • var sample_stats = sample_var.divide(sample_mean.multiply(sample_mean));
  • //# Calculate localNoiseVariance
  • var sigmaV = sample_stats.toArray().arraySort().arraySlice(0,0,5).arrayReduce(ee.Reducer.mean(), [0]);
  • //# Set up the 7*7 kernels for directional statistics
  • var rect_weights = ee.List.repeat(ee.List.repeat(0,7),3).cat(ee.List.repeat(ee.List.repeat(1,7),4));
  • var diag_weights = ee.List([[1,0,0,0,0,0,0], [1,1,0,0,0,0,0], [1,1,1,0,0,0,0],
  • [1,1,1,1,0,0,0], [1,1,1,1,1,0,0], [1,1,1,1,1,1,0], [1,1,1,1,1,1,1]]);
  • var rect_kernel = ee.Kernel.fixed(7,7, rect_weights, 3, 3, false);
  • var diag_kernel = ee.Kernel.fixed(7,7, diag_weights, 3, 3, false);
  • //# Create stacks for mean and variance using the original kernels. Mask with relevant direction.
  • var dir_mean = img.reduceNeighborhood(ee.Reducer.mean(), rect_kernel).updateMask(directions.eq(1));
  • var dir_var = img.reduceNeighborhood(ee.Reducer.variance(), rect_kernel).updateMask(directions.eq(1));
  • var dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)));
  • var dir_var= dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), diag_kernel).updateMask(directions.eq(2)));
  • //# and add the bands for rotated kernels
  • for (var i=1; i<4; i++) {
  • dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)));
  • dir_var = dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)));
  • dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)));
  • dir_var = dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)));
  • }
  • //# "collapse" the stack into a single band image (due to masking, each pixel has just one value in it's directional band, and is otherwise masked)
  • dir_mean = dir_mean.reduce(ee.Reducer.sum());
  • dir_var = dir_var.reduce(ee.Reducer.sum());
  • //# A finally generate the filtered value
  • var varX = dir_var.subtract(dir_mean.multiply(dir_mean).multiply(sigmaV)).divide(sigmaV.add(1.0));
  • var b = varX.divide(dir_var)
  • var result = dir_mean.add(b.multiply(img.subtract(dir_mean))).arrayFlatten([['sum']]).float();
  • return ee.Image(toDB(result)).rename(bandNames).copyProperties(img)
  • }
  • //Convert the new mosaic to natural color, apply fine plum spot filtering, and then convert back to db
  • var before_filtered = ee.Image(toDB(RefinedLee(toNatural(before_image))));
  • var after_filtered = ee.Image(toDB(RefinedLee(toNatural(after_image))));
  • //The threshold value (pixel) depends on the region of interest
  • var flood = before_filtered.gt(-20).and(after_filtered.lt(-20));
  • //Show value equal to 1 as flood
  • var flood_mask = flood.updateMask(flood.eq(1));
  • //This kind of water mask is used to distinguish between water and flooded areas
  • var water = before_filtered.lt(-20).and(after_filtered.lt(-20));
  • var water_mask = water.updateMask(water.eq(1));
  • //Open the center of the map at the area of ​​interest
  • Map.centerObject(square);
  • //Load different layers
  • Map.addLayer(before_filtered, { min: -25, max: 0 }, 'Before_Filt');
  • Map.addLayer(after_filtered, { min: -25, max: 0 }, 'After_Filt');
  • Map.addLayer(flood_mask, { palette: ['Yellow'] }, 'Flood_Inundation');
  • Map.addLayer(water_mask, { palette: ['Blue'] }, 'Water');
  • //Calculate the area and flooded area and print it on the screen
  • var stats = flood_mask.multiply(ee.Image.pixelArea()).reduceRegion({
  • reducer: ee.Reducer.sum(),
  • geometry: square,
  • scale: 10,
  • maxPixels: 1e13,
  • tileScale: 16
  • });
  • print(stats)
  • var flood_area = ee.Number(stats.get('constant')).divide(10000).round();
  • print('Flooded Area (Ha)', flood_area);