This filter retains image edges, fine details, and thin points
Code was implemented from here.
Find the earth engine code here
/*Copyright (c) 2021 SERVIR-Mekong
Permission is hereby granted, free of charge, to any person obtaining a copy
of the data and associated documentation files, to deal in the data
without restriction, including without limitation the rights to use, copy, modify,
merge, publish, distribute, sublicense, and/or sell copies, and to permit persons
to whom the data is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE DATA IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.*/
// Import Sentinel-1 Collection
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filterBounds(geometry)
.filterDate("2020-10-01","2020-10-31");
var firstNoTerrainCorrection = ee.Image(s1.first());
Map.addLayer(firstNoTerrainCorrection,{min:-25,max:20},"no terrain correction");
s1 = s1.map(terrainCorrection);
s1_lee_sigma = s1.map(leeSigma);
var firstTerrainCorrection = ee.Image(s1.first());
var s1_lee_sigma = ee.Image(s1_lee_sigma.first());
Map.addLayer(firstTerrainCorrection,{min:-25,max:20},"Terrain corrected");
Map.addLayer(s1_lee_sigma,{min:-25,max:20},"lee sigma");
// Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
function terrainCorrection(image) {
var imgGeom = image.geometry();
var srtm = ee.Image('USGS/SRTMGL1_003').clip(imgGeom); // 30m srtm
var sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0));
// Article ( numbers relate to chapters)
// 2.1.1 Radar geometry
var theta_i = image.select('angle');
var phi_i = ee.Terrain.aspect(theta_i)
.reduceRegion(ee.Reducer.mean(), theta_i.get('system:footprint'), 1000)
.get('aspect');
// 2.1.2 Terrain geometry
var alpha_s = ee.Terrain.slope(srtm).select('slope');
var phi_s = ee.Terrain.aspect(srtm).select('aspect');
// 2.1.3 Model geometry
// reduce to 3 angle
var phi_r = ee.Image.constant(phi_i).subtract(phi_s);
// convert all to radians
var phi_rRad = phi_r.multiply(Math.PI / 180);
var alpha_sRad = alpha_s.multiply(Math.PI / 180);
var theta_iRad = theta_i.multiply(Math.PI / 180);
var ninetyRad = ee.Image.constant(90).multiply(Math.PI / 180);
// slope steepness in range (eq. 2)
var alpha_r = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan();
// slope steepness in azimuth (eq 3)
var alpha_az = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan();
// local incidence angle (eq. 4)
var theta_lia = (alpha_az.cos().multiply((theta_iRad.subtract(alpha_r)).cos())).acos();
var theta_liaDeg = theta_lia.multiply(180 / Math.PI);
// 2.2
// Gamma_nought_flat
var gamma0 = sigma0Pow.divide(theta_iRad.cos());
var gamma0dB = ee.Image.constant(10).multiply(gamma0.log10());
var ratio_1 = gamma0dB.select('VV').subtract(gamma0dB.select('VH'));
// Volumetric Model
var nominator = (ninetyRad.subtract(theta_iRad).add(alpha_r)).tan();
var denominator = (ninetyRad.subtract(theta_iRad)).tan();
var volModel = (nominator.divide(denominator)).abs();
// apply model
var gamma0_Volume = gamma0.divide(volModel);
var gamma0_VolumeDB = ee.Image.constant(10).multiply(gamma0_Volume.log10());
// we add a layover/shadow maskto the original implmentation
// layover, where slope > radar viewing angle
var alpha_rDeg = alpha_r.multiply(180 / Math.PI);
var layover = alpha_rDeg.lt(theta_i);
// shadow where LIA > 90
var shadow = theta_liaDeg.lt(85);
// calculate the ratio for RGB vis
var ratio = gamma0_VolumeDB.select('VV').subtract(gamma0_VolumeDB.select('VH'));
var output = gamma0_VolumeDB.addBands(ratio).addBands(alpha_r).addBands(phi_s).addBands(theta_iRad)
.addBands(layover).addBands(shadow).addBands(gamma0dB).addBands(ratio_1);
return image.addBands(
output.select(['VV', 'VH'], ['VV', 'VH']),
null,
true
);
}
function powerToDb(img){
return ee.Image(10).multiply(img.log10());
}
function dbToPower(img){
return ee.Image(10).pow(img.divide(10));
}
function leeSigma(img){
var meanReducer = ee.Reducer.mean(),
varianceReducer = ee.Reducer.variance(),
sumReducer = ee.Reducer.sum();
var bandNames = img.bandNames();
var pow = dbToPower(img);
var kernelWeights = ee.List.repeat(ee.List.repeat(1,9),9);
var kernel = ee.Kernel.fixed(9,9,kernelWeights,5,5);
var targetWeights = ee.List.repeat(ee.List.repeat(1,3),3);
var targetKernel = ee.Kernel.fixed(3,3,targetWeights,1,1);
var Tk = ee.Image(7);
var sigma = '0.9';
var looks = '4';
// Lookup table for range and eta values for intensity
var sigmaLookup = ee.Dictionary({
1: ee.Dictionary({
0.5: ee.Dictionary({
'A1': 0.436,
'A2': 1.92,
'η': 0.4057
}),
0.6: ee.Dictionary({
'A1': 0.343,
'A2': 2.21,
'η': 0.4954
}),
0.7: ee.Dictionary({
'A1': 0.254,
'A2': 2.582,
'η': 0.5911
}),
0.8: ee.Dictionary({
'A1': 0.168,
'A2': 3.094,
'η': 0.6966
}),
0.9: ee.Dictionary({
'A1': 0.084,
'A2': 3.941,
'η': 0.8191
}),
0.95: ee.Dictionary({
'A1': 0.043,
'A2': 4.840,
'η': 0.8599
})
}),
2: ee.Dictionary({
0.5: ee.Dictionary({
'A1': 0.582,
'A2': 1.584,
'η': 0.2763
}),
0.6: ee.Dictionary({
'A1': 0.501,
'A2': 1.755,
'η': 0.3388
}),
0.7: ee.Dictionary({
'A1': 0.418,
'A2': 1.972,
'η': 0.4062
}),
0.8: ee.Dictionary({
'A1': 0.327,
'A2': 2.260,
'η': 0.4819
}),
0.9: ee.Dictionary({
'A1': 0.221,
'A2': 2.744,
'η': 0.5699
}),
0.95: ee.Dictionary({
'A1': 0.152,
'A2': 3.206,
'η': 0.6254
}),
}),
3: ee.Dictionary({
0.5: ee.Dictionary({
'A1': 0.652,
'A2': 1.458,
'η': 0.2222
}),
0.6: ee.Dictionary({
'A1': 0.580,
'A2': 1.586,
'η': 0.2736
}),
0.7: ee.Dictionary({
'A1': 0.505,
'A2': 1.751,
'η': 0.3280
}),
0.8: ee.Dictionary({
'A1': 0.419,
'A2': 1.865,
'η': 0.3892
}),
0.9: ee.Dictionary({
'A1': 0.313,
'A2': 2.320,
'η': 0.4624
}),
0.95: ee.Dictionary({
'A1': 0.238,
'A2': 2.656,
'η': 0.5084
}),
}),
4: ee.Dictionary({
0.5: ee.Dictionary({
'A1': 0.694,
'A2': 1.385,
'η': 0.1921
}),
0.6: ee.Dictionary({
'A1': 0.630,
'A2': 1.495,
'η': 0.2348
}),
0.7: ee.Dictionary({
'A1': 0.560,
'A2': 1.627,
'η': 0.2825
}),
0.8: ee.Dictionary({
'A1': 0.480,
'A2': 1.804,
'η': 0.3354
}),
0.9: ee.Dictionary({
'A1': 0.378,
'A2': 2.094,
'η': 0.3991
}),
0.95: ee.Dictionary({
'A1': 0.302,
'A2': 2.360,
'η': 0.4391
}),
})
});
// extract data from lookup
var looksDict = ee.Dictionary(sigmaLookup.get(ee.String(looks)));
var sigmaImage = ee.Dictionary(looksDict.get(ee.String(sigma))).toImage();
var a1 = sigmaImage.select('A1');
var a2 = sigmaImage.select('A2');
var aRange = a2.subtract(a1);
var eta = sigmaImage.select('η').pow(2);
// MMSE estimator
var mmseMask = pow.gte(a1).or(pow.lte(a2));
var mmseIn = pow.updateMask(mmseMask);
var oneImg = ee.Image(1);
var z = mmseIn.reduceNeighborhood(meanReducer,kernel,null,true);
var varz = mmseIn.reduceNeighborhood(varianceReducer,kernel);
var varx = (varz.subtract(z.abs().pow(2).multiply(eta))).divide(oneImg.add(eta));
var b = varx.divide(varz);
var mmse = oneImg.subtract(b).multiply(z.abs()).add(b.multiply(mmseIn));
// workflow
var z99 = ee.Dictionary(pow.reduceRegion({
reducer: ee.Reducer.percentile([99],null,255,0.001,1e6),
geometry: img.geometry(),
scale:10,
bestEffort:true
})).toImage();
var overThresh = pow.gte(z99);
var K = overThresh.reduceNeighborhood(sumReducer,targetKernel,null,true);
var retainPixel = K.gte(Tk);
var xHat = powerToDb(pow.updateMask(retainPixel).unmask(mmse));
return ee.Image(xHat).clip(img.geometry()).rename(bandNames)
.copyProperties(img,null,null)
.set('system:time_start',img.get('system:time_start'));
}
Open Geo Blog 