drainageBasins {SAiVE}

Watershed/basin delineation

Description

Hydro-processes a DEM, creating flow accumulation, direction, and streams rasters, and (optionally) delineates watersheds above one or more points using Whitebox Tools. To facilitate this task in areas with poor quality/low resolution DEMs, can "burn-in" a stream network to the DEM to ensure proper stream placement (see details). Many time-consuming raster operations are performed, so the function will attempt to use existing rasters if they are present in the same path as the base DEM and named according to the function's naming conventions. In practice, this means that only the first run of the function needs to be very time consuming. See details for more information.

NOTE 1: This tool can be slow to execute and will use a lot of memory. Be patient, it might take several hours with a large DEM.

NOTE 2: ESRI shapefiles, on which the Whitebox Tools functions depend, truncate column names to 10 characters. You may want to save and re-assign column names to the output terra object after this function has run.

NOTE 3: If you are have already run this tool and are using a DEM in the same directory as last time, you only need to specify the DEM and the points (and, optionally, a projection for the points output). Operations using the optional streams shapefile and generating flow accumulation direction, and the artificial streams raster do not need to be repeated unless you want to use a different DEM or streams shapefile.

NOTE 4: This function is very memory (RAM) intensive. You'll want at least 16GB of RAM, and to ensure that most of it is free. If you get an error such as 'cannot allocate xxxxx bytes', you probably don't have the resources to run the tool. All rasters are un-compressed and converted to 64-bit float type before starting work, and there needs to be room to store more than twice that uncompressed raster size in memory.

Usage

drainageBasins(
  DEM,
  streams = NULL,
  breach_dist = 10000,
  threshold = 500,
  overwrite = FALSE,
  projection = NULL,
  points = NULL,
  points_name_col = NULL,
  save_path = NULL,
  snap = "nearest",
  snap_dist = 200,
  burn_dist = 10,
  n.cores = NULL,
  force_update_wbt = FALSE,
  silent_wbt = TRUE
)

Arguments

DEM

The path to a DEM including extension from which to delineate watersheds/catchments. Must be in .tif format. Derived layers such as flow accumulation, flow direction, and streams will inherit the DEM coordinate reference system.

streams

Optionally, the path to the polylines shapefile/geopackage containing lines, which can be used to improve accuracy when using poor quality DEMs. If this shapefile is the only input parameter being modified from previous runs (i.e. you've found a new/better streams shapefile but the DEM is unchanged) then specify a shapefile or geopackage lines file here and overwrite = TRUE.

breach_dist

The max radius (in raster cells) for which to search for a path to breach depressions, passed to whitebox::wbt_breach_depressions_least_cost(). This value should be high to ensure all depressions are breached. Note that the DEM is not breached in order of lowest elevation to greatest, nor is it breached sequentially (order is unknown, but the raster is presumably searched in some grid pattern for depressions). This means that flow paths may need to cross multiple depressions, especially in low relief areas.

threshold

The accumulation threshold in DEM cells necessary to start defining a stream. This streams raster is necessary to snap pout points to, so make sure not to make this number too great!

overwrite

If applicable, should rasters present in the same directory as the DEM be overwritten? This will also force the recalculation of derived layers.

projection

Optionally, a projection string in the form "epsg:3579" (find them here). The derived watersheds and point output layers will use this projection. If NULL the projection of the points will be used.

points

The path to the points shapefile (extension .shp) containing the points from which to build watersheds. The attribute of each point will be attached to the newly-created drainage polygons. Leave NULL (along with related parameters) to only process the DEM without defining watersheds.

points_name_col

The name of the column in the points shapefile containing names to assign to the watersheds. Duplicates are allowed, and are labelled with the suffix _duplicate and a number for duplicates 2 +.

save_path

The directory where you want the output shapefiles saved.

snap

Snap to the "nearest" derived (calculated) stream, or to the "greatest" flow accumulation cell within the snap distance? Beware that "greatest" will move the point downstream by up to the 'snap_dist' specified, while nearest might snap to the wrong stream.

snap_dist

The search radius within which to snap points to streams. Snapping method depends on 'snap' parameter. Note that distance units will match the projection, so probably best to work on a meter grid.

burn_dist

If specifying a streams layer polyline layer, the number of DEM units to depress the DEM along the stream trace.

n.cores

The maximum number of cores to use. Leave NULL to use all cores minus 1.

force_update_wbt

Whitebox Tools is by default only downloaded if it cannot be found on the computer, and no check are performed to ensure the local version is current. Set to TRUE if you know that there is a new version and you would like to use it.

silent_wbt

Should Whitebox tools messages be suppressed? This function prints messages to the console already but these messages can be useful if you need to do some debugging.

Details

This function uses software from the Whitebox geospatial analysis package, built by Prof. John Lindsay. Refer to this link for more information.

Creating derived raster layers without defining watersheds

This function can be run without having any specific point above which to define a watershed. This can come in handy if you need to know where the synthetic streams raster will end up to ensure that your defined watershed pour points do not end up on the wrong stream branch, or if you simply want to front-load work while you work on defining the watershed pour points. To do this, leave the parameter points and associated parameters as NULL.

Explanation of process:

Starting from a supplied DEM, the function will fill single-cell pits, burn-in a stream network depression if requested (ensuring that flow accumulations happen in the correct location), breach depressions in the digital elevation model using a least-cost algorithm (i.e. using the pathway resulting in minimal changes to the DEM considering distance and elevation) then calculate flow accumulation and direction rasters. Then, a raster of streams is created where flow accumulation is greatest. The points provided by the user are then snapped to the derived streams raster and watersheds are computed using the flow direction rasters. Finally, the watershed/drainage basin polygons are saved to the specified save path along with the provided points and the snapped pour points.

Using a streams shapefile to burn-in depressions to the DEM:

Be aware that this part of the function should ideally be used with a "simplified" streams shapefile. In particular, avoid or pre-process stream shapefiles that represent side-channels, as these will burn-in several parallel tracks to the DEM. ESRI has a tool called "simplify hydrology lines" which is great if you can ever get it to work, and WhiteboxTools has functions whitebox::wbt_remove_short_streams() to trim the streams raster, and whitebox::wbt_repair_stream_vector_topology() to help in converting a corrected streams vector to raster in the first place.

Value

A list of terra objects. If points are specified: delineated drainages, pour points as provided, snapped pour points, and the derived streams network. If no points: flow accumulation and direction rasters, and the derived streams network. If points specified, also saved to disk: an ESRI shapefile for each drainage basin, plus the associated snapped pour point and the point as provided and a shapefiles for all basins/points together. In all cases the created or discovered rasters will be in the same folder as the DEM.

Author(s)

Ghislain de Laplante ([e-mail] or [e-mail])

Examples

# Must be run with file paths as well as a save_path

# Interim raster are created in the same path as the DEM

file.copy(system.file("extdata/basin_rast.tif", package = "SAiVE"),
  paste0(tempdir(), "/basin_rast.tif"))

basins <- drainageBasins(save_path = tempdir(),
  DEM = paste0(tempdir(), "/basin_rast.tif"),
  streams = system.file("extdata/streams.gpkg", package = "SAiVE"),
  points = system.file("extdata/basin_pts.gpkg", package = "SAiVE"),
  points_name_col = "ID",
  breach_dist = 500,
  n.cores = 2)

terra::plot(basins$delineated_basins)