samsammurphy/gee-atmcorr-s2

Atmospheric correction of a (single) Sentinel 2 image

Introduction

Atmospheric correction of Sentinel 2 imagery in Google Earth Engine using Py6S.

Installation

Recommended: Docker

The following docker container has all dependencies to run the code in this repository

docker pull samsammurphy/ee-python3-jupyter-atmcorr:v1.2

You can also build a docker container from scratch using this Dockerfile from the gee-community repo.

Alternative: Manual installation

This repo has the following prerequisites

• Python 3.x
• Google Earth Engine Python API
• Jupyter
• Py6S

Usage

To run the docker container with access to a web browser

$ docker run -i -t -p 8888:8888 samsammurphy/ee-python3-jupyter-atmcorr:v1.2

Once inside the container, authenticate Earth Engine

earthengine authenticate

then grab the source code for this repo

git clone https://github.com/samsammurphy/gee-atmcorr-S2

and run the example jupyter notebook

cd gee-atmcorr-S2/jupyer_notebooks/
jupyter-notebook sentinel2_atmospheric_correction.ipynb --ip='*' --port=8888 --allow-root

this will print out a URL that you can copy/paste into your web browser to run the code.

If the URL is http://(something_in_parentheses) then you will need to change the parentheses and its contents for localhost. A valid URL should look something like..

http://localhost:8888/?token=...

Saving authentication

After authenticating with earthengine and cloning the repository you can save the changes you've made to the container with docker commit. Docker commit will create a new image from a running containers current state.

With the container still running, open a new terminal to get the container's ID:

docker ps

copy the ID to clipboard and run

docker commit [ID] gee-atmcorr-s2:myauth

to commit the image. Now if you run

docker images

your newly committed image should be at the top of the list.

You can now start the note book with

docker run -i -t -p 8888:8888 gee-atmcorr-s2:myauth 

jupyter-notebook /gee-atmcorr-S2/jupyer_notebooks/sentinel2_atmospheric_correction.ipynb --ip='*' --port=8888 --allow-root

gee-atmcorr-s2's Issues

how to run this program

Hi, Dr. Murphy, I am a student from China, I had download this program, but I can not run it correct, please help me how to run it correct, thank you very much

Problem in function spectralResponseFunction()

We need to change the function to the following. Otherwise, an error will be reported 'PredefinedWavelengths' has no attribute 'S2A_MSI_13, because there is no B13 in S2A_MSI SRF (see the Py6S source code)

def spectralResponseFunction(bandname):

bandSelect = {
    'B1':PredefinedWavelengths.S2A_MSI_01,
    'B2':PredefinedWavelengths.S2A_MSI_02,
    'B3':PredefinedWavelengths.S2A_MSI_03,
    'B4':PredefinedWavelengths.S2A_MSI_04,
    'B5':PredefinedWavelengths.S2A_MSI_05,
    'B6':PredefinedWavelengths.S2A_MSI_06,
    'B7':PredefinedWavelengths.S2A_MSI_07,
    'B8':PredefinedWavelengths.S2A_MSI_08,
    'B8A':PredefinedWavelengths.S2A_MSI_8A,
    'B9':PredefinedWavelengths.S2A_MSI_09,
    'B10':PredefinedWavelengths.S2A_MSI_10,
    'B11':PredefinedWavelengths.S2A_MSI_11,
    'B12':PredefinedWavelengths.S2A_MSI_12,
    }
return Wavelength(bandSelect[bandname])

Where to find atmospheric library?

Hi author:

You import "atmospheric" library and use it to obtain aot parameters. But where to find this library?

Thanks,
Hanze

Sentinel 3

Although the repo title may answer this question but I wanted to ask just to make sure. Is this applicable to Sentinel 3? Also what is the current condition of this? Is this properly working (properly doing atmospheric correction of S2 in gee)?

I am unable to try this yet due to not being successful in installing the Python API of gee. Thanks.

add docker instructions

move current installation instructions to 'manual installation' wiki

update docker image

needs updated version of Earth Engine Python API and testing

Want to participate in subsequent development

Want to participate in subsequent development, anyone want to be together?

docker commit to save authentication

write a wiki on how to use docker commit to save earthengine authentication details

update docker image

needs updated version of Earth Engine Python API and testing

Correction AOI

Hi,
I see in your code that you calculate some atmospheric values (ozone, water vapor, etc.) at the centroid of the area of interest (AOI). Does it mean that the corrected image will be accurate only for centroid? What if the AOI is as big as the whole Sentinel-2 scene? Does the P6yS algorithm plays a role here?

Best Regards
Shahriar

Layer error: Image.constant: Parameter 'value' is required.

Hi Sam,
Thank you for providing amazing repo for image calibration in GEE. All was good till now but with reduced image collections i.e. mean median or mosaic image, google earth engine shows this error.
"Layer error: Image.constant: Parameter 'value' is required."

