Skip to content
Extraits de code Groupes Projets
Valider ea1c9e5e rédigé par Dries De Bièvre's avatar Dries De Bièvre
Parcourir les fichiers

correct imgs

parent 1e21e85e
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-success"> <div class="alert alert-success">
<h4 class="alert-heading"><strong> Preprocessing Sentinel-2 images <h4 class="alert-heading"><strong> Preprocessing Sentinel-2 images
</strong></h4> </strong></h4>
<hr> <hr>
### **OBJECTIVES** ### **OBJECTIVES**
- Resample 20m spatial resolutions bands to 10m (if you want to work with 20m bands) - 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) - Cropping images to the extent of Region of Interest (ROI)
- Apply Scene Classification map (SCL) on reflectance images to mask invalid pixels - Apply Scene Classification map (SCL) on reflectance images to mask invalid pixels
- Calculate the fraction of valid pixels in your ROI - 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. 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> </div>
<figure class="image"> <figure class="image">
<img src="prepro_S2.png" alt="S2 preprocessing" width="1100"> <img src="../../../imgs/prepro_S2.png" alt="S2 preprocessing" width="1100">
</figure> </figure>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import glob, os import glob, os
import numpy as np import numpy as np
import geopandas as gpd import geopandas as gpd
import rasterio import rasterio
import rasterio.mask import rasterio.mask
import rasterio.plot import rasterio.plot
from rasterio.enums import Resampling from rasterio.enums import Resampling
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from pathlib import Path from pathlib import Path
print('All libraries successfully imported!') print('All libraries successfully imported!')
print(f'Rasterio : {rasterio.__version__}') print(f'Rasterio : {rasterio.__version__}')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Set directory** **Set directory**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
computer_path = 'U:/' computer_path = 'U:/'
grp_nb = '30' grp_nb = '30'
work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/' # Directory for all your work files work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/' # Directory for all your work files
roi_path = f'{work_path}TP/ROI/' roi_path = f'{work_path}TP/ROI/'
dwl_path = f'{work_path}DOWNLOAD/' dwl_path = f'{work_path}DOWNLOAD/'
# For each step of the preprocessing, # For each step of the preprocessing,
# a folder will be created to store # a folder will be created to store
# the intermediary files. # the intermediary files.
resampled_path = f'{work_path}1_L2A_RESAMPLED/' resampled_path = f'{work_path}1_L2A_RESAMPLED/'
clipped_path = f'{work_path}2_L2A_CLIPPED/' clipped_path = f'{work_path}2_L2A_CLIPPED/'
masked_path = f'{work_path}3_L2A_MASKED/' masked_path = f'{work_path}3_L2A_MASKED/'
Path(resampled_path).mkdir(parents=True, exist_ok=True) Path(resampled_path).mkdir(parents=True, exist_ok=True)
Path(clipped_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) Path(masked_path).mkdir(parents=True, exist_ok=True)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Set region of interest** **Set region of interest**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
## change these variables to your own crs and roi name ## change these variables to your own crs and roi name
crs = 'EPSG:32631' crs = 'EPSG:32631'
roi_name = 'roi' roi_name = 'roi'
roi_filename = f'{roi_name}_{crs[5:]}.shp' roi_filename = f'{roi_name}_{crs[5:]}.shp'
roi_file = f'{roi_path}{roi_filename}' roi_file = f'{roi_path}{roi_filename}'
roi_gdf = gpd.read_file(roi_file) roi_gdf = gpd.read_file(roi_file)
display(roi_gdf) display(roi_gdf)
# print(f'ROI shapefile : {roi_file}') # print(f'ROI shapefile : {roi_file}')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Select reflectance bands you want to work with Select reflectance bands you want to work with
%% Cell type:markdown id: tags:
<figure class="image"> <figure class="image">
<img src="s2_bands.png" alt="S2 bands" width="500"> <img src="../../../imgs/s2_bands.png" alt="S2 bands" width="500">
</figure> </figure>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
band_10m_list = ['B04','B08'] band_10m_list = ['B04','B08']
band_20m_list = ['B11'] band_20m_list = ['B11']
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Set parameters for resampling** **Set parameters for resampling**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# We will upscale the image by 2: 20m --> 10m # We will upscale the image by 2: 20m --> 10m
upscale_factor = 2 upscale_factor = 2
# As SCL is categorical data, we MUST use "nearest neighbor" resampling method # As SCL is categorical data, we MUST use "nearest neighbor" resampling method
resampling_method_categorical = Resampling.nearest resampling_method_categorical = Resampling.nearest
# As BOA is continuous data, we can use other resampling methods : nearest, bilinear, cubic # As BOA is continuous data, we can use other resampling methods : nearest, bilinear, cubic
resampling_method_continuous = Resampling.bilinear resampling_method_continuous = Resampling.bilinear
nodata_val = -10000 nodata_val = -10000
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-warning"> <div class="alert alert-warning">
<h4 class="alert-heading"><strong>First exercise of the session</strong></h4> <h4 class="alert-heading"><strong>First exercise of the session</strong></h4>
<hr> <hr>
<p class="mb-0">List all your images using the glob library </p> <p class="mb-0">List all your images using the glob library </p>
</div> </div>
# #
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Using `glob.glob` in Python ## 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. 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: ### Common Wildcards:
- `*` → Matches any number of characters (`*.txt` finds all `.txt` files). - `*` → Matches any number of characters (`*.txt` finds all `.txt` files).
- `?` → Matches a single character (`file?.csv` matches `file1.csv`, `file2.csv`, etc.). - `?` → 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`). - `[abc]` → Matches any character inside the brackets (`file[12].txt` matches `file1.txt` and `file2.txt`).
### Example: ### Example:
```python ```python
import glob import glob
# Get all `.txt` files in the current directory # Get all `.txt` files in the current directory
txt_files = glob.glob("*.txt") txt_files = glob.glob("*.txt")
# Get all '.tif' files in the C:/Download/ directory # Get all '.tif' files in the C:/Download/ directory
tif_files = glob.glob("C:/Download/*.tif") tif_files = glob.glob("C:/Download/*.tif")
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-warning" > <div class="alert alert-warning" >
### Exercise: Find `.SAFE` Files in Your Downloads Folder ### 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. 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 The list should be stored in the `list_L2A` variable
#### Steps: #### Steps:
1. Use `glob.glob()` with the correct path pattern to find `.SAFE` files in your download path. 1. Use `glob.glob()` with the correct path pattern to find `.SAFE` files in your download path.
2. Print the resulting list of images. 2. Print the resulting list of images.
**Hint 1 :** your images are located in "X:/STUDENTS/GROUP_{your group number}/DOWNLOAD/", Adjust accordingly. **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 **Hint 2 :** We stored the path to your downloaded images in a variable a the top of the notebook
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
list_L2A = list_L2A =
print(f'{len(list_L2A)} L2A will be pre-processed \n') print(f'{len(list_L2A)} L2A will be pre-processed \n')
for L2A_safe in list_L2A: for L2A_safe in list_L2A:
L2A_name = os.path.basename(L2A_safe) L2A_name = os.path.basename(L2A_safe)
print(L2A_name) print(L2A_name)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## 1. Resample images at 20m resolution to 10m ## 1. Resample images at 20m resolution to 10m
| Methods | Type of data | How it works | | 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>| | 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>| | 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> | | 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 ### 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! > **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: %% Cell type:markdown id: tags:
<div class = 'alert alert-success'> <div class = 'alert alert-success'>
### Function: `resample_SCL(list_L2A, resampled_path, upscale_factor)` ### Function: `resample_SCL(list_L2A, resampled_path, upscale_factor)`
This function **resamples the Scene Classification Layer (SCL)** from 20m to 10m resolution. This function **resamples the Scene Classification Layer (SCL)** from 20m to 10m resolution.
#### **Arguments:** #### **Arguments:**
1. **`list_L2A`** *(list of str)* → A list of Sentinel-2 L2A `.SAFE` directories. 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. 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). 3. **`upscale_factor`** *(int or float)* → The factor by which the image resolution will be increased (e.g., `2` for 20m → 10m).
#### **Process:** #### **Process:**
- Finds the `*_SCL_20m.jp2` file in each L2A folder. - Finds the `*_SCL_20m.jp2` file in each L2A folder.
- Reads the image using `rasterio`. - Reads the image using `rasterio`.
- Resamples it using **nearest neighbor interpolation** (to preserve categorical values). - Resamples it using **nearest neighbor interpolation** (to preserve categorical values).
- Writes the new **GeoTIFF (.tif)** image at 10m resolution. - Writes the new **GeoTIFF (.tif)** image at 10m resolution.
- Saves it to `resampled_path`. - Saves it to `resampled_path`.
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from s2_func_prepro import resample_SCL from s2_func_prepro import resample_SCL
resample_SCL(list_L2A, resampled_path, upscale_factor) resample_SCL(list_L2A, resampled_path, upscale_factor)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### 1.2 Resample Bottom-Of-Atmosphere corrected reflectance ### 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. 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. 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: %% Cell type:markdown id: tags:
<div class="alert alert-success"> <div class="alert alert-success">
### Function: `resample_bands(list_L2A, resampled_path, upscale_factor, bands_20m_list, resampling_method=Resampling.bilinear)` ### 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. This function **resamples selected Sentinel-2 bands** from 20m to 10m resolution.
#### **Arguments:** #### **Arguments:**
1. **`list_L2A`** *(list of str)* → A list of Sentinel-2 L2A `.SAFE` directories. 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. 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. 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']`). 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**. 5. **`resampling_method`** *(Resampling method, optional)* → Defines how pixel values are interpolated during resampling. Default: **Bilinear interpolation**.
#### **Process:** #### **Process:**
- Iterates through the specified bands in `bands_20m_list`. - Iterates through the specified bands in `bands_20m_list`.
- Finds the corresponding `*_20m.jp2` file. - Finds the corresponding `*_20m.jp2` file.
- Reads the image using `rasterio`. - Reads the image using `rasterio`.
- Resamples it using the chosen interpolation method (**bilinear by default** for continuous reflectance values). - Resamples it using the chosen interpolation method (**bilinear by default** for continuous reflectance values).
- Writes the new **GeoTIFF (.tif)** image at 10m resolution. - Writes the new **GeoTIFF (.tif)** image at 10m resolution.
- Saves it to `resampled_path`. - Saves it to `resampled_path`.
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from s2_func_prepro import resample_bands from s2_func_prepro import resample_bands
resample_bands(list_L2A, resampled_path, upscale_factor, band_20m_list, resampling_method_continuous) resample_bands(list_L2A, resampled_path, upscale_factor, band_20m_list, resampling_method_continuous)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## 2. Clip images to the extent of Region of Interest ## 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. 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: %% Cell type:markdown id: tags:
Get list of all images to clip Get list of all images to clip
> **WARNING**: All images must be at same spatial resolution ! > **WARNING**: All images must be at same spatial resolution !
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-warning"> <div class="alert alert-warning">
<h4 class="alert-heading"><strong>Second exercise of the session</strong></h4> <h4 class="alert-heading"><strong>Second exercise of the session</strong></h4>
<hr> <hr>
### Exercise: Find All 10m Resolution Images for Clipping ### 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**. 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. The list should be stored in the `list_im_to_clip` variable.
#### Steps: #### Steps:
1. Use `glob.glob()` to find all `*_10m.tif` images in your **resampled folder**. 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**. 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`. 3. Append the found images to `list_im_to_clip`.
4. Print the total number of images found. 4. Print the total number of images found.
**Hint 1:** **Hint 1:**
- Your **resampled images** are in `"X:/STUDENTS/GROUP_{your group number}/RESAMPLED/"`. - 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/"`. - Your **original images** are in `"X:/STUDENTS/GROUP_{your group number}/DOWNLOAD/*.SAFE/GRANULE/*/IMG_DATA/R10m/"`.
**Hint 2:** **Hint 2:**
- The path to your **resampled images** is stored in `resampled_path`. - The path to your **resampled images** is stored in `resampled_path`.
- The path to your **downloaded images** is stored in `dwl_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 - You have to loop through the spectral bands and find all the jp2 files corresponding to these spectral bands
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
list_im_to_clip = list_im_to_clip =
for band in band_10m_list: for band in band_10m_list:
list_im_to_clip += glob.glob() list_im_to_clip += glob.glob()
print(f'') print(f'')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Check if the ROI is located in one of the image to clip** **Check if the ROI is located in one of the image to clip**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, ax = plt.subplots(1, 1, figsize=(10,10)) fig, ax = plt.subplots(1, 1, figsize=(10,10))
# Plot image # Plot image
# First image of the list # First image of the list
im_file = list_im_to_clip[0] im_file = list_im_to_clip[0]
# A cloudy image # A cloudy image
#im_file = [s for s in list_im_to_clip if all(xs in s for xs in ['20200621','B04'])][0] #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") src = rasterio.open(im_file, "r")
im = src.read(1) im = src.read(1)
im = im.astype(float) im = im.astype(float)
im = np.where(im == nodata_val, np.nan, im) im = np.where(im == nodata_val, np.nan, im)
p5 = np.nanpercentile(im, 5) p5 = np.nanpercentile(im, 5)
p95 = np.nanpercentile(im, 95) p95 = np.nanpercentile(im, 95)
rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax) rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax)
# Plot vector # Plot vector
roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 4, ax=ax) roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 4, ax=ax)
plt.box(False) plt.box(False)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## **Clip your images** ## **Clip your images**
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-success"> <div class="alert alert-success">
### Function: `clip_imgs(list_im_to_clip, clipped_path, roi_gdf)` ### 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). This function **clips a list of raster images** to the extent of a given region of interest (ROI).
#### **Arguments:** #### **Arguments:**
1. **`list_im_to_clip`** *(list of str)* → A list of file paths to the images that need to be clipped. 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. 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. 3. **`roi_gdf`** *(GeoDataFrame)* → A GeoDataFrame containing the region of interest (ROI) geometry used for clipping.
#### **Process:** #### **Process:**
- Iterates through each image file in `list_im_to_clip`. - Iterates through each image file in `list_im_to_clip`.
- Generates the output filename by appending `_ROI.tif` to the original filename. - Generates the output filename by appending `_ROI.tif` to the original filename.
- Opens the raster file using `rasterio`. - Opens the raster file using `rasterio`.
- Clips the raster to the extent of the `roi_gdf` geometry using `rasterio.mask.mask()`. - 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. - Updates the metadata to match the new clipped image dimensions and transformation.
- Writes the clipped image to a new **GeoTIFF (.tif)** file. - Writes the clipped image to a new **GeoTIFF (.tif)** file.
- Saves the clipped image to `clipped_path`. - Saves the clipped image to `clipped_path`.
</div> </div>
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-warning"> <div class="alert alert-warning">
<h4 class="alert-heading"><strong>Third exercise of the session</strong></h4> <h4 class="alert-heading"><strong>Third exercise of the session</strong></h4>
<hr> <hr>
### Exercise: Clip your images using the clip_imgs function ### 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)**. 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: #### Steps:
1. Verify that all variables have already been defined above (your ROI GeoDataFrame, your list of images to clip, and the output path) 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 2. Run the functions with the correct arguments
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from s2_func_prepro import clip_imgs from s2_func_prepro import clip_imgs
clip_imgs() clip_imgs()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Check if clipped image is located inside the ROI** **Check if clipped image is located inside the ROI**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, ax = plt.subplots(1, 1, figsize=(10,10)) fig, ax = plt.subplots(1, 1, figsize=(10,10))
# Plot image # Plot image
# First image of the list # First image of the list
im_file = glob.glob(f'{clipped_path}*_B*.tif')[0] im_file = glob.glob(f'{clipped_path}*_B*.tif')[0]
print(im_file) print(im_file)
src = rasterio.open(im_file, "r") src = rasterio.open(im_file, "r")
im = src.read(1) im = src.read(1)
im = im.astype(float) im = im.astype(float)
im = np.where(im == nodata_val, np.nan, im) im = np.where(im == nodata_val, np.nan, im)
p5 = np.nanpercentile(im, 5) p5 = np.nanpercentile(im, 5)
p95 = np.nanpercentile(im, 95) p95 = np.nanpercentile(im, 95)
rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax) rasterio.plot.show(src, cmap='Greys_r', vmin=p5, vmax=p95, ax=ax)
# Plot vector # Plot vector
roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 6, ax=ax) roi_gdf.plot(facecolor='none', edgecolor='red', linewidth = 6, ax=ax)
plt.box(False) plt.box(False)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## 3. Atmospheric correction & Cloud screening ## 3. Atmospheric correction & Cloud screening
### 3.1 Atmospheric correction ### 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. 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 ### 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. 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.) Scene Classification map (SCL) aims at providing a pixel classification map (cloud, cloud shadows, vegetation, soils/deserts, water, snow, etc.)
The SC algorithm enables: 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 - 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. - provision of associated quality indicators corresponding to a map of cloud probability and a map of snow probability.
SCL class| Description SCL class| Description
:---------:|:----------: :---------:|:----------:
0| No data 0| No data
1| Saturated or defective 1| Saturated or defective
2| Dark area pixels 2| Dark area pixels
3| Cloud shadows 3| Cloud shadows
4| Vegetation 4| Vegetation
5| Not vegetated 5| Not vegetated
6| Water 6| Water
7| Unclassified 7| Unclassified
8| Cloud medium probability 8| Cloud medium probability
9| Cloud high probability 9| Cloud high probability
10| Thin cirrus 10| Thin cirrus
11| Snow 11| Snow
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Get a list with all reflectance bands clipped on the ROI Get a list with all reflectance bands clipped on the ROI
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
list_im_ROI = glob.glob(f'{clipped_path}*_B*_ROI.tif') 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') print(f'There are {len(list_im_ROI)} images where we have to apply SCL map')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Apply SCL Apply SCL
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from s2_func_prepro import apply_SCL from s2_func_prepro import apply_SCL
apply_SCL(list_im_ROI, clipped_path, masked_path, nodata_val) apply_SCL(list_im_ROI, clipped_path, masked_path, nodata_val)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## **Check if the SCL was correctly applied** ## **Check if the SCL was correctly applied**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
date_to_display = '20250218' date_to_display = '20250218'
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20,10)) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20,10))
# Plot reflectance image before applying SCL # 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}*{date_to_display}*B04*.tif')[0]
#im_file_reflect = glob.glob(f'{clipped_path}*31UFR*B04*.tif')[0] #im_file_reflect = glob.glob(f'{clipped_path}*31UFR*B04*.tif')[0]
src = rasterio.open(im_file_reflect, "r") src = rasterio.open(im_file_reflect, "r")
im = src.read(1) im = src.read(1)
im = im.astype(float) im = im.astype(float)
im = np.where(im == nodata_val, np.nan, im) im = np.where(im == nodata_val, np.nan, im)
p5 = np.nanpercentile(im, 5) p5 = np.nanpercentile(im, 5)
p95 = np.nanpercentile(im, 95) p95 = np.nanpercentile(im, 95)
rasterio.plot.show(src, ax=ax1, cmap='Greys_r', vmin=p5, vmax=p95, title="Before SCL") rasterio.plot.show(src, ax=ax1, cmap='Greys_r', vmin=p5, vmax=p95, title="Before SCL")
# Plot reflectance image after applying 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}*{date_to_display}*B04*.tif')[0]
#im_file_scl = glob.glob(f'{masked_path}*31UFR*B04*.tif')[0] #im_file_scl = glob.glob(f'{masked_path}*31UFR*B04*.tif')[0]
src = rasterio.open(im_file_scl, "r") src = rasterio.open(im_file_scl, "r")
color_map = plt.cm.get_cmap("Greys_r") color_map = plt.cm.get_cmap("Greys_r")
color_map.set_bad(color='red') color_map.set_bad(color='red')
rasterio.plot.show(src, ax=ax2, cmap=color_map, vmin=p5, vmax=p95, title="After SCL") rasterio.plot.show(src, ax=ax2, cmap=color_map, vmin=p5, vmax=p95, title="After SCL")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-warning"> <div class="alert alert-warning">
<h4 class="alert-heading"><strong>Fourth exercise of the session</strong></h4> <h4 class="alert-heading"><strong>Fourth exercise of the session</strong></h4>
<hr> <hr>
### Exercise: Sum Cloud Pixels and Plot the Results ### 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**. 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: #### Steps:
1. Use `glob.glob` to list all the SCL masks paths 1. Use `glob.glob` to list all the SCL masks paths
2. Iterate through the list with a `for` loop 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 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**. 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. 5. Plot the **cloud pixel count** on the y-axis and **image dates** on the x-axis.
#### **Hints:** #### **Hints:**
- Use numpy to sum the number of pixels with categorical values equal to clouds - 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 - 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 - For plotting, use `matplotlib` or another plotting library
SCL class| Description SCL class| Description
:---------:|:----------: :---------:|:----------:
0| No data 0| No data
1| Saturated or defective 1| Saturated or defective
2| Dark area pixels 2| Dark area pixels
3| Cloud shadows 3| Cloud shadows
4| Vegetation 4| Vegetation
5| Not vegetated 5| Not vegetated
6| Water 6| Water
7| Unclassified 7| Unclassified
8| Cloud medium probability 8| Cloud medium probability
9| Cloud high probability 9| Cloud high probability
10| Thin cirrus 10| Thin cirrus
11| Snow 11| Snow
</div> </div>
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-info"> <div class="alert alert-info">
### How to Open an Image with `rasterio` and Convert it to a NumPy Array ### 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: 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()`. 1. **Open the raster file** using `rasterio.open()`.
2. **Read the image data** into a NumPy array using the `read()` method. 2. **Read the image data** into a NumPy array using the `read()` method.
Here’s an example: Here’s an example:
```python ```python
import rasterio import rasterio
# Open the image file # Open the image file
with rasterio.open('path_to_your_image.tif') as src: with rasterio.open('path_to_your_image.tif') as src:
# Convert the image to a NumPy array # Convert the image to a NumPy array
image_array = src.read(1) # '1' refers to the first band (change if needed) 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 # Now 'image_array' contains the raster data as a NumPy array
......
0% Chargement en cours ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter