16  Using WhiteboxTools in R

16.1 Topics Covered

1. Access external geospatial software from R
2. Calculate terrain derivatives using WhiteboxTools
3. Perform vector geoprocessing using WhiteboxTools
4. Process lidar point clouds using WhiteboxTools

16.2 Introduction

R has a large and growing collection of native packages for geospatial data handling and analysis, such as sf for vector data, terra for raster data, and lidR for lidar point clouds. However, many mature geospatial software tools have been developed outside of R, and some of them implement algorithms that are not available, or are less efficient, in native R packages. Instead of re-implementing these tools, a number of R packages act as bridges that allow R to call functions from external geospatial software. This allows the user to combine the analytical, statistical, and visualization capabilities of R with the specialized tools available in other environments, all within a single, reproducible script.

Below are some commonly used packages that link R to other geospatial software:

• whitebox: provides an R interface to WhiteboxTools, an open-source geospatial analysis platform with more than 450 tools for terrain analysis, hydrology, image processing, vector geoprocessing, and lidar processing. This is the package that we explore in this chapter.
• rgrass: provides an interface to GRASS GIS, allowing R to run GRASS modules and transfer raster and vector data between R and a GRASS database.
• qgisprocess: provides access to the QGIS processing framework. Since QGIS processing includes native QGIS algorithms as well as providers such as GDAL, GRASS, and SAGA, this package makes a very large number of tools available from R.
• Rsagacmd and RSAGA: provide interfaces to SAGA GIS, which is especially well known for terrain analysis tools.
• link2GI: helps locate and configure installations of external software, such as GRASS, SAGA, and Orfeo Toolbox, so that they can be called from R. The OTBsegm package provides an R interface specifically for Orfeo Toolbox segmentation algorithms.
• gdalUtilities and the gdal_utils() function from sf: provide access to the GDAL command line utilities, such as gdalwarp, gdal_translate, and ogr2ogr.
• arcgisbinding: part of the R-ArcGIS Bridge, which allows R to read and write data from ArcGIS Pro and to be used for building geoprocessing tools within ArcGIS.
• rgee: provides an R interface to the Google Earth Engine cloud computing platform.
• reticulate: although not geospatial-specific, this package allows R to call Python code and libraries, which opens up the large Python geospatial ecosystem (e.g., ArcPy, rasterio, and geopandas).

Most of these packages work by passing commands to the external software, which then reads input from, and writes output to, files on disk. As a result, a common workflow is to (1) write data to disk from R, (2) run the external tool, and (3) read the results back into R for further analysis or visualization. We will follow this pattern throughout this chapter.

Using external software requires that it be installed on your machine and that R can find it. Some tools, such as GRASS GIS, QGIS, and SAGA GIS, must be installed separately. WhiteboxTools is especially easy to set up since the whitebox package can download the WhiteboxTools executable for you.

Since this is a short chapter meant to simply introduce how R can interface with other geospatial software tools, it does not include any questions or exercises.

16.3 Setup

The terra package is used for reading and processing raster data while sf is used for vector data. Maps are generated with tmap and tmaptools. The tidyverse is loaded for general data handling; specifically, we use the str_glue() function from stringr to define file paths. Lastly, the whitebox package provides access to WhiteboxTools.

library(terra)
library(tidyverse)
library(tmap)
library(tmaptools)
library(sf)
library(whitebox)

The whitebox package does not include the WhiteboxTools software itself. If this is your first time using the package, you will need to download the WhiteboxTools executable using install_whitebox(). This only needs to be done once. Once installed, wbt_init() initializes WhiteboxTools by locating the executable and confirming that it can be called from R.

# Only needs to be run once
install_whitebox()

wbt_init()

All WhiteboxTools functions in the whitebox package are prefixed with wbt_. A complete list of available tools and their parameters are available in the WhiteboxTools User Manual. You can also use wbt_list_tools() to list the available tools and wbt_tool_help() to view documentation for a specific tool. It can also be helpful to reference the whitebox package documentation, which is available here.

16.4 Terrain Analysis

We begin with a terrain analysis example. First, we read in a digital terrain model (DTM) using terra. To reduce processing time, we aggregate the DTM to a coarser cell size by a factor of 5 using the mean of the cells being combined. Since WhiteboxTools operates on files rather than R objects, we must save the aggregated DTM to disk using writeRaster(). Note that writeRaster() also returns the raster, so we can save it to a variable at the same time.

