diff --git a/DESCRIPTION b/DESCRIPTION index 2c210f9..eedeea9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: watershed Title: Tools for Watershed Delineation -Version: 0.1.0.9007 +Version: 0.2.0 Authors@R: person(given = "Matthew", family = "Talluto", diff --git a/NAMESPACE b/NAMESPACE index 66a8808..a287be4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(delineate) +export(extract_reach) export(fill_dem) export(find_grass) export(pixel_topology) diff --git a/R/grass.r b/R/grass.r index 161d75c..d14ff15 100644 --- a/R/grass.r +++ b/R/grass.r @@ -42,7 +42,7 @@ #' @name read_rasters -#' @rdname .read_rasters +#' @rdname read_rasters #' @title Read files from grass #' Read and format raster and vector layers from a grass session #' @param layers A vector of names of rasters to read @@ -63,7 +63,7 @@ ras } -#' @rdname .read_rasters +#' @rdname read_rasters #' @keywords internal .read_vector = function(x) { v = rgrass7::readVECT(x) diff --git a/R/stream.r b/R/stream.r index c7edb45..66b82ff 100644 --- a/R/stream.r +++ b/R/stream.r @@ -4,13 +4,15 @@ #' @param pretty Should flat areas be made prettier? Corresponds to r.watershed -b flag #' @param file The file name of the raster to be returned, see `details`. #' @param outlet The location of the outlet of the target stream, see `details`. +#' @param reach_len Optional; if provided, reaches will be resized to meet this target length +#' @param ... If `reach_len` is specified, additional parameters to be passed to [resize_reaches()] #' @details This is a wrapper for [r.watershed](https://grass.osgeo.org/grass76/manuals/r.watershed.html). #' #' The threshold parameter controls the level of detail in the delineated streams. Higher thresholds result in faster #' computation and fewer streams in the output. For finer control, see [extract_stream()]. #' -#' If the `outlet` parameter is missing (the default), then all streams in the DEM will be returned. If `outlet` is `NA`, -#' then the largest stream in the area will be set as the outlet, and only streams in that watershed will be used. +#' IIf `outlet` is `NA` (the default), then the largest stream in the area will be set as the outlet, +#' and only streams in that watershed will be used. #' If a smaller stream is the focus, then outlet can be a pair of x-y coordinates which will determine the farthest #' downstream point to use. #' @@ -20,8 +22,13 @@ #' file format; e.g., .tif, .grd). If not specified, a temp file will be created and will be #' lost at the end of the R session. #' @return A [raster::stack()], containing the drainage map, the flow accumulation map, and the delineated streams. +#' @examples +#' \donttest{ +#' data(kamp_dem) +#' x = delineate(kamp_dem) +#' } #' @export -delineate = function(dem, threshold = 1e6, pretty = FALSE, file, outlet) { +delineate = function(dem, threshold = 1e6, pretty = FALSE, file, outlet=NA, reach_len, ...) { cell_area = prod(raster::res(dem)) threshold_cell = floor(threshold/cell_area) @@ -47,32 +54,39 @@ delineate = function(dem, threshold = 1e6, pretty = FALSE, file, outlet) { ws_env$rasters = c(ws_env$rasters, accum, drainage, stream) ## gather data - res = .read_rasters(c(accum, drainage, stream), file) - - ## Clip to a single network if desired - if(!missing(outlet)) { - if(length(outlet) == 1 && is.na(outlet)) - outlet = raster::coordinates(res)[which.max(abs(raster::values(res$accum))),] - - rgrass7::execGRASS("r.water.outlet", flags=c("overwrite", "quiet"), input=drainage, output = "catchment", - coordinates = outlet) - ws_env$rasters = c(ws_env$rasters, "catchment") - catchment = .read_rasters("catchment") - catchment[catchment == 0] = NA - stream_clip = raster::mask(res$stream, catchment) - res = raster::dropLayer(res, which(names(res) == stream)) - res = raster::addLayer(res, stream = stream_clip, catchment = catchment) - if(!missing(file)) - res = raster::writeRaster(res, filname = file) - } - - ## renumber reaches to go from 1:nreaches - res[['stream']] = raster::match(res[['stream']], + res = .read_rasters(c(accum, drainage, stream)) + + ## Clip to a single network + if(length(outlet) == 1 && is.na(outlet)) + outlet = raster::coordinates(res)[which.max(abs(raster::values(res$accum))),] + + rgrass7::execGRASS("r.water.outlet", flags=c("overwrite", "quiet"), input=drainage, output = "catchment", + coordinates = outlet) + ws_env$rasters = c(ws_env$rasters, "catchment") + catchment = .read_rasters("catchment") + catchment[catchment == 0] = NA + catchment = raster::trim(catchment) + res = raster::crop(res, catchment) + stream_clip = raster::mask(res$stream, catchment) + res = raster::dropLayer(res, which(names(res) == stream)) + res = raster::addLayer(res, stream = stream_clip, catchment = catchment) + + ## renumber reaches to go from 1:nreaches, make sure streams are integers + res[['stream']] = raster::match(res[['stream']], raster::unique(res[['stream']])) + storage.mode(res[['stream']][]) = 'integer' + + if(!missing(reach_len)) { + Tp = pixel_topology(res) + res[['stream']] = resize_reaches(res[['stream']], Tp, len = reach_len, ...) + } ## clean up files .clean_grass() + if(!missing(file)) + res = raster::writeRaster(res, filname = file) + res } @@ -93,22 +107,22 @@ vectorise_stream = function(x) { inpname = "stream" outname = "stream_thin" vectname = "stream_vect" - + .start_grass(x, inpname) flags = c("overwrite", "quiet") rgrass7::execGRASS("r.thin", flags=flags, input=inpname, output = outname) ws_env$rasters = c(ws_env$rasters, outname) - + rgrass7::execGRASS("r.to.vect", flags=flags, input=outname, output = vectname, type="line", column="reach_id") ws_env$vectors = c(ws_env$vectors, vectname) vect = .read_vector(vectname) - + ## todo: v.clean .clean_grass() - vect + vect } #' Resize reaches from a delineated stream @@ -116,43 +130,59 @@ vectorise_stream = function(x) { #' @param len The target length of resized reaches (in map units) #' @param min_len How small is the smallest acceptable reach? #' @param trim Logical; remove headwater reaches smaller than min_length? -#' @param Tp Optional; a pixel topology; will be computed if missing -#' @param Tr Optional; a reach topology; will be computed if missing +#' @param Tp A pixel topology #' @return A raster similar to `x` with reaches resized #' @examples #' \donttest{ #' data(kamp_dem) #' x = delineate(kamp_dem) -#' x[['stream']] = resize_reaches(x[['stream']], 1000, 600) +#' Tp = pixel_topology(x) +#' x[['stream']] = resize_reaches(x[['stream']], Tp, 1000, 600) #' } #' @export -resize_reaches = function(x, len, min_len = 0.5 * len, trim = TRUE, Tp, Tr) { - +resize_reaches = function(x, Tp, len, min_len = 0.5 * len, trim = TRUE) { + + if(!is.null(getOption("mc.cores")) && getOption("mc.cores") > 1) { + lapplfun = parallel::mclapply + } else { + lapplfun = lapply + } + + ids = raster::unique(x) + pids = lapplfun(ids, extract_reach, stream=x, Tp = Tp, sorted = TRUE) + new_reaches = unlist(lapplfun(pids, .resize_r, Tp = Tp, trim = trim, len = len, min_len = min_len), + recursive = FALSE) + x = x * NA + for(i in seq_along(new_reaches)) + x[new_reaches[[i]]] = i + + storage.mode(x[]) = 'integer' + x } -# split a single reach -# returns a list of vectors -# each vector is a new reach -# pixels in the old reach that are not in the new reach should be set to NA by the calling function -.resize_r = function(Tp) { - +#' Split a single reach +#' @param ids A vector of pixel ids in the reach, sorted from upstream to downstream +#' @param Tp A pixel topology +#' @param trim Logical, should we trim small headwater reaches? +#' @param len Target reach length +#' @param min_len Minimum acceptable reach length, only used if `trim` is `TRUE` +#' @return A list of vectors, each vector is a single new reach. +#' @keywords internal +.resize_r = function(ids, Tp, trim, len, min_len) { + ## only trim headwater reaches - trim = trim && .is_headwater(Tr, i) - - ## extract the IDs of the reach of interest - ## then find the reach outlet - pids = which(raster::values(stream) == i) - Tp_r = Tp[pids, pids, drop=FALSE] - bottom = pids[.outlet(Tp_r)] - top = pids[.headwater(Tp_r)] - - pix_lens = Matrix::rowSums(Tp) - - ## number of reaches - rlen = sum(pix_lens[pids]) + trim = trim && .is_headwater(Tp)[ids[1]] + + # we will resize from bottom to top, then reverse again at the end + ids = rev(ids) + + pix_lens = Matrix::rowSums(Tp[ids,]) + + ## compute target number of reaches + rlen = sum(pix_lens) nr = ceiling(rlen/len) - - ## if not trimming, decide on new reach length + + ## if not trimming, decide on new reach length to get reaches mostly even lengths ## try the computed number of reaches +/- 1 ## pick the one that gets us closest to the target length if(!trim) { @@ -161,42 +191,26 @@ resize_reaches = function(x, len, min_len = 0.5 * len, trim = TRUE, Tp, Tr) { len = cand[which.min(ldiff)] nr = nr + (-1:1)[which.min(ldiff)] } - - ## biggest possible reach, in terms of the number of pixels - maxpix = ceiling(len / min(Tp@x)) - - - ## traverse, accumulating each reach as we go upstream - ## start with memory pre-allocated, will drop excess 0s later - reaches = lapply(rep(maxpix, nr), numeric) - cpix = bottom - r = 1 - while(cpix != top) { - p = 1 - rl = 0 - while(rl < len) { - rl = rl + pix_lens[cpix] - reaches[[r]][p] = cpix - p = p+1 - cpix = .upstream(Tp, cpix) - if(cpix == top) - break() - } - r = r+1 - } - - ## drop 0s - reaches = lapply(reaches, function(r) r[r>0]) - if(length(reaches[[nr]]) == 0) { - reaches[[nr]] = NULL - nr = nr - 1 + + new_rids = ceiling(cumsum(pix_lens) / len) + + ## outlet of the whole river network will end up with an id of 0; fix it + new_rids[new_rids == 0] = 1 + + # any pixels given a new id > nr are short headwater pixels that should be dropped + new_rids = new_rids[new_rids <= nr] + + new_reaches = lapply(unique(new_rids), function(i) ids[new_rids == i]) + + ## if we are trimming, drop the headwater if it is too short + if(trim) { + r = new_reaches[[length(new_reaches)]] + hwlen = sum(pix_lens[match(r, ids)]) + if(hwlen < min_len) + new_reaches[[length(new_reaches)]] = NULL } - rlens = lapply(reaches, function(x) sum(pix_lens[x])) - - ## trim the top if it is a headwater - if(trim && rlens[nr] < min_len) - reaches[[nr]] = NULL - - reaches + + # reorder the list so that the reaches are sorted from headwater to outlet + rev(new_reaches) } diff --git a/R/topology.r b/R/topology.r index e7c5e1b..5e937c3 100644 --- a/R/topology.r +++ b/R/topology.r @@ -14,8 +14,9 @@ #' @examples #' \donttest{ #' data(kamp_dem) -#' Tp = pixel_topology(kamp_dem) -#' Tr = reach_topology(kamp_dem, Tp) +#' kamp = delineate(kamp_dem, outlet = NA) +#' Tp = pixel_topology(kamp) +#' Tr = reach_topology(kamp, Tp) #' } NULL @@ -44,7 +45,7 @@ pixel_topology = function(x, drainage, stream) { xy = .flowto(vals, xmax = nc, ymax = nr) # res = xy[,c('from_id', 'drainage')] # colnames(res)[1] = 'id' - r = raster::res(x) + r = raster::res(stream) len = ifelse(xy[, 'drainage'] %in% c(1,3,5,7), sqrt(r[1]^2 + r[2]^2), ifelse(xy[, 'drainage'] %in% c(2, 6), r[2], r[1])) res = Matrix::sparseMatrix(i = xy[,'from_id'], j = xy[,'to_id'], dims=rep(raster::ncell(drainage), 2), x = len) @@ -55,6 +56,9 @@ pixel_topology = function(x, drainage, stream) { #' @rdname topologies #' @export reach_topology = function(x, Tp, stream) { + if(missing(Tp)) + stop("A topology is required for this function") + if(!requireNamespace("Matrix")) stop("This functionality requires the Matrix package; please install it with install.packages('Matrix') and try again.") @@ -74,12 +78,48 @@ reach_topology = function(x, Tp, stream) { adj = lapplfun(rids, .upstream_r, stream = stream, Tx = Tp) adj = do.call(rbind, adj) - Matrix::sparseMatrix(adj[,2], adj[,1], dims=rep(max(rids), 2), - dimnames = list(rids, rids)) + + ## compute the length of each reach as the sum of all pixels + pids = lapplfun(rids, extract_reach, stream=stream, Tp = Tp) + lens = unlist(lapplfun(pids, function(x) sum(Tp[x,,drop=FALSE]))) + + # add a virtual node for the outlet to flow to + outlet = which(!(rids %in% adj[,'from'])) + adj = rbind(adj, c(max(rids)+1, outlet)) + rids = c(rids, max(rids)+1) + + Matrix::sparseMatrix(adj[,'from'], adj[,'to'], dims=rep(max(rids), 2), + dimnames = list(rids, rids), x = lens) } +#' Find pixel IDs matching reach i +#' @param i Reach number to extract +#' @param stream Raster stream layer +#' @param Tp optional Topology, only needed for sorting +#' @param sorted Logical, should the ids be returned sorted from upstream to down? +#' @return A vector of pixel ids +#' #' @examples +#' \donttest{ +#' data(kamp_dem) +#' kamp = delineate(kamp_dem, outlet = NA) +#' Tp = pixel_topology(kamp) +#' extract_reach(5, kamp$stream, Tp) +#' } +#' @export +extract_reach = function(i, stream, Tp, sorted = FALSE) { + pids = which(raster::values(stream) == i) + if(sorted) { + pid_s = numeric(length(pids)) + Tp_r = Tp[pids, pids, drop=FALSE] + pid_s[1] = .headwater(Tp_r) + for(i in 2:length(pid_s)) + pid_s[i] = .downstream(Tp_r, pid_s[i-1]) + pids = pids[pid_s] + } + pids +} @@ -112,8 +152,7 @@ NULL #' @rdname upstream #' @keywords internal .headwater = function(Tx) { - ## check for rowsums as well because we exclude nodes that are connected to nothing - which(Matrix::colSums(Tx) == 0 & Matrix::rowSums(Tx) != 0) + which(.is_headwater(Tx)) } #' @rdname upstream @@ -122,13 +161,20 @@ NULL which(Matrix::rowSums(Tx) == 0 & Matrix::colSums(Tx) != 0) } +#' @rdname upstream +#' @keywords internal +.is_headwater = function(Tx) { + ## check for rowsums as well because we exclude nodes that are connected to nothing + Matrix::colSums(Tx) == 0 & Matrix::rowSums(Tx) != 0 +} + #' @rdname upstream #' @return For `.upstream_r`, A named2-column matrix, 'to' is the downstream reach (i.e., i) and #' 'from' contains ids of the reach(es) upstream of i; or NULL if there are none #' @keywords internal .upstream_r = function(i, stream, Tx) { - pids = which(raster::values(stream) == i) + pids = extract_reach(i, stream) Tp_r = Tx[pids, pids, drop=FALSE] r_top = pids[.headwater(Tp_r)] r_nb = .upstream(Tx, r_top) diff --git a/inst/testdata/kamp_sm.rds b/inst/testdata/kamp_sm.rds index 16b1d9d..6cf7cb6 100644 Binary files a/inst/testdata/kamp_sm.rds and b/inst/testdata/kamp_sm.rds differ diff --git a/man/delineate.Rd b/man/delineate.Rd index fd234e3..d20bddb 100644 --- a/man/delineate.Rd +++ b/man/delineate.Rd @@ -4,7 +4,15 @@ \alias{delineate} \title{Delineate streams from an elevation model} \usage{ -delineate(dem, threshold = 1e+06, pretty = FALSE, file, outlet) +delineate( + dem, + threshold = 1e+06, + pretty = FALSE, + file, + outlet = NA, + reach_len, + ... +) } \arguments{ \item{dem}{A \link[raster:raster]{raster::raster}; a digital elevation model.} @@ -16,6 +24,10 @@ delineate(dem, threshold = 1e+06, pretty = FALSE, file, outlet) \item{file}{The file name of the raster to be returned, see \code{details}.} \item{outlet}{The location of the outlet of the target stream, see \code{details}.} + +\item{reach_len}{Optional; if provided, reaches will be resized to meet this target length} + +\item{...}{If \code{reach_len} is specified, additional parameters to be passed to \code{\link[=resize_reaches]{resize_reaches()}}} } \value{ A \code{\link[raster:stack]{raster::stack()}}, containing the drainage map, the flow accumulation map, and the delineated streams. @@ -29,8 +41,8 @@ This is a wrapper for \href{https://grass.osgeo.org/grass76/manuals/r.watershed. The threshold parameter controls the level of detail in the delineated streams. Higher thresholds result in faster computation and fewer streams in the output. For finer control, see \code{\link[=extract_stream]{extract_stream()}}. -If the \code{outlet} parameter is missing (the default), then all streams in the DEM will be returned. If \code{outlet} is \code{NA}, -then the largest stream in the area will be set as the outlet, and only streams in that watershed will be used. +IIf \code{outlet} is \code{NA} (the default), then the largest stream in the area will be set as the outlet, +and only streams in that watershed will be used. If a smaller stream is the focus, then outlet can be a pair of x-y coordinates which will determine the farthest downstream point to use. @@ -40,3 +52,9 @@ It is recommended to specify the \code{file} parameter (including the extension file format; e.g., .tif, .grd). If not specified, a temp file will be created and will be lost at the end of the R session. } +\examples{ +\donttest{ + data(kamp_dem) + x = delineate(kamp_dem) +} +} diff --git a/man/dot-resize_r.Rd b/man/dot-resize_r.Rd new file mode 100644 index 0000000..d6572fb --- /dev/null +++ b/man/dot-resize_r.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stream.r +\name{.resize_r} +\alias{.resize_r} +\title{Split a single reach} +\usage{ +.resize_r(ids, Tp, trim, len, min_len) +} +\arguments{ +\item{ids}{A vector of pixel ids in the reach, sorted from upstream to downstream} + +\item{Tp}{A pixel topology} + +\item{trim}{Logical, should we trim small headwater reaches?} + +\item{len}{Target reach length} + +\item{min_len}{Minimum acceptable reach length, only used if \code{trim} is \code{TRUE}} +} +\value{ +A list of vectors, each vector is a single new reach. +} +\description{ +Split a single reach +} +\keyword{internal} diff --git a/man/extract_reach.Rd b/man/extract_reach.Rd new file mode 100644 index 0000000..c73c9e3 --- /dev/null +++ b/man/extract_reach.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/topology.r +\name{extract_reach} +\alias{extract_reach} +\title{Find pixel IDs matching reach i} +\usage{ +extract_reach(i, stream, Tp, sorted = FALSE) +} +\arguments{ +\item{i}{Reach number to extract} + +\item{stream}{Raster stream layer} + +\item{Tp}{optional Topology, only needed for sorting} + +\item{sorted}{Logical, should the ids be returned sorted from upstream to down?} +} +\value{ +A vector of pixel ids +#' @examples +\donttest{ + data(kamp_dem) + kamp = delineate(kamp_dem, outlet = NA) + Tp = pixel_topology(kamp) + extract_reach(5, kamp$stream, Tp) +} +} +\description{ +Find pixel IDs matching reach i +} diff --git a/man/read_rasters.Rd b/man/read_rasters.Rd new file mode 100644 index 0000000..59a3ae1 --- /dev/null +++ b/man/read_rasters.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/grass.r +\name{read_rasters} +\alias{read_rasters} +\alias{.read_rasters} +\alias{.read_vector} +\title{Read files from grass +Read and format raster and vector layers from a grass session} +\usage{ +.read_rasters(layers, file) + +.read_vector(x) +} +\arguments{ +\item{layers}{A vector of names of rasters to read} + +\item{file}{The file name to save the raster} + +\item{x}{A single vector layer name to read} +} +\description{ +Read files from grass +Read and format raster and vector layers from a grass session +} +\keyword{internal} diff --git a/man/resize_reaches.Rd b/man/resize_reaches.Rd index 8298c3f..1b8159e 100644 --- a/man/resize_reaches.Rd +++ b/man/resize_reaches.Rd @@ -4,20 +4,18 @@ \alias{resize_reaches} \title{Resize reaches from a delineated stream} \usage{ -resize_reaches(x, len, min_len = 0.5 * len, trim = TRUE, Tp, Tr) +resize_reaches(x, Tp, len, min_len = 0.5 * len, trim = TRUE) } \arguments{ \item{x}{A \link{raster} stream map, such as one created by \code{\link[=delineate]{delineate()}}.} +\item{Tp}{A pixel topology} + \item{len}{The target length of resized reaches (in map units)} \item{min_len}{How small is the smallest acceptable reach?} \item{trim}{Logical; remove headwater reaches smaller than min_length?} - -\item{Tp}{Optional; a pixel topology; will be computed if missing} - -\item{Tr}{Optional; a reach topology; will be computed if missing} } \value{ A raster similar to \code{x} with reaches resized @@ -29,6 +27,7 @@ Resize reaches from a delineated stream \donttest{ data(kamp_dem) x = delineate(kamp_dem) - x[['stream']] = resize_reaches(x[['stream']], 1000, 600) + Tp = pixel_topology(x) + x[['stream']] = resize_reaches(x[['stream']], Tp, 1000, 600) } } diff --git a/man/topologies.Rd b/man/topologies.Rd index 9baa9c1..4baa5f6 100644 --- a/man/topologies.Rd +++ b/man/topologies.Rd @@ -36,7 +36,8 @@ upstream pixel/reach. \examples{ \donttest{ data(kamp_dem) - Tp = pixel_topology(kamp_dem) - Tr = reach_topology(kamp_dem, Tp) + kamp = delineate(kamp_dem, outlet = NA) + Tp = pixel_topology(kamp) + Tr = reach_topology(kamp, Tp) } } diff --git a/man/upstream.Rd b/man/upstream.Rd index cdab6aa..13f0eb4 100644 --- a/man/upstream.Rd +++ b/man/upstream.Rd @@ -6,6 +6,7 @@ \alias{.downstream} \alias{.headwater} \alias{.outlet} +\alias{.is_headwater} \alias{.upstream_r} \title{Simple topology analysis Simple functions for extracting information from river topologies} @@ -18,6 +19,8 @@ Simple functions for extracting information from river topologies} .outlet(Tx) +.is_headwater(Tx) + .upstream_r(i, stream, Tx) } \arguments{ diff --git a/vignettes/watershed.Rmd b/vignettes/watershed.Rmd index 2888429..71cb784 100644 --- a/vignettes/watershed.Rmd +++ b/vignettes/watershed.Rmd @@ -36,51 +36,24 @@ Many workflows advise either using a "burned-in" DEM or a "filled" DEM to avoid The pared-down workflow offered by `watershed` provides a relatively quick way to go from a dem to a stream map. We will make a first-draft map using the `kamp_dem` dataset provided by this package, for the Kamp river in Austria. The DEM is a `raster`, so we load this package as well. -```{r dem} -library(raster) -library(sf) -data(kamp_dem) -plot(kamp_dem, col=terrain.colors(20), axes = FALSE) -``` +The `delineate()` function is the main workhorse here. The resulting raster stack contains 3 layers: `accum` is the flow accumulation, `drainage` is the drainage direction, and `stream` is the stream map. We also produce a vector layer of the stream map for plotting. Some notes about our use of `delineate` here: -The `delineate()` function is the main workhorse here. The resulting raster stack contains 3 layers: `accum` is the flow accumulation, `drainage` is the drainage direction, and `stream` is the stream map. We also produce a vector layer of the stream map for plotting. +* We do not set the `threshold` parameter, which is the minimum catchment area before something is marked as a stream. The default of 1e6 $m^2$ (10 $km^2$) works well; increasing the threshold will result in more streams, decreasing it will take longer and make more streams. +* We restrict the output to a single river basin with `outlet = NA` (this chooses the largest basin). To select a different basin, we could instead use `outlet = c(4704588, 2847762)`; the coordinates must be exactly on the stream for this to work. +* Because `outlet` is provided, the output rasters also have their extent reduced to the catchment of our stream. ```{r delineate, cache = TRUE} -kamp = delineate(kamp_dem) -kv = vectorise_stream(kamp[["stream"]]) -plot(kamp_dem, col=terrain.colors(20), axes = FALSE) -plot(kv, col='blue', add = TRUE) -``` - -The defaults produce pretty good results in this case. Perhaps we know from on-the-ground work, or aerial photo inspection, that there are more streams that what is produced here. We can improve this with the `threshold` parameter, which sets the default minimum basin size. This defaults to 1e6 m^2, or 1 km^2. Here we reduce it by half; note that this gives us a warning about possible long computation times; it's not such a problem here, because the Ybbs is a relatively small area, but with large DEMs this can make the computation very slow. - -```{r delineate2, cache = TRUE} -kamp2 = delineate(kamp_dem, threshold = 0.5e5) -plot(kamp_dem, col=terrain.colors(20), axes = FALSE) -kv = vectorise_stream(kamp2[["stream"]]) -plot(kv, col='blue', add = TRUE) - -rm(kamp2) -``` - -This produces a lot of streams! From our fieldwork in the Kamp, we know that the default produces a pretty good network, so we will keep the first one. - -## 2. Selecting a single basin - -We would like to restrict our stream map to a single catchment. `delineate` offers this option; here we set `outlet = NA`, which will use the largest stream (with a catchment completely contained within the DEM) as the outlet for the desired basin. Alternatively, we could zoom in on the layer and choose the x-y coordinates of the outlet we wanted to use, and provide these as a vector in outlet. - - -```{r crop, cache = TRUE} +library(raster) +library(sf) +data(kamp_dem) kamp = delineate(kamp_dem, outlet = NA) -## alternative - note that the coordinates must be exactly on the delineated stream -## for this to work -# kamp = delineate(kamp_dem, outlet = c(4704588, 2847762)) kv = vectorise_stream(kamp[["stream"]]) plot(kamp_dem, col=terrain.colors(20), axes = FALSE) -plot(kv, col='blue', add = TRUE) +plot(st_geometry(kv), col='blue', add = TRUE) ``` -## 3. Generating topologies + +## 2. Generating topologies Many later features require a topology; this simply a square matrix `T`, with one row/column for each pixel or reach in the river network. `T[i,j] == 1` indicates that pixel/reach `i` is directly upstream of pixel/reach `j` (in other words, `j` gets water from `i` without that water flowing through any other intermediate pixels/reaches). @@ -89,10 +62,6 @@ Generating the topology can take some time for large and/or complex river networ Our first step is to crop the rasters from the previous steps to only include our catchment of interest. Then we load the Matrix package to handle the `sparseMatrix` returned by `pixel_topology`. Finally, we create two topologies; one for the pixels, and one for reaches. ```{r topol, cache=TRUE, warning=TRUE} -catch = kamp$catchment -catch = trim(catch) -kamp = crop(kamp, catch) - library(Matrix) kamp_Tp = pixel_topology(kamp) kamp_Tr = reach_topology(kamp, kamp_Tp) @@ -101,17 +70,32 @@ kamp_Tr = reach_topology(kamp, kamp_Tp) Note the warning, indicating that there is a confluence where more than 2 streams flow into the same pixel. This is not a serious problem, so we will ignore it for now. +## 3. Producing smaller reaches -## 4. Revising the delineation - +The default output from `delineate` is to assign each pixel to a reach, where a reach is a stretch of river between confluences. For many applications, it is useful to split these into smaller sections. `watershed` includes a function to split reaches into fixed (approximately) equal-length sections. -## 4. Producing smaller reaches - -The default output from `delineate` is to assign each pixel to a reach, where a reach is a stretch of river between confluences. For many applications, it is useful to split these into smaller sections. `watershed` includes a function to split reaches into fixed (approximately) equal-length sections. Unequal-length reaches is a feature that may be added if we receive enough requests. +* Note that resizing the reaches invalidates the old topologies, and a new ones must be made. +* Alternatively, we could define the reach lengths directly in `delineate` by setting the `reach_len` parameter. +* Setting `trim=TRUE` will remove small (by default, length < target_length/2) headwater reaches from the network ```{r reach_split, cache=TRUE} - - +new_reaches = resize_reaches(kamp$stream, kamp_Tp, 800, trim=TRUE) +new_v = vectorise_stream(new_reaches) +new_Tp = pixel_topology(drainage = kamp$drainage, stream = new_reaches) +new_Tr = reach_topology(Tp = new_Tp, stream = new_reaches) + +## how many reaches, and what range of lengths? +c(nrow(kamp_Tr), nrow(new_Tr)) +library(ggplot2) +pl = ggplot(data.frame(len = c(kamp_Tr@x, new_Tr@x), + reach = c(rep("old", length(kamp_Tr@x)), rep("new", length(new_Tr@x)))), aes(x=len, fill=reach)) +pl + geom_histogram(aes(y=..density..), position="dodge", binwidth = 200, colour="gray") + theme_minimal() + xlab("Reach length (m)") +## update the raster stack to use the new streams +kamp[["stream"]] = new_reaches +kv = new_v +kamp_Tp = new_Tp +kamp_Tp = new_Tr +rm(list=c("new_reaches", "new_v", "new_Tp", "new_Tr")) ```