A modular land cover system part 4: brdf correction
Step 1: create a new script and call it brdf_module
Step 2: copy the code below and store it in the script
var PI = ee.Number(3.14159265359);
var MAX_SATELLITE_ZENITH = 7.5;
var MAX_DISTANCE = 1000000;
var UPPER_LEFT = 0;
var LOWER_LEFT = 1;
var LOWER_RIGHT = 2;
var UPPER_RIGHT = 3;
/* Creates a collection of mosaics with a given temporal interval.
*
* collection - the collection from which to make composites.
* start - the date of the first composite (either a string or an ee.Date)
* count - the number of composites to make
* interval - The time between composites, in units of "units".
* units - The units of step (day, week, month, year; see ee ee.Date.advance)
*/
exports.brdfS2 = function(collection) {
collection = collection.map(applyBRDF);
return collection;
function applyBRDF(image){
var date = image.date();
var footprint = ee.List(image.geometry().bounds().bounds().coordinates().get(0));
var angles = getsunAngles(date, footprint);
var sunAz = angles[0];
var sunZen = angles[1];
var viewAz = azimuth(footprint);
var viewZen = zenith(footprint);
var kval = _kvol(sunAz, sunZen, viewAz, viewZen);
var kvol = kval[0];
var kvol0 = kval[1];
var result = _apply(image, kvol.multiply(PI), kvol0.multiply(PI));
return result;}
/* Get sunAngles from the map given the data.
*
* date: ee.date object
* footprint: geometry of the image
*/
function getsunAngles(date, footprint){
var jdp = date.getFraction('year');
var seconds_in_hour = 3600;
var hourGMT = ee.Number(date.getRelative('second', 'day')).divide(seconds_in_hour);
var latRad = ee.Image.pixelLonLat().select('latitude').multiply(PI.divide(180));
var longDeg = ee.Image.pixelLonLat().select('longitude');
// Julian day proportion in radians
var jdpr = jdp.multiply(PI).multiply(2);
var a = ee.List([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]);
var meanSolarTime = longDeg.divide(15.0).add(ee.Number(hourGMT));
var localSolarDiff1 = value(a, 0)
.add(value(a, 1).multiply(jdpr.cos()))
.subtract(value(a, 2).multiply(jdpr.sin()))
.subtract(value(a, 3).multiply(jdpr.multiply(2).cos()))
.subtract(value(a, 4).multiply(jdpr.multiply(2).sin()));
var localSolarDiff2 = localSolarDiff1.multiply(12 * 60);
var localSolarDiff = localSolarDiff2.divide(PI);
var trueSolarTime = meanSolarTime
.add(localSolarDiff.divide(60))
.subtract(12.0);
// Hour as an angle;
var ah = trueSolarTime.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2).multiply(PI.divide(180))) ;
var b = ee.List([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]);
var delta = value(b, 0)
.subtract(value(b, 1).multiply(jdpr.cos()))
.add(value(b, 2).multiply(jdpr.sin()))
.subtract(value(b, 3).multiply(jdpr.multiply(2).cos()))
.add(value(b, 4).multiply(jdpr.multiply(2).sin()))
.subtract(value(b, 5).multiply(jdpr.multiply(3).cos()))
.add(value(b, 6).multiply(jdpr.multiply(3).sin()));
var cosSunZen = latRad.sin().multiply(delta.sin())
.add(latRad.cos().multiply(ah.cos()).multiply(delta.cos()));
var sunZen = cosSunZen.acos();
// sun azimuth from south, turning west
var sinSunAzSW = ah.sin().multiply(delta.cos()).divide(sunZen.sin());
sinSunAzSW = sinSunAzSW.clamp(-1.0, 1.0);
var cosSunAzSW = (latRad.cos().multiply(-1).multiply(delta.sin())
.add(latRad.sin().multiply(delta.cos()).multiply(ah.cos())))
.divide(sunZen.sin());
var sunAzSW = sinSunAzSW.asin();
sunAzSW = where(cosSunAzSW.lte(0), sunAzSW.multiply(-1).add(PI), sunAzSW);
sunAzSW = where(cosSunAzSW.gt(0).and(sinSunAzSW.lte(0)), sunAzSW.add(PI.multiply(2)), sunAzSW);
var sunAz = sunAzSW.add(PI);
// # Keep within [0, 2pi] range
sunAz = where(sunAz.gt(PI.multiply(2)), sunAz.subtract(PI.multiply(2)), sunAz);
var footprint_polygon = ee.Geometry.Polygon(footprint);
sunAz = sunAz.clip(footprint_polygon);
sunAz = sunAz.rename(['sunAz']);
sunZen = sunZen.clip(footprint_polygon).rename(['sunZen']);
return [sunAz, sunZen];
}
/* Get azimuth.
*
*
* footprint: geometry of the image
*/
function azimuth(footprint){
function x(point){return ee.Number(ee.List(point).get(0))}
function y(point){return ee.Number(ee.List(point).get(1))}
var upperCenter = line_from_coords(footprint, UPPER_LEFT, UPPER_RIGHT).centroid().coordinates();
var lowerCenter = line_from_coords(footprint, LOWER_LEFT, LOWER_RIGHT).centroid().coordinates();
var slope = ((y(lowerCenter)).subtract(y(upperCenter))).divide((x(lowerCenter)).subtract(x(upperCenter)));
var slopePerp = ee.Number(-1).divide(slope);
var azimuthLeft = ee.Image(PI.divide(2).subtract((slopePerp).atan()));
return azimuthLeft.rename(['viewAz']);
}
/* Get zenith.
*
*
* footprint: geometry of the image
*/
function zenith(footprint){
var leftLine = line_from_coords(footprint, UPPER_LEFT, LOWER_LEFT);
var rightLine = line_from_coords(footprint, UPPER_RIGHT, LOWER_RIGHT);
var leftDistance = ee.FeatureCollection(leftLine).distance(MAX_DISTANCE);
var rightDistance = ee.FeatureCollection(rightLine).distance(MAX_DISTANCE);
var viewZenith = rightDistance.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2))
.divide(rightDistance.add(leftDistance))
.subtract(ee.Number(MAX_SATELLITE_ZENITH))
.clip(ee.Geometry.Polygon(footprint))
.rename(['viewZen']);
return viewZenith.multiply(PI.divide(180));
}
/* apply function to all bands
*
* http://www.mdpi.com/2072-4292/9/12/1325/htm#sec3dot2-remotesensing-09-01325
* https://www.sciencedirect.com/science/article/pii/S0034425717302791
*
* image : the image to apply the function to
* kvol:
* kvol0
*
*/
function _apply(image, kvol, kvol0){
var f_iso = 0;
var f_geo = 0;
var f_vol = 0;
var blue = _correct_band(image, 'blue', kvol, kvol0, f_iso=0.0774, f_geo=0.0079, f_vol=0.0372);
var green = _correct_band(image, 'green', kvol, kvol0, f_iso=0.1306, f_geo=0.0178, f_vol=0.0580);
var red = _correct_band(image, 'red', kvol, kvol0, f_iso=0.1690, f_geo=0.0227, f_vol=0.0574);
var re1 = _correct_band(image, 're1', kvol, kvol0, f_iso=0.2085, f_geo=0.0256, f_vol=0.0845);
var re2 = _correct_band(image, 're2', kvol, kvol0, f_iso=0.2316, f_geo=0.0273, f_vol=0.1003);
var re3 = _correct_band(image, 're3', kvol, kvol0, f_iso=0.2599, f_geo=0.0294, f_vol=0.1197);
var nir = _correct_band(image, 'nir', kvol, kvol0, f_iso=0.3093, f_geo=0.0330, f_vol=0.1535);
var re4 = _correct_band(image, 're4', kvol, kvol0, f_iso=0.2907, f_geo=0.0410, f_vol=0.1611);
var swir1 = _correct_band(image, 'swir1', kvol, kvol0, f_iso=0.3430, f_geo=0.0453, f_vol=0.1154);
var swir2 = _correct_band(image, 'swir2', kvol, kvol0, f_iso=0.2658, f_geo=0.0387, f_vol=0.0639);
return image.select([]).addBands([blue, green, red, nir,re1,re2,re3,nir,re4,swir1, swir2]);
}
/* correct band function
*
*
* image : the image to apply the function to
* band_name
* kvol
* kvol0
* f_iso
* f_geo
* f_vol
*
*/
function _correct_band(image, band_name, kvol, kvol0, f_iso, f_geo, f_vol){
//"""fiso + fvol * kvol + fgeo * kgeo"""
var iso = ee.Image(f_iso);
var geo = ee.Image(f_geo);
var vol = ee.Image(f_vol);
var pred = vol.multiply(kvol).add(geo.multiply(kvol)).add(iso).rename(['pred']);
var pred0 = vol.multiply(kvol0).add(geo.multiply(kvol0)).add(iso).rename(['pred0']);
var cfac = pred0.divide(pred).rename(['cfac']);
var corr = image.select(band_name).multiply(cfac).rename([band_name]);
return corr;
}
/* calculate kvol and kvol0
*
* sunAZ
* sunZen
* viewAz
* viewZen
*
*/
function _kvol(sunAz, sunZen, viewAz, viewZen){
//"""Calculate kvol kernel.
//From Lucht et al. 2000
//Phase angle = cos(solar zenith) cos(view zenith) + sin(solar zenith) sin(view zenith) cos(relative azimuth)"""
var relative_azimuth = sunAz.subtract(viewAz).rename(['relAz']);
var pa1 = viewZen.cos().multiply(sunZen.cos());
var pa2 = viewZen.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
var phase_angle1 = pa1.add(pa2);
var phase_angle = phase_angle1.acos();
var p1 = ee.Image(PI.divide(2)).subtract(phase_angle);
var p2 = p1.multiply(phase_angle1);
var p3 = p2.add(phase_angle.sin());
var p4 = sunZen.cos().add(viewZen.cos());
var p5 = ee.Image(PI.divide(4));
var kvol = p3.divide(p4).subtract(p5).rename(['kvol']);
var viewZen0 = ee.Image(0);
var pa10 = viewZen0.cos().multiply(sunZen.cos());
var pa20 = viewZen0.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
var phase_angle10 = pa10.add(pa20);
var phase_angle0 = phase_angle10.acos();
var p10 = ee.Image(PI.divide(2)).subtract(phase_angle0);
var p20 = p10.multiply(phase_angle10);
var p30 = p20.add(phase_angle0.sin());
var p40 = sunZen.cos().add(viewZen0.cos());
var p50 = ee.Image(PI.divide(4));
var kvol0 = p30.divide(p40).subtract(p50).rename(['kvol0']);
return [kvol, kvol0]}
/* helper function
*
*
*
*/
function line_from_coords(coordinates, fromIndex, toIndex){
return ee.Geometry.LineString(ee.List([
coordinates.get(fromIndex),
coordinates.get(toIndex)]));
}
function where(condition, trueValue, falseValue){
var trueMasked = trueValue.mask(condition);
var falseMasked = falseValue.mask(invertMask(condition));
return trueMasked.unmask(falseMasked);
}
function invertMask(mask){
return mask.multiply(-1).add(1);
}
function value(list,index){
return ee.Number(list.get(index));
}
};
/* Creates a collection of mosaics with a given temporal interval.
*
* collection - the collection from which to make composites.
* start - the date of the first composite (either a string or an ee.Date)
* count - the number of composites to make
* interval - The time between composites, in units of "units".
* units - The units of step (day, week, month, year; see ee ee.Date.advance)
*/
exports.brdfL8 = function(image) {
var date = image.date();
var footprint = ee.List(image.geometry().bounds().bounds().coordinates().get(0));
var angles = getsunAngles(date, footprint);
var sunAz = angles[0];
var sunZen = angles[1];
var viewAz = azimuth(footprint);
var viewZen = zenith(footprint);
var kval = _kvol(sunAz, sunZen, viewAz, viewZen);
var kvol = kval[0];
var kvol0 = kval[1];
var result = _apply(image, kvol.multiply(PI), kvol0.multiply(PI));
return result;
/* Get sunAngles from the map given the data.
*
* date: ee.date object
* footprint: geometry of the image
*/
function getsunAngles(date, footprint){
var jdp = date.getFraction('year');
var seconds_in_hour = 3600;
var hourGMT = ee.Number(date.getRelative('second', 'day')).divide(seconds_in_hour);
var latRad = ee.Image.pixelLonLat().select('latitude').multiply(PI.divide(180));
var longDeg = ee.Image.pixelLonLat().select('longitude');
// Julian day proportion in radians
var jdpr = jdp.multiply(PI).multiply(2);
var a = ee.List([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]);
var meanSolarTime = longDeg.divide(15.0).add(ee.Number(hourGMT));
var localSolarDiff1 = value(a, 0)
.add(value(a, 1).multiply(jdpr.cos()))
.subtract(value(a, 2).multiply(jdpr.sin()))
.subtract(value(a, 3).multiply(jdpr.multiply(2).cos()))
.subtract(value(a, 4).multiply(jdpr.multiply(2).sin()));
var localSolarDiff2 = localSolarDiff1.multiply(12 * 60);
var localSolarDiff = localSolarDiff2.divide(PI);
var trueSolarTime = meanSolarTime
.add(localSolarDiff.divide(60))
.subtract(12.0);
// Hour as an angle;
var ah = trueSolarTime.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2).multiply(PI.divide(180))) ;
var b = ee.List([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]);
var delta = value(b, 0)
.subtract(value(b, 1).multiply(jdpr.cos()))
.add(value(b, 2).multiply(jdpr.sin()))
.subtract(value(b, 3).multiply(jdpr.multiply(2).cos()))
.add(value(b, 4).multiply(jdpr.multiply(2).sin()))
.subtract(value(b, 5).multiply(jdpr.multiply(3).cos()))
.add(value(b, 6).multiply(jdpr.multiply(3).sin()));
var cosSunZen = latRad.sin().multiply(delta.sin())
.add(latRad.cos().multiply(ah.cos()).multiply(delta.cos()));
var sunZen = cosSunZen.acos();
// sun azimuth from south, turning west
var sinSunAzSW = ah.sin().multiply(delta.cos()).divide(sunZen.sin());
sinSunAzSW = sinSunAzSW.clamp(-1.0, 1.0);
var cosSunAzSW = (latRad.cos().multiply(-1).multiply(delta.sin())
.add(latRad.sin().multiply(delta.cos()).multiply(ah.cos())))
.divide(sunZen.sin());
var sunAzSW = sinSunAzSW.asin();
sunAzSW = where(cosSunAzSW.lte(0), sunAzSW.multiply(-1).add(PI), sunAzSW);
sunAzSW = where(cosSunAzSW.gt(0).and(sinSunAzSW.lte(0)), sunAzSW.add(PI.multiply(2)), sunAzSW);
var sunAz = sunAzSW.add(PI);
// # Keep within [0, 2pi] range
sunAz = where(sunAz.gt(PI.multiply(2)), sunAz.subtract(PI.multiply(2)), sunAz);
var footprint_polygon = ee.Geometry.Polygon(footprint);
sunAz = sunAz.clip(footprint_polygon);
sunAz = sunAz.rename(['sunAz']);
sunZen = sunZen.clip(footprint_polygon).rename(['sunZen']);
return [sunAz, sunZen];
}
/* Get azimuth.
*
*
* footprint: geometry of the image
*/
function azimuth(footprint){
function x(point){return ee.Number(ee.List(point).get(0))}
function y(point){return ee.Number(ee.List(point).get(1))}
var upperCenter = line_from_coords(footprint, UPPER_LEFT, UPPER_RIGHT).centroid().coordinates();
var lowerCenter = line_from_coords(footprint, LOWER_LEFT, LOWER_RIGHT).centroid().coordinates();
var slope = ((y(lowerCenter)).subtract(y(upperCenter))).divide((x(lowerCenter)).subtract(x(upperCenter)));
var slopePerp = ee.Number(-1).divide(slope);
var azimuthLeft = ee.Image(PI.divide(2).subtract((slopePerp).atan()));
return azimuthLeft.rename(['viewAz']);
}
/* Get zenith.
*
*
* footprint: geometry of the image
*/
function zenith(footprint){
var leftLine = line_from_coords(footprint, UPPER_LEFT, LOWER_LEFT);
var rightLine = line_from_coords(footprint, UPPER_RIGHT, LOWER_RIGHT);
var leftDistance = ee.FeatureCollection(leftLine).distance(MAX_DISTANCE);
var rightDistance = ee.FeatureCollection(rightLine).distance(MAX_DISTANCE);
var viewZenith = rightDistance.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2))
.divide(rightDistance.add(leftDistance))
.subtract(ee.Number(MAX_SATELLITE_ZENITH))
.clip(ee.Geometry.Polygon(footprint))
.rename(['viewZen']);
return viewZenith.multiply(PI.divide(180));
}
/* apply function to all bands
*
* http://www.mdpi.com/2072-4292/9/12/1325/htm#sec3dot2-remotesensing-09-01325
* https://www.sciencedirect.com/science/article/pii/S0034425717302791
*
* image : the image to apply the function to
* kvol:
* kvol0
*
*/
function _apply(image, kvol, kvol0){
var f_iso = 0;
var f_geo = 0;
var f_vol = 0;
var blue = _correct_band(image, 'blue', kvol, kvol0, f_iso=0.0774, f_geo=0.0079, f_vol=0.0372);
var green = _correct_band(image, 'green', kvol, kvol0, f_iso=0.1306, f_geo=0.0178, f_vol=0.0580);
var red = _correct_band(image, 'red', kvol, kvol0, f_iso=0.1690, f_geo=0.0227, f_vol=0.0574);
var nir = _correct_band(image, 'nir', kvol, kvol0, f_iso=0.3093, f_geo=0.0330, f_vol=0.1535);
var swir1 = _correct_band(image, 'swir1', kvol, kvol0, f_iso=0.3430, f_geo=0.0453, f_vol=0.1154);
var swir2 = _correct_band(image, 'swir2', kvol, kvol0, f_iso=0.2658, f_geo=0.0387, f_vol=0.0639);
return image.select([]).addBands([blue, green, red, nir, swir1, swir2]);
}
/* correct band function
*
*
* image : the image to apply the function to
* band_name
* kvol
* kvol0
* f_iso
* f_geo
* f_vol
*
*/
function _correct_band(image, band_name, kvol, kvol0, f_iso, f_geo, f_vol){
//"""fiso + fvol * kvol + fgeo * kgeo"""
var iso = ee.Image(f_iso);
var geo = ee.Image(f_geo);
var vol = ee.Image(f_vol);
var pred = vol.multiply(kvol).add(geo.multiply(kvol)).add(iso).rename(['pred']);
var pred0 = vol.multiply(kvol0).add(geo.multiply(kvol0)).add(iso).rename(['pred0']);
var cfac = pred0.divide(pred).rename(['cfac']);
var corr = image.select(band_name).multiply(cfac).rename([band_name]);
return corr;
}
/* calculate kvol and kvol0
*
* sunAZ
* sunZen
* viewAz
* viewZen
*
*/
function _kvol(sunAz, sunZen, viewAz, viewZen){
//"""Calculate kvol kernel.
//From Lucht et al. 2000
//Phase angle = cos(solar zenith) cos(view zenith) + sin(solar zenith) sin(view zenith) cos(relative azimuth)"""
var relative_azimuth = sunAz.subtract(viewAz).rename(['relAz']);
var pa1 = viewZen.cos().multiply(sunZen.cos());
var pa2 = viewZen.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
var phase_angle1 = pa1.add(pa2);
var phase_angle = phase_angle1.acos();
var p1 = ee.Image(PI.divide(2)).subtract(phase_angle);
var p2 = p1.multiply(phase_angle1);
var p3 = p2.add(phase_angle.sin());
var p4 = sunZen.cos().add(viewZen.cos());
var p5 = ee.Image(PI.divide(4));
var kvol = p3.divide(p4).subtract(p5).rename(['kvol']);
var viewZen0 = ee.Image(0);
var pa10 = viewZen0.cos().multiply(sunZen.cos());
var pa20 = viewZen0.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
var phase_angle10 = pa10.add(pa20);
var phase_angle0 = phase_angle10.acos();
var p10 = ee.Image(PI.divide(2)).subtract(phase_angle0);
var p20 = p10.multiply(phase_angle10);
var p30 = p20.add(phase_angle0.sin());
var p40 = sunZen.cos().add(viewZen0.cos());
var p50 = ee.Image(PI.divide(4));
var kvol0 = p30.divide(p40).subtract(p50).rename(['kvol0']);
return [kvol, kvol0]}
/* helper function
*
*
*
*/
function line_from_coords(coordinates, fromIndex, toIndex){
return ee.Geometry.LineString(ee.List([
coordinates.get(fromIndex),
coordinates.get(toIndex)]));
}
function where(condition, trueValue, falseValue){
var trueMasked = trueValue.mask(condition);
var falseMasked = falseValue.mask(invertMask(condition));
return trueMasked.unmask(falseMasked);
}
function invertMask(mask){
return mask.multiply(-1).add(1);
}
function value(list,index){
return ee.Number(list.get(index));
}
};
Step 4: Open the main script and add the lines below
var brdf = require("users/.../mrc:brdf_module");
print("applying sentinel brdf");
s2 = brdf.brdfS2(s2);
Map.addLayer(ee.Image(s2.??()),{min:0,max:0.??,bands:"swir2,nir,red"},"after brdf");
print(ee.Image(s2.first()));
Step 5: Run the script
Open Geo Blog 
Appreciate for this great script! Could this script be applied to Landsat 5/7/8?
LikeLike
yes!
LikeLike