fldrPth <- "gslrData/chpt16/data/"
outFld <- "gslrData/chpt16/output/"

dtm <- rast(str_glue("{fldrPth}dtm.tif")) |>
  aggregate(fact=5, fun=mean) |>
  writeRaster(str_glue("{outFld}dtm5.tif"), overwrite=TRUE)

|---------|---------|---------|---------|
=========================================
                                          

The aggregated DTM is then mapped using tmap with a reversed grey color ramp so that higher elevations are displayed with lighter shades.

tm_shape(dtm) +
  tm_raster(col.scale = tm_scale_continuous(values="-brewer.greys"))

Topographic slope is calculated using wbt_slope(). Note that the input and output are provided as file paths as opposed to R objects. We specify that we want slope to be calculated in degrees as opposed to percent or radians. The function writes the result to disk, so we must read it back into R using rast() before mapping it.

wbt_slope(dem=str_glue("{outFld}dtm5.tif"),
          output=str_glue("{outFld}slp.tif"),
          units = "degrees")

slp <- rast(str_glue("{outFld}slp.tif"))

tm_shape(slp) +
  tm_raster(col.scale = tm_scale_continuous(values="brewer.yl_or_rd"))

Next, we calculate the terrain ruggedness index (TRI) after Riley et al. (1999) using wbt_ruggedness_index(). This measure characterizes local topographic heterogeneity by quantifying the elevation differences between each cell and its eight neighboring cells. Larger values indicate more rugged terrain.

wbt_ruggedness_index(dem=str_glue("{outFld}dtm5.tif"),
                     output=str_glue("{outFld}rug.tif"))

rug <- rast(str_glue("{outFld}rug.tif"))

tm_shape(rug) +
  tm_raster(col.scale = tm_scale_continuous(values="brewer.yl_or_br"))

WhiteboxTools provides many additional terrain derivatives that are useful as predictor variables in machine learning and deep learning workflows. Examples include topographic position index (TPI), surface curvatures, topographic openness, hillshades, and hydrologic variables such as flow accumulation and topographic wetness index.

16.5 Vector Geoprocessing

WhiteboxTools also provides tools for working with vector data. However, it only supports vector data in the shapefile format. So, we first read layers from a GeoPackage using sf and then save them to shapefiles using st_write(). We read in cities in th United States (point features), interstates (line features), and state boundaries (polygon features). Similar to writeRaster(), st_write() returns the data invisibly, so we can save the result to a variable and write it to disk in a single step. The delete_dsn argument allows existing files to be overwritten if the code is re-run.

us_cities <- st_read(dsn=str_glue("{fldrPth}chpt16.gpkg"),
                     "citiesUS",
                     quiet=TRUE) |>
  st_write(str_glue("{outFld}cities.shp"), delete_dsn=TRUE, quiet=TRUE)

us_inter <- st_read(dsn=str_glue("{fldrPth}chpt16.gpkg"),
                    "interstatesUS",
                    quiet=TRUE) |>
  st_write(str_glue("{outFld}inter.shp"), delete_dsn=TRUE, quiet=TRUE)

states <- st_read(dsn=str_glue("{fldrPth}chpt16.gpkg"),
                  "statesUS",
                  quiet=TRUE) |>
  st_write(str_glue("{outFld}states.shp"), delete_dsn=TRUE, quiet=TRUE)

The three layers are then mapped together: state boundaries as polygons, interstates as lines, and cities as points.

tm_shape(states) +
  tm_polygons(fill="tan")+
  tm_shape(us_inter)+
  tm_lines(lwd=1,
           col="black")+
  tm_shape(us_cities)+
  tm_bubbles(fill="red",
             size=.5)

The wbt_dissolve() function merges, or dissolves, polygons that share the same value in a specified attribute field. Here, we dissolve the state boundaries based on the “SUB_REGION” field to generate a layer of US sub-regions. The result is then read back into R and mapped.

wbt_dissolve(input=str_glue("{outFld}states.shp"),
             output=str_glue("{outFld}subregions.shp"),
             field = "SUB_REGION")

subReg <- st_read(dsn=str_glue("{outFld}subregions.shp"),
                  quiet=TRUE)

tm_shape(subReg)+
  tm_polygons(fill="SUB_REGION")+
  tm_layout(legend.outside = TRUE)

Next, we calculate the centroid of each sub-region using wbt_centroid_vector(). When we map the results, you will notice that some sub-regions have more than one centroid. This is because the dissolve operation produced a separate, single-part polygon for each spatially disconnected piece of a sub-region (for example, islands or areas separated by water bodies). A centroid is then calculated for each of these separate features.

