Terrain Algorithms from Scratch

There are plenty of tools to calculate slope, aspect, and hillshading from elevation data, but if you've ever been curious about how they're calculated, this post goes through the process of implementing those algorithms from scratch in Python using just Numpy.

Elevation Data

Terrain algorithms won't do much good without elevation data to apply them to, so the first thing we'll do is to grab some of that. Specifically, we'll download a digital elevation model (DEM) over Mount St. Helens from Microsoft Planetary Computer. Using STAC (via pystac-client) will allow us to quickly find the asset we need, and the stackstac package will make it easy to download and clip it into a 2D Numpy array.

python

from pystac_client import Client import planetary_computer as pc import stackstac # Load the Planetary Computer catalog catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1") # Define the bounding box around Mount St. Helens bbox = (-122.259265, 46.149409, -122.112323, 46.250639) # Search the catalog for Copernicus GLO-30 elevation data over Mount St. Helens search = catalog.search( collections=["cop-dem-glo-30"], bbox=bbox, ) # Sign the search results with Planetary Computer to allow asset access items = pc.sign(search) # Download the data, squeeze out extra dimensions, and grab the Numpy array dem = stackstac.stack(items, bounds_latlon=bbox).squeeze().values dem.shape
(366, 530)

Now that we have a 2D array of elevation values, let's take a look using matplotlib.

python

import matplotlib.pyplot as plt plt.figure(figsize=(10, 6)) plt.imshow(dem, cmap="terrain") plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7f7e7d125580>

png

Slope

The concept of slope is simple: How much does elevation change within an area? Steep areas have lots of elevation change over a short distance while flat areas have very little. However, applying that concept to a 2D array raises the question, how do we define change within an area?

There have been many different slope algorithms created to answer that question, and each produces slightly different results. We're going to use the Horn 1981 algorithm since it is widely used (supported by GDAL, GRASS, Whitebox Tools, ESRI, etc.) but if you want to dive into the rabbit hole of slope algorithm comparisons, check out this site or this paper.

Horn's algorithm calculates slope for each pixel using the following equation, where dzdx\frac{dz}{dx} is the east-west gradient of neighboring pixels and dzdy\frac{dz}{dy} is the north-south gradient.

slopepercent=dzdx2+dzdy2slope _{percent} = \sqrt {\frac{dz}{dx}^2 + \frac{dz}{dy}^2}

Let's break that equation down by calculating the slope of a single pixel.

One Pixel at a Time

Imagine we want to calculate slope for the center pixel of a 3x3 pixel neighborhood ww with the following elevations:

python

import numpy as np w = np.array([ [52, 74, 73], [63, 98, 89], [72, 73, 75], ]) plt.imshow(w, cmap="terrain") for (y, x), z in np.ndenumerate(w): plt.text(x, y, z, ha="center", va="center")

png

East-West Gradient

The first term in the Horn algorithm, the east-west gradient dzdx\frac{dz}{dx}, describes how elevation changes between the east and west side of the pixel neighborhood. To solve it, we'll break it down into the change in elevation dzdz over the horizontal distance dxdx.

The vertical distance between pixels, dzdz, is calculated with the following equation, where the northwest pixel in the neighborhood is labeled znwz_{nw}, the southwest pixel is labeled zswz_{sw}, and so on.

dz=(znw+zsw+2zw)(zne+zse+2ze)8dz = \frac{(z_{nw} + z_{sw} + 2z_{w}) - (z_{ne} + z_{se} + 2z_{e})}{8}

There are a few things to notice in the equation above.

  1. We're calculating the difference between the sum of the western and eastern pixels.
  2. The directly west and east pixels, zwz_{w} and zez_{e}, are multiplied by 2 to increase their weight over the diagonal pixels.
  3. The result is divided by 8 to normalize the value based on the input weights.

Here's dzdz in code:

python

dz = ((w[0][0] + w[1][0] * 2 + w[2][0]) - (w[0][2] + w[1][2] * 2 + w[2][2])) / 8

The horizontal distance between pixels, dxdx, is simply the raster resolution (for this example, 30).

