Reach resizing, and other features and bugfixes

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",

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)

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)

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,66 +107,82 @@ 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
 #' @param x A [raster] stream map, such as one created by [delineate()].
 #' @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)
 }