Py6S_atmCorrectionFINAL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#######################################################################\n",
    "############### ATMOSPHERIC CORRECTION USING PY6S #####################\n",
    "#######################################################################\n",
    "## Daniel Borini Alves; Dhemerson Conciani; and Swanni Alvarado\n",
    "\n",
    "## --  adapted  based   on  Minh Nguyen  script   <https://github.com/ndminhhus/geeguide/blob/master/02.Atm-correction.md>  and \n",
    "## the original script developed by Sam Murphy <https://github.com/samsammurphy/gee-atmcorr-S2>, using PY6s python API v1.8.0,\n",
    "## <https://py6s.readthedocs.io/en/latest/> (Robin Wilson). See all references in the end of the script.\n",
    "\n",
    "## Inputs: Top of Atmosphere (TOA) reflectance images (Sentinel 2 and Landsat collections)\n",
    "## Output: Bottom of Atmosphere (BOA) reflectance of each collection, exported in Google Drive\n",
    "\n",
    "## ACTIVATE PACKAGES\n",
    "from Py6S import * \n",
    "import datetime\n",
    "import math\n",
    "import os\n",
    "import sys\n",
    "sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))\n",
    "from atmospheric import Atmospheric # custom-defined class by Sam Murphy -- check link above\n",
    "import ee\n",
    "\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available images in the period:\n",
      "S2A images =  2\n",
      "S2B images =  3\n",
      "L8 images =  1\n",
      "L7 images =  2\n"
     ]
    }
   ],
   "source": [
    "## FILTER IMAGE COLLECTIONS BY SENSORS/PLATFORMS\n",
    "\n",
    "# Define your time interval, centroid of study area and area of interest\n",
    "startDate = ee.Date('2019-09-01')\n",
    "endDate   = ee.Date('2019-10-01') #long periods can spend a lot of processing time...\n",
    "geom = ee.Geometry.Point(-61.5775286643,-8.4852883834)\n",
    "#region = geom.buffer(20000).bounds().getInfo()['coordinates']\n",
    "region = ee.FeatureCollection(\"users/danielborini/parcelas_PNCA/Recorte_ZonasAeB_Sent2\")\n",
    "\n",
    "print ('Available images in the period:')\n",
    "# Sentinel 2A collection images\n",
    "S2A_col = ee.ImageCollection('COPERNICUS/S2').filterBounds(geom).filterDate(startDate,endDate).filterMetadata('SPACECRAFT_NAME', 'equals', 'Sentinel-2A')\n",
    "col_length_S2A = S2A_col.size().getInfo()\n",
    "print ('S2A images = ', col_length_S2A)\n",
    "\n",
    "# Sentinel 2B collection images\n",
    "S2B_col = ee.ImageCollection('COPERNICUS/S2').filterBounds(geom).filterDate(startDate,endDate).filterMetadata('SPACECRAFT_NAME', 'equals', 'Sentinel-2B')\n",
    "col_length_S2B = S2B_col.size().getInfo()\n",
    "print ('S2B images = ', col_length_S2B)\n",
    "\n",
    "# Landsat OLI collection images\n",
    "L8_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterBounds(geom).filterDate(startDate,endDate)\n",
    "col_length_L8 = L8_col.size().getInfo()\n",
    "print ('L8 images = ', col_length_L8)\n",
    "\n",
    "# Landsat ETM+ collection images\n",
    "L7_col = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA').filterBounds(geom).filterDate(startDate,endDate)\n",
    "col_length_L7 = L7_col.size().getInfo()\n",
    "print ('L7 images = ', col_length_L7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sentinel 2A atmopheric correction\n",
      "exporting 20190903_S2A--->done\n",
      "exporting 20190913_S2A--->done\n"
     ]
    }
   ],
   "source": [
    "### SENTINEL 2A COLLECTION ATMOSPHERIC CORRECTION\n",
    "\n",
    "print('Sentinel 2A atmopheric correction')\n",
    "\n",
    "for i in range (0,col_length_S2A):\n",
    "    # OPEN IMAGE AND EXTRACT METADATA ATMOSPHERIC PARAMETERS\n",
    "    \n",
    "    S2_list = S2A_col.toList(col_length_S2A)\n",
    "    img = ee.Image(S2_list.get(i))\n",
    "    \n",
    "    info = img.getInfo()['properties']\n",
    "    scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)\n",
    "    solar_z = img.getInfo()['properties']['MEAN_SOLAR_ZENITH_ANGLE']\n",
    "    h2o = Atmospheric.water(geom,img.date()).getInfo()\n",
    "    o3 = Atmospheric.ozone(geom,img.date()).getInfo()\n",
    "    aot = Atmospheric.aerosol(geom,img.date()).getInfo()\n",
    "    SRTM = ee.Image('CGIAR/SRTM90_V4')  # Shuttle Radar Topography mission covers *most* of the Earth\n",
    "    alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo()\n",
    "    km = alt/1000                       # i.e. Py6S uses units of kilometers\n",
    "    \n",
    "    # APPLY CORRECTION\n",
    "    \n",
    "    s = SixS()\n",
    "    # Atmospheric constituents\n",
    "    s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3)\n",
    "    s.aero_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical) #altered from- s.aero_profile = AeroProfile.Continental #\n",
    "    s.aot550 = aot\n",
    "    # Earth-Sun-satellite geometry\n",
    "    s.geometry = Geometry.User()\n",
    "    s.geometry.view_z = 0               # always NADIR (I think..)\n",
    "    s.geometry.solar_z = solar_z        # solar zenith angle\n",
    "    s.geometry.month = scene_date.month # month and day used for Earth-Sun distance\n",
    "    s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance\n",
    "    s.altitudes.set_sensor_satellite_level()\n",
    "    s.altitudes.set_target_custom_altitude(km)\n",
    "    \n",
    "    def spectralResponseFunction(bandname):        \n",
    "        \n",
    "        bandSelect = {\n",
    "            'B1':PredefinedWavelengths.S2A_MSI_01,\n",
    "            'B2':PredefinedWavelengths.S2A_MSI_02,\n",
    "            'B3':PredefinedWavelengths.S2A_MSI_03,\n",
    "            'B4':PredefinedWavelengths.S2A_MSI_04,\n",
    "            'B5':PredefinedWavelengths.S2A_MSI_05,\n",
    "            'B6':PredefinedWavelengths.S2A_MSI_06,\n",
    "            'B7':PredefinedWavelengths.S2A_MSI_07,\n",
    "            'B8':PredefinedWavelengths.S2A_MSI_08,\n",
    "            'B8A':PredefinedWavelengths.S2A_MSI_8A,\n",
    "            'B9':PredefinedWavelengths.S2A_MSI_09,\n",
    "            'B10':PredefinedWavelengths.S2A_MSI_10,\n",
    "            'B11':PredefinedWavelengths.S2A_MSI_11,\n",
    "            'B12':PredefinedWavelengths.S2A_MSI_12,\n",
    "                        }\n",
    "        return Wavelength(bandSelect[bandname])\n",
    "    def toa_to_rad(bandname):\n",
    "        ESUN = info['SOLAR_IRRADIANCE_'+bandname]\n",
    "        solar_angle_correction = math.cos(math.radians(solar_z))\n",
    "        # Earth-Sun distance (from day of year)\n",
    "        doy = scene_date.timetuple().tm_yday\n",
    "        d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year\n",
    "        # conversion factor\n",
    "        multiplier = ESUN*solar_angle_correction/(math.pi*d**2)\n",
    "        # at-sensor radiance\n",
    "        rad = img.select(bandname).multiply(multiplier)\n",
    "        return rad\n",
    "    def surface_reflectance(bandname):\n",
    "        # run 6S for this waveband\n",
    "        s.wavelength = spectralResponseFunction(bandname)\n",
    "        s.run()\n",
    "\n",
    "        # extract 6S outputs\n",
    "        Edir = s.outputs.direct_solar_irradiance             #direct solar irradiance\n",
    "        Edif = s.outputs.diffuse_solar_irradiance            #diffuse solar irradiance\n",
    "        Lp   = s.outputs.atmospheric_intrinsic_radiance      #path radiance\n",
    "        absorb  = s.outputs.trans['global_gas'].upward       #absorption transmissivity\n",
    "        scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity\n",
    "        tau2 = absorb*scatter                                #total transmissivity\n",
    "\n",
    "        # radiance to surface reflectance\n",
    "        rad = toa_to_rad(bandname)\n",
    "        ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif))\n",
    "        return ref\n",
    "    \n",
    "    #dblue  = surface_reflectance('B1') #costal aerosol \n",
    "    blue = surface_reflectance('B2')\n",
    "    green = surface_reflectance('B3')\n",
    "    red = surface_reflectance('B4')\n",
    "    #redge1 = surface_reflectance('B5')\n",
    "    #redge2 = surface_reflectance('B6')\n",
    "    #redge3 = surface_reflectance('B7')\n",
    "    #nir0 = surface_reflectance('B8')\n",
    "    nir = surface_reflectance('B8A') # used as NIR due to better compatibility with Landsat <doi.org/10.1016/j.rse.2018.04.031>\n",
    "    #wvap = surface_reflectance('B9')\n",
    "    #cirrus = surface_reflectance('B10')\n",
    "    swir1 = surface_reflectance('B11')\n",
    "    swir2 = surface_reflectance('B12')\n",
    "    #qa = img.select('QA60')\n",
    "    \n",
    "    # to use all bands (-- activate alsdo the bands in the list above)\n",
    "    #boa = dblue.addBands(blue).addBands(green).addBands(red).addBands(redge1).addBands(redge2).addBands(redge3).addBands(nir0).addBands(nir).addBands(wvap).addBands(cirrus).addBands(swir1).addBands(swir2) #.addBands(qa)\n",
    "    # selected bands\n",
    "    boa = blue.addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2)\n",
    "    \n",
    "    # EXPORT FINAL IMAGE\n",
    "    \n",
    "    fname = img.get('DATASTRIP_ID').getInfo()\n",
    "    fname = str(fname [25:33]) + '_' + str(fname [0:3])\n",
    "    #fname = ee.String(img.get('system:index')).getInfo()\n",
    "    task = ee.batch.Export.image.toDrive(**{\n",
    "        'image': boa,\n",
    "        'description': fname + '_BOAzz',\n",
    "        'folder':'P6s_export',\n",
    "        'scale': 20,\n",
    "        'fileFormat': 'GeoTIFF',\n",
    "        'region': region.geometry().bounds(),\n",
    "        'maxPixels': 1e13\n",
    "    })\n",
    "\n",
    "    task.start()\n",
    "    print('exporting ' + fname + '--->done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sentinel 2B atmopheric correction\n",
      "exporting 20190101_S2B--->done\n",
      "exporting 20190111_S2B--->done\n",
      "exporting 20190121_S2B--->done\n",
      "exporting 20190131_S2B--->done\n",
      "exporting 20190210_S2B--->done\n",
      "exporting 20190220_S2B--->done\n",
      "exporting 20190302_S2B--->done\n",
      "exporting 20190312_S2B--->done\n",
      "exporting 20190322_S2B--->done\n",
      "exporting 20190401_S2B--->done\n",
      "exporting 20190411_S2B--->done\n",
      "exporting 20190421_S2B--->done\n",
      "exporting 20190501_S2B--->done\n",
      "exporting 20190511_S2B--->done\n",
      "exporting 20190521_S2B--->done\n",
      "exporting 20190531_S2B--->done\n",
      "exporting 20190610_S2B--->done\n",
      "exporting 20190620_S2B--->done\n",
      "exporting 20190630_S2B--->done\n",
      "exporting 20190710_S2B--->done\n",
      "exporting 20190720_S2B--->done\n",
      "exporting 20190730_S2B--->done\n",
      "exporting 20190809_S2B--->done\n",
      "exporting 20190819_S2B--->done\n",
      "exporting 20190829_S2B--->done\n",
      "exporting 20190908_S2B--->done\n",
      "exporting 20190918_S2B--->done\n",
      "exporting 20190928_S2B--->done\n",
      "exporting 20191008_S2B--->done\n",
      "exporting 20191018_S2B--->done\n",
      "exporting 20191028_S2B--->done\n",
      "exporting 20191107_S2B--->done\n",
      "exporting 20191117_S2B--->done\n",
      "exporting 20191127_S2B--->done\n",
      "exporting 20191207_S2B--->done\n",
      "exporting 20191217_S2B--->done\n",
      "exporting 20191227_S2B--->done\n"
     ]
    }
   ],
   "source": [
    "### SENTINEL 2B COLLECTION ATMOSPHERIC CORRECTION\n",
    "\n",
    "print('Sentinel 2B atmopheric correction')\n",
    "\n",
    "for i in range (0,col_length_S2B):\n",
    "    # OPEN IMAGE AND EXTRACT METADATA ATMOSPHERIC PARAMETERS\n",
    "    \n",
    "    S2_list = S2B_col.toList(col_length_S2B)\n",
    "    img = ee.Image(S2_list.get(i))\n",
    "    \n",
    "    info = img.getInfo()['properties']\n",
    "    scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)\n",
    "    solar_z = img.getInfo()['properties']['MEAN_SOLAR_ZENITH_ANGLE']\n",
    "    h2o = Atmospheric.water(geom,img.date()).getInfo()\n",
    "    o3 = Atmospheric.ozone(geom,img.date()).getInfo()\n",
    "    aot = Atmospheric.aerosol(geom,img.date()).getInfo()\n",
    "    SRTM = ee.Image('CGIAR/SRTM90_V4')  # Shuttle Radar Topography mission covers *most* of the Earth\n",
    "    alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo()\n",
    "    km = alt/1000                       # i.e. Py6S uses units of kilometers\n",
    "    \n",
    "    # APPLY CORRECTION\n",
    "    \n",
    "    s = SixS()\n",
    "    # Atmospheric constituents\n",
    "    s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3)\n",
    "    s.aero_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)  #altered from-s.aero_profile = AeroProfile.Continental\n",
    "    s.aot550 = aot\n",
    "    # Earth-Sun-satellite geometry\n",
    "    s.geometry = Geometry.User()\n",
    "    s.geometry.view_z = 0               # always NADIR (I think..)\n",
    "    s.geometry.solar_z = solar_z        # solar zenith angle\n",
    "    s.geometry.month = scene_date.month # month and day used for Earth-Sun distance\n",
    "    s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance\n",
    "    s.altitudes.set_sensor_satellite_level()\n",
    "    s.altitudes.set_target_custom_altitude(km)\n",
    "    \n",
    "    def spectralResponseFunction(bandname):        \n",
    "        \n",
    "        bandSelect = {\n",
    "            'B1':PredefinedWavelengths.S2B_MSI_01,\n",
    "            'B2':PredefinedWavelengths.S2B_MSI_02,\n",
    "            'B3':PredefinedWavelengths.S2B_MSI_03,\n",
    "            'B4':PredefinedWavelengths.S2B_MSI_04,\n",
    "            'B5':PredefinedWavelengths.S2B_MSI_05,\n",
    "            'B6':PredefinedWavelengths.S2B_MSI_06,\n",
    "            'B7':PredefinedWavelengths.S2B_MSI_07,\n",
    "            'B8':PredefinedWavelengths.S2B_MSI_08,\n",
    "            'B8A':PredefinedWavelengths.S2B_MSI_8A,\n",
    "            'B9':PredefinedWavelengths.S2B_MSI_09,\n",
    "            'B10':PredefinedWavelengths.S2B_MSI_10,\n",
    "            'B11':PredefinedWavelengths.S2B_MSI_11,\n",
    "            'B12':PredefinedWavelengths.S2B_MSI_12,\n",
    "                        }\n",
    "        return Wavelength(bandSelect[bandname])\n",
    "    def toa_to_rad(bandname):\n",
    "        ESUN = info['SOLAR_IRRADIANCE_'+bandname]\n",
    "        solar_angle_correction = math.cos(math.radians(solar_z))\n",
    "        # Earth-Sun distance (from day of year)\n",
    "        doy = scene_date.timetuple().tm_yday\n",
    "        d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year\n",
    "        # conversion factor\n",
    "        multiplier = ESUN*solar_angle_correction/(math.pi*d**2)\n",
    "        # at-sensor radiance\n",
    "        rad = img.select(bandname).multiply(multiplier)\n",
    "        return rad\n",
    "    def surface_reflectance(bandname):\n",
    "        # run 6S for this waveband\n",
    "        s.wavelength = spectralResponseFunction(bandname)\n",
    "        s.run()\n",
    "\n",
    "        # extract 6S outputs\n",
    "        Edir = s.outputs.direct_solar_irradiance             #direct solar irradiance\n",
    "        Edif = s.outputs.diffuse_solar_irradiance            #diffuse solar irradiance\n",
    "        Lp   = s.outputs.atmospheric_intrinsic_radiance      #path radiance\n",
    "        absorb  = s.outputs.trans['global_gas'].upward       #absorption transmissivity\n",
    "        scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity\n",
    "        tau2 = absorb*scatter                                #total transmissivity\n",
    "\n",
    "        # radiance to surface reflectance\n",
    "        rad = toa_to_rad(bandname)\n",
    "        ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif))\n",
    "        return ref\n",
    "    \n",
    "    #dblue  = surface_reflectance('B1') #costal aerosol \n",
    "    blue = surface_reflectance('B2')\n",
    "    green = surface_reflectance('B3')\n",
    "    red = surface_reflectance('B4')\n",
    "    #redge1 = surface_reflectance('B5')\n",
    "    #redge2 = surface_reflectance('B6')\n",
    "    #redge3 = surface_reflectance('B7')\n",
    "    #nir0 = surface_reflectance('B8')\n",
    "    nir = surface_reflectance('B8A') # used as NIR due to better compatibility with Landsat <doi.org/10.1016/j.rse.2018.04.031>\n",
    "    #wvap = surface_reflectance('B9')\n",
    "    #cirrus = surface_reflectance('B10')\n",
    "    swir1 = surface_reflectance('B11')\n",
    "    swir2 = surface_reflectance('B12')\n",
    "    #qa = img.select('QA60')\n",
    "    \n",
    "    # to use all bands (-- activate alsdo the bands in the list above)\n",
    "    #boa = dblue.addBands(blue).addBands(green).addBands(red).addBands(redge1).addBands(redge2).addBands(redge3).addBands(nir0).addBands(nir).addBands(wvap).addBands(cirrus).addBands(swir1).addBands(swir2) #.addBands(qa)\n",
    "    # selected bands\n",
    "    boa = blue.addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2)\n",
    "    \n",
    "    # EXPORT FINAL IMAGE\n",
    "    \n",
    "    fname = img.get('DATASTRIP_ID').getInfo()\n",
    "    fname = str(fname [25:33]) + '_' + str(fname [0:3])\n",
    "    #fname = ee.String(img.get('system:index')).getInfo() ### adicionado\n",
    "    task = ee.batch.Export.image.toDrive(**{\n",
    "        'image': boa,\n",
    "        'description': fname + '_BOA',\n",
    "        'folder':'P6s_export',\n",
    "        'scale': 20,\n",
    "        'fileFormat': 'GeoTIFF',\n",
    "        'region': region.geometry().bounds(),\n",
    "        'maxPixels': 1e13\n",
    "    })\n",
    "\n",
    "    task.start()\n",
    "    print('exporting ' + fname + '--->done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Landsat 8 atmopheric correction\n",
      "exporting 20190125_L08--->done\n",
      "exporting 20190314_L08--->done\n",
      "exporting 20190501_L08--->done\n",
      "exporting 20190517_L08--->done\n",
      "exporting 20190618_L08--->done\n",
      "exporting 20190704_L08--->done\n",
      "exporting 20190720_L08--->done\n",
      "exporting 20190805_L08--->done\n",
      "exporting 20190821_L08--->done\n",
      "exporting 20190906_L08--->done\n",
      "exporting 20191008_L08--->done\n",
      "exporting 20191109_L08--->done\n",
      "exporting 20191125_L08--->done\n",
      "exporting 20191227_L08--->done\n"
     ]
    }
   ],
   "source": [
    "### LANDSAT 8 COLLECTION ATMOSPHERIC CORRECTION\n",
    "\n",
    "print('Landsat 8 atmopheric correction')\n",
    "\n",
    "for i in range (0,col_length_L8):\n",
    "    # OPEN IMAGE AND EXTRACT METADATA ATMOSPHERIC PARAMETERS\n",
    "    \n",
    "    L8_list = L8_col.toList(col_length_L8)\n",
    "    img = ee.Image(L8_list.get(i))\n",
    "\n",
    "    info = img.getInfo()['properties']\n",
    "    scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)\n",
    "    solar_elv = img.getInfo()['properties']['SUN_ELEVATION']\n",
    "    solar_z = ee.Number(90).subtract(solar_elv).getInfo()\n",
    "    h2o = Atmospheric.water(geom,img.date()).getInfo()\n",
    "    o3 = Atmospheric.ozone(geom,img.date()).getInfo()\n",
    "    aot = Atmospheric.aerosol(geom,img.date()).getInfo()\n",
    "    SRTM = ee.Image('CGIAR/SRTM90_V4') # Shuttle Radar Topography mission covers *most* of the Earth\n",
    "    alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo()\n",
    "    km = alt/1000                      # i.e. Py6S uses units of kilometers\n",
    "    \n",
    "    # APPLY CORRECTION\n",
    "    \n",
    "    s = SixS()\n",
    "    # Atmospheric constituents\n",
    "    s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3)\n",
    "    s.aero_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical) #altered from- s.aero_profile = AeroProfile.Continental\n",
    "    s.aot550 = aot\n",
    "    # Earth-Sun-satellite geometry\n",
    "    s.geometry = Geometry.User()\n",
    "    s.geometry.view_z = 0               # always NADIR (I think..)\n",
    "    s.geometry.solar_z = solar_z        # solar zenith angle\n",
    "    s.geometry.month = scene_date.month # month and day used for Earth-Sun distance\n",
    "    s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance\n",
    "    s.altitudes.set_sensor_satellite_level()\n",
    "    s.altitudes.set_target_custom_altitude(km)\n",
    "    \n",
    "    def spectralResponseFunction(bandname):        \n",
    "        bandSelect = {\n",
    "            'B1':PredefinedWavelengths.LANDSAT_OLI_B1,\n",
    "            'B2':PredefinedWavelengths.LANDSAT_OLI_B2,\n",
    "            'B3':PredefinedWavelengths.LANDSAT_OLI_B3,\n",
    "            'B4':PredefinedWavelengths.LANDSAT_OLI_B4,\n",
    "            'B5':PredefinedWavelengths.LANDSAT_OLI_B5,\n",
    "            'B6':PredefinedWavelengths.LANDSAT_OLI_B6,\n",
    "            'B7':PredefinedWavelengths.LANDSAT_OLI_B7,\n",
    "            'B8':PredefinedWavelengths.LANDSAT_OLI_B8,\n",
    "            'B9':PredefinedWavelengths.LANDSAT_OLI_B9,\n",
    "                    }\n",
    "        return Wavelength(bandSelect[bandname])\n",
    "    def toa_to_rad(bandname):\n",
    "        ESUN_L8 = [1895.33, 2004.57, 1820.75, 1549.49, 951.76, 247.55, 85.46, 1723.8, 366.97]\n",
    "        ESUN_BAND = {\n",
    "            'B1':ESUN_L8[0],\n",
    "            'B2':ESUN_L8[1],\n",
    "            'B3':ESUN_L8[2],\n",
    "            'B4':ESUN_L8[3],\n",
    "            'B5':ESUN_L8[4],\n",
    "            'B6':ESUN_L8[5],\n",
    "            'B7':ESUN_L8[6],\n",
    "            'B8':ESUN_L8[7],\n",
    "            'B9':ESUN_L8[8],\n",
    "            }\n",
    "        solar_angle_correction = math.cos(math.radians(solar_z))\n",
    "        # Earth-Sun distance (from day of year)\n",
    "        doy = scene_date.timetuple().tm_yday\n",
    "        d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year\n",
    "        # conversion factor\n",
    "        multiplier = ESUN_BAND[bandname]*solar_angle_correction/(math.pi*d**2)\n",
    "        # at-sensor radiance\n",
    "        rad = img.select(bandname).multiply(multiplier)\n",
    "        return rad\n",
    "    def surface_reflectance(bandname):\n",
    "        # run 6S for this waveband\n",
    "        s.wavelength = spectralResponseFunction(bandname)\n",
    "        s.run()\n",
    "\n",
    "        # extract 6S outputs\n",
    "        Edir = s.outputs.direct_solar_irradiance             #direct solar irradiance\n",
    "        Edif = s.outputs.diffuse_solar_irradiance            #diffuse solar irradiance\n",
    "        Lp   = s.outputs.atmospheric_intrinsic_radiance      #path radiance\n",
    "        absorb  = s.outputs.trans['global_gas'].upward       #absorption transmissivity\n",
    "        scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity\n",
    "        tau2 = absorb*scatter                                #total transmissivity\n",
    "\n",
    "        # radiance to surface reflectance\n",
    "        rad = toa_to_rad(bandname)\n",
    "        ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif))\n",
    "        return ref\n",
    "    \n",
    "    #dblue  = surface_reflectance('B1') #costal aerosol \n",
    "    blue = surface_reflectance('B2')\n",
    "    green = surface_reflectance('B3')\n",
    "    red = surface_reflectance('B4')\n",
    "    nir = surface_reflectance('B5')\n",
    "    swir1 = surface_reflectance('B6')\n",
    "    swir2 = surface_reflectance('B7')\n",
    "    #pan = surface_reflectance('B8')\n",
    "    #cirrus = surface_reflectance('B9')\n",
    "    #qa = img.select('BQA')\n",
    "    \n",
    "    # to use all bands (-- activate alsdo the bands in the list above)\n",
    "    #boa = dblue.addBands(blue).addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2) #.addBands(qa)\n",
    "    # selected bands\n",
    "    boa = blue.addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2)\n",
    "    \n",
    "    # EXPORT FINAL IMAGE\n",
    "    \n",
    "    fname = ee.String(img.get('system:index')).getInfo()\n",
    "    fname = str(fname [12:20]) + '_L08'\n",
    "\n",
    "    task = ee.batch.Export.image.toDrive(**{\n",
    "        'image': boa,\n",
    "        'description': fname + '_BOA',\n",
    "        'folder':'P6s_export',\n",
    "        'scale': 30,\n",
    "        'fileFormat': 'GeoTIFF',\n",
    "        'region': region.geometry().bounds(),\n",
    "        'maxPixels': 1e13\n",
    "    })\n",
    "\n",
    "    task.start()\n",
    "    print('exporting ' + fname + '--->done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Landsat 7 atmopheric correction\n",
      "exporting 20190117_L07--->done\n",
      "exporting 20190306_L07--->done\n",
      "exporting 20190322_L07--->done\n",
      "exporting 20190509_L07--->done\n",
      "exporting 20190525_L07--->done\n",
      "exporting 20190610_L07--->done\n",
      "exporting 20190626_L07--->done\n",
      "exporting 20190712_L07--->done\n",
      "exporting 20190728_L07--->done\n",
      "exporting 20190813_L07--->done\n",
      "exporting 20190829_L07--->done\n",
      "exporting 20190914_L07--->done\n",
      "exporting 20190930_L07--->done\n",
      "exporting 20191101_L07--->done\n",
      "exporting 20191203_L07--->done\n"
     ]
    }
   ],
   "source": [
    "### LANDSAT 7 COLLECTION ATMOSPHERIC CORRECTION\n",
    "\n",
    "print('Landsat 7 atmopheric correction')\n",
    "\n",
    "for i in range (0,col_length_L7):\n",
    "    # OPEN IMAGE AND EXTRACT METADATA ATMOSPHERIC PARAMETERS\n",
    "    \n",
    "    L7_list = L7_col.toList(col_length_L7)\n",
    "    img = ee.Image(L7_list.get(i))\n",
    "\n",
    "    info = img.getInfo()['properties']\n",
    "    scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)\n",
    "    solar_elv = img.getInfo()['properties']['SUN_ELEVATION']\n",
    "    solar_z = ee.Number(90).subtract(solar_elv).getInfo()\n",
    "    h2o = Atmospheric.water(geom,img.date()).getInfo()\n",
    "    o3 = Atmospheric.ozone(geom,img.date()).getInfo()\n",
    "    aot = Atmospheric.aerosol(geom,img.date()).getInfo()\n",
    "    SRTM = ee.Image('CGIAR/SRTM90_V4')  # Shuttle Radar Topography mission covers *most* of the Earth\n",
    "    alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo()\n",
    "    km = alt/1000                       # i.e. Py6S uses units of kilometers\n",
    "    \n",
    "    # APPLY CORRECTION\n",
    "    \n",
    "    s = SixS()\n",
    "    # Atmospheric constituents\n",
    "    s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3)\n",
    "    s.aero_profile = AeroProfile.Continental # altered from: s.aero_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)\n",
    "    s.aot550 = aot\n",
    "    # Earth-Sun-satellite geometry\n",
    "    s.geometry = Geometry.User()\n",
    "    s.geometry.view_z = 0               # always NADIR (I think..)\n",
    "    s.geometry.solar_z = solar_z        # solar zenith angle\n",
    "    s.geometry.month = scene_date.month # month and day used for Earth-Sun distance\n",
    "    s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance\n",
    "    s.altitudes.set_sensor_satellite_level()\n",
    "    s.altitudes.set_target_custom_altitude(km)\n",
    "    \n",
    "    def spectralResponseFunction(bandname):        \n",
    "        bandSelect = {\n",
    "            'B1':PredefinedWavelengths.LANDSAT_ETM_B1,\n",
    "            'B2':PredefinedWavelengths.LANDSAT_ETM_B2,\n",
    "            'B3':PredefinedWavelengths.LANDSAT_ETM_B3,\n",
    "            'B4':PredefinedWavelengths.LANDSAT_ETM_B4,\n",
    "            'B5':PredefinedWavelengths.LANDSAT_ETM_B5,\n",
    "            'B7':PredefinedWavelengths.LANDSAT_ETM_B7,\n",
    "                    }\n",
    "        return Wavelength(bandSelect[bandname])\n",
    "    def toa_to_rad(bandname):\n",
    "        ESUN_L8 = [1895.33, 2004.57, 1820.75, 1549.49, 951.76, 247.55, 85.46, 1723.8, 366.97]\n",
    "        ESUN_L7 = [1997, 1812, 1533, 1039, 230.8, 84.9] # PAN =  1362 (removed to match Py6S)\n",
    "        ESUN_BAND = {\n",
    "            'B1':ESUN_L7[0],\n",
    "            'B2':ESUN_L7[1],\n",
    "            'B3':ESUN_L7[2],\n",
    "            'B4':ESUN_L7[3],\n",
    "            'B5':ESUN_L7[4],            \n",
    "            'B7':ESUN_L7[5],\n",
    "           }\n",
    "        solar_angle_correction = math.cos(math.radians(solar_z))\n",
    "        # Earth-Sun distance (from day of year)\n",
    "        doy = scene_date.timetuple().tm_yday\n",
    "        d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year\n",
    "        # conversion factor\n",
    "        multiplier = ESUN_BAND[bandname]*solar_angle_correction/(math.pi*d**2)\n",
    "        # at-sensor radiance\n",
    "        rad = img.select(bandname).multiply(multiplier)\n",
    "        return rad\n",
    "    def surface_reflectance(bandname):\n",
    "        # run 6S for this waveband\n",
    "        s.wavelength = spectralResponseFunction(bandname)\n",
    "        s.run()\n",
    "\n",
    "        # extract 6S outputs\n",
    "        Edir = s.outputs.direct_solar_irradiance             #direct solar irradiance\n",
    "        Edif = s.outputs.diffuse_solar_irradiance            #diffuse solar irradiance\n",
    "        Lp   = s.outputs.atmospheric_intrinsic_radiance      #path radiance\n",
    "        absorb  = s.outputs.trans['global_gas'].upward       #absorption transmissivity\n",
    "        scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity\n",
    "        tau2 = absorb*scatter                                #total transmissivity\n",
    "\n",
    "        # radiance to surface reflectance\n",
    "        rad = toa_to_rad(bandname)\n",
    "        ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif))\n",
    "        return ref\n",
    "    \n",
    "    blue = surface_reflectance('B1')\n",
    "    green = surface_reflectance('B2')\n",
    "    red = surface_reflectance('B3')\n",
    "    nir = surface_reflectance('B4')    \n",
    "    swir1 = surface_reflectance('B5')\n",
    "    swir2 = surface_reflectance('B7')\n",
    "        \n",
    "    boa = blue.addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2) ###\n",
    "    \n",
    "    \n",
    "    # EXPORT FINAL IMAGE\n",
    "    fname = ee.String(img.get('system:index')).getInfo()\n",
    "    fname = str(fname [12:20]) + '_L07'\n",
    "    \n",
    "    task = ee.batch.Export.image.toDrive(**{\n",
    "        'image': boa,\n",
    "        'description': fname + '_BOA',\n",
    "        'folder':'P6s_export',\n",
    "        'scale': 30,\n",
    "        'fileFormat': 'GeoTIFF',\n",
    "        'region': region.geometry().bounds(),\n",
    "        'maxPixels': 1e13\n",
    "    })\n",
    "\n",
    "    task.start()\n",
    "    print('exporting ' + fname + '--->done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "## REFERENCES ##\n",
    "\n",
    "# -- PY6s Python API\n",
    "# WILSON, R. T., Py6S: A Python interface to the 6S radiative transfer model, Computers and Geosciences, 51, p.166-171, 2013.\n",
    "\n",
    "# -- Second Simulation of the Satellite Signal in the Solar Spectrum (6S) method \n",
    "# VERMOTE, E. F.; TANRÉ, D.; DEUZÉ, J. L.; HERMAN, M.; MORCRETTE, J.-J. Second simulation of the satellite signal in the solar spectrum, 6S:\n",
    "# an overview. IEEE Transactions on Geoscience and Remote Sensing, v. 35, n. 3, p. 675–686, 1997. \n",
    "# VERMOTE, E. F.; EL SALEOUS, N.; JUSTICE, C. O.; KAUFMAN, Y. J.; PRIVETTE, J. L.; REMER, L.; ROGER, J. C.; TANRÉ, D. Atmospheric \n",
    "# correction of visible to middle-infrared EOS-MODIS data over land surfaces: Background, operational algorithm and validation. Journal \n",
    "# of Geophysical Research, v. 102, p. 17131, 1997. \n",
    "\n",
    "# -- Harmonization of Landsat-Sentinel2 example, in the context of crop monitoring (GEE scripts available, used as reference here)\n",
    "# NGUYEN, M. D.; BAEZ-VILLANUEVA, O. M.; BUI, D. D.; NGUYEN, P. T.; RIBBE, L. Harmonization of Landsat and Sentinel 2 for crop monitoring \n",
    "# in drought prone areas: Case studies of Ninh Thuan (Vietnam) and Bekaa (Lebanon). Remote Sensing, v. 12, n. 2, p. 1–18, 2020."
   ]
  }
 ],
 "metadata": {
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}