python

dx = 30

Now we can calculate dzdx\frac {dz}{dx}, the change in elevation over the xx dimension, by simply dividing the two terms.

python

dz_dx = dz / dx

North-South Gradient

The second term in the Horn algorithm, the north-south gradient dzdy\frac{dz}{dy}, describes how elevation changes between the north and south side of the pixel neighborhood. It's solution is nearly identical to the east-west gradient, after swapping in the appropriate pixels.

dz=(zsw+zse+2zs)+(znw+zne+2zn)8dz = \frac{(z_{sw} + z_{se} + 2z_{s}) + (z_{nw} + z_{ne} + 2z_{n})}{8}

And in code:

python

dz = ((w[2][0] + w[2][1] * 2 + w[2][2]) - (w[0][0] + w[0][1] * 2 + w[0][2])) / 8

Assuming our image has square pixels, the distance between pixels in the north-south dimension dydy will be the same as dxdx.

python

dy = 30

The last step in calculating the north-south gradient is to divide the change in elevation over the yy dimension by the horizontal distance between cells.

python

dz_dy = dz / dy

Putting It Together

With the east-west and north-south gradients calculated, dzdx\frac{dz}{dx} and dzdy\frac{dz}{dy} respectiely, solving the Horn algorithm is straightforward. Just plug the two solved terms into the original equation.

slopepercent=dzdx2+dzdy2slope _{percent} = \sqrt {\frac{dz}{dx}^2 + \frac{dz}{dy}^2}

python

slope_pct = np.sqrt(dz_dx ** 2 + dz_dy ** 2)

With a little more work, we can convert the percent slope into more familiar degrees of slope.

slopedegrees=arctan(slopepercent)(180π)slope _{degrees} = \arctan \left(slope _{percent}\right) * \left(\frac {180}{\pi}\right)

python

slope = np.arctan(slope_pct) * (180 / np.pi) slope
18.130932944052347

For convenience, let's simplify the code above and package it into a function that will calculate the slope of a single pixel given it's 3x3 neighborhood of pixels.

python

def pixel_slope(w, resolution): dz_dx = ((w[0][0] + w[1][0] * 2 + w[2][0]) - (w[0][2] + w[1][2] * 2 + w[2][2])) / (8 * resolution) dz_dy = ((w[2][0] + w[2][1] * 2 + w[2][2]) - (w[0][0] + w[0][1] * 2 + w[0][2])) / (8 * resolution) return np.arctan(np.sqrt(dz_dx ** 2 + dz_dy ** 2)) * (180 / np.pi) pixel_slope(w, 30)
18.130932944052347

With the fundamentals of Horn's algorithm down, the challenge now is simply to calculate it for each pixel.

Scaling Up

To calculate slope from our elevation data, we'll iterate over each row and column in the DEM, grab the 3x3 window of neighboring pixels, use the pixel_slope function to calculate the slope of the center pixel, and store the result in an empty slope array.

First, we'll create the empty array to store slope values. We'll make it two pixels smaller than the DEM in the x and y dimensions to account for the fact that edge pixels don't have the eight required neighbors.

python

slope = np.empty((dem.shape[0] - 2, dem.shape[1] - 2)) slope.shape
(364, 528)

Now we'll iterate over rows and columns (dropping one pixel from each side to account for edges) and calculate slope for each neighborhood of pixels.

python

for row in range(1, dem.shape[0] - 1): for col in range(1, dem.shape[1] - 1): w = dem[row-1:row+2, col-1:col+2] slope[row-1][col-1] = pixel_slope(w, 30)

Finally, let's see what our slope map looks like, with flat areas in blue and steep areas in red.

python

plt.figure(figsize=(10, 6)) plt.imshow(slope, cmap="RdBu_r") plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7f7e6f70d4f0>

png

Aspect

Aspect is closely related to slope, describing orientation rather than steepness. Since we've already implemented the Horn slope algorithm, we'll use that for calculating aspect as well, with the following equation.

