In jmigueldelgado/a5udes: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(sf,warn.conflicts=FALSE)
library(a5udes,warn.conflicts=FALSE)
library(dplyr,warn.conflicts=FALSE)
library(igraph,warn.conflicts=FALSE)
library(knitr,warn.conflicts=FALSE)
library(ggplot2,warn.conflicts=FALSE)

Global datasets of rivers and reservoirs

knitr::include_graphics('/home/delgado/proj/a5udes/screenshot.png')

Rivers

We use HydroSHEDS as a consistent global river dataset. The resolution is about 90 m around the equator (3 arc-seconds).

Lehner, B., Verdin, K., Jarvis, A. (2008): New global hydrography derived from spaceborne elevation data. Eos, Transactions, AGU, 89(10): 93-94.

For example in Ceará:

plot(river_geometry["UP_CELLS"])
# ggplot(river_geometry) + geom_sf(aes(color=UP_CELLS))

Disjoint components of the graph

Since the river network dataset of HydroSHEDS has valid topology, we can operate on it with package igraph, using graph theory. We can for example obtain all disjoint components of the graph:

riv_split = split_river_network(river_geometry)
plot(riv_split["membership"])

split_river_network

Selecting a river system by a reach id, eg in the Jaguaribe river basin:

reach_id=140877
riv_jagua = select_disjoint_river(reach_id,riv_split)
plot(riv_jagua["UP_CELLS"])

Reservoirs

We use the JRC global surface water dataset.

Jean-Francois Pekel, Andrew Cottam, Noel Gorelick, Alan S. Belward, High-resolution mapping of global surface water and its long-term changes. Nature 540, 418-422 (2016). (doi:10.1038/nature20584)

Resolution is 30 m. We run

allocate_reservoir_to_river(riv_jagua), with e.g. riv_jagua obtained in the previous step; and
build_reservoir_topology(reservoir_geometry_subset,riv_jagua,riv_graph).

in advance in order to create a graph. Then we can obtain any subgraph of choice.

head(reservoir_graph)

If you wish to pass attributes to the reservoir graph using tidygraph makes it very easy. For example passing the attribute area_max to the graph:

library(tidygraph)
res_area=reservoir_geometry %>%
  st_set_geometry(NULL) %>%
  select(id_jrc,area_max) %>%
  rename(name=id_jrc)

reservoir_tidygraph = reservoir_tidygraph %>%
  activate(nodes) %>%
  mutate(name=as.integer(name)) %>%
  left_join(res_area)

Small example with reservoir 34560

Example 1: Identify adjacent reservoirs upstream of 34560:

id='34560'
adj_subgraph=make_ego_graph(reservoir_graph, order = 1, nodes = id, mode = "in", mindist = 0)
plot(adj_subgraph[[1]])

...and plot them on a georreferenced map:

adj=neighbors(reservoir_graph,id,"in")
res_adj=filter(reservoir_geometry,id_jrc %in% adj$name)
res_id=filter(reservoir_geometry,id_jrc %in% id)



ind=st_combine(res_adj) %>% st_convex_hull %>% st_intersects(river_geometry)
river_subset=river_geometry[ind[[1]],]

ggplot() +
  geom_sf(data=res_adj,color='darkblue',fill='darkblue') +
  geom_sf(data=res_id,color='red',fill='red') +
  geom_sf(data=river_subset,color='black',size=0.1)

Example 2: Identify all reservoirs upstream of 34560:

sub=all_simple_paths(reservoir_graph,from=id,mode='in') %>%
  unlist %>%
  unique

Vsubgraph = induced_subgraph(reservoir_graph, sub ,impl='auto')

plot(Vsubgraph)

...and plot them on a georreferenced map:

vertices = V(Vsubgraph)  

res_subset = filter(reservoir_geometry,reservoir_geometry$id_jrc %in% vertices$name)

reach_id = filter(res_subset,`id_jrc`==id) %>%
  pull(`nearest river`)
riv_upstr=river_upstream(reach_id,river_geometry,river_graph)

res_subset %>%
  mutate(`nearest river`=as.factor(`nearest river`)) %>%
  ggplot(.) +
    geom_sf(aes(color=`nearest river`)) +
    geom_sf(data=riv_upstr)