feature request: `terrain(x, v = "flowacc")`

[dimfalk]

Hi Robert,

as already briefly mentioned in #1180, I'd very much appreciate the possibility to estimate flow accumulation using {terra} based on flow direction estimates coming from terrain(x, v = "flowdir").

Using another illustrated example published in ESRI's online help for reference, I thought I'd at least be able to provide a working implementation in R, which could have been ported to C++ but I have to admit that I was not able to resolve some minor flaws in the end (and currently I don't really have the time to investigate further) to fully reproduce the exampe given - but you'll probably solve this in a more elegant way anyway, I guess... :smirk:

However, at least I'll leave this elevation raster from upper reference for convenience if you wanted to test whether the algorithm is working as expected as input for terrain():

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)

library(terra)
#> terra 1.7.29

dem <- rast(elev)

Thank you very much in advance!

[kadyb]

Wouldn't it be easier to just use the existing tools in SAGA or GRASS from R?

[dimfalk]

Thanks for pointing to these resources!

I had something similar in mind at second thought but since Robert seemed to be interested in adding flow accumulation to terra, I wanted to leave this little reminder...

However, I feel like I should have a look at rgrass nonetheless.

[ecor]

Hi all, I inserted a flow accumulation method, called flowAccu2 within my fork (ecor/terra) of this repository (rspatial/terra). This is experimental but it looks working. Please feel free to test it on your cases and do not hesitatate to send me your feedback. Thank you best Emanuele @ecor

> remotes::install_github("ecor/terra")
> library(terra)
terra 1.7.62
> help(flowAccu2) ## See/run examples

[chris-english]

I wrote a naive flow direction lines function that returns a flow direction lines whose intermediate products to creating the geom_mtx can be used to 'rank' accumulation and write new dem values therefrom. It works on terra 'elev', 'volcano', but does not scale usefully thereafter. naive_descent, just has .txt file with function and command line example. Look forward to testing flowAccu2.

I think that some users might expect the vect product as part of flow accumulation.

[chris-english]

@ecor your NIDP2 and FlowAccu2 work in all cases I tested. Look forward to testing more.

[ecor]

@chris-english , Thank you

[ecor]

Hi @dimfalk and others,

here is a snippet of code :

###


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)

library(terra)  ## remotes::install_github("ecor/terra")


dem <- rast(elev)

flowdir <- terrain(dem,"flowdir")
flowacc <- flowAccu2(flowdir)

Any feedback is appreciated. I would like to create a pull request onto the main fork/branch in the future. Thank you Best @ecor

[dimfalk]

Hi @ecor,

sorry for the late reply, I have been a little bit busy lately, unfortuantely... Anyway, your implementation looks quite promising - thank you very much! - but it seems like your resulting SpatRaster object, from a numerical perspective, is by one larger compared to the example in ESRI's online help, c.f. my opening post?

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)

# devtools::install_github("ecor/terra")
library("terra")
#> terra 1.7.73

dem <- rast(elev)

flowdir <- terrain(dem, "flowdir")

flowacc <- flowAccu2(flowdir)
flowacc
#> class       : SpatRaster 
#> dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 6, 0, 6  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : flowdir 
#> min value   :       1 
#> max value   :      36

plot(flowacc)
text(flowacc)

To put this a little naively, it looks like your "background value" is not correct here - at least in comparison to the reference. Based on the definiton used ("The result of Flow Accumulation is a raster of accumulated flow to each cell [...]"), the first row and column of the resulting flow accumulation raster based on flowdir input should indeed be zero, from my point of view, since no other cells are draining into these cells. The flow accumulation at the outlet should be ncell(dem) - 1 in this case, right?

Of course this can be fixed quite pragmatically by considering flowacc - 1 to be the result, but maybe this can be fixed natively?

Or am I just falling for some ESRI-specific implementations here and this interpretation/approach is not common in other GIS environments like SAGA & GRASS?

[ecor]

Thanks for the feedback, @dimfalk and sorry too for my delay. I'm not so familiar with these operations in QGIS or GRASS. I did this kind of job with jgrasstools / the Horton Machine many years ago. In my personal opinion, the pixel should belong to the upstream watershed and so it must be considered, so that the area should be always >=1. Anyway, I will check in QGIS and GRASS and related tools in order to standardize the operation. I think both options might be taken into account, they differ only by 1 . Thank you for the suggestion.

[dimfalk]

No worries! :smirk:

My personal and careful interpretation (based on ESRI's deinition, c.f. above) would be that these grid cells definitely belong to the watershed, yes, and are also generally considered in watershed deliniation, but in this specific case, where we are talking about flow accumulation only, they are still 0 - which does not mean they are not considered, from my point of view. There are just no additional cells further upstream (because the region of interest is limited to a subwatershed, or we are effectively located at the watershed boundary) draining into these cells. But I guess this is just a matter of definition and the underlying algorithms used.

However, yes, let's have a look how this is implemented in other software using the elev raster definied initially. I have only worked with ArcGIS so far when it comes to watershed operations, to be honest, so I also have to do some digging first.

I'm not yet familiar with {rgrass} and {RSAGA}, so I apologise for the following information not being reproducible:

• Using r.flow from GRASS 8.3.1 the cells at the upper and left border were indeed also 0. The max value at the outlet was deviating from ESRI's approach and also the spatial pattern was not completely identical, but I assume this might be the result of using a different flow direction algorithm instead of D8?

_Supplement material confirms this_: r.flow uses an original vector-grid algorithm which uses an infinite number of directions between 0.0000... and 360.0000... and traces the flow as a line (vector) in the direction of gradient (rather than from cell to cell in one of the 8 directions = D-infinity algorithm). They are traced in any direction using aspect (so there is no limitation to 8 directions here). The D8 algorithm produces zig-zag lines. The value in the outlet is very similar for r.flow algorithm (because it is essentially the watershed area), however the spatial distribution of flow, especially on hillslopes is quite different. It is still a 1D flow routing so the dispersal flow is not accurately described, but still better than D8.

https://grass.osgeo.org/grass83/manuals/r.flow.html

• saga_cmd ta_hydrology 29 from SAGA 9.2.0 seems to behave differently, at least per default. Here the grid cells of interest are 1, but I'm missing additional information in the docs and currently don't have the time to dig into the referenced mentioned, unfortunately.

https://saga-gis.sourceforge.io/saga_tool_doc/9.2.0/ta_hydrology_29.html

[ecor]

Thanks @dimfalk for the references, if the pour point cells must not be included, one can use flowAccu2 using weight arguments a raster map with values 1 for each cell except the pour point which is 0. Otherwise , putting the result in the R session memory, one can operates with -1 at the end. Otherwise I have to modify the internal C/C++ code a little bit. If @rhijmans and the other authors agree , I can give this contributions through a pull request from the @ecor fork to the @rspatial one. Thank you

[ecor]

I just created a pull request : https://github.com/rspatial/terra/pull/1477