{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Random Forest Classification\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", "\n", "Random Forest (RF), an improved implementation of Decision Trees (DT), is an ensemble-learning algorithm that combines multiple classifications of the same data to produce higher classification accuracies than other forms of DT. RF works by fitting many DT-based classifications to a data set, and then uses a rule-based approach to combine the predictions from all the trees. During this process, individual trees are grown from differing subsets of training data using a process called “bagging”. Bagging involves the random subsampling (with replacement) of the original data for growing each tree. Generally, for each tree grown, two thirds of the training data are used to grow the tree, while the remaining one third are left unused (out-of-bag, or OOB) for later error assessment. A classification is then fit to each bootstrap sample; however, at each node (split), only a small number of randomly selected predictor variables are used in the binary partitioning. The splitting process continues until further subdivision no longer reduces the Gini index. Each tree contributes to the assignment of the most frequent class to the input data with a single vote. The predicted class of an observation is calculated by the majority vote for that observation, with ties split randomly.\n", "\n", "\n", "> **The random forest combines hundreds or thousands of decision trees, trains each one on a slightly different set of the observations, splitting nodes in each tree considering a limited number of the 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. Scikit-learn is an amazing machine learning library that provides easy and consistent interfaces to many of the most popular machine learning algorithms. It is built on top of the pre-existing scientific Python libraries, including NumPy, SciPy, and matplotlib, which makes it very easy to incorporate into your workflow. The number of available methods for accomplishing any task contained within the library is its real strength.\n", "\n", "\n", "
\n", " \"Image\n", "
\n", "\n", "\n", "\n", "
\n", " \"Image\n", "
Namur, 2020 (NDVI & monthly composites S1 backscattering VV)
\n", "
\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": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "All libraries successfully imported!\n", "Scikit-learn : 1.2.0\n" ] } ], "source": [ "import glob, os, time, math\n", "import numpy as np\n", "import pandas as pd\n", "import geopandas as gpd\n", "import rasterio\n", "import rasterio.plot\n", "from rasterio import features\n", "\n", "import sklearn\n", "from sklearn.ensemble import RandomForestClassifier\n", "\n", "from pathlib import Path\n", "from IPython.display import display\n", "\n", "print('All libraries successfully imported!')\n", "print(f'Scikit-learn : {sklearn.__version__}')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Set directory**" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Classification path is set to : /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/CLASSIF/\n" ] } ], "source": [ "computer_path = 'X:/'\n", "grp_nb = '99'\n", "\n", "data_path = f'{computer_path}data/' # Directory with data shared by the assistant\n", "work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/' # Directory for all work files\n", "\n", "# Input directories\n", "in_situ_path = f'{work_path}C_IN_SITU_CAL_VAL/'\n", "s2_path = f'{work_path}3_L2A_MASKED/'\n", "ndvi_path = f'{work_path}NDVI/'\n", "s1_path = f'{data_path}S1_GRD/'\n", "lut_path = f'{data_path}LUT/'\n", "\n", "# Output directory\n", "classif_path = f'{work_path}CLASSIF/'\n", "\n", "Path(classif_path).mkdir(parents=True, exist_ok=True)\n", "\n", "print(f'Classification path is set to : {classif_path}')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Set parameters**" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [], "source": [ "site = 'NAMUR'\n", "year = '2020'\n", "\n", "no_data = -999\n", "\n", "ws = 3 # Must be a odd number 3x3, 5x5, 7x7, ... (filtering post classification)\n", "\n", "# Field used for classification\n", "field_classif_code = 'grp_1_nb'\n", "field_classif_name = 'grp_1'\n", "\n", "# Field used for reclassification\n", "field_reclassif_code = 'grp_A_nb'\n", "field_reclassif_name = 'grp_A'\n", "\n", "# Group of features used in classification\n", "feat_nb = 2\n", "\n", "if feat_nb == 1:\n", " feat_name = ['NDVI']\n", "elif feat_nb == 2:\n", " feat_name = ['NDVI','S1_monthly_mean']" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Set filenames**" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Input files" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [], "source": [ "in_situ_cal_shp = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_CAL.shp'\n", "\n", "s4s_lut_xlsx = f'{lut_path}crop_dictionary.xlsx'" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Output files" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [], "source": [ "in_situ_cal_tif = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_CAL.tif'\n", "\n", "classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}.tif'\n", "reclassif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}_reclassify_{field_reclassif_name}.tif'\n", "reclassif_filter_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}_reclassify_{field_reclassif_name}_filter_ws_{ws}x{ws}.tif'" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Prepare classification features associated to *in situ* data" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Rasterize *in situ* data calibration shapefile" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Raster template file : /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/3_L2A_MASKED/T31UFS_20200116T105309_B02_10m_ROI_SCL.tif\n", "The CRS of in situ data is : 32631\n", "The CRS of raster template is : 32631\n", "CRS are the same\n", "Rasterize starts : /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/C_IN_SITU_CAL_VAL/NAMUR_2020_IN_SITU_ROI_CAL.shp\n", "Rasterize is done : /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/C_IN_SITU_CAL_VAL/NAMUR_2020_IN_SITU_ROI_CAL.tif\n" ] } ], "source": [ "# Open the calibration polygons with GeoPandas\n", "in_situ_gdf = gpd.read_file(in_situ_cal_shp)\n", "\n", "# Open the raster file you want to use as a template for rasterize\n", "img_temp_tif = glob.glob(f'{s2_path}*.tif')[0]\n", "\n", "print(f'Raster template file : {img_temp_tif}')\n", "\n", "src = rasterio.open(img_temp_tif, \"r\")\n", "\n", "# Update metadata\n", "\n", "out_meta = src.meta\n", "out_meta.update(nodata=no_data)\n", "\n", "crs_shp = str(in_situ_gdf.crs).split(\":\",1)[1]\n", "crs_tif = str(src.crs).split(\":\",1)[1]\n", "\n", "print(f'The CRS of in situ data is : {crs_shp}')\n", "print(f'The CRS of raster template is : {crs_tif}')\n", "\n", "if crs_shp == crs_tif:\n", " print(\"CRS are the same\")\n", "\n", " print(f'Rasterize starts : {in_situ_cal_shp}')\n", "\n", " # Burn the features into the raster and write it out\n", "\n", " dst = rasterio.open(in_situ_cal_tif, 'w+', **out_meta)\n", " dst_arr = dst.read(1)\n", "\n", " # This is where we create a generator of geom, value pairs to use in rasterizing\n", "\n", " geom_col = in_situ_gdf.geometry\n", " code_col = in_situ_gdf[field_classif_code].astype(int)\n", "\n", " shapes = ((geom,value) for geom, value in zip(geom_col, code_col))\n", "\n", " in_situ_arr = features.rasterize(shapes=shapes,\n", " fill=no_data,\n", " out=dst_arr,\n", " transform=dst.transform)\n", "\n", " dst.write_band(1, in_situ_arr)\n", "\n", " print(f'Rasterize is done : {in_situ_cal_tif}')\n", "\n", " # Close rasterio objects\n", " src.close()\n", " dst.close()\n", "\n", "else:\n", " print('CRS are different --> repoject in-situ data shapefile with \"to_crs\"')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 List all the classification features" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Create an empty list to append all feature rasters one by one\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "list_src_arr = []" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "1 NDVI image per month" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Shape of features : (570, 986)\n", "Number of features : 12\n" ] } ], "source": [ "if 'NDVI' in feat_name:\n", "\n", " list_im = sorted(glob.glob(f'{ndvi_path}*.tif'))\n", "\n", " for im_file in list_im:\n", "\n", " src = rasterio.open(im_file, \"r\")\n", " im = src.read(1)\n", " list_src_arr.append(im)\n", " src.close()\n", " \n", " print(f'Shape of features : {im.shape}')\n", " print(f'Number of features : {len(list_src_arr)}')\n", "\n", "else:\n", " print(\"No NDVI in the set of features\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "S1 monthly mean composite (obtained with Google Earth Engine)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Shape of features : (570, 986)\n", "Number of features : 24\n" ] } ], "source": [ "if 'S1_monthly_mean' in feat_name:\n", "\n", " s1_montlhy_mean_tif = f'{s1_path}monthly_mean_{site}_{year}.tif'\n", "\n", " src = rasterio.open(s1_montlhy_mean_tif, \"r\")\n", " im = src.read()\n", " src.close()\n", "\n", " for i in range(len(im)):\n", " band = im[i]\n", " list_src_arr.append(band)\n", "\n", " print(f'Shape of features : {band.shape}')\n", " print(f'Number of features : {len(list_src_arr)}')\n", "\n", "else:\n", " print(\"No S1 monthly mean in the set of features\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Merge all the 2D matrices from the list into one 3D matrix\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(570, 986, 24)\n", "There are 24 features\n", "The features type is : float32\n" ] } ], "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}')\n", "\n", "#feat_arr_1 = np.stack(list_src_arr, axis=0)\n", "#print(feat_arr_1.shape)" ] }, { "attachments": {}, "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": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "We have 34909 samples (= calibration pixels)\n" ] } ], "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)')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "What are our classification labels?" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The training data include 19 classes: [ 3 21 22 69 81 84 121 1111 1121 1152 1171 1192 1435 1511\n", " 1771 1811 1911 1923 9212]\n" ] } ], "source": [ "labels = np.unique(cal_arr[cal_arr != no_data])\n", "\n", "print(f'The training data include {labels.size} classes: {labels}')" ] }, { "attachments": {}, "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": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our X matrix is sized: (34909, 24)\n", "Our y array is sized: (34909,)\n" ] } ], "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}')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Train the Random Forest" ] }, { "attachments": {}, "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 web page to find the usage of RandomForestClassifier from scikit-learn." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Random Forest training : 00:00:13.95\n" ] } ], "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))" ] }, { "attachments": {}, "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": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our OOB prediction of accuracy is: 99.49%\n" ] } ], "source": [ "print(f'Our OOB prediction of accuracy is: {round(rf.oob_score_ * 100,2)}%')" ] }, { "attachments": {}, "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": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
feat_bandfeat_nameGini
67NDVI July0.0994
78NDVI August0.0885
45NDVI May0.0871
89NDVI September0.0784
34NDVI April0.0722
23NDVI March0.0702
01NDVI January0.0555
1718S1 mean composite - June0.0479
910NDVI October0.0476
1819S1 mean composite - July0.0447
1011NDVI November0.0427
1112NDVI December0.0320
1920S1 mean composite - August0.0318
1617S1 mean composite - May0.0305
1516S1 mean composite - April0.0302
56NDVI June0.0285
12NDVI February0.0213
2021S1 mean composite - September0.0194
1213S1 mean composite - January0.0173
1314S1 mean composite - February0.0128
2324S1 mean composite - December0.0118
1415S1 mean composite - March0.0113
2122S1 mean composite - October0.0101
2223S1 mean composite - November0.0088
\n", "
" ], "text/plain": [ " feat_band feat_name Gini\n", "6 7 NDVI July 0.0994\n", "7 8 NDVI August 0.0885\n", "4 5 NDVI May 0.0871\n", "8 9 NDVI September 0.0784\n", "3 4 NDVI April 0.0722\n", "2 3 NDVI March 0.0702\n", "0 1 NDVI January 0.0555\n", "17 18 S1 mean composite - June 0.0479\n", "9 10 NDVI October 0.0476\n", "18 19 S1 mean composite - July 0.0447\n", "10 11 NDVI November 0.0427\n", "11 12 NDVI December 0.0320\n", "19 20 S1 mean composite - August 0.0318\n", "16 17 S1 mean composite - May 0.0305\n", "15 16 S1 mean composite - April 0.0302\n", "5 6 NDVI June 0.0285\n", "1 2 NDVI February 0.0213\n", "20 21 S1 mean composite - September 0.0194\n", "12 13 S1 mean composite - January 0.0173\n", "13 14 S1 mean composite - February 0.0128\n", "23 24 S1 mean composite - December 0.0118\n", "14 15 S1 mean composite - March 0.0113\n", "21 22 S1 mean composite - October 0.0101\n", "22 23 S1 mean composite - November 0.0088" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "feat_band_list = []\n", "gini_list = []\n", "\n", "for band_nb, imp in enumerate(rf.feature_importances_, start=1):\n", "\n", " feat_band_list.append(band_nb)\n", " gini_list.append(imp)\n", "\n", "gini_dict = {'feat_band':feat_band_list,'Gini':gini_list} \n", "\n", "gini_df = pd.DataFrame(gini_dict).round(4)\n", "\n", "feat_name_df = pd.read_excel(f'{data_path}feature_name.xlsx')\n", "\n", "gini_df = feat_name_df.merge(gini_df, on='feat_band').sort_values(by='Gini', ascending=False)\n", "\n", "gini_df" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at a crosstabulation to see the class confusion" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
predict32122698184121111111211152117111921435151117711811191119239212All
truth
359710000000000000000005971
2101551000000000000000001551
22001240000000000000000124
6900022960000000000000002296
8100001338000000000000001338
84000001830000000000000183
121000000566000000000000566
111100000007403000000000007403
112100000000213200000000002132
115200000000016630000000001663
1171000000000053400000000534
119200000000000147700000001477
143500000000000011760000001176
151100000000000001577000001577
177100000000000000245400002454
181100000000000000018500001850
1911000000000000000040300403
192300000000000000000205902059
9212000000000000000000152152
All597115511242296133818356674032132166353414771176157724541850403205915234909
\n", "
" ], "text/plain": [ "predict 3 21 22 69 81 84 121 1111 1121 1152 1171 1192 \\\n", "truth \n", "3 5971 0 0 0 0 0 0 0 0 0 0 0 \n", "21 0 1551 0 0 0 0 0 0 0 0 0 0 \n", "22 0 0 124 0 0 0 0 0 0 0 0 0 \n", "69 0 0 0 2296 0 0 0 0 0 0 0 0 \n", "81 0 0 0 0 1338 0 0 0 0 0 0 0 \n", "84 0 0 0 0 0 183 0 0 0 0 0 0 \n", "121 0 0 0 0 0 0 566 0 0 0 0 0 \n", "1111 0 0 0 0 0 0 0 7403 0 0 0 0 \n", "1121 0 0 0 0 0 0 0 0 2132 0 0 0 \n", "1152 0 0 0 0 0 0 0 0 0 1663 0 0 \n", "1171 0 0 0 0 0 0 0 0 0 0 534 0 \n", "1192 0 0 0 0 0 0 0 0 0 0 0 1477 \n", "1435 0 0 0 0 0 0 0 0 0 0 0 0 \n", "1511 0 0 0 0 0 0 0 0 0 0 0 0 \n", "1771 0 0 0 0 0 0 0 0 0 0 0 0 \n", "1811 0 0 0 0 0 0 0 0 0 0 0 0 \n", "1911 0 0 0 0 0 0 0 0 0 0 0 0 \n", "1923 0 0 0 0 0 0 0 0 0 0 0 0 \n", "9212 0 0 0 0 0 0 0 0 0 0 0 0 \n", "All 5971 1551 124 2296 1338 183 566 7403 2132 1663 534 1477 \n", "\n", "predict 1435 1511 1771 1811 1911 1923 9212 All \n", "truth \n", "3 0 0 0 0 0 0 0 5971 \n", "21 0 0 0 0 0 0 0 1551 \n", "22 0 0 0 0 0 0 0 124 \n", "69 0 0 0 0 0 0 0 2296 \n", "81 0 0 0 0 0 0 0 1338 \n", "84 0 0 0 0 0 0 0 183 \n", "121 0 0 0 0 0 0 0 566 \n", "1111 0 0 0 0 0 0 0 7403 \n", "1121 0 0 0 0 0 0 0 2132 \n", "1152 0 0 0 0 0 0 0 1663 \n", "1171 0 0 0 0 0 0 0 534 \n", "1192 0 0 0 0 0 0 0 1477 \n", "1435 1176 0 0 0 0 0 0 1176 \n", "1511 0 1577 0 0 0 0 0 1577 \n", "1771 0 0 2454 0 0 0 0 2454 \n", "1811 0 0 0 1850 0 0 0 1850 \n", "1911 0 0 0 0 403 0 0 403 \n", "1923 0 0 0 0 0 2059 0 2059 \n", "9212 0 0 0 0 0 0 152 152 \n", "All 1176 1577 2454 1850 403 2059 152 34909 " ] }, "metadata": {}, "output_type": "display_data" } ], "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" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Unbelievable? I highly doubt the real confusion matrix will be 100% accuracy. What is likely going on is that we used a large number of trees within a machine learning algorithm to best figure out the pattern in our training data. 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." ] }, { "attachments": {}, "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": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 6.60263658e-01 4.29129243e-01 6.68883979e-01 ... -1.23333731e+01\n", " -1.23078318e+01 -1.07777815e+01]\n", " [ 5.07605612e-01 4.32377875e-01 5.90189874e-01 ... -9.39306545e+00\n", " -1.21026096e+01 -1.09584541e+01]\n", " [ 4.00820792e-01 3.75474095e-01 4.71971065e-01 ... -8.66642475e+00\n", " -1.05296297e+01 -1.01824131e+01]\n", " ...\n", " [-9.99000000e+02 -9.99000000e+02 -9.99000000e+02 ... -7.79056311e+00\n", " -8.99753094e+00 -9.10786438e+00]\n", " [-9.99000000e+02 -9.99000000e+02 -9.99000000e+02 ... -8.75600719e+00\n", " -8.76734066e+00 -9.15415192e+00]\n", " [-9.99000000e+02 -9.99000000e+02 -9.99000000e+02 ... -1.01634216e+01\n", " -9.73724937e+00 -9.64400387e+00]]\n", "Reshaped from (570, 986, 24) to (562020, 24)\n", "Random Forest training : 00:00:09.11\n", "[[ 3 21 81 ... 3 3 3]\n", " [ 3 21 21 ... 3 3 3]\n", " [ 3 3 21 ... 3 3 1121]\n", " ...\n", " [ 69 69 69 ... 1111 1111 1111]\n", " [ 69 69 69 ... 1111 1111 1111]\n", " [ 81 81 81 ... 81 9212 81]]\n" ] } ], "source": [ "# Take our full image and reshape into long 2d array (nrow * ncol, nband) for classification\n", "\n", "img = feat_arr\n", "\n", "img = np.nan_to_num(img, nan=no_data)\n", "\n", "new_shape = (img.shape[0] * img.shape[1], img.shape[2])\n", "\n", "img_as_array = img[:, :, :].reshape(new_shape)\n", "\n", "print(img_as_array)\n", "\n", "print(f'Reshaped from {img.shape} to {img_as_array.shape}')\n", "\n", "start_classification = time.time()\n", "\n", "# Now predict for each pixel\n", "class_prediction = rf.predict(img_as_array)\n", "\n", "# Reshape our classification map\n", "class_prediction = class_prediction.reshape(img[:, :, 0].shape)\n", "\n", "end_classification = time.time()\n", "\n", "hours, rem = divmod(end_classification-start_classification, 3600)\n", "minutes, seconds = divmod(rem, 60)\n", "print(\"Random Forest training : {:0>2}:{:0>2}:{:05.2f}\".format(int(hours),int(minutes),seconds))\n", "\n", "print(class_prediction)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Reclassify classification" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 4.1 Open LUT and sort values" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
grp_1_nbgrp_1grp_A_nbgrp_A
00Remove0Remove
710Remove0Remove
750Remove0Remove
820Remove0Remove
860Remove0Remove
\n", "
" ], "text/plain": [ " grp_1_nb grp_1 grp_A_nb grp_A\n", "0 0 Remove 0 Remove\n", "71 0 Remove 0 Remove\n", "75 0 Remove 0 Remove\n", "82 0 Remove 0 Remove\n", "86 0 Remove 0 Remove" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lut_df = pd.read_excel(s4s_lut_xlsx)\n", "\n", "lut_df = lut_df.sort_values(by=field_classif_code, ascending=True)\n", "\n", "display(lut_df[[field_classif_code, field_classif_name, field_reclassif_code, field_reclassif_name]].head())" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 4.2 Reclassify prediction" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Classification : \n", " [[ 3 21 81 ... 3 3 3]\n", " [ 3 21 21 ... 3 3 3]\n", " [ 3 3 21 ... 3 3 1121]\n", " ...\n", " [ 69 69 69 ... 1111 1111 1111]\n", " [ 69 69 69 ... 1111 1111 1111]\n", " [ 81 81 81 ... 81 9212 81]]\n", "Re-classification : \n", " [[ 3 21 8 ... 3 3 3]\n", " [ 3 21 21 ... 3 3 3]\n", " [ 3 3 21 ... 3 3 112]\n", " ...\n", " [ 6 6 6 ... 111 111 111]\n", " [ 6 6 6 ... 111 111 111]\n", " [ 8 8 8 ... 8 9 8]]\n" ] } ], "source": [ "reclass_prediction = np.copy(class_prediction)\n", "\n", "for i, row in lut_df.iterrows():\n", " \n", " old_class = row[field_classif_code]\n", " new_class = row[field_reclassif_code]\n", "\n", " #print(f'{old_class} --> {new_class}')\n", "\n", " #array[np.where(array == old_class)] = new_class\n", "\n", " reclass_prediction[reclass_prediction == old_class] = new_class\n", "\n", "print(f'Classification : \\n {class_prediction}')\n", "print(f'Re-classification : \\n {reclass_prediction}')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Filter classification with moving window" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Re-classification : \n", " [[ 3 21 8 ... 3 3 3]\n", " [ 3 21 21 ... 3 3 3]\n", " [ 3 3 21 ... 3 3 112]\n", " ...\n", " [ 6 6 6 ... 111 111 111]\n", " [ 6 6 6 ... 111 111 111]\n", " [ 8 8 8 ... 8 9 8]]\n", "Re-classification with filter : \n", " [[[ 3 21 8 ... 3 3 3]\n", " [ 3 3 21 ... 3 3 3]\n", " [ 3 3 3 ... 3 3 112]\n", " ...\n", " [ 6 6 6 ... 111 111 111]\n", " [ 6 6 6 ... 111 111 111]\n", " [ 8 8 8 ... 8 8 8]]]\n" ] } ], "source": [ "sizey = reclass_prediction.shape[0]\n", "sizex = reclass_prediction.shape[1]\n", "\n", "rad_wind = math.floor(ws/2)\n", "\n", "X = np.pad(reclass_prediction, ((rad_wind,rad_wind),(rad_wind,rad_wind)), 'edge')\n", "\n", "majority = np.empty((sizey,sizex), dtype='int16')\n", "\n", "for i in range(sizey):\n", " for j in range(sizex):\n", " window = X[i:i+ws , j:j+ws]\n", " window = window.flatten()\n", " counts = np.bincount(window)\n", " maj = np.argmax(counts)\n", " majority[i,j]= maj\n", "\n", "\n", "majority = majority.reshape((1,sizey,sizex))\n", "\n", "print(f'Re-classification : \\n {reclass_prediction}')\n", "print(f'Re-classification with filter : \\n {majority}')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Write classification products into GeoTIFF files" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Open template image to get metadata" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -10000.0, 'width': 986, 'height': 570, 'count': 1, 'crs': CRS.from_epsg(32631), 'transform': Affine(10.0, 0.0, 627260.0,\n", " 0.0, -10.0, 5596180.0), 'blockysize': 4, 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with rasterio.open(img_temp_tif) as src:\n", " profile = src.profile\n", "\n", "profile" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Write classification**" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "with rasterio.open(classif_tif, \"w\", **profile) as dst:\n", " dst.write(class_prediction, 1)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Write re-classification**" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "with rasterio.open(reclassif_tif, \"w\", **profile) as dst:\n", " dst.write(reclass_prediction, 1)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Write re-classification with moving window filtering**" ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [], "source": [ "with rasterio.open(reclassif_filter_tif, \"w\", **profile) as dst:\n", " dst.write(majority)" ] } ], "metadata": { "interpreter": { "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" }, "kernelspec": { "display_name": "Python 3.6.10 64-bit ('python_3_6': conda)", "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.8" }, "orig_nbformat": 2 }, "nbformat": 4, "nbformat_minor": 2 }