Creating a Relative Elevation Model (REM) of the Waimakariri River

1 Introduction

Recently, I created a Relative Elevation Model (REM) for the Waimakariri River, covering the stretch from Gorge Bridge to the coast. A REM is similar to a Digital Elevation Model (DEM) but differs in that it standardises each cell’s height relative to a reference point—in this case, the river centreline. This method highlights relative variations in elevation across the riverbed that may not be as evident in a traditional DEM.

1.1 Why use a REM?

The Waimakariri River’s braided structure creates a complex topography. By using a REM, we can highlight how the riverbed varies relative to its centreline, providing insight into fluvial geomorphology, erosion, sediment deposition, and flood management. This is particularly useful for understanding and managing dynamic river systems like braided rivers. REMs can easily discern relic channels, terraces, and oxbow lakes.

1.2 Why not use a DEM?

While a Digital Elevation Model (DEM) provides the base information about the absolute elevation of the landscape, it doesn’t always capture the nuances needed for specific analyses. A DEM gives a broad overview of the topography, but it might not reveal subtle variations in elevation that are critical for understanding river geomorphology.

Because rivers always flow downhill, DEMs of rivers will always have a decreasing trend in elevation from the source to the mouth. Added to this are the smaller-scale elevation variations within the river corridor. These variations are often lost in a DEM but can be highlighted in a REM.

Comparison of DEM (left) to a REM (right) of the same area. From: https://dancoecarto.com/creating-rems-in-qgis-the-idw-method

2 Creating the REM

To create the REM, I used a combination of LiDAR data and river centreline information. The LiDAR data provided detailed elevation information, while the centreline served as the reference point for each cell.

Detrending a DEM to generate a Relative Elevation Model (REM) involves three key steps:

1. Sample Elevation Values Along the River Centreline: Begin by extracting elevation values along the river’s centreline, which serves as the reference point.

2. Interpolate River Elevations Using Inverse Distance Weighting (IDW): Instead of simply using the nearest elevation sample, IDW is applied to calculate a weighted average of multiple nearby river elevations. This technique smooths the output surface by considering the influence of several points, resulting in a more accurate representation of local elevation relative to the river.

3. Subtract River Elevations from the DEM: Finally, subtract these interpolated river elevations from the original DEM. The resulting REM assigns a value of 0 to areas at the same elevation as the nearest point on the river channel, with positive values indicating the height above the river channel.

Fortunately, the RiverREM Python package takes care of most of the heavy lifting. The package is designed specifically for creating REMs of river systems and streamlining the process for researchers and practitioners. An equivalent package isn’t available for R, so it’s create to be able to mix and match R and Python code in an analysis.

import rasterio
import rasterio.merge
import glob
import os

# Set the directory containing the GeoTIFF files
input_directory = 'data/dem'
output_directory = 'out/'
output_file = os.path.join(output_directory, 'merged_output.tif')

# Ensure the output directory exists
os.makedirs(output_directory, exist_ok=True)

# Find all GeoTIFF files in the directory
geotiff_files = glob.glob(os.path.join(input_directory, '*.tif'))

# Debugging: Print the list of found files
print("Found GeoTIFF files:", geotiff_files)

# Open all the GeoTIFF files and add them to a list
src_files_to_mosaic = []
for fp in geotiff_files:
    src = rasterio.open(fp)
    src_files_to_mosaic.append(src)

# Merge all GeoTIFFs into a single mosaic
if src_files_to_mosaic:
    mosaic, out_trans = rasterio.merge.merge(src_files_to_mosaic)
    
    # Copy the metadata from one of the input files
    out_meta = src.meta.copy()
    
    # Update the metadata to reflect the number of layers and the transform
    out_meta.update({
        "driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": out_trans,
        "count": mosaic.shape[0]  # Number of bands
    })
    
    # Write the merged raster to disk
    with rasterio.open(output_file, "w", **out_meta) as dest:
        dest.write(mosaic)
    
    print(f"Merged GeoTIFF saved to {output_file}")
