Atmospheric Correction of Satellite images using Python
This article will help you write the python code to extract the reflectance from DN(Digital Number) of the Satellite images. For this tutorial, we will use images Landsat of 7 and Landsat 8 both. You can get the whole code here.
Why can’t you use the Satellite images Directly?
Irradiance originates from two sources :
- Directly reflected sunlight
- Diffused Skylight
- The relative dominance of sunlight versus skylight in any given image
is strongly dependent on weather conditions (e.g., sunny vs. hazy vs.
Cloudy). - Irradiance varies with the seasonal changes in solar elevation angle and
the changing distance between the earth and sun.
- With a sensor positioned close to the earth’s surface, the path radiance will generally be small or negligible.
- Imagery from satellite systems will be more strongly affected by the path radiance due to the longer atmospheric path between the earth’s surface and the spacecraft.
Due to the thickness of the atmosphere between the earth’s surface and
the satellite’s position above the atmosphere, this second spectral response
pattern shows an elevated signal at short wavelengths due to the
extraneous path radiance. In raw form, the near-surface measurement from the eld spectro-radiometer could not be directly compared to the satellite measurement because one observes surface reflectance while the other is observing the so-called top of atmosphere (TOA) reflectance. Before such a comparison, the satellite image needs to be processed for atmospheric correction. The raw spectral data are modified to compensate for the expected effects of atmospheric scattering and absorption.
How to Get the relevant data?
Well, for this tutorial, I got you covered, but if you really want to do your own experiment, You are free to explore USGS Dataset.
First of all, add this zip file to your google drive inside a folder that says “Geoinformatics Now Open Google Colab and create a New Notebook. Now from the left-hand side panel, ‘mount your Google Drive.’
Once it is done, you will see a folder named ‘drive’ there. Now go and locate the folder where you have saved the dataset in the folder.
Now once you have located the Data.zip file, right-click on it and copy the path. Now run the following command to unzip the dataset.
# unzipping the data.zip file!unzip /content/drive/MyDrive/Geoinformatics/Lab/Data.zip
Now you will find 3 unzipped files in the root directory. One is the Bhopal_Mask.zip, and the other two are the folder that contains satellite images from Landsat 7 and Landsat 8. Now we will extract these file too
#creating the folder named L7 and extracting the landsat7 data into it.!mkdir L7!tar -xvf /content/LE07_L1TP_145044_20021210_20170127_01_T1.tar.gz -C /content/L7#creating the folder named L8 and extracting the landsat8 data into it.
!mkdir L8
!tar -xvf /content/LC08_L1TP_146044_20201226_20201226_01_RT.tar.gz -C /content/L8# Unzipping the bhopal mask.zip
!unzip Bhopal_Mask.zip
By now Data Collection part is done.
What are the Python Libraries we gonna use?
Install the required Libraries with the following code
# installing the all the required library!pip install elevation
!pip install richdem
!pip install pysheds
!pip install plantcv
!apt install imagemagick
!pip install fiona
!pip install rasterio
!pip install shapely
!pip install pyproj
!pip install geopandas
!pip install rioxarray
!pip install earthpy
Then import the libraries:
# importing the required librariesfrom osgeo import gdal
import matplotlib
import elevation
from matplotlib.pyplot import figure
import os
import earthpy as et
import fiona
import rasterio
import rasterio.mask
import matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
from rasterio.crs import CRS
import rioxarray as rxr
import earthpy as et
import geopandas as gpd
import xarray as xr
import rioxarray as rxr
import earthpy.spatial as es
import earthpy.plot as ep
Working With Landsat 7 images
Atmospheric Correction has two steps involved:
- DN to Radiance
- Radiance to Reflectance
So let's dive into DN to Radiance Conversion. We have two methods for that.
- Gain and Bias Method
- Spectral Radiance Scaling
We will see both techniques. Also, both techniques give the same result.
DN to Radiance
Gain And Bias method:
Now the main question that arises is, is that how to get the G and B values?
So for each band of the different satellites, you have another value for G and B., and all these values can be found in the meta-datasheet of the satellite. SO if you look in the L7 folder that we created, you will find a “…..MTL.txt” file inside it. That is a metadata file. In that file, G is stored as “RADIANCE_MULT_BAND_<number>,” and B is stored as “RADIANCE_ADD_BAND_<number>.” SO if you want, you can set the value of each band and each satellite manually, or you could be smarter like me and do it in a python way.
This will create a dictionary of that metadata file. If you followed me this far, it seems like you are interested in this article.
Why do masking?
Masking is required because the images we are using are enormous and. If you operate on these images, you will probably end up with a Memory overflow error. SO make our code crash-free and fast. We have to apply masking. One more reason for using masking is when you only want to study a particular area in the image. Then we mask the other part of the picture.
Radiance to Reflectance
Now you have everything you need now. If you plot the corrected images, you will get something like this.
Spectral Radiance method
Similar to the gain and bias method, we will get the values of Lmin, Lmax…. form the meta-datafile.
Now let's plot the corrected images using the Spectral radiance method.
The thing to notice is that images created by both methods are the same, so it does not matter which way you choose. Both will give you the same results. Also, for radiance to reflectance, we have only one method, and for DN to radiance, we have two ways.
Now, don’t you want to see the RGB composite image?
False-color composite image:
Now, if you perform similar steps for Landsat8 images, you can get the results very quickly.
Your final images will be something like this:
RGB image:
False-color Composition:
Get the whole code from following GitHub or colab notebook.
