Flow accumulation on grid in R

Question

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

Answer 1

Short answer:

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

Long answer:

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() %>% 
  head
#      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).