examples/user_guide/3_Zonal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xarray-spatial\n",
    "### User Guide: Zonal\n",
    "-----\n",
    "\n",
    "Xarray-spatial's zonal functions provide an easy way to generate statistics for zones within a raster aggregate. It's set up with a default set of calculations, or you can input any set of custom calculations you'd like to perform.\n",
    "\n",
    "[Generate terrain](#Generate-Terrain-Data)\n",
    "[Zonal Statistics](#Zonal-Statistics)\n",
    "\n",
    "-----------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's use datashader to render our images...\n",
    "\n",
    "To get started we'll bring in numpy and some functions from datashader to easily render our images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Terrain Data\n",
    "\n",
    "The rest of the geo-related functions focus on raster data or data that's been aggregates into the row-column grid of cells for an image raster. To demonstrate using these raster-based functions, we'll first use xarray-spatial's generate_terrain to generate a fake elevation terrain raster. We use datashader's Canvas as a quick base to set up a new raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import generate_terrain\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "template_terrain = xr.DataArray(np.zeros((H, W)))\n",
    "x_range=(-20, 20)\n",
    "y_range=(-20, 20)\n",
    "\n",
    "terrain = generate_terrain(template_terrain, x_range=x_range, y_range=y_range)\n",
    "\n",
    "shade(terrain, cmap=['black', 'white'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also apply datashader's Elevation colormap imported above to give a more intuitive terrain image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shade(terrain, cmap=Elevation, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zonal Statistics\n",
    "\n",
    "Zonal statistics calculates summary statistics for specific areas or *zones* within an xarray.DataArray aggregate. Specific zones within an aggregate are defined by creating a corresponding aggregate of the same shape and setting the value at each cell to a unique non-zero integer representing a unique zone id. \n",
    "\n",
    "For example, if we set all the values in the top row of the zones aggregate to 3 and apply this to the original values aggregate, zonal stats will calculate statisitics for all the values in the corresponding top row of the values aggregate and return the results as stats for zone \\#3. \n",
    "\n",
    "The output of zonal stats is in the form of a pandas DataFrame, with a row for each zone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set up an example.\n",
    "\n",
    "Imagine you go on a six-day hike.\n",
    "- We can represent the area with a terrain raster.\n",
    "- In that terrain, we can represent each day's path as a line segment from your start to finish point.\n",
    "- We can set this up with a pandas dataframe containing the start and finish points, which we then aggregate with Canvas.line. \n",
    "\n",
    "Let's take a look at these line zones overlayed on the fake terrain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import hillshade\n",
    "from datashader.colors import Set1\n",
    "\n",
    "terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')\n",
    "\n",
    "illuminated = hillshade(terrain)\n",
    "illuminated_shaded = shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')\n",
    "\n",
    "zone_df = pd.DataFrame({\n",
    "   'x': [-11, -5, 4, 12, 14, 18, 19],\n",
    "   'y': [-5, 4, 10, 13, 13, 13, 10],\n",
    "   'trail_segement_id': [11, 12, 13, 14, 15, 16, 17]\n",
    "})\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=x_range, y_range=y_range)\n",
    "\n",
    "zones_agg = cvs.line(zone_df, 'x', 'y', ds.sum('trail_segement_id'))\n",
    "zones_shaded = dynspread(shade(zones_agg, cmap=Set1), max_px=5, threshold=1)\n",
    "\n",
    "stack(illuminated_shaded, terrain_shaded, zones_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can apply zonal stats, after quickly correcting for nan values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import zonal_stats\n",
    "\n",
    "zones_agg.values = np.nan_to_num(zones_agg.values, copy=False).astype(int)\n",
    "zonal_stats(zones_agg, terrain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate custom stats for each zone:\n",
    "\n",
    "We can also put in our own set of stats calculations to perform instead of the default ones above. \n",
    "- We set up a dict with our desired functions and input that as the third argument to `zonal_stats`.\n",
    "- Below, we try out a range function and min and max functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_stats = dict(elevation_change=lambda zone: zone.max() - zone.min(),\n",
    "                    elevation_min=np.min,\n",
    "                    elevation_max=np.max)\n",
    "\n",
    "zonal_stats(zones=zones_agg, values=terrain, stats_funcs=custom_stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the zones are defined by line segments, but they can be any spatial pattern or, more specifically, any region computable as a Datashader aggregate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zonal crosstab\n",
    "\n",
    "Zonal crosstab function can be used when we want to calculate cross-tabulated (categorical stats) areas between two datasets. As an example, assume we have gender data (male/female), and we want to see how all the gender groups distributed over all the zones defined above.\n",
    "\n",
    "First, let's create a mock dataset for gender: 1 for male, and 2 for female. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a mask where 0s are water, and 1s are land areas\n",
    "mask = terrain.data.astype(bool).astype(int)\n",
    "\n",
    "# gender data where 0s are nodata, 1s are male, and 2s are female\n",
    "# assume that the population only takes 10% overall, and the probability of male and female are equal.\n",
    "genders = np.random.choice([0, 1, 2], p=[0.9, 0.05, 0.05], size=zones_agg.shape) * mask\n",
    "genders_agg = xr.DataArray(genders, coords=terrain.coords, dims=terrain.dims)\n",
    "\n",
    "# visualize the results\n",
    "genders_shaded = shade(genders_agg, cmap=['white', 'blue', 'red'])\n",
    "stack(genders_shaded, terrain_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's calculate the stats between the 2 datasets of genders and zones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import zonal_crosstab\n",
    "\n",
    "zonal_crosstab(zones=zones_agg, values=genders_agg, zone_ids=[11, 12, 13, 14, 15, 16], cat_ids=[1, 2])"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}