Digital elevation models (DEM) provide a representaion of surface topography (elevation) in two dimensional space. DEMs are a 3D representaion of a terrain’s surface such as the Earth. Typically, DEM data can be represented as a raster which is most easily expressed as being a 2D array with each individual cell having an elevation associated with it.
With a DEM, you can analyze certain derived quantities using prebuild software. This tutorial will go over how to compute the slope and aspect of an area using DEM data.
- Read in DEM data
- Compute the slope for the entirety of the data
- Compute the aspect for the entirety of the data
- Plot the slope of the terrain
- Plot the aspect of the terrain
import os import numpy as np import matplotlib import matplotlib.pyplot as plt import elevation import richdem as rd
matplotlib.rcParams['figure.figsize'] = (8, 5.5)
Download a digital elevation model
elevation module provides access to elevation data from the NASA SRTM mission, with areas of interest specified via a bounding box. Here you draw a bounding box around Mt. Shasta in northern California, and use the clip funciton to output a GeoTIFF.
dem_path = os.path.join(os.getcwd(), 'Shasta-30m-DEM.tif') elevation.clip(bounds=(-122.4, 41.2, -122.1, 41.5), output=dem_path)
Load the DEM with richdem
richdem package has a
LoadGDAL function for loading DEMs that you can use here.
shasta_dem = rd.LoadGDAL(dem_path)
plt.imshow(shasta_dem, interpolation='none') plt.colorbar() plt.show()
Compute and plot slope and aspect
You can use the
rd.TerrainAttribute function to compute slope and aspect for each pixel. Note that there are multiple ways to represent the slope values. Read the richdem docs for more options.
To visualize slope and aspect, you can use the
slope = rd.TerrainAttribute(shasta_dem, attrib='slope_riserun') rd.rdShow(slope, axes=False, cmap='magma', figsize=(8, 5.5)) plt.show()
aspect = rd.TerrainAttribute(shasta_dem, attrib='aspect') rd.rdShow(aspect, axes=False, cmap='jet', figsize=(8, 5.5)) plt.show()