examples/user_guide/5_Classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xarray-spatial\n",
    "### User Guide: Classification tools\n",
    "-----\n",
    "\n",
    "The classification tools let you reclassify the values in an xarray DataArray into a new set of values based on set bins.\n",
    "\n",
    "[Generate terrain](#Generate-Terrain-Data)\n",
    "[Reclassify](#Reclassify)\n",
    "\n",
    "-----------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's use datashader to render our images...\n",
    "\n",
    "First, we'll import numpy and some datashader functions to help with rendering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "import datashader as ds\n",
    "from datashader.transfer_functions import shade\n",
    "from datashader.transfer_functions import stack\n",
    "from datashader.transfer_functions import dynspread\n",
    "from datashader.transfer_functions import set_background\n",
    "from datashader.colors import Elevation\n",
    "\n",
    "import xrspatial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Terrain Data\n",
    "\n",
    "To test out the classification functions, we'll need some rasterized data. We can generate an artificial digital elevation model (*dem*), or terrain, with xarray-spatial's `generate_terrain`, with the help of datashader's Canvas for the aggregation of values into a raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import generate_terrain\n",
    "from xrspatial import hillshade\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "template_terrain = xr.DataArray(np.zeros((H, W)))\n",
    "x_range=(-20e6, 20e6)\n",
    "y_range=(-20e6, 20e6)\n",
    "\n",
    "terrain = generate_terrain(\n",
    "    template_terrain, x_range=x_range, y_range=y_range, seed=1, zfactor=1000\n",
    ")\n",
    "stack(shade(hillshade(terrain), cmap=['grey', 'white']), shade(terrain, cmap=Elevation, alpha=128))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reclassify"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantile Reclassify\n",
    "\n",
    "One method of reclassification is by *quantile*. In this method, the set of all values is divided into bins such that each bin contains the same number of cell values. This method works better for non-evenly distributed value sets.\n",
    "\n",
    "Let's try it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import hillshade\n",
    "from xrspatial import quantile\n",
    "\n",
    "qcut_agg = quantile(terrain, k=15)\n",
    "\n",
    "stack(shade(hillshade(qcut_agg), cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(qcut_agg,     cmap=Elevation,         alpha=128, how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equal Interval Reclassify\n",
    "\n",
    "Another method of reclassification is equal interval. This simply divides the entire values range by the given number of bins and assigns the vlaues into those bins based on where the value lies.\n",
    "\n",
    "As you can see below, for our terrain, this flattens out a lot of our values since our data is not evenly distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial.classify import equal_interval\n",
    "from xrspatial import hillshade\n",
    "\n",
    "equal_interval_agg = equal_interval(terrain, k=15)\n",
    "\n",
    "stack(shade(hillshade(equal_interval_agg), cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(equal_interval_agg, cmap=Elevation, alpha=128, how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Natural Breaks (Jenks) Reclassify\n",
    "\n",
    "This is another non-linear classification method that is best for non-evenly distributed data that does not skew only towards the high or low range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial.classify import natural_breaks\n",
    "from xrspatial import hillshade\n",
    "\n",
    "natural_breaks_agg = natural_breaks(terrain, num_sample=1000, k=15)\n",
    "\n",
    "stack(shade(hillshade(natural_breaks_agg), cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(natural_breaks_agg, cmap=Elevation, alpha=128, how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regions: Groupby Pixel-Value Connectivity\n",
    "\n",
    "Xarray-spatial's `regions` function creates a raster with unique regions based on connected pixel areas with the same value. Regions assigns each such area a unique integer value.\n",
    "\n",
    "For our reclassified rasters, this can easily be applied to separate out all of the equal-elevation 'plateaus' and name each one.\n",
    "\n",
    "Connectivity can be set to either 4 or 8-pixel neighborhoods, with the default being 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial.zonal import regions\n",
    "from xrspatial import hillshade\n",
    "from datashader.colors import Set1, Set3\n",
    "\n",
    "regions_agg = regions(equal_interval_agg, neighborhood=4)\n",
    "\n",
    "stack(shade(hillshade(regions_agg), cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(regions_agg, cmap=Set1, alpha=128, how='eq_hist'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "### References\n",
    "- Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), p. 406.\n",
    "- Making Maps with Noise Functions: https://www.redblobgames.com/maps/terrain-from-noise/\n",
    "- How Aspect Works: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-aspect-works.htm#ESRI_SECTION1_4198691F8852475A9F4BC71246579FAA"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}