Relative elevation models using lidar, R, and Python

Pretending I can make beautiful maps!

Michael Schramm (Texas Water Resources Institute)

Daniel Coe has been making wonderful relative elevation maps of meandering rivers using LIDAR data for quite some time:

60 river images in 60 seconds! Now all available in 4k resolution (3840x2160px) for non-commercial use—great for all your screen wallpaper needs. #lidar #geomorphology #geography

— Daniel Coe (@geo_coe) July 21, 2022

We even have a print hanging in our living room and I’ve been itching to try to make one myself using R. Recently, OpenTopography released the RiverREM package for Python which is kind of an easy mode button for produce a relative elevation model. I can use Python, I can use R. Now I just need to learn how to use LIDAR data…


The lidR has a host of functions for reading in lidar files or catalogs, creating digital terrain models, digital surface models, and a bunch of other stuff that is currently over my head! I primarily was interested in reading in the data, creating a high resolution digital terrain model, then getting it into RiverREM. The instruction manual linked above is a great resources and it took me about a day to become comfortable with reading in and manipulating data and making some really neat 3d maps. I’m actually excited about incorporating LIDAR into some of our hydrology and water quality projects now.


This example isn’t going to be fully reproducible because (as far as I know) downloading LIDAR data requires point and click web interfaces for the most part, and I created polygons for clipping data in ArcGIS.

library(lidR) # for lidar stuff
library(sf) # for vector spatial files
library(rayshader) # for fun 3d rendering
library(terra) # for raster data
library(reticulate) # need this to talk to python
library(viridis) # pretty colors
library(pracma) # for just one function, trust me

The general workflow is:

  1. read the LIDAR catalog (this does not load the entire LIDAR data into memory but keeps it on disk because of size)

  2. read in the polygon that has the area of interest

  3. clip the LIDAR data to the area of interest

  4. create the digital terrain model, which is the rasterized representation of the land and water surfaces in the LIDAR data.

  5. load the raster into riverREM and make a relative elevation model

  6. get the relative elevation model into rayshader and render it with desired colors and hillshading.

For reference, the LIDAR data was obtained from the TNRIS web downloader tool. I created a polygon in ArcGIS at the confluence of the Brazos and Navasota rivers for this example.

## Load/create the lidar catalog
brazos <- readLAScatalog("data/brazos/")

## this sets the temp files path and names used for files associated with this catalog
opt_output_files(brazos) <- paste0(tempdir(), "/{XCENTER}_{YCENTER}_{ID}")

## read the polygon shapefile of the area of interest
## and transform to the CRS used by the LIDAR data
conf <-st_read("data/navasota/confluence.shp") |> 
  st_transform(crs = 6343)

## clip polygon
conf_lidr <- clip_roi(brazos, conf)

## set path/filenames for conf_lidr so we don't
## overwrite brazos files
opt_output_files(conf_lidr) <- paste0(tempdir(), "/{*}_dtm")

## create the dtm, this might take a minute
terrain <- rasterize_terrain(conf_lidr, pkg = "terra")

## save the raster to a file
writeRaster(terrain, "confluence_raw_1.tif", overwrite = TRUE)

Now we have a high resolution digital terrain model raster created from the LIDAR data. This is basically a high res DEM that we can treat like any other raster data for spatial analysis. The next step is to get this into the RiverREM pacakage in Python. I’m not going into how to setup RStudio and Anaconda so you can use Python in RStudio. Plenty of tutorials and it varies slightly by OS. But RStudio has done a great job getting Python integrated into R workflows.

The RiverREM package loads the raster, looks up the hydrologic lines using the OpenStreeetMap API, points are sampled along the hydrologic feature and interpolated using inverse distance weighting across the elevation model, then the river elevation is detrended. More detailed info on the process is available on the OpenTopography blog.

## use the conda environment I have setup for this

## import riverrem
remmaker <- import("riverrem.REMMaker")

## initiate the REMmaker module
rem_maker <- remmaker$REMMaker(dem = "confluence_raw_1.tif", 
                               out_dir = "rem",
                               interp_pts = 10000, k = 100)

## create the REM

## you can create a hillshaded viz using riverREM
## if you want