GEE案例——利用GEE中的Sentinel-1哨兵 1 号卫星数据绘制洪水（水体）淹没图

//从哥白尼那里收集资料，并根据不同方面进行筛选

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')

));

//根据洪水期前后的不同日期进行筛选。我们使用了 2021 年的数据

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);

//对采集的内容进行马赛克拼接，并裁剪到感兴趣的区域，本例中 “正方形 ”就是感兴趣的区域

var before_image = before.select('VH').mosaic().clip(square);

var after_image = after.select('VH').mosaic().clip(square);

//函数

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)

}

//将新马赛克转换为自然色，应用精细李斑点滤波，然后再转换回 db

var before_filtered = ee.Image(toDB(RefinedLee(toNatural(before_image))));

var after_filtered = ee.Image(toDB(RefinedLee(toNatural(after_image))));

//阈值（像素）值取决于感兴趣的区域

var flood = before_filtered.gt(-20).and(after_filtered.lt(-20));

//将等于 1 的值显示为洪水

var flood_mask = flood.updateMask(flood.eq(1));

//这种水面具用于区分水体和淹没区

var water = before_filtered.lt(-20).and(after_filtered.lt(-20));

var water_mask = water.updateMask(water.eq(1));

//将地图中心对准感兴趣的区域

Map.centerObject(square);

//加载不同图层

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');

//计算区域面积和淹没面积并打印在屏幕上

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);