examples/user_guide/1_Surface.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xarray-spatial\n",
    "### User Guide: Surface tools\n",
    "-----\n",
    "With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.\n",
    "\n",
    "Starting with a raster elevation surface, represented as an Xarray DataArray, these tools can help you identify some specific patterns that may not be readily apparent in the original surface. The return of each function is also an Xarray DataArray.\n",
    "\n",
    "[Hillshade](#Hillshade): Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.\n",
    "\n",
    "[Slope](#Slope): Identifies the slope for each cell of a raster.\n",
    "\n",
    "[Curvature](#Curvature): Calculates the curvature of a raster surface.\n",
    "\n",
    "[Aspect](#Aspect): Derives the aspect for each cell of a raster surface.\n",
    "\n",
    "[Viewshed](#Viewshed): Determines visible locations in the input raster surface from a viewpoint with an optional observer height.\n",
    "\n",
    "-----------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's use datashader to render our images..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll need the basic Numpy and Pandas, as well as datashader,  \n",
    "a data rasterization package highly compatible with Xarray-spatial.  \n",
    "Along with the base package, we'll import several nested functions (shade, stack...)  \n",
    "including Elevation, which we'll use below."
   ]
  },
  {
   "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",
    "\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",
    "The rest of the geo-related functions focus on raster data, i.e. data that's been aggregated into the row-column grid of cells in a raster image. Datashader's Canvas object provides a convenient frame to set up a new raster, so we'll use that with our `generate_terrain` function to generate some fake terrain as an elevation raster. Once we have that, we'll use datashader's shade for easy visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import generate_terrain\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "terrain = xr.DataArray(np.zeros((H, W)))\n",
    "terrain = generate_terrain(terrain)\n",
    "\n",
    "shade(terrain, cmap=['black', 'white'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grayscale values in the image above show elevation, scaled linearly in black-to-white color intensity (with the large black areas indicating low elevation). This shows the data, but it would look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops. Let's try the Elevation colormap we imported above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shade(terrain, cmap=Elevation, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hillshade\n",
    "\n",
    "[Hillshade](https://en.wikipedia.org/wiki/Terrain_cartography) is a technique used to visualize terrain as shaded relief by illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which can be calculated from slope and aspect.\n",
    "\n",
    "Let's apply Hillshade to our terrain and visualize the result with shade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import hillshade\n",
    "\n",
    "illuminated = hillshade(terrain)\n",
    "\n",
    "hillshade_gray_white = shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')\n",
    "hillshade_gray_white"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying hillshade reveals a lot of detail in the 3D shape of the terrain.\n",
    "\n",
    "To add even more detail, we can add the Elevation colormapped terrain from earlier and combine it with the hillshade terrain using datashader's stack function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "terrain_elevation = shade(terrain, cmap=Elevation, alpha=128, how='linear')\n",
    "stack(hillshade_gray_white, terrain_elevation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope\n",
    "[Slope](https://en.wikipedia.org/wiki/Slope) is the inclination of a surface. \n",
    "In geography, *slope* is the amount of change in elevation for an area in a terrain relative to its surroundings.\n",
    "\n",
    "Xarray-spatial's slope function returns the slope at each cell in degrees.\n",
    "Because Xarray-spatial is integrated with Xarray and Numpy, we can apply standard Numpy filters. For example, we can highlight only slopes in the [avalanche risk](https://www.gravityprotection.co.uk/blog/slope-steepness-avalanche-risk.html) range of 25 - 50 degrees. (Note the use of risky.data since these are DataArrays).\n",
    "Stacking the resulting raster with the hillshaded and plain terrain ones from above gives an image with areas of avalanche risk neatly highlighted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import slope\n",
    "\n",
    "risky = slope(terrain)\n",
    "risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)\n",
    "\n",
    "stack(shade(terrain,     cmap=['black', 'white'], how='linear'),\n",
    "      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(risky,       cmap='red',              how='linear', alpha=200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Curvature\n",
    "[Curvature](https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/curvature.htm) is the second derivative of a surface's elevation, or the *slope-of-the-slope*; in other words, how fast the slope is increasing or decreasing as we move along a surface.\n",
    "- A positive curvature means the surface is curving up (upwardly convex) at that cell. \n",
    "- A negative curvature means the surface is curving down (downwardly convex) at that cell. \n",
    "- A curvature of 0 means the surface is striaght and constant in whatever angle it's sloped towards.\n",
    "\n",
    "The Xarray-spatial curvature function returns a raster in units one hundredth (1/100) of the z-factor, or scaling factor (which you can set explicitly in generate _terrain as \"zfactor\").  \n",
    "Reasonably expected values in the curvature raster for a hilly area (moderate relief) would be between -0.5 and 0.5, while for steep, rugged mountains (extreme relief) these can range as far as -4 and 4. For certain raster surfaces it is possible to go even larger than that.\n",
    "\n",
    "Let's generate a terrain with an appropriate z-factor and apply the curvature function to it. Then, we can apply some Numpy filtering (remember, we have access to all those functions) to highlight steeper and gentler curves in the slopes.\n",
    "Stacking these with the hillshaded and plain terrains gives us a fuller picture of the slopes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import curvature\n",
    "\n",
    "terrain_z_one = xr.DataArray(np.zeros((H, W)))\n",
    "terrain_z_one = generate_terrain(terrain_z_one, zfactor=1)\n",
    "curv = curvature(terrain_z_one)\n",
    "curv_hi, curv_low = curv.copy(), curv.copy()\n",
    "curv_hi.data = np.where(np.logical_and(curv_hi.data > 1, curv_hi.data < 4), 1, np.nan)\n",
    "curv_low.data = np.where(np.logical_and(curv_low.data > 0.5, curv_low.data < 1), 1, np.nan)\n",
    "\n",
    "stack(shade(terrain,     cmap=['black', 'white'], how='linear'),\n",
    "      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(curv_hi,     cmap='red',              how='log',    alpha=200),\n",
    "      shade(curv_low,    cmap='green',            how='log',    alpha=200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aspect\n",
    "\n",
    "[Aspect](https://en.wikipedia.org/wiki/Aspect_(geography)) is the orientation of a slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.\n",
    "\n",
    "The Xarray-spatial aspect function returns the aspect in degrees for each cell in an elevation terrain.\n",
    "\n",
    "We can apply aspect to our terrain, then use Numpy to filter out only slopes facing close to North. Then, we can stack that with the hillshaded and plain terrains.\n",
    "(Note: the printout images are from a North point-of-view.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import aspect\n",
    "\n",
    "north_faces = aspect(terrain)\n",
    "north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,\n",
    "                                          north_faces.data < 10), 1, np.nan)\n",
    "\n",
    "stack(shade(terrain,     cmap=['black', 'white'], how='linear'),\n",
    "      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(north_faces, cmap=['aqua'],           how='linear', alpha=100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewshed\n",
    "\n",
    "The `xrspatial.viewshed` function operates on a given aggregate to calculate the viewshed (the visible cells in the raster) for a given viewpoint, or *observer location*.  \n",
    "\n",
    "The visibility model is as follows: Two cells are visible to each other if the line of sight that connects their centers is not blocked at any point by another part of the terrain. If the line of sight does not pass through the cell center, elevation is determined using bilinear interpolation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Simple Viewshed Example\n",
    "\n",
    "- The example below creates a datashader aggregate from a 2d normal distribution.\n",
    "- To calculate the viewshed, we need an observer location so we'll set up an aggregate for that as well.\n",
    "- Then, we can visualize all of that with hillshade, shade, and stack.\n",
    "- The observer location is indicated by the orange point in the upper-left of the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import viewshed\n",
    "\n",
    "OBSERVER_X = -12.5\n",
    "OBSERVER_Y = 10\n",
    "\n",
    "canvas = ds.Canvas(plot_width=W, plot_height=H,\n",
    "                   x_range=(-20, 20), y_range=(-20, 20))\n",
    "\n",
    "normal_df = pd.DataFrame({\n",
    "   'x': np.random.normal(.5, 1, 10000000),\n",
    "   'y': np.random.normal(.5, 1, 10000000)\n",
    "})\n",
    "normal_agg = canvas.points(normal_df, 'x', 'y')\n",
    "normal_agg.values = normal_agg.values.astype(\"float64\")\n",
    "normal_shaded = shade(normal_agg)\n",
    "\n",
    "observer_df = pd.DataFrame({'x': [OBSERVER_X], 'y': [OBSERVER_Y]})\n",
    "observer_agg = canvas.points(observer_df, 'x', 'y')\n",
    "observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),\n",
    "                            threshold=1, max_px=4)\n",
    "\n",
    "normal_illuminated = hillshade(normal_agg)\n",
    "normal_illuminated_shaded = shade(normal_illuminated, cmap=['black', 'white'], \n",
    "                                  alpha=128, how='linear')\n",
    "\n",
    "stack(normal_illuminated_shaded, observer_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Calculate viewshed using the observer location\n",
    "\n",
    "Now we can apply viewshed to the normal_agg, with the observer_agg for the viewpoint. We can then visualize it and stack it with the hillshade and observer rasters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Will take some time to run...\n",
    "%time view = viewshed(normal_agg, x=OBSERVER_X, y=OBSERVER_Y)\n",
    "\n",
    "view_shaded = shade(view, cmap=['white', 'red'], alpha=128, how='linear')\n",
    "\n",
    "stack(normal_illuminated_shaded, observer_shaded, view_shaded)                         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the image highlights in red all points visible from the observer location marked with the orange dot. As one might expect, the areas behind the normal distribution *mountain* are blocked from the viewer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Viewshed on Terrain\n",
    "Now we can try using viewshed on our more complicated terrain.\n",
    "\n",
    "- We'll set up our terrain aggregate and apply hillshade and shade for easy visualization.\n",
    "- We'll also set up an observer location aggregate, setting the location to the center, at (x, y) = (0, 0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrspatial import viewshed\n",
    "\n",
    "x_range=(-20e6, 20e6)\n",
    "y_range=(-20e6, 20e6)\n",
    "terrain = xr.DataArray(np.zeros((H, W)))\n",
    "terrain = generate_terrain(terrain, x_range=x_range, y_range=y_range)\n",
    "terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')\n",
    "\n",
    "illuminated = hillshade(terrain)\n",
    "\n",
    "OBSERVER_X = 0.0\n",
    "OBSERVER_Y = 0.0\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=x_range, y_range=y_range)\n",
    "observer_df = pd.DataFrame({'x': [OBSERVER_X],'y': [OBSERVER_Y]})\n",
    "observer_agg = cvs.points(observer_df, 'x', 'y')\n",
    "observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),\n",
    "                            threshold=1, max_px=4)\n",
    "\n",
    "stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      terrain_shaded,\n",
    "      observer_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can apply viewshed.\n",
    "- Notice the use of the `observer_elev` argument, which is the height of the observer above the terrain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time view = viewshed(terrain, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=100)\n",
    "\n",
    "view_shaded = shade(view, cmap='fuchsia', how='linear')\n",
    "stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      terrain_shaded,\n",
    "      view_shaded,\n",
    "      observer_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fuchsia areas are those visible to an observer of the given height at the indicated orange location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "### References\n",
    "- An overview of the Surface toolset: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/an-overview-of-the-surface-tools.htm\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
}