3D catchment model for the Upper Ruamahanga

1 Introduction

In a previous post, I offered to create an STL file for 3D printing for any catchment groups interested in visualising their catchments. Chris Hollis, a geoscience educator and member of the Upper Ruamahanga River Management Advisory Group (URRMAC) established by Greater Wellington Regional Council, was the first to get in touch. He expressed interest in the Upper Ruamahanga, particularly in helping to understand how the active faults in the region interact with rivers.

The Ruamahanga Water Management Zone is a vast area of 3,555 km², encompassing the entire Ruamahanga River catchment and it occupies nearly half the land area of the Wellington Region. The catchment features diverse land uses, including urban townships, intensive agriculture and horticulture, native bush, and wetlands.

Chris also mentions:

A significant feature of this catchment is how the rivers flowing from the Tararua Range have been affected by movement along active faults linked to the region’s setting alongside the Hikurangi Subduction Zone.

Show the code

library(sf)
library(elevatr)
library(geodata)
library(rayshader)
library(dplyr)
library(leaflet)

2 Preparing the data

Show the code

catchment_file <-  "data/upr_ruamahanga_catchment/upr_ruamahanga_catchment.shp"
catchment_name <-  "ruamahanga"

# import catchment boundary
catchment <- st_read(catchment_file, crs = 2193, quiet = TRUE)

# get elevation data using elevatr package, from global dataset
catchment_elev <- get_elev_raster(catchment, # location of interest
                                z = 11, # zoom level (detail)
                                clip = "location", # clip to catchment boundary
                                prj = 2193) # inherit projection from catchment boundary

# convert to matrix
catchment_mat <- raster_to_matrix(catchment_elev)

# active faults database
afd_raw <- st_read("data/NZAFD/Shapefile/NZAFD_250K_Oct_2023_NZTM.shp", crs = 2193, quiet = TRUE)
afd <- st_crop(afd_raw, catchment)

# rivers
rivers_raw <- st_read("data/lds-nz-river-centrelines-topo-1500k-SHP/nz-river-centrelines-topo-1500k.shp", crs = 2193, quiet = TRUE)
rivers <- st_crop(rivers_raw, catchment)

Catchment boundary
The first step involved preparing a boundary to delineate the catchment area of interest. For the Upper Ruamahanga, this was defined as all land upstream of the confluence of the Ruamahanga and Waiohine Rivers near Greytown.

To delineate the catchment boundary, I used the Digital River Network (also known as REC). The steps were as follows:

Identify the NZSEGMENT ID of the furthest downstream river segment.
Trace the river network upstream using code to find all river segments above this point, stopping at the headwaters.
Merge and dissolve the watersheds of all contributing river segments to form a single contiguous boundary.

The code for these steps is quite involved, so I have not included it here. However, if there is enough interest, I can cover it in a future blog post. If you already have a sea-draining catchment or another pre-computed catchment boundary, you can skip this step.

Here is the catchment boundary we are working with:

Show the code

# Create leaflet map
leaflet() |>
  addProviderTiles(providers$Esri.WorldImagery, group = "World Imagery") |>
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB Positron") |>
  addPolygons(
    data = catchment |> st_transform(4326),
    color = "blue",
    weight = 2,
    opacity = 1,
    fillOpacity = 0.5,
    group = "Catchment Layer"
  ) |>
  addLayersControl(
    baseGroups = c("World Imagery", "CartoDB Positron"),
    overlayGroups = c("Catchment Layer"),
    options = layersControlOptions(collapsed = FALSE)
  )

Elevation data
Next, I used the elevatr package to download a high-resolution digital elevation model (DEM) for the catchment area. The get_elev_raster() function downloads a raster file from the Mapzen Terrain Tiles API and clips it to the catchment boundary. The raster_to_matrix() function converts the raster to a matrix that can be used by the rayshader package.

Overlays
I also imported two additional datasets to overlay on the 3D model. The first is the active faults database from GNS Science, which I cropped to the catchment boundary. The second is the NZ river centrelines dataset from the LINZ Data Service, which I similarly cropped to the catchment boundary.

3 3D rendering

The first step in rendering was to process all the data in Rayshader to create the actual mesh object. The rendered media can be utilised in various digital formats, including websites, presentations, print, and social media.

