# 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: sp
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(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:

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_wd

ch_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_wd

ch_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_fillsinks

ch_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_removesinks

ch_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_carea

ch_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_catchment

ch_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:

  1. a shape file containing the drainage network,
  2. a grid containing flow directions, and
  3. 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_channels

ch_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_checkcatchment

ch_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_checkchannels

Application - 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':
## 
##     extract
library(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-1
library(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, union
library(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 demNULL
ch_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^2
channel_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:    NA
channels = 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()