classification_2025.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. In-situ data\n",
    "\n",
    "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.\n",
    "\n",
    "In-situ data can be obtained via\n",
    "\n",
    "- Photo-interpretation of satellite images\n",
    "    * This can be done manually in GIS software or with the help of tools (e.g. [AcATaMa](https://plugins.qgis.org/plugins/AcATaMa/))\n",
    "- Existing databases\n",
    "    * 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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio import features\n",
    "import rasterio.plot\n",
    "from IPython.display import display\n",
    "\n",
    "import plotly.express as px\n",
    "import plotly.offline\n",
    "plotly.offline.init_notebook_mode()\n",
    "from pathlib import Path\n",
    "\n",
    "import glob, os, time, math\n",
    "\n",
    "import sklearn\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "print('All libraries successfully imported!')\n",
    "print(f'Pandas    : {pd.__version__}')\n",
    "print(f'GeoPandas : {gpd.__version__}')\n",
    "print(f'Scikit-learn : {sklearn.__version__}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set directory**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "computer_path = 'X:/'\n",
    "grp_nb        = '30'\n",
    "\n",
    "# Directory for all work files\n",
    "work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/'\n",
    "data_path = f'{work_path}DATA/'\n",
    "\n",
    "reflectance_path = f'{work_path}3_L2A_MASKED/'\n",
    "\n",
    "in_situ_path = f'{data_path}A_IN_SITU/'\n",
    "\n",
    "Path(in_situ_path).mkdir(parents=True, exist_ok=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "<h4 class=\"alert-heading\"><strong>First Exercise of the Session</strong></h4>  \n",
    "<hr>  \n",
    "\n",
    "### **Exercise: Create Labeled Training Data as a Shapefile in QGIS**  \n",
    "\n",
    "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:  \n",
    "\n",
    "- **`id`** *(int)*: A unique identifier for each sample.  \n",
    "- **`label`** *(str)*: The class name (e.g., \"Water\", \"Forest\", \"Urban\").  \n",
    "- **`labelID`** *(int)*: A numerical code corresponding to the class label.  \n",
    "\n",
    "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.\n",
    "\n",
    "#### **Steps:**  \n",
    "\n",
    "1. **Open QGIS and Load a Basemap:**  \n",
    "   - Open QGIS and add a relevant satellite image or basemap as a reference.  \n",
    "\n",
    "2. **Create a random sample of points using the AcaTama plug-in**\n",
    "   - Go to **Plugins**\n",
    "   - Open the AcaTama plugin (first install if needed)\n",
    "   - Select your clipped image as the thematic map\n",
    "   - Generate random points in the sampling design window\n",
    "\n",
    "2. **Create a New Shapefile Layer:**  \n",
    "   - Go to **Layer** → **Create Layer** → **New Shapefile Layer**.  \n",
    "   - Save the layer as **`in_situ.shp`** in the folder **`DATA/A_IN_SITU`** \n",
    "   - Select **\"Polygon\"** as the geometry type (or \"Point\" if appropriate).  \n",
    "   - Set the **CRS (Coordinate Reference System)** to match your satellite image.  \n",
    "\n",
    "3. **Create the following fields:**\n",
    "    - `id` (Integer)  \n",
    "    - `label` (Text/String)  \n",
    "    - `labelID` (Integer)  \n",
    "\n",
    "4. **Digitize Training Samples:**  \n",
    "   - Activate **Editing Mode** (click the pencil icon).  \n",
    "   - Use the **\"Add Feature\"** tool to draw polygons (or points) over different land cover types.  \n",
    "   - In the pop-up window, assign a **unique id, a class label, and a numerical labelID**.  \n",
    "\n",
    "5. **Save and Export:**  \n",
    "   - Stop editing and save changes.  \n",
    "   - Ensure the shapefile is properly stored and named **`in_situ.shp`**.  \n",
    "\n",
    "#### **Hints:**  \n",
    "\n",
    "- Be consistent with labeling (e.g., `labelID = 1` for \"Water\", `labelID = 2` for \"Forest\", etc.).  \n",
    "- Ensure your samples cover a variety of land cover types to improve classification accuracy.  \n",
    "- You can later use this shapefile as ground truth data for supervised classification.  \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Read & print in situ data\n",
    "in_situ = gpd.read_file(f'{in_situ_path}in_situ.shp')\n",
    "in_situ.head()\n",
    "\n",
    "### Split in training & validation data\n",
    "cal_gdf, val_gdf = train_test_split(in_situ, test_size=0.25, random_state=4916)\n",
    "\n",
    "print(f\"Training data: {len(cal_gdf)} samples\")\n",
    "print(f\"Validation data: {len(val_gdf)} samples\")\n",
    "\n",
    "cal_gdf.to_file(f\"{in_situ_path}in_situ_cal.shp\")\n",
    "val_gdf.to_file(f\"{in_situ_path}in_situ_val.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supervised classification with the Random Forest algorithm\n",
    "\n",
    "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:\n",
    "- the **supervised type**, which uses a training data set to calibrate the algorithm a priori;\n",
    "- 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.\n",
    "\n",
    "> **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.**\n",
    "\n",
    "[Watch this video if Decision Trees are not clear for you !](https://www.youtube.com/watch?v=7VeUPuFGJHk)\n",
    "\n",
    "[Watch this video if Random Forest are not clear for you !](https://www.youtube.com/watch?v=J4Wdy0Wc_xQ)\n",
    "\n",
    "\n",
    "In this chapter we will see how to use the Random Forest implementation provided by the `scikit-learn` library.\n",
    "\n",
    "<figure class=\"image\">\n",
    "  <img src=\"im_classif_pixel.png\" alt=\"Image classification\" width=\"600\">\n",
    "</figure>\n",
    "\n",
    "\n",
    "\n",
    "<figure class=\"image\">\n",
    "  <img src=\"classif_namur_2020.png\" alt=\"Image classification\" width=\"800\">\n",
    "  <figcaption>Namur, 2020 (NDVI & monthly composites S1 backscattering VV)</figcaption>\n",
    "</figure>\n",
    "\n",
    "---\n",
    "\n",
    "[Handbook on remote sensing for agricultural statistics](https://nicolasdeffense.github.io/eo-toolbox/docs/Remote_Sensing_for_Agricultural_Statistics.pdf)\n",
    "\n",
    "[Chris Holden's tutorial](https://ceholden.github.io/open-geo-tutorial/python/chapter_5_classification.html)\n",
    "\n",
    "[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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Input directories\n",
    "in_situ_path = f'{data_path}A_IN_SITU/'\n",
    "s2_path      = f'{work_path}3_L2A_MASKED/'\n",
    "\n",
    "# Output directory\n",
    "classif_path = f'{work_path}CLASSIF/'\n",
    "\n",
    "# no data value\n",
    "no_data = -10000\n",
    "\n",
    "Path(classif_path).mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "print(f'Classification path is set to : {classif_path}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Prepare classification features associated to *in situ* data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Rasterize *in situ* data calibration shapefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "### Function: `rasterize(in_situ_gdf, in_situ_cal_tif, img_temp_path, no_data=-10000)`\n",
    "\n",
    "This function **rasterizes a vector dataset** by burning geometries from an input GeoDataFrame into a raster template.  \n",
    "\n",
    "#### **Arguments:**  \n",
    "1. **`in_situ_gdf`** *(GeoDataFrame)* → A GeoDataFrame containing the vector data to be rasterized, including geometries and class labels.  \n",
    "2. **`in_situ_cal_tif`** *(str)* → The file path where the rasterized output will be saved.  \n",
    "3. **`img_temp_path`** *(str)* → The file path to an existing raster file, used as a template for spatial reference and resolution.  \n",
    "4. **`no_data`** *(int, optional)* → The NoData value for the output raster. Defaults to `-10000`.  \n",
    "\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the calibration polygons with GeoPandas\n",
    "in_situ_cal_shp = f\"{in_situ_path}in_situ_cal.shp\"\n",
    "in_situ_gdf = gpd.read_file(in_situ_cal_shp)\n",
    "\n",
    "in_situ_cal_tif = f'{in_situ_path}IN_SITU_ROI_CAL.tif'\n",
    "\n",
    "# Open the raster file you want to use as a template for rasterize\n",
    "img_temp_path = glob.glob(f'{s2_path}*.tif')[0]\n",
    "\n",
    "from func_classif import rasterize\n",
    "\n",
    "rasterize(in_situ_gdf, in_situ_cal_tif, img_temp_path, no_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 List all the classification features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "<h4 class=\"alert-heading\"><strong>Second exercise: Load and Store Feature Rasters</strong></h4>  \n",
    "<hr>  \n",
    "\n",
    "### **Task:**  \n",
    "\n",
    "Write a Python script to load satellite image features stored as **GeoTIFF files** and save them as arrays for further processing.  \n",
    "\n",
    "#### **Steps:**  \n",
    "\n",
    "1. **Define a List of Features**  \n",
    "   - Create a list named `features` with the feature names you want to load (e.g., `\"NDVI\"`, `\"EVI\"`).  \n",
    "\n",
    "2. **Initialize an Empty List**  \n",
    "   - Define `list_src_arr = []` to store raster arrays.  \n",
    "\n",
    "3. **Search and Load Raster Files**  \n",
    "   - Use `glob.glob()` to find `.tif` files for each feature.  \n",
    "   - Adapt the search path according to your file structure.  \n",
    "   - Open each file with `rasterio`, read band **1**, append it to `list_src_arr`, and close the file.  \n",
    "\n",
    "4. **Print Summary**  \n",
    "   - Print the shape of the last loaded feature and the total number of features stored.  \n",
    "\n",
    "#### **Hints:**  \n",
    "\n",
    "- Modify the search path in `glob.glob()` to match your dataset structure.  \n",
    "- You can easily extend the script by adding more features to `features`.  \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = [\"NDVI\", \"NBR\"]\n",
    "\n",
    "list_src_arr = [] # Leave empty, this will be used to store the arrays\n",
    "feature_names = [] # Leave empty, this will be used to store the feature names\n",
    "\n",
    "for feature in features:\n",
    "    list_im = glob.glob(f\"{work_path}/{feature}/*2023*.tif\")\n",
    "    feature_names += [path.split(\"\\\\\")[-1] for path in list_im]\n",
    "\n",
    "    for im_file in list_im:\n",
    "        print(im_file)\n",
    "        with rasterio.open(im_file, 'r') as src:\n",
    "            im = src.read(1)\n",
    "        list_src_arr.append(im)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge all the 2D matrices from the list into one 3D matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat_arr = np.dstack(list_src_arr).astype(np.float32)\n",
    "\n",
    "print(feat_arr.shape)\n",
    "print(f'There are {feat_arr.shape[2]} features')\n",
    "print(f'The features type is : {feat_arr.dtype}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Pairing *in situ* data (Y) with EO classification features (X)\n",
    "\n",
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open in-situ used for calibration\n",
    "\n",
    "src = rasterio.open(in_situ_cal_tif, \"r\")\n",
    "cal_arr = src.read(1)\n",
    "src.close()\n",
    "\n",
    "# Find how many labeled entries we have -- i.e. how many training data samples?\n",
    "n_samples = (cal_arr != no_data).sum()\n",
    "\n",
    "print(f'We have {n_samples} samples (= calibration pixels)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are our classification labels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = np.unique(cal_arr[cal_arr != no_data])\n",
    "\n",
    "print(f'The training data include {labels.size} classes: {labels}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need :\n",
    "- **\"X\" 2D matrix** containing classification features\n",
    "- **\"y\" 1D matrix** containing our labels\n",
    "\n",
    "These will have `n_samples` rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = feat_arr[cal_arr != no_data, :]\n",
    "y = cal_arr[cal_arr != no_data]\n",
    "\n",
    "# Replace NaN in classification features by the no_data value\n",
    "X = np.nan_to_num(X, nan=no_data)\n",
    "\n",
    "print(f'Our X matrix is sized: {X.shape}')\n",
    "print(f'Our y array is sized: {y.shape}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Train the Random Forest model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our X 2D-matrix of feature inputs and our y 1D-matrix containing the labels, we can train our model.\n",
    "\n",
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_training = time.time()\n",
    "\n",
    "# Initialize our model\n",
    "rf = RandomForestClassifier(n_estimators=100, # The number of trees in the forest.\n",
    "                            bootstrap=True,   # Whether bootstrap samples are used when building trees. If False, the whole dataset is used to build each tree.\n",
    "                            oob_score=True)   # Whether to use out-of-bag samples to estimate the generalization score. Only available if bootstrap=True.\n",
    "\n",
    "# Fit our model to training data\n",
    "rf = rf.fit(X, y)\n",
    "\n",
    "end_training = time.time()\n",
    "\n",
    "# Get time elapsed during the Random Forest training\n",
    "hours, rem = divmod(end_training-start_training, 3600)\n",
    "minutes, seconds = divmod(rem, 60)\n",
    "print(\"Random Forest training : {:0>2}:{:0>2}:{:05.2f}\".format(int(hours),int(minutes),seconds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our Random Forest model fit, we can check out the \"Out-of-Bag\" (OOB) prediction score.\n",
    "\n",
    "> Score of the training dataset obtained using an out-of-bag estimate. This attribute exists only when oob_score is True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Our OOB prediction of accuracy is: {round(rf.oob_score_ * 100,2)}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To help us get an idea of which features bands were important, we can look at the feature importance scores.\n",
    "\n",
    "> 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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importance_df = pd.DataFrame({\"Importance\": rf.feature_importances_, \"Feature\": feature_names})\n",
    "importance_df = importance_df.sort_values(by=\"Importance\", ascending=False)\n",
    "importance_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a crosstabulation to see the class confusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup a dataframe\n",
    "df = pd.DataFrame()\n",
    "\n",
    "df['truth'] = y\n",
    "df['predict'] = rf.predict(X)\n",
    "\n",
    "# Cross-tabulate predictions\n",
    "\n",
    "cross_tab = pd.crosstab(df['truth'], df['predict'], margins=True)\n",
    "display(cross_tab)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unbelievable?\n",
    "\n",
    "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.\n",
    "\n",
    "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",
   "metadata": {},
   "source": [
    "## 3. Predict the rest of the image\n",
    "\n",
    "With our Random Forest classifier fit, we can now proceed by trying to classify the entire image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "### Function: `apply_rf_to_raster(rf, img, no_data=-10000)`\n",
    "\n",
    "This function **applies a trained Random Forest (RF) model** to classify each pixel in a raster image.  \n",
    "\n",
    "#### **Arguments:**  \n",
    "1. **`rf`** *(RandomForestClassifier)* → A pre-trained Random Forest classifier used for pixel-wise classification.  \n",
    "2. **`img`** *(numpy array)* → A 3D NumPy array representing the raster image, where dimensions are `(rows, columns, bands)`.  \n",
    "3. **`no_data`** *(int, optional)* → The NoData value indicating missing or invalid pixels. Defaults to `-10000`.  \n",
    "\n",
    "#### **Process:**  \n",
    "- Reshapes the input raster from a 3D shape `(rows, columns, bands)` to a 2D array `(rows * columns, bands)`, making it suitable for classification.  \n",
    "- Uses the trained Random Forest model (`rf`) to predict the class label for each pixel.  \n",
    "- Replaces predictions with `no_data` for pixels where any input feature has the `no_data` value.  \n",
    "- Reshapes the classification output back to the original 2D raster shape `(rows, columns)`.  \n",
    "- Measures and prints the time taken for classification.  \n",
    "- Returns the classified raster as a NumPy array.  \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from func_classif import apply_rf_to_raster\n",
    "\n",
    "class_prediction = apply_rf_to_raster(rf, feat_arr, no_data)\n",
    "\n",
    "print(class_prediction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Filter classification with moving window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "### Function: `smooth_classification(class_prediction, window_size, no_data=-10000)`\n",
    "\n",
    "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.  \n",
    "\n",
    "#### **Arguments:**  \n",
    "1. **`class_prediction`** *(numpy array)* → A 2D NumPy array representing the classified raster, where each pixel has a class label.  \n",
    "2. **`window_size`** *(int)* → The size of the moving window (must be an odd integer).  \n",
    "3. **`no_data`** *(int, optional)* → The NoData value indicating missing or invalid pixels. Defaults to `-10000`.  \n",
    "\n",
    "#### **Process:**  \n",
    "- Determines the image dimensions (`sizey`, `sizex`).  \n",
    "- Computes the window radius (`rad_wind`) as `floor(window_size / 2)`.  \n",
    "- Pads the input raster using edge padding to handle boundary pixels.  \n",
    "- Initializes an empty array `majority` to store the smoothed classification results.  \n",
    "- Iterates over each pixel in the image:  \n",
    "  - If the pixel itself is `no_data`, it remains `no_data`.  \n",
    "  - Extracts a properly centered `window_size × window_size` neighborhood around the pixel.  \n",
    "  - Flattens the window and removes `no_data` values.  \n",
    "  - Determines the most frequent class within the window and assigns it to the pixel.  \n",
    "- Reshapes the output array to match the original image dimensions.  \n",
    "- Returns the smoothed classification as a NumPy array with shape `(1, sizey, sizex)`.  \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = 3 # window size for smoothing, use an uneven number so the pixel of interest is in the center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from func_classif import smooth_classification\n",
    "\n",
    "majority = smooth_classification(class_prediction, ws)\n",
    "\n",
    "print(f'Classification : \\n {class_prediction}')\n",
    "print(f'Classification with filter : \\n {majority}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Write classification products into GeoTIFF files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open template image to get metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(img_temp_path) as src:\n",
    "    profile = src.profile\n",
    "\n",
    "profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Write classification**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classif_tif = f'{classif_path}CLASSIF_RF.tif'\n",
    "\n",
    "with rasterio.open(classif_tif, \"w\", **profile) as dst:\n",
    "    dst.write(class_prediction, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Write classification with moving window filtering**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reclassif_filter_tif = f'{classif_path}CLASSIF_RF_FILTER.tif'\n",
    "\n",
    "with rasterio.open(reclassif_filter_tif, \"w\", **profile) as dst:\n",
    "    dst.write(majority)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "<h4 class=\"alert-heading\"><strong> Third exercise: Visualizing and Improving Classifications in QGIS</strong></h4>  \n",
    "<hr>  \n",
    "\n",
    "### **Task:**  \n",
    "\n",
    "Open your classified raster in **QGIS**, visualize the results with distinct colors for each class, and assess how to improve the classification.  \n",
    "\n",
    "1. **Load the Classification Raster**  \n",
    "   - Open **QGIS** and add your classified raster layer.  \n",
    "\n",
    "2. **Apply Unique Colors to Each Class**  \n",
    "   - Open the **Layer Properties** → **Symbology** tab.  \n",
    "   - Set **Render Type** to **\"Unique Values\" (Paletted/Unique Values)**.  \n",
    "   - Click **Classify** to assign colors to each class.  \n",
    "\n",
    "3. **Suggest Improvements**  \n",
    "   - Consider adding other spectral indices/ other images/composites to the features.\n",
    "   - Consider adding training samples and/or change the classes (e.g. distinguish between croplands and grasslands)\n",
    "   - Consider changing the window size for smoothing the classification.\n",
    "\n",
    "#### **Hints:**  \n",
    "\n",
    "- Zoom into different regions to identify classification errors.  \n",
    "- Document misclassified areas and suggest possible fixes.  \n",
    "\n",
    "</div>\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rf-env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}