Whitebox Workflows for Python (WbW) Tutorial 2: Geomorphometric Analysis

Introduction

This tutorial is part of a series that includes:

• Tutorial 1: Hydrological Analysis
• Tutorial 2: Geomorphometric Analysis
• Tutorial 3: Mapping Building Footprints from LiDAR

This tutorial is intended to demonstrate how to use the WbW geospatial analysis library to perform a geomorphometric analysis. Geomorphometry is the field focused on understanding landscape processes using digital topography, i.e., digital elevation models (DEMs). This tutorial will not cover all of the functionality related to geomorphometry contained within WbW. For more information, you may refer to the user manual. You may download a copy of the raw Jupyter Notebooks file (*.ipynb) from here.

Setting up WbW

If you haven't already done so, install the Whitebox Worfklows for Python pip package.

pip install whitebox-workflows

We need to import the whitebox_workflows library into our script and set up our floating license ID, which will be used by WbW. Once you register your WbW license, you will be emailed your unique floating license ID, which will likely be a randomly generated three-word phrase, involving an animal. The license below was used for the Geomorphometry 2023 conference in Iasi, Romania and will be valid until August 2023. After that point, you will need to purchase a license for eithter WbW (about 10USD) or WbW-Pro (about 350USD) to use the scripts below. Licenses can be purchased from Whitebox Geospatial Inc.

import whitebox_workflows

license_id = 'geomorphometry-2023' # Update this value with your own license ID when this one expires.

wbe = whitebox_workflows.WbEnvironment(license_id)
wbe.verbose = True # Let each of the function calls output to stdout.

print(wbe.version()) # Let's see what version of WbW we're working with

Working with DEM data

Let's create a new script to download some sample data. Here we'll grab the 'peterborough_drumlins' DEM. The script below will download the data for us, assign the directory to which these data are downloaded to the WbEnvironment working directory and lastly print this location so we can know where the data are being stored. Notice that it may take a few minutes to download the data. In the event that the download takes more than a few minutes, the connection may timeout and you will receive an error. If this should happen, you may download the dataset directly from here but you will need to update the wbe.working_directory to your download folder.

# Download a sample dataset and set the working directory to the location of these data
wbe.working_directory = whitebox_workflows.download_sample_data('Grand_Junction')
print(f'Data have been stored in: {wbe.working_directory}')

Now let's read in the DEM file and generate a multidirectional hillshade image for visualization. When the script below is complete, open the DEM and the newly created hillshade image in QGIS to familiarize yourself with the data set.

# Read in the DEM file.
dem = wbe.read_raster('DEM.tif') # This DEM file is contained in the downloaded data folder.

# There are some NoData holes in the DEM that we should fill in.
# Notice that we can overwrite the 'dem' object
dem = wbe.fill_missing_data(dem, filter_size=35, exclude_edge_nodata=True)

# create a hillshade image for visualization.
hs = wbe.multidirectional_hillshade(dem, full_360_mode=True)
wbe.write_raster(hs, 'hillshade.tif', compress=True)

What are the characteristics of this DEM raster?

print(f"Rows: {dem.configs.rows}")
print(f"Columns: {dem.configs.columns}")
print(f"Resolution (x direction): {dem.configs.resolution_x}")
print(f"Resolution (y direction): {dem.configs.resolution_y}")
print(f"NoData value: {dem.configs.nodata}")
dem.update_min_max() # Find the raster min/max values
print(f"Min. value: {dem.configs.minimum}")
print(f"Max. value: {dem.configs.maximum}")
print(f"Data type: {dem.configs.data_type}")
print(f"Projection: {dem.configs.projection}")

Extracting land-surface parameters (LSPs)

Now let's extract some common land-surface parameters (LSPs), the basic building blocks of a geomorphometric analysis.

# Slope and aspect are two of the most common LSPs. Notice that we're combining the writing of the
# output raster and the running of the function in one line. We won't reuse the raster objects
# created by each function and are only saving them to file and so this makes sense.
wbe.write_raster(wbe.slope(dem, units="degrees"), 'slope.tif', compress=True)
wbe.write_raster(wbe.aspect(dem), 'aspect.tif', compress=True)

# Surface curvatures describe surface shape
wbe.write_raster(wbe.profile_curvature(dem, log_transform=True), 'prof_curv.tif', compress=True)
wbe.write_raster(wbe.tangential_curvature(dem, log_transform=True), 'tan_curv.tif', compress=True)
wbe.write_raster(wbe.plan_curvature(dem, log_transform=True), 'plan_curv.tif', compress=True)
wbe.write_raster(wbe.minimal_curvature(dem, log_transform=True), 'min_curv.tif', compress=True)
wbe.write_raster(wbe.maximal_curvature(dem, log_transform=True), 'max_curv.tif', compress=True)
wbe.write_raster(wbe.mean_curvature(dem, log_transform=True), 'mean_curv.tif', compress=True)
wbe.write_raster(wbe.gaussian_curvature(dem, log_transform=True), 'gauss_curv.tif', compress=True)
wbe.write_raster(wbe.total_curvature(dem, log_transform=True), 'total_curv.tif', compress=True)

