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

classif 2025

parent fdde0106
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
Masquer les modifications d'espaces
En ligne Côte à côte
docs/source/classification/classification_2025.ipynb 0 → 100644
+ 777
− 0
Voir le fichier @ cbbbede1
%% Cell type:markdown id: tags:
## 1. In-situ data
In situ data refers to ground-truth measurements collected to train and validate satellite-based classification models. It can consist of points or polygons with an assigned class for a specific location and/or time.
In-situ data can be obtained via
- Photo-interpretation of satellite images
* This can be done manually in GIS software or with the help of tools (e.g. [AcATaMa](https://plugins.qgis.org/plugins/AcATaMa/))
- Existing databases
* e.g. agricultural census data (e.g. [Wallonia](https://geoportail.wallonie.be/catalogue/414cdf16-c697-4244-8f63-0ad6f1770400.html)), land cover databases, ...
%% Cell type:code id: tags:
``` python
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from rasterio import features
import rasterio.plot
from IPython.display import display
import plotly.express as px
import plotly.offline
plotly.offline.init_notebook_mode()
from pathlib import Path
import glob, os, time, math
import sklearn
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
print('All libraries successfully imported!')
print(f'Pandas : {pd.__version__}')
print(f'GeoPandas : {gpd.__version__}')
print(f'Scikit-learn : {sklearn.__version__}')
```
%% Cell type:markdown id: tags:
**Set directory**
%% Cell type:code id: tags:
``` python
computer_path = 'X:/'
grp_nb = '30'
# Directory for all work files
work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/'
data_path = f'{work_path}DATA/'
reflectance_path = f'{work_path}3_L2A_MASKED/'
in_situ_path = f'{data_path}A_IN_SITU/'
Path(in_situ_path).mkdir(parents=True, exist_ok=True)
```
%% Cell type:markdown id: tags:
<div class="alert alert-warning">
<h4 class="alert-heading"><strong>First Exercise of the Session</strong></h4>
<hr>
### **Exercise: Create Labeled Training Data as a Shapefile in QGIS**
Your task is to create a **shapefile** in QGIS containing labeled training data for satellite image classification. The shapefile must include an **attribute table** with the following variables:
- **`id`** *(int)*: A unique identifier for each sample.
- **`label`** *(str)*: The class name (e.g., "Water", "Forest", "Urban").
- **`labelID`** *(int)*: A numerical code corresponding to the class label.
Use the [AcaTama plug-in](https://plugins.qgis.org/plugins/AcATaMa/) to generate a random set of points inside your image. Use these points to draw the polygons.
#### **Steps:**
1. **Open QGIS and Load a Basemap:**
- Open QGIS and add a relevant satellite image or basemap as a reference.
2. **Create a random sample of points using the AcaTama plug-in**
- Go to **Plugins**
- Open the AcaTama plugin (first install if needed)
- Select your clipped image as the thematic map
- Generate random points in the sampling design window
2. **Create a New Shapefile Layer:**
- Go to **Layer** → **Create Layer** → **New Shapefile Layer**.
- Save the layer as **`in_situ.shp`** in the folder **`DATA/A_IN_SITU`**
- Select **"Polygon"** as the geometry type (or "Point" if appropriate).
- Set the **CRS (Coordinate Reference System)** to match your satellite image.
3. **Create the following fields:**
- `id` (Integer)
- `label` (Text/String)
- `labelID` (Integer)
4. **Digitize Training Samples:**
- Activate **Editing Mode** (click the pencil icon).
- Use the **"Add Feature"** tool to draw polygons (or points) over different land cover types.
- In the pop-up window, assign a **unique id, a class label, and a numerical labelID**.
5. **Save and Export:**
- Stop editing and save changes.
- Ensure the shapefile is properly stored and named **`in_situ.shp`**.
#### **Hints:**
- Be consistent with labeling (e.g., `labelID = 1` for "Water", `labelID = 2` for "Forest", etc.).
- Ensure your samples cover a variety of land cover types to improve classification accuracy.
- You can later use this shapefile as ground truth data for supervised classification.
</div>
%% Cell type:code id: tags:
``` python
### Read & print in situ data
in_situ = gpd.read_file(f'{in_situ_path}in_situ.shp')
in_situ.head()
### Split in training & validation data
cal_gdf, val_gdf = train_test_split(in_situ, test_size=0.25, random_state=4916)
print(f"Training data: {len(cal_gdf)} samples")
print(f"Validation data: {len(val_gdf)} samples")
cal_gdf.to_file(f"{in_situ_path}in_situ_cal.shp")
val_gdf.to_file(f"{in_situ_path}in_situ_val.shp")
```
%% Cell type:markdown id: tags:
## Supervised classification with the Random Forest algorithm
The classification step consists in one or many numerical processes to finally allocate every pixel or object to one of the classes of the land cover typology. The vast diversity of classification algorithms can be split into two main types:
- the **supervised type**, which uses a training data set to calibrate the algorithm a priori;
- and **the unsupervised type**, which produces clusters of pixels to be labelled a posteriori as land cover class in light of in situ or ancillary information.
> **The random forest is a model consisting hundreds of decision trees. Each decision tree is trained on a subset of observations and features. The final predictions of the random forest are made by averaging the predictions of each individual tree.**
[Watch this video if Decision Trees are not clear for you !](https://www.youtube.com/watch?v=7VeUPuFGJHk)
[Watch this video if Random Forest are not clear for you !](https://www.youtube.com/watch?v=J4Wdy0Wc_xQ)
In this chapter we will see how to use the Random Forest implementation provided by the `scikit-learn` library.
<figure class="image">
<img src="im_classif_pixel.png" alt="Image classification" width="600">
</figure>
<figure class="image">
<img src="classif_namur_2020.png" alt="Image classification" width="800">
<figcaption>Namur, 2020 (NDVI & monthly composites S1 backscattering VV)</figcaption>
</figure>
---
[Handbook on remote sensing for agricultural statistics](https://nicolasdeffense.github.io/eo-toolbox/docs/Remote_Sensing_for_Agricultural_Statistics.pdf)
[Chris Holden's tutorial](https://ceholden.github.io/open-geo-tutorial/python/chapter_5_classification.html)
[An Implementation and Explanation of the Random Forest in Python](https://towardsdatascience.com/an-implementation-and-explanation-of-the-random-forest-in-python-77bf308a9b76)
%% Cell type:code id: tags:
``` python
# Input directories
in_situ_path = f'{data_path}A_IN_SITU/'
s2_path = f'{work_path}3_L2A_MASKED/'
# Output directory
classif_path = f'{work_path}CLASSIF/'
# no data value
no_data = -10000
Path(classif_path).mkdir(parents=True, exist_ok=True)
print(f'Classification path is set to : {classif_path}')
```
%% Cell type:markdown id: tags:
## 1. Prepare classification features associated to *in situ* data
%% Cell type:markdown id: tags:
### 1.1 Rasterize *in situ* data calibration shapefile
%% Cell type:markdown id: tags:
<div class="alert alert-success">
### Function: `rasterize(in_situ_gdf, in_situ_cal_tif, img_temp_path, no_data=-10000)`
This function **rasterizes a vector dataset** by burning geometries from an input GeoDataFrame into a raster template.
#### **Arguments:**
1. **`in_situ_gdf`** *(GeoDataFrame)* → A GeoDataFrame containing the vector data to be rasterized, including geometries and class labels.
2. **`in_situ_cal_tif`** *(str)* → The file path where the rasterized output will be saved.
3. **`img_temp_path`** *(str)* → The file path to an existing raster file, used as a template for spatial reference and resolution.
4. **`no_data`** *(int, optional)* → The NoData value for the output raster. Defaults to `-10000`.
</div>
%% Cell type:code id: tags:
``` python
# Open the calibration polygons with GeoPandas
in_situ_cal_shp = f"{in_situ_path}in_situ_cal.shp"
in_situ_gdf = gpd.read_file(in_situ_cal_shp)
in_situ_cal_tif = f'{in_situ_path}IN_SITU_ROI_CAL.tif'
# Open the raster file you want to use as a template for rasterize
img_temp_path = glob.glob(f'{s2_path}*.tif')[0]
from func_classif import rasterize
rasterize(in_situ_gdf, in_situ_cal_tif, img_temp_path, no_data)
```
%% Cell type:markdown id: tags:
### 1.2 List all the classification features
%% Cell type:markdown id: tags:
<div class="alert alert-warning">
<h4 class="alert-heading"><strong>Second exercise: Load and Store Feature Rasters</strong></h4>
<hr>
### **Task:**
Write a Python script to load satellite image features stored as **GeoTIFF files** and save them as arrays for further processing.
#### **Steps:**
1. **Define a List of Features**
- Create a list named `features` with the feature names you want to load (e.g., `"NDVI"`, `"EVI"`).
2. **Initialize an Empty List**
- Define `list_src_arr = []` to store raster arrays.
3. **Search and Load Raster Files**
- Use `glob.glob()` to find `.tif` files for each feature.
- Adapt the search path according to your file structure.
- Open each file with `rasterio`, read band **1**, append it to `list_src_arr`, and close the file.
4. **Print Summary**
- Print the shape of the last loaded feature and the total number of features stored.
#### **Hints:**
- Modify the search path in `glob.glob()` to match your dataset structure.
- You can easily extend the script by adding more features to `features`.
</div>
%% Cell type:code id: tags:
``` python
features = ["NDVI", "NBR"]
list_src_arr = [] # Leave empty, this will be used to store the arrays
feature_names = [] # Leave empty, this will be used to store the feature names
for feature in features:
list_im = glob.glob(f"{work_path}/{feature}/*2023*.tif")
feature_names += [path.split("\\")[-1] for path in list_im]
for im_file in list_im:
print(im_file)
with rasterio.open(im_file, 'r') as src:
im = src.read(1)
list_src_arr.append(im)
```
%% Cell type:markdown id: tags:
Merge all the 2D matrices from the list into one 3D matrix
%% Cell type:code id: tags:
``` python
feat_arr = np.dstack(list_src_arr).astype(np.float32)
print(feat_arr.shape)
print(f'There are {feat_arr.shape[2]} features')
print(f'The features type is : {feat_arr.dtype}')
```
%% Cell type:markdown id: tags:
### 1.3 Pairing *in situ* data (Y) with EO classification features (X)
Now that we have the image we want to classify (our X feature inputs), and the ROI with the land cover labels (our Y labeled data), we need to pair them up in NumPy arrays so we may feed them to Random Forest.
%% Cell type:code id: tags:
``` python
# Open in-situ used for calibration
src = rasterio.open(in_situ_cal_tif, "r")
cal_arr = src.read(1)
src.close()
# Find how many labeled entries we have -- i.e. how many training data samples?
n_samples = (cal_arr != no_data).sum()
print(f'We have {n_samples} samples (= calibration pixels)')
```
%% Cell type:markdown id: tags:
What are our classification labels?
%% Cell type:code id: tags:
``` python
labels = np.unique(cal_arr[cal_arr != no_data])
print(f'The training data include {labels.size} classes: {labels}')
```
%% Cell type:markdown id: tags:
We need :
- **"X" 2D matrix** containing classification features
- **"y" 1D matrix** containing our labels
These will have `n_samples` rows.
%% Cell type:code id: tags:
``` python
X = feat_arr[cal_arr != no_data, :]
y = cal_arr[cal_arr != no_data]
# Replace NaN in classification features by the no_data value
X = np.nan_to_num(X, nan=no_data)
print(f'Our X matrix is sized: {X.shape}')
print(f'Our y array is sized: {y.shape}')
```
%% Cell type:markdown id: tags:
## 2. Train the Random Forest model
%% Cell type:markdown id: tags:
Now that we have our X 2D-matrix of feature inputs and our y 1D-matrix containing the labels, we can train our model.
Visit this <a href="https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html" target="_blank">web page</a> to find the usage of RandomForestClassifier from scikit-learn.
%% Cell type:code id: tags:
``` python
start_training = time.time()
# Initialize our model
rf = RandomForestClassifier(n_estimators=100, # The number of trees in the forest.
bootstrap=True, # Whether bootstrap samples are used when building trees. If False, the whole dataset is used to build each tree.
oob_score=True) # Whether to use out-of-bag samples to estimate the generalization score. Only available if bootstrap=True.
# Fit our model to training data
rf = rf.fit(X, y)
end_training = time.time()
# Get time elapsed during the Random Forest training
hours, rem = divmod(end_training-start_training, 3600)
minutes, seconds = divmod(rem, 60)
print("Random Forest training : {:0>2}:{:0>2}:{:05.2f}".format(int(hours),int(minutes),seconds))
```
%% Cell type:markdown id: tags:
With our Random Forest model fit, we can check out the "Out-of-Bag" (OOB) prediction score.
> Score of the training dataset obtained using an out-of-bag estimate. This attribute exists only when oob_score is True.
%% Cell type:code id: tags:
``` python
print(f'Our OOB prediction of accuracy is: {round(rf.oob_score_ * 100,2)}%')
```
%% Cell type:markdown id: tags:
To help us get an idea of which features bands were important, we can look at the feature importance scores.
> The impurity-based feature importances. The higher, the more important the feature. The importance of a feature is computed as the (normalized) total reduction of the criterion brought by that feature. It is also known as the Gini importance.
%% Cell type:code id: tags:
``` python
importance_df = pd.DataFrame({"Importance": rf.feature_importances_, "Feature": feature_names})
importance_df = importance_df.sort_values(by="Importance", ascending=False)
importance_df
```
%% Cell type:markdown id: tags:
Let's look at a crosstabulation to see the class confusion
%% Cell type:code id: tags:
``` python
# Setup a dataframe
df = pd.DataFrame()
df['truth'] = y
df['predict'] = rf.predict(X)
# Cross-tabulate predictions
cross_tab = pd.crosstab(df['truth'], df['predict'], margins=True)
display(cross_tab)
```
%% Cell type:markdown id: tags:
Unbelievable?
Given enough information and effort, this algorithm precisely learned what we gave it. Asking to validate a machine learning algorithm on the training data is a useless exercise that will overinflate the accuracy.
Instead, we could have done a crossvalidation approach where we train on a subset the dataset, and then predict and assess the accuracy using the sections we didn't train it on.
%% Cell type:markdown id: tags:
## 3. Predict the rest of the image
With our Random Forest classifier fit, we can now proceed by trying to classify the entire image.
%% Cell type:markdown id: tags:
<div class="alert alert-success">
### Function: `apply_rf_to_raster(rf, img, no_data=-10000)`
This function **applies a trained Random Forest (RF) model** to classify each pixel in a raster image.
#### **Arguments:**
1. **`rf`** *(RandomForestClassifier)* → A pre-trained Random Forest classifier used for pixel-wise classification.
2. **`img`** *(numpy array)* → A 3D NumPy array representing the raster image, where dimensions are `(rows, columns, bands)`.
3. **`no_data`** *(int, optional)* → The NoData value indicating missing or invalid pixels. Defaults to `-10000`.
#### **Process:**
- Reshapes the input raster from a 3D shape `(rows, columns, bands)` to a 2D array `(rows * columns, bands)`, making it suitable for classification.
- Uses the trained Random Forest model (`rf`) to predict the class label for each pixel.
- Replaces predictions with `no_data` for pixels where any input feature has the `no_data` value.
- Reshapes the classification output back to the original 2D raster shape `(rows, columns)`.
- Measures and prints the time taken for classification.
- Returns the classified raster as a NumPy array.
</div>
%% Cell type:code id: tags:
``` python
from func_classif import apply_rf_to_raster
class_prediction = apply_rf_to_raster(rf, feat_arr, no_data)
print(class_prediction)
```
%% Cell type:markdown id: tags:
## 4. Filter classification with moving window
%% Cell type:markdown id: tags:
<div class="alert alert-success">
### Function: `smooth_classification(class_prediction, window_size, no_data=-10000)`
This function **applies a majority filter** to smooth a classified raster by replacing each pixel with the most frequent class within a moving window, ensuring proper centering.
#### **Arguments:**
1. **`class_prediction`** *(numpy array)* → A 2D NumPy array representing the classified raster, where each pixel has a class label.
2. **`window_size`** *(int)* → The size of the moving window (must be an odd integer).
3. **`no_data`** *(int, optional)* → The NoData value indicating missing or invalid pixels. Defaults to `-10000`.
#### **Process:**
- Determines the image dimensions (`sizey`, `sizex`).
- Computes the window radius (`rad_wind`) as `floor(window_size / 2)`.
- Pads the input raster using edge padding to handle boundary pixels.
- Initializes an empty array `majority` to store the smoothed classification results.
- Iterates over each pixel in the image:
- If the pixel itself is `no_data`, it remains `no_data`.
- Extracts a properly centered `window_size × window_size` neighborhood around the pixel.
- Flattens the window and removes `no_data` values.
- Determines the most frequent class within the window and assigns it to the pixel.
- Reshapes the output array to match the original image dimensions.
- Returns the smoothed classification as a NumPy array with shape `(1, sizey, sizex)`.
</div>
%% Cell type:code id: tags:
``` python
ws = 3 # window size for smoothing, use an uneven number so the pixel of interest is in the center
```
%% Cell type:code id: tags:
``` python
from func_classif import smooth_classification
majority = smooth_classification(class_prediction, ws)
print(f'Classification : \n {class_prediction}')
print(f'Classification with filter : \n {majority}')
```
%% Cell type:markdown id: tags:
## 5. Write classification products into GeoTIFF files
%% Cell type:markdown id: tags:
Open template image to get metadata
%% Cell type:code id: tags:
``` python
with rasterio.open(img_temp_path) as src:
profile = src.profile
profile
```
%% Cell type:markdown id: tags:
**Write classification**
%% Cell type:code id: tags:
``` python
classif_tif = f'{classif_path}CLASSIF_RF.tif'
with rasterio.open(classif_tif, "w", **profile) as dst:
dst.write(class_prediction, 1)
```
%% Cell type:markdown id: tags:
**Write classification with moving window filtering**
%% Cell type:code id: tags:
``` python
reclassif_filter_tif = f'{classif_path}CLASSIF_RF_FILTER.tif'
with rasterio.open(reclassif_filter_tif, "w", **profile) as dst:
dst.write(majority)
```
%% Cell type:markdown id: tags:
<div class="alert alert-warning">
<h4 class="alert-heading"><strong> Third exercise: Visualizing and Improving Classifications in QGIS</strong></h4>
<hr>
### **Task:**
Open your classified raster in **QGIS**, visualize the results with distinct colors for each class, and assess how to improve the classification.
1. **Load the Classification Raster**
- Open **QGIS** and add your classified raster layer.
2. **Apply Unique Colors to Each Class**
- Open the **Layer Properties** → **Symbology** tab.
- Set **Render Type** to **"Unique Values" (Paletted/Unique Values)**.
- Click **Classify** to assign colors to each class.
3. **Suggest Improvements**
- Consider adding other spectral indices/ other images/composites to the features.
- Consider adding training samples and/or change the classes (e.g. distinguish between croplands and grasslands)
- Consider changing the window size for smoothing the classification.
#### **Hints:**
- Zoom into different regions to identify classification errors.
- Document misclassified areas and suggest possible fixes.
</div>
docs/source/classification/func_classif.py 0 → 100644
+ 113
− 0
Voir le fichier @ cbbbede1
import numpy as np
import math, time
import rasterio
from rasterio import features
def rasterize(in_situ_gdf, in_situ_cal_tif, img_temp_path, no_data = -10000):
print(f'Raster template file : {img_temp_path}')
src = rasterio.open(img_temp_path, "r")
# Update metadata
out_meta = src.meta
out_meta.update(nodata=no_data)
crs_shp = str(in_situ_gdf.crs).split(":",1)[1]
crs_tif = str(src.crs).split(":",1)[1]
print(f'The CRS of in situ data is : {crs_shp}')
print(f'The CRS of raster template is : {crs_tif}')
if crs_shp == crs_tif:
print("CRS are the same")
print(f'Rasterize starts')
# Burn the features into the raster and write it out
dst = rasterio.open(in_situ_cal_tif, 'w+', **out_meta)
dst_arr = dst.read(1)
# This is where we create a generator of geom, value pairs to use in rasterizing
geom_col = in_situ_gdf.geometry
code_col = in_situ_gdf["labelID"].astype(int)
shapes = ((geom,value) for geom, value in zip(geom_col, code_col))
in_situ_arr = features.rasterize(shapes=shapes,
fill=no_data,
out=dst_arr,
transform=dst.transform)
dst.write_band(1, in_situ_arr)
print(f'Rasterize is done : {in_situ_cal_tif}')
# Close rasterio objects
src.close()
dst.close()
else:
print('CRS are different --> repoject in-situ data shapefile with "to_crs"')
def apply_rf_to_raster(rf, img, no_data = -10000):
# Take our full image and reshape into long 2d array (nrow * ncol, nband) for classification
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :].reshape(new_shape)
print(f'Reshaped from {img.shape} to {img_as_array.shape}')
start_classification = time.time()
# Now predict for each pixel
class_prediction = rf.predict(img_as_array)
# predictions with no_data in any of the input features are replaced by the no_data value
class_prediction[np.any(img_as_array == no_data, axis=1)] = no_data
# Reshape our classification map
class_prediction = class_prediction.reshape(img[:, :, 0].shape)
end_classification = time.time()
hours, rem = divmod(end_classification-start_classification, 3600)
minutes, seconds = divmod(rem, 60)
print("Random Forest prediction : {:0>2}:{:0>2}:{:05.2f}".format(int(hours),int(minutes),seconds))
return class_prediction
def smooth_classification(class_prediction, window_size, no_data = -10000):
sizey = class_prediction.shape[0]
sizex = class_prediction.shape[1]
rad_wind = math.floor(window_size/2)
X = np.pad(class_prediction, ((rad_wind,rad_wind),(rad_wind,rad_wind)), 'edge')
majority = np.empty((sizey,sizex), dtype='int16')
for i in range(sizey):
for j in range(sizex):
window = X[i:i+window_size , j:j+window_size]
window = window.flatten()
window = window[window != no_data]
if len(window) == 0:
majority[i,j] = no_data
else:
counts = np.bincount(window)
maj = np.argmax(counts)
majority[i,j]= maj
majority = majority.reshape((1,sizey,sizex))
return majority
\ No newline at end of file
0% 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