Skip to content

Commit

Permalink
add supplement scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@spatialthoughts
spatialthoughts committed Feb 17, 2024
1 parent 2d0f302 commit 0df435e
Show file tree
Hide file tree
Showing 6 changed files with 938 additions and 1,023 deletions.
171 changes: 35 additions & 136 deletions code/gee_water_resources_management/Supplement/SAR/SAR_Water_Detection
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,144 +7,43 @@ var geometry = basin.geometry();
Map.addLayer(basin, {color: 'grey'}, 'Arkavathy Sub Basin')
Map.centerObject(geometry, 12);

var s1 = ee.ImageCollection("COPERNICUS/S1_GRD")
var filtered = s1
// Filter to get images with VV and VH dual polarization.
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
// Change the pass to ASCENDING depending on your location
.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
.filterDate('2019-01-01', '2020-01-01')
.filterBounds(geometry)

// Apply Speckle Filter
var speckleFiltered = filtered.map(refinedLee);
// Use Sentinel-1 ARD Pre-Processing
// See https://github.com/adugnag/gee_s1_ard/
var wrapper = require('users/adugnagirma/gee_s1_ard:wrapper');
var helper = require('users/adugnagirma/gee_s1_ard:utilities');


var parameters = {
// Input
START_DATE: '2019-01-01',
STOP_DATE: '2020-01-01',
POLARIZATION:'VVVH',
ORBIT : 'DESCENDING',
GEOMETRY: geometry,

// Speckle filter
APPLY_SPECKLE_FILTERING: true,
SPECKLE_FILTER_FRAMEWORK: 'MONO',
SPECKLE_FILTER: 'REFINED LEE',

// Output
FORMAT : 'DB',
CLIP_TO_ROI: true,
SAVE_ASSETS: false
};

// Preprocess the S1 collection
var output = wrapper.s1_preproc(parameters);
var s1 = output[0];
var s1_preprocess = output[1];

// Convert to DB and Add VV/VH ratio band
var speckleFiltered = s1_preprocess
.map(helper.lin_to_db);

// Mean is preferred for SAR data
var sarComposite = speckleFiltered.mean().select('VH')
var sarComposite = speckleFiltered.mean().select('VH');
var threshold = -25;
var sarWater = sarComposite.lt(threshold);
Map.addLayer(sarWater.clip(geometry).selfMask(), {min:0, max:1, palette: ['blue']}, 'Water')


// ***************************
// Speckle Filtering functions
// Credit: SERVIR-Mekong, adapted from
// https://mygeoblog.com/2021/01/21/sentinel-1-speckle-filter-refined-lee/
function powerToDb(img){
return ee.Image(10).multiply(img.log10());
}

function dbToPower(img){
return ee.Image(10).pow(img.divide(10));
}

// The RL speckle filter
function refinedLee(image) {
var date = image.date();
var bandNames = image.bandNames();
image = dbToPower(image);

var result = ee.ImageCollection(bandNames.map(function(b){
var img = image.select([b]);

// img must be in natural units, i.e. not in dB!
// Set up 3x3 kernels
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();
gradients = gradients.addBands(sample_mean.select(6).subtract(sample_mean.select(2)).abs());
gradients = gradients.addBands(sample_mean.select(3).subtract(sample_mean.select(5)).abs());
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
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
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)
directions = directions.reduce(ee.Reducer.sum());

//var pal = ['ffffff','ff0000','ffff00', '00ff00', '00ffff', '0000ff', 'ff00ff', '000000'];
//Map.addLayer(directions.reduce(ee.Reducer.sum()), {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));

dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)));
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);

return dir_mean.add(b.multiply(img.subtract(dir_mean)))
.arrayProject([0])
// Get a multi-band image bands.
.arrayFlatten([['sum']])
.float();
})).toBands().rename(bandNames);
var resultImage = powerToDb(ee.Image(result));
return resultImage.set('system:time_start', date.millis());
}

23 changes: 8 additions & 15 deletions code/gee_water_resources_management/Supplement/Surface_Water/Detecting_Algae_Blooms
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,3 @@
/**** Start of imports. If edited, may not auto-convert in the playground. ****/
var geometry =
/* color: #98ff00 */
/* shown: false */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[73.55833843987628, 21.30221866314927],
[73.55833843987628, 21.16269103558375],
[73.80278424065753, 21.16269103558375],
[73.80278424065753, 21.30221866314927]]], null, false);
/***** End of imports. If edited, may not auto-convert in the playground. *****/
// Example script showing how to compute
// Normalized Difference Chlorophyll Index (NDCI) to detect
// Algal Blooms in inland waters
Expand All @@ -21,6 +6,14 @@ var geometry =
// Ukai Dam Reservoir, Gujarat, India
// during February, 2022
// More info at https://x.com/rajbhagatt/status/1497212432195809280

var geometry = ee.Geometry.Polygon([[
[73.5583, 21.3022],
[73.5583, 21.1626],
[73.8027, 21.1626],
[73.8027, 21.3022]
]]);

var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');

var startDate = ee.Date.fromYMD(2022, 2, 15);
Expand Down
171 changes: 35 additions & 136 deletions docs/code/gee_water_resources_management/Supplement/SAR/SAR_Water_Detection
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,144 +7,43 @@ var geometry = basin.geometry();
Map.addLayer(basin, {color: 'grey'}, 'Arkavathy Sub Basin')
Map.centerObject(geometry, 12);

var s1 = ee.ImageCollection("COPERNICUS/S1_GRD")
var filtered = s1
// Filter to get images with VV and VH dual polarization.
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
// Change the pass to ASCENDING depending on your location
.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
.filterDate('2019-01-01', '2020-01-01')
.filterBounds(geometry)

// Apply Speckle Filter
var speckleFiltered = filtered.map(refinedLee);
// Use Sentinel-1 ARD Pre-Processing
// See https://github.com/adugnag/gee_s1_ard/
var wrapper = require('users/adugnagirma/gee_s1_ard:wrapper');
var helper = require('users/adugnagirma/gee_s1_ard:utilities');


var parameters = {
// Input
START_DATE: '2019-01-01',
STOP_DATE: '2020-01-01',
POLARIZATION:'VVVH',
ORBIT : 'DESCENDING',
GEOMETRY: geometry,

// Speckle filter
APPLY_SPECKLE_FILTERING: true,
SPECKLE_FILTER_FRAMEWORK: 'MONO',
SPECKLE_FILTER: 'REFINED LEE',

// Output
FORMAT : 'DB',
CLIP_TO_ROI: true,
SAVE_ASSETS: false
};

// Preprocess the S1 collection
var output = wrapper.s1_preproc(parameters);
var s1 = output[0];
var s1_preprocess = output[1];

// Convert to DB and Add VV/VH ratio band
var speckleFiltered = s1_preprocess
.map(helper.lin_to_db);

// Mean is preferred for SAR data
var sarComposite = speckleFiltered.mean().select('VH')
var sarComposite = speckleFiltered.mean().select('VH');
var threshold = -25;
var sarWater = sarComposite.lt(threshold);
Map.addLayer(sarWater.clip(geometry).selfMask(), {min:0, max:1, palette: ['blue']}, 'Water')


// ***************************
// Speckle Filtering functions
// Credit: SERVIR-Mekong, adapted from
// https://mygeoblog.com/2021/01/21/sentinel-1-speckle-filter-refined-lee/
function powerToDb(img){
return ee.Image(10).multiply(img.log10());
}

function dbToPower(img){
return ee.Image(10).pow(img.divide(10));
}

// The RL speckle filter
function refinedLee(image) {
var date = image.date();
var bandNames = image.bandNames();
image = dbToPower(image);

var result = ee.ImageCollection(bandNames.map(function(b){
var img = image.select([b]);

// img must be in natural units, i.e. not in dB!
// Set up 3x3 kernels
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();
gradients = gradients.addBands(sample_mean.select(6).subtract(sample_mean.select(2)).abs());
gradients = gradients.addBands(sample_mean.select(3).subtract(sample_mean.select(5)).abs());
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
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
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)
directions = directions.reduce(ee.Reducer.sum());

//var pal = ['ffffff','ff0000','ffff00', '00ff00', '00ffff', '0000ff', 'ff00ff', '000000'];
//Map.addLayer(directions.reduce(ee.Reducer.sum()), {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));

dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)));
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);

return dir_mean.add(b.multiply(img.subtract(dir_mean)))
.arrayProject([0])
// Get a multi-band image bands.
.arrayFlatten([['sum']])
.float();
})).toBands().rename(bandNames);
var resultImage = powerToDb(ee.Image(result));
return resultImage.set('system:time_start', date.millis());
}

Loading

0 comments on commit 0df435e

Please sign in to comment.