jupyter_notebooks/atmcorr_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Atmospheric Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In remote sensing we often need to calculate surface reflectance (ρ) from radiance (L) measure by a given sensor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ρ = π(L -L<sub>p</sub>) / τ(E<sub>dir</sub> + E<sub>dif</sub>)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "\n",
    "* ρ = surface reflectance\n",
    "* L = at-sensor radiance\n",
    "* L<sub>p</sub> = path radiance\n",
    "* τ = transmissivity (from surface to satellite)\n",
    "* E<sub>dir</sub> = direct solar irradiance \n",
    "* E<sub>dif</sub> = diffuse solar irradiance\n",
    "* π = 3.1415\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say a satellite sensor measures a radiance of 120"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = 120"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 4 unknowns remaining and they on i) atmospheric conditions and ii) Earth-Sun-Satellite geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Atmospheric Conditions** <br />\n",
    "let's say we have measured values of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "H2O = 1    # water vapour (g cm-2)\n",
    "O3 = 0.4   # ozone (atm-cm)\n",
    "AOT = 0.3  # aerosol optical thickness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Earth-Sun-Satellite Geometry** <br />\n",
    "and some additional measurements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alt = 0      # target altitude (km)\n",
    "solar_z = 20 # solar zenith angle (degrees)\n",
    "view_z = 0   # view zentith angle (degrees)\n",
    "doy = 4      # Earth-Sun distance (Astronomical Units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Potential data sources**\n",
    "\n",
    "* Water vapour: [NCEP/NCAR](http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%281996%29077%3C0437%3ATNYRP%3E2.0.CO%3B2)\n",
    "* Ozone: [TOMS/OMI](http://ozoneaq.gsfc.nasa.gov/missions). \n",
    "* Aerosol optical thickness: [MODIS Aerosol Product](http://modis-atmos.gsfc.nasa.gov/MOD04_L2/index.html) or in-scene techniques\n",
    "* Geometry and day-of-year: satellite image metadata\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6S emulator**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# loading dependencies\n",
    "import os\n",
    "import sys\n",
    "sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))\n",
    "from interpolated_LUTs import Interpolated_LUTs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# instantiate interpolated look up table class\n",
    "iLUTs = Interpolated_LUTs('COPERNICUS/S2') # i.e. Sentinel 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# download look-up tables\n",
    "iLUTs.download_LUTs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# interpolate look-up tables\n",
    "iLUTs.interpolate_LUTs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iLUTs_all_wavebands = iLUTs.get()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# for example let's look at band 1\n",
    "iLUT_B1 = iLUTs_all_wavebands['B1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get atmospheric correction coefficients (a, b) for *perihelion*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a, b = iLUT_B1(solar_z,H2O,O3,AOT,alt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and corrected for Earth's [elliptical orbit](https://github.com/samsammurphy/6S_LUT/wiki/Elliptical-Orbit-Correction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "elliptical_orbit_correction = 0.03275104*math.cos(doy/59.66638337) + 0.96804905\n",
    "a *= elliptical_orbit_correction\n",
    "b *= elliptical_orbit_correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to get the correction coefficients we need to calculate surface reflectance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ρ = (L-a)/b\n",
    "\n",
    "print('Surface Reflectance = {:.3f}'.format(ρ))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}