I am working on a script that needs to do the following; 1. smooth a Landsat time series 2a. apply a filter to remove noise from sentinel-1 using refined Lee 2b. smooth a sentinel-1 time series 3. Join collections using a ee.Join.saveAll 4. create a series of functions to create yearly mosaics and export on a yearly basis.
This post has been edited to include a new issue in the same script. I have added the new function I am incorporating which was referenced from this page. The issue is that after the refined lee filtering is applied to S-1, I have the issue of missing properties. I attempted to fix that with a copyProperties() function, in order to retrieve the system time start property which i need for the following data fusion. The code error is that when I try to retrieve that property, the variable var = img is not defined althought it lives within the function.Here is my reproducible example,
/*1A.Sentinel-1 Functions////////////////////////////
Code adapted from :
-https://mygeoblog.com/2021/01/21/sentinel-1-speckle-filter-refined-lee/
-https://groups.google.com/g/google-earth-engine-developers/c/ExepnAmP-hQ/m/7e5DnjXXAQAJ
- 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
);
}
// Functions to convert from/to dB
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 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);
return powerToDb(ee.Image(result)).copyProperties(img, 'system:time_start');
}
//Indices Funciton
function addIndices(image){
var ndvi = image.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI')
var ndmi = image.normalizedDifference(['SR_B4', 'SR_B5']).rename('NDMI')
var ndwi = image.normalizedDifference(['SR_B2', 'SR_B4']).rename('NDWI')
return image.addBands(ndvi).addBands(ndmi).addBands(ndwi)
}
//Function for acquiring Landsat SR image collection
function getLandsatImageCollection(studyArea,startDate,endDate,startJulian,endJulian){
var ls;var l5SR;var l7SR;var l8SR;var l9SR;var out;
var sensorBandDictLandsatSR = ee.Dictionary({L9: ee.List([0,1,2,3,4,5,6,8,17,18]),
L8 : ee.List([0,1,2,3,4,5,6,8,17,18]),
L7 : ee.List([0,1,2,3,4,5,8,17,18]),
L5 : ee.List([0,1,2,3,4,5,8,17,18]),
L4 : ee.List([0,1,2,3,4,5,8,17,18])
});
var bandNamesLandsatSR = ee.Dictionary({L89: ee.List(['SR_B1_UB', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'ST_B6', 'QA_PIXEL', 'QA_RADSAT']),
L57: ee.List(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'ST_B6', 'QA_PIXEL', 'QA_RADSAT'])
})
//ee.List(['blue','green','red','nir','swir1',
// 'swir2','TEMP','QA_PIXEL', 'QA_RADSAT']);
// l4SR = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2")
// .filterDate(startDate,endDate)
// .filter(ee.Filter.calendarRange(startJulian,endJulian))
// .filterBounds(studyArea)
// .select(sensorBandDictLandsatSR.get('L4'), bandNamesLandsatSR.get('L57'));
l5SR = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
.filterDate(startDate,endDate)
.filter(ee.Filter.calendarRange(startJulian,endJulian))
.filterBounds(studyArea)
.select(sensorBandDictLandsatSR.get('L5'), bandNamesLandsatSR.get('L57'));
l7SR = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
.filterDate(startDate,endDate)
.filter(ee.Filter.calendarRange(startJulian,endJulian))
.filterBounds(studyArea)
.select(sensorBandDictLandsatSR.get('L7'), bandNamesLandsatSR.get('L57'));
l8SR = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterDate(startDate,endDate)
.filter(ee.Filter.calendarRange(startJulian,endJulian))
.filterBounds(studyArea)
.select(sensorBandDictLandsatSR.get('L8'), bandNamesLandsatSR.get('L89'));
l9SR = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
.filterDate(startDate,endDate)
.filter(ee.Filter.calendarRange(startJulian,endJulian))
.filterBounds(studyArea)
.select(sensorBandDictLandsatSR.get('L9'), bandNamesLandsatSR.get('L89'));
ls = ee.ImageCollection(l5SR.merge(l7SR).merge(l8SR).merge(l9SR))
.map(function(i){
var qaMask = i.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
var saturationMask = i.select('QA_RADSAT').eq(0);
var optical = i.select('^SR_B.*').multiply(0.0000275).add(-0.2);
var thermal = i.select('^ST_B.*').multiply(0.00341802).add(149.0);
return i.addBands(optical, null, true).addBands(thermal, null, true).updateMask(qaMask).updateMask(saturationMask).clip(studyArea);
});
out = ls.set('system:time_start', ls.get('system:time_start')) ;
return out.select(['^SR_B.*','^ST_B.*']);
}
function gets1(studyArea,startDate,endDate,startJulian,endJulian){
var s1 = ee.ImageCollection("COPERNICUS/S1_GRD")
.filterDate(startDate,endDate)
.filter(ee.Filter.calendarRange(startJulian,endJulian))
.filterBounds(studyArea)
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.map(function(image) {
var edge = image.lt(-30.0);
var maskedImage = image.mask().and(edge.not());
return image.updateMask(maskedImage);
});
return ee.ImageCollection(s1).select(['VV','VH']);
}
/////////////////////////////////////////////////////
//1B. USER INPUTS
var AssetIDRoot = 'projects/ee-pbaltezar91/assets/MIT2023/RASTER/BENIN_LSR30M_S1_2000_2022_Version01'////A. CHANGE THIS TO YOUR PREFERRED ASSET ID ROOT
var countryName = 'BEN'//B. CHANGE TO THE COUNTRY NAME YOU ARE WORKING WITH
var version = 'Version.01'//C. CHANGE BASED ON DATA PRODUCT VERSIONS
var author = 'BaltezarP_20230721'//D. ADD YOUR NAME AND DATE: LASTFIRST _YDM
var roi = roi //E. Change roi to your own digitized OR IMPORTED roi geometry
var libraryScript = "https://code.earthengine.google.com/bcd86a5d2f6c54a1a771f24f732cb28c"//F.IMPORTANT: MUST SAVE THIS SCRIPT LINK IN YOUR ASSETS, OTHERWISE WONT WORK
var script = "https://code.earthengine.google.com/0f921cb6a0af9a2129316adee42275c2"//G. CHANGE TO THE MOST CURRENT SCRPT BEFORE EXPORTING DATA PRODUCT
var studyArea = ee.FeatureCollection("WWF/HydroATLAS/v1/Basins/level12").filterBounds(roi)//H. CHANGE TO YOUR DESIRED SATELLITE IMAGERY FILTERING GEOMETRY AREA, CURRENTLY USING HYDROLOGY OF BENIN
Map.addLayer(ee.Image().byte().paint(studyArea), {palette:'white', opacity:0.5}, "Study Area")
Map.addLayer(ee.Image().byte().paint(studyArea.geometry().bounds(), 3,3), {palette:'white'}, "Bounding Area")
var startYear = 2015 //I. CHANGE TO YOUR YEAR OF INTEREST FOR START AND END YEAR AND JULIAN ; the earliest year available for Sentinel-1 is 2015
var endYear = 2022 //
///Gather imagery that is from September to December, peak season for max hyacynth extent
///althought peak season is slated till November, greenkeeper claims there is still extensive
///harvesting
var startJulian = 244//first day of september julian day
var endJulian = 364//last day of december julian day
if(startJulian > endJulian){endJulian = endJulian + 365}
var startDate = ee.Date.fromYMD(startYear,1,1).advance(startJulian,'day');
var endDate = ee.Date.fromYMD(endYear,1,1).advance(endJulian,'day');
print(startDate, endDate);
var exportName = countryName;
var crs = 'EPSG:4326';//J. CHANGE TO YOUR DESIRED PROJECTION, BUT GEE HANDLES PROJECTS AS WELL
Map.centerObject(studyArea, 12);
Map.setOptions('satellite');
//==================================================================================
//2. PROCESSING FUNCTIONS
//K. IMPORTANT!!!!:CHANGE THIS TO WHERE YOUR LIBRARY SCRIPT LIVES (SEE VAR LIBRARY SCRIPT)
//EXAMPLE: 'users/YOUR_USER_NAME/YOUR_PARENT_FOLDER:YOUR_LIBRARY_SCRIPT_NAME'
///THE CODE WILL NOT WORK IF YOU BREAK FUNCTIONS THAT EXIST IN LIBRARY OR IF YOU DO NOT DO THIS
//External Library Functions
//var constant = require('users/pbaltezar91/MIT_Spring2023:A_LandsatTimeSeries_Library');
var Sentinel_1_notworking = gets1(studyArea,startDate,endDate,startJulian,endJulian).map(RefinedLee)//Refined Lee did not work, need to trouble shoot in following version
print(Sentinel_1_notworking,'s-1 not copying property')
var Sentinel_1 = gets1(studyArea,startDate,endDate,startJulian,endJulian)//.map(RefinedLee)//Refined Lee did not work, need to trouble shoot in following version
Map.addLayer(Sentinel_1, {gamma:1.00, min:-24, max:10}, "Sentinel-1 Collection", false);
print(Sentinel_1, "Sentinel-1 Collection")
//Landsat Collection 2 Tier 1 SR Time Series
var Landsat = getLandsatImageCollection(studyArea,startDate,endDate,startJulian,endJulian).map(addIndices)
print(Landsat.first(), "Landsat Image Collection ")
// /////////////join sentinel-1 to Landsat 2015 - 2022
// Moving-Window Smoothing function
var days = 15
// Specify the time-window
// Convert to milliseconds
var millis = ee.Number(days).multiply(1000*60*60*24);
// We use a 'save-all join' to find all images
// that are within the time-window
// This filter will match all images that are captured
// within the specified day of the source image
var diffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
});
// Each image in the joined collection will contain
// matching images in the 'images' property
// Extract and return the mean of matched images
var extractAndComputeMedian = function(image) {
var matchingImages = ee.ImageCollection.fromImages(image.get('images'));
var medianImage = matchingImages.reduce(
ee.Reducer.median().setOutputs(['moving_median']))
return ee.Image(image).addBands(medianImage)
}
var smoothedCollection = ee.ImageCollection(
ee.Join.saveAll({matchesKey:'images'}).apply({
primary: Landsat,
secondary: Landsat,
condition: diffFilter
}).map(extractAndComputeMedian)).select('.*_moving_median$');
var smoothedCollection = ee.ImageCollection(
ee.Join.saveAll({matchesKey:'images'}).apply({
primary: smoothedCollection,
secondary: Sentinel_1,
condition: diffFilter
}).map(extractAndComputeMedian))
Map.addLayer(smoothedCollection, {bands:['nir_moving_median', 'swir1_moving_median', 'red_moving_median'], gamma:1.00,min:0.015,max:0.35}, "Landsat Collection + Sentinel-1", false)
/////Make a list of years to generate composites for.
var yearList = ee.List.sequence(startYear, endYear);
var fracToms = function(frac){
var fyear = ee.Number(frac);
var year = fyear.floor();
var d = fyear.subtract(year).multiply(365);
var day_one = ee.Date.fromYMD(year, 1, 1);
return day_one.advance(d, 'day').millis();
};
var count_mask = function(image) {
var MOD = image.select('nir_moving_median').neq(0)
return image.select('nir_moving_median').mask( MOD )//.multiply(10000).toShort();
};
print(smoothedCollection,"Collection Before Mapping")
// Map over the list of years to generate a composite for each year.
var yearCompList = yearList.map(function(year){
var yearCol = smoothedCollection.filter(ee.Filter.calendarRange(year, year, 'year'));
var n_images = yearCol.size();
var imgList = yearCol.aggregate_array('LANDSAT_SCENE_ID');
//var yearComp = yearCol.reduce(ee.Reducer.median())
var yearComp = yearCol.mosaic().addBands(yearCol.map(count_mask).count().divide(10000).rename('count'))
//metadata
var nBands = yearComp.bandNames().size();
return yearComp
.set({
'year': year,
'system:time_start': fracToms(year),
'image_list': imgList,
'n_bands': nBands,
'n_images': n_images,
'version': version,
'author_name':author,
'scriptLink':script
})
})
print(yearCompCol, "Image collection after year mapping, where it fails 1")
// Convert the annual composite image list to an ImageCollection
var yearCompCol = ee.ImageCollection.fromImages(yearCompList);
// Filter out years with no bands.
//(can happen if there were no images to composite)
yearCompCol = yearCompCol.filter(ee.Filter.gt('n_images', 0));
print(yearCompCol, "Image collection after year mapping, where it fails 2")
yearCompCol = yearCompCol.map(function(img){
return img
//.clip(clipGeometry)
//.multiply(10000).toShort()//Save space rather than float(), convert back to float post-processing?
.set('year', img.get('year'))
.set('n_images', img.get('n_images'))
.set('image_list', img.get('image_list'))
.set('version', img.get('version'))
.set('script', img.get('scriptLink'))
.set('author', img.get('author_name'));
});
//Export each yearly median composite using all bands
yearCompCol
.aggregate_array('year')
.evaluate(function (systemIndexes) {
systemIndexes.forEach(function (systemIndex, i) {
var image = yearCompCol
.filterMetadata('year', 'equals', systemIndex)
.first();
Export.image.toAsset({
image: image,
description: countryName +'_'+ systemIndex,
assetId: AssetIDRoot +countryName +'_' + systemIndex,
region: exportGeometry,
scale: 30,
crs: 'EPSG:4326',
maxPixels: 4e10
});
});
});
but the final scripts will have additional processing. Note the printed errors where it shows the collections were not properly constructed hence them being undefined. My final goal is to export to the final Image Collection asset in my folders.
This question is similar to this post I created, but this question seemed different.
