Rayleigh Correction for Sentinel data using SNAP

Removing atmospheric effects from satellite images is tricky, since the theoretical foundations are not completed understood, it is difficult to model the physical conditions, and also because we generally lack many information regarding the atmospheric condition at the time and location of satellite image acquisition.

One procedure that attempts to partially remove atmospheric effects is implemented in the SNAP software application, that is released by ESA and mainly devoted to the processing of Sentinel data, both optical and radar.

This method is named "Rayleigh Correction" and attempts to mitigate gaseous and Rayleigh effects. It does however do not consider aerosol effects.

In the screenshot below you see how to open the method window.

The GUI allows to select the bands to be corrected (B2-4 in the example below), besides other parameters (e.g., sea-level pressure). For this example the output resolution was changed to 10 meters, to match that of input bands (B2-4). The other parameters were left unchanged.

The two windows below display a an example of Sentinel 2 data (S2A MSI, 2017-05-27) pre- and post- Rayleigh correction. The displayed zone is comprised between Chivasso to the West and Verolengo to the East, and the Po River to the South and the Canale Cavour to the North (Piedmont, Italy). 

In the upper window the original image is displayed as a 4-3-2 composite. In the lower window the same zone is represented using a 4-3-2 composite of Rayleigh-corrected bands.

The original bands are darker and with a low contrast, while the corrected bands present a larger chromatic variance. Note that the lightest pixels in the original band tend to saturate to white, somewhat loosing information.

Water, vegetation and urban constructions are more vivid and looks like more "natural", with more chromatic variations than in the original bands.

Quick calculation of band indices from satellite images with GDAL

Band indices of satellite images, such as the NDVI, can be routinely calculated from within applications, such as ENVI, SNAP and so on.

Sometimes it is preferable to calculate band indices using a programming language: it could be faster than using applications, we could not dispose of the required applications or we need to automate a set of processing by using a programming language.

Having GDAL installed, it is quick and easy to calculate band indices (and a lot more) directly from the prompt/shell. I was drawn to this useful solution by the answer provided by

The used command is gdal_calc.py and we can specify one, two or more input bands to perform calculations, for instance, for calculating the Normalized Difference Water Index (NDWI) from a spatial subset of a Sentinel 2 image saved from SNAP as BEAM-Dimap, we use:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))"

Obviously the shell working directory is the folder containing the source image files. Morevover we do not have to define the used band within the input .img file since there is only one band.
The complete description of the command is available at the page 'gdal_calc.py'.


An example of calculation of the Normalized Difference Water Index is here presented. Three Sentinel 2 SLC image of the Po river in the Piedmont area between Chivasso and Casale Monferrato (Northern Italy) were spatially subset in the range lat 45.1 -> 45.2 and long 7.8 -> 8.6 from within SNAP (Fig. 1).

 Fig. 1. Map of the studied area in the central Piedmont region (Northern Italy). Basemap: Google Earth. Created with QGIS.

The NDWI was calculated with the command line as listed above:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))" 

From the analysis of the histogram of the 2017 image, a value of 10 was chosen of a threshold for masking the non-water pixels:

gdal_calc.py -A ndwi.tif --outfile=ndwi_masked.tif --calc="A>10"

The result is represented in Fig. 2A-C, created with QGIS. The result quality is more than acceptable (for the 2017 image the spotted fields to the North of the Po river could possibly correspond to rice fields, but it has to be better investigated).

A - 2017-05-27

B - 2021-03-30

C- 2023-04-04

 Fig. 2. NDWI > 10 for 2017 (a), 2021 (b) and 2023 (c). Basemap: Google Earth. Created with QGIS.

By comparing the results at different times, it is possible to study the temporal evolution of the rivers (in this case, the Po river, Orco torrent and Canale Cavour near the city of Chivasso (Fig. 3).

 Fig. 2. NDWI > 10 for 2017 (red), 2021 (green) and 2023 (blue) in correspondence of the city of Chivasso. The water bodies are the Po river, the Orco torrent and the artificial Canale Cavour. Basemap: Google Earth. Created with QGIS.

 
If we want to automate the task, when we have for instance a lot of input images or we have to perform routinely the same task on new datasets, the run command can be created on-the-fly and executed from within a script created with languages such as Go, Rust, Ocaml and so on, perhaps even concurrently.