wbt_centroid_vector(input=str_glue("{outFld}subregions.shp"),
                    output=str_glue("{outFld}subregionsC.shp"))

subRegC <- st_read(dsn=str_glue("{outFld}subregionsC.shp"),
                   quiet=TRUE)

tm_shape(subReg)+
  tm_polygons(fill="SUB_REGION")+
  tm_layout(legend.outside = TRUE)+
  tm_shape(subRegC)+
  tm_bubbles(fill="red",
             size=.5)

To obtain a single centroid per sub-region, we first convert the single-part polygons to multi-part polygons using wbt_single_part_to_multi_part(). All of the parts that share the same “SUB_REGION” value are combined into a single multi-part feature. We then re-calculate the centroids.

wbt_single_part_to_multi_part(input=str_glue("{outFld}subregions.shp"),
                              output=str_glue("{outFld}subregions2.shp"),
                              field = "SUB_REGION")

wbt_centroid_vector(input=str_glue("{outFld}subregions2.shp"),
                    output=str_glue("{outFld}subregionsC2.shp"))

subRegC2 <- st_read(dsn=str_glue("{outFld}subregionsC2.shp"),
                    quiet=TRUE)

Now, each sub-region is represented by a single centroid.

tm_shape(subReg)+
  tm_polygons(fill="SUB_REGION")+
  tm_layout(legend.outside = TRUE)+
  tm_shape(subRegC2)+
  tm_bubbles(fill="red",
             size=.5)

Since the centroid of a multi-part feature is calculated using all of its parts, it may not fall within any of the polygons that make up the feature. This is also true for single polygons with irregular, concave shapes.

16.6 Lidar Processing

Lastly, WhiteboxTools provides a variety of tools for processing lidar point clouds stored in the LAS format. Here, we use an example point cloud for an area around Madison, West Virginia, USA (“madison.las”), which is provided in the data folder for the chapter.

The wbt_lidar_digital_surface_model() function creates a digital surface model (DSM) from the point cloud. Unlike a DTM, which represents the bare earth surface, a DSM represents the elevation of the highest surface, including vegetation canopy and buildings. The DSM is generated using the highest points within the local neighborhood of each cell. We specify an output cell size of 5 m.

wbt_lidar_digital_surface_model(input=str_glue("{fldrPth}madison.las"),
                                output=str_glue("{outFld}dsm.tif"),
                                resolution=5)

dsm <- rast(str_glue("{outFld}dsm.tif"))

tm_shape(dsm) +
  tm_raster(col.scale = tm_scale_continuous(values="brewer.purples"))

The wbt_lidar_point_density() function calculates the density of lidar returns per unit area. Here, we include all returns, generate an output with a 5 m cell size, and use a 10 m search radius around each cell center to count points. Mapping point density is useful for assessing the quality and consistency of a lidar dataset. For example, low density areas may be associated with water bodies, which absorb the laser pulses, or gaps between flight lines, while high density areas may occur where flight lines overlap.

wbt_lidar_point_density(input=str_glue("{fldrPth}madison.las"),
                        output=str_glue("{outFld}pntDen.tif"),
                        returns = "all",
                        resolution = 5,
                        radius = 10)

pntDen <- rast(str_glue("{outFld}pntDen.tif"))

tm_shape(pntDen) +
  tm_raster(col.scale = tm_scale_continuous(values="brewer.reds"))

As discussed in Chapter 14, the lidR package provides native R tools for lidar processing, including generating DTMs, DSMs, and canopy height models and segmenting individual trees. WhiteboxTools is a good alternative when processing large datasets, since its tools are implemented in Rust and are generally very fast.

16.7 Concluding Remarks

In this chapter, we demonstrated how to use the whitebox package to access WhiteboxTools from R for terrain analysis, vector geoprocessing, and lidar processing. The general workflow is the same regardless of the tool being used: write the input data to disk, run the WhiteboxTools function with file paths as input and output, and read the result back into R. This same pattern applies to many of the other packages discussed in the introduction that link R to external geospatial software. Being able to access these tools greatly expands the set of analyses that can be performed within an R-based workflow.

16.8 References

Riley, S.J., DeGloria, S.D. and Elliot, R., 1999. A terrain ruggedness index that quantifies topographic heterogeneity. Intermountain Journal of Sciences, 5(1-4), pp.23-27.