Watershed Analysis Methods

.Rbuildignore

@@ -20,6 +20,8 @@ README.md
 ^\.github$
 ^.*\.Rproj$
 ^\.Rproj\.user$
+
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
+

.gitignore

@@ -12,7 +12,9 @@ src/Makevars.ucrt.rej
 *.dll
 *.zip
 .Rproj.user
-
+af_acc*
+.Rbuildignore
+.RData
 revdep/
 python/
 cmdapp/terra

NAMESPACE

@@ -21,3 +21,24 @@ S3method(str, SpatGraticule)
 
 export(add_box, add_legend, add_grid, focalMat, extractAlong, gdal, getGDALconfig, graticule, halo, setGDALconfig, map.pal, north, sbar, terraOptions, tmpFiles, makeVRT, mem_info, free_RAM, same.crs, shade, gdalCache, fileBlocksize, vector_layers, vrt_tiles, names, round)
 
+## EC 20210702 ecor
+exportMethods(watershed2)
+## END EC 20210702 ecor
+## EC 20220809
+exportMethods(pitfinder2)
+## END EC 20220809 ecor
+
+## EC 20231031 ecor
+exportMethods(NIDP2)
+## END EC 20231031
+
+## EC 20231031 ecor
+exportMethods(flowAccu2)
+## END EC 20231031
+
+## EC 20231114 ecor
+exportMethods(flowAccu2_weight)
+## END EC 2023114
+
+
+

R/Agenerics.R