else:
    print("No GeoTIFF files found in the directory. Please check the directory path and file extensions.")

# Close all the source files
for src in src_files_to_mosaic:
    src.close()

from riverrem.REMMaker import REMMaker
# provide the DEM file path and desired output directory
rem_maker = REMMaker(dem='out/merged_output.tif', out_dir='out/')
# create an REM
rem_maker.make_rem()
# create an REM visualization with the given colormap
rem_maker.make_rem_viz(cmap='mako_r')

3 Visualisations

The REM can be visualised in various ways to highlight different aspects of the riverbed. For example, a colour gradient can be used to show the relative elevation of each cell, with darker colours indicating higher elevations and cooler colours indicating lower elevations. This can help identify areas of erosion and deposition, as well as changes in the river channel over time.

Remember that in a REM, every height is by definition relative to the river centreline. And apologies for the blocky outline of the below images! This is due to the outline of the captured LiDAR imagery, which just followed the river.

3.1 Gorge Bridge

In this section, you can observe the Waimakariri River being squeezed through a narrow constriction at the Gorge Bridge and then expanding outwards as it moves downstream. This constriction-expansion sequence is characteristic of braided rivers, which transport huge quantities of sediment (especially gravel) down from high-energy mountainous regions and distribute it across floodplains further downstream.

Interestingly, this is also the finish line of the kayak leg of the Coast to Coast race, making it a notable point both geographically and recreationally!

3.2 Courtenay shallowing

Just downstream of Courtenay, the riverbed becomes noticeably shallower across its entire width, a feature clearly indicated by the relatively darker colour in the main riverbed on the REM. This section highlights the natural variability in the river’s depth. Additionally, you can observe the stopbanks on the southern side of the river, strategically placed to prevent the river from spilling over onto the Old West Coast Road.

3.3 West Melton flood defences

In this section, the extensive flood engineering works are clearly visible on both sides of the river, particularly in the form of stop-banks and groynes. These hard engineering structures are designed to control the river’s flow, manage sediment deposition, and protect surrounding areas from flooding. While highly effective in preventing immediate flood risks, they exemplify the traditional hard engineering approach, which often requires significant maintenance and can alter the natural dynamics of the river. Time will also tell how resilient these structures are in the face of climate change.

There is growing interest in complementing these hard engineering with nature-based solutions. Integrating natural floodplains, constructing wetlands, or re-establishing vegetation could enhance a river’s natural ability to manage floodwaters, offering a more sustainable and resilient approach to flood management in the long term. Combining these methods could provide the necessary protection while also supporting the river’s ecological health.

4 Interactive map

The above images were cropped, but here you can explore the Waimakariri REM along its entire length in an interactive map. Zoom in and pan around to see the detailed elevation changes and how these structures interact with the landscape.

Due to file size restrictions, the map has been downsampled to approximately 15MB, which lowers the resolution.

# Load the required libraries
library(terra)
library(leaflet)
library(leaflet.extras)

# Suppress progress bar
terraOptions(progress = 0)

# Load the raster
raster_path <- "data/merged_output_hillshade-color.tif"
r <- rast(raster_path)

desired_cells <- 200 * 1e6 / 4  # Each cell in a raster typically takes 4 bytes

# Calculate the aggregation factor to achieve the desired number of cells
agg_factor <- round(sqrt(ncell(r) / desired_cells))

# Resample the raster to reduce its size
r_resampled <- terra::aggregate(r, fact=agg_factor, method="bilinear")

# Create a leaflet map
map <- leaflet() %>% 
  addProviderTiles('Esri.WorldImagery', group = "Base Map") |> 
  addRasterImage(r_resampled, project = TRUE, opacity = 1, group = "REM", maxBytes = 30 * 1024 * 1024) |> 
  addLayersControl(
    baseGroups = c("Base Map"),
    overlayGroups = c("REM"),
    options = layersControlOptions(collapsed = FALSE)
  ) |> 
  addFullscreenControl()

# Display the map
map