examples/user_guide/2_Proximity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xarray-spatial\n",
    "### User Guide: Proximity (Distance) tools\n",
    "-----\n",
    "The Proximity tools let you perform analyses using one of three distance metrics: Euclidean (straight-line), Manhattan (with path obstacles, i.e., taxicab distance), and Great Circle (distance on the surface of a sphere).\n",
    "\n",
    "[Proximity Distance](#Proximity-Distance): Calculates the distance to the nearest of a set of target - or source - points for each point in the input raster.\n",
    "\n",
    "[Proximity Allocation](#Proximity-Allocation): Identifies the nearest source (target) point (the 'allocation' point) for each cell in the input raster.\n",
    "\n",
    "[Proximity Direction](#Proximity-Direction): Returns the direction to the nearest source point (the 'allocation') for each cell in the input raster.\n",
    "\n",
    "-----------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's use datashader to render our images...\n",
    "\n",
    "We'll start by importing some basic packages, including several functions from datashader that will help us render our images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\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": [
    "## Proximity Distance\n",
    "\n",
    "The `xrspatial.proximity` function operates on a given values aggregate to produce a new distance aggregate based on a set of target values and using a distance metric. The metric is used to calculate the distance at each each array cell (pixel) in the values aggregate from the nearest of the target points. For `proximity`, this smallest distance is set as each cell's value in the returned aggregate.\n",
    "\n",
    "A powerful feature of `proximity` is that you can target specific points, the `target_values` parameter, in the aggregate for distance calculation and ignore the rest. If this isn't set, then all non-zero pixels are set as targets. Feel free to play with the parameter and see the difference between using `target_values=[1,2,3,4]` vs. the default option. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set up a values aggregate raster and try out the proximity function.\n",
    "\n",
    "- We can set up a regular cell-grid aggregate with datashader's `Canvas.points`.\n",
    "- We'll set this one up in reverse to make it easier.\n",
    "- We define the target points in a pandas DataFrame\n",
    "- We aggregate those with Canvas.points into an aggregates raster, which embeds the target points into a large grid of nan's.\n",
    "- We change all the non-target points from nan to 0.\n",
    "\n",
    "Now we have an aggregate raster full of zeros but with a few non-zero target/source points at our desired locations and each with a unique 'id' value.\n",
    "\n",
    "We'll set a quick shade and background to visualize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import proximity\n",
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame({\n",
    "   'x': [-13, -11, -5, 4, 9, 11, 18, 6],\n",
    "   'y': [-13, -5, 0, 10, 7, 2, 5, -5],\n",
    "   'id': [1, 2, 3, 4, 5, 6, 7, 8]\n",
    "})\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H,\n",
    "                x_range=(-20, 20), y_range=(-20, 20))\n",
    "\n",
    "points_agg = cvs.points(df, x='x', y='y', agg=ds.min('id'))\n",
    "points_agg.data[~np.isfinite(points_agg.data)] = 0\n",
    "\n",
    "points_shaded = dynspread(shade(points_agg, cmap='salmon', min_alpha=0, span=(0,1), how='linear'),\n",
    "                          threshold=1,\n",
    "                          max_px=5)\n",
    "set_background(points_shaded, 'black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply Proximity:\n",
    "Now we can apply the proximity function to our points aggregate.\n",
    "   \n",
    "##### Set Targets:\n",
    "We can set the target_values parameter to a subset of the points, for example, \\[1, 2, 3, 4\\], to only measure distance everywhere from the nearest of these four points.\n",
    "\n",
    "Let's do that, and then shade it for visualization.\n",
    "\n",
    "Note: feel free to remove or add to the targets from among the digits 1-8 and see how the image changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = [1, 2, 3, 4]\n",
    "target_proximity_agg = proximity(points_agg, target_values=targets, distance_metric=\"GREAT_CIRCLE\")\n",
    "\n",
    "stack(shade(target_proximity_agg, cmap=['#bbeb9e', 'black'], how='eq_hist'),\n",
    "      points_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the brighter areas nearer the target points and the black lines along areas that are equal distance to two points. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Proximity with Default Targets: All Non-zero Points\n",
    "Leaving out the target_values parameter entirely defaults it to all non-zero points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proximity_agg = proximity(points_agg, distance_metric=\"GREAT_CIRCLE\")\n",
    "\n",
    "stack(shade(proximity_agg, cmap=['#bbeb9e', 'black'], how='eq_hist'),\n",
    "      points_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply proximity to a line:\n",
    "Proximity can be applied to any shape appregate, so applying it to a raster aggregate of lines works as well. \n",
    "\n",
    "For a quick example we can once again aggregate the points DataFrame from above, but, using Canvas.line, into lines rather than isolated points.\n",
    "\n",
    "Take a look below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_agg = cvs.line(df, x='x', y='y')\n",
    "line_shaded = dynspread(shade(line_agg, cmap=['salmon',  'salmon']),\n",
    "                          threshold=1,\n",
    "                          max_px=2)\n",
    "set_background(line_shaded, 'black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Apply `proximity`:\n",
    "Now we can apply proximity and it gives us the smallest distance at every cell from any point on the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_proximity = proximity(line_agg)\n",
    "stack(shade(line_proximity, cmap=['#bbeb9e', 'black'], how='linear'),\n",
    "      line_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Apply Proximity with transformations:\n",
    "Since Xarray-spatial is built on Xarray and Numpy, we can also use proximity and then apply some nifty transformations from the wide range of transformations available through the xarray DataArray API.\n",
    "\n",
    "Here's an example using `DataArray.where()`:\n",
    "- We'll bring in the line_proximity we calculated above.\n",
    "- Then we'll put a minimum and maximum distance clause into `DataArray.where()`.\n",
    "\n",
    "The result is shaded and visualized below. As you can see, this provides a neat ouline around the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "where_clause = (line_proximity > 1) & (line_proximity < 1.1)\n",
    "proximity_shaded = shade(line_proximity.where(where_clause), cmap=['darkturquoise', 'darkturquoise'])\n",
    "proximity_shaded = set_background(proximity_shaded, 'black')\n",
    "stack(proximity_shaded, line_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proximity Allocation\n",
    "\n",
    "Similar to `xrspatial.proximity`, the `xrspatial.allocation` function operates on an aggregate and finds the smallest distance from each cell to any one of the target, or source, points. But, instead of the distance, it returns the value at that source point and sets the current cell to that value. So, the result is a raster with block areas, with all cells in each area set to the value of that area's nearest source.\n",
    "\n",
    "In the same way as for `proximity`, we can set the target points with the `target_values` parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply allocation to our aggregate from above, leaving the default all target points active."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import allocation\n",
    "\n",
    "allocation_agg = allocation(points_agg, distance_metric=\"GREAT_CIRCLE\")\n",
    "                            \n",
    "stack(shade(allocation_agg, cmap=['darkturquoise', 'black'], how='linear'),\n",
    "      points_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the blocks: each of the differently shaded blocks contains all of the points that have the target point in the center as their nearest target point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proximity Direction\n",
    "\n",
    "Also similar to `xrspatial.proximity`, `xrspatial.direction` also starts by calculating the smallest distance for each cell to a source point. But, instead of distance, it returns the direction in degrees to that nearest source point. The output values range from 0 to 360:\n",
    "- 0 is for the source cell itself\n",
    "- 90 is the east\n",
    "- 180 is the south\n",
    "- 270 is the west\n",
    "- 360 is the north\n",
    "\n",
    "Once again, you can target specific values with the `target_values` parameter.\n",
    "\n",
    "Let's apply it to our aggrega from above, leaving all the default target points active.\n",
    "\n",
    "(Note: the image is from a north point-of-view; i.e. north is down.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import direction\n",
    "\n",
    "direction_agg = direction(points_agg, distance_metric=\"GREAT_CIRCLE\")\n",
    "                            \n",
    "stack(shade(direction_agg, cmap=['darkturquoise', 'black'], how='linear'),\n",
    "      points_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "### References\n",
    "An overview of the Distance toolset: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/an-overview-of-the-distance-tools.htm"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}