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.

An example of atmospheric effect in differential interferograms in the Piedmont Po Plain area (Vercelli-Trino, Northern Italy)?

With differential interferometry from Interferometric Synthetic Aperture Radar (InSAR) it is possible to monitor subsidence or uplift along line-of-sight. One popular satellite constellation for interferometry is Sentinel-1, operated by ESA. It is composed by two satellites, Sentinel-1A and Sentinel-1B, carrying a C band SAR sensor.

Sentinel data are freely available from the Copernicus Open Access Hub (https://scihub.copernicus.eu/).

I downloaded a few Sentinel-1 images of the Piedmont Po Plain (Vercelli-Trino Vercellese and surroundings) to try deriving some interferometric results to interpret.

Fig. 1: map of the studied are in the Piedmont Po Plain (Northern Italy). Google satellite basemap. Created with QGIS.

To see the temporal evolution of the results, a set of four images for the following 2020 periods were used:

• 05/01/2020: S1B_IW_SLC__1SDV_20200501T171426_20200501T171453_021391_0289B0_577E
  
• 05/13/2020: S1B_IW_SLC__1SDV_20200513T171427_20200513T171454_021566_028F1C_D01D
  
• 05/25/2020: S1B_IW_SLC__1SDV_20200525T171428_20200525T171455_021741_029439_E72E
  
• 06/06/2020: S1B_IW_SLC__1SDV_20200606T171428_20200606T171455_021916_02997C_907D
  

Three differential interferometric results were obtained, using SNAP software by ESA and following the indications presented in tutorials and videos published by ESA. The calculated LOS displacements were masked to hide areas were results have low coherence.

The displacements are presented in the following Fig. 2A-C, where Fig. 2A refers to the 05/01->13 time interval, Fig. 2B to the 05/13 -> 25 interval and  Fig. 2C refers to the final 05/25->06/06 time interval.

 (A)

(B)

(C)

Fig. 2: masked displacements along Line-Of-Sight for the three time intervals.  A) 05/01->13 time interval; B) 05/13->25 time interval; C) 05/25->06/06 time interval. The central E-W line is the profile trace used to derive Fig. 3. Google satellite basemap. Created with QGIS.


We can see that in the first time interval (05/01->13) the regional displacements are mainly negative, in the second (05/13->25) positive with a clear E-W gradient, while the situation is more nuanced in the third one (05/13->06/06), without a clear prevalence of positive or negative values.

A profile of the three displacement rasters, along a E-W direction (red line in Fig. 2), was created using the 'gst' Python module (https://gitlab.com/alberese/gst), as described in the 'Trino' Jupyter notebook (https://gitlab.com/alberese/gst/-/blob/master/docs/interferometry/Trino.ipynb).

Fig. 3: Profiles of the masked displacement rasters along the trace in Figs 2A-C. The first period (I) is in red, the second (II) in green and the third (III) in blue. Created with Python modules matplotlib and gst.

The profiles depict the same situation as visible in map, and evidence, near to rightmost profile end, in correspondence of Vigevano town, a local maximum for the first period while the second period presents a local minimum and in the third period there is a flat (Fig. 3).  

It is difficult to consider these variations, both regional and local (Figs. 2 and 3) as related to topographic surface movements, since no constant trend is observed in the three periods, but the displacements during the second period somewhat reverse those of the first periods. The observed results are therefore interpreted as the results of atmospheric effects, particularly strong in the 05/13 image, that impacted with opposite results the displacements of the first (05/01->13) and second period (05/13->25). The results of the third period are more similar to those expected for a no-movement differential interferogram and possibly reflect only local atmospheric effects for both the 05/25/2020 and 06/06/2020 images.