In [1]:
import ee
ee.Initialize()
In [2]:
import os
import sys, math
sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))
from sixs_emulator_ee_sentinel2_batch import SixS_emulator
from atmcorr_input import Atmcorr_input
from atmospheric_correction import atmospheric_correction
from radiance import radiance_from_TOA
from interpolated_LUTs import Interpolated_LUTs
In [3]:
# a place and a mission
geom = ee.Geometry.Point(103.16, 21.32)
mission = 'COPERNICUS/S2'
In [4]:
# 6S emulator
se = SixS_emulator(mission)
In [5]:
# image collection
ic = ee.ImageCollection(mission)\
.filterBounds(geom)\
.filterDate('2017-02-01','2017-04-01')\
.sort('CLOUD_COVERAGE_ASSESSMENT')
In [6]:
# iLUTs instance
iLUTs = Interpolated_LUTs('COPERNICUS/S2') # i.e. Sentinel2
In [7]:
# if this is first time you might have to download the look up tables
iLUTs.download_LUTs()
In [8]:
# and run the interpolation
iLUTs.interpolate_LUTs()
If you are running this notebook in a docker container then you can save these interpolated look-up tables (and your Earth Engine authentication) for later using a docker commit. This will save them to memory so that you only have to do the set-up once.
In [9]:
# otherwise, you can just load iLUTs from file
se.iLUTs = iLUTs.get()
In [10]:
# extract atmcorr inputs as feature collection
Atmcorr_input.geom = geom # specify target location (would use image centroid otherwise)
atmcorr_inputs = ic.map(Atmcorr_input.extractor).getInfo()
features = atmcorr_inputs['features']
In [11]:
# atmospherically correct image collection
imgCollection = []
for feature in features:
# at-sensor radiance
toa = ee.Image(mission+'/'+feature['properties']['imgID'])
rad = radiance_from_TOA(toa, feature)
# 6S emulator
cc = se.run(feature['properties']['atmcorr_inputs'])
# atmospheric correction
SR = atmospheric_correction(rad, cc)
imgCollection.append(SR)
In [12]:
# Earth Engine image collection
imgCollection = ee.ImageCollection(imgCollection)
In [13]:
# test
SR = ee.Image(imgCollection.first())
#print(SR.reduceRegion(ee.Reducer.mean(),geom).getInfo())
In [14]:
# this is a nice idea!
def projectShadows(cloudMask, sunAzimuth, offset):
# ^
# |
# shouldn't this offet depend on the solar elevation angle?
azimuth = ee.Number(sunAzimuth).multiply(math.pi).divide(180.0).add(ee.Number(0.5).multiply(math.pi))
x = azimuth.cos().multiply(10.0).round()
y = azimuth.sin().multiply(10.0).round()
shadow = cloudMask.changeProj(cloudMask.projection(), cloudMask.projection().translate(x.multiply(ee.Number(offset)), y.multiply(ee.Number(offset))))
return shadow
In [15]:
def cloudMaskS2 (image, SunAzimuth):
# is the multiplication of band 9 SQUARED on purpose?
# | |
# V V
# cld = image.select('B2').multiply(image.select('B3')).multiply(image.select('B9')).multiply(image.select('B9')).divide(1e11)
# cloud = cld.lt(210).focal_min(10)
# shadow = projectShadows(cloud, SunAzimuth, 6.5)
# shouldn't use mask() to set a mask
# |
# V
#cldsh = image.mask().reduce('min').And(cloud).And(sh)
#return image.mask(cldsh)
# let's use this test cloud mask for now
cloud = image.select('B2').gt(0.5)
shadow = projectShadows(cloud, SunAzimuth, 6.5)
cloudAndShadow = cloud.Or(shadow)
return image.updateMask(cloudAndShadow.eq(0)).unmask()
# ^
# |
# this sets masked pixels to zero so that they show up in this display
In [16]:
# original image (top of atmosphere)
# TOA = ee.Image(ic.sort('CLOUD_COVERAGE_ASSESSMENT', False).map(cloudMaskS2).mosaic()).divide(10000)
# ^ ^
# | |
#----------------- did you want to display the whole image collection in one go? ------------------------
# 1) same images
TOA = ee.Image(ic.first()).divide(10000)
BOA = ee.Image(SR)
# when you divide the TOA by 10,000 earth engine will not preserve all the original metadata of the input image
# (because it will often be really inefficent) this means you have to explicitly .set() it again if you need it
# or save it into a new variable, e.g.
SunAzimuth = ee.Image(ic.first()).get('MEAN_SOLAR_AZIMUTH_ANGLE')
# 2) same cloud mask
TOA_masked = cloudMaskS2(TOA, SunAzimuth)
BOA_masked = cloudMaskS2(BOA, SunAzimuth)
In [17]:
from IPython.display import display, Image
region = geom.buffer(10000).bounds().getInfo()['coordinates']
before = Image(url=TOA_masked.select(['B4','B3','B2']).getThumbUrl({
'region':region,
'min':0,
'max':0.25,
'gamma':1.5
}))
after = Image(url=BOA_masked.select(['B4','B3','B2']).getThumbUrl({
'region':region,
'min':0,
'max':0.25,
'gamma':1.5
}))
display(before, after)
In [ ]: