CLA_3_tree_based.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3 Tree based classification\n",
    "-------------------\n",
    "\n",
    "**Decision Tree** is the most widely used **additive** model for **supervised learning** that combines simple classification nodes for larger scale tasks of both **classification** and **regression**. In most cases, the basic element in the **decision tree** is a binary node that split the data further into two parts. In a tree, many of the nodes would split the data incrementally along the tree, and bring the final classes at a node (usually referred as the leaf) where the data would not be split any more. There are also cases that the data can be split into more than two classes as opposed to binary classes, such as the [**CHAID**](https://en.wikipedia.org/wiki/Chi-square_automatic_interaction_detection), which will not be covered in this session, but please do feel free to explore further through the link.\n",
    "\n",
    "The logic of **additivity** can also be applied to **decision trees** to form an **ensemble model** call the **random forest** for **ensemble learning**. In **ensemble learning**, each tree is applied to the data and trained in same manner. The results of all the trees would then be aggregated, e.g. through averaging, to produce the final result. The **ensembled trees** will be covered later in this course, but not in this session. \n",
    "\n",
    "In this session, we will be focusing on **decision tree** and its application to the real GIS and Remote Sensing data that we have used in the previous sessions. The structure will be:\n",
    "\n",
    "- 3.0 Decision tree based classification of real world GIS and Remote Sensing datasets\n",
    "- 3.1 Further inspection: features and model configurations\n",
    "\n",
    "<font color='blue'> ***Throughout the notebook, you will encounter Questions, highlighted in blue. Please feel free to discuss them with your classmates and advisors***</font>\n",
    "\n",
    "<font color='green'> ***Throughout the notebook, you will also encounter shorter and longer Exercises, highlighted in green.***</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.0 Decision tree based classification of real world GIS and Remote Sensing datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we will use the same dataset of the satellite images along with the labelled LULC types of the Netherlands. The labelled LULC types in the AOIs will be split into **training** and **test** datasets. The trained model will be applied to the larger area covering almost the entire Netherlands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline  # In order to plot figures inline in Jupyter Notebook, we need to run this. But please ignore this in Colab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Before reading the data we need to first clone the data on Github to our Colab workspace\n",
    "!git clone https://github.com/jonwangio/uu_ml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same real world dataset as you have already encountered in the previous sessions.\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "% pip install geopandas\n",
    "import geopandas as gpd\n",
    "% pip install rasterio\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "from sklearn.svm import SVC\n",
    "\n",
    "# The data contains sample LULC areas around dutch provinces North Holland and Utrecht.\n",
    "aoi = gpd.read_file('uu_ml/data/aoi_NL_5_classes.shp')\n",
    "\n",
    "print (aoi.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Below is a visualization of the sample LULC areas superimposed on a satellite image of part of the Netherlands\n",
    "file_location = 'uu_ml/data/b5_2015.TIF'\n",
    "b5_2020 = rasterio.open(file_location, nodata=0)\n",
    "\n",
    "# We also prepare the color codes for visualization\n",
    "colors = [(257, 71, 27), (98, 93, 78), (14, 79, 58), (26, 0, 50), (75, 90, 85), (347, 72, 60), (246, 79, 60)]\n",
    "cols = []\n",
    "for col in colors:\n",
    "    pal = sns.light_palette(col, input=\"husl\", n_colors=4)\n",
    "    for rgb in pal[1:]:\n",
    "        cols.append(rgb)\n",
    "\n",
    "# A preview of color codes. Please delete the triple quotation marks to run the code.\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(20, 5))\n",
    "for i, c in enumerate(cols):\n",
    "    ax.add_artist(plt.Circle((i, 0), 0.4, color=c))\n",
    "    plt.text(i, -1, i, horizontalalignment='center')\n",
    "    ax.set_axis_off()\n",
    "    ax.set_aspect(1)\n",
    "    ax.autoscale()\n",
    "    plt.xlim(-1.25,43.25)\n",
    "    plt.ylim(-1,1)\n",
    "\n",
    "\n",
    "# Assign color codes to LULC types \n",
    "symbology = {'Agriculture': cols[5],\n",
    "             'Clear water': cols[20],\n",
    "             'Deciduous forest': cols[13],\n",
    "             'Residential': cols[17],\n",
    "             'Sand': cols[11]}\n",
    "\n",
    "# Visualize\n",
    "fig,ax = plt.subplots(1,1, figsize=(10,10))\n",
    "show(b5_2020, ax=ax, cmap='gray', alpha=0.25)\n",
    "aoi.plot(ax=ax, column='land_cover', legend=True, color=aoi['land_cover'].map(symbology))\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_points = [Line2D([0], [0], marker=\"o\", linestyle=\"none\", markersize=5, color=color) for color in symbology.values()]\n",
    "leg_points = ax.legend(custom_points, symbology.keys(), loc='upper right', frameon=False)\n",
    "ax.add_artist(leg_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you may recall, we will prepare the **training** and **test** datasets from the sample LULC types in AOIs. We will use the two bands of satellite images as the inputs and the manually delineated LULC types in the AOIs as the output labels. Did you still remember how we stacked the data together and extract the input band pixels along with labels within the AOIs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary modules and functionalities\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib.pyplot import imread\n",
    "from sklearn.svm import SVC\n",
    "\n",
    "# As before, we start with loading and stacking the image bands.\n",
    "# Again, we start with 2 bands that appeared to be useful for the small dataset.\n",
    "file_list = ['uu_ml/data/b5_2015.TIF', 'uu_ml/data/b6_2015.TIF']  # List to store file names\n",
    "\n",
    "# You can print to see how the file_list looks like\n",
    "print(file_list)\n",
    "\n",
    "# Read the files and stack them together by calling their names\n",
    "# Use the 'for' loop to iterate over the names to read files\n",
    "stack = np.array([])  # Empty array to store the stacked images\n",
    "for file in file_list:\n",
    "    img = imread(file)  # Read each image file\n",
    "    print(img.shape)  # Each time, also check the size of the image\n",
    "    \n",
    "    # In order to do clustering, image should be reshaped into a single column\n",
    "    img_col = img.reshape(-1, 1)\n",
    "    \n",
    "    # Each time put the reshaped image into the stack\n",
    "    stack = np.hstack((stack,img_col)) if stack.size else img_col\n",
    "    # Also to check the size of the stack\n",
    "    print(stack.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We also need to rasterize our manually delineated LULC types in the AOIs as we did in previous sessions\n",
    "# Use the rasterio again to rasterize the *.shp file\n",
    "\n",
    "from rasterio import features\n",
    "import pandas as pd\n",
    "\n",
    "# Labels from the AOIs\n",
    "aoi = gpd.read_file('uu_ml/data/aoi_NL_5_classes.shp')\n",
    "aoi['aoi_cat'] = pd.Categorical(aoi['class'])\n",
    "\n",
    "# Rasterize\n",
    "rst = rasterio.open('uu_ml/data/b5_2015.TIF')  # Base image to rasterize the *.shp\n",
    "meta = rst.meta.copy()  # Copy metadata from the base image\n",
    "meta.update(compress='lzw')\n",
    "\n",
    "# Burn the AOIs *.shp file into raster and save it\n",
    "out_rst = 'uu_ml/data/aoi_rasterized.tif'\n",
    "with rasterio.open(out_rst, 'w+', **meta) as out:\n",
    "    out_arr = out.read(1)\n",
    "\n",
    "    # Create a generator of geom, value pairs to use in rasterizing\n",
    "    shapes = ((geom,value) for geom, value in zip(aoi.geometry, aoi.aoi_cat))\n",
    "\n",
    "    burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)\n",
    "    out.write_band(1, burned)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will further stack the rasterized labels from the AOIs with the input bands, and only focus on the pixels within the AOIs.\n",
    "# The training and test datasets will be prepared from the pixels and labels within the AOIs\n",
    "\n",
    "# Load the rasterized LULC types in the AOI and concatenate it together with the images\n",
    "\n",
    "aoi_rst = rasterio.open('uu_ml/data/aoi_rasterized.tif').read(1)\n",
    "\n",
    "# Stack the label with the input bands\n",
    "data = np.c_[stack, aoi_rst.reshape(-1,)]\n",
    "\n",
    "# Of course, we are only interested in pixels with LULC type labelled\n",
    "data = data[np.where(data[:,data.shape[1]-1]!=0)]\n",
    "\n",
    "# Recall our function for preparing training and test datasets.\n",
    "# This time we re-write it a little bit to let the users of the function to split the data into training and test sets.\n",
    "\n",
    "def trainTestSplit(x, y, training_proportion):\n",
    "    data = np.c_[x, y]\n",
    "    np.random.shuffle(data)  # Shuffle the data so that LULC types can spread over training and test sets\n",
    "    x_train = data[:int(training_proportion*len(data)), :2]  # 70% of data for training\n",
    "    x_test = data[int(training_proportion*len(data)):, :2]  # 30% for testing\n",
    "    y_train = data[:int(training_proportion*len(data)), 2:].reshape(-1,)  # 70% of data for training\n",
    "    y_test = data[int(training_proportion*len(data)):, 2:].reshape(-1,)  # 30% for testing\n",
    "    return x_train, y_train, x_test, y_test\n",
    "\n",
    "# This time, use a very small proportion of the data for training, say, 30%.\n",
    "X_train, Y_train, X_test, Y_test = trainTestSplit(data[:,:-1], data[:,-1], 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the training and test datasets\n",
    "\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# Assign color codes to LULC types \n",
    "symbology2 = {11: cols[20],\n",
    "              21: cols[11],\n",
    "              31: cols[17],\n",
    "              52: cols[5],\n",
    "              61: cols[13]}\n",
    "\n",
    "cm = ListedColormap(symbology2.values())\n",
    "imin = min(symbology2)  # Colormap range\n",
    "imax = max(symbology2)\n",
    "\n",
    "# Visualize\n",
    "classes = ['Clear water', 'Sand', 'Residential', 'Deciduous forest', 'Agriculture']\n",
    "\n",
    "fig1,(ax1, ax2) = plt.subplots(1,2, figsize=(20,10))\n",
    "\n",
    "scatter1 = ax1.scatter(X_train[:, 0], X_train[:, 1], c=Y_train, cmap=cm, vmin=imin, vmax=imax, label='LULC types')\n",
    "ax1.set_title('Training data')\n",
    "ax1.legend(handles=scatter1.legend_elements()[0], labels=classes)\n",
    "\n",
    "ax2.scatter(X_test[:, 0], X_test[:, 1], c=Y_test, cmap=cm, vmin=imin, vmax=imax, label='LULC types')\n",
    "ax2.set_title('Test data')\n",
    "ax2.legend(handles=scatter1.legend_elements()[0], labels=classes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have both the **training** and **test** datasets, where the inputs are pixel values from the two bands (2-dimensional data points), and the outputs are categorical LULC types (encoded as numerical numbers).\n",
    "\n",
    "Given the fact that, as you already saw in previous sessions, the dataset we have is relatively simple to classify in the 2-dimensional feature space formed by the two bands, we will first try to inspect the **decision tree** with individual band, and see the model performance with only one band.\n",
    "\n",
    "Try to interpret the visualization in terms of values for splits, number of samples in each split, how the split compared to the patterns in the scatterplots in the 2-dimensional feature space above? And where are the problems?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import sklearn along with all necessary modules\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn import tree\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.metrics import r2_score,mean_squared_error\n",
    "\n",
    "# Initiate a tree model\n",
    "tree_depth = 3\n",
    "model_tree = DecisionTreeClassifier(max_depth=tree_depth)\n",
    "\n",
    "# Fit the model to your data.\n",
    "# Please note that the output of this fitting is a model with several parameters that are configurable, so far you only configured \"max_depth\" while training/fitting.\n",
    "model_tree.fit(X_train[:,0].reshape(-1,1), Y_train)\n",
    "\n",
    "# Visualize split\n",
    "fig = plt.figure(figsize=(tree_depth*4,tree_depth*3))\n",
    "tree.plot_tree(model_tree, filled=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='blue'>***Question 0.1: Please interpret the information at the nodes. For instance, what does the gini mean? What are the values in shown in the squared brackets? How the split sample numbers correspond to the data points shown in the feature space above?***</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further inspect the predictive performance by using the **confusion matrix** along with more detailed accuracy metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_pred = pd.Series(list(model_tree.predict(X_test[:,0].reshape(-1,1))), name='DT prediction')  # Store the predicted value in Y_pred\n",
    "Y_actu = pd.Series(list(Y_test), name='Manual delineation')\n",
    "\n",
    "# Map the LULC codes to the actual name of LULC types\n",
    "\n",
    "# First we need a mapping from the LULC codes to the actual LULC type name.\n",
    "code_lulc = { 52: 'Agriculture',\n",
    "              11: 'Clear water',\n",
    "              61: 'Deciduous forest',\n",
    "              31: 'Residential',\n",
    "              21: 'Sand'}\n",
    "\n",
    "# Now replace the non-intuitive numbers with actual LULC type names and store them into new variables\n",
    "Y_actu2 = Y_actu.replace(code_lulc)\n",
    "Y_pred2 = Y_pred.replace(code_lulc)\n",
    "\n",
    "# Show the LULC coded confusion matrix\n",
    "df_confusion2 = pd.crosstab(Y_actu2, Y_pred2)\n",
    "df_confusion2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='blue'>***Question 0.2: Now, do the accuracy metrics capture your visual comparison between the decision tree and the scatterplots in the feature space above. Which classes suffer the most from misclassification? And why?***</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print out more detailed accuracy assessment report\n",
    "\n",
    "from sklearn.metrics import classification_report\n",
    "\n",
    "print(classification_report(Y_actu2, Y_pred2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Further inspection: features and model configurations\n",
    "\n",
    "Now you may want to try both bands in the sample datasets, and experiment with model configurations. Please pay attention to how the splits are shown with two features (bands). And please do modify the model configuration such as the ***tree depth***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import tree\n",
    "\n",
    "# Initiate a tree model\n",
    "tree_depth = 3\n",
    "model_tree = DecisionTreeClassifier(max_depth=tree_depth)\n",
    "\n",
    "# Fit the model to your data.\n",
    "# Please note that the output of this fitting is a model with several parameters that are configurable, so far you only configured \"max_depth\" while training/fitting.\n",
    "model_tree.fit(X_train, Y_train)\n",
    "\n",
    "# Visualize split\n",
    "fig = plt.figure(figsize=(tree_depth*4,tree_depth*3))\n",
    "tree.plot_tree(model_tree, filled=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='blue'>***Question 0.3: You may notice that at each node, the data has been split by using one of the two features scripted as either [0] or [1]. Does the split value corresponds to the pattern in the scatterplot above?***</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Again, make predictions and compare with the known labels.\n",
    "\n",
    "Y_pred = pd.Series(list(model_tree.predict(X_test)), name='DT prediction')  # Store the predicted value in Y_pred\n",
    "Y_actu = pd.Series(list(Y_test), name='Manual delineation')\n",
    "\n",
    "# Now replace the non-intuitive numbers with actual LULC type names and store them into new variables\n",
    "Y_actu2 = Y_actu.replace(code_lulc)\n",
    "Y_pred2 = Y_pred.replace(code_lulc)\n",
    "\n",
    "# Show the LULC coded confusion matrix\n",
    "df_confusion2 = pd.crosstab(Y_actu2, Y_pred2)\n",
    "df_confusion2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print out more detailed accuracy assessment report\n",
    "\n",
    "from sklearn.metrics import classification_report\n",
    "\n",
    "report=classification_report(Y_actu2, Y_pred2,output_dict=True)\n",
    "\n",
    "rp = pd.DataFrame(report).transpose()\n",
    "rp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='green'>***Exercise 0.1: Please try attach precison and recall as extra row and column to the confusion matrix above to create a more coherent summary of the accuracy assessment.***</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do have a lot of **features** of the input *x* in the form of *(x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>,..., x<sub>n</sub>)*, you may soon get lost about how the data is split in the feature space and what the values mean in the tree nodes. However, there are many options to visualize the data splits along the tree. \n",
    "\n",
    "At the same time, you have already seen that some data value are better split along specific feature/dimension. It means that certain data point of a dimension or feature can be used best for capture data variance. Hence it is possible to rank feature importance in classification. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The tree can also be visualized in a different way to see how data has been splitted.\n",
    "# Let's plot how the features split the data\n",
    "\n",
    "# Plot parameters\n",
    "plot_colors = \"ryb\"\n",
    "plot_step = 100\n",
    "\n",
    "# Plot the decision boundary\n",
    "plt.figure(figsize=(6, 6))\n",
    "\n",
    "x_min, x_max = X_test[:, 0].min() - 1000, X_test[:, 0].max() + 1000\n",
    "y_min, y_max = X_test[:, 1].min() - 1000, X_test[:, 1].max() + 1000\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),\n",
    "                     np.arange(y_min, y_max, plot_step))\n",
    "plt.tight_layout(h_pad=0.5, w_pad=0.5, pad=2.5)\n",
    "\n",
    "Z = model_tree.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "Z = Z.reshape(xx.shape)\n",
    "cs = plt.contourf(xx, yy, Z, cmap=plt.cm.Blues)\n",
    "\n",
    "plt.xlabel('feature_1')\n",
    "plt.ylabel('feature_2')\n",
    "\n",
    "plt.scatter(X_test[:, 0], X_test[:, 1], c=Y_test.reshape(Y_test.shape[0]), cmap='Oranges', edgecolor='black', s=45)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='blue'>***Question 0.4: How will you interpret the usefulness of the features in classifying the data?***</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='green'>***Exercise 0.2: How is the classification? Which way of visualization do you prefer? Please try to configure your model with different depths of trees and inspect how that impact the data split by using one of the above visualization option.***</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get feature importance\n",
    "\n",
    "importance = model_tree.feature_importances_\n",
    "# summarize feature importance\n",
    "for i,v in enumerate(importance):\n",
    "    print('Feature: %0d, Score: %.5f' % (i,v))\n",
    "    \n",
    "# Plot feature importance\n",
    "plt.bar([ind for ind in range(len(importance))], importance)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='blue'>***Question 0.5: Does the ranking capture the usefulness of the feature shown above?***</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, once you are satisfied with the **training** and **testing**, use the model to predict on the larger area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict on the image stack\n",
    "# In this case, the more powerful non-linear kernel model is used\n",
    "\n",
    "Y_pred_all = model_tree.predict(stack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign color codes to LULC types \n",
    "symbology = {'Agriculture': cols[5],\n",
    "             'Clear water': cols[20],\n",
    "             'Deciduous forest': cols[3],\n",
    "             'Residential': cols[17],\n",
    "             'Sand': cols[11]}\n",
    "\n",
    "# Visualize\n",
    "fig1,(ax1, ax2) = plt.subplots(1,2, figsize=(20,10))\n",
    "show(b5_2020, ax=ax1, cmap='gray', alpha=0.25)\n",
    "aoi.plot(ax=ax1, column='land_cover', legend=True, color=aoi['land_cover'].map(symbology))\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_points = [Line2D([0], [0], marker=\"o\", linestyle=\"none\", markersize=5, color=color) for color in symbology.values()]\n",
    "leg_points = ax1.legend(custom_points, symbology.keys(), loc='upper right', frameon=False)\n",
    "ax1.add_artist(leg_points)\n",
    "\n",
    "# Assign color codes to LULC types \n",
    "symbology2 = {31: cols[17],\n",
    "              52: cols[5],\n",
    "              11: cols[20],\n",
    "              21: cols[11],\n",
    "              61: cols[3]}\n",
    "\n",
    "# Visualize\n",
    "# Because the predicted labels are still in one column, you need to reshape it back to original image shape\n",
    "row, col = img.shape  # Get the original dimensions of the image\n",
    "imin = min(symbology2)  # Colormap range\n",
    "imax = max(symbology2)\n",
    "\n",
    "print('Printing large image takes time...')\n",
    "ax2.imshow(Y_pred_all.reshape(row, col), cmap=cm, interpolation='none', vmin=imin, vmax=imax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can further save the predicted raster map into a georeferenced TIFF, so that we can inspect it along with other data or maps. Please feel free to drag the save TIFF file into other free and open-source software (FOSS) for further inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Again we use our image data for georeferencing information\n",
    "rst = rasterio.open('uu_ml/data/b5_2015.TIF')  # Base image to rasterize the *.shp\n",
    "meta = rst.meta.copy()  # Copy metadata from the base image\n",
    "meta.update(compress='lzw')\n",
    "\n",
    "# Burn the AOIs *.shp file into raster and save it\n",
    "out_rst = 'uu_ml/data/tree_prediction.tif'\n",
    "out_file = rasterio.open(\n",
    "    out_rst,\n",
    "    'w',\n",
    "    driver='GTiff',\n",
    "    height=row,\n",
    "    width=col,\n",
    "    count=1,\n",
    "    dtype=Y_pred_all.dtype,\n",
    "    crs=rst.crs,\n",
    "    transform=rst.transform)\n",
    "\n",
    "out_file.write(Y_pred_all.reshape(row, col),1)\n",
    "out_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='green'>***Exercise 0.3: As always, please compile your own program in a compact way, where you can deal with different input data for training, and also staying flexible with model configuration. Please use the .csv file containing multiple bands and corresponding labels to test your model performance on different input bands.***</font>\n",
    "\n",
    "*Start to load training data:*\n",
    "\n",
    "```\n",
    "import pandas as pd\n",
    "lulc = pd.read_csv('uu_ml/data/stack_aoi_2015.csv')\n",
    "\n",
    "# View some sample rows of the data\n",
    "lulc.head()\n",
    "```\n",
    "\n",
    "*You are also encourage to create your own AOI labels, load it, rasterize it and concantenate it together with the images:*\n",
    "\n",
    "```\n",
    "aoi = gpd.read_file('...')\n",
    "\n",
    "rst = rasterio.open('...')  # Base image to rasterize the *.shp\n",
    "meta = rst.meta.copy()  # Copy metadata from the base image\n",
    "meta.update(compress='lzw')\n",
    "\n",
    "...\n",
    "```\n",
    "\n",
    "*Train your model and inspect the training processes*\n",
    "\n",
    "```\n",
    "from sklearn... import ...\n",
    "    ...\n",
    "    ...\n",
    "tree_depth = ##    \n",
    "band_list = [##, ##, ..., ##]\n",
    "    \n",
    "```\n",
    "\n",
    "*Evaluate the accuracy with different metrics*\n",
    "\n",
    "```\n",
    "...\n",
    "...\n",
    "```\n",
    "\n",
    "*Generalize your prediction from the trained model to larger areas*\n",
    "\n",
    "```\n",
    "...\n",
    "...\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}