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).
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.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.