TD23: Creating indexes from a Sentinel scene.¶

Vincent GODARD - V1- 31/01/2024

GIS and Remote Sensing - PhD curriculum¶

Debre Berhan University and University of Paris 8¶

Freely inspired by :

  • Maurício Cordeiro, 2021, Python for Geosciences: Satellite Image Analysis (step by step): https://medium.com/analytics-vidhya/python-for-geosciences-satellite-image-analysis-step-by-step-6d49b1ad567

  • Henrikki Tenkanen & Vuokko Heikinheimo, 2021, Raster map algebra - Calculating NDVI: https://autogis-site.readthedocs.io/en/2019/notebooks/Raster/raster-map-algebra.html

Sources :

  • Data ESA (https://www.esa.int/ESA_Multimedia/Images/2017/03/Sentinel-2_global_coverage);

  • The Level-2A product (https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/product-types/level-2a);

  • Sentinel-2 Products Specification Document (https://sentinel.esa.int/documents/247904/685211/sentinel-2-products-specification-document);

  • The Copernicus Browser (https://dataspace.copernicus.eu/browser).

Download required documents :

Compressed file to download => here.


The dataset was downloaded in a previous tutorial (tdpy21semdoc).

For this tutorial, a window has been cut over the north of Addis Ababa. The data are from 15/12/2023 (T10SFG_20231215T185801_B0?_10m.jp2).

To make this tutorial independent, these first two chapters are identical to those in tutorial "tdpy22semdoc".

1. Load libraries and data¶

TD21 continuation¶

So far we have only downloaded the Sentinel-2 satellite bands for our AOI point.

If you haven't done so, return to TD21.

1.1. Loading data¶

Upload the following file from your storage area to the data directory (TD23/data) under this notebook is located (.ipynb):

○ "T10SFG_20231215T185801_B02_10m.jp2";

○ "T10SFG_20231215T185801_B03_10m.jp2";

○ "T10SFG_20231215T185801_B04_10m.jp2";

○ "T10SFG_20231215T185801_B08_10m.jp2".

This noteboock and the "data" directory will be at the same level in the tree, both under TD23.

1.2. Loading Libraries¶

In [ ]:
## Loading Libraries

# To install Skimage, delete the # (before the "!").
# !pip install rasterio

# !pip install numpy
# !pip install matplotlib
In [ ]:
## Importing libraries
 
import numpy as np
import matplotlib.pyplot as plt
import rasterio as io
from pathlib import Path

2. Download the images¶

Select four images (B02_10m, B03_10m, B04_10m and B08_10m). The size of the pixel is 10 by 10 m.

This four images must be read to produce indexes. Indexes amplify the differences between surface conditions, making it easier to read the landscape.

2.1. Open images in a long but simple way¶

In [ ]:
## Opening the 4 files one by one

B02 = io.open('data/T10SFG_20231215T185801_B02_10m.jp2')
B03 = io.open('data/T10SFG_20231215T185801_B03_10m.jp2')
B04 = io.open('data/T10SFG_20231215T185801_B04_10m.jp2')
B08 = io.open('data/T10SFG_20231215T185801_B08_10m.jp2')

# To test the opening
B08

2.2 Automate the opening with a fonction¶

A function is a piece of code that receives arguments as input, performs a calculation, and then generates output. It can be used several times in the same script or called in other scripts.

In [ ]:
# A few comments :
# Creation of a function which contains a loop for reading certain files (which contain an ident and end with jp2!)
# As Sentinel images contain "T37PDL_.jp2", only they will be selected!
# The letter "f" at the beginning of the string "next(path.glob(f'T10SFG_*_{band}_10m.jp2'))" indicates that it is a formatted string
# this allows Python commands to be included in parentheses to: display a result.

def load_sentinel_image(img_folder, bands): # function name
    image = {}
    path = Path(img_folder)
    for band in bands: # start of loop            
        file = next(path.glob(f'T10SFG_*_{band}_10m.jp2')) # here the formatted string
        print(f'Opening file {file}')
        ds = io.open(file) # creating a dataset with rasterio.open()
        image.update({band: ds.read(1)})

    return image
In [ ]:
# A Python dictionary is an object that can store pairs of keys and values.
# Keys can match any "internal" Python data type (i.e. string, numbers, or tuples)
# Values can be of any arbitrary type/object (here an image).
# Using the previous function to complete the dictionary of "keys" with the loaded bands of our image

img = load_sentinel_image('data', ['B02', 'B03', 'B04', 'B08'])

img.keys()

Just for information

In [ ]:
## Querying the row x column size of our images

img['B02'].shape, img['B03'].shape, img['B04'].shape, img['B08'].shape
In [ ]:
## Description of our variables

type(img), type(img['B04'])
In [ ]:
## Control by displaying of our variables

fig, ax = plt.subplots(1, 4, figsize=(20, 5)) # subplot => to display side by side
ax[0].imshow(img['B02'], cmap='Blues_r') # Append _r to the name of any built-in colormap to get the reversed version;
ax[1].imshow(img['B03'], cmap='Greens_r')
ax[2].imshow(img['B04'], cmap='Reds_r')
ax[3].imshow(img['B08'], cmap='gray_r')

Our four images are plotted with the primary colors in matplotlib library (see https://matplotlib.org/2.0.2/api/colors_api.html)

3. Creation of normalized difference indexes¶

For example, the objective is to find an index to separate vegetated spaces from mineralized spaces or submerged land from emerged land!

3.1. Elimination of zeros and NaN¶

In [ ]:
# As there is a risk of division by "0" if the denominator pixel has this value:
# The first two lines of code check this condition!
# "==" checks if the pixel value is equal to "0"
# "&" means "AND", which must be strictly equivalent
# Similarly, if the same pixel is equal to "0" in both bands then it is not a number (NaN, not a number).

def normalized_difference(img, B1, B2, eps=0.0001):
    Band1 = np.where((img[B1]==0) & (img[B1]==0), np.nan, img[B1])
    Band2 = np.where((img[B1]==0) & (img[B2]==0), np.nan, img[B2])
    
    return (Band1 - Band2) / (Band1 + Band2)

3.2. Creation of the standardized vegetation index (IVN)¶

In [ ]:
## Calculation of the normalized vegetation index (NVI)

ndvi = normalized_difference(img, 'B08', 'B04')
In [ ]:
## Display index

# magic function that renders the figure in a notebook (https://www.data-blogger.com/2017/11/15/python-matplotlib-pyplot-a-perfect-combination/)
# No comment after '%matplotlib inline'
%matplotlib inline

# Plot the NDVI

fig = plt.figure(figsize=(10, 10)) # Figure size in inches (default). Change values if necessary (columns, lines).
plt.imshow(ndvi, cmap='terrain')

# Add colorbar to show the index

plt.colorbar(label='Index', shrink = 0.8)

3.3. Creating the Modified Normalized Difference Water Index (MNDWI)¶

In [ ]:
## Calculation of MNDVI (Modified Normalized Difference Water Index)

mndwi = normalized_difference(img, 'B03', 'B08')
In [ ]:
## Display index

# magic function that renders the figure in a notebook (https://www.data-blogger.com/2017/11/15/python-matplotlib-pyplot-a-perfect-combination/)
# No comment after '%matplotlib inline'
%matplotlib inline

# Plot the MNDVI

fig = plt.figure(figsize=(10, 10)) # Figure size in inches (default). Change values if necessary (columns, lines).
plt.imshow(mndwi, cmap='gray_r')

# Add colorbar to show the index

plt.colorbar(label='Index', shrink = 0.8)

How could you interpret the MNDWI?

3.4. Displaying index images¶

In [ ]:
## Image control

fig, ax = plt.subplots(1, 2, figsize=(18, 9))
ax[0].imshow(ndvi, cmap='gray_r')
ax[1].imshow(mndwi, cmap='gray_r')

Why is it usefull to use the same colormap for the both images?

4. Another index?¶

4.1. Would you like to test your favorite index?¶

In [ ]: