

In [1]:

    
import ee
ee.Initialize()



In [2]:

    
import os
import sys
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









    



---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-2-b2793274f3f3> in <module>()
      3 sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))
      4 from sixs_emulator_ee_sentinel2_batch import SixS_emulator
----> 5 from atmcorr_input import Atmcorr_input
      6 from atmospheric_correction import atmospheric_correction
      7 from radiance import radiance_from_TOA

ModuleNotFoundError: No module named 'atmcorr_input'



In [3]:

    
# a place and a mission
geom = ee.Geometry.Point(-157.816222, 21.297481)
mission = 'COPERNICUS/S2'



In [4]:

    
# 6S emulator
se = SixS_emulator(mission)



In [5]:

    
# image collection
ic = ee.ImageCollection(mission)\
    .filterBounds(geom)\
    .filterDate('2017-01-01','2017-02-01')\
    .filter(ee.Filter.lt('MEAN_SOLAR_ZENITH_ANGLE',75))



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()









    



Downloading look up table (LUT) zip file..
Extracting zip file..
Done: LUT files available locally



In [8]:

    
# and run the interpolation
iLUTs.interpolate_LUTs()









    



running n-dimensional interpolation may take a few minutes...
iLUT file already exists (skipping interpolation): S2A_MSI_01.lut
iLUT file already exists (skipping interpolation): S2A_MSI_02.lut
iLUT file already exists (skipping interpolation): S2A_MSI_03.lut
iLUT file already exists (skipping interpolation): S2A_MSI_04.lut
iLUT file already exists (skipping interpolation): S2A_MSI_05.lut
iLUT file already exists (skipping interpolation): S2A_MSI_06.lut
iLUT file already exists (skipping interpolation): S2A_MSI_07.lut
iLUT file already exists (skipping interpolation): S2A_MSI_08.lut
iLUT file already exists (skipping interpolation): S2A_MSI_09.lut
iLUT file already exists (skipping interpolation): S2A_MSI_10.lut
iLUT file already exists (skipping interpolation): S2A_MSI_11.lut
iLUT file already exists (skipping interpolation): S2A_MSI_12.lut
iLUT file already exists (skipping interpolation): S2A_MSI_13.lut

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
ic_atmospherically_corrected = []

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)
    ic_atmospherically_corrected.append(SR)



In [12]:

    
# Earth Engine image collection
ic_atmospherically_corrected = ee.ImageCollection(ic_atmospherically_corrected)



In [13]:

    
# test
SR = ee.Image(ic_atmospherically_corrected.first())
print(SR.reduceRegion(ee.Reducer.mean(),geom).getInfo())









    



{'B1': 0.15776946933165464, 'B10': 1.0800332271125255, 'B11': 0.26621009885365277, 'B12': 0.1998808958940913, 'B2': 0.23538158062086803, 'B3': 0.25703131294111686, 'B4': 0.28050916213457583, 'B5': 0.28684050502120506, 'B6': 0.32420061980635617, 'B7': 0.32387646618908117, 'B8': 0.3293928233008832, 'B8A': 0.32055806565533307, 'B9': 0.44896412979808153}



In [14]:

    
# original image (top of atmosphere)
TOA = ee.Image(ic.first()).divide(10000)

# surface reflectance (bottom of atmosphere)
SR = ee.Image(ic_atmospherically_corrected.first())

from IPython.display import display, Image

region = geom.buffer(5000).bounds().getInfo()['coordinates']

before = Image(url=TOA.select(['B4','B3','B2']).getThumbUrl({
                'region':region,
                'min':0,
                'max':0.25,
                'gamma':1.5
                }))

after = Image(url=SR.select(['B4','B3','B2']).getThumbUrl({
                'region':region,
                'min':0,
                'max':0.25,
                'gamma':1.5
                }))

display(before, after)