Focal#

Importing Packages#

To get started, we’ll import numpy and xarray-spatial, along with matplotlib for rendering images.

[1]:

import numpy as np
import xarray as xr

import xrspatial

[2]:

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap, Normalize
import numpy as np

# Rendering helpers (replace datashader shade/stack).
# _shade maps a 2D array through a colormap to an (H, W, 4) uint8 RGBA,
# with NaN pixels transparent. `alpha` (0-255) overrides finite-pixel alpha.
def _shade(arr, cmap=None, alpha=None, how='linear'):
    arr = np.asarray(arr)
    finite = np.isfinite(arr)
    if finite.any():
        lo = np.nanmin(arr)
        hi = np.nanmax(arr)
    else:
        lo = hi = 0.0
    norm = np.zeros_like(arr, dtype=np.float32)
    if hi > lo:
        norm = np.where(finite, (arr - lo) / (hi - lo), 0.0).astype(np.float32)
    if isinstance(cmap, (list, tuple)):
        cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
    elif cmap is None:
        cmap = plt.get_cmap('terrain')
    rgba = cmap(norm, bytes=True)
    a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
    rgba[..., 3] = a
    return rgba

# _stack alpha-composites a sequence of RGBA layers bottom-to-top.
def _stack(*layers):
    base = layers[0].astype(np.float64)
    for top in layers[1:]:
        a = (top[..., 3:4] / 255.0)
        rgb = top[..., :3] * a + base[..., :3] * (1 - a)
        out_a = np.maximum(top[..., 3], base[..., 3])
        base = np.concatenate([rgb, out_a[..., None]], axis=-1)
    return np.clip(base, 0, 255).astype(np.uint8)

def _show(rgba, figsize=(8, 6)):
    plt.figure(figsize=figsize)
    plt.imshow(rgba)
    plt.axis('off')
    plt.show()

# Elevation-style colormap (replaces datashader.colors.Elevation).
Elevation = LinearSegmentedColormap.from_list(
    'elevation',
    [(0.0, (0.0, 0.3, 0.6)), (0.2, (0.2, 0.6, 0.85)), (0.3, (0.55, 0.78, 0.92)),
     (0.32, (0.4, 0.73, 0.45)), (0.45, (0.3, 0.55, 0.25)), (0.6, (0.5, 0.5, 0.3)),
     (0.75, (0.65, 0.55, 0.4)), (0.9, (0.8, 0.8, 0.78)), (1.0, (1.0, 1.0, 1.0))],
    N=256)

# matplotlib's Set1 (replaces datashader.colors.Set1).
Set1 = plt.get_cmap('Set1')

Generate Terrain Data#

The rest of the geo-related functions work with raster data, i.e. data that’s been aggregated into the regular row-column grid pattern of a raster. In the code below, we build a raster grid directly as 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):

[3]:

from xrspatial import generate_terrain

W = 800
H = 600

template_terrain = xr.DataArray(np.zeros((H, W)))
x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)

terrain = generate_terrain(template_terrain, x_range=x_range, y_range=y_range)
terrain.attrs["unit"] = "meter"

_show(_shade(terrain, cmap=["black", "white"]))
../_images/user_guide_focal_7_0.png

The grayscale values above show the elevation, scaled linearly in intensity from black 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.

  • The matplotlib Set1 colormap 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.

[4]:

from xrspatial import hillshade

terrain_shaded = _shade(terrain, cmap=Elevation, alpha=128)

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

_show(_stack(illuminated_shaded, terrain_shaded))
../_images/user_guide_focal_9_0.png

Focal Statistics and Convolutions#

Similar to zonal statistics, focal statistics are also used to calculate statistics locally, but relative to a focal point and a neighborhood around that rather than an inflexible zone. The neighborhood is defined with a kernel representing the neighboring cells which should be used in the calculations for each cell. Currently, only circle and annulus kernels are implemented, but any custom kernel can be used, as long as:

  • The kernel is a numpy array

  • The kernel’s dimensions are odd: for example, 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 (or elevation) relative to nearby neighbors. The TPI is scale-dependent and will vary based on the relative sizes of the inner and outer radii in the annulus kernel used, so one TPI does not define all TPIs. Once calculated, a TPI can be used to classify slope positions and landforms within landscapes, and can also be used as a numeric feature for model inputs.

We calculate a TPI for our terrain as follows:

  • 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.

[5]:

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 = str(cellsize_x * 30) + "m"
inner_radius = str(cellsize_x * 25) + "m"
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)
_show(_stack(illuminated_shaded, tpi_terrain_shaded))
../_images/user_guide_focal_13_0.png

Convolutions#

The focal.apply function can be computationally expensive depending on the sizes of the kernel and image. Additionally, we’d like to extend focal’s capabilities to the use of custom convolution kernels. This is where the convolution functions come in.

Let’s try an example with kernels from image processing. The Sobel operator is a crude, but computationally inexpensive way to do edge-detection. Let’s try setting one up and applying it to our terrain image.

For our example, we’ll set up a horizontal Sobel operator, which calculates an approximation of the derivative in the horizontal dimension to get horizontal edges.

(Note: By default, the convolution module will use the local CUDA-enabled GPU unit if it is available. To use only the CPU, you can pass use_cuda=False to convolution.convolve_2d.)

  • First, we’ll set up a horizontal Sobel kernel manually as a 2D numpy array. Notice the vertical column of zeros and the symmetry on either side of this column.

  • Next, we’ll apply this kernel to the terrain values to generate a corresponding array of Sobel values.

  • Finally, we’ll set those values into a DataArray raster with the proper coordinates, dimensions, and attributes from the terrain.

  • To visualize all this, we’ll apply hillshade and shade to the sobel terrain and stack it in an image with the original.

Notice all the emphasized horizontal edges. (You can also edit the kernel numbers below to make it vertical and see the vertical edges emphasized, instead.)

[6]:

from xrspatial import convolution
from xarray import DataArray

# Use Sobel operator
kernel = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]])
print("Horizontal Sobel Kernel:")
print(kernel)

sobel_values = convolution.convolve_2d(terrain.values, kernel)
sobel = DataArray(
    sobel_values,
    coords=terrain.coords,
    dims=terrain.dims,
    attrs=terrain.attrs,
)
sobel_terrain = hillshade(sobel)
sobel_terrain_shaded = _shade(sobel_terrain, cmap=["white", "black"], alpha=255)
_show(_stack(illuminated_shaded, sobel_terrain_shaded))

Horizontal Sobel Kernel:
[[ 1  0 -1]
 [ 2  0 -2]
 [ 1  0 -1]]
../_images/user_guide_focal_17_1.png