Surface

With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.

Starting with a raster elevation surface that represented as an Xarray DataArray, these tools help you in identifying some specific patterns that were not readily apparent in the original surface. Return of each function is also an Xarray DataArray.

• Hillshade: Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.

• Slope: Identifies the slope from each cell of a raster.

• Curvature: Calculates the curvature of a raster surface.

• Aspect: Derives the aspect from each cell of a raster surface.

• Viewshed: Determines visible locations in the input raster surface from a viewpoint with some optional observer features.

Importing Packages

import numpy as np
import pandas as pd
import xarray as xr

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

from xrspatial import generate_terrain

W = 800
H = 600

terrain = xr.DataArray(np.zeros((H, W)))
terrain = generate_terrain(terrain)

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

/home/giancastro/miniconda3/envs/xrspatial/lib/python3.8/site-packages/numba/np/ufunc/parallel.py:365: NumbaWarning: The TBB threading layer requires TBB version 2019.5 or later i.e., TBB_INTERFACE_VERSION >= 11005. Found TBB_INTERFACE_VERSION = 9107. The TBB threading layer is disabled.
  warnings.warn(problem)

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:

shade(terrain, cmap=Elevation, how="linear")

Hillshade

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

Let’s apply Hillshade to our terrain and visualize the result with shade.

from xrspatial import hillshade

illuminated = hillshade(terrain)

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

Applying hillshade reveals a lot of detail in the 3D shape of the terrain.

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.

terrain_elevation = shade(terrain, cmap=Elevation, alpha=128, how="linear")
stack(hillshade_gray_white, terrain_elevation)