The following advanced curvatures are found in WbW-Pro. To run the script below, you'll need a valid WbW-Pro license, otherwise you will receive an error.

wbe.write_raster(wbe.accumulation_curvature(dem, log_transform=True), 'accum_curv.tif', compress=True)
wbe.write_raster(wbe.curvedness(dem, log_transform=True), 'curvedness.tif', compress=True)
wbe.write_raster(wbe.difference_curvature(dem, log_transform=True), 'diff_curv.tif', compress=True)
wbe.write_raster(wbe.generating_function(dem, log_transform=True), 'generating_function.tif', compress=True)
wbe.write_raster(wbe.horizontal_excess_curvature(dem, log_transform=True), 'horizontal_excess_curv.tif', compress=True)
wbe.write_raster(wbe.ring_curvature(dem, log_transform=True), 'ring_curv.tif', compress=True)
wbe.write_raster(wbe.rotor(dem, log_transform=True), 'rotor.tif', compress=True)
wbe.write_raster(wbe.shape_index(dem), 'shape_index.tif', compress=True)
wbe.write_raster(wbe.unsphericity(dem, log_transform=True), 'unsphericity.tif', compress=True)
wbe.write_raster(wbe.vertical_excess_curvature(dem, log_transform=True), 'vertical_excess_curv.tif', compress=True)

The following LSPs can be used to characterize surface roughness and complexity.

wbe.write_raster(wbe.circular_variance_of_aspect(dem, filter_size = 21), 'circular_variance_of_aspect.tif', compress=True)
wbe.write_raster(wbe.edge_density(dem, filter_size=21, normal_diff_threshold=5.0), 'edge_density.tif', compress=True)
wbe.write_raster(wbe.spherical_std_dev_of_normals(dem, filter_size = 21), 'spherical_sd_norms.tif', compress=True)
wbe.write_raster(wbe.standard_deviation_of_slope(dem, filter_size = 21), 'stdev_slope.tif', compress=True)
wbe.write_raster(wbe.surface_area_ratio(dem), 'surface_area_ratio.tif', compress=True)
wbe.write_raster(wbe.ruggedness_index(dem), 'ruggedness_index.tif', compress=True)

Most of the measures of surface roughness and complexity above are measured for local neighbourhoods of a specified size (filter_size = 21). You may modify the filter_size parameter to see the impact of changing the scale of analysis on the output spatial distributions. Note that the grid resolution of the DEM is 5 m. Later we'll explore an alternative method for evaluating the multiscale nature of LSPs.

Measures of local topographic position measure how elevated or low-lying a site is relative to it's neighbouring landscape.

wbe.write_raster(wbe.deviation_from_mean_elevation(dem, filter_size_x=21, filter_size_y=21), 'dev.tif', compress=True)
wbe.write_raster(wbe.difference_from_mean_elevation(dem, filter_size_x=21, filter_size_y=21), 'diff.tif', compress=True)
wbe.write_raster(wbe.elevation_percentile(dem, filter_size_x=21, filter_size_y=21), 'ep.tif', compress=True)
wbe.write_raster(wbe.percent_elev_range(dem, filter_size_x=21, filter_size_y=21), 'percent_elev_range.tif', compress=True)

It can be insightful to evaluate the relations between slope and aspect and slope and elevation, as well as to perform a hypsometric analysis (i.e. area-elevation relation).

wbe.slope_vs_aspect_plot(dem, output_html_file='slope_v_aspect.html', aspect_bin_size=2.0, min_slope=0.1)
wbe.slope_vs_elev_plot(dem_rasters=[dem], output_html_file='slope_v_z.html')
wbe.hypsometric_analysis(dem_rasters=[dem], output_html_file='hypsomeric_analysis.html')

Geomorphons land classes can be useful for interpreting landscape structure. The geomorphons classes are as follows:

Value	Landform Type
1	Flat
2	Peak (summit)
3	Ridge
4	Shoulder
5	Spur (convex)
6	Slope
7	Hollow (concave)
8	Footslope
9	Valley
10	Pit (depression)

wbe.write_raster(
    wbe.geomorphons(
        dem, 
        search_distance=100, 
        flatness_threshold=0.0, 
        flatness_distance=0, 
        skip_distance=0, 
        output_forms=True, 
        analyze_residuals=False
    ),
    'geomorphons.tif',
    compress=True
)

Multiscale geomorphometric analysis

