5 - Classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "76bcbe41-3588-4ccb-8fe6-3a15cbab839a",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<u><strong>Authors:</strong></u> <b>Alberto Vavassori</b> (alberto.vavassori@polimi.it), <b>Emanuele Capizzi</b> (emanuele.capizzi@polimi.it) - DICA - Politecnico di Milano - GIS GEOLab <br>\n",
    "Developed within the LCZ-ODC project, funded by the Italian Space Agency (agreement n. 2022-30-HH.0).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30ee2ff6-9718-4c22-98c0-59c78d59fbee",
   "metadata": {},
   "source": [
    "# Classification of PRISMA imagery into LCZ\n",
    "\n",
    "<a id='TOC_TOP'></a>\n",
    "Notebook structure:  <br>\n",
    " 1. [Import training samples](#sec1)  \n",
    " 2. [Rasterize training samples](#sec2)\n",
    " 3. [Import useful layers and prepare data for training](#sec3)\n",
    " 4. [Classifier training](#sec4)\n",
    " 5. [Accuracy assessment](#sec5)\n",
    " 6. [Classified image filtering and export](#sec6)\n",
    " 7. [Plot classified image](#sec7)\n",
    "<hr>\n",
    "\n",
    "This Notebook allows the user to perform the classification of PRISMA imagery in LCZ.\n",
    "The study area is the Metropolitan City of Milan.<br>\n",
    "The Notebook employs machine learning algorithms to perform the classification, exploiting built-in functions of the *[`scikit-learn`](https://ceholden.github.io/open-geo-tutorial/python/chapter_5_classification.html)* library. It provides options for different classification methods, including [Random Forest](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html#sklearn.ensemble.RandomForestClassifier) (RF), [AdaBoost](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.AdaBoostClassifier.html) (AB), [Gradient Boosting](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html#sklearn.ensemble.GradientBoostingClassifier) (GB), and [XGBoost](https://xgboost.readthedocs.io/en/stable/) (XGB)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bce83a1b-98c6-4b6d-9022-bc301d9918c6",
   "metadata": {},
   "source": [
    "### Additional resources\n",
    "\n",
    "* Hyperparameter tuning with [*scikit-learn*](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html): the [`GridSearchCV`](https://towardsdatascience.com/gridsearchcv-for-beginners-db48a90114ee) model\n",
    "* [Random Forest](https://stackabuse.com/random-forest-algorithm-with-python-and-scikit-learn/) and [parameters of Random Forest](https://www.analyticsvidhya.com/blog/2015/06/tuning-random-forest-model/)\n",
    "* [Gini index](https://medium.com/analytics-steps/understanding-the-gini-index-and-information-gain-in-decision-trees-ab4720518ba8) and [Entropy](https://towardsdatascience.com/entropy-how-decision-trees-make-decisions-2946b9c18c8) in the Random Forest: a [comparison](https://quantdare.com/decision-trees-gini-vs-entropy/)\n",
    "* [AdaBoost](https://www.analyticsvidhya.com/blog/2021/09/adaboost-algorithm-a-complete-guide-for-beginners/)\n",
    "* [Gradient Boosting](https://www.analyticsvidhya.com/blog/2021/09/gradient-boosting-algorithm-a-complete-guide-for-beginners/)\n",
    "* [XGBoost](https://towardsdatascience.com/beginners-guide-to-xgboost-for-classification-problems-50f75aac5390)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b790f9f3-2818-485c-8a4f-a9737f5e53e2",
   "metadata": {},
   "source": [
    "### Import libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ae32ed-b34c-4cf5-a725-26c8057a3e3e",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\">\n",
    "<span>&#9888;</span>\n",
    "Note: the Notebook relies on the <a href='https://gdal.org/' target='_blank'><em>gdal</em></a> Python library; make sure you have it installed in your environment.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbb9620a-737f-4d93-8e08-d3d58a9adc54",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import os\n",
    "import geopandas as gpd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "from rasterio.mask import mask\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from osgeo import gdal, ogr, gdalconst, gdal_array, osr\n",
    "from scipy.ndimage import median_filter\n",
    "import ipywidgets as widgets\n",
    "from shapely.geometry import box\n",
    "\n",
    "# ML libraries\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, classification_report\n",
    "from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier, GradientBoostingClassifier\n",
    "from sklearn.model_selection import GridSearchCV, train_test_split\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "#import xgboost as xgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "336a23bb-b002-4aed-b812-f2aadd78db5e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import functions and set auto-reload\n",
    "from functions import *\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39755524-7756-42da-9eb6-77604dc795bc",
   "metadata": {},
   "source": [
    "## 1. <a id='sec1'></a> Import training samples\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69c7d250-fc81-4ea7-9626-00cd1591923d",
   "metadata": {},
   "source": [
    "First, select the PRISMA image date:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3660241-06b9-4f23-8457-bb005df237c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "date_prisma_w = widgets.Dropdown(\n",
    "    options=['2023-02-09', '2023-03-22', '2023-04-08', '2023-06-17', '2023-07-10', '2023-08-08'],\n",
    "    value='2023-02-09',\n",
    "    description='PRISMA date:',\n",
    "    disabled=False,\n",
    "    layout={'width': 'max-content'},\n",
    "    style = {'description_width': 'initial'}\n",
    ")\n",
    "date_prisma_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86ca03b5-2f47-4e3e-82ea-0edfda63166d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sel_prisma_date = date_prisma_w.value\n",
    "selected_prisma_image = 'PRISMA_outputs/coregistered/PR_'+ sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "print(f\"The selected date is --> PRISMA: {sel_prisma_date}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "409883b8-b905-4f88-9c41-eb7823be850e",
   "metadata": {},
   "source": [
    "Set the legend that will be used henceforth for the plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe544258-f833-4805-86d1-53d548eebe4f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "legend = {\n",
    "    2: ['Compact mid-rise', '#D10000'],\n",
    "    3: ['Compact low-rise', '#CD0000'],\n",
    "    5: ['Open mid-rise', '#FF6600'],\n",
    "    6: ['Open low-rise', '#FF9955'],\n",
    "    8: ['Large low-rise', '#BCBCBC'],\n",
    "    101: ['Dense trees', '#006A00'],\n",
    "    102: ['Scattered trees', '#00AA00'],\n",
    "    104: ['Low plants', '#B9DB79'],\n",
    "    105: ['Bare rock or paved', '#545454'],\n",
    "    106: ['Bare soil or sand', '#FBF7AF'],\n",
    "    107: ['Water', '#6A6AFF']\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ec0e6ce-fb34-4546-a3d4-bdcd0aea5da8",
   "metadata": {},
   "source": [
    "The following function imports the training samples and computes the area of each LCZ class. Indeed, it is important to keep data balanced for the next classification steps (this is relevant expecially for urban classes, while natural classed usually are more easily classified). The function outputs a plot with the total area of each LCZ class as well as the path to the vector layer that will be used in the following."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1c0f09d-97ae-436d-87c7-3884cfd6c01b",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\">\n",
    "<span>&#9888;</span>\n",
    "<a id='warning'></a> \n",
    "\n",
    "**Note**:\n",
    "Training samples must be stored in `.gpkg` format as vector multi-polygons, and the file must contain a column `LCZ` with integer values corresponding to the LCZ class, as reported in the dictionary `legend`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164e6ba8-168f-4be9-921a-2b545212ef87",
   "metadata": {},
   "source": [
    "</div>\n",
    "<div class=\"alert alert-warning\" role=\"alert\">\n",
    "<span>&#9888;</span>\n",
    "<a id='warning'></a> Training samples are time-dependent, especially for the natural classes; they must be updated if there are changes (e.g. in the land cover).  The following function imports the training samples specific to the selected PRISMA acquisition date.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74570145-a26e-4e35-ac15-828bdca2d356",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vector_LCZ_path = training_area(sel_prisma_date, legend)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65806579-09b6-4547-ae2e-9ad721d2a15b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. <a id='sec2'></a>Rasterize training samples\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcdbd9bd-03b6-4570-aa1e-a1f5034eb66b",
   "metadata": {},
   "source": [
    "In order to perform the classification, it is necessary to convert the training set (provided in vector format as Geopackage) in raster format using the following functions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b250a06-d927-47d0-99a6-3e082291ea78",
   "metadata": {},
   "source": [
    "You need to set the following variables: (1) `raster_reference`, image to be classified (GeoTIFF) that will be used as reference image; metadata (including the spatial resolution) of the rasterized training set is obtained from the reference image; (2) `attribute`: field to be used for creating the raster; in this case, the `LCZ` field number of the vector file is used to set the value of each pixel in the rasterized training set; (3) `output`: output location and name of the complete rasterized training set (GeoTIFF); (4) `projection`: projection of the rasterized image; in this case, `EPSG:32632` is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76ae4f6d-d5b1-488f-86e7-6b2111e657ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_reference = 'PCs/PCs_'+ sel_prisma_date.replace('-', '') +'_30m.tif'\n",
    "output = './layers/training_samples/training_set_'+ sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "attribute = 'LCZ'\n",
    "projection = 32632"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5072c39-a192-4a62-a4d8-114099b1837a",
   "metadata": {},
   "source": [
    "The following function is meant to rasterize the training sample file. <br>\n",
    "In this case we used the image of 9th February 2023 as reference image for rasterization, since it covers the whole study area. This is done in order to get the complete rasterized training samples over the whole area (without clipped parts that may be generated using images that are not exacly located over the study area).<br>\n",
    "However, in principle, it is possible to pass to the function the specified `raster_reference`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d762bef2-b126-46d7-86a2-1714493ffe6e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "rasterized_result = rasterize_training('PCs/PCs_20230209_30m.tif', vector_LCZ_path, output, attribute, projection)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0dceadf-f5fa-47ea-b93b-d3f2003ed954",
   "metadata": {},
   "source": [
    "Optionally, you can use the following function to plot the rasterized training samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0a0544-b7b0-4817-a720-82fd59a9cbb3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#plot_raster_training(output, legend)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8503ea02-94ca-41f1-b3ca-d2abd97980e6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. <a id='sec3'></a>Import useful layers and prepare data for training\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b34e564-27a4-424e-a21a-3a728047983a",
   "metadata": {},
   "source": [
    "Before performing the classification, it is necessary to import the useful layers and to prepare them for training."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be8cc559-f796-4615-95a1-9bea36dc11f3",
   "metadata": {},
   "source": [
    "First, the bounding box for clipping the images to the correct extent must be obtained. The following function generates a Geopackage file with the extent of the selected PRISMA image footprint. The Geopackage is saved in the `study_area` folder which is automatically created if not yet existing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e69dbd4f-82d1-4364-9acf-d7246430c7ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bbox = prisma_bbox(selected_prisma_image, sel_prisma_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "521c1eb1-31dd-4c0e-9b0f-4ce85bacbcbf",
   "metadata": {},
   "source": [
    "Then, it is necessary to create a mask defining the extent of the area where PRISMA data is available.<br>\n",
    "The mask will be applied to all the layers used in the classification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9ddc4b2-e962-4d2b-ad23-580a7543e950",
   "metadata": {},
   "source": [
    "</div>\n",
    "<div class=\"alert alert-warning\" role=\"alert\">\n",
    "<span>&#9888;</span>\n",
    "<a id='warning'></a>\n",
    "\n",
    "**Note**:\n",
    "The mask contains `1` where PRISMA data is available and `0` at the edges, since `nan` is not supported by the `scikit-learn` library for classification purposes.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f09273ae-931e-4149-ab6f-a9531907d70c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mask_prisma = mask_prisma_image(selected_prisma_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7348135c-1aef-47a5-8ba8-00703684c46c",
   "metadata": {},
   "source": [
    "At this point, the user must import the layers that will be used for the classification. In this case, the following layers (GeoTIFF) are imported:\n",
    "* Urban Canopy Layers (*UCL*) including:\n",
    "    - *imperviousness layer* retrieved from [Copernicus Imperviousness Density 2018](https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness/status-maps/imperviousness-density-2018);\n",
    "    - *building heights* retrieved from [Geo-Topographic Database (DBGT) of Regione Lombardia](https://www.regione.lombardia.it/wps/portal/istituzionale/HP/DettaglioServizio/servizi-e-informazioni/Enti-e-Operatori/Territorio/sistema-informativo-territoriale-sit/database-topografico-regionale/database-topografico-regionale);\n",
    "    - *sky view factor* computed from DBGT using the homonymous built-in function of [SAGA-GIS](https://saga-gis.sourceforge.io/en/index.html);\n",
    "    - *building percentage* for each pixel, calculated from DBGT;\n",
    "    - *tree canopy height* retrieved from [ETH Global Sentinel-2 10m Canopy Height (2020)](https://gee-community-catalog.org/projects/canopy/);\n",
    "* PRISMA image Principal Components (*PCs*).\n",
    "\n",
    "The rasterized *training samples*, which will be used to train the classifier, have already been computed and imported in the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2d2990e-8afa-4133-8e22-196400ad6efa",
   "metadata": {},
   "source": [
    "Import the UCL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a82cbbe7-5b1a-4453-b7b8-a5acf205d138",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imperv_new_path = 'layers/ucp/imperviousness_' + sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "perc_build_new_path = 'layers/ucp/percentage_buildings_' + sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "svf_new_path = 'layers/ucp/SVF_' + sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "canopy_height_new_path = 'layers/ucp/canopy_height_ETH_' + sel_prisma_date.replace('-', '') + '_30m.tif'\n",
    "buildings_new_path = 'layers/ucp/buildings_' + sel_prisma_date.replace('-', '') + '_30m.tif'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f41325e-e1a2-4cea-b912-8dff4124076e",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following function clips the UCL to the extent of the PRISMA image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de2c978c-500a-4aa8-953b-b579f4f4f8bc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load the GeoPackage vector file representing the bounding box\n",
    "study_area = gpd.read_file(bbox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbb76eeb-fc2e-4736-83f3-62bfb0357c97",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "clip_image_study_area('layers/ucp/imperviousness_30m.tif', imperv_new_path, study_area)\n",
    "clip_image_study_area('layers/ucp/percentage_buildings_30m.tif', perc_build_new_path, study_area)\n",
    "clip_image_study_area('layers/ucp/SVF_30m.tif', svf_new_path, study_area)\n",
    "clip_image_study_area('layers/ucp/canopy_height_ETH_30m.tif', canopy_height_new_path, study_area)\n",
    "clip_image_study_area('layers/ucp/buildings_30m.tif', buildings_new_path, study_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f4dfc28-5245-4e18-b8e3-ea9f53c7d3bf",
   "metadata": {},
   "source": [
    "With the following function, the UCL are rescaled so that their values are between 0 and 1, and the PRISMA image mask is applied to each of them in order to have area matching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14f2141a-c137-4b42-a693-f4be79e18a38",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imperv = open_layer(imperv_new_path, mask_prisma)\n",
    "perc_build = open_layer(perc_build_new_path, mask_prisma)\n",
    "svf = open_layer(svf_new_path, mask_prisma)\n",
    "canopy_height = open_layer(canopy_height_new_path, mask_prisma)\n",
    "buildings = open_layer(buildings_new_path, mask_prisma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de82c228-ed6b-4ae6-bb07-3207c5929ef1",
   "metadata": {},
   "source": [
    "The following function displays the UCL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33e2e809-16c1-4dd1-a599-93c9bcf45f06",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_ucl(imperv, perc_build, svf, canopy_height, buildings)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb811579-98d1-4128-a1f2-0fcc77107f70",
   "metadata": {},
   "source": [
    "Make also sure that the training sample is clipped to the correct extent with the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daf11df4-0e55-493d-add1-e5176d681b14",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img, roi = clip_training_sample(raster_reference, output, sel_prisma_date, study_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bab018a-83b9-4d6d-a3dc-bdcf073cc992",
   "metadata": {},
   "source": [
    "Check the dimension of the layers used for the classification and expand if necessary. This is done to match layers dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "791a9bca-97c4-4a8b-961f-822c0b6d769e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imperv, perc_build, svf, canopy_height, buildings, roi = check_layers_dimension(imperv, perc_build, svf, canopy_height, buildings, roi, img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6228ba16-223b-48a1-9cb4-6bd0a24a7309",
   "metadata": {},
   "source": [
    "Stack all the layers (PRISMA PCs and the previously shown layers):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2690963-813f-4a4c-9761-1b7dbeb46a2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img = np.dstack((imperv, perc_build, svf, canopy_height, buildings, img))\n",
    "print(f\"The stacked array shape is --> {img.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00ce7943-0aaa-4a94-bd32-669a2c9b30f2",
   "metadata": {},
   "source": [
    "Set the classification labels, and the training samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7e58a0a-836f-4f3d-9d3b-757d850d31a9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "labels = np.unique(roi[roi > 0])\n",
    "print(f'The training data include {labels.size} classes: {labels}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b71be35-21ef-429a-9a45-ea68190097be",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "X = img[roi > 0, :] \n",
    "y = roi[roi > 0]\n",
    "print(f'X matrix size: {X.shape}')\n",
    "print(f'y array size: {y.shape}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fab372b5-2919-4ecd-8cd4-554b62bbd192",
   "metadata": {},
   "source": [
    "## 4. <a id='sec4'></a>Classifier training\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f3f2809-b2cb-4620-9221-38c3f73cd46f",
   "metadata": {},
   "source": [
    "This section performs the classification. First, training samples are split into training and testing sets. <br>\n",
    "Then, the user can select the classification method. Available methods are **Random Forest** (RF), XGBoost (XGB), AdaBoosting (AD) and GradientBoosting (BG). <br>\n",
    "Hyperparameters tuning is performed so that the combination of parameters yielding the best accuracy on the testing set is exploited."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2ec5601-15f7-4a3c-a9b2-981f183fb8d5",
   "metadata": {},
   "source": [
    "A remark about the classification methods:<br>\n",
    "RF: The Random Forest classifier is created using `RandomForestClassifier()`. It uses an ensemble of decision trees for classification. The classifier is configured with different hyperparameters, such as the number of estimators and maximum features, which are specified in the `param_grid` dictionary. The algorithm aims to find the best combination of hyperparameters using `GridSearchCV`.<br>\n",
    "XGB: The XGBoost classifier is created using `xgb.XGBClassifier()`. XGBoost is an optimized gradient boosting library that provides fast and accurate implementations of gradient boosting algorithms. Similar to the previous classifiers, hyperparameters such as the number of estimators and learning rate are tuned using `GridSearchCV`.<br>\n",
    "AB: The AdaBoost classifier is created using `AdaBoostClassifier()`. It combines multiple weak classifiers to create a strong classifier. Similar to Random Forest, the classifier is tuned using `GridSearchC` to optimize the number of estimators and learning rate.<br>\n",
    "GB: The Gradient Boosting classifier is created using `GradientBoostingClassifier()`. It builds an ensemble of decision trees in a sequential manner. The classifier is tuned using `GridSearchCV` to optimize hyperparameters such as learning rate, maximum features, criterion, and number of estimators."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e455617-c123-4ea2-8372-17afa8af5cf9",
   "metadata": {},
   "source": [
    "First, it's needed to split train and test data (80% for training and 20% for testing):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "225bd163-c20d-469c-aa10-b4a7f746ac40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Split data into training and test sets\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6ff1ae1-46ea-4a80-b3fc-840f1dd9e1ee",
   "metadata": {},
   "source": [
    "Select the classification method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3abfefde-5b6b-4653-a45f-3a4e4b4a8404",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "classification_method_w = widgets.RadioButtons(\n",
    "    options=['RandomForest', 'XGBoost', 'AdaBoost', 'GradientBoost'],\n",
    "    value='RandomForest', # Defaults to 'pineapple'\n",
    "    layout={'width': 'max-content'}, # If the items' names are long\n",
    "    description='Classifier:',\n",
    "    disabled=False\n",
    ")\n",
    "classification_method_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e01b3bc6-73d8-477f-a275-31738e8198bf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "classification_method = classification_method_w.value\n",
    "print(f'Selected classification method: {classification_method}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2165150-29e0-46b9-a2b6-c0877420417d",
   "metadata": {},
   "source": [
    "The following function performs hyperparameter tuning for the selected classification algorithm, and returns the best hyperparameter. The selection is done by computing the accuracy score:\n",
    "\n",
    "$ {Accuracy Score} = {(TP+TN)\\over(TP+FN+TN+FP)}$\n",
    "\n",
    "where $TP$ and $TN$ are true positives/negatives, $FP$ abd $FN$ are false positives/negatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb6c8e1b-591a-4784-8e1d-e057ca963b2e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "best_params = parameter_tuning(classification_method, X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3644a95-7e88-4a26-aa44-71da40b2ed88",
   "metadata": {},
   "source": [
    "The following function creates classifier objects based on the selected classification method. The classifier is created with the hyperparameters yielding the best performance based on the accuracy score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef3bcb3c-7d7d-4137-9a95-50d9c7e04290",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "y_pred, clc = classification(classification_method, best_params, X_train, y_train, X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "567848c8-0dd4-48a9-bde2-0b7bc1ff3eaa",
   "metadata": {},
   "source": [
    "## 5. <a id='sec5'></a> Accuracy assessment\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9080091f-97e4-4320-9d97-bb45aa665529",
   "metadata": {},
   "source": [
    "Calculate accuracy, confusion matrix and statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28028e4b-6f2a-473c-85c8-fef5efba3d23",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculate the accuracy of the best model on the test set\n",
    "accuracy = accuracy_score(y_test, y_pred)\n",
    "print(f\"Accuracy: {accuracy:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b709ed64-ab08-4a25-b32e-0f3ae827aa7b",
   "metadata": {},
   "source": [
    "https://github.com/jrkreiger/random-forest-trees/blob/master/random-forest-for-trees.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f481d086-6780-47f8-9601-b281e3ed95bb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print_metrics(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3c22da8-6147-42b4-863e-187bcb09851d",
   "metadata": {},
   "source": [
    "## 6. <a id='sec6'></a> Classified image filtering and export\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a83fe3df-9ca2-4d95-9070-db8053dcf240",
   "metadata": {},
   "source": [
    "The classified map is exported with the following function. The path is directly set within the function depending on the name of the original PRISMA image.<br>\n",
    "The function also applies a post-processing to the classified image using the [scipy median filter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.median_filter.html). This is done to remove isolated pixels obtained from the classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f94cbe79-8c79-425a-bfb9-8d0d5c306714",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "export_classified_map(img, clc, X, mask_prisma, selected_prisma_image, classification_method, sel_prisma_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f1a5beb-0e28-409e-86e8-79afe2da38ee",
   "metadata": {},
   "source": [
    "## 7. <a id='sec7'></a>Plot classified image\n",
    "[Back to top](#TOC_TOP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4822ef66-4a67-4fdc-a72d-55703ca7c23c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#def plot_classified_image(classification_method, sel_prisma_date):\n",
    "    \n",
    "#    # Open the image file using rasterio\n",
    "#    folder_path = 'classified_images'\n",
    "#    with rasterio.open(folder_path + '/classified_' + classification_method + '_' + sel_prisma_date.replace('-', '') + '_medianfilter_30m.tif') as dataset:\n",
    "#        # Read the image data\n",
    "#        image_data = dataset.read()\n",
    "    \n",
    "#    # Plot the raster data with a custom figure size\n",
    "#    fig = plt.figure(figsize=(14, 14))\n",
    "#    ax = fig.add_subplot(1, 1, 1)\n",
    "#    plt.title(f\"Classified LCZ image with {classification_method}\")\n",
    "#    show(image_data, cmap=cmap, ax=ax, interpolation='none',norm=norm);\n",
    "\n",
    "#    cmap = colors.ListedColormap([colors_dict[key] for key in colors_dict])\n",
    "\n",
    "#    # Add a legend with the correct class colors\n",
    "#    labels = list(colors_dict.keys())\n",
    "#    handles = [plt.Rectangle((0, 0), 1, 1, color=colors_dict[label]) for label in labels]\n",
    "#    legend = plt.legend(handles, labels, title='LCZ Classes', bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "#    plt.setp(legend.get_title(), fontsize='12')  # Adjust the font size of the legend title\n",
    "\n",
    "#    # Adjust the plot layout to accommodate the legend outside\n",
    "#    plt.tight_layout(rect=[0, 0, 0.85, 1])  # Increase the left margin to make space for the legend\n",
    "\n",
    "#    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0dfed4f-5676-43e7-b54d-f314da561e58",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a186af8a-41d0-4758-9c13-165baadcf1b4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}