diff --git a/sentinel2.py b/sentinel2.py new file mode 100644 index 0000000..abde8e5 --- /dev/null +++ b/sentinel2.py @@ -0,0 +1,915 @@ +#sentinel-2 package + +import ee +from Py6S import * +import math +import datetime +import os, sys +from utils import * +sys.path.append("/gee-atmcorr-S2/bin/") +from atmospheric import Atmospheric +import sun_angles +import view_angles +import time + +class env(object): + + def __init__(self): + """Initialize the environment.""" + + # Initialize the Earth Engine object, using the authentication credentials. + ee.Initialize() + + self.dem = ee.Image("JAXA/ALOS/AW3D30_V1_1").select(["AVE"]) + self.epsg = "EPSG:4326" + self.feature = 0 + ########################################## + # variable for the getSentinel algorithm # + ########################################## + + self.metadataCloudCoverMax = 30 + + + ########################################## + # variable for the shadowMask algorithm # + ########################################## + + # zScoreThresh: Threshold for cloud shadow masking- lower number masks out + # less. Between -0.8 and -1.2 generally works well + self.zScoreThresh = -1 + + # shadowSumThresh: Sum of IR bands to include as shadows within TDOM and the + # shadow shift method (lower number masks out less) + self.shadowSumThresh = 0.10; + + # contractPixels: The radius of the number of pixels to contract (negative buffer) clouds and cloud shadows by. Intended to eliminate smaller cloud + # patches that are likely errors (1.5 results in a -1 pixel buffer)(0.5 results in a -0 pixel buffer) + # (1.5 or 2.5 generally is sufficient) + self.contractPixels = 1.5; + + # dilatePixels: The radius of the number of pixels to dilate (buffer) clouds + # and cloud shadows by. Intended to include edges of clouds/cloud shadows + # that are often missed (1.5 results in a 1 pixel buffer)(0.5 results in a 0 pixel buffer) + # (2.5 or 3.5 generally is sufficient) + self.dilatePixels = 2.5; + + + ########################################## + # variable for cloudScore algorithm # + ########################################## + + # 9. Cloud and cloud shadow masking parameters. + # If cloudScoreTDOM is chosen + # cloudScoreThresh: If using the cloudScoreTDOMShift method-Threshold for cloud + # masking (lower number masks more clouds. Between 10 and 30 generally works best) + self.cloudScoreThresh = 35; + + # Percentile of cloud score to pull from time series to represent a minimum for + # the cloud score over time for a given pixel. Reduces commission errors over + # cool bright surfaces. Generally between 5 and 10 works well. 0 generally is a bit noisy + self.cloudScorePctl = 5; + + ########################################## + # variable for terrain algorithm # + ########################################## + + self.terrainScale = 3000 + + ########################################## + # Export variables # + ########################################## + + self.assetId ="users/servirmekong/rubberCollection/" + self.name = "1S2_2018_winter_final" + self.exportScale = 20 + + ########################################## + # variable band selection # + ########################################## + + self.s2BandsIn = ee.List(['QA60','B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12']) + self.s2BandsOut = ee.List(['QA60','cb','blue','green','red','re1','re2','re3','nir1','nir2','waterVapor','cirrus','swir1','swir2']) + self.divideBands = ee.List(['blue','green','red','re1','re2','re3','nir1','nir2','cb','cirrus','swir1','swir2','waterVapor']) + self.medianIncludeBands = ee.List(['blue','green','red','re1','re2','re3','nir1','nir2','cb','cirrus','swir1','swir2','waterVapor']) + self.medoidBands = ee.List(['blue','green','red','re1','re2','re3','nir1','nir2','cb','cirrus','swir1','swir2','waterVapor']) + + ########################################## + # enable / disable modules # + ########################################## + self.calcSR = True + self.brdf = False + self.QAcloudMask = True + self.cloudMask = True + self.shadowMask = True + self.terrainCorrection = False + self.percentiles = [20,80] + +class functions(): + def __init__(self): + """Initialize the Surfrace Reflectance app.""" + + # get the environment + self.env = env() + + + def main(self,studyArea,startDate,endDate,startDay,endDay,week): + + self.env.startDate = startDate + self.env.endDate = endDate + + self.env.startDoy = startDay + self.env.endDoy = endDay + + s2 = self.getSentinel2(startDate,endDate,studyArea); + + print(s2.size().getInfo()) + + if s2.size().getInfo() > 0: + + print(self.env.startDate.getInfo()) + print(self.env.endDate.getInfo()) + + print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + + s2 = s2.map(self.scaleS2) + print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + # masking the shadows + print("Masking shadows..") + if self.env.shadowMask == True: + s2 = self.maskShadows(s2,studyArea) + + + #print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + self.collectionMeta = s2.getInfo()['features'] + print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + + #print("scaling bands..") + #s2 = s2.map(self.scaleS2) + #print(ee.Image(s2.first()).get('system:footprint').getInfo()) + + # applying the atmospheric correction + #if self.env.calcSR == True: + # print("applying atmospheric correction") + # s2 = s2.map(self.TOAtoSR) #.select(self.env.s2BandsIn,self.env.s2BandsOut) + #print(ee.Image(s2.first()).get('system:time_start').getInfo()) + s2 = s2.select(self.env.s2BandsIn,self.env.s2BandsOut) + + + if self.env.QAcloudMask == True: + print("use QA band for cloud Masking") + s2 = s2.map(self.QAMaskCloud) + print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + + + if self.env.cloudMask == True: + print("sentinel cloud score...") + s2 = s2.map(self.sentinelCloudScore) + s2 = self.cloudMasking(s2) + print(ee.Image(s2.first()).get('SOLAR_IRRADIANCE_B1').getInfo()) + + + s2 = s2.map(self.removeHaze) + s2 = s2.select(self.env.s2BandsOut,self.env.s2BandsIn) + + #s2 = s2.map(self.pixelArea) + #s2 = s2.filter(ee.Filter.gt("pixelArea",7000)) + print("-------") + #print(s2.size().getInfo()) + + # applying the atmospheric correction + if self.env.calcSR == True: + print("applying atmospheric correction") + s2 = s2.map(self.TOAtoSR) #.select(self.env.s2BandsIn,self.env.s2BandsOut) + #print(ee.Image(s2.first()).get('system:time_start').getInfo()) + s2 = s2.select(self.env.s2BandsIn,self.env.s2BandsOut) + + + if self.env.brdf == True: + print("apply brdf correction..") + s2 = s2.map(self.brdf) + + + if self.env.terrainCorrection == True: + print("apply terrain correction..") + s2 = s2.map(self.terrain) + + print("calculating medoid") + img = self.medoidMosaic(s2) + + medoid = self.medoidMosaic(s2) + medoidDown = ee.Image(self.medoidMosaicPercentiles(s2,self.env.percentiles[0])) + medoidUp = self.medoidMosaicPercentiles(s2,self.env.percentiles[1]) + stdevBands = self.addSTDdev(s2) + + #med = s2.median() + #up = s2.reduce(ee.Reducer.percentile([20])) + #down = s2.reduce(ee.Reducer.percentile([80])) + #img = med.addBands(down).addBands(up).addBands(stdevBands) + + img = img.addBands(medoidDown).addBands(medoidUp).addBands(stdevBands) + + print("rescale") + img = self.reScaleS2(img) + + print("set MetaData") + img = self.setMetaData(img) + + print("exporting composite") + self.exportMap(img,studyArea,week) + + + def addSTDdev(self,collection): + + def addSTDdevIndices(img): + """ Function to add common (and less common) spectral indices to an image. + Includes the Normalized Difference Spectral Vector from (Angiuli and Trianni, 2014) """ + img = img.addBands(img.normalizedDifference(['green','swir1']).rename(['ND_green_swir1'])); # NDSI, MNDWI + img = img.addBands(img.normalizedDifference(['nir1','red']).rename(['ND_nir_red'])); # NDVI + img = img.addBands(img.normalizedDifference(['nir1','swir2']).rename(['ND_nir_swir2'])); # NBR, MNDVI + + return img; + + + + blue_stdDev = collection.select(["blue"]).reduce(ee.Reducer.stdDev()).rename(['blue_stdDev']) + red_stdDev = collection.select(["red"]).reduce(ee.Reducer.stdDev()).rename(['red_stdDev']) + green_stdDev = collection.select(["green"]).reduce(ee.Reducer.stdDev()).rename(['green_stdDev']) + nir_stdDev = collection.select(["nir1"]).reduce(ee.Reducer.stdDev()).rename(['nir_stdDev']) + swir1_stdDev = collection.select(["swir1"]).reduce(ee.Reducer.stdDev()).rename(['swir1_stdDev']) + swir2_stdDev = collection.select(["swir2"]).reduce(ee.Reducer.stdDev()).rename(['swir2_stdDev']) + + col = collection.map(addSTDdevIndices) + + ND_green_swir1 = col.select(['ND_green_swir1']).reduce(ee.Reducer.stdDev()).rename(['ND_green_swir1_stdDev']); + ND_nir_red = col.select(['ND_nir_red']).reduce(ee.Reducer.stdDev()).rename(['ND_nir_red_stdDev']); + ND_nir_swir2 = col.select(['ND_nir_swir2']).reduce(ee.Reducer.stdDev()).rename(['ND_nir_swir2_stdDev']); + + stdevBands = ee.Image(blue_stdDev.addBands(red_stdDev).addBands(green_stdDev).addBands(nir_stdDev).addBands(swir1_stdDev).addBands(swir2_stdDev)\ + .addBands(ND_green_swir1).addBands(ND_nir_red).addBands(ND_nir_swir2)) + + return stdevBands + + def getSentinel2(self,start,end,studyArea): + + s2s = ee.ImageCollection('COPERNICUS/S2').filterDate(start,end) \ + .filterBounds(studyArea) \ + .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',self.env.metadataCloudCoverMax)) \ + .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',self.env.metadataCloudCoverMax)) + + s2s = s2s.map(self.pixelArea) + #s2s = s2s.filter(ee.Filter.gt("pixelArea",2000)) + return s2s + + + def maskShadows(self,collection,studyArea): + + def TDOM(image): + zScore = image.select(shadowSumBands).subtract(irMean).divide(irStdDev) + irSum = image.select(shadowSumBands).reduce(ee.Reducer.sum()) + TDOMMask = zScore.lt(self.env.zScoreThresh).reduce(ee.Reducer.sum()).eq(2)\ + .And(irSum.lt(self.env.shadowSumThresh)).Not() + TDOMMask = TDOMMask.focal_min(self.env.dilatePixels) + return image.addBands(TDOMMask.rename(['TDOMMask'])) + + def mask(image): + return image.updateMask(image.select(['TDOMMask'])) + + shadowSumBands = ['B8','B11'] + + allCollection = ee.ImageCollection('COPERNICUS/S2').filterBounds(studyArea).filterDate("2017-01-01","2018-01-01").map(self.QAMaskCloud) + + # Get some pixel-wise stats for the time series + irStdDev = allCollection.select(shadowSumBands).reduce(ee.Reducer.stdDev()) + irMean = allCollection.select(shadowSumBands).reduce(ee.Reducer.mean()) + + # Mask out dark dark outliers + collection_tdom = collection.map(TDOM) + + return collection_tdom.map(mask) + + def getTopo(self,img): + ''' funtion to filter for areas with terrain and areas without''' + dem = self.env.dem.unmask(0) + geom = ee.Geometry(img.get('system:footprint')).bounds() + slp_rad = ee.Terrain.slope(dem).clip(geom); + + slope = slp_rad.reduceRegion(reducer= ee.Reducer.percentile([80]), \ + geometry= geom,\ + scale= 100 ,\ + maxPixels=10000000) + return img.set('slope',slope.get('slope')) + + def scaleS2(self,img): + + divideBands = ['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12'] + bandNames = img.bandNames() + otherBands = bandNames.removeAll(divideBands) + + others = img.select(otherBands) + out = img.select(divideBands).divide(10000) + return out.addBands(others).set(img.toDictionary()).copyProperties(img,['system:time_start','system:footprint','MEAN_SOLAR_ZENITH_ANGLE','MEAN_SOLAR_AZIMUTH_ANGLE']).set("centroid",img.geometry().centroid()); + + + + def reScaleS2(self,img): + + #bandNames = img.bandNames() + #otherBands = bandNames.removeAll(self.env.divideBands) + #others = img.select(otherBands) + + # = img.select(self.env.divideBands); + # = t.multiply(10000) + + #ut = ee.Image(t.copyProperties(img).copyProperties(img,['system:time_start'])).addBands(others).int16() + out = img.multiply(10000).int16() + return out; + + + def pixelArea(self,img): + geom = ee.Geometry(img.get('system:footprint')).bounds() + + area = img.select(['B4']).gt(0).reduceRegion(reducer= ee.Reducer.sum(),\ + geometry= geom,\ + scale= 1000,\ + maxPixels=10000000) + + return img.set("pixelArea",area.get("B4")) + + def removeHaze(self,img): + ndsi = img.normalizedDifference(['green','swir1']) + ndvi = img.normalizedDifference(['nir1','red']) + ndwi = img.normalizedDifference(['green','nir1']) + + meanVis = img.expression('(blue + green + red) / 3',{ + 'blue' : img.select("blue"), + 'red' : img.select("red"), + 'green' : img.select("green")}) + + whiteness = img.expression('(abs(blue - meanVis) + abs(green - meanVis) + abs(red - meanVis)) / meanVis', + { 'blue' : img.select("blue"), + 'red' : img.select("red"), + 'green' : img.select("green"), + 'meanVis' : meanVis}) + + variabilityProbability = img.expression('1 - max(max(abs(ndvi), abs(ndsi)), whiteness)', + { 'ndvi' : ndvi, + 'ndsi' : ndsi, + 'whiteness' : whiteness + }) + + mask = variabilityProbability.lt(-0.5) + variabilityProbability = variabilityProbability.max(0.1) + + + cirrusCloudProbability = img.select("cirrus").divide(0.04) + + hazeProbability = img.expression('min(50 * max(blue - 0.5 * red - 0.06, 0), 1)', + {'blue' : img.select("blue"), + 'red' : img.select("red")}) + + + aerosolProbability = img.expression('max(100 * pow(aerosol - 0.15, 2), 0)', + { 'aerosol' : img.select("cb")}) + + hazeScore = img.expression('max(1 - (hazeProbability + variabilityProbability + cirrusCloudProbability + aerosolProbability), 0)', + {'hazeProbability' : hazeProbability, + 'variabilityProbability' :variabilityProbability, + 'cirrusCloudProbability' : cirrusCloudProbability, + 'aerosolProbability' : aerosolProbability}) + + return img.updateMask(hazeScore.neq(0)) + + + def TOAtoSR(self,img): + + #TDOMMask = img.select(['TDOMMask']) + info = self.collectionMeta[self.env.feature]['properties'] + scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)# i.e. Python uses seconds, EE uses milliseconds + solar_z = info['MEAN_SOLAR_ZENITH_ANGLE'] + + geom = ee.Geometry(info['system:footprint']).centroid() + #geom = ee.Geometry.Point([info['centroid']['coordinates'][0],info['centroid']['coordinates'][1]]) + date = ee.Date.fromYMD(scene_date.year,scene_date.month,scene_date.day) + + h2o = Atmospheric.water(geom,date).getInfo() + o3 = Atmospheric.ozone(geom,date).getInfo() + aot = Atmospheric.aerosol(geom,date).getInfo() + + SRTM = ee.Image('CGIAR/SRTM90_V4')# Shuttle Radar Topography mission covers *most* of the Earth + alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom).get('elevation').getInfo() + + if alt: + km = alt/1000 # i.e. Py6S uses units of kilometers + + else: + km = 0 + # Instantiate + s = SixS() + + # Atmospheric constituents + s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3) + s.aero_profile = AeroProfile.Continental + s.aot550 = aot + + # Earth-Sun-satellite geometry + s.geometry = Geometry.User() + s.geometry.view_z = 0 # always NADIR (I think..) + s.geometry.solar_z = solar_z # solar zenith angle + s.geometry.month = scene_date.month # month and day used for Earth-Sun distance + s.geometry.day = scene_date.day # month and day used for Earth-Sun distance + s.altitudes.set_sensor_satellite_level() + s.altitudes.set_target_custom_altitude(km) + + def spectralResponseFunction(bandname): + """ + Extract spectral response function for given band name + """ + + bandSelect = { + 'B1':PredefinedWavelengths.S2A_MSI_01, + 'B2':PredefinedWavelengths.S2A_MSI_02, + 'B3':PredefinedWavelengths.S2A_MSI_03, + 'B4':PredefinedWavelengths.S2A_MSI_04, + 'B5':PredefinedWavelengths.S2A_MSI_05, + 'B6':PredefinedWavelengths.S2A_MSI_06, + 'B7':PredefinedWavelengths.S2A_MSI_07, + 'B8':PredefinedWavelengths.S2A_MSI_08, + 'B8A':PredefinedWavelengths.S2A_MSI_09, + 'B9':PredefinedWavelengths.S2A_MSI_10, + 'B10':PredefinedWavelengths.S2A_MSI_11, + 'B11':PredefinedWavelengths.S2A_MSI_12, + 'B12':PredefinedWavelengths.S2A_MSI_13} + + return Wavelength(bandSelect[bandname]) + + def toa_to_rad(bandname): + """ + Converts top of atmosphere reflectance to at-sensor radiance + """ + + # solar exoatmospheric spectral irradiance + ESUN = info['SOLAR_IRRADIANCE_'+bandname] + solar_angle_correction = math.cos(math.radians(solar_z)) + + # Earth-Sun distance (from day of year) + doy = scene_date.timetuple().tm_yday + d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year + + # conversion factor + multiplier = ESUN*solar_angle_correction/(math.pi*d**2) + + # at-sensor radiance + rad = img.select(bandname).multiply(multiplier) + return rad + + + def surface_reflectance(bandname): + """ + Calculate surface reflectance from at-sensor radiance given waveband name + """ + + # run 6S for this waveband + s.wavelength = spectralResponseFunction(bandname) + s.run() + + # extract 6S outputs + Edir = s.outputs.direct_solar_irradiance #direct solar irradiance + Edif = s.outputs.diffuse_solar_irradiance #diffuse solar irradiance + Lp = s.outputs.atmospheric_intrinsic_radiance #path radiance + absorb = s.outputs.trans['global_gas'].upward #absorption transmissivity + scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity + tau2 = absorb*scatter #total transmissivity + + # radiance to surface reflectance + rad = toa_to_rad(bandname) + + ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif)) + + return ref + + # all wavebands + output = img.select('QA60') + for band in ['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12']: + output = output.addBands(surface_reflectance(band)) + + self.env.feature +=1 + + return output.copyProperties(img,['system:time_start','system:footprint','MEAN_SOLAR_ZENITH_ANGLE','MEAN_SOLAR_AZIMUTH_ANGLE']) + + + # Function to mask clouds using the Sentinel-2 QA band. + def QAMaskCloud(self,img): + + bandNames = img.bandNames() + #otherBands = bandNames.removeAll(self.env.divideBands) + #others = img.select(otherBands) + + + qa = img.select('QA60').int16().unmask(0); + + #img = img.select(self.env.divideBands) + + + # Bits 10 and 11 are clouds and cirrus, respectively. + cloudBitMask = int(math.pow(2, 10)); + cirrusBitMask = int(math.pow(2, 11)); + + # Both flags should be set to zero, indicating clear conditions. + mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(qa.bitwiseAnd(cirrusBitMask).eq(0)); + + img = img.updateMask(mask)#.addBands(others) + + # Return the masked and scaled data. + return img + + def sentinelCloudScore(self,img): + """ + Computes spectral indices of cloudyness and take the minimum of them. + + Each spectral index is fairly lenient because the group minimum + is a somewhat stringent comparison policy. side note -> this seems like a job for machine learning :) + + originally written by Matt Hancher for Landsat imagery + adapted to Sentinel by Chris Hewig and Ian Housman + """ + + + def rescale(img, thresholds): + """ + Linear stretch of image between two threshold values. + """ + return img.subtract(thresholds[0]).divide(thresholds[1] - thresholds[0]) + + # cloud until proven otherwise + score = ee.Image(1) + blueCirrusScore = ee.Image(0) + + # clouds are reasonably bright + blueCirrusScore = blueCirrusScore.max(rescale(img.select(['blue']), [0.1, 0.5])) + blueCirrusScore = blueCirrusScore.max(rescale(img.select(['cb']), [0.1, 0.5])) + blueCirrusScore = blueCirrusScore.max(rescale(img.select(['cirrus']), [0.1, 0.3])) + score = score.min(blueCirrusScore) + + score = score.min(rescale(img.select(['red']).add(img.select(['green'])).add(img.select('blue')), [0.2, 0.8])) + score = score.min(rescale(img.select(['nir1']).add(img.select(['swir1'])).add(img.select('swir2')), [0.3, 0.8])) + + # clouds are moist + ndsi = img.normalizedDifference(['green','swir1']) + score=score.min(rescale(ndsi, [0.8, 0.6])) + score = score.multiply(100).byte(); + score = score.clamp(0,100); + + return img.addBands(score.rename(['cloudScore'])) + + + def cloudMasking(self,collection): + + def maskClouds(img): + + cloudMask = img.select(['cloudScore']).lt(self.env.cloudScoreThresh)\ + .focal_min(self.env.dilatePixels) \ + .focal_max(self.env.contractPixels) \ + .rename(['cloudMask']) + + bandNames = img.bandNames() + otherBands = bandNames.removeAll(self.env.divideBands) + others = img.select(otherBands) + + img = img.select(self.env.divideBands).updateMask(cloudMask) + + return img.addBands(cloudMask).addBands(others); + + + # Find low cloud score pctl for each pixel to avoid comission errors + minCloudScore = collection.select(['cloudScore']).reduce(ee.Reducer.percentile([self.env.cloudScorePctl])); + + collection = collection.map(maskClouds) + + return collection + + def brdf(self,img): + + + + def _apply(image, kvol, kvol0): + blue = _correct_band(image, 'blue', kvol, kvol0, f_iso=0.0774, f_geo=0.0079, f_vol=0.0372) + green = _correct_band(image, 'green', kvol, kvol0, f_iso=0.1306, f_geo=0.0178, f_vol=0.0580) + red = _correct_band(image, 'red', kvol, kvol0, f_iso=0.1690, f_geo=0.0227, f_vol=0.0574) + re1 = _correct_band(image, 're1', kvol, kvol0, f_iso=0.2085, f_geo=0.0256, f_vol=0.0845) + re2 = _correct_band(image, 're2', kvol, kvol0, f_iso=0.2316, f_geo=0.0273, f_vol=0.1003) + re3 = _correct_band(image, 're3', kvol, kvol0, f_iso=0.2599, f_geo=0.0294, f_vol=0.1197) + nir1 = _correct_band(image, 'nir1', kvol, kvol0, f_iso=0.3093, f_geo=0.0330, f_vol=0.1535) + nir2 = _correct_band(image, 'nir2', kvol, kvol0, f_iso=0.2907, f_geo=0.0410, f_vol=0.1611) + swir1 = _correct_band(image, 'swir1', kvol, kvol0, f_iso=0.3430, f_geo=0.0453, f_vol=0.1154) + swir2 = _correct_band(image, 'swir2', kvol, kvol0, f_iso=0.2658, f_geo=0.0387, f_vol=0.0639) + return replace_bands(image, [blue, green, red,re1,re2,re3, nir1,nir2, swir1, swir2]) + + + def _correct_band(image, band_name, kvol, kvol0, f_iso, f_geo, f_vol): + """fiso + fvol * kvol + fgeo * kgeo""" + iso = ee.Image(f_iso) + geo = ee.Image(f_geo) + vol = ee.Image(f_vol) + pred = vol.multiply(kvol).add(geo.multiply(kvol)).add(iso).rename(['pred']) + pred0 = vol.multiply(kvol0).add(geo.multiply(kvol0)).add(iso).rename(['pred0']) + cfac = pred0.divide(pred).rename(['cfac']) + corr = image.select(band_name).multiply(cfac).rename([band_name]) + return corr + + + def _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)""" + + relative_azimuth = sunAz.subtract(viewAz).rename(['relAz']) + pa1 = viewZen.cos() \ + .multiply(sunZen.cos()) + pa2 = viewZen.sin() \ + .multiply(sunZen.sin()) \ + .multiply(relative_azimuth.cos()) + phase_angle1 = pa1.add(pa2) + phase_angle = phase_angle1.acos() + p1 = ee.Image(PI().divide(2)).subtract(phase_angle) + p2 = p1.multiply(phase_angle1) + p3 = p2.add(phase_angle.sin()) + p4 = sunZen.cos().add(viewZen.cos()) + p5 = ee.Image(PI().divide(4)) + + kvol = p3.divide(p4).subtract(p5).rename(['kvol']) + + viewZen0 = ee.Image(0) + pa10 = viewZen0.cos() \ + .multiply(sunZen.cos()) + pa20 = viewZen0.sin() \ + .multiply(sunZen.sin()) \ + .multiply(relative_azimuth.cos()) + phase_angle10 = pa10.add(pa20) + phase_angle0 = phase_angle10.acos() + p10 = ee.Image(PI().divide(2)).subtract(phase_angle0) + p20 = p10.multiply(phase_angle10) + p30 = p20.add(phase_angle0.sin()) + p40 = sunZen.cos().add(viewZen0.cos()) + p50 = ee.Image(PI().divide(4)) + + kvol0 = p30.divide(p40).subtract(p50).rename(['kvol0']) + + return (kvol, kvol0) + + date = img.date() + footprint = determine_footprint(img) + footprint = ee.List(img.geometry().bounds().bounds().coordinates().get(0)); + (sunAz, sunZen) = sun_angles.create(date, footprint) + (viewAz, viewZen) = view_angles.create(footprint) + (kvol, kvol0) = _kvol(sunAz, sunZen, viewAz, viewZen) + + bandNames = img.bandNames() + otherBands = bandNames.removeAll(self.env.divideBands) + others = img.select(otherBands) + + img = ee.Image(_apply(img, kvol.multiply(PI()), kvol0.multiply(PI()))) + + return img + + def terrain(self,img): + degree2radian = 0.01745; + + geom = ee.Geometry(img.get('system:footprint')).bounds() #.buffer(10000) + dem = self.env.dem + + bandNames = img.bandNames() + otherBands = bandNames.removeAll(self.env.divideBands) + others = img.select(otherBands) + + bandList = ['blue','green','red','re1','re2','re3','nir1','nir2','cb','cirrus','swir1','swir2','waterVapor'] + + def topoCorr_IC(img): + + # Extract image metadata about solar position + SZ_rad = ee.Image.constant(ee.Number(img.get('MEAN_SOLAR_ZENITH_ANGLE'))).multiply(degree2radian).clip(geom); + SA_rad = ee.Image.constant(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE'))).multiply(degree2radian).clip(geom); + + + # Creat terrain layers + slp = ee.Terrain.slope(dem).clip(geom); + slp_rad = ee.Terrain.slope(dem).multiply(degree2radian).clip(geom); + asp_rad = ee.Terrain.aspect(dem).multiply(degree2radian).clip(geom); + + + # Calculate the Illumination Condition (IC) + # slope part of the illumination condition + cosZ = SZ_rad.cos(); + cosS = slp_rad.cos(); + slope_illumination = cosS.expression("cosZ * cosS", \ + {'cosZ': cosZ, 'cosS': cosS.select('slope')}); + + # aspect part of the illumination condition + sinZ = SZ_rad.sin(); + sinS = slp_rad.sin(); + cosAziDiff = (SA_rad.subtract(asp_rad)).cos(); + aspect_illumination = sinZ.expression("sinZ * sinS * cosAziDiff", \ + {'sinZ': sinZ, \ + 'sinS': sinS, \ + 'cosAziDiff': cosAziDiff}); + + # full illumination condition (IC) + ic = slope_illumination.add(aspect_illumination); + + # Add IC to original image + img_plus_ic = ee.Image(img.addBands(ic.rename(['IC'])).addBands(cosZ.rename(['cosZ'])).addBands(cosS.rename(['cosS'])).addBands(slp.rename(['slope']))); + + return ee.Image(img_plus_ic); + + + def topoCorr_SCSc(img): + img_plus_ic = img; + mask1 = img_plus_ic.select('nir1').gt(-0.1); + mask2 = img_plus_ic.select('slope').gte(5) \ + .And(img_plus_ic.select('IC').gte(0)) \ + .And(img_plus_ic.select('nir1').gt(-0.1)); + + img_plus_ic_mask2 = ee.Image(img_plus_ic.updateMask(mask2)); + + + + def apply_SCSccorr(band): + method = 'SCSc'; + + out = ee.Image(1).addBands(img_plus_ic_mask2.select('IC', band)).reduceRegion(reducer= ee.Reducer.linearRegression(2,1), \ + geometry= ee.Geometry(img.geometry()), \ + scale= 50, \ + bestEffort =True, + maxPixels=1e10) + + + fit = out.combine({"coefficients": ee.Array([[1],[1]])}, False); + + #Get the coefficients as a nested list, + #cast it to an array, and get just the selected column + out_a = (ee.Array(fit.get('coefficients')).get([0,0])); + out_b = (ee.Array(fit.get('coefficients')).get([1,0])); + out_c = out_a.divide(out_b) + + + # apply the SCSc correction + SCSc_output = img_plus_ic_mask2.expression("((image * (cosB * cosZ + cvalue)) / (ic + cvalue))", { + 'image': img_plus_ic_mask2.select([band]), + 'ic': img_plus_ic_mask2.select('IC'), + 'cosB': img_plus_ic_mask2.select('cosS'), + 'cosZ': img_plus_ic_mask2.select('cosZ'), + 'cvalue': out_c }); + + return ee.Image(SCSc_output); + + img_SCSccorr = ee.Image([apply_SCSccorr(band) for band in bandList]).addBands(img_plus_ic.select('IC')); + + bandList_IC = ee.List([bandList, 'IC']).flatten(); + + img_SCSccorr = img_SCSccorr.unmask(img_plus_ic.select(bandList_IC)).select(bandList); + + return img_SCSccorr.unmask(img_plus_ic.select(bandList)) + + img = topoCorr_IC(img) + img = topoCorr_SCSc(img).addBands(others ) + + return img + + + def medoidMosaicPercentiles(self,inCollection,p): + ' calculate the medoid of a percentile' + + + inCollection = inCollection.select(self.env.medoidBands) + + p1 = p + p2 = 100 -p + + #med1 = self.medoidPercentiles(inCollection,p1).select(["green","re1","re2","re3","nir1","nir2"]) + #med2 = self.medoidPercentiles(inCollection,p2).select(["blue","red","swir1","swir2","cb","cirrus","waterVapor"]) + + med1 = inCollection.select(["green","re1","re2","re3","nir1","nir2"]).reduce(ee.Reducer.percentile([p1])) #.select(["green","re1","re2","re3","nir1","nir2"]) + med2 = inCollection.select(["blue","red","swir1","swir2","cb","cirrus","waterVapor"]).reduce(ee.Reducer.percentile([p2])) #.select(["blue","red","swir1","swir2","cb","cirrus","waterVapor"]) + + #medoidP = self.renameBands(ee.Image(med1).addBands(med2),str("p")+str(p)) + medoidP = ee.Image(med1).addBands(med2) + return medoidP + + def medoidPercentiles(self,inCollection,p): + + # Find band names in first image + bandNumbers = ee.List.sequence(1,self.env.medoidBands.length()); + + # Find the median + percentile = inCollection.select(self.env.medoidBands).reduce(ee.Reducer.percentile([p])); + + def subtractPercentile(img): + diff = ee.Image(img).subtract(percentile).pow(ee.Image.constant(2)); + return diff.reduce('sum').addBands(img); + + percentile = inCollection.map(subtractPercentile) + + percentile = ee.ImageCollection(percentile).reduce(ee.Reducer.min(self.env.medoidBands.length().add(1))).select(bandNumbers,self.env.medoidBands); + + return percentile; + + def medoidMosaic(self,collection): + """ medoid composite with equal weight among indices """ + + bandNames = ee.Image(collection.first()).bandNames() + otherBands = bandNames.removeAll(self.env.divideBands) + + others = collection.select(otherBands).reduce(ee.Reducer.mean()).rename(otherBands); + + collection = collection.select(self.env.divideBands) + + bandNumbers = ee.List.sequence(1,self.env.divideBands.length()); + + median = ee.ImageCollection(collection).median() + + def subtractmedian(img): + diff = ee.Image(img).subtract(median).pow(ee.Image.constant(2)); + return diff.reduce('sum').addBands(img); + + medoid = collection.map(subtractmedian) + + medoid = ee.ImageCollection(medoid).reduce(ee.Reducer.min(self.env.divideBands.length().add(1))).select(bandNumbers,self.env.divideBands); + + return medoid.addBands(others); + + + def medianMosaic(self,collection): + + """ median composite """ + median = collection.select(medianIncludeBands).median(); + othersBands = bandNames.removeAll(medianIncludeBands); + others = collection.select(otherBands).mean(); + + return median.addBands(others) + + + def setMetaData(self,img): + """ add metadata to image """ + + img = ee.Image(img).set({'system:time_start':ee.Date(self.env.startDate).millis(), \ + 'startDOY':str(self.env.startDoy), \ + 'endDOY':str(self.env.endDoy), \ + 'useCloudScore':str(self.env.cloudMask), \ + 'useTDOM':str(self.env.shadowMask), \ + 'useQAmask':str(self.env.QAcloudMask), \ + 'useCloudProject':str(self.env.cloudMask), \ + 'terrain':str(self.env.terrainCorrection), \ + 'surfaceReflectance':str(self.env.calcSR), \ + 'cloudScoreThresh':str(self.env.cloudScoreThresh), \ + 'cloudScorePctl':str(self.env.cloudScorePctl), \ + 'zScoreThresh':self.env.zScoreThresh, \ + 'shadowSumThresh':self.env.shadowSumThresh, \ + 'contractPixels':str(self.env.contractPixels), \ + 'crs':self.env.epsg, \ + 'cloudFilter':str(self.env.metadataCloudCoverMax),\ + 'dilatePixels':self.env.dilatePixels}) + + return img + + def exportMap(self,img,studyArea,week): + + geom = studyArea.bounds().getInfo(); + + task_ordered= ee.batch.Export.image.toAsset(image=img.clip(studyArea.buffer(10000)), + description = self.env.name, #+ str(week), + assetId= self.env.assetId + self.env.name + str(week).zfill(3), + region=geom['coordinates'], + maxPixels=1e13, + crs=self.env.epsg, + scale=self.env.exportScale) + + task_ordered.start() + + + def renameBands(self,image,prefix): + 'rename bands with prefix' + + bandnames = image.bandNames(); + + def mapBands(band): + band = ee.String(prefix).cat('_').cat(band); + return band; + + bandnames = bandnames.map(mapBands) + + image = image.rename(bandnames); + + return image; + + +if __name__ == "__main__": + + ee.Initialize() + + startDate = year.advance(startDay,"day") + endDate = year.advance(endDay,"day") + + studyArea = na.merge(th).geometry().bounds() + studyArea = ee.FeatureCollection("users/KarisTenneson/MyanmarRubberPlantations/StudyRegion").geometry().bounds() + + functions().main(studyArea,startDate,endDate,startDay,endDay,startWeek) + diff --git a/sun_angles.py b/sun_angles.py new file mode 100644 index 0000000..eefa4a8 --- /dev/null +++ b/sun_angles.py @@ -0,0 +1,66 @@ + +from utils import * + +def create(date, footprint): + jdp = date.getFraction('year') + seconds_in_hour = 3600 + hourGMT = ee.Number(date.getRelative('second', 'day')) \ + .divide(seconds_in_hour) + + latRad = degToRad(ee.Image.pixelLonLat().select('latitude')) + longDeg = ee.Image.pixelLonLat().select('longitude') + + # Julian day proportion in radians + jdpr = jdp.multiply(PI()).multiply(2) + + a = ee.List([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]) + meanSolarTime = longDeg.divide(15.0).add(ee.Number(hourGMT)) + 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())) + + localSolarDiff2 = localSolarDiff1.multiply(12 * 60) + + localSolarDiff = localSolarDiff2.divide(PI()) + trueSolarTime = meanSolarTime \ + .add(localSolarDiff.divide(60)) \ + .subtract(12.0) + + # Hour as an angle + ah = trueSolarTime.multiply(degToRad(ee.Number(MAX_SATELLITE_ZENITH * 2))) + b = ee.List([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]) + 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())) + cosSunZen = latRad.sin().multiply(delta.sin()) \ + .add(latRad.cos().multiply(ah.cos()).multiply(delta.cos())) + sunZen = cosSunZen.acos() + + # sun azimuth from south, turning west + sinSunAzSW = ah.sin().multiply(delta.cos()).divide(sunZen.sin()) + sinSunAzSW = sinSunAzSW.clamp(-1.0, 1.0) + + cosSunAzSW = (latRad.cos().multiply(-1).multiply(delta.sin()) + .add(latRad.sin().multiply(delta.cos()).multiply(ah.cos()))) \ + .divide(sunZen.sin()) + 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) + + sunAz = sunAzSW.add(PI()) + # Keep within [0, 2pi] range + sunAz = where(sunAz.gt(PI().multiply(2)), sunAz.subtract(PI().multiply(2)), sunAz) + + 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) diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..b4b2b36 --- /dev/null +++ b/utils.py @@ -0,0 +1,83 @@ +# gee helpers + +import ee +import math + +UPPER_LEFT = 0 +LOWER_LEFT = 1 +LOWER_RIGHT = 2 +UPPER_RIGHT = 3 +PI = lambda: ee.Number(math.pi) +MAX_SATELLITE_ZENITH = 7.5 + +def line_from_coords(coordinates, fromIndex, toIndex): + return ee.Geometry.LineString(ee.List([ + coordinates.get(fromIndex), + coordinates.get(toIndex)])) + + +def line(start, end): + return ee.Geometry.LineString(ee.List([start, end])) + + +def degToRad(deg): + return deg.multiply(PI().divide(180)) + + +def value(list, index): + return ee.Number(list.get(index)) + + +def radToDeg(rad): + return rad.multiply(180).divide(PI()) + + +def where(condition, trueValue, falseValue): + trueMasked = trueValue.mask(condition) + falseMasked = falseValue.mask(invertMask(condition)) + return trueMasked.unmask(falseMasked) + + +def invertMask(mask): + return mask.multiply(-1).add(1) + + +def x(point): + return ee.Number(ee.List(point).get(0)) + + +def y(point): + return ee.Number(ee.List(point).get(1)) + + +def determine_footprint(image): + footprint = ee.Geometry(image.get('system:footprint')) + bounds = ee.List(footprint.bounds().coordinates().get(0)) + coords = footprint.coordinates() + + xs = coords.map(lambda item: x(item)) + ys = coords.map(lambda item: y(item)) + + def findCorner(targetValue, values): + diff = values.map(lambda value: ee.Number(value).subtract(targetValue).abs()) + minValue = diff.reduce(ee.Reducer.min()) + idx = diff.indexOf(minValue) + return coords.get(idx) + + lowerLeft = findCorner(x(bounds.get(0)), xs) + lowerRight = findCorner(y(bounds.get(1)), ys) + upperRight = findCorner(x(bounds.get(2)), xs) + upperLeft = findCorner(y(bounds.get(3)), ys) + + return ee.List([upperLeft, lowerLeft, lowerRight, upperRight, upperLeft]) + + +def replace_bands(image, bands): + result = image + for band in bands: + result = result.addBands(band, overwrite=True) + return result + + + + diff --git a/view_angles.py b/view_angles.py new file mode 100644 index 0000000..29a3285 --- /dev/null +++ b/view_angles.py @@ -0,0 +1,30 @@ +import ee + +from utils import * + +MAX_DISTANCE = 1000000 + +def create(footprint): + return (azimuth(footprint), zenith(footprint)) + + +def azimuth(footprint): + upperCenter = line_from_coords(footprint, UPPER_LEFT, UPPER_RIGHT).centroid().coordinates() + lowerCenter = line_from_coords(footprint, LOWER_LEFT, LOWER_RIGHT).centroid().coordinates() + slope = ((y(lowerCenter)).subtract(y(upperCenter))).divide((x(lowerCenter)).subtract(x(upperCenter))) + slopePerp = ee.Number(-1).divide(slope) + azimuthLeft = ee.Image(PI().divide(2).subtract((slopePerp).atan())) + return azimuthLeft.rename(['viewAz']) + + +def zenith(footprint): + leftLine = line_from_coords(footprint, UPPER_LEFT, LOWER_LEFT) + rightLine = line_from_coords(footprint, UPPER_RIGHT, LOWER_RIGHT) + leftDistance = ee.FeatureCollection(leftLine).distance(MAX_DISTANCE) + rightDistance = ee.FeatureCollection(rightLine).distance(MAX_DISTANCE) + 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 degToRad(viewZenith)