# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Atmospheric Correction

# In remote sensing we often need to calculate surface reflectance (ρ) from radiance (L) measure by a given sensor

# ### ρ = π(L -L<sub>p</sub>) / τ(E<sub>dir</sub> + E<sub>dif</sub>)

# where
#
# * ρ = surface reflectance
# * L = at-sensor radiance
# * L<sub>p</sub> = path radiance
# * τ = transmissivity (from surface to satellite)
# * E<sub>dir</sub> = direct solar irradiance 
# * E<sub>dif</sub> = diffuse solar irradiance
# * π = 3.1415
#
#

# Let's say a satellite sensor measures a radiance of 120

L = 120

# There are 4 unknowns remaining and they on i) atmospheric conditions and ii) Earth-Sun-Satellite geometry

# **Atmospheric Conditions** <br />
# let's say we have measured values of:

H2O = 1    # water vapour (g cm-2)
O3 = 0.4   # ozone (atm-cm)
AOT = 0.3  # aerosol optical thickness

# **Earth-Sun-Satellite Geometry** <br />
# and some additional measurements:

alt = 0      # target altitude (km)
solar_z = 20 # solar zenith angle (degrees)
view_z = 0   # view zentith angle (degrees)
doy = 4      # Earth-Sun distance (Astronomical Units)

# **Potential data sources**
#
# * Water vapour: [NCEP/NCAR](http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%281996%29077%3C0437%3ATNYRP%3E2.0.CO%3B2)
# * Ozone: [TOMS/OMI](http://ozoneaq.gsfc.nasa.gov/missions). 
# * Aerosol optical thickness: [MODIS Aerosol Product](http://modis-atmos.gsfc.nasa.gov/MOD04_L2/index.html) or in-scene techniques
# * Geometry and day-of-year: satellite image metadata
#

# **6S emulator**

# loading dependencies
import os
import sys
sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))
from interpolated_LUTs import Interpolated_LUTs

# The 6S emulator is **100x** faster than the radiative transfer code. This speed increased is acheived by using interpolated look-up tables (iLUTs). This trades set-up time for execution time.

# instantiate interpolated look up table class
iLUTs = Interpolated_LUTs('COPERNICUS/S2') # i.e. Sentinel 2

# download look-up tables
iLUTs.download_LUTs()

# interpolate look-up tables
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](https://github.com/samsammurphy/6S_emulator/wiki/docker-commits). This will save them to memory so that you only have to do the set-up once.

iLUTs_all_wavebands = iLUTs.get()

# for example let's look at band 1
iLUT_B1 = iLUTs_all_wavebands['B1']

# We can get atmospheric correction coefficients (a, b) for *perihelion*:

a, b = iLUT_B1(solar_z,H2O,O3,AOT,alt)

# and corrected for Earth's [elliptical orbit](https://github.com/samsammurphy/6S_LUT/wiki/Elliptical-Orbit-Correction):

# +
import math

elliptical_orbit_correction = 0.03275104*math.cos(doy/59.66638337) + 0.96804905
a *= elliptical_orbit_correction
b *= elliptical_orbit_correction
# -

# to get the correction coefficients we need to calculate surface reflectance.

# +
ρ = (L-a)/b

print('Surface Reflectance = {:.3f}'.format(ρ))

StarcoderdataJupyterScriptsDedupFiltered