Prelim work on reach topology

Author: Matthew Talluto

R/topology.r

@@ -25,7 +25,7 @@ reach_topology = function(x, Tp, stream) {
 	if(!all.equal(rids, 1:length(rids)))
 		stop("ReachIDs not in strict ascending order, they must be renumbered")
 
-	adj = lapplfun(reaches, function(r) .reach_upstream(ws, r))
+	adj = lapplfun(rids, .upstream_r, stream = stream, Tp = Tp)
 	adj = do.call(rbind, adj)
 	stop("check the below; I think columns are reversed")
 	adjspMat = Matrix::sparseMatrix(adjMat[,1], adjMat[,2], dims=rep(max(ws$data$reachID), 2),
@@ -34,30 +34,59 @@ reach_topology = function(x, Tp, stream) {
 
 }
 
-#' Find all reaches directly upstream from a given reach
+#' Find all reaches directly upstream from a given reach using a pixel topology
+#' @param stream A stream raster (as produced by [delineate()]) with reachIDs as the values
+#' @param Tp Pixel topology for the stream raster
+#' @param rid The focal reach ID
+#' @return A named2-column matrix, 'to' is the downstream reach (i.e., rid) and
+#' 'from' contains ids of the reach(es) upstream of rid, or NULL if there are none
 #' @keywords internal
-.reach_upstream <- function(Tp, rid) {
-
-	stop("update")
-	ids <- which(ws$data$reachID == rch)
-	reachAdj <- ws$adjacency[ids,ids, drop = FALSE]
-	mostUpstream <- ids[which(Matrix::rowSums(reachAdj) == 0)]
-	adjMatUp <- as.matrix(ws$adjacency[mostUpstream,])
-	upRch <- which(adjMatUp == 1)
-	if(length(upRch) > 0) {
-		upRch <- ws$data$reachID[upRch]
-		return(cbind(rch, upRch))
-	} else return(NULL)
+.upstream_r = function(rid, stream, Tp) {
+	pids = which(values(stream) == rid)
+	Tp_r = Tp[pids, pids, drop=FALSE]
+	r_top = .headwater(Tp_r)
+	r_nb = .upstream(Tp, r_top)
+	if(length(r_nb) > 0) {
+		from_ids = stream[r_nb]
+		res = cbind(to = rid, from = from_ids)
+	} else {
+		res = NULL
+	}
+	res
+}
+
+#' @keywords internal
+.headwater = function(Tx) {
+	## check for rowsums as well because we exclude nodes that are connected to nothing
+	which(colSums(Tx) == 0 & rowSums(Tx) != 0)
 }
 
+#' Simple topological functions
+#' @details upstream/downstream take a single node as input, and return the up- or downstream node.
+#' Outlets/headwaters analyzes the entire toplogy and finds nodes that have no downstream or upstream neighbours.
+#' @param Tx A (pixel or reach) topology
+#' @param i,j Focal node
+#' @return The id (corresponding to dims of Tx) of the desired node(s)
+#' @keywords internal
+.upstream = function(Tx, j) {
+	which(Tx[,j] != 0)
+}
 
+#' @keywords internal
+.downstream = function(Tx, i) {
+	which(Tx[,j] != 0)
+}
 
 
-#' Buid a topology for each pixel in a delineated stream
+#' Buid a topology for a delineated stream
 #' @param x A [raster::stack], such as one created by [delineate()], or specify layers separately with `drain` and `stream`.
 #' @param drainage Optional, ignored if `x` is provided; a drainage direction raster
 #' @param stream Optional, ignored if `x` is provided; a delineated stream raster, all non-stream cells must be `NA`
-#' @return A [Matrix::sparseMatrix] giving the pixel topology
+#' @details The topology is a square matrix showing immediate adjacency between pixels/reaches. Each row/column
+#' in the topology is a single pixel/reach. Zero entries indicate no adjacency between two nodes. A non-zero entry
+#' for row `i` column `j` indicates that `i` is upstream of `j`. The value in the entry is the length of the 
+#' upstream pixel/reach.
+#' @return A [Matrix::sparseMatrix] giving the pixel or reach topology
 #' @export
 pixel_topology = function(x, drainage, stream) {
 	if(!requireNamespace("Matrix"))

R/zzz.r

@@ -2,7 +2,7 @@
 	Sys.setenv("GRASS_VERBOSE"=0)
 	rgrass7::use_sp()
 
-	if(requireNamespace("parallel") & !grepl("[Ww]indows", Sys.info()['sysname'])) {
+	if(is.null(getOption("mc.cores")) && !grepl("[Ww]indows", Sys.info()['sysname'])) {
 		message("For faster performance of some fucntions on a machine with lots of RAM, you can set")
 		message("options(mc.cores = parallel::detectCores())")
 	}

README.md

@@ -12,7 +12,7 @@ Before installing this package, you should install [GRASS GIS](https://grass.osg
 
 **Mac**: `watershed` can auto-detect [grass binary](http://grassmac.wikidot.com) installations (which are usually installed in `/Applications/`). Command-line installations may also work; try following the linux instructions.
 
-**Linux/Unix**: Install grass however you like. For auto-detection to work, make sure the grass binary (`grass74`, `grass76`, `grass78`) is in your `$PATH` and visible from within R; try `system2('grass76', args = c("--config", "path"))` to see if this is working, it should output the location of your grass installation. If you want a different version of grass, you can alias `grass78` to whatever binary you like, or see **Changing the default grass* below.
+**Linux/Unix**: Install grass however you like. For auto-detection to work, make sure the grass binary (`grass74`, `grass76`, `grass78`) is in your `$PATH` and visible from within R; try, e.g., `system2('grass76', args = c("--config", "path"))` to see if this is working, it should output the location of your grass installation. If you want a different version of grass, you can alias `grass78` to whatever binary you like, or see **Changing the default grass* below.
 
 ## Installing `watershed`.