It only works with a single image which is filtered by ID or .first() method, not the mean median or mosaic image. :(

Please check this out
https://code.earthengine.google.com/e74efdb28503f7d942d90d831c337e28

Help needed for OLCI sentinel 3

I hope you are doing well. I am currently working with OLCI in Google Earth Engine and have a few questions:
My study area near coast. does not cover srtm dem.
I have set up the following link https://py6s.readthedocs.io/en/latest/installation.html
I am trying to convert radiance to reflectance, here is the code, please review it.

import ee
from Py6S import *
import datetime
import math
import os
import sys
import pytz
from pysolar.solar import get_altitude

# Initialize Earth Engine
ee.Initialize()

# Define the dates to download
datesstart = ['2024-03-10', '2024-03-11', '2024-03-12', '2024-03-13', '2024-06-24', '2024-06-25', '2024-06-26']

# Define the study area
study_area = ee.FeatureCollection('projects/ee-wfahydrologymsu100/assets/WMS')
geometry = study_area.geometry()

# Function to calculate solar zenith angle
def calculate_solar_z(geometry, scene_date):
   latitude = geometry.centroid().coordinates().get(1).getInfo()
   longitude = geometry.centroid().coordinates().get(0).getInfo()
   solar_altitude = get_altitude(latitude, longitude, scene_date)
   solar_z = 90 - solar_altitude
   return solar_z

# Function to calculate surface reflectance
def surface_reflectance(bandname, image):
   # Instantiate the 6S model
   s = SixS(r"C:\Users\hafez\build\6SV\1.1\6SV1.1\sixsV1.1")
   
   # Extract image metadata
   info = image.getInfo()['properties']
   scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start'] / 1000)
   scene_date = scene_date.replace(tzinfo=pytz.UTC)
   
   # Calculate solar zenith angle
   solar_z = calculate_solar_z(geometry, scene_date)
   
   # Get atmospheric data
   h2o = Atmospheric.water(geometry, ee.Date(scene_date)).getInfo()
   o3 = Atmospheric.ozone(geometry, ee.Date(scene_date)).getInfo()
   aot = Atmospheric.aerosol(geometry, ee.Date(scene_date)).getInfo()
   
   # Set the atmospheric constituents
   s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o, o3)
   s.aero_profile = AeroProfile.Continental
   s.aot550 = aot
   
   # Set the Earth-Sun-satellite geometry
   s.geometry = Geometry.User()
   s.geometry.view_z = 0               # always NADIR
   s.geometry.solar_z = solar_z        # solar zenith angle
   s.geometry.month = scene_date.month # month and day used for Earth-Sun distance
   s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance
   s.altitudes.set_sensor_satellite_level()
   
   # Define the spectral response function for Sentinel-3 OLCI
   def spectralResponseFunction(bandname):
       bandSelect = {
           'Oa02_radiance': PredefinedWavelengths.S3A_OLCI_02,
           'Oa03_radiance': PredefinedWavelengths.S3A_OLCI_03,
           'Oa04_radiance': PredefinedWavelengths.S3A_OLCI_04,
           'Oa05_radiance': PredefinedWavelengths.S3A_OLCI_05,
           'Oa06_radiance': PredefinedWavelengths.S3A_OLCI_06,
           'Oa07_radiance': PredefinedWavelengths.S3A_OLCI_07,
           'Oa08_radiance': PredefinedWavelengths.S3A_OLCI_08,
           'Oa09_radiance': PredefinedWavelengths.S3A_OLCI_09,
           'Oa11_radiance': PredefinedWavelengths.S3A_OLCI_11,
           'Oa12_radiance': PredefinedWavelengths.S3A_OLCI_12,
           'Oa13_radiance': PredefinedWavelengths.S3A_OLCI_13,
       }
       return Wavelength(bandSelect[bandname])

   def toa_to_rad(bandname, image):
       rad = image.select(bandname)
       return rad

   # Run 6S for this waveband
   s.wavelength = spectralResponseFunction(bandname)
   s.run()
   
   # Extract  Edir 
   Edir = s.outputs.direct_solar_irradiance  # Direct solar irradiance
   Edif = s.outputs.diffuse_solar_irradiance  # Diffuse solar irradiance
   Lp = s.outputs.atmospheric_intrinsic_radiance  # Path radiance
   absorb = s.outputs.trans['global_gas'].upward  # Absorption transmissivity
   scatter = s.outputs.trans['total_scattering'].upward  # Scattering transmissivity
   tau2 = absorb * scatter  # Total transmissivity
   
   # Radiance to surface reflectance
   rad = toa_to_rad(bandname, image)
   ref = rad.subtract(Lp).multiply(math.pi).divide(tau2 * (Edir + Edif))
   
   return ref

# Process each date and download the image
for date in datesstart:
   date_obj = ee.Date(date)
   asset_name = 'S3OLCI_' + date_obj.format('YYYY_MM_dd').getInfo()
   
   # Get the Sentinel-3 OLCI image
   S3_image = ee.Image(
       ee.ImageCollection('COPERNICUS/S3/OLCI')
       .filterBounds(geometry)
       .filterDate(date_obj, date_obj.advance(1, 'day'))
       .sort('system:time_start')
       .first()
   )
   
   # Calculate surface reflectance for each band
   bands = ['Oa02_radiance', 'Oa03_radiance', 'Oa04_radiance', 'Oa05_radiance', 'Oa06_radiance', 'Oa07_radiance',
            'Oa08_radiance', 'Oa09_radiance', 'Oa11_radiance', 'Oa12_radiance', 'Oa13_radiance']
   
   surface_reflectance_images = [surface_reflectance(band, S3_image) for band in bands]
   
   # Combine all surface reflectance bands into a single image
   combined_reflectance = ee.Image.cat(surface_reflectance_images)
   
   # Clip the image to the study area geometry
   clipped_reflectance = combined_reflectance.clip(geometry)
   
   # Define the export parameters
   export_params = {
       'image': clipped_reflectance.toFloat(),
       'description': asset_name,
       'scale': 10,  # Adjust the scale as needed
       'region': geometry,
       'fileFormat': 'GeoTIFF'
   }
   
   # Export the image to Google Drive
   task = ee.batch.Export.image.toDrive(**export_params)
   task.start()
   
   # Print the export task status
   print(f'Exporting {asset_name}...')

print('All exports queued successfully.')

`


Best,
Hafez