qq2898 / 6s_emulator Goto Github PK

The 6S emulator is an open-source atmospheric correction tool. It is based on the 6S radiative transfer model but it runs 100x faster with minimal additional error (i.e. < 0.5 %).

This speed increase is acheived by building interpolated look-up tables. This trades set-up time for execution time. The look-up tables take a long time (i.e. hours) to build, here are some prebuilts for: Sentinel 2, Landsat 8, Landsat 7 & Landsat 4 and 5. You only need to build (or download) a look-up table once.

Interpolated look-up tables are the core of the 6S emulator. Essentially, they are used to calculate atmospheric correction coefficients (a, b) which convert at-sensor radiance (L) to surface reflectance (ρ) as follows:

ρ = (L - a) / b

The installation instructions (below) are for building look-up tables. To use a pre-existing look-up table, all that is required are python3.x, numpy and scipy.

We interact with 6S through an excellent Python wrapper called Py6S and share the same dependencies.

The recommended installation method is to use the conda package and environment manager.

$ conda create -n py6s-env -c conda-forge py6s

This will create a new environment that needs to be activated.

You could permanently add the conda-forge channel if you prefer to avoid (de)activating environments.

$ conda config --add channels conda-forge

$ conda install py6s

You could optionally run the following docker container instead, which has all dependencies pre-installed

$ docker run -it samsammurphy/ee-python3-jupyter-atmcorr:v1.0

then clone this repository into the container

# git clone https://github.com/samsammurphy/6S_emulator

Here are a list of all dependencies for manual installation.

See the jupyter notebook for a quick start example of atmospheric correction.

It is much more bandwidth efficient to send and receive look-up tables, and then interpolate them locally, which is why building and interpolating are handled by separate modules. To see a more complete list of examples of how to build a look-up table (for any satellite mission) see this wiki. Here is a short example.

$ python3 LUT_build.py --wavelength 0.42

which will build a look-up table for a wavelength of 0.42 microns, it can be interpolated as follows

$ python3 LUT_interpolate.py path/to/LUT_directory

where the 'path/to/LUT_directory' is the full path to the look-up table files ('.lut').

An interpolated look-up tables is a pickle file of a scipy linear n-dimensional interpolator. It can be loaded like this:

import pickle

fpath = 'path/to/interpolated_lookup_table.ilut'

with open(fpath,"rb") as ilut_file:
    iLUT = pickle.load(ilut_file)

An interpolated look-up table requires the following input variables (in order) to provide atmospheric correction coefficients:

1. solar zentith [degrees] (0 - 75)
2. water vapour [g/m2] (0 - 8.5)
3. ozone [cm-atm] (0 - 0.8)
4. aerosol optical thickness [unitless] (0 - 3)
5. surface altitude [km] (0 - 7.75)

In code it might look something like this

a, b = iLUT(solar_z, h2o, o3, aot, km)

where a and b are the atmospheric correction coefficients at perihelion. The look-up tables are built at perihelion (i.e. January 4th) to save space because Earth's elliptical orbit can be corrected as follows:

import math

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

Surface reflectance can then be calculated from at-sensor radiance:

surface_reflectance = (L - a) / b