Vincent GODARD - V1- 31/01/2024
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).
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".
So far we have only downloaded the Sentinel-2 satellite bands for our AOI point.
If you haven't done so, return to TD21.
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.
## Loading Libraries
# To install Skimage, delete the # (before the "!").
# !pip install rasterio
# !pip install numpy
# !pip install matplotlib
## Importing libraries
import numpy as np
import matplotlib.pyplot as plt
import rasterio as io
from pathlib import Path
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.
## 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
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.
# 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
# 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
## Querying the row x column size of our images
img['B02'].shape, img['B03'].shape, img['B04'].shape, img['B08'].shape
## Description of our variables
type(img), type(img['B04'])
## 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)
For example, the objective is to find an index to separate vegetated spaces from mineralized spaces or submerged land from emerged land!
# 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)
## Calculation of the normalized vegetation index (NVI)
ndvi = normalized_difference(img, 'B08', 'B04')
## 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)
## Calculation of MNDVI (Modified Normalized Difference Water Index)
mndwi = normalized_difference(img, 'B03', 'B08')
## 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?
## 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?