CUU_burn_extent_pixel_box_plots.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CUU_burn_extent_pixel_box_plots.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "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/burn-mapping/blob/master/CUU_burn_extent_pixel_box_plots.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"
      },
      "source": [
        "##Script to extract pixel values from a raster relating to a set of training polygons and then generate box plots\n",
        "Inputs are\n",
        "\n",
        "1)  indices generated from sentinel 2 imagery (generated from previous script) \n",
        "\n",
        "2)\tpolygon dataset of training samples . Needs Poly_ID_ field containing polygon ID (short integer) and 'Class_labe' field containing the class label.\n",
        "\n",
        "Main steps:\n",
        "•\textracts pixel values for polygons\n",
        "• generates box plots using seaborn package\n",
        "\n",
        "Duncan Blake NatureScot 11 September 2020"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bQ4BbT9b5SGT"
      },
      "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": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZqbFH-cc48L1"
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import os\n",
        "import json\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 seaborn as sns"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YUp4Uh2R5sZh"
      },
      "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": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AMsQ8v3Ro1ns"
      },
      "source": [
        "Set working drive, file paths and dataset name parameters\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MYjD-LPt48Ml"
      },
      "source": [
        "\n",
        "wd = '/content/drive/My Drive'\n",
        "image_dir = os.path.join(wd, 'Imagery')\n",
        "training_dir= os.path.join(wd, 'Training_Data')\n",
        "plots_dir = os.path.join(wd, 'Plots')\n",
        "# parameters for Skye\n",
        "# image = 'CUU_burn_skye_indices.tif'\n",
        "# training_polys = 'Box_plot_polygons_Skye.shp'\n",
        "# filename_diff = ('Skye_difference_indices_boxplot.png')\n",
        "# filename = ('Skye_postburn_indices_boxplot.png')\n",
        "\n",
        "# parameters for east cairngorms\n",
        "image = 'CUU_burn_cairngorms_indices.tif'\n",
        "training_polys = 'Box_plot_polygons_Cairngorms.shp'\n",
        "filename_diff = ('East_cairngorms_difference_indices_boxplot.png')\n",
        "filename = ('East_cairngorms_postburn_indices_boxplot.png')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVIFMdeuwhIu"
      },
      "source": [
        "###Open input image and polygon dataset"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yI995WqY48Nw"
      },
      "source": [
        "#Read image\n",
        "s2 = rasterio.open(os.path.join(image_dir,image))\n",
        "#Print number of bands\n",
        "B = s2.count\n",
        "print(B)\n",
        "# read polygons\n",
        "TaggedPolys= gpd.read_file(os.path.join(training_dir, training_polys))\n",
        "print(TaggedPolys.count)\n",
        "# show column headings\n",
        "for col in TaggedPolys.columns:\n",
        "  print(col)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1vnLpetG1oXs"
      },
      "source": [
        "##Extract all data for all pixels where the centroid falls within polygons, add polygon id and text class label to final 2 columns\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Dtuhh2un48Od"
      },
      "source": [
        "\n",
        "def getPixels(image, poly, indexInput, polygons,  target):\n",
        "    global B\n",
        "    shape=[mapping(poly)] # create list of coordinates in gjson format from gpandas?\n",
        "    #print((shape))\n",
        "\n",
        "    #reduce imagery to pixels overlapping polygon, noData values outside the polygon\n",
        "    outImage, out_transform = mask(image, shape, crop=True, nodata=np.nan)\n",
        "\n",
        "    # reshape output array to rows equal to number of bands, and number of columns to match input (-1) .T switches rows and columns\n",
        "    outList = outImage.reshape((B, -1)).T \n",
        "    # get current polygon ID\n",
        "    currentPolyID = polygons.loc[indexInput,\"Poly_ID_\"]\n",
        "    # creates 1D array of polyID, size equal to number of pixels (shape returns rows, columns)\n",
        "    currentPolyIDarr= np.repeat(currentPolyID, outList.shape[0])\n",
        "    # creates 2D array, 1column\n",
        "    currentPolyIDarr= currentPolyIDarr.reshape((outList.shape[0],1))\n",
        "    currentCategory =polygons.loc[indexInput,\"Class_labe\"]\n",
        "    currentCategoryarr= np.repeat(currentCategory, outList.shape[0])\n",
        "    # create 2D array of current class / categor\n",
        "    currentCategoryarr= currentCategoryarr.reshape((outList.shape[0],1))\n",
        "    # convert to dataframe to cope with string attribute\n",
        "    currentCategorydf = pd.DataFrame(currentCategoryarr)\n",
        "    currentCategorydf = currentCategorydf.astype(str)\n",
        "    #print(currentCategorydf.shape)\n",
        "\n",
        "    # add poly ID to pixel values\n",
        "    outList = np.concatenate((outList,currentPolyIDarr), axis = 1)\n",
        "    # convert to pandas dataframe\n",
        "    outList = pd.DataFrame(outList)\n",
        "    # add class to pixel values by appending two dataframes together\n",
        "    outList = np.append(outList,currentCategorydf, axis=1) \n",
        "    # drop noData values from the df\n",
        "    outList = pd.DataFrame(outList).dropna()\n",
        "    #print(outList.shape)\n",
        "    #print(outList.head())\n",
        "    #print((\"Number of pixels = \" + str(outList.shape[0])))\n",
        "    return np.append(target, outList, axis=0)\n",
        "\n",
        "\n",
        "def extractAllPolygons(image, featuresgeom, features):\n",
        "    global B \n",
        "    # number of columns in extracted pixel dataset\n",
        "    finalcolno = B+2 \n",
        "    # create empty dataset with number of colums set and datatype set to float\n",
        "    flatten = np.array([]).reshape(0,finalcolno).astype(float)\n",
        "    # iterate through each polygon\n",
        "    for index, f in enumerate(featuresgeom): \n",
        "      indexInput= index\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()\n",
        "\n",
        "# run both of above functions, input is imagery, geometry part of pandas arraym, gp dataframe\n",
        "totValues = extractAllPolygons(s2,TaggedPolys.geometry.values, TaggedPolys)\n",
        "print(totValues.head())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fXigZlPJcqs1"
      },
      "source": [
        "## Generate box plots"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2yC9Iy-jcwRF"
      },
      "source": [
        "# set the 23 column names\n",
        "totValues.columns = ['prenbr','postnbr','dnbr','prenbr2','postnbr2','dnbr2','prendvi','postndvi','dndvi','prenmirbi','postnmirbi','dnmirbi',\n",
        "                     'precsi','postcsi','dcsi','prebais2','postbais2','dbais2','presave','postsavi','dsavi','Poly_ID_','Class']\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xNKOZ11WE_3i"
      },
      "source": [
        "### Run loop to create box plots for difference indices which should be run with the change classes only"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6bJZqBk-FU3H"
      },
      "source": [
        "# Subset to contain difference indices only\n",
        "df_differences = totValues[['dnbr','dnbr2','dndvi','dnmirbi','dcsi','dbais2','dsavi','Class']]\n",
        "print(df_differences.head())\n",
        "\n",
        "# next select only the change classes we want to compare\n",
        "dClassList = [\"Burns\", \"Field to bare field\", \"Trees to no trees\", \"Cloud to no cloud\",\"No cloud to cloud\",\"Shadow to no shadow\",\n",
        "              \"No shadow to shadow\", \"Snow to no snow\"]\n",
        "\n",
        "# define class list as the order we want to display the classes in the boxplots.  This will include any categories with no data though.\n",
        "orderList = dClassList\n",
        "\n",
        "df_differences_classes = df_differences[df_differences.Class.isin(dClassList)]\n",
        "print(df_differences_classes.head())\n",
        "\n",
        "# and select the list of difference indices we want to create box plots for\n",
        "dList = ['dnbr','dnbr2','dndvi','dnmirbi','dcsi','dbais2','dsavi']\n",
        "\n",
        "fig, axs = plt.subplots(4,2, figsize=(12,15))\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[0] , data = df_differences_classes, ax=axs[0][0], order=orderList, showfliers=False).set_title(dList[0])\n",
        "_ = axs[0][0].get_xaxis().set_visible(False) \n",
        "_ = plt.ylabel(dList[0] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[1] , data = df_differences_classes, ax=axs[0][1], order=orderList, showfliers=False).set_title(dList[1])\n",
        "_ = axs[0][1].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[1] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[2] , data = df_differences_classes, ax=axs[1][0], order=orderList, showfliers=False).set_title(dList[2])\n",
        "_ = axs[1][0].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[2] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[3] , data = df_differences_classes, ax=axs[1][1], order=orderList, showfliers=False).set_title(dList[3])\n",
        "_ = axs[1][1].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[3] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[4] , data = df_differences_classes, ax=axs[2][0], order=orderList,  showfliers=False).set_title(dList[4])\n",
        "_ = axs[2][0].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[4] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[5] , data = df_differences_classes, ax=axs[2][1], order=orderList, showfliers=False).set_title(dList[5])\n",
        "_ = axs[2][1].set_xticklabels(axs[2][1].get_xticklabels(), rotation=45, horizontalalignment='right')\n",
        "_ = plt.ylabel(dList[5] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[6] , data = df_differences_classes, ax=axs[3][0], order=orderList, showfliers=False).set_title(dList[6])\n",
        "_ = axs[3][0].set_xticklabels(axs[3][0].get_xticklabels(), rotation=45, horizontalalignment='right')\n",
        "_ = plt.ylabel(dList[6] + ' values')\n",
        "\n",
        "_ = axs[3][1].set_visible(False)\n",
        "\n",
        "# export the figure - needs to happen before displaying - reduce whitespace with the bbox_inches command\n",
        "plt.savefig(f\"{plots_dir}/{filename_diff}\", bbox_inches='tight')\n",
        "# display the figure\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2OlG4vwRKH9-"
      },
      "source": [
        "### Run loop to create box plots for post burn indices which should be run with the land cover classes only"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j2l18FnyKH-E"
      },
      "source": [
        "# Subset to contain post burn indices only\n",
        "df_differences = totValues[['postnbr','postnbr2','postndvi','postnmirbi','postcsi','postbais2','postsavi','Class']]\n",
        "print(df_differences.head())\n",
        "\n",
        "# next select only the change classes we want to compare\n",
        "dClassList = [\"Burns\", \"Bare field\", \"Bare peat\", \"Felled\", \"Rock\", \n",
        "              \"Shadow\", \"Vegetation\", \"Water\"]\n",
        "\n",
        "# define class list as the order we want to display the classes in the boxplots.  This will include any categories with no data though.\n",
        "orderList = dClassList\n",
        "df_differences_classes = df_differences[df_differences.Class.isin(dClassList)]\n",
        "print(df_differences_classes.head())\n",
        "\n",
        "# and select the list of indices we want to create box plots for\n",
        "dList = ['postnbr','postnbr2','postndvi','postnmirbi','postcsi','postbais2','postsavi']\n",
        "\n",
        "fig, axs = plt.subplots(4,2, figsize=(12,15))\n",
        "\n",
        "# x axis ordered and outliers are removed because some were quite extreme and the plots were hard to interpret\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[0] , data = df_differences_classes, ax=axs[0][0], order=orderList, showfliers=False).set_title(dList[0])\n",
        "_ = axs[0][0].get_xaxis().set_visible(False) \n",
        "_ = plt.ylabel(dList[0] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[1] , data = df_differences_classes, ax=axs[0][1], order=orderList, showfliers=False).set_title(dList[1])\n",
        "_ = axs[0][1].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[1] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[2] , data = df_differences_classes, ax=axs[1][0], order=orderList, showfliers=False).set_title(dList[2])\n",
        "_ = axs[1][0].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[2] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[3] , data = df_differences_classes, ax=axs[1][1], order=orderList, showfliers=False).set_title(dList[3])\n",
        "_ = axs[1][1].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[3] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[4] , data = df_differences_classes, ax=axs[2][0], order=orderList, showfliers=False).set_title(dList[4])\n",
        "_ = axs[2][0].get_xaxis().set_visible(False)\n",
        "_ = plt.ylabel(dList[4] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[5] , data = df_differences_classes, ax=axs[2][1], order=orderList, showfliers=False).set_title(dList[5])\n",
        "_ = axs[2][1].set_xticklabels(axs[2][1].get_xticklabels(), rotation=45, horizontalalignment='right')\n",
        "_ = plt.ylabel(dList[5] + ' values')\n",
        "\n",
        "_ = sns.set_style(\"whitegrid\") \n",
        "_ = sns.boxplot(x = 'Class', y = dList[6] , data = df_differences_classes, ax=axs[3][0], order=orderList, showfliers=False).set_title(dList[6])\n",
        "_ = axs[3][0].set_xticklabels(axs[3][0].get_xticklabels(), rotation=45, horizontalalignment='right')\n",
        "_ = plt.ylabel(dList[6] + ' values')\n",
        "\n",
        "_ = axs[3][1].set_visible(False)\n",
        "\n",
        "# export the figure - needs to happen before displaying - reduce whitespace with the bbox_inches command\n",
        "plt.savefig(f\"{plots_dir}/{filename}\", bbox_inches='tight')\n",
        "# display the figure\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}