R/ch_wbt_catchment.R

In CSHShydRology: Canadian Hydrological Analyses

#' Delineate catchment boundaries
#'
#' @param fn_pp_snap Name of file containing snapped pour points
#' @param fn_flowdir Name of file containing flow accumulations.
#' @param fn_catchment_ras Raster file to contain delineated catchment.
#' @param fn_catchment_vec Vector file to contain delineated catchment.
#' @param return_vector If \code{TRUE} (the default) a vector of the catchment will be returned.
#' 
#' @author Dan Moore and Kevin Shook
#' @importFrom raster raster
#' @importFrom whitebox wbt_watershed wbt_raster_to_vector_polygons
#' @importFrom sf st_crs write_sf st_read
#' @importFrom magrittr %>%
#' @return If \code{return_vector == TRUE} a vector of the catchment is returned. Otherwise
#' nothing is returned.
#' @export
#' @seealso  \code{\link{ch_wbt_catchment_onestep}}
#' @examples 
#' # Only proceed if Whitebox executable is installed
#' library(whitebox)
#' if (check_whitebox_binary()){
#'   library(raster)
#'   test_raster <- ch_volcano_raster()
#'   dem_raster_file <- tempfile(fileext = ".tif")
#'   no_sink_raster_file <- tempfile("no_sinks", fileext = ".tif")
#' 
#'   # write test raster to file
#'   writeRaster(test_raster, dem_raster_file, format = "GTiff")
#' 
#'   # remove sinks
#'   removed_sinks <- ch_wbt_removesinks(dem_raster_file, no_sink_raster_file, method = "fill")
#' 
#'   # get flow accumulations
#'   flow_acc_file <- tempfile("flow_acc", fileext = ".tif")
#'   flow_acc <- ch_wbt_flow_accumulation(no_sink_raster_file, flow_acc_file)
#' 
#'   # get pour points
#'   pourpoint_file <- tempfile("volcano_pourpoints", fileext = ".shp")
#'   pourpoints <- ch_volcano_pourpoints(pourpoint_file)
#'   snapped_pourpoint_file <- tempfile("snapped_pourpoints", fileext = ".shp")
#'   snapped_pourpoints <- ch_wbt_pourpoints(pourpoints, flow_acc_file, pourpoint_file,
#'   snapped_pourpoint_file, snap_dist = 10)
#' 
#'   # get flow directions
#'   flow_dir_file <- tempfile("flow_dir", fileext = ".tif")
#'   flow_dir <- ch_wbt_flow_direction(no_sink_raster_file, flow_dir_file)
#'   fn_catchment_ras <- tempfile("catchment", fileext = ".tif")
#'   fn_catchment_vec <- tempfile("catchment", fileext = ".shp")
#'   catchments <- ch_wbt_catchment(snapped_pourpoint_file, flow_dir_file, 
#'   fn_catchment_ras, fn_catchment_vec)
#' } else {
#'   message("Examples not run as Whitebox executable not found")
#' }

ch_wbt_catchment <- function(fn_pp_snap, fn_flowdir, fn_catchment_ras, 
                             fn_catchment_vec, return_vector = TRUE) {
  ch_wbt_check_whitebox()
  if (!file.exists(fn_pp_snap)) {
    stop("Error: file containing snapped pour points does not exist")
  }
  if (!file.exists(fn_flowdir)) {
    stop("Error: input flow direction file does not exist")
  }
  
  message("ch_wbt: Delineating catchment boundaries")
  crs_pp <- sf::st_crs(st_read(fn_pp_snap))
  crs_fd <- sf::st_crs(raster(fn_flowdir))
  if (crs_pp != crs_fd) {
    stop("Error: pour points and flow direction grid have different crs")
  }
  wbt_watershed(d8_pntr = fn_flowdir, pour_pts = fn_pp_snap, 
                          output = fn_catchment_ras)
  wbt_raster_to_vector_polygons(fn_catchment_ras, fn_catchment_vec)
  catchment_vec <- st_read(fn_catchment_vec) %>% st_as_sf()
  if (is.na(st_crs(catchment_vec))) {
    sf::st_crs(catchment_vec) <- st_crs(raster(fn_catchment_ras))
    write_sf(catchment_vec, fn_catchment_vec)
  }
  if (return_vector) {
    return(st_read(fn_catchment_vec))
  } else {
    return()
  }
}