Importing Packages

To get started, we’ll import numpy and xarray-spatial, along with datashader and a set of its functions to help us quickly render images.

import numpy as np
import datashader as ds
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.transfer_functions import dynspread
from datashader.transfer_functions import set_background
from datashader.colors import Elevation

import xrspatial

Generate Terrain Data

The rest of the geo-related functions work with raster data - or rasterized data, i.e. data that’s been aggregated into the regular row-column grid pattern of a raster. In the code below, the datashader Canvas object helps us do this aggregation and returns a raster in the format of an xarray DataArray.

To demonstrate using these raster-based functions, let’s generate some fake terrain as an elevation raster (or digital elevation model (dem):

from xrspatial import generate_terrain

W = 800
H = 600

cvs = ds.Canvas(
    plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6)
terrain = generate_terrain(canvas=cvs)

shade(terrain, cmap=["black", "white"], how="linear")

The grayscale values above show the elevation, scaled linearly in intensity from dark to bright white (with the large black areas being low elevation). This gives us a glimpse of a fair amount of detail, but we could make it more intuitive by shading it like a landscape.

  • Datashader’s Set1 maps low values to colors representing water and high ones to colors representing mountaintops, with a range of landscape color in between.

  • We’ll apply this, but first we’ll also apply xarray-spatial’s hillshade to give an illuminated representation.

from xrspatial import hillshade
from datashader.colors import Set1
import pandas as pd

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))

terrain = generate_terrain(canvas=cvs)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how="linear")

illuminated = hillshade(terrain)
illuminated_shaded = shade(illuminated, cmap=["gray", "white"], alpha=255, how="linear")

stack(illuminated_shaded, terrain_shaded)

Focal Statistics and Convolutions

Similar to zonal statistics, focal statistics are used to calculate local statistics, but relative to a focal point and a given neighborhood. The neighborhood is a user-defined kernel representing the cells which should be used in calculations. Currently, only circle and annulus kernels are implemented, but any custom kernel can be used given that: - The kernel is a numpy array - The kernel’s dimensions are odd, so a 3x1 kernel is valid while a 3x2 is not. This is required for symmetry around the focal point in the current implementation.

The following example uses focal statistics to calculate the topographic position index (TPI), a measure of local topographic position relative to nearby neighbors. The TPI is scale-dependent and will vary based on the inner and outer radii chosen for the annulus kernel used for this calculation. Once calculated, the TPI can be used to classify slope positions and landforms within landscapes, and can also be used as a numeric feature for model inputs.

  • First, we use focal’s calc_cellsize to get the scaling right.

  • Next, we multiply by a larger and smaller integer for the outer and inner annulus radii, respectively.

  • Then, we generate the kernel with annulus_kernel.

  • And, finally, we use that kernel with focal.apply to calculate the TPI at each point in the terrain.

  • To visualize all this, we can shade the TPI and stack it in a composite image with the original terrain.

from xrspatial import convolution
from xrspatial import focal

cellsize_x, cellsize_y = convolution.calc_cellsize(terrain)
# Use an annulus kernel with a ring at a distance from 25-30 cells away from focal point
outer_radius = cellsize_x * 30
inner_radius = cellsize_x * 25
kernel = convolution.annulus_kernel(cellsize_x, cellsize_y, outer_radius, inner_radius)
tpi = terrain - focal.apply(terrain, kernel)
tpi_terrain = hillshade(terrain - focal.apply(terrain, kernel))
tpi_terrain_shaded = shade(
    tpi_terrain, cmap=["white", "black"], alpha=255, how="linear"
stack(illuminated_shaded, tpi_terrain_shaded)