@@ -382,3 +382,6 @@ if (!isGeneric("xmax<-")) { setGeneric("xmax<-", function(x, ..., value) standar
 if (!isGeneric("ymin<-")) { setGeneric("ymin<-", function(x, ..., value) standardGeneric("ymin<-"))}
 if (!isGeneric("ymax<-")) { setGeneric("ymax<-", function(x, ..., value) standardGeneric("ymax<-"))}
 if (!isGeneric("zoom")) {setGeneric("zoom", function(x, ...)standardGeneric("zoom"))}
+
+
+

R/Agenerics_watershed2_extesion.R

@@ -0,0 +1,13 @@
+
+## ADDED BY ecor
+
+if (!isGeneric("watershed2")) {setGeneric("watershed2", function(p,pp_offset,...) standardGeneric("watershed2"))}
+
+if (!isGeneric("pitfinder2")) {setGeneric("pitfinder2", function(p, ...) standardGeneric("pitfinder2"))}
+
+
+if (!isGeneric("NIDP2")) {setGeneric("NIDP2", function(p, ...) standardGeneric("NIDP2"))}
+
+if (!isGeneric("flowAccu2")) {setGeneric("flowAccu2", function(p, ...) standardGeneric("flowAccu2"))}
+
+if (!isGeneric("flowAccu2_weight")) {setGeneric("flowAccu2_weight", function(p,weight,...) standardGeneric("flowAccu2_weight"))}

R/generics.R

@@ -1058,11 +1058,17 @@ setMethod("t", signature(x="SpatVector"),
 
 setMethod("terrain", signature(x="SpatRaster"),
 	function(x, v="slope", neighbors=8, unit="degrees", filename="", ...) {
-		unit <- match.arg(unit, c("degrees", "radians"))
-		opt <- spatOptions(filename, ...)
-		seed <- ifelse("flowdir" %in% v, .seed(), 0)
-		x@ptr <- x@ptr$terrain(v, neighbors[1], unit=="degrees", seed, opt)
-		messages(x, "terrain")
+
+	  if (v=="flowacc") {
+	       flowAccu2(x,filename=filename,...)
+
+	  } else {
+		  unit <- match.arg(unit, c("degrees", "radians"))
+		  opt <- spatOptions(filename, ...)
+		  seed <- ifelse("flowdir" %in% v, .seed(), 0)
+		  x@ptr <- x@ptr$terrain(v, neighbors[1], unit=="degrees", seed, opt)
+		  messages(x, "terrain")
+	  }
 	}
 )
 
@@ -1239,3 +1245,25 @@ setMethod("sort", signature(x="data.frame"),
 		x[i, ]
 	}
 )
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

R/generics_watershed2_extension.R

@@ -0,0 +1,89 @@
+# Author: Emanuele Cordano
+# Date : October 2023
+# Version 1.0
+# License GPL v3
+
+setMethod("watershed2", signature(p="SpatRaster",pp_offset="integer"), 
+          function(p,pp_offset,filename="", ...) { 
+
+
+            opt <- spatOptions(filename, ...)
+           ## p@ptr <- p@ptr$watershed2(as.integer(pp_offset-1),opt)
+            p@ptr <- p@ptr$watershed2(as.integer(pp_offset-1),opt)
+            messages(p, "watershed2") ## EC 20210318
+            return(p)
+            ## p@ptr <- uu
+            ##messages(p, "watershed2")
+          }
+
+)
+
+#### END EC 20210702
+
+#### EC 20220809
+setMethod("pitfinder2", signature(p="SpatRaster"), 
+          function(p,filename="", ...) { 
+
+            opt <- spatOptions(filename, ...)
+           # p@ptr <- p@ptr$pitfinder2(opt)# EC 20231026
+            p@ptr <- p@ptr$pitfinder2(opt)
+            messages(p, "pitfinder2") ## EC 20210318
+            return(p)
+            ## p@ptr <- uu
+            ##messages(p, "watershed2")
+          }
+
+)
+#### END EC 20220809
+
+
+#### EC 20231031
+setMethod("NIDP2", signature(p="SpatRaster"), 
+          function(p,filename="", ...) { 
+
+            opt <- spatOptions(filename, ...)
+
+            p@ptr <- p@ptr$NIDP2(opt)
+            messages(p, "NIDP2") ## EC 20231031
+            return(p)
+
+          }
+
+)
+#### END EC 20231031
+
+#### EC 20231104
+setMethod("flowAccu2", signature(p="SpatRaster"), 
+          function(p,filename="", ...) { 
+
+            opt <- spatOptions(filename, ...)
+
+            p@ptr <- p@ptr$flowAccu2(opt)
+            messages(p, "flowAccu2") ## EC 20231104
+            return(p)
+
+          }
+
+)
+#### END EC 20231104
+
+
+#### EC 2023114
+setMethod("flowAccu2_weight", signature(p="SpatRaster",weight="SpatRaster"), 
+          function(p,weight,filename="", ...) { 
+
+            opt <- terra:::spatOptions(filename,...)
+            p@ptr <- p@ptr$flowAccu2_weight(weight@ptr,opt)
+            messages(p, "flowAccu2_weight") ## EC 20231104
+            return(p)
+
+          }
+
+)
+#### END EC 20231104
+
+
+
+
+
+

inst/tinytest/test_flowacc2.R

@@ -0,0 +1,34 @@
+###
+library(terra)
+#> terra 1.7.29
+
+elev <- c(78, 72, 69, 71, 58, 49,
+          74, 67, 56, 49, 46, 50,
+          69, 53, 44, 37, 38, 48,
+          64, 58, 55, 22, 31, 24,
+          68, 61, 47, 21, 16, 19,
+          74, 53, 34, 12, 11, 12) |> matrix(ncol = 6, byrow = TRUE) |> rast()
+
+
+flowdir0 <- c(002,002,002,004,004,008,
+             002,002,002,004,004,008,
+             001,001,002,004,008,004,
+             128,128,001,002,004,008,
+             002,002,001,004,004,004,
+             001,001,001,001,000,016) |> matrix(ncol = 6, byrow = TRUE) |> rast()
+
+flowacc0 <- c(001,001,001,001,001,001,
+              001,002,002,003,003,001,
+              001,004,008,006,005,001,
+              001,001,001,021,001,002,
+              001,001,001,002,025,001,
+              001,003,005,008,036,002) |> matrix(ncol = 6, byrow = TRUE) |> rast()
+
+
+
+flowdir1 <- terrain(elev,"flowdir")
+flowacc1 <- flowAccu2(flowdir1)
+
+result <- (flowacc1==flowacc0) & (flowdir1==flowdir0)
+
+expect_equal(all(result[]),TRUE)

inst/tinytest/test_pitfinder2.R

@@ -0,0 +1,53 @@
+
+## Creation of a Digital Elevation Model 
+
+elev <- array(NA,c(9,9))
+dx <- 1
+dy <- 1 
+for (r in 1:nrow(elev)) {
+  x <- (r-5)*dx
+  for (c in 1:ncol(elev)) {
+
+    y <- (c-5)*dy
+    elev[r,c] <- 10+5*(x^2+y^2)
+  }
+} 
+
+elev <- cbind(elev,elev) 
+elev <- rbind(elev,elev) 
+elev <- rast(elev)
+
+
+
+
+
+## Flow Directions
+
+flowdir<- terrain(elev,v="flowdir")
+
+
+#t(array(flowdir[],rev(dim(flowdir)[1:2])))
+
+## Pit Detect
+
+pits1 <- pitfinder2(flowdir)
+
+xypit <- as.data.frame(pits1,xy=TRUE)
+names(xypit) <- c("x","y","pit")
+xypit$icell <- 1:nrow(xypit)
+xypit[which(xypit$pit!=0),]
+
+
+xypit2 <-  xypit[which(xypit$pit!=0),]
+# > xypit2
+# x    y pit icell
+# 77   4.5 13.5   1    77
+# 86  13.5 13.5   2    86
+# 239  4.5  4.5   3   239
+# 248 13.5  4.5   4   248
+pits0 <- elev*0
+pits0[c(77,86,239,248)] <- c(1,2,3,4)
+
+result <- (pits1==pits0) 
+
+expect_equal(all(result[]),TRUE)

man/NIPD2.Rd

@@ -0,0 +1,101 @@
+\docType{methods}
+\name{NIDP2}
+
+\alias{NIDP2}
+\alias{NIDP2,SpatRaster-method}
+
+
+\title{Number of immediate adjacent cells flowing into each cell (Watershed Analysis)}
+
+\description{
+It computes the number of immediate adjacent cells flowing into each cell (Watershed Analysis)
+}
+
+\usage{
+\S4method{NIDP2}{SpatRaster}(p,filename="",...) 
+}
+
+\arguments{
+  \item{p}{Flow Direction raster map, a \code{\link{SpatRaster-class}} object, see also \code{\link{terrain}}. }
+  \item{filename,...}{list. Options for writing files as in \code{\link{writeRaster}}}
+}
+
+\value{ 
+
+A \code{\link{SpatRaster-class}} (raster) map containing the value for each cell. 
+}
+
+\details{
+NDIP is computed to preliminarly to compute Flow Accumulation area in the algorithm proposed by Zhouo at al,2019. 
+
+
+}
+
+\seealso{\code{\link{flowAccu2}},\code{\link{flowAccu2_weight}}}
+\author{
+
+Emanuele Cordano (implementation of the algorithm within the package),Guiyun Zhou, Wenyan Dong and Hongqiang Wei
+}
+
+\examples{
+
+
+elev1 <- array(NA,c(9,9))
+elev2 <- elev1
+dx <- 1
+dy <- 1 
+for (r in 1:nrow(elev1)) {
+  y <- (r-5)*dx
+  for (c in 1:ncol(elev1)) {
+
+    x <- (c-5)*dy
+    elev1[r,c] <- 5*(x^2+y^2)
+    elev2[r,c] <- 10+5*(abs(x))-0.001*y ### 5*(x^2+y^2)
+  }
+} 
+
+
+## Elevation Raster Maps
+elev1 <- rast(elev1)
+elev2 <- rast(elev2)
+
+t(array(elev1[],rev(dim(elev1)[1:2])))
+t(array(elev2[],rev(dim(elev2)[1:2])))
+
+plot(elev1)
+plot(elev2)
+
+## Flow Direction Raster Maps
+flowdir1<- terrain(elev1,v="flowdir")
+flowdir2<- terrain(elev2,v="flowdir")
+
+
+t(array(flowdir1[],rev(dim(flowdir1)[1:2])))
+t(array(flowdir2[],rev(dim(flowdir2)[1:2])))
+
+plot(flowdir1)
+plot(flowdir2)
+
+## 
+nidp1 <- NIDP2((flowdir1))
+nidp2 <- NIDP2((flowdir2))
+
+t(array(nidp1[],rev(dim(nidp1)[1:2])))
+t(array(nidp2[],rev(dim(nidp2)[1:2])))
+
+plot(nidp1)
+plot(nidp2)
+
+}
+
+
+
+
+\references{
+Zhou, G., Wei, H. & Fu, S. A fast and simple algorithm for calculating flow accumulation matrices from raster digital elevation. Front. Earth Sci. 13, 317–326 (2019). https://doi.org/10.1007/s11707-018-0725-9
+\url{https://link.springer.com/article/10.1007/s11707-018-0725-9}
+}
+
+
+
+\keyword{spatial}