aspect=arctan2(dzdx,dzdy)aspect = \arctan2 \left( \frac{dz}{dx} , \frac{dz}{dy} \right)

The east-west and north-south gradients, dzdx\frac{dz}{dx} and dzdy\frac{dz}{dy} respectiely, are calculated identically to slope. The only difference is that instead of taking the square root of their sum to get the overall slope, we use the arctangent to calculate the angle between them.

Here's that equation in code, plus conversion to degrees and rescaling to compass bearings:

python

def pixel_aspect(w, resolution): """Calculate the aspect of a pixel in degrees given its 3x3 neighborhood `w` and cell resolution.""" dz_dx = ((w[0][0] + w[1][0] * 2 + w[2][0]) - (w[0][2] + w[1][2] * 2 + w[2][2])) / (8 * resolution) dz_dy = ((w[2][0] + w[2][1] * 2 + w[2][2]) - (w[0][0] + w[0][1] * 2 + w[0][2])) / (8 * resolution) aspect = np.arctan2(dz_dy, dz_dx) * (180 / np.pi) # Convert to compass bearings 0 - 360 aspect = 450 - aspect if aspect > 90 else 90 - aspect return aspect pixel_aspect(w, 30)
184.12761236173162

We calculate aspect for each pixel the same way we did slope, by iteratively filling an empty 2D array with aspects calculated from each pixel's 3x3 neighborhood.

python

aspect = np.empty((dem.shape[0] - 2, dem.shape[1] - 2)) for row in range(1, dem.shape[0] - 1): for col in range(1, dem.shape[1] - 1): w = dem[row-1:row+2, col-1:col+2] aspect[row-1][col-1] = pixel_aspect(w, 30) plt.figure(figsize=(10, 6)) plt.imshow(aspect, cmap="magma") plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7f7e6f650ac0>

png

Hillshading

With slope and aspect calculated, generating hillshading to visualize the terrain is simple.

The formula to calculate hillshading is below, with all units in radians. The zenith and azimuth parameters describe the position of the simulated light source, and can be tuned to adjust the hillshading effect.

hillshade=cos(zenith)cos(slope)+sin(zenith)sin(slope)cos(azimuthaspect)hillshade = \cos(zenith) * \cos(slope) + \sin(zenith) * \sin(slope) * \cos(azimuth - aspect)

Here's that equation in code, using the slope and aspect arrays we calculated earlier.

python

azimuth = 315 altitude = 45 # Convert solar altitude to zenith and convert everything to radians zenith_rad = 90 - altitude * np.pi / 180 azimuth_rad = azimuth * np.pi / 180 slope_rad = slope * np.pi / 180 aspect_rad = aspect * np.pi / 180 # Calculate hillshade and scale to 8-bit hs = 255 * (np.cos(zenith_rad) * np.cos(slope_rad) + np.sin(zenith_rad) * np.sin(slope_rad) * np.cos(azimuth_rad - aspect_rad)) # Clip out-of-bounds values hs = np.clip(hs, 0, 255) plt.figure(figsize=(10, 6)) plt.imshow(hs, cmap="Greys_r")
<matplotlib.image.AxesImage at 0x7f7e7d0339d0>

png

Wrapping Up

Now that we know how to implement terrain algorithms from scratch, the next step is to uninstall QGIS, WhiteboxTools, GDAL, and any other geospatial tools we no longer need!

Okay, probably not. There are faster and more convenient ways to generate terrain data than rolling your own implementations, but getting a glimpse at the underlying algorithms does provide some interesting insights into how they work.

python

fig, axs = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 6)) axs[0, 0].set_title("Elevation") axs[0, 0].imshow(dem, cmap="terrain") axs[0, 1].set_title("Hillshade") axs[0, 1].imshow(hs, cmap="Greys_r") axs[1, 0].set_title("Slope") axs[1, 0].imshow(slope, cmap="RdBu_r") axs[1, 1].set_title("Aspect") axs[1, 1].imshow(aspect, cmap="magma")
<matplotlib.image.AxesImage at 0x7f7e6f560d90>

png



© 2023 Aaron Zuspan