ch_saga
Digital terrain analysis tools for CSHShydRology based on SAGA GIS
# ch_saga
#
# set of functions for hydrological terrain analysis to comlement those
# currently available in raster/sf/sp; uses SAGA GIS via RSAGA as a
# computational engine
#
# 2019-Sept-04 RDM
#
# updated 2019-Sept-24 RDM
# - moved ch_clear_wd() to end to avoid error when trying to access
# raster data stored on disk
###################################################################################
rm(list = ls())
library(raster)## Loading required package: splibrary(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(ggplot2)Overview
The raster and sf packages provide a good range of tools for analysing and visualizing spatial data for applications in hydrology and watershed science. However, a handful of critical tools are currently lacking, including the following:
- removal of sinks prior to flowpath modelling and catchment delineation
- flowpath modelling and determination of contributing area grids
- catchment delineation
- determination of drainage networks
At present, to accomplish these tasks, hydrologists must resort to using GIS tools such as SAGA GIS, QGIS, GRASS and ArcInfo. However, use of these tools constrains the development of integrated, reproducible workflows.
The objective of the spatial_hydrology theme within CSHShydrology is to develop functions that provide these capabilities and thus facilitate the development of integrated reproducible workflows within the R universe.
While an ultimate aim is to develop functions to be compiled as native R functions, the interim approach here is to develop tools based on RSAGA, which is an interface that allows SAGA GIS commands to be executed from within an R script. One feature of SAGA is that it stores function inputs and outputs in files within a directory/folder on the user’s hard drive. Therefore, one needs to read the inputs from SAGA formats into R objects in order to use them within R scripts.
The functions described below are essentially wrappers around functions in the RSAGA package. In particular, they handle the reading and writing of data between R objects and SAGA files in a manner that is mostly unseen by the user.
Some of the SAGA functions support many options. In these functions, we have opted to support only those functions that should be most important for routine hydrological analyses, particularly in relation to catchment delineation. Users that desire access to the full range of options within SAGA would need to work directly with RSAGA functions.
One disadvantage of using RSAGA is that we have heard reports that some users have experienced difficulty in running RSAGA on Mac OS. Therefore, a next step will be to develop a parallel set of functions based on GRASS GIS using the RGRASS package. GRASS is well supported for Windows, Linux/UNIX and Mac OS.
Function definitions
ch_create_wd
This function simply creates a working directory for reading and writing SAGA files. The directory name is specified by the user. It will output a warning if the directory already exists.
#######################################################################################
# Function definitions
#######################################################################################
## ch_create_wd()
ch_create_wd <- function(wd) {
# creates working directory for saga files
if (dir.exists(wd)) {
# print warning if directory already exists
print(paste("Warning: a directory named", wd, "exists"))
} else {
# create directory
dir.create(wd)
}
} # end function ch_create_wdch_clear_wd
Should the user not wish to retain the directory containing the SAGA files, this function will empty the directory and remove it.
Note: the data for raster layers read in as sdat files (e.g., carea) are held on disk rather than in memory; if you clear the working directory you will get an error if you try to use the raster layer.
## ch_clear_wd
ch_clear_wd <- function(wd) {
# empty and remove saga working directory
# note: the data for raster layers read in as sdat files
# (e.g., carea) are held on disk rather than in memory; if you
# clear the working directory you will get an error if you try
# to use the raster layer
filelist = list.files(wd)
file.remove(paste0(wd, "/", filelist))
unlink(wd, recursive = TRUE)
} # end function ch_clear_wdch_saga_fillsinks
This function removes sinks by filling them. To ensure continuous flow through the DEM, a minimum slope can be imposed on cells. Although there are other options, the most imporant are method and minslope; these are implemented in the function. For a full list of available options type ?rsaga.fill.sinks at the console.
The main input is a DEM as a raster object. The main output will be a sink-filled DEM with the same topology and coordinate reference system as the input raster.
## ch_saga_fillsinks
ch_saga_fillsinks <- function(dem_raw, saga_wd,
sinkmethod = "planchon.darboux.2001",
minslope = 0.1,
saga.env = RSAGA::rsaga.env()) {
#
# removes sinks in a dem by filling; required prior to modelling
# flow accumulation and catchment delineation; see also
# ch_saga_removesinks
#
# arguments
# ---------
# dem_raw Raster object of your raw dem in the desired projection
# saga_wd working directory to write and read saga files
# sinkmethod character. other options are "xxl.wang.liu.2006" and "wang.liu.2006"
# minslope minimum slope to impose on filled cells to force drainage
# saga.env Saga environment object. Default is to let saga find it on its own.
#
# outputs
# -------
# dem_ns - filled dem as a raster object
library(RSAGA)
library(raster)
# error trap - saga_wd does not exist
if (!dir.exists(saga_wd)) {
print("saga_wd does not exist")
return(NA)
}
# store the input dem in a file in the working directory
raster::writeRaster(dem_raw, paste0(saga_wd, "/dem_raw.sdat"), format = "SAGA",
NAflag = -9999, overwrite = TRUE)
# fill sinks
RSAGA::rsaga.fill.sinks(in.dem = paste0(saga_wd, "/dem_raw.sgrd"),
out.dem = paste0(saga_wd, '/dem_ns.sgrd'),
method = sinkmethod, minslope = minslope,
env = saga.env)
# create filled dem as a raster object
ns_file <- paste0(saga_wd, "/dem_ns.sdat")
dem_ns <- raster::raster(ns_file, format = "SAGA")
crs(dem_ns) <- crs(dem_raw)
return(dem_ns)
} # end function ch_saga_fillsinksch_saga_removesinks
This function removes sinks by deepening the drainage network.
The main input is a DEM as a raster object. The main output will be a sinkless DEM with the same topology and coordinate reference system as the input raster.
The function is still in progress, pending research to determine the options that are available and implementing them (if deemed appropriate).
## ch_saga_removesinks
ch_saga_removesinks <- function(dem_raw, saga_wd,
saga.env = RSAGA::rsaga.env()) {
#
# removes sinks by deepening drainage network; alternative to
# ch_saga_fillsinks
#
# arguments
# ---------
# raw_dem Raster object of your raw dem in the desired projection
# saga_wd working directory to write and read saga files
# saga.env Saga environment object. Default is to let saga find it on its own.
#
# outputs
# -------
# dem_ns processed dem as a raster object
#
# TODO
# - find out what options are available and implement them
library(RSAGA)
library(raster)
# error trap - saga_wd does not exist
if (!dir.exists(saga_wd)) {
print("saga_wd does not exist")
return(NA)
}
# store the input dem in a file in the working directory
raster::writeRaster(dem_raw, paste0(saga_wd, "/dem_raw.sdat"), format = "SAGA",
NAflag = -9999, overwrite = TRUE)
# remove sinks
RSAGA::rsaga.sink.removal(in.dem = paste0(saga_wd, "/dem_raw.sgrd"),
out.dem = paste0(saga_wd, '/dem_ns.sgrd'),
method = "deepen drainage route",
env = saga.env)
# create filled dem as a raster object
ns_file <- paste0(saga_wd, "/dem_ns.sdat")
dem_ns <- raster::raster(ns_file, format = "SAGA")
crs(dem_ns) <- crs(dem_raw)
return(dem_ns)
} # end function ch_saga_removesinksch_saga_carea
This function generates a raster of contributing areas based on a DEM, which should have had sinks removed in a pre-processing step.
For more information, see ?rsaga.topdown.processing.
## ch_saga_carea
ch_saga_carea <- function(dem, saga_wd,
method = 4,
linear_threshold = Inf,
saga.env = RSAGA::rsaga.env()) {
#
# generates a grid of contributing area for each dem grid cell
#
# arguments
# ---------
# dem
# - Raster object of your dem in the desired projection - should
# have had sinks removed
# saga.env
# - Saga environment object. Default is to let saga find it on its own.
#
# outputs
# -------
# carea raster of contributing areas for each pixel
#
# by-product: carea.sdat in saga_wd
#
library(RSAGA)
library(raster)
library(sf)
# error trap - saga_wd does not exist
if (!dir.exists(saga_wd)) {
print("saga_wd does not exist")
return(NA)
}
# store the dem object in the working directory
raster::writeRaster(dem, paste0(saga_wd, "/dem.sdat"), format = "SAGA",
NAflag = -9999, overwrite = TRUE)
# generate contributing area grid
RSAGA::rsaga.topdown.processing(paste0(saga_wd, "/dem.sgrd"),
out.carea = paste0(saga_wd, "/carea.sgrd"),
method = method,
linear.threshold = linear_threshold,
env = saga.env)
# read the catchment area grid into a raster object and return
carea <- raster::raster(paste0(saga_wd, '/carea.sdat'))
crs(carea) <- crs(dem)
return(carea)
} # end function ch_saga_careach_saga_catchment
This function generates catchment boundaries for a set of specified points of interest along a stream network. The output is a sf object containing the polygons.
The function has two gridded inputs: a DEM (with sinks removed as appropriate) and a grid of contributing areas. Other inputs include outlet, a sf object containing the coordinates of the catchment outlets, and outlet_label, an optional vector of character strings that serve as labels or descriptions of the outlet points.
The DEM is input as a raster object via the dem argument.
There are three ways to provide the contributing area grid, which can be accomplished using the carea and carea_flag arguments, as summarized below.
| carea | carea_flag | Explanation |
|---|---|---|
| NULL | 0 | compute from DEM |
| NULL | 1 | read in existing carea.sdat |
| raster_name | NULL | input as raster object |
The first option (carea = NULL and carea_flag = 0) would be appropriate if you have not generated a contributing area grid prior to application of ch_saga_catchment. If you use this option, the contributing area grid is generated using the standard SAGA defaults, with no options for the user to change them.
The second option (carea = NULL and carea_flag = 1) would be appropriate if you have run the ch_saga_carea function and have not cleared the working directory.
The third option would be appropriate if you have generated a contributing area grid as a raster object using either ch_saga_carea or another function.
The second and third options are particularly useful if the user either (a) wants to use a different approach than the SAGA defaults to generate the contributing area grid or (b) wants to use the contributing area grid in other analyses (e.g., computing a topographic wetness index).
Note that only one of the carea and carea_flag arguments need be specified.
## ch_saga_catchment
ch_saga_catchment <- function(dem, saga_wd, outlet,
carea = NULL, carea_flag = 0,
outlet_label = NULL,
buffsize = 100,
saga.env = RSAGA::rsaga.env()) {
#
# determines catchment polygons for a set of points of interest
#
# arguments
# ---------
# dem
# - Raster object of your dem in the desired projection - should
# have had sinks removed
# carea
# - raster object containing contributing areas (default none provided)
# carea_flag
# - if carea = NULL, 0 = create carea from dem; 1 = read in carea.sdat
# sagawd
# - name of working directory
# outlet
# - sf point object with coordinates of catchment outlets in
# the same projection as the dem
# outlet_label
# - character vector of labels; if "NULL", numbers are assigned
# buffsize
# - numeric; the buffer radius (m) around catchment outlet to
# find location on digital stream network.
# saga.env
# - Saga environment object. Default is to let saga find it on its own.
#
# outputs
# -------
# catchment_sf sf object containing catchment polygons
#
library(RSAGA)
library(raster)
library(sf)
# error trap - saga_wd does not exist
if (!dir.exists(saga_wd)) {
print("saga_wd does not exist")
return(NA)
}
# store the dem in the working directory
raster::writeRaster(dem, paste0(saga_wd, "/dem.sdat"), format = "SAGA",
NAflag = -9999, overwrite = TRUE)
# if a contributing area raster was not included, create one
if (is.null(carea)) {
if (carea_flag == 0) {
# generate contributing area grid - use all defaults
RSAGA::rsaga.topdown.processing(in.dem = paste0(saga_wd, "/dem.sgrd"),
out.carea = paste0(saga_wd, "/carea.sgrd"),
env = saga.env)
}
# read the contributing area grid into a raster object
carea <- raster::raster(paste0(saga_wd, '/carea.sdat'))
}
# extract coordinates into a matrix
xy = st_coordinates(outlet)
nc = nrow(outlets)
polygon_sfc <- sf::st_sfc()
# print catchment area with units
if (is.null(outlet_label)) {
labels <- as.character(1:nc)
} else {
labels <- outlet_label
}
for (i in 1:nc) {
# check that outlet lies within dem bbox
bb <- bbox(dem)
if (xy[i, 1] < bb[1, 1] | xy[i, 1] > bb[1,2] |
xy[i, 2] < bb[2, 1] | xy[i, 2] > bb[2,2]) {
print(cat("outlet", i, "does not lie within dem"))
polygon_sfc[i] <- NULL
} else {
# extract a window around around the outlet
buffer <- raster::extract(carea, outlets[i,], buffer = buffsize,
cellnumbers = TRUE)[[1]] %>%
as.data.frame
# find the location of the maximum catchment area on the grid, given as the id from the raster
snap_id <- buffer$cell[which.max(buffer$value)]
# get the xy coordinates at that max location, which is now going to be the location of the outlet
snap_xy <- raster::xyFromCell(carea, snap_id)
# make catchment as grid
RSAGA::rsaga.geoprocessor(lib = 'ta_hydrology', 4,
param = list(TARGET_PT_X = snap_xy[1,1],
TARGET_PT_Y = snap_xy[1,2],
ELEVATION = paste0(saga_wd, '/dem.sgrd'),
AREA = paste0(saga_wd, '/bounds.sgrd'),
METHOD = 0),
env = saga.env)
# convert grid to vector shapefile
RSAGA::rsaga.geoprocessor(lib = 'shapes_grid', 6,
param = list(GRID = paste0(saga_wd, '/bounds.sgrd'),
POLYGONS = paste0(saga_wd, '/bounds.shp'),
CLASS_ALL = 0,
CLASS_ID = 100,
SPLIT = 0),
env = saga.env)
# read polygon shapefile and save in sf object
polygon_sfc[i] <- st_geometry(sf::st_read(paste0(saga_wd, '/bounds.shp')))
}
}
# create data frame
catchment_sf <- data.frame(label = outlet_label,
outlet_x = xy[, 1], outlet_y = xy[, 2]) %>%
cbind(polygon_sfc) %>%
st_as_sf()
# set crs
sf::st_crs(catchment_sf) <- raster::crs(dem)@projargs
return(catchment_sf)
} # end function ch_saga_catchmentch_saga_channels
This function generates a drainage network using the SAGA ta_channels module with option 0.
The function requires a DEM and a contributing area grid. Like ch_saga_catchment, the user can provide the contributing area grid either as a raster object or by reading in a SAGA file; alternatively, the function can compute the contributing area grid from the DEM.
The SAGA function generates three outputs:
- a shape file containing the drainage network,
- a grid containing flow directions, and
- a grid in which, if a cell is part of a channel, its value equals the channel order; otherwise the cell is marked as no-data.
The output is a list, which will have elements that will include the channel network shape as an sf object (if out_shp = TRUE), the flow direction grid as a raster object (if out_route = TRUE) and the network grid as a raster object (if out_ntwrk = TRUE)
For more details, see the following:
## ch_saga_channels
ch_saga_channels <- function(dem, saga_wd, carea = NULL, carea_flag = 0,
out_shp = TRUE,
out_ntwrk = FALSE, out_route = FALSE,
initmethod = 2, initvalue = 0,
divcells = 5, minlen = 10,
saga.env = RSAGA::rsaga.env()) {
#
# generates a channel network as a sf vector object
#
# arguments
# ---------
# dem
# - Raster object of your dem in the desired projection - must
# have had sinks removed
# saga_wd
# - working directory
# carea
# - contributing area grid as a raster object; default is none provided
# carea_flag
# - if carea not provided (= NULL)
# 0 => create one from dem
# 1 => read one named 'carea.sdat' from saga_wd
# out_shp
# - if TRUE, return channel network as sf object
# out_ntwrk, out_route
# - if TRUE, return specified object as a raster
# initvalue, initmethod, divcells, minlen
# - parameters for SAGA function
# saga.env
# - Saga environment object. Default is to let saga find it on its own.
#
# outputs
# -------
# out_list
# - list containing (as specified) channels (shp), ntwrk (raster)
# and route (raster)
library(RSAGA)
library(raster)
library(sf)
# error trap - saga_wd does not exist
if (!dir.exists(saga_wd)) {
print("saga_wd does not exist")
return(NA)
}
# make a temporary directory within the working directory and store the dem there
raster::writeRaster(dem, paste0(saga_wd, "/dem.sdat"), format = "SAGA",
NAflag = -9999, overwrite = TRUE)
# if a contributing area raster was not included, create one
if (is.null(carea)) {
if (carea_flag == 0) {
# generate contributing area grid - use all defaults
RSAGA::rsaga.topdown.processing(in.dem = paste0(saga_wd, "/dem.sgrd"),
out.carea = paste0(saga_wd, "/carea.sgrd"),
env = saga.env)
}
carea <- raster::raster(paste0(saga_wd, '/carea.sdat'))
}
# generate channel network
RSAGA::rsaga.geoprocessor(lib = 'ta_channels', 0,
param = list(ELEVATION = paste0(saga_wd, '/dem_ns.sgrd'),
INIT_GRID = paste0(saga_wd, '/carea.sgrd'),
INIT_METHOD = initmethod,
INIT_VALUE = initvalue,
DIV_CELLS = divcells,
MINLEN = minlen,
CHNLNTWRK = paste0(saga_wd, '/ntwrk.sdat'),
CHNLROUTE = paste0(saga_wd, '/route.sdat'),
SHAPES = paste0(saga_wd, '/channels.shp')),
env = saga.env)
out_list <- list(channels = NULL, ntwrk = NULL, route = NULL)
if (out_shp) {
channels_sf <- st_geometry(sf::st_read(paste0(saga_wd, '/channels.shp')))
sf::st_crs(channels_sf) <- raster::crs(dem)@projargs
out_list$channels <- channels_sf
}
if (out_ntwrk) {
ntwrk <- raster::raster(paste0(saga_wd, '/ntwrk.sdat'))
crs(ntwrk) <- crs(dem)
out_list$ntwrk <- ntwrk
}
if (out_route) {
route <- raster::raster(paste0(saga_wd, '/route.sdat'))
crs(route) <- crs(dem)
out_list$route <- route
}
return(out_list)
} # end ch_channelsch_contours
This function generates contour lines from a DEM, which are returned as an sf object. The user can either provide a vector of elevation values by specifying the z_levels argument, or by supplying the minimum and maximum elevations (zmin and zmax) and the number of contour lines (n_levels).
## ch_contours
ch_contours <- function(dem,
zmin = NULL, zmax = NULL,
n_levels = 10,
z_levels = NULL) {
# arguments
# ---------
# dem
# - Raster object of your dem in the desired projection
# zmin, zmax, n_levels
# - min and max values for contours and number of contours
# z_levels
# - levels at which to plot contours
#
# outputs
# -------
# contours_sf sf object containing contours
#
library(raster)
library(sf)
# determine contour levels
if (is.null(z_levels)) {
z = raster::getValues(dem)
if (is.null(zmin)) zmin <- min(z, na.rm = TRUE)
if (is.null(zmax)) zmax <- max(z, na.rm = TRUE)
z_levels <- seq(zmin, zmax, length.out = n_levels)
}
# if dem includes sea level, start contours at 0.1 m to mimic coastline
if (z_levels[1] <= 0) z_levels[1] <- 0.1
# generate contours as a sf object
contours_sf <- raster::rasterToContour(dem, levels = z_levels) %>%
sf::st_as_sf()
st_crs(contours_sf) <- crs(dem)
return(contours_sf)
} # end function ch_contours()ch_checkcatchment
This function generates a simple map to allow a visual assessment of the catchment boundaries relative to the elevation contours. It also generates a table summarizing the catchments, including the coordinates of the outlet point and the catchment area.
## ch_checkcatchment
ch_checkcatchment <- function(dem, catchment, outlet, outlet_label = NULL) {
# arguments
# ---------
# dem raster. dem catchment was generated from
# catchment sf. catchment polygon
# outlet sf. location of catchment outlet
#
# outputs
# -------
# generates contour map with catchment polygon and outlet
#
library(sf)
library(raster)
library(ggplot2)
library(ggspatial)
# create contours and get bounding box to set map limits
contours = ch_contours(dem)
bb = st_bbox(contours)
# generate map
check_map <- ggplot2::ggplot(data = contours) +
geom_sf(data = contours, color = "grey") +
geom_sf(data = outlet, pch = 21, bg = "blue") +
geom_sf(data = sf::st_geometry(catchment), fill = NA, color = "red") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
print(check_map)
nc <- length(catchment)
# print catchment area with units
if (is.null(outlet_label)) {
labels <- as.character(1:nc)
} else {
labels <- outlet_label
}
area <- sf::st_area(catchment)
units <- rep(paste0(attr(area, "units")$numerator[1], "^2"), length(area))
value <- round(as.numeric(area))
area_df <- outlet %>%
mutate(label = labels, area = value, units = units)
print(area_df)
} # end function ch_checkcatchmentch_checkchannels
This function generates a simple map of the drainage network plotted over the contours to allow a visual assessment.
## ch_checkchannels
ch_checkchannels <- function(dem, channels, outlet) {
# arguments
# ---------
# dem raster. dem catchment was generated from
# catchment sf. catchment polygon
# outlet sf. location of catchment outlet
#
# outputs
# -------
# generates contour map with catchment polygon and outlet
#
library(sf)
library(raster)
library(ggplot2)
library(ggspatial)
contours = ch_contours(dem)
# get bounding box for contours to set map limits
bb = st_bbox(contours)
# generate map
check_map <- ggplot2::ggplot(data = contours) +
geom_sf(data = contours, color = "grey") +
geom_sf(data = outlet, pch = 21, bg = "black") +
geom_sf(data = sf::st_geometry(channels), color = "blue") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
print(check_map)
} # end function ch_checkchannelsApplication - Shannon Falls Creek, BC
To illustrate the functions, they will be applied to a DEM containing the area near the head of Howe Sound, including the drainage area for Shannon Falls Creek. If you want to run this script, you can access the DEM via the following link:
Note that the name of the DEM in the zipped archive is different than than the name of the DEM on my hard drive, so that will have to be changed in the script.
Read in DEM (in longlat) and reproject to UTM
#######################################################################################
# application - Shannon Falls Creek, BC
#######################################################################################
library(magrittr)##
## Attaching package: 'magrittr'## The following object is masked from 'package:raster':
##
## extractlibrary(raster)
library(sf)
library(rgdal)## rgdal: version: 1.4-4, (SVN revision 833)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/dan moore/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/dan moore/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggspatial)
setwd("c:/Extension/R/CSHShydRology/catchment")
# read in digital elevation model and project to UTM
dem <- raster::raster("c:/Research/Shannon Falls/Map/cdem_dem_151108_134502_tif/cdem_dem_151108_134502.tif")
dem.utm <- raster::projectRaster(dem, crs = "+proj=utm +zone=10 +ellps=WGS84", res = 25, method = "bilinear") Define outlets and specify arguments
# catchment outlet in long, lat, converted to sf object, then projected to same crs as dem.utm
long <- c(-123.16139, -123.151173, -123.152317, -123.16)
lat <- c(49.67105, 49.675904, 49.666367, 50.1)
label <- c("Shannon Creek 1", "Chief Creek", "Shannon Creek 2", "not in dem")
outlet_ll <- data.frame(label, long, lat) %>%
mutate(label = label) %>%
sf::st_as_sf(coords = c("long", "lat"))
sf::st_crs(outlet_ll) <- 4326
# reproject to utm
(outlets <- sf::st_transform(outlet_ll, crs = raster::crs(dem.utm)@projargs))## Simple feature collection with 4 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 488354.9 ymin: 5501548 xmax: 489093.2 ymax: 5549762
## epsg (SRID): 32610
## proj4string: +proj=utm +zone=10 +ellps=WGS84 +units=m +no_defs
## label geometry
## 1 Shannon Creek 1 POINT (488354.9 5502070)
## 2 Chief Creek POINT (489093.2 5502608)
## 3 Shannon Creek 2 POINT (489008.5 5501548)
## 4 not in dem POINT (488557.1 5549762)# specify arguments
sinkmethod <- "xxl.wang.liu.2006"
minslope <- 0.1
buffsize <- 50
saga_wd <- paste0(getwd(), "/scratch")Apply functions
1. Filling sinks and running ch_carea
The workflow below removes sinks in the DEM by filling them. It then generates a contributing area raster using the carea() function. When the ch_catchments function is run, the contributing area grid is read in from the working directory.
# apply functions - first by filling sinks and creating carea
ch_create_wd(saga_wd)## [1] "Warning: a directory named C:/Extension/R/CSHShydRology/ch_saga/scratch exists"dem_filled <- ch_saga_fillsinks(dem_raw = dem.utm, saga_wd = saga_wd,
sinkmethod = sinkmethod,
minslope = minslope)## Loading required package: gstat## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Loading required package: shapefiles## Loading required package: foreign##
## Attaching package: 'shapefiles'## The following objects are masked from 'package:foreign':
##
## read.dbf, write.dbf## Loading required package: plyr## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_preprocessor
## library : ta_preprocessor
## tool : Fill Sinks XXL (Wang & Liu)
## identifier : 5
## author : Copyrights (c) 2007 by Volker Wichmann
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_raw.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## DEM: cdem_dem_151108_134502
## Filled DEM: Filled DEM
## Minimum Slope [Degree]: 0.100000
##
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_ns.sgrd...carea <- ch_saga_carea(dem_filled, saga_wd)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Flow Accumulation (Top-Down)
## identifier : 0
## author : O.Conrad (c) 2001-2016, Portions by T.Grabs (c) 2010
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Elevation: dem_ns
## Sink Routes: <not set>
## Weights: <not set>
## Flow Accumulation: Flow Accumulation
## Input for Mean over Catchment: <not set>
## Material for Accumulation: <not set>
## Step: 1
## Flow Accumulation Unit: cell area
## Flow Path Length: <not set>
## Channel Direction: <not set>
## Method: Multiple Flow Direction
## Thresholded Linear Flow: no
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/carea.sgrd...catchments <- ch_saga_catchment(dem_filled, saga_wd = saga_wd,
outlet = outlets,
carea_flag = 1,
outlet_label = label, buffsize = buffsize)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 488344.659487
## Target Y coordinate: 5502095.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 488307.2 ymin: 5497983 xmax: 494607.2 ymax: 5502458
## epsg (SRID): NA
## proj4string: NA
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 489119.659487
## Target Y coordinate: 5502620.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 489107.2 ymin: 5501508 xmax: 492832.2 ymax: 5504908
## epsg (SRID): NA
## proj4string: NA
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 488994.659487
## Target Y coordinate: 5501570.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 488657.2 ymin: 5497983 xmax: 494607.2 ymax: 5502458
## epsg (SRID): NA
## proj4string: NA
## outlet 4 does not lie within demNULLch_checkcatchment(catchment = catchments, dem = dem.utm, outlet = outlets,
outlet_label = label)
## Simple feature collection with 4 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 488354.9 ymin: 5501548 xmax: 489093.2 ymax: 5549762
## epsg (SRID): 32610
## proj4string: +proj=utm +zone=10 +ellps=WGS84 +units=m +no_defs
## label geometry area units
## 1 Shannon Creek 1 POINT (488354.9 5502070) 13620625 m^2
## 2 Chief Creek POINT (489093.2 5502608) 6596875 m^2
## 3 Shannon Creek 2 POINT (489008.5 5501548) 13096875 m^2
## 4 not in dem POINT (488557.1 5549762) 0 m^2channel_list <- ch_saga_channels(dem_filled, saga_wd, carea_flag = 1,
initmethod = 2, initvalue = 1e5,
divcells = 5, minlen = 10)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_channels
## library : ta_channels
## tool : Channel Network
## identifier : 0
## author : (c) 2001 by O.Conrad
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_ns.sgrd...
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/carea.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Elevation: dem_ns
## Flow Direction: <not set>
## Channel Network: Channel Network
## Channel Direction: Channel Direction
## Channel Network: Channel Network
## Initiation Grid: carea
## Initiation Type: Greater than
## Initiation Threshold: 100000.000000
## Divergence: <not set>
## Tracing: Max. Divergence: 5
## Tracing: Weight: <not set>
## Min. Segment Length: 10
##
## Create index: dem_ns
## Channel Network: Pass 1
## Channel Network: Pass 2
## Channel Network: Pass 3
## Channel Network: Pass 4
## Channel Network: Pass 5
## Channel Network: Pass 6
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/ntwrk.sdat...
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/route.sdat...
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/channels.shp...
## Reading layer `channels' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\channels.shp' using driver `ESRI Shapefile'
## Simple feature collection with 993 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 482669.7 ymin: 5496396 xmax: 495669.7 ymax: 5516421
## epsg (SRID): NA
## proj4string: NAchannels = channel_list$channels
ch_checkchannels(dem.utm, channels, outlet = outlets)
contours <- ch_contours(dem.utm)
bb <- st_bbox(contours)The code chunk below generates maps by combining several of the layers.
# map with several layers, including log10(contributing area),
# elevation contours, outlets and catchment boundaries
ggplot() +
layer_spatial(log10(carea)) +
scale_fill_continuous(expression(log[10](italic(CA))),
low = "lightgrey", high = "black") +
geom_sf(data = contours, color = "yellow") +
geom_sf(data = outlets, pch = 21, bg = "red") +
geom_sf(data = sf::st_geometry(catchments), fill = NA, color = "red") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
# map with elevation contours, channels, outlets and catchment boundaries
ggplot() +
geom_sf(data = contours, color = "grey") +
geom_sf(data = outlet_ll, pch = 21, bg = "red") +
geom_sf(data = sf::st_geometry(catchments), fill = NA, color = "red") +
geom_sf(data = sf::st_geometry(channels), color = "blue") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
ch_clear_wd(saga_wd)2. Remove sinks by deepening the drainage network and have ch_saga_catchment generate carea from the DEM
In the script below, ch_carea is not run. Therefore, ch_saga_catchment is run with carea_flag = 0, so that the contributing area grid is generated within the function and is stored as a file named carea.sdat in the working directory. When ch_saga_channels is run, carea_flag is set to 1, so that the contributing area grid is read in from carea.sdat.
It is interesting to see the contrast between the channel network generated by deepening the drainage network versus that generated by filling sinks.
# now remove sinks by deepening the drainage network and
# have ch_saga_catchment generate carea from dem
ch_create_wd(saga_wd)
dem_ns <- ch_saga_removesinks(dem = dem.utm, saga_wd = saga_wd)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_preprocessor
## library : ta_preprocessor
## tool : Sink Removal
## identifier : 2
## author : O. Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_raw.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## DEM: cdem_dem_151108_134502
## Sink Route: <not set>
## Preprocessed DEM: Preprocessed DEM
## Method: Fill Sinks
## Threshold: no
## Threshold Height: 100.000000
##
## Create index: cdem_dem_151108_134502 [no sinks]
## Find Pits
## Find Outlets
## Routing
## Finalize
## number of processed sinks: 1938
## Initializing direction matrix...
## I'm fillin'...
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_ns.sgrd...catchments <- ch_saga_catchment(dem_ns, saga_wd = saga_wd,
outlet = outlets,
carea_flag = 0,
outlet_label = label, buffsize = buffsize)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Flow Accumulation (Top-Down)
## identifier : 0
## author : O.Conrad (c) 2001-2016, Portions by T.Grabs (c) 2010
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Elevation: dem_ns
## Sink Routes: <not set>
## Weights: <not set>
## Flow Accumulation: Flow Accumulation
## Input for Mean over Catchment: <not set>
## Material for Accumulation: <not set>
## Step: 1
## Flow Accumulation Unit: cell area
## Flow Path Length: <not set>
## Channel Direction: <not set>
## Method: Multiple Flow Direction
## Thresholded Linear Flow: no
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/carea.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 488344.659487
## Target Y coordinate: 5502095.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 488307.2 ymin: 5497983 xmax: 494607.2 ymax: 5502458
## epsg (SRID): NA
## proj4string: NA
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 489119.659487
## Target Y coordinate: 5502620.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 489107.2 ymin: 5501508 xmax: 492832.2 ymax: 5504908
## epsg (SRID): NA
## proj4string: NA
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Upslope Area
## identifier : 4
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Target Area: <not set>
## Target X coordinate: 488994.659487
## Target Y coordinate: 5501570.907952
## Elevation: dem_ns
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Deterministic 8
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: shapes_grid
## library : shapes_grid
## tool : Vectorising Grid Classes
## identifier : 6
## author : O.Conrad (c) 2008
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Grid: bounds
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## edge detection
## edge collection
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/bounds.shp...
## Reading layer `bounds' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\bounds.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 488657.2 ymin: 5497983 xmax: 494607.2 ymax: 5502458
## epsg (SRID): NA
## proj4string: NA
## outlet 4 does not lie within demNULLch_checkcatchment(catchment = catchments, dem = dem.utm, outlet = outlets,
outlet_label = label)
## Simple feature collection with 4 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 488354.9 ymin: 5501548 xmax: 489093.2 ymax: 5549762
## epsg (SRID): 32610
## proj4string: +proj=utm +zone=10 +ellps=WGS84 +units=m +no_defs
## label geometry area units
## 1 Shannon Creek 1 POINT (488354.9 5502070) 13617500 m^2
## 2 Chief Creek POINT (489093.2 5502608) 6543125 m^2
## 3 Shannon Creek 2 POINT (489008.5 5501548) 13093750 m^2
## 4 not in dem POINT (488557.1 5549762) 0 m^2channel_list <- ch_saga_channels(dem_ns, saga_wd, carea_flag = 1,
initmethod = 2, initvalue = 5e4,
divcells = 5, minlen = 10)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_channels
## library : ta_channels
## tool : Channel Network
## identifier : 0
## author : (c) 2001 by O.Conrad
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem_ns.sgrd...
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/carea.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Elevation: dem_ns
## Flow Direction: <not set>
## Channel Network: Channel Network
## Channel Direction: Channel Direction
## Channel Network: Channel Network
## Initiation Grid: carea
## Initiation Type: Greater than
## Initiation Threshold: 50000.000000
## Divergence: <not set>
## Tracing: Max. Divergence: 5
## Tracing: Weight: <not set>
## Min. Segment Length: 10
##
## Create index: dem_ns
## Channel Network: Pass 1
## Channel Network: Pass 2
## Channel Network: Pass 3
## Channel Network: Pass 4
## Channel Network: Pass 5
## Channel Network: Pass 6
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/ntwrk.sdat...
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/route.sdat...
## Save shapes: C:/Extension/R/CSHShydRology/ch_saga/scratch/channels.shp...
## Reading layer `channels' from data source `C:\Extension\R\CSHShydRology\ch_saga\scratch\channels.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1184 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 482669.7 ymin: 5496396 xmax: 495669.7 ymax: 5516421
## epsg (SRID): NA
## proj4string: NAchannels <- channel_list$channels
ch_checkchannels(dem_ns, channels, outlet = outlets)
contours <- ch_contours(dem.utm)
bb <- st_bbox(contours)# map with several layers, including log10(contributing area),
# elevation contours, outlets and catchment boundaries
ggplot() +
layer_spatial(log10(carea)) +
scale_fill_continuous(expression(log[10](italic(CA))),
low = "lightgrey", high = "black") +
geom_sf(data = contours, color = "yellow") +
geom_sf(data = outlets, pch = 21, bg = "red") +
geom_sf(data = sf::st_geometry(catchments), fill = NA, color = "red") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
# map with elevation contours, channels, outlets and catchment boundaries
ggplot() +
geom_sf(data = contours, color = "grey") +
geom_sf(data = outlet_ll, pch = 21, bg = "red") +
geom_sf(data = sf::st_geometry(catchments), fill = NA, color = "red") +
geom_sf(data = sf::st_geometry(channels), color = "blue") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
3. Compute topographic wetness index (combines ch_saga and raster functions)
The topographic wetness index was classically defined as follows:
\[ TWI = ln(a/tan \beta) \]
where \(a\) is the specific contributing area, computed as
\[ a = A/w_c \]
where \(A\) is the contributing area for a pixel (m2), \(w_c\) is the contour width (m), and \(\beta\) is slope angle (rad). The contour width is a function of flow direction, equal to the cell width (\(dx\), m) for flow in the N, E, S or W directions, or \(\sqrt{2} \cdot dx\) for flow in the diagonal directions.
The code chunk below computes \(TWI\) using the classical definition by combining the ch_saga_carea function with the terrain function from the raster package.
# compute and map topographic wetness index
dx = xres(dem_filled)
carea <- ch_saga_carea(dem_filled, saga_wd)
slope <- terrain(dem_filled, opt = 'slope', units = 'radians', neighbours = 8)
flowdir <- terrain(dem_filled, opt = 'flowdir')
wc <- calc(flowdir, function(x) ifelse(x %in% c(1, 4, 16, 64), dx, dx*sqrt(2)))
twi <- log((carea/wc)/tan(slope))
ggplot() +
layer_spatial(twi) +
scale_fill_continuous(expression(italic(TWI)),
low = "lightgrey", high = "black") +
geom_sf(data = contours, color = "yellow") +
geom_sf(data = outlets, pch = 21, bg = "red") +
geom_sf(data = sf::st_geometry(catchments), fill = NA, color = "red") +
ggspatial::annotation_north_arrow(style = north_arrow_fancy_orienteering,
location = "tr",
pad_x = unit(4, "mm"),
pad_y = unit(6.5, "mm")) +
ggspatial::annotation_scale() +
coord_sf(xlim = c(bb[1], bb[3]), ylim = c(bb[2], bb[4])) +
theme_bw()
ch_clear_wd(saga_wd)## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 6.4.0 (64 bit)
##
## library path: C:\Progra~1\SAGA-GIS\tools\
## library name: ta_hydrology
## library : ta_hydrology
## tool : Flow Accumulation (Top-Down)
## identifier : 0
## author : O.Conrad (c) 2001-2016, Portions by T.Grabs (c) 2010
## processors : 4 [4]
##
## Loading grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/dem.sgrd...
##
##
## Parameters
##
##
## Grid system: 25; 531x 812y; 482544.659487x 5496270.907952y
## Elevation: dem_ns
## Sink Routes: <not set>
## Weights: <not set>
## Flow Accumulation: Flow Accumulation
## Input for Mean over Catchment: <not set>
## Material for Accumulation: <not set>
## Step: 1
## Flow Accumulation Unit: cell area
## Flow Path Length: <not set>
## Channel Direction: <not set>
## Method: Multiple Flow Direction
## Thresholded Linear Flow: no
## Convergence: 1.100000
##
## Create index: dem_ns
## Saving grid: C:/Extension/R/CSHShydRology/ch_saga/scratch/carea.sgrd...