Skip to contents

Orchestrates the full VCA pipeline: slope thresholding, cost-distance analysis, flood surface modelling, and morphological cleanup to identify unconfined valley bottoms.

Usage

fl_valley_confine(
  dem,
  streams,
  field = "channel_width",
  slope = NULL,
  slope_threshold = 9,
  max_width = 2000,
  cost_threshold = 2500,
  flood_factor = 6,
  precip = 1,
  waterbodies = NULL,
  channel_buffer = NULL,
  size_threshold = 5000,
  hole_threshold = 2500
)

Arguments

dem

A SpatRaster of elevation.

streams

An sf linestring object or a SpatRaster of rasterized streams. If sf, it is rasterized using field.

field

Character. Column name for fl_stream_rasterize() when streams is sf. Default "channel_width".

slope

A SpatRaster of percent slope. If NULL, derived from dem.

slope_threshold

Numeric. Maximum percent slope for valley floor. Default 9.

max_width

Numeric. Maximum valley width in map units (metres). Default 2000.

cost_threshold

Numeric. Maximum accumulated cost distance. Default 2500.

flood_factor

Numeric. Multiplier on bankfull depth. Default 6.

precip

A SpatRaster or numeric scalar of precipitation. Default 1.

waterbodies

An sf polygon object of lakes and/or wetlands, or NULL (default). Waterbody polygons are rasterized onto the valley grid and added to the output after morphological cleanup. No buffer is applied — a lake or wetland in the valley is part of the flood system as-is. No spatial filtering is applied — all polygons are rasterized. Pre-filter to valley-bottom features before calling if headwater waterbodies are not wanted.

channel_buffer

Logical. Buffer streams by their channel_width attribute and add to the valley output. Default TRUE when streams is an sf object with a channel_width column, FALSE otherwise. The stream channel is floodplain but can be sub-pixel at coarse DEM resolution.

size_threshold

Numeric. Minimum valley patch area (m²). Default 5000.

hole_threshold

Numeric. Maximum hole area to fill (m²). Default 2500.

Value

A SpatRaster with binary values: 1 = unconfined valley, 0 = confined / hillslope, NA = outside analysis extent.

Details

The algorithm combines four criteria via intersection (AND):

  1. Slope mask — cells with slope <= slope_threshold

  2. Distance mask — cells within max_width / 2 of a stream

  3. Cost distance mask — cells with accumulated cost < cost_threshold

  4. Flood mask — cells identified as flooded by bankfull regression

The combined mask then undergoes morphological cleanup:

  • Closing filter (3x3) to bridge small gaps

  • Fill small holes (< hole_threshold)

  • Remove small patches (< size_threshold)

  • Majority filter (3x3) to smooth edges

After cleanup, optional features are added via logical OR:

  • Channel buffer — streams buffered by channel_width (DEM correction)

  • Waterbodies — user-supplied lake/wetland polygons rasterized as-is

Adapted from the USDA Valley Confinement Algorithm Toolbox (BlueGeo implementation by Devin Cairns, MIT license) and bcfishpass lateral habitat assembly (Simon Norris, Apache 2.0).

Performance

Several internal operations (focal filters, distance calculations, raster math) support multi-threading via terra::terraOptions(). Set threads before calling this function to speed up processing on large rasters:

terra::terraOptions(threads = 12)

On an Apple M4 Max (16 cores), 12 threads reduced runtime from ~3.5 minutes to ~1 minute for a 27M-cell raster (~2,700 km² at 10 m).

Examples

dem <- terra::rast(system.file("testdata/dem.tif", package = "flooded"))
streams <- sf::st_read(
  system.file("testdata/streams.gpkg", package = "flooded"),
  quiet = TRUE
)
precip_r <- fl_stream_rasterize(streams, dem, field = "map_upstream")

# Basic VCA (channel buffer auto-detected from streams$channel_width)
valleys <- fl_valley_confine(
  dem, streams,
  field = "upstream_area_ha",
  precip = precip_r
)
terra::plot(valleys, col = c("grey90", "darkgreen"), main = "Unconfined valleys")


# With waterbodies — fills lake/wetland donut holes
waterbodies <- sf::st_read(
  system.file("testdata/waterbodies.gpkg", package = "flooded"),
  quiet = TRUE
)
valleys_wb <- fl_valley_confine(
  dem, streams,
  field = "upstream_area_ha",
  precip = precip_r,
  waterbodies = waterbodies
)
terra::plot(valleys_wb, col = c("grey90", "darkgreen"), main = "With waterbodies")