# Flow accumulation on grid in R

I would like to calculate the number of cells from gridded elevation data that drain into a particular cell. First, I calculate the flow direction by selecting the cell with the steepest slope from one out of eight adjacent cells.

``````library(raster)
library(gdistance) # package for adjacent function

# make the blank raster, the extents and resolution of which are up to you.
# the cells are numbered 1 to ncell (starting in the top left in R)
r <- raster(xmn=0,xmx=6,ymn=0,ymx=6,res=c(2,2))
r[] <- 1:ncell(r)

# flow direction of digital elevation model
flow <- terrain(subelevation, opt="flowdir", neighbors=8)
# transform to matrix
dim(flow@data@values) <- c(sqrt(length(flow@data@values)), sqrt(length(flow@data@values)))
# yields values from 1 to 128, see Figure 1
flow@data@values
``````

Figure 1

Figure 2 plots `flow@data@values` transformed to an `igraph` object to illustrate the problem. For any chosen cell, I would like to calculate the number of incoming cells connected at any distance. For example if the cell labelled A is chosen, it should have the value 3 as B, C and D drain into it. The more efficient the solution, the better as the gridded dataset is large.

Figure 2

• Why not do this in a GIS that has more robust algorithms than D8? If you would like to stay in R you can use the package interfaces to GRASS (`rgrass7`), SAGA (`RSAGA`), or indirectly QGIS (`RQGIS`). Computationally, this will be much more efficient as well as providing flow direction algorithms such as, D-infinity. Commented Mar 26, 2020 at 19:21
• I have been mostly using `R` which is why I would prefer this program. There must be a neat way of doing this. `igraph` allows finding clusters of nodes. Maybe, I can achieve this by identifying clusters of nodes. Commented Mar 27, 2020 at 0:08

`whitebox::wbt_d8_flow_accumulation()` calculates a flow accumulation grid.

I was trying to find a method for doing this in R, for stream network and hydrology analysis. It was surprisingly hard to find an answer, with most suggestions involving bridging to another GIS software which I wasn't overly excited about, but I have found a very simple way!

