undo

9adf863c · Dries De Bièvre · ea1c9e5e · 9adf863c
--- a/docs/source/sentinel_2/s2_prepro_2025.ipynb
+++ b/docs/source/sentinel_2/s2_prepro_2025.ipynb
@@ -23,7 +23,7 @@
    "</div>\n",
    "\n",
    "<figure class=\"image\">\n",
-    "  <img src=\"../../../imgs/prepro_S2.png\" alt=\"S2 preprocessing\" width=\"1100\">\n",
+    "  <img src=\"prepro_S2.png\" alt=\"S2 preprocessing\" width=\"1100\">\n",
    "</figure>"
   ]
  },
@@ -118,7 +118,7 @@
    "Select reflectance bands you want to work with\n",
    "\n",
    "<figure class=\"image\">\n",
-    "  <img src=\"../../../imgs/s2_bands.png\" alt=\"S2 bands\" width=\"500\">\n",
+    "  <img src=\"s2_bands.png\" alt=\"S2 bands\" width=\"500\">\n",
    "</figure>"
   ]
  },
@@ -780,7 +780,7 @@
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
-   "version": "3.11.10"
+   "version": "3.11.4"
  },
  "metadata": {
   "interpreter": {

 %% Cell type:markdown id: tags:


 <div class="alert alert-success">
 <h4 class="alert-heading"><strong> Preprocessing Sentinel-2 images
 </strong></h4>
 <hr>

 ### **OBJECTIVES**

 - Resample 20m spatial resolutions bands to 10m (if you want to work with 20m bands)
 - Cropping images to the extent of Region of Interest (ROI)
 - Apply Scene Classification map (SCL) on reflectance images to mask invalid pixels
 - Calculate the fraction of valid pixels in your ROI

 If your ROI spans several S2 tiles, start by preprocessing each tile one by one using this notebook and then use the "Mosaic Sentinel-2 tiles" notebook to mosaic the entire ROI.

 </div>

 <figure class="image">
-  <img src="../../../imgs/prepro_S2.png" alt="S2 preprocessing" width="1100">
+  <img src="prepro_S2.png" alt="S2 preprocessing" width="1100">
 </figure>

 %% Cell type:code id: tags:

 ``` python
 import glob, os
 import numpy as np
 import geopandas as gpd
 import rasterio
 import rasterio.mask
 import rasterio.plot
 from rasterio.enums import Resampling
 import matplotlib.pyplot as plt
 from pathlib import Path

 print('All libraries successfully imported!')
 print(f'Rasterio : {rasterio.__version__}')
 ```

 %% Cell type:markdown id: tags:

 **Set directory**

 %% Cell type:code id: tags:

 ``` python
 computer_path = 'U:/'
 grp_nb        = '30'

 work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/'  # Directory for all your work files

 roi_path = f'{work_path}TP/ROI/'
 dwl_path = f'{work_path}DOWNLOAD/'


 # For each step of the preprocessing,
 # a folder will be created to store
 # the intermediary files.

 resampled_path = f'{work_path}1_L2A_RESAMPLED/'
 clipped_path   = f'{work_path}2_L2A_CLIPPED/'
 masked_path    = f'{work_path}3_L2A_MASKED/'

 Path(resampled_path).mkdir(parents=True, exist_ok=True)
 Path(clipped_path).mkdir(parents=True, exist_ok=True)
 Path(masked_path).mkdir(parents=True, exist_ok=True)
 ```

 %% Cell type:markdown id: tags:

 **Set region of interest**

 %% Cell type:code id: tags:

 ``` python
 ## change these variables to your own crs and roi name
 crs = 'EPSG:32631'
 roi_name = 'roi'


 roi_filename = f'{roi_name}_{crs[5:]}.shp'
 roi_file     = f'{roi_path}{roi_filename}'

 roi_gdf = gpd.read_file(roi_file)
 display(roi_gdf)
 # print(f'ROI shapefile : {roi_file}')
 ```

 %% Cell type:markdown id: tags:

 Select reflectance bands you want to work with

 <figure class="image">
-  <img src="../../../imgs/s2_bands.png" alt="S2 bands" width="500">
+  <img src="s2_bands.png" alt="S2 bands" width="500">
 </figure>

 %% Cell type:code id: tags:

 ``` python
 band_10m_list = ['B04','B08']

 band_20m_list = ['B11']
 ```

 %% Cell type:markdown id: tags:

 **Set parameters for resampling**

 %% Cell type:code id: tags:

 ``` python
 # We will upscale the image by 2: 20m --> 10m
 upscale_factor = 2

 # As SCL is categorical data, we MUST use "nearest neighbor" resampling method
 resampling_method_categorical = Resampling.nearest

 # As BOA is continuous data, we can use other resampling methods : nearest, bilinear, cubic
 resampling_method_continuous = Resampling.bilinear

 nodata_val = -10000
 ```

 %% Cell type:markdown id: tags:


 <div class="alert alert-warning">

 <h4 class="alert-heading"><strong>First exercise of the session</strong></h4>
 <hr>

 <p class="mb-0">List all your images using the glob library </p>

 </div>


 #

 %% Cell type:markdown id: tags:

 ## Using `glob.glob` in Python

 The `glob` module helps find file paths matching a specified pattern. It is useful for listing files of a specific type or format.

 ### Common Wildcards:
 - `*` → Matches any number of characters (`*.txt` finds all `.txt` files).
 - `?` → Matches a single character (`file?.csv` matches `file1.csv`, `file2.csv`, etc.).
 - `[abc]` → Matches any character inside the brackets (`file[12].txt` matches `file1.txt` and `file2.txt`).

 ### Example:
 ```python
 import glob

 # Get all `.txt` files in the current directory
 txt_files = glob.glob("*.txt")

 # Get all '.tif' files in the C:/Download/ directory

 tif_files = glob.glob("C:/Download/*.tif")

 %% Cell type:markdown id: tags:

 <div class="alert alert-warning" >

 ### Exercise: Find `.SAFE` Files in Your Downloads Folder

 Your task is to use the `glob` module to create a list of all `.SAFE` files in your **Downloads** folder.

 The list should be stored in the `list_L2A`  variable

 #### Steps:
 1. Use `glob.glob()` with the correct path pattern to find `.SAFE` files in your download path.
 2. Print the resulting list of images.

 **Hint 1 :** your images are located in "X:/STUDENTS/GROUP_{your group number}/DOWNLOAD/", Adjust accordingly.

  **Hint 2 :** We stored the path to your downloaded images in a variable a the top of the notebook

 </div>

 %% Cell type:code id: tags:

 ``` python
 list_L2A =

 print(f'{len(list_L2A)} L2A will be pre-processed \n')

 for L2A_safe in list_L2A:
    L2A_name = os.path.basename(L2A_safe)
    print(L2A_name)
 ```

 %% Cell type:markdown id: tags:

 ## 1. Resample images at 20m resolution to 10m

 | Methods | Type of data | How it works |
 |:---------:|:----------:| ---- |
 | Nearest Neighbor | categorical| <font size="2">The value of of the output cell is determined by the nearest cell center on the input grid </font>|
 | Bilinear Interpolation | continuous | <font size="2"> Weighted average of the four nearest cell centers. <br/> The closer an input cell center is to the output cell center, the higher the influence of its value is on the output cell value. The output value could be different than the nearest input but is always within the same range of values as the input. </font>|
 | Cubic Convolution | continuous | <font size="2">Looks at the 16 nearest cell centers to the output and fits a smooth curve through the points to find the value. <br/>Not only does this change the values of the input but it could also cause the output value to be outside of the range of input values (imagine a sink or a peak occurring on a surface). </font> |


 ### 1.1 Resample Scene Classification map

 > **WARNING**: Only if you are planning to work at 10m resolution. If you are planning to work at 20m resolution, you can skip this step!

 %% Cell type:markdown id: tags:

 <div class = 'alert alert-success'>

 ### Function: `resample_SCL(list_L2A, resampled_path, upscale_factor)`

 This function **resamples the Scene Classification Layer (SCL)** from 20m to 10m resolution.

 #### **Arguments:**
 1. **`list_L2A`** *(list of str)* → A list of Sentinel-2 L2A `.SAFE` directories.
 2. **`resampled_path`** *(str)* → The directory where the resampled SCL images will be saved.
 3. **`upscale_factor`** *(int or float)* → The factor by which the image resolution will be increased (e.g., `2` for 20m → 10m).

 #### **Process:**
 - Finds the `*_SCL_20m.jp2` file in each L2A folder.
 - Reads the image using `rasterio`.
 - Resamples it using **nearest neighbor interpolation** (to preserve categorical values).
 - Writes the new **GeoTIFF (.tif)** image at 10m resolution.
 - Saves it to `resampled_path`.

 </div>

 %% Cell type:code id: tags:

 ``` python
 from s2_func_prepro import resample_SCL

 resample_SCL(list_L2A, resampled_path, upscale_factor)
 ```

 %% Cell type:markdown id: tags:

 ### 1.2 Resample Bottom-Of-Atmosphere corrected reflectance

 If you work only with Blue (B02), Green (B03), Red (B04) and NIR (B08) bands, you don't have to do the resampling step because theses reflectances are already available at 10m resolution.

 This step is only necessary if you work with bands 5,6,7,8A,11,12 which are only available at 20m resolution.

 %% Cell type:markdown id: tags:

 <div class="alert alert-success">

 ### Function: `resample_bands(list_L2A, resampled_path, upscale_factor, bands_20m_list, resampling_method=Resampling.bilinear)`

 This function **resamples selected Sentinel-2 bands** from 20m to 10m resolution.

 #### **Arguments:**
 1. **`list_L2A`** *(list of str)* → A list of Sentinel-2 L2A `.SAFE` directories.
 2. **`resampled_path`** *(str)* → The directory where the resampled band images will be saved.
 3. **`upscale_factor`** *(int or float)* → The factor by which the image resolution will be increased.
 4. **`bands_20m_list`** *(list of str)* → List of band names to be resampled (e.g., `['B05', 'B06', 'B07']`).
 5. **`resampling_method`** *(Resampling method, optional)* → Defines how pixel values are interpolated during resampling. Default: **Bilinear interpolation**.

 #### **Process:**
 - Iterates through the specified bands in `bands_20m_list`.
 - Finds the corresponding `*_20m.jp2` file.
 - Reads the image using `rasterio`.
 - Resamples it using the chosen interpolation method (**bilinear by default** for continuous reflectance values).
 - Writes the new **GeoTIFF (.tif)** image at 10m resolution.
 - Saves it to `resampled_path`.

 </div>

 %% Cell type:code id: tags:

 ``` python
 from s2_func_prepro import resample_bands

 resample_bands(list_L2A, resampled_path, upscale_factor, band_20m_list, resampling_method_continuous)
 ```

 %% Cell type:markdown id: tags:

 ## 2. Clip images to the extent of Region of Interest

 We can use our ROI (vector) to clip satellite images (raster) into a smaller area to reduce image storage and speed up further processing.

 %% Cell type:markdown id: tags:

 Get list of all images to clip

 > **WARNING**: All images must be at same spatial resolution !

 %% Cell type:markdown id: tags:

 <div class="alert alert-warning">

 <h4 class="alert-heading"><strong>Second exercise of the session</strong></h4>
 <hr>

 ### Exercise: Find All 10m Resolution Images for Clipping

 Your task is to use the `glob` module to create a list of **all 10m resolution images** from both the **resampled folder** and the **original downloaded SAFE files**.

 The list should be stored in the `list_im_to_clip` variable.

 #### Steps:
 1. Use `glob.glob()` to find all `*_10m.tif` images in your **resampled folder**.
 2. Loop through each band in `band_10m_list` and use `glob.glob()` to find `*_10m.jp2` images in the **original `.SAFE` directories**.
 3. Append the found images to `list_im_to_clip`.
 4. Print the total number of images found.

 **Hint 1:**
 - Your **resampled images** are in `"X:/STUDENTS/GROUP_{your group number}/RESAMPLED/"`.
 - Your **original images** are in `"X:/STUDENTS/GROUP_{your group number}/DOWNLOAD/*.SAFE/GRANULE/*/IMG_DATA/R10m/"`.

 **Hint 2:**
 - The path to your **resampled images** is stored in `resampled_path`.
 - The path to your **downloaded images** is stored in `dwl_path`.
 - You have to loop through the spectral bands and find all the jp2 files corresponding to these spectral bands

 </div>

 %% Cell type:code id: tags:

 ``` python
 list_im_to_clip =

 for band in band_10m_list:

    list_im_to_clip += glob.glob()

 print(f'')
 ```

 %% Cell type:markdown id: tags:

 **Check if the ROI is located in one of the image to clip**

 %% Cell type:code id: tags:

 ``` python
 fig, ax = plt.subplots(1, 1, figsize=(10,10))

 # Plot image

 # First image of the list
 im_file = list_im_to_clip[0]

 # A cloudy image
 #im_file = [s for s in list_im_to_clip if all(xs in s for xs in ['20200621','B04'])][0]

 src = rasterio.open(im_file, "r")
 im = src.read(1)

 im = im.astype(float)
 im = np.where(im == nodata_val, np.nan, im)
 p5 = np.nanpercentile(im, 5)
 p95 = np.nanpercentile(im, 95)
 rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax)

 # Plot vector
 roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 4, ax=ax)

 plt.box(False)
 ```

 %% Cell type:markdown id: tags:

 ## **Clip your images**

 %% Cell type:markdown id: tags:

 <div class="alert alert-success">

 ### Function: `clip_imgs(list_im_to_clip, clipped_path, roi_gdf)`

 This function **clips a list of raster images** to the extent of a given region of interest (ROI).

 #### **Arguments:**
 1. **`list_im_to_clip`** *(list of str)* → A list of file paths to the images that need to be clipped.
 2. **`clipped_path`** *(str)* → The directory where the clipped images will be saved.
 3. **`roi_gdf`** *(GeoDataFrame)* → A GeoDataFrame containing the region of interest (ROI) geometry used for clipping.

 #### **Process:**
 - Iterates through each image file in `list_im_to_clip`.
 - Generates the output filename by appending `_ROI.tif` to the original filename.
 - Opens the raster file using `rasterio`.
 - Clips the raster to the extent of the `roi_gdf` geometry using `rasterio.mask.mask()`.
 - Updates the metadata to match the new clipped image dimensions and transformation.
 - Writes the clipped image to a new **GeoTIFF (.tif)** file.
 - Saves the clipped image to `clipped_path`.

 </div>

 %% Cell type:markdown id: tags:

 <div class="alert alert-warning">

 <h4 class="alert-heading"><strong>Third exercise of the session</strong></h4>
 <hr>

 ### Exercise: Clip your images using the clip_imgs function

 Your task is to use the `clip_imgs` function to clip all **10m resolution images** in `list_im_to_clip` to the extent of the provided **region of interest (ROI)**.


 #### Steps:

 1. Verify that all variables have already been defined above (your ROI GeoDataFrame, your list of images to clip, and the output path)
 2. Run the functions with the correct arguments


 </div>

 %% Cell type:code id: tags:

 ``` python
 from s2_func_prepro import clip_imgs

 clip_imgs()
 ```

 %% Cell type:markdown id: tags:

 **Check if clipped image is located inside the ROI**

 %% Cell type:code id: tags:

 ``` python
 fig, ax = plt.subplots(1, 1, figsize=(10,10))

 # Plot image

 # First image of the list
 im_file = glob.glob(f'{clipped_path}*_B*.tif')[0]
 print(im_file)
 src = rasterio.open(im_file, "r")
 im = src.read(1)

 im = im.astype(float)
 im = np.where(im == nodata_val, np.nan, im)
 p5 = np.nanpercentile(im, 5)
 p95 = np.nanpercentile(im, 95)

 rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax)

 # Plot vector
 roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 6, ax=ax)

 plt.box(False)
 ```

 %% Cell type:markdown id: tags:

 ## 3. Atmospheric correction & Cloud screening

 ### 3.1 Atmospheric correction

 Atmospheric correction allows calculation of Bottom Of Atmosphere (BOA) reflectance from Top Of Atmosphere (TOA) reflectance images available in L1C products. For Sentinel-2 images, we will directly use the L2A data as available from the ground segment.

 ### 3.2 Cloud screening

 Reliable identification of clouds and cloud shadows are necessary for any optical remote sensing image analysis, especially in operational and fully automatic setups. The cloud screening can be achieved using different algorithms (Sen2Cor, MAJA, Fmask). In this course we will work with Sen2Cor as the mask is already present in L2A products.

 Scene Classification map (SCL) aims at providing a pixel classification map (cloud, cloud shadows, vegetation, soils/deserts, water, snow, etc.)

 The SC algorithm enables:
 - generation of a classification map which includes four different classes for clouds (including cirrus) and six different classifications for shadows, cloud shadows, vegetation, soils/deserts, water and snow
 - provision of associated quality indicators corresponding to a map of cloud probability and a map of snow probability.

 SCL class| Description
 :---------:|:----------:
 0| No data
 1| Saturated or defective
 2|	Dark area pixels
 3|	Cloud shadows
 4|	Vegetation
 5|	Not vegetated
 6|	Water
 7|	Unclassified
 8|	Cloud medium probability
 9|	Cloud high probability
 10|	Thin cirrus
 11|	Snow

 %% Cell type:markdown id: tags:

 Get a list with all reflectance bands clipped on the ROI

 %% Cell type:code id: tags:

 ``` python
 list_im_ROI = glob.glob(f'{clipped_path}*_B*_ROI.tif')

 print(f'There are {len(list_im_ROI)} images where we have to apply SCL map')
 ```

 %% Cell type:markdown id: tags:

 Apply SCL

 %% Cell type:code id: tags:

 ``` python
 from s2_func_prepro import apply_SCL

 apply_SCL(list_im_ROI, clipped_path, masked_path, nodata_val)
 ```

 %% Cell type:markdown id: tags:

 ## **Check if the SCL was correctly applied**

 %% Cell type:code id: tags:

 ``` python
 date_to_display = '20250218'
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20,10))


 # Plot reflectance image before applying SCL
 # ------------------------------------------

 im_file_reflect = glob.glob(f'{clipped_path}*{date_to_display}*B04*.tif')[0]
 #im_file_reflect = glob.glob(f'{clipped_path}*31UFR*B04*.tif')[0]


 src = rasterio.open(im_file_reflect, "r")
 im = src.read(1)

 im = im.astype(float)
 im = np.where(im == nodata_val, np.nan, im)
 p5 = np.nanpercentile(im, 5)
 p95 = np.nanpercentile(im, 95)

 rasterio.plot.show(src, ax=ax1, cmap='Greys_r', vmin=p5, vmax=p95, title="Before SCL")


 # Plot reflectance image after applying SCL
 # -----------------------------------------

 im_file_scl = glob.glob(f'{masked_path}*{date_to_display}*B04*.tif')[0]
 #im_file_scl = glob.glob(f'{masked_path}*31UFR*B04*.tif')[0]


 src = rasterio.open(im_file_scl, "r")

 color_map = plt.cm.get_cmap("Greys_r")
 color_map.set_bad(color='red')

 rasterio.plot.show(src, ax=ax2, cmap=color_map, vmin=p5, vmax=p95, title="After SCL")
 ```

 %% Cell type:markdown id: tags:

 <div class="alert alert-warning">

 <h4 class="alert-heading"><strong>Fourth exercise of the session</strong></h4>
 <hr>

 ### Exercise: Sum Cloud Pixels and Plot the Results

 Your task is to sum all cloud pixels for each image based on the SCL cloud mask and plot the total number of cloud pixels against the **image dates**.

 #### Steps:
 1. Use `glob.glob` to list all the SCL masks paths
 2. Iterate through the list with a `for` loop
 3. In the loop, open each image with Rasterio (explained below) and read it as a Numpy array
 4. Store the **sum of cloud pixels** and the corresponding **image date**.
 5. Plot the **cloud pixel count** on the y-axis and **image dates** on the x-axis.

 #### **Hints:**
 - Use numpy to sum the number of pixels with categorical values equal to clouds
 - You can extract the **image date** from the image filename or metadata. For example, the filename might contain the date in `YYYYMMDD` format
 - For plotting, use `matplotlib` or another plotting library

 SCL class| Description
 :---------:|:----------:
 0| No data
 1| Saturated or defective
 2|	Dark area pixels
 3|	Cloud shadows
 4|	Vegetation
 5|	Not vegetated
 6|	Water
 7|	Unclassified
 8|	Cloud medium probability
 9|	Cloud high probability
 10|	Thin cirrus
 11|	Snow

 </div>

 %% Cell type:markdown id: tags:

 <div class="alert alert-info">

 ### How to Open an Image with `rasterio` and Convert it to a NumPy Array

 To open a raster image and convert it to a NumPy array, you can use the following steps:

 1. **Open the raster file** using `rasterio.open()`.
 2. **Read the image data** into a NumPy array using the `read()` method.

 Here’s an example:

 ```python
 import rasterio

 # Open the image file
 with rasterio.open('path_to_your_image.tif') as src:
    # Convert the image to a NumPy array
    image_array = src.read(1)  # '1' refers to the first band (change if needed)

 # Now 'image_array' contains the raster data as a NumPy array