Data_Split_Classification_.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Data_Split_Classification_for_github.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "machine_shape": "hm",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/duncansnh/Bare-peat/blob/master/Data_Split_Classification_for_github.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AzhcWTCF4xvQ",
        "colab_type": "text"
      },
      "source": [
        "##Script to perform supervised classification. Applied to bare peat classification on nationwide scale using sentinel 2.  Joblog 92279 1November 2019  \n",
        "Inputs are\n",
        "•  indices generated from sentinel 2 imagery (23 layers are generated from previous script,  but other numbers of layers may be applied) \n",
        "•\tpolygon dataset of training samples (known classes  of bare peat = 10, other = 20, rock/ stone=30, shadow =40, water=50.). ID field must contain the class label (short integer), Poly_ID field contains unique polygon ID (short integer).\n",
        "Main steps:\n",
        "•\textracts pixel values for polygons\n",
        "• resamples pixels , default is 300 pixel per class (total number pixels required from each polygon are estimated based on polygon size, sampling is random where total number of pixels required is less than that in polygon, if number of pixels required is greater than number in polygon then all are sampled and addtional pixels are randomly sampled.)\n",
        "•\tsplits into training/ validation data (does not take into account source polygon), scales values based on training data\n",
        "•\tshows accuracy matrices for 6 classifiers : random forest, KNN, logistic regression, xg boost, linear discriminant analysis, ensemble classifier\n",
        "•\tneeds user input to choose most accurate classifier\n",
        "•\tclassifies input image and writes out:\n",
        "  o\tProbability for class == bare peat.\n",
        "\n",
        "  Currently based on running in Google co- lab, with folders in google drive labelled:\n",
        "  Imagery, Training_Data\n",
        "\n",
        "Adapted from original script from Tom Wilson, Foresty Research, which was used with sentinel 1 data to classify forestry as felled/ mature trees. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G9KdEVlBTyNl",
        "colab_type": "text"
      },
      "source": [
        "##To use a different number of bands:\n",
        " \n",
        "\n",
        "*   Update the script to the new number in the the 9th code block ('extract all data for all pixels'), where code is highlighted with lines of ######, and select which bands should be included (again highlighted in this code block with lines of #########).\n",
        "*   Update the final code block with band to be excluded (highlighted with lines of ####)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bQ4BbT9b5SGT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#This is only required if running in colab notebook to install the libraries\n",
        "#If running Python code elsewhere need to make sure below libraries are installed\n",
        "! pip install geopandas\n",
        "! pip install descartes\n",
        "! pip install rasterio\n",
        "! pip install rasterstats"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZqbFH-cc48L1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import os\n",
        "import json\n",
        "from sklearn.model_selection import train_test_split\n",
        "import geopandas as gpd\n",
        "import descartes\n",
        "import rasterio\n",
        "from rasterio.mask import mask\n",
        "from rasterio.features import geometry_mask\n",
        "from shapely.geometry import mapping\n",
        "from rasterstats import zonal_stats\n",
        "import datetime\n",
        "import math"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YUp4Uh2R5sZh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#Only if running in Google Colab, in which case input image, training polygons and output results need to be in Google Drive.\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive/', force_remount=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AMsQ8v3Ro1ns",
        "colab_type": "text"
      },
      "source": [
        "set working drive, iteration (for naming outputs), number of pixels per class"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MYjD-LPt48Ml",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "wd = '/content/drive/My Drive'\n",
        "iteration = 'ML13_v3'\n",
        "samplesize = 300 #number of pixels to be sampled per class (approximate)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3BaBVh-E48M-",
        "colab_type": "text"
      },
      "source": [
        "###Set image directory, training directory"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x-Kn8iYc48NE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "image_dir = os.path.join(wd, 'Imagery')\n",
        "training_dir= os.path.join(wd, 'Training_Data')\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVIFMdeuwhIu",
        "colab_type": "text"
      },
      "source": [
        "###'open' input image"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yI995WqY48Nw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#Read image\n",
        "s2 = rasterio.open(os.path.join(image_dir,'ML13_23_indices.tif'))\n",
        "#Print number of bands\n",
        "B = s2.count\n",
        "print(B)\n",
        "print(s2.shape)\n",
        "#Copy raster profile to for later output\n",
        "s2prof = s2.profile.copy()\n",
        "s2prof.update(count = 1, nodata=None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BaQJNmawfgIW",
        "colab_type": "text"
      },
      "source": [
        "OPTIONAL - get max and min values from each band ( if raster is small) "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TeTrI7gXffgY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "array = s2.read()\n",
        "\n",
        "# Calculate statistics for each band\n",
        "stats = []\n",
        "\n",
        "for band in array:\n",
        "  stats.append({\n",
        "    'min': band.min(),\n",
        "#   'mean': band.mean(),\n",
        "#   'median': np.median(band),\n",
        "    'max': band.max()})\n",
        "\n",
        "del(array)\n",
        "print((stats))\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WcaqvPsW48N_",
        "colab_type": "text"
      },
      "source": [
        "#open polygon dataset"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DjMY7BlelSPj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "TaggedPolys= gpd.read_file(os.path.join(training_dir, 'Training_ML13.shp'))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1vnLpetG1oXs",
        "colab_type": "text"
      },
      "source": [
        "##extract all data for all pixels where the centroid falls within polygons, add polygon id and category id to final 2 columns\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Dtuhh2un48Od",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "#####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "#ALTER B if changing the number of bands\n",
        "#B = 22\n",
        "B=23\n",
        "\n",
        "#####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "#####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "\n",
        "\n",
        "def getPixels(image, poly, indexInput, polygons,  target):\n",
        "    global B\n",
        "    shape=[mapping(poly)] \n",
        "    print((shape))\n",
        "    print((\"got to -1\"))\n",
        "\n",
        "\n",
        "    ######################## USE THE 2 LINES BELOW TO UPDATE BAND LIST IF NEEDING TO EXCLUDE BAND, COMMENT OUT FIRST LINE##############################\n",
        "    #####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "    outImage, out_transform = mask(image, shape, crop=True, nodata=np.nan)#reduce imagery to pixels overlapping polygon\n",
        "    #bandlist = [1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20, 21,22,23] # EXCLUDES DARKNESS\n",
        "    #outImage, out_transform = mask(image, shape, crop=True, nodata=np.nan, indexes = bandlist)# optionally reduce number of bands with indexes, CHANGE B if this is the case\n",
        "    \n",
        "    \n",
        "    #####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "    #####################################XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX###########################################\n",
        "\n",
        "\n",
        "    # if nodata is set to a figure pixels where centroid is outwith polygon are included and not excluded with drop.na\n",
        "    outList = outImage.reshape((B, -1)).T# reshape output array to rows equal to number of bands, and number of columns to match input (-1) \n",
        "    \n",
        "    currentPolyID = polygons.loc[indexInput,\"Poly_ID\"]# get current polygon ID\n",
        "   \n",
        "    currentPolyIDarr= np.repeat(currentPolyID, outList.shape[0])# creates 1D array of polyID, size equal to number of pixels (shape returns rows, columns)\n",
        "    currentPolyIDarr= currentPolyIDarr.reshape((outList.shape[0],1))# creates 2D array, 1column\n",
        "    currentCategory =polygons.loc[indexInput,\"ID\"]\n",
        "    currentCategoryarr= np.repeat(currentCategory, outList.shape[0])\n",
        "    currentCategoryarr= currentCategoryarr.reshape((outList.shape[0],1))# create 2D array of current class / category\n",
        "    \n",
        "    outList = np.concatenate((outList,currentPolyIDarr), axis = 1)# add poly ID to pixel values\n",
        "    outList = np.append(outList,currentCategoryarr, axis=1)# add class to pixel values\n",
        "    outList = pd.DataFrame(outList).dropna()\n",
        " \n",
        "    return np.append(target, outList, axis=0)\n",
        "\n",
        "\n",
        "def extractAllPolygons(image, featuresgeom, features):\n",
        "    global B # number of bands in input imagt\n",
        "    finalcolno = B+2 # number of columns in extracted pixel dataset\n",
        "    flatten = np.array([]).reshape(0,finalcolno).astype(float)# empty dataset with number of colums set and datatype set to float\n",
        "    for index, f in enumerate(featuresgeom): #iterate through each polygon\n",
        "      indexInput= index# iteration number\n",
        "      flatten = getPixels(image,f,indexInput, features, flatten)\n",
        "    flattenArr = np.ma.masked_array(flatten, mask=(flatten == np.nan))\n",
        "    return pd.DataFrame(flattenArr).dropna()# remove any na - machine learning models can't deal with them \n",
        "\n",
        "\n",
        "totValues = extractAllPolygons(s2,TaggedPolys.geometry.values, TaggedPolys)# run both of above functions, input is imagery, geometry part of pandas arraym, gp dataframe\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T93RO0nYdN0g",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print((totValues.shape))\n",
        "print((totValues.size))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GQzHCJla48Qe",
        "colab_type": "text"
      },
      "source": [
        "#### select appropriate number of pixels for each polygon, based on size of polygon and total size of area for each class. Select all data and resample where required number of pixels exceeds total pixels in polygon, where number of pixels required is less than that in polygon take random selection.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_Gccj-C_48Qh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#get unique class IDs\n",
        "ClassColumnIndex = B + 1 # get column index for class ID\n",
        "PolygonColumnIndex = B # get column index for class ID\n",
        "Classes = totValues.iloc[:,ClassColumnIndex].unique()\n",
        "\n",
        "\n",
        "#iloc is index rather than name\n",
        "FinValues = pd.DataFrame()#create empty pandas dataframe\n",
        "for Class in Classes: #iterate through each class\n",
        "  ClassValues = totValues[totValues.iloc[:,ClassColumnIndex]==Class]# select all rows from totValues for this class \n",
        "  totpixels = ClassValues.shape[0] # get number of rows (equals number of pixels) for this class from all training polygons\n",
        "  print((' class :{}  tot pixels = {}'. format(Class,totpixels) ))\n",
        "\n",
        "  ClassPolygonIDs = ClassValues.iloc[:,B].unique() #obtain polyIDs for this class\n",
        "  for polyID in ClassPolygonIDs: #iterate through polygons, sample from each polygon based on size\n",
        "    ClassPolygonValues = ClassValues[ClassValues.iloc[:,B]==polyID]# pixel values for this class , this polygon\n",
        "    PolySize = ClassPolygonValues.shape[0] # get number of rows (equals number of pixels) for this polygon\n",
        "    ReqPixels = ((int(math.ceil((PolySize/totpixels)*samplesize)))) #number of pixels required from this polygon, taking into account overall number in training data and total pixels required for class\n",
        "    print((\"Total pixels for poly = {}\".format(PolySize)))\n",
        "    print((\"Number of pixels for class {}, polyId {} is :{}\".format(Class, polyID,ReqPixels )))\n",
        "\n",
        "    if PolySize>= ReqPixels:# random selection if number of pixels in polygon is greater than number required\n",
        "      SelectedValues = ClassPolygonValues.sample(ReqPixels, replace=False)\n",
        "      \n",
        "      FinValues = FinValues.append(SelectedValues)\n",
        "      \n",
        "    elif PolySize < ReqPixels:# select all pixels and then resample with replacement if number of pixels in polygon is lower than number required\n",
        "      Extrasamplesize = ReqPixels - PolySize\n",
        "      ExtraValues = ClassPolygonValues.sample(Extrasamplesize, replace=True)\n",
        "      print(type(ExtraValues))\n",
        "      SelectedValues = pd.concat((ExtraValues,ClassPolygonValues), axis = 0)\n",
        "      FinValues = FinValues.append(SelectedValues)\n",
        "      print((FinValues.shape))\n",
        "\n",
        "print((type(Classes)))\n",
        "print((Classes))\n",
        "\n",
        "\n",
        "\n",
        "print(FinValues.shape)\n",
        "print((FinValues.iloc[0,:]))\n",
        "\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-8hKDbX0d_g1",
        "colab_type": "text"
      },
      "source": [
        "##Split into training and validation datasets, (not grouped by polygon so pixels from the same polygon may appear in training AND test samples - ensures fraction of split is honoured and may give better output classification as uses data from wider range of pixels )"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z468uCind4oT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#iterate through each class\n",
        "\n",
        "#unique class IDs from previous step are in Classes\n",
        "\n",
        "trainX =pd.DataFrame(columns=range(B))\n",
        "testX = pd.DataFrame(columns=range(B))\n",
        "trainy = pd.DataFrame(columns=range(1))\n",
        "testy = pd.DataFrame(columns=range(1))\n",
        "print((testy.shape))\n",
        "\n",
        "\n",
        "\n",
        "for Class in Classes:\n",
        "  print((Class))\n",
        "  ClassSampledValues = FinValues[FinValues.iloc[:,ClassColumnIndex]==Class ] #select current class from all pixel values\n",
        "\n",
        "\n",
        "  ClassSampledValuesTrain, ClassSampledValuesTest = train_test_split(ClassSampledValues, test_size = 0.2, random_state = 999 )\n",
        "\n",
        "  trainXClass = ClassSampledValuesTrain.iloc[:,0:B]#select training data colums\n",
        "  #selects columns from index 0 up to but not including the number of bands - so if you start counting from 1 it selects up to  column B\n",
        "  trainyClass = ClassSampledValuesTrain.iloc[:, ClassColumnIndex].values.reshape(-1,1)# select class column only \n",
        "\n",
        "  testXClass = ClassSampledValuesTest.iloc[:,0:B]#select training data colums\n",
        "  testyClass = ClassSampledValuesTest.iloc[:, ClassColumnIndex].values.reshape(-1,1)# select class column only \n",
        " \n",
        "  trainX =  np.append(trainX, trainXClass, axis = 0)\n",
        "  testX = np.append(testX, testXClass, axis = 0)\n",
        "  trainy = np.append(trainy,trainyClass, axis = 0)\n",
        "  testy = np.append(testy, testyClass, axis = 0)\n",
        "  \n",
        "\n",
        "  \n",
        "\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Sbw2ak3WQsSS",
        "colab_type": "text"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RD4T6LKC48O3",
        "colab_type": "text"
      },
      "source": [
        "### scale DFs based on training data\n",
        "May need to Check no missing, 0 or repeated fields in dataframe or something has gone wrong in image creation process "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ygy58cG448O6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.preprocessing import StandardScaler\n",
        "sc = StandardScaler()\n",
        "\n",
        "trainX = sc.fit_transform(trainX)# scale training data\n",
        "testX = sc.transform(testX)# scale test data based on training data \n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bXbSLdBX48PI",
        "colab_type": "text"
      },
      "source": [
        "### Supervised ML classifiers. Scale training, testing matrices, train different models, show different accuracy metrics\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "liHUgL8m48PV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.metrics import confusion_matrix, accuracy_score, f1_score\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LfOa2YvN48Pe",
        "colab_type": "text"
      },
      "source": [
        "Random Forest"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mfaRef9p48Ph",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "trainy= trainy.astype(int)#  class format needs to be particular format - int works - for classifer\n",
        "testy= testy.astype(int)#  class format needs to be particular format - int works - for classifer\n",
        "\n",
        "from sklearn.ensemble import RandomForestClassifier\n",
        "\n",
        "rf = RandomForestClassifier( random_state = 99, max_features=2, n_estimators=1000)\n",
        "                            \n",
        "# note that higher values of  n_estimators (n trees in the forest) may increase time to generate final output \n",
        "#max_features = max number of features considered for splitting a node\n",
        "rf.fit(trainX, trainy)# train model\n",
        "pred_y = rf.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y, average = None  )\n",
        "\n",
        "print (\"RF:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*99, 3),(f1)))\n",
        "\n",
        "print((\"Parameters currently in use : {}\".format((rf.get_params()))))\n",
        "\n",
        "print((\"classes: {}\".format(rf.classes_)))\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FutlH7MxbUBR",
        "colab_type": "text"
      },
      "source": [
        "#print feature importance table and plot"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jXQjtgBibbJU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "from sklearn.ensemble import ExtraTreesClassifier\n",
        "\n",
        "#get importances\n",
        "importances = rf.feature_importances_\n",
        "\n",
        "\n",
        "std = np.std([tree.feature_importances_ for tree in rf.estimators_],\n",
        "             axis=0)\n",
        "indices = np.argsort(importances)[::-1] # invert sorted array\n",
        "\n",
        "# Print the feature rank\n",
        "print(\"Feature ranking:\")\n",
        "\n",
        "for f in range(trainX.shape[1]):\n",
        "    print(\"{}. feature {} ({})\" .format(f + 1, indices[f]+1, importances[indices[f]]))\n",
        "\n",
        "# Plot the feature importances of the forest\n",
        "plt.figure()\n",
        "plt.title(\"Feature importances\")\n",
        "plt.bar(range(trainX.shape[1]), importances[indices],\n",
        "       color=\"r\", yerr=std[indices], align=\"center\")\n",
        "plt.xticks(range(trainX.shape[1]), indices+1)\n",
        "plt.xlim([-1, trainX.shape[1]])\n",
        "plt.show()\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qbtn14QlnzXZ",
        "colab_type": "text"
      },
      "source": [
        "#Optional  - fine tune RF paramters - note that it is worth trying default values too , as not all combinations are tested in the fine tune . Note that fine tune for random selection of training data may not be the same when applied to all data combined"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I2oaGiUrn3Xy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "\n",
        "import sklearn.ensemble  \n",
        "import sklearn.metrics  \n",
        "from sklearn.model_selection import GridSearchCV\n",
        " \n",
        "\n",
        "\n",
        "rftune = sklearn.ensemble.RandomForestClassifier()  \n",
        "\n",
        "param_grid = {  \n",
        "           \"n_estimators\" : [5, 10, 100,500, 1000],  \n",
        "           \"max_features\" : [\"auto\", \"sqrt\", \"log2\", 2,3,4,10], #add more parameters here if required - but this will increase time to run\n",
        "           \"random_state\" : [99]\n",
        "           }  \n",
        "\n",
        "flattrainy = np.ndarray.flatten(trainy)\n",
        "CV_rf = GridSearchCV(estimator=rftune, param_grid=param_grid)  \n",
        "CV_rf.fit(trainX, flattrainy)  \n",
        "print(CV_rf.best_params_)  "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R04JE3p_nzNF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O6Y9cFYt48Po",
        "colab_type": "text"
      },
      "source": [
        "KNN"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R2t0oUxL48Pp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.neighbors import KNeighborsClassifier\n",
        "knn = KNeighborsClassifier(n_neighbors = 2, metric = 'minkowski', p = 2)\n",
        "knn.fit(trainX, trainy)\n",
        "pred_y = knn.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y, average=None)\n",
        "print (\"KNN:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*100, 3),(f1)))\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k9NMYocm48Py",
        "colab_type": "text"
      },
      "source": [
        "Logistic Regression"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tOCdhCb248P0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.linear_model import LogisticRegression\n",
        "lgr = LogisticRegression(random_state = 99, max_iter=5000)\n",
        "lgr.fit(trainX, trainy)\n",
        "pred_y = lgr.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y, average=None)\n",
        "print (\"logistic regression:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*100, 3),(f1)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7d-VLwUh48P_",
        "colab_type": "text"
      },
      "source": [
        "XG Boost"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "przVHjKt48QD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from xgboost import XGBClassifier\n",
        "xg = XGBClassifier()\n",
        "xg.fit(trainX, trainy)\n",
        "pred_y = xg.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y,average=None)\n",
        "print (\"XGBoost:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*100, 3),(f1)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SWuRp5fm48QK",
        "colab_type": "text"
      },
      "source": [
        "Linear Discriminant Analysis"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EzO7IZD048QM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\n",
        "ld = LDA()\n",
        "ld.fit(trainX, trainy)\n",
        "pred_y = ld.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y, average=None)\n",
        "print (\"LDA:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*100, 3),(f1)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4AicvOK848QV",
        "colab_type": "text"
      },
      "source": [
        "Ensemble classifier - try different combinations of models to try and get zero error on class of interest (bare peat is class 10, first row/ column) - it may be possible to 'balance' errors out eg if one model has class 10 ommission error and another has class 10 commission error for the problematic class - usually 20)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YNr073IS48QZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.ensemble import VotingClassifier\n",
        "#voting_clf = VotingClassifier(estimators=[('lr', lgr),('rf', rf),('knn',knn),('xg',xg)],voting='soft')#voting='hard'\n",
        "voting_clf = VotingClassifier(estimators=[('rf',rf) ,('knn',knn)],voting='soft')#voting='hard'\n",
        "voting_clf.fit(trainX, trainy)\n",
        "pred_y = voting_clf.predict(testX)\n",
        "cm = confusion_matrix(testy, pred_y)\n",
        "acc = accuracy_score(testy, pred_y)\n",
        "f1 = f1_score(testy, pred_y, average=None)\n",
        "print (\"Ensemble:\\n{0}\\nOverall: {1}%\\nF1: {2}\".format(cm,round(acc*100, 3),(f1)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FIJXI_Zmd_WT",
        "colab_type": "text"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xljON-5751r1",
        "colab_type": "text"
      },
      "source": [
        "Optional - run on all data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6JDuZV2U54e8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "alltrainX = np.append(testX, trainX, axis = 0)\n",
        "alltrainy = np.append(testy, trainy, axis = 0)\n",
        "alltrainy = np.ndarray.flatten(alltrainy)\n",
        "#knn.fit(alltrainX, alltrainy)\n",
        "#rf.fit(alltrainX, alltrainy)\n",
        "\n",
        "\n",
        "voting_clf.fit(alltrainX, alltrainy)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lq1BkByL48Qo",
        "colab_type": "text"
      },
      "source": [
        "Output probability image - this stage takes longer for ensemble classifier and time may also depend on parameters of random forest "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aFSMUg7448Qq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "\n",
        "#Select model of choice\n",
        "model = voting_clf #xg, rf, knn, lgr ,voting_clf\n",
        "\n",
        "\n",
        "\n",
        "##NOTE - unable to predict probability for ensemble classifer when voting = hard.\n",
        "##NOTE - knn probability output only has a few categories\n",
        "s2prof.update(count=1, nodata=None, dtype=np.float32)\n",
        "dst = rasterio.open(os.path.join(image_dir,'{}_RF_KNN_alldata1000trees.tif'.format(iteration)), 'w', **s2prof)\n",
        "\n",
        "for block_index, window in s2.block_windows(1):\n",
        "    s2_block = s2.read(window=window, masked=True)\n",
        "    ######################################################################################\n",
        "    #SELECTION OF BANDS##########################################\n",
        "    # delete the band that needs excluding (band number minus 1, eg darkness is index 13)\n",
        "    #s2_block = np.delete(s2_block,13, 0)\n",
        "    ################################################################################\n",
        "    v= s2_block.shape\n",
        "    s2_block = s2_block.reshape(B, -1).T\n",
        "    s2_block = sc.transform(s2_block)\n",
        "    s2_block[s2_block<-3.4e+35]=9999\n",
        "    result_block = model.predict_proba(s2_block).astype('float32')\n",
        "    #select probabilies for class 1 only\n",
        "    result_block = result_block[:,0]\n",
        "    result_block = result_block.reshape(1,v[1],v[2])\n",
        "    dst.write(result_block, window=window)\n",
        "dst.close()"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}