I’m quite pleased with the results. This was my first time draping vector polylines (or polygons) over the rendering, which turned out to be quite effective. It is fascinating to observe the interactions between the landforms, faults (shown in white), and river centrelines (shown in blue).

Show the code

catchment_mat |> 
  sphere_shade(texture = "imhof1", sunangle = 45) |> # creates colour gradient
  add_overlay(generate_line_overlay(afd,
                                    heightmap = catchment_mat,
                                    extent = st_bbox(catchment_elev),
                                    color = "white",
                                    linewidth = 3,
                                    offset = c(0,0))
                                    ) |> 
  add_overlay(generate_line_overlay(rivers,
                                    heightmap = catchment_mat,
                                    extent = st_bbox(catchment_elev),
                                    color = "blue",
                                    linewidth = 4,
                                    offset = c(0,0)), alphalayer = 0.7) |> 
  add_shadow(ray_shade(catchment_mat), 0.5) |> # adds shadow
  add_shadow(ambient_shade(catchment_mat), 0) |>  # adds ambient light (looks nice but very slow)
  plot_3d(catchment_mat, zscale = 28, zoom = 0.4, theta = 0, phi = 30, windowsize = c(1920, 1080)) # define plot parameters

render_highquality(paste0(catchment_name, "_static.png"),
                   parallel = TRUE,
                   samples = 128 * 2,
                   lightdirection = 75,
                   lightaltitude = 80,
                   lightintensity = 1000,
                   width = 2000,
                   height = 2000)

# Create spinning animation
render_camera(theta = seq(0, 360, length.out = 360), phi = 30, zoom = 0.5) # define camera parameters, spin around y-axis
render_movie(paste0(catchment_name, "_spinning.mp4"), fps = 30, frames = 360, zoom = 0.5) # create spinning animation

# export
save_3dprint(paste0("data/", catchment_name, ".stl"),
             maxwidth = 200,
             rotate = TRUE)

4 Export STL file for 3D printing

Show the code

# export
save_3dprint(paste0("data/", catchment_name, ".stl"),
             maxwidth = 200,
             rotate = TRUE)

The next step is to convert the data into a format suitable for 3D printing (or CNC routing, laser cutting, etc.). For this, we need an .STL file, which is accepted by slicing software used for 3D printing. The slicing software will then convert the .STL file into G-code, providing precise instructions for the 3D printer to construct the object layer by layer. The G-code details include the path the printer’s nozzle should follow, the speed at which it should move, and the amount of material to be extruded.

Much of the information included in the rendering is not needed by the slicer, such as colours, shadows, and polylines. Only the elevation data is required. This elevation data provides the necessary topographical information to the slicer, which uses it to create the appropriate toolpaths and determine the layer heights for the 3D printer.

If you want to download this model, you can find it on Thingiverse here.

4.1 3D printing at Wellington Libraries

3D printers are becoming increasingly accessible through various means. Many public libraries, universities, and community makerspaces now offer 3D printing services, providing individuals with the opportunity to use these advanced machines without needing to own one. Additionally, 3D printing service providers allow users to upload their digital models online and have the printed objects shipped to them.

Wellington Libraries are one such organisation which provide access to 3D printers, through the Tūhura HIVE makerspace at Johnsonville Library. They have four 3D printers in the HIVE, which can print up to 250mm x 210mm x 210mm in grey PLA at a cost of $0.10 per gram. You can even submit 3D printer jobs to them over email by sending them the .STL file.

5 Finished product

I’m lucky to have a 3D printer at home, so here’s a timelapse video of the model printing out on my Bambu Labs X1 Carbon where you can see the layers being extruded line by line. This took around 6.5 hours to print, but the timelapse is sped up to 6 seconds.

And here’s photos of the finished model. It measures approximately 200mm x 170mm x 15mm and is printed in green PLA.

6 Conclusion

The creation of a 3D catchment model for the Upper Ruamahanga region showcases the power of R, Rayshader, and 3D printers. The model offers a unique perspective on active faults, landforms, and river networks, serving as both an educational tool and a resource for stakeholders. This integration of geographic data and modern technology highlights the potential for enhanced visualisation and analysis in environmental management.