It uses the R package `whitebox`, which admittedly is an R-frontend for Whitebox Tools software; see here for details on the package, and here for info on Whitebox Tools software and detailed functionalities. Don't be put off by the fact it isn't 100% R (I didn't even realise this until I started looking at other functions), because there is an almost endless list of awesome functions inside it.

I discovered this package through JP Gannon's "Hydroinfomatics at VT" bookdown, which has some great examples of stream network analysis. Here, I have basically copied from those examples (for the sake of providing an example rather than just the link), but modified to include different data. Specifically, I am taking the workflow from Chapters 15 and 16 - do read these as it explains the process in better detail than the following minimal example.

First up, load the packages we will be using:

``````library(terra)
library(tmap)
library(tidyverse)

#install.packages("whitebox", repos="http://R-Forge.R-project.org")
library(whitebox)
``````

For this example we use the `srtm` elevation raster of Zion National Park, from the `spDataLarge` package. We attach it and then save it locally in `my_folder` (this is important), and then import and plot it.

``````# load elevation data from spDataLarge package (elevation data for Zion National Park)
# save sample data locally
system.file("raster/srtm.tif", package = "spDataLarge") %>%
rast %>%
writeRaster('my_folder/my_rast.tif', overwrite = TRUE)

# import dem
rast('my_folder/my_rast.tif') %>%
{. ->> my_rast}

tmap_mode('view')

tm_shape(my_rast)+
tm_raster(style = 'cont', palette = 'viridis')
``````

Now we crack into `whitebox`. Without repeating the original example, basically `whitebox` functions look a little different to what you're probably used to - instead of R objects as inputs and outputs, you work with filepaths (this is why we saved our raster locally in our own folder). I assume this has something to do with `whitebox` package being a frontend for the standalone Whitebox Tools software.

The first thing we do is to prepare the raster for analysis. This involves 'filling in' single-cell holes/pits/depressions as this can cause drainage paths to terminate unexpectedly. Likewise, we 'breach' larger depressions (see the original examples for more details). You can see here we supply the functions with filepaths for (input) `dem` and `output` objects.

``````# fill and breach
wbt_fill_single_cell_pits(
dem = 'my_folder/my_rast.tif',
output = 'my_folder/my_rast_filled.tif'
)

wbt_breach_depressions_least_cost(
dem = 'my_folder/my_rast_filled.tif',
output = 'my_folder/my_rast_filled_breached.tif',
dist = 5
)
``````

Next step we calculate a (D8) flow accumulation grid - cell values in the output represent the number of cells that flow into it. This can be changed to catchment area or specific catchment area - see the original example.

This is essentially the answer to the original question, so once it is created we also import it into R.

``````# D8 flow accumulation
wbt_d8_flow_accumulation(
input = 'my_folder/my_rast_filled_breached.tif',
output = 'my_folder/my_rast_d8fa.tif',
out_type = 'cells'
)

# now, we can import the flow accumulation raster
rast('my_folder/my_rast_d8fa.tif') %>%
{. ->> my_rast_d8fa}

my_rast_d8fa

# class       : SpatRaster
# dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
# resolution  : 0.0008333333, 0.0008333333  (x, y)
# extent      : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs
# source      : my_rast_d8fa.tif
# name        : my_rast_d8fa
``````

And, we plot it:

``````tm_shape(my_rast_d8fa)+
tm_raster(palette = 'viridis')
``````

Due to massive variation in values (some cells have 10-100 cells draining into them, others at the bottom of the landscape have 100,000-200,000), this is hard to visualise. Rather, we plot log-transformed values:

``````tm_shape(my_rast_d8fa %>% log10)+
tm_raster(palette = 'viridis')
``````

Unsurprisingly, this looks like a stream network!

So, once you've got the raster in R, you can pull out the values for each cell whichever way suits:

``````my_rast_d8fa %>%
raster::as.data.frame(xy = TRUE) %>%
as_tibble %>%
rename(
n_cells = my_rast_d8fa
)
# # A tibble: 212,404 x 3
#        x     y n_cells
#    <dbl> <dbl>   <dbl>
# 1  -113.  37.5       1
# 2  -113.  37.5       3
# 3  -113.  37.5       4
# 4  -113.  37.5       5
# 5  -113.  37.5       6
# 6  -113.  37.5       1
# 7  -113.  37.5       1
# 8  -113.  37.5       1
# 9  -113.  37.5       1
# 10 -113.  37.5       2
# # ... with 212,394 more rows

my_rast_d8fa %>%
terra::values() %>%
#      my_rast_d8fa
# [1,]            1
# [2,]            3
# [3,]            4
# [4,]            5
# [5,]            6
# [6,]            1
``````

Bonus - If you're doing stream analysis, you can also extract streams from the flow accumulation raster. You just need to provide a threshold for flow accumulation values to distinguish streams from non-streams:

``````# extract streams
wbt_extract_streams(
flow_accum = 'my_folder/my_rast_d8fa.tif',
output = 'my_folder/raster_streams.tif',
threshold = 5000
)

# import and inspect streams raster
rast('my_folder/raster_streams.tif') %>%
{. ->> my_rast_streams}

# plot streams raster
tm_shape(my_rast)+
tm_raster(style = 'cont', palette = 'viridis')+
tm_shape(my_rast_streams)+
tm_raster(palette = c('white'), legend.show = FALSE)
``````

#### Notes:

For me, the `wbt_fill_*` and `wbt_breach_*` functions didn't cooperate when I had the `dem` in the working directory, but was fine when I put `dem` in a subfolder (`my_folder`).

D8 and D-Infinity flow accumulation grids can both be created - the example explains the difference but which one you want will depend on your exact situation I guess. D-Infinity looks like it would generally be a more 'realistic' representation of overland water flow, but I don't know enough about this type of analysis to say which one is more appropriate (from what I have read there is no right/wrong).

• I have also had issues when there is whitespace in the file path or file name. Replacing whitespace with underscores has helped me a couple of times. Commented Oct 11, 2021 at 23:36