Many of the LSPs calculated above are scale dependent and require us to specify a measurement scale, usually in the form of a filter_size parameter. This parameter allows us to calculate these parameters at a single, uniform spatial scale. For many LSPs, Whitebox also allows us to calculate scale mosiacs. An LSP scale mosaic are calculated by estimating the normalized value of the LSP across a range of spatial scales, known as a scale space or stack, and then to identify the scale at which each grid cell is most expressive. This is known as the characteristic scale or key scale and it can be different for different locations. An LSP scale mosaic therefore represents the value of that LSP at the key scale of each individual grid cell. In comparison to the uniform scale approach, scale mosaics provide locally scale optimized representations of the LSP being measured.

Let's begin by downloading a sample DEM for an area near Peterborough Ontario, Canada, that contains numerous drumlin. Drumlins are streamlined hills formed by glaciation. Open the DEM (peterborough_drumlins.tif) and our newly created hillshade image in QGIS once the following script has completed to visualize the landscape.

wbe.working_directory = whitebox_workflows.download_sample_data('peterborough_drumlins')
print(f'Data have been stored in: {wbe.working_directory}')

# Read in the DEM file and create a hillshade image for visualization.
dem = wbe.read_raster('peterborough_drumlins.tif')
hs = wbe.multidirectional_hillshade(dem, full_360_mode=True)
wbe.write_raster(hs, 'hillshade.tif', compress=True)

Now let's create scale mosaics of deviation from mean elevation (DEV), a measure of relative topographic position, at three broadly defined scale ranges, including local, intermediate (meso), and broad scale ranges.

# Let's turn off the verbose mode because these tools are pretty chatty
wbe.verbose = False

print('Calculating the local scale range...')
dev_local, key_scales = wbe.max_elevation_deviation(dem, min_scale=1, max_scale=25, step_size=1)
wbe.write_raster(dev_local, 'dev_multiscale_local.tif')
wbe.write_raster(key_scales, 'key_scales_local.tif')

print('Calculating the intermediate scale range...')
dev_meso, key_scales = wbe.max_elevation_deviation(dem, min_scale=30, max_scale=100, step_size=1)
wbe.write_raster(dev_meso, 'dev_multiscale_meso.tif')
wbe.write_raster(key_scales, 'key_scales_meso.tif')

print('Calculating the broad scale range...')
dev_broad, key_scales = wbe.max_elevation_deviation(dem, min_scale=500, max_scale=1000, step_size=10)
wbe.write_raster(dev_broad, 'dev_multiscale_broad.tif')
wbe.write_raster(key_scales, 'key_scales_broad.tif')

print('Calculating the multiscale topographic position image...')
mstp = wbe.multiscale_topographic_position_image(dev_local, dev_meso, dev_broad)
wbe.write_raster(mstp, 'mstp.tif')

print('Calculating the full scale range...')
dev_broad, key_scales = wbe.max_elevation_deviation(dem, min_scale=1, max_scale=1000, step_size=2)
wbe.write_raster(dev_broad, 'dev_multiscale_full.tif')
wbe.write_raster(key_scales, 'key_scales_full.tif')

wbe.verbose = True

The multiscale topographic position (MSTP) image is really only meant for visualization purposes, but can be very effective for interpreting landscapes. In effect, with a MSTP image, you are using colour to represent spatial scales. Pixels that are blue are most deviated (either elevated or low-lying) at the local scale, green pixels are most deviated at the intermediate scale, and redish pixels are most deviated at the broadest tested scale range. Of course, you can also have combinations therein, e.g. a yellow pixel is one that is deviated at the intermediate and broad scales, but not particularly deviated at the local scale.

While the MSTP image is useful for visualization and interpretation, ultimately it is the dev_multiscale images that are most useful as modelling inputs (predictors). For example, compare the density of information about landforms contained within the dev_multiscale_full.tif to that of the uniform scale dev.tif generated above. Notice how much more information is contain within this DEV scale mosaic.

Whitebox contains functions for calculating other multiscale LSPs as well, including multiscale_elevation_percentile, multiscale_roughness, multiscale_std_dev_normals, max_anisotropy_dev, gaussian_scale_space, and multiscale_curvatures. Below we use the multiscale_curvatures function (requires a WbW-Pro license) to estimate a multiscale version of mean curvature. Compare this output with that of the mean_curv.tif raster generated previously.

Also notice that in addition to the scale mosaic, each of the multiscale LSP functions also output a second raster for the key scales (i.e. the characteristic scale at which the 'optimal' LSP value is calculated). You may wish to display some of our generated key scale rasters and explore the information contained within these ancillary rasters.

curv_mosaic, key_scales = wbe.multiscale_curvatures(
    dem, 
    curv_type = 'MeanCurvature', 
    min_scale = 5, 
    step_size = 1, 
    num_steps = 10, 
    step_nonlinearity = 1.0,
    log_transform = True
  )
wbe.write_raster(curv_mosaic, 'ms_mean_curv.tif')
wbe.write_raster(key_scales, 'ms_mean_curv_key_scales.tif')

Wrapping things up

Don't forget to check your license in after you're done using it. If you skip this step, your checked-out license won't be returned to the license pool for others to use for one hour. Be kind to others!

print(wbe.check_in_license(license_id))