1.1 Packages needed

This guide relies on several packages for data processing (the links direct to useful resources and tutorials):

• sf, lwgeom, rgdal, raster, elevatr, RANN packages for obtaining and processing spatial data;
• igraph package for performing network analyses;
• tidyverse packages for efficient and easy-to-read data pipeline (see here and here for the most frequently used functions);
• ggnetwork for producing easy, ggplot-like, plot of class igraph objects and ggspatial for easy, ggplot-like plot, of spatial data;
• corrmorant to plot nice correlation plots;
  
• riverconn to calculate the fragmentation indices for the riverscape.

Make sure those libraries are installed, updated, and loaded.

library(tidyverse)
library(sf)
library(raster)
library(ggspatial)
library(viridis)
library(igraph)
library(riverconn)
library(elevatr)
library(gridExtra)
library(ggnetwork)
library(lwgeom)
library(gridExtra)
library(corrmorant)
library(RANN)
library(ggpubr)
library(cowplot)

Most of the geospatial processing in this tutorial occurs on sf objects. Thanks to their data.frame-like structure, spatial operations (e.g., selecting and updating columns and rows, join with other data frames) are “easier” on such objects. Most spatial objects are converted to this class as soon as possible in the script.

1.2 Workflow steps

The workflow suggested in this document follows several steps:

• data retrieval from existing datasets;
  
• data pre-processing;
  
• network creation;
  
• indices calculation and comparison.

1.3 Sources of the spatial data

The construction of the igraph object requires at least some information on the spatial geometry of the river system and the position of the barriers.

The river network shapefile can be retrieved from HydroSHEDS, a global hydrographic atlas with geo-referenced datasets derived from spaceborn elevation data. HydroSHEDS is widely used for global to regional analyses (Lehner et al., 2008). River network data (HydroRIVERS) can be downloaded from here. Basins and catchments shapefiles (HydroBASINS) are also available with 12 detail levels (Lehner et al., 2013). The Ebro catchment border can be selected from the polygons available in the HydroBASINS layer (use level 5 basins shapefile for the identification) and the corresponding network can be extracted from the HydroSHEDS layer.

Other options are available for catchments and rivers. The European Catchments and Rivers network system (ECRINS) dataset provides good quality data covering most of Europe. The ECRINS data for rivers are structured with segments and nodes, making the creation of a network quite straight-forward. ECRINS has a higher resolution than HydroRIVERS. The recent Hydrography90m dataset can also be used. Local sources (e.g. hydrography layers from local authorities) are also generally available, but particular care should be taken in preprocessing to correct the presence of disconnected segments and loops in the shapefile. Different HydroSHEDS and HydroRIVERS versions are available. This tutorial uses HydroRIVERS v1.1: note that the previous versions were slighthly different in terms of attributes (for further details see the HydroRIVERS technical documentation).

The shapefile with the geographic position of the dams was obtained from the Ebro river public Geographic Information System CHE website.

The AMBER dataset described in Belletti et al. (2020) contains barriers data compiled from different national and international sources. Metadata on barriers type (weir, dam, culvert, etc) and height are available but most barriers are unclassified.

The corresponding shapefiles are available here and can be imported in R (note that the river shapefile “Ebre_riv_modif.shp” has been sligtly modified from the original hydrosheds shapefile -see next sections on the checks carried out on the shapefiles. The original shapefile is saved in the Ebro_shape_original/ folder for comparison).

shape_river <- st_read("Ebro_shape_corrected/Ebro_rivers.shp")
## Reading layer `Ebro_rivers' from data source 
##   `C:\Users\Damiano\Dropbox\riverconn_tutorial\Ebro_shape_corrected\Ebro_rivers.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4985 features and 14 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -4.3625 ymin: 40.3625 xmax: 2.116667 ymax: 43.15417
## Geodetic CRS:  WGS 84
shape_basin <- st_read("Ebro_catchment/Ebro_mask.shp")
## Reading layer `Ebro_mask' from data source 
##   `C:\Users\Damiano\Dropbox\riverconn_tutorial\Ebro_catchment\Ebro_mask.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4.4 ymin: 40.34167 xmax: 2.154167 ymax: 43.17917
## Geodetic CRS:  WGS 84
shape_dams <- st_read("Ebro_dams/Embalses/Embalses.shp")
## Reading layer `Embalses' from data source 
##   `C:\Users\Damiano\Dropbox\riverconn_tutorial\Ebro_dams\Embalses\Embalses.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 224 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 408868.8 ymin: 4505766 xmax: 903207.1 ymax: 4766629
## Projected CRS: ETRS89 / UTM zone 30N

1.4 Functions to support data processing

Some custom functions were developed and were used in the data processing. The source code and the documentation for these functions is found in Appendix B. Before starting, make sure those functions are loaded in the global environment: get_outlet, create_network, snap_to_river, headwaters_dam, multiple_confluences, gaps_detection, check_components. Detailed information on the use of these functions is presented in the next sections.

2.1 Shapefiles processing

Shapefiles derived from the HydroRIVERS layer have a high detail level: a network with a high detail level will have a large amount of nodes and links, making the algorithms for network analysis slow. Thus it makes sense to reduce the size of the network by pruning the most upstream reaches. The information on the upstream area (suquare kilometers) is stored in the field ‘UPLAND_SKM’.

Note that different HydroRIVERS versions are available. This tutorial uses HydroRIVERS v1.1, wuere different fields are available for pruning the river network based on the size (e.g. ORD_STRA, ORD_CLAS, ORD_FLOW). Previous versions have different attributes. For instance, in HydroRIVERS v1.0 the field UP_CELLS (number of upstream cells) can be used to prune the river network. Check the dataset documentation for further details.

# Set a threshold of 80 square kilometers
threshold = 80

# Prune HydroRIVERS network based on upstream area
shape_river_small <- shape_river[as.numeric(shape_river$UPLAND_SKM) > threshold,]

Next, the dams shapefile can be converted to a shapefile and a unique id can be given to each point. As the dams dataset is a polygon shapefile the centroid of each poligon has to be calculated to obtain a “point” layer. The reference system of the dams layer is changed to WGS84, the same reference system of HydroSHEDS and HydroRIVERS.

dams_to_points <- shape_dams %>%
  st_as_sf %>%
  st_centroid %>%
  mutate(id = 1:nrow(.))%>% 
  dplyr::select(id)%>% 
  st_transform(crs = "+proj=longlat +datum=WGS84")

The shapefiles can be plotted to check there are no errors in the geolocalization of the data.

ggplot() +
  coord_fixed() +
  theme_minimal() +
  ggspatial::layer_spatial(shape_basin, fill = NA, color = "gray90") +
  ggspatial::layer_spatial(shape_river_small, aes(color = log10(UPLAND_SKM)) )+
  ggspatial::layer_spatial(dams_to_points) +
  scale_color_viridis(direction = -1, name= "upstream area \n(log10[Km^2])") +
  theme(legend.position = "bottom") +
  ggspatial::annotation_scale(location = "bl", style = "ticks") +
  ggspatial::annotation_north_arrow(location = "br")+
  labs(caption = "Black dots are the position of the dams")

2.2 Confluences processing

Every river network contains confluences, i.e. points where two reaches merge into one. In the river graph defined above, confluences are edges (links) between nodes. To define confluences the pruned river shapefile is first simplified and its geometry is casted to “POINT”. Note that for large shapefiles the function rgeos::gLineMerge might be more efficient than sf::st_union.

# Simplify river shapefile
shape_river_simple <- shape_river_small %>%
  st_as_sf %>%
  st_union()

# Convert shapefile to point object
river_to_points <- shape_river_simple %>%
  st_as_sf %>%
  st_cast("POINT") %>%
  mutate(id = 1:nrow(.))

Since points are created at the extremes of every segment in the original shapefile, a confluence is delineated where three (or more) of such nodes are overlapping.

# Check points that overlap
joins_selection <- river_to_points %>%
  st_equals() %>%
  Filter(function(x){length(x) > 2}, .) %>%
  lapply(., FUN = min) %>%
  unlist() %>%
  unique()

# Filter original point shapefile to retain only confluences
river_joins <- river_to_points %>% 
  filter(id %in% joins_selection)

The last step is the splitting of the simplified shapefile where it intersects a confluence. This allows to identify unequivocally the river sections that are comprised between two confluences. Such elements will then became vertices (nodes) of the graph. The function lwgeom::st_split is used. This function splits an object of class ‘sf’ with ‘MULTILINESTRING’ geometry into and object of class ‘sf’ containing multiple ‘LINESTRING’ geometries records based on the position of the points. Note that for this function to work properly, the point object must perfectly overlap with the polyline (this can be checked with the function sf::st_distance applied to the points and the polyline, which must return a vector of zeroes).

# Split polyline
shape_river_simplified <- lwgeom::st_split(shape_river_simple, river_joins) %>%
  st_collection_extract(.,"LINESTRING") %>%
  data.frame(id = 1:length(.), geometry = .) %>%
  st_as_sf() %>%
  mutate(length = st_length(.))

Results can be checked with a plot.

ggplot() +
  coord_fixed() +
  ggspatial::layer_spatial(shape_river_simplified, aes(color = id))+
  scale_color_viridis(direction = -1, name= "Reach ID") +
  ggspatial::layer_spatial(river_joins, shape = 1)+
  theme_minimal() +
  theme(legend.position = "bottom")+
  ggspatial::annotation_scale(location = "bl", style = "ticks") +
  ggspatial::annotation_north_arrow(location = "br")+
  labs(caption = "Hollow points are the position of the junctions")

2.3 Dams processing

Dams are processed in a way similar to confluences. The general idea is to obtain a shapefile that is sliced in points where dams or confluences are located. However, dams shapefiles require special care in the processing of the data. First, dams need to be snapped to the river network. Second, dams that are too far away from the river network need to be removed from the analysis because they are probably located in smaller streams that were removed beforehand.

The function snap_to_river (see Appendix B) is used to snap the dams points to the river network based on a tolerance threshold (here, 1000 m were used). After this function, the distance of snapped dams from the river network can be checked to make sure the algorithm worked properly.

# Snap dams
dams_snapped <- snap_to_river(dams_to_points,
                              shape_river_simple %>% st_sf(),
                              max_dist = 1000)
# Retain dams that were snapped
dams_snapped_reduced <-
  dams_snapped[st_contains(shape_river_simple %>% st_sf(), dams_snapped, prepared = FALSE)[[1]],]

# Check if dams were snapped properly to the network (all the distances should be zero)
st_distance(dams_snapped_reduced, shape_river_simple) %>% sum
## 0 [m]

The final step is to check if two dams were snapped in the same geographic point, and to assign passablity values to the obtained snapped points. Here, uniform (not dam-dependent) passabilities were assigned: 0.8 for downstream passability, and 0.1 for upstream passability. An id is assigned to each of the geographic points. Finally, each dam is assigned the value of the upstream area with a spatial join with the original river shapefile (HydroSheds expresses the upstream area as number of cells).

dams_snapped_reduced_joined <- dams_snapped_reduced %>%
  mutate(cluster =
           st_equals(dams_snapped_reduced, dams_snapped_reduced) %>%
           sapply(., FUN = function(x) paste0(x, collapse = "_"))) %>%
  group_by(cluster) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(id_dam = as.character(1:nrow(.)), pass_u = 0.1, pass_d = 0.8) %>%
  as.data.frame %>%
  st_as_sf() %>%
  st_join(., shape_river_small, join = st_is_within_distance, dist = 10 ) %>% 
  group_by(id) %>%
  slice(which.max(UPLAND_SKM)) %>% 
  ungroup()
  

The function headwaters_dam (see Appendix B) can identify dams that are located in the headwaters, meaning that they do not have any river network segment upstream. When no upstream segment is present, the network can not be created (because the dam would be an edge without an upstream vertex). When calling the flag_headwater you should get a column with only “FALSE” statements.

headwaters_checking <- headwaters_dam(dams_snapped_reduced_joined, shape_river_simple)
head(headwaters_checking$flag_headwater)
##       [,1]
## [1,] FALSE
## [2,] FALSE
## [3,] FALSE
## [4,] FALSE
## [5,] FALSE
## [6,] FALSE

Results can be checked with a plot.

ggplot() +
  coord_fixed() +
  ggspatial::layer_spatial(shape_river_simple, color = "gray70")+
  ggspatial::layer_spatial(dams_snapped_reduced, shape = 1) +
  theme_minimal() +
  theme(legend.position = "bottom")+
  ggspatial::annotation_north_arrow(location = "br")+
  ggspatial::annotation_scale(location = "bl", style = "ticks") +
  labs(caption = "Hollow points are the position of the dams")

2.4 Shapefile slicing

The last step is to slice the river shapefile based on both the position of dams and confluences. First a combined ‘junctions’ dataset is created (the network_links object), then the function lwgeom::st_split is applied again. note that the confluences are assigned NA values to the fields present in the dams shapefile. Note also that in the river network shapefile, the field length is calculated. Finally, a spatial join is performed with the original river shapefile to assign to each reach the corresponding upstream area.

The field NodeID is created in the river_net_simplified object. This field will be used later for linking graph properties back to the shapefile for further plotting.

# Create junction point shapefile
network_links <- rbind(
  dams_snapped_reduced_joined %>% 
    mutate(type = "dam", id_barrier = id_dam) %>%
    dplyr::select(type, id_barrier, pass_u, pass_d),
  river_joins %>% mutate(type = "joint") %>%
    dplyr::select(type) %>%
    mutate(id_barrier = NA, pass_u = NA, pass_d = NA) %>%
    rename(geometry = x)) %>%
  mutate(id_links = 1:nrow(.))

# Split river network
river_net_simplified <- lwgeom::st_split(shape_river_simple, network_links) %>%
  st_collection_extract(.,"LINESTRING") %>%
  data.frame(NodeID = 1:length(.), geometry = .) %>%
  st_as_sf() %>%
  mutate(length = st_length(.)) %>%
  st_join(., shape_river_small, join = st_is_within_distance, dist = 0.01 ) %>% 
  group_by(NodeID) %>%
  slice(which.max(UPLAND_SKM)) %>% 
  ungroup()

A plot can reveal if the process was successful.

ggplot() +
  coord_fixed() +
  ggspatial::layer_spatial(river_net_simplified, color = "gray70")+
  ggspatial::layer_spatial(network_links, aes(shape = type))+
  scale_shape_manual(name = "Splitting points", values=c("dam" =17,"joint" = 23))+
  theme_minimal() +
  theme(legend.position = "bottom")+
  ggspatial::annotation_north_arrow(location = "br")+
  ggspatial::annotation_scale(location = "bl", style = "ticks")

Sometimes, the HydroSHEDS-derived shapefiles might have confluences where four or more rivers join. Such type of confluences can cause problems in the next steps. The function multiple_confluences (Appendix B) identifies such problematic confluences so that they can be corrected in a GIS environment. The correction consists in editing the river shapefile to move re-locate one of the reaches (either upstream or downstream) on the intersections that have more than three reaches (this can be easily achieved once these intersections are located by the following function).

Once edited, re-run the tutorial and check again that there are no more four reaches confluences.

confluences <- multiple_confluences(river_net_simplified)
head(confluences)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -4.047917 ymin: 42.78958 xmax: -1.364583 ymax: 42.97292
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 4
##   NodeID             geometry n_confluences flag_confluences
##    <int>          <POINT [°]>         <int> <lgl>           
## 1      1  (-3.46875 42.95208)             3 FALSE           
## 2      2 (-4.047917 42.97292)             3 FALSE           
## 3      5 (-3.552083 42.96875)             3 FALSE           
## 4      7 (-1.364583 42.81042)             3 FALSE           
## 5     10 (-1.835417 42.90625)             3 FALSE           
## 6     11 (-3.352083 42.78958)             3 FALSE

Another common problem in HydroSHEDS-derived shapefiles is the presence of gaps between segments. This issue can cause problems in the network_creation function, which expects a shapefile with no gaps. The functions gaps_detection and check_components (Appendix B) can be used to detect such problematic points in the shapefile. In particular, the function check_components returns an sf file that is a copy of the original river_net_simplified with the additional field component that contains a progressive number identifying multiple components. A correct shapefile should contain only ‘1’ values in this field. A plot can help identifying the location of those components and a manual correction in a gis environment can fix the problem.

Once edited, re-run the tutorial and check again that the shapefile has only one component.

shp_check <- check_components(network_links, river_net_simplified)
head(shp_check)
## Simple feature collection with 6 features and 2 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -4.277083 ymin: 42.93958 xmax: -3.464583 ymax: 43.11458
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 3
##   NodeID                                                        geometry compo…¹
##    <int>                                                <LINESTRING [°]>   <dbl>
## 1      1 (-3.58125 43.11458, -3.572917 43.10625, -3.572917 43.10208, -3…       1
## 2      2 (-4.277083 43.01042, -4.272917 43.01458, -4.264583 43.01458, -…       1
## 3      3 (-3.95625 43.00625, -3.960417 43.00625, -3.964583 43.00208, -3…       1
## 4      4 (-3.927083 43.01458, -3.94375 43.01458, -3.947917 43.01042, -3…       1
## 5      5 (-3.602083 43.01458, -3.597917 43.01042, -3.589583 43.01042, -…       1
## 6      6 (-3.74375 43.00208, -3.739583 42.99792, -3.739583 42.98958, -3…       1
## # … with abbreviated variable name ¹​component

2.5 Adding additional attributes to the river shapefile

Additional attributes useful for the network calculation can be added to the river shapefile. Elevation is an example. Here we use the function elevatr::get_elev_raster to get the elevation raster. Then, the raster is converted to data frame, and elevation is joined with the shapefile based on the nearest point. This procedure is approximated, but is fine for most of the large-scale analyses.

# get DEM and transform to data frame with coordinates
elevation <- get_elev_raster(shape_basin, z = 8)
catchment_DEM <- raster::as.data.frame(elevation, xy = TRUE)

# Get coordinates of the river network segments
river_net_simplified_centroids <- river_net_simplified %>%
  st_as_sf() %>%
  st_centroid()

# Get coordinates of both elements for joining
Coord_Edges <- st_coordinates(river_net_simplified_centroids) #Coordinates of the joins
Coord_DEM <- catchment_DEM[,1:2] #Coordinates of the dams

# Matching each centroid with its closer altittude point to later obtain the altitudes
matching_altitudes <- RANN::nn2(data=Coord_DEM, query = Coord_Edges, k=1, radius = 1)[[1]] 

# Get values and add to the river shapefile
catchment_DEM <- catchment_DEM[matching_altitudes,3]
river_net_simplified <- river_net_simplified %>% 
  mutate(alt = catchment_DEM)

# To avoid issues in the network creation process, retain only the important columns
river_net_simplified <- river_net_simplified %>% 
  dplyr::select(NodeID, length, alt, DIST_DN_KM, UPLAND_SKM)

Plotting can show weird behaviors of the join operation.

ggplot() +
  coord_fixed() +
  ggspatial::layer_spatial(river_net_simplified, aes(color = alt))+
  scale_color_viridis(name = "Elevation")+
  theme_minimal() +
  theme(legend.position = "bottom")+
  ggspatial::annotation_north_arrow(location = "br")+
  ggspatial::annotation_scale(location = "bl", style = "ticks")

3.1 Editing the igraph object

The edges and verices attributes can be further edited. For instance, the ‘length’ vertices attributes is expressed in meters. The value is changed to tens of kilometers for the following calculations.

# update length attribute
V(river_graph)$length <- V(river_graph)$length / 10000
hist(V(river_graph)$length)

Next, the name attribute is converted from numeric to charachter to avoid issues with the index_calculation function.

# update length attribute
V(river_graph)$name <- as.character(V(river_graph)$name)

Another step is to add some fields with habitat suitability information. Two different habitat curves are implemented to link elevation and habitat preference, one for high-altitude organisms and one for low-altitude organisms. Relative reach-specific habitat suitability indices (HSI) are added as attributes of the vertices. Weighted usable length (WUL) is also calculated as the product of length and HSI as the reach-specific fraction of length that is available for organisms.

# Function for organism that prefers high elevation
suit_fun_high <- function(x){dnorm(x, mean = 1500, sd = 500)*410*sqrt(3*pi)}

# Function for organism that prefers low elevation
suit_fun_low <- function(x){exp(- 0.001*x)}

# Calculate HSI for the igraph nodes
V(river_graph)$HSI_low <- suit_fun_low(V(river_graph)$alt)
V(river_graph)$HSI_high <- suit_fun_high(V(river_graph)$alt)

# Calculate weighted usable length for igraph nodes
V(river_graph)$WUL_low <- V(river_graph)$HSI_low * V(river_graph)$length
V(river_graph)$WUL_high <- V(river_graph)$HSI_high * V(river_graph)$length

# Plot the two response functions
1:10:2500 %>%
  data.frame("Elevation" = .,
             "Low" = suit_fun_low(.), 
             "High" = suit_fun_high(.)) %>%
  pivot_longer(c("Low", "High"), 
               names_to = "Type", values_to = "HSI") %>%
  ggplot() + 
  geom_line(aes(x = Elevation, y = HSI, color = Type))+
  theme_bw() + xlab("Elevation (m)") + ylab("Habitat Suitability Index (HSI)")

3.2 Plotting the igraph object

The ggnetwork package can be used to fortify igraph objects for ggplot-friendly plots. For a nice network layout that resembles the geometry of the starting network, the centroids of the reaches are extracted. The habitat suitabilities are plotted.

# Extract reaches centroids
river_net_simplified_centroids <- river_net_simplified %>%
  st_as_sf() %>%
  st_centroid() 

# get the centroids coordinates
coordinates <- river_net_simplified_centroids %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,1], lon = sf::st_coordinates(.)[,2]) %>%
  dplyr::select(lat, lon) %>%
  st_set_geometry( NULL)

# fortify the igraph object
gg0 <- ggnetwork(river_graph, layout = coordinates %>% as.matrix(), scale = FALSE)

grid.arrange(
ggplot(gg0, aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_fixed() +
  geom_nodes(aes(color = HSI_high)) +
  geom_edges(alpha = 0.5) +
  scale_color_viridis()+
  theme_blank()+
  ggtitle("High-altitude organism") +
  labs(caption = "Network directionality not shown"), 

ggplot(gg0, aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_fixed() +
  geom_nodes(aes(color = WUL_high)) +
  geom_edges(alpha = 0.5) +
  scale_color_viridis()+
  theme_blank()+
  ggtitle("High-altitude organism") +
  labs(caption = "Network directionality not shown"), 

ggplot(gg0, aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_fixed() +
  geom_nodes(aes(color = HSI_low)) +
  geom_edges(alpha = 0.5) +
  scale_color_viridis()+
  theme_blank()+
  ggtitle("Low-altitude organism") +
  labs(caption = "Network directionality not shown"),

ggplot(gg0, aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_fixed() +
  geom_nodes(aes(color = WUL_low)) +
  geom_edges(alpha = 0.5) +
  scale_color_viridis()+
  theme_blank()+
  ggtitle("Low-altitude organism") +
  labs(caption = "Network directionality not shown"),
  ncol=2, nrow=2)

4.1 Reach-scale indices

The funtion index_calculation from the riverconn package is used to calculate reach-scale connectivity indices (RCI). Refer to the package vignette for details on the calculation of the indices. For comparison purposes, different indices are calculated: (A) Dendritic connectivity index with symmetric dams passabilities; (B) Dendritic connectivity index with asymmetric dams passabilities; (C) Integral Index of connectivity (dams are not passable and a threshold on dispersal distance is used); (D) Integral Index of connectivity with uniform weights; (E) lowland-preferring organisms with active aquatic dispersal (e.g. fish), (F) upland organisms with active aquatic dispersal (e.g. fish), (G) lowland-preferring organism with passive aquatic dispersal (e.g., invertebrate larval stage), (H) upland organism with passive aquatic dispersal (e.g., invertebrate larval stage, bivalve). Note that for all those indices, the option index_mode = "from" was used, meaning paths exiting from the node are considered for calculating the coincidence probabilities. Thus, such indices can be interpreted as the potential for each node to support dispersal of organisms to the whole river network. An alternative could be to specify index_mode = "to". In this case, only incoming links are considered for the calculation of the index, yielding a prioritization based on the node potential for recolonization. Note that when symmetric dispersal and symmetric barriers passability are used, the two indices should overlap.

Before calculating indices, you should consider carefully the size of the network and its complexity. The index_calculation function might fail for large networks on comuters with limited available RAM because R would not be able to allocate enough memory to efficiently store the \(n * n\) matrices generated in the process. For example, the function was failing with a network size of 10^5 nodes on a computer with 8 GB RAM, but was succesfull on a computer with 32 GB RAM. Few tests might be needed to ballance complexity and computational resources.

# Initialize list where to store all the index calculation outputs
index <- list()
lab_index <- list()
letter_index <- list()

# 1: Symmetric Dendritic Connectivity Index (no biotic effects)
lab_index[[1]] <- "Symmetric DCI"
letter_index[[1]] <- "A"
index[[1]] <- index_calculation(graph = river_graph,
                                weight = "length",
                                B_ij_flag = FALSE,
                                index_type = "reach",
                                index_mode = "from")

# 2: Asymmetric Dendritic Connectivity Index (no biotic effects)
lab_index[[2]] <- "Asymmetric DCI"
letter_index[[2]] <- "B"
index[[2]] <- index_calculation(graph = river_graph,
                                weight = "length",
                                B_ij_flag = FALSE,
                                dir_fragmentation_type = "asymmetric",
                                index_type = "reach",
                                index_mode = "from")

## Before calculating IIC a binary passability has to be defined
E(river_graph)$pass_u_bin <- ifelse(is.na(E(river_graph)$pass_u), NA, 0)
E(river_graph)$pass_d_bin <- ifelse(is.na(E(river_graph)$pass_d), NA, 0)

# 3: Symmetric Integral Index of Connectivity
lab_index[[3]] <- "IIC"
letter_index[[3]] <- "C"
index[[3]] <- index_calculation(graph = river_graph,
                                weight = "length",
                                pass_u = "pass_u_bin",
                                pass_d = "pass_d_bin",
                                param = 3,
                                disp_type = "threshold",
                                index_type = "reach",
                                index_mode = "from")

## Defining uniform weights for IIC
V(river_graph)$unif_w <- 1

# 4: Integral Index of Connectivity with uniform weights
lab_index[[4]] <- "IIC with uniform weights"
letter_index[[4]] <- "D"
index[[4]] <- index_calculation(graph = river_graph,
                                weight = "unif_w",
                                dir_fragmentation_type = "asymmetric",
                                pass_u = "pass_u_bin",
                                pass_d = "pass_d_bin",
                                param = 3,
                                disp_type = "threshold",
                                index_type = "reach",
                                index_mode = "from")

# 5: Population Connectivity Index for lowland fish
lab_index[[5]] <- "PCI lowland fish"
letter_index[[5]] <- "E"
index[[5]] <- index_calculation(graph = river_graph,
                                weight = "WUL_low", 
                                param = 0.8,
                                index_type = "reach",
                                index_mode = "from") 

# 6: Population Connectivity Index for upland fish
lab_index[[6]] <- "PCI upland fish"
letter_index[[6]] <- "F"
index[[6]] <- index_calculation(graph = river_graph,
                                weight = "WUL_high", 
                                param = 0.8,
                                index_type = "reach",
                                index_mode = "from") 


# 7: Population Connectivity Index for lowland passive drifter
## note that units of param_d must be consistent with length field (10s of km)
lab_index[[7]] <- "PCI lowland passive"
letter_index[[7]] <- "G"
index[[7]] <- index_calculation(graph = river_graph,
                                weight = "WUL_low", 
                                dir_distance_type  = "asymmetric",
                                disp_type = "threshold", 
                                param_u = 0, 
                                param_d = 3,
                                index_type = "reach",
                                index_mode = "from")

# 8: Population Connectivity Index for upland passive drifter
## note that units of param_d must be consistent with length field (10s of km)
lab_index[[8]] <- "PCI upland passive"
letter_index[[8]] <- "H"
index[[8]] <- index_calculation(graph = river_graph,
                                weight = "WUL_high", 
                                dir_distance_type  = "asymmetric",
                                disp_type = "threshold",
                                param_u = 0, 
                                param_d = 3,
                                index_type = "reach",
                                index_mode = "from")

A quick check to the variation range.

data.frame("index" = do.call(rbind, lab_index),
           "min" = t(sapply(index, sapply, min))[,4], 
           "mean" = t(sapply(index, sapply, mean))[,4],
           "max" = t(sapply(index, sapply, max))[,4] )
##                      index                  min        mean                max
## 1            Symmetric DCI 0.000136068138828111 0.089094209  0.183348672040763
## 2           Asymmetric DCI 0.000198840008657542 0.243700325  0.575848057363519
## 3                      IIC 4.21802688098917e-05 0.007121115 0.0209559901505482
## 4 IIC with uniform weights  0.00154320987654321 0.008921087 0.0308641975308642
## 5         PCI lowland fish 5.49335271641571e-05 0.016012598 0.0536956831924345
## 6          PCI upland fish 0.000182545268578881 0.009541135 0.0410567699684218
## 7      PCI lowland passive 3.79438573745419e-05 0.003693739  0.020852867904853
## 8       PCI upland passive 1.66202468527347e-05 0.003000132 0.0353297075947445

Resulting maps can be plotted to visually compare the results. For comparisons between different indices, reach ranks are plotted. The lowest rank is the total number of reaches (ca. 650). Reaches with higher rankings offer higher potential for being dispersal sources.

# Initialize empty list
plot_list <- list()

# iterate through list length
for (i in 1:length(index)) {
  
  # Join d_i information with dams shapefile
  river_plot <- river_net_simplified %>%
    mutate(name = as.character(NodeID)) %>%
    left_join(index[[i]]) %>%
    mutate(rank = rank(desc(index)))

  # plot
  plot_iter <- ggplot() +
    coord_fixed() +
    ggspatial::layer_spatial(river_plot, aes(color = rank))+
    scale_color_viridis(name = "Ranking", direction = -1)+
    theme_void()+
    guides(size = FALSE) +
    theme(legend.position = "bottom")+
    ggtitle(paste0(letter_index[[i]], ") ", lab_index[[i]]))
  
    legend <- cowplot::get_legend(plot_iter)
  
  plot_list[[i]] <- plot_iter + 
    theme(legend.position = "none")

}

plot_list[[length(index)+1]] <- legend

plot_final <- ggpubr::ggarrange( plotlist = plot_list)
plot_final

ggsave(plot = plot_final, "Figures/habitat_prioritization.jpeg",  
       width = 18, height = 12, units = "cm")

A correlation plot helps identifying the correlated indices

d_index_chart <- data.frame("d_i_1" = index[[1]]$index, 
                            "d_i_2" = index[[2]]$index,
                            "d_i_3" = index[[3]]$index,
                            "d_i_4" = index[[4]]$index,
                            "d_i_5" = index[[5]]$index,
                            "d_i_6" = index[[6]]$index,
                            "d_i_7" = index[[7]]$index,
                            "d_i_8" = index[[8]]$index)

colnames(d_index_chart) <- c(letter_index[[1]], letter_index[[2]], letter_index[[3]], 
                             letter_index[[4]], letter_index[[5]], letter_index[[6]],
                             letter_index[[7]], letter_index[[8]])

# Spearman corrlations between rankings calculated with different priorities
cor(d_index_chart, method = "spearman")
##            A           B         C          D          E           F         G
## A  1.0000000  0.91224397 0.5542838  0.4941878  0.8789122  0.19375360 0.3691311
## B  0.9122440  1.00000000 0.4955071  0.4322901  0.8570346 -0.08416734 0.3894249
## C  0.5542838  0.49550709 1.0000000  0.8277318  0.7113399  0.41582814 0.5582716
## D  0.4941878  0.43229005 0.8277318  1.0000000  0.7109633  0.30121569 0.4021395
## E  0.8789122  0.85703457 0.7113399  0.7109633  1.0000000  0.13442831 0.5122386
## F  0.1937536 -0.08416734 0.4158281  0.3012157  0.1344283  1.00000000 0.1011033
## G  0.3691311  0.38942489 0.5582716  0.4021395  0.5122386  0.10110333 1.0000000
## H -0.2105778 -0.34538754 0.1183752 -0.1075774 -0.2894423  0.55415117 0.4008371
##            H
## A -0.2105778
## B -0.3453875
## C  0.1183752
## D -0.1075774
## E -0.2894423
## F  0.5541512
## G  0.4008371
## H  1.0000000
corrmorant::corrmorant(d_index_chart, corr_method = "spearman", style = "binned")+
  theme(legend.position = "bottom") +
  scale_color_viridis(name = "Correlation", option = "E", direction = -1, limits = c(-1, 1)) +
  theme(axis.text = element_text(size=5),
        axis.text.x = element_text(angle = 45) )

ggsave("Figures/ranking_habitat_comparison.jpeg",  width = 15, height = 15, units = "cm")

4.2 Barriers prioritization

The function d_index_calculation from the riverconn package is used to identify barriers whose passability improvement can best benefit the catchment connectivity. Refer to the package vignette for details on the calculation of the indices. Check the package vignette for details on the indices used to this end. Here we use the same setups described above (section 4.1).

First dams metadata are to be defined. The object is saved as a dataframe that stores the id of the dams to be individually tested in a column, and the relative improvements in upstream and downstream passabilies (here improvement is assumed to be uniformly 1).

barriers_metadata <- data.frame("id_barrier" =  E(river_graph)$id_barrier[!is.na(E(river_graph)$id_barrier)],
                            "pass_u_updated" = 1,
                            "pass_d_updated" = 1)
head(barriers_metadata)
##   id_barrier pass_u_updated pass_d_updated
## 1         76              1              1
## 2         72              1              1
## 3          7              1              1
## 4          6              1              1
## 5         78              1              1
## 6         51              1              1

Next, the dams prioritization can be executed. The calculation may take a while depending on the computer where it is run. There is the option to parallelize the process (marking -> parallel=TRUE and the number of cores devoted -> ncores= ). It takes 4 minutes with parallel=FALSE for this example in a computer with this procesor: 11th Gen Intel Core i7, 3.30 GHz and 32 GbRAM. On the other hand it takes 1.5 minutes when parallelizing with 7 cores (parallel = TRUE, ncores = 7). To parallelize calculations, the doParallel package is required.

# Initialize list where to store all the index calculation outputs
d_index <- list()
lab_d_index <- list()
letter_d_index <- list()

# 1: Symmetric Dendritic Connectivity Index (no biotic effects)
lab_d_index[[1]] <- "Symmetric DCI"
letter_d_index[[1]] <- "A"
d_index[[1]] <- d_index_calculation(graph = river_graph, 
                                    barriers_metadata = barriers_metadata,
                                    B_ij_flag = FALSE, 
                                    parallel = FALSE)

# 2: Asymmetric Dendritic Connectivity Index (no biotic effects)
lab_d_index[[2]] <- "Asymmetric DCI"
letter_d_index[[2]] <- "B"
d_index[[2]] <- d_index_calculation(graph = river_graph, 
                                    barriers_metadata = barriers_metadata, 
                                    B_ij_flag = FALSE, 
                                    parallel = FALSE, 
                                    dir_fragmentation_type = "asymmetric")

## Before calculating IIC a binary passability has to be defined
E(river_graph)$pass_u_bin <- ifelse(is.na(E(river_graph)$pass_u), NA, 0)
E(river_graph)$pass_d_bin <- ifelse(is.na(E(river_graph)$pass_d), NA, 0)

# 3: Symmetric Integral Index of Connectivity
lab_d_index[[3]] <- "Symmetric IIC"
letter_d_index[[3]] <- "C"
d_index[[3]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "WUL_low",
                                    pass_u = "pass_u_bin",
                                    pass_d = "pass_d_bin",
                                    param = 3,
                                    disp_type = "threshold",
                                    parallel = FALSE)

## Defining uniform weights for IIC
V(river_graph)$unif_w <- 1

# 4: Integral Index of Connectivity with uniform weights
lab_d_index[[4]] <- "IIC with uniform weights"
letter_d_index[[4]] <- "D"
d_index[[4]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "unif_w",
                                    dir_fragmentation_type = "asymmetric",
                                    pass_u = "pass_u_bin",
                                    pass_d = "pass_d_bin",
                                    param = 3,
                                    disp_type = "threshold",
                                    parallel = FALSE)

# 5: Population Connectivity Index for lowland fish
lab_d_index[[5]] <- "PCI lowland fish"
letter_d_index[[5]] <- "E"
d_index[[5]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "WUL_low", 
                                    param = 0.8,
                                    parallel = FALSE) 

# 6: Population Connectivity Index for upland fish
lab_d_index[[6]] <- "PCI upland fish"
letter_d_index[[6]] <- "F"
d_index[[6]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "WUL_high", 
                                    param = 0.8,
                                    parallel = FALSE )


# 7: Population Connectivity Index for lowland passive drifter
## note that units of param_d must be consistent with length field (10s of km)
lab_d_index[[7]] <- "PCI lowland passive"
letter_d_index[[7]] <- "G"
d_index[[7]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "WUL_low", 
                                    dir_distance_type  = "asymmetric",
                                    disp_type = "threshold", 
                                    param_u = 0, 
                                    param_d = 3,
                                    parallel = FALSE)

# 8: Population Connectivity Index for upland passive drifter
## note that units of param_d must be consistent with length field (10s of km)
lab_d_index[[8]] <- "PCI upland passive"
letter_d_index[[8]] <- "H"
d_index[[8]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "WUL_high", 
                                    dir_distance_type  = "asymmetric",
                                    disp_type = "threshold",
                                    param_u = 0, 
                                    param_d = 3,
                                    parallel = FALSE)

# 9: Catchment Area Fragmentation Index
## note that units of param_d must be consistent with length field (10s of km)
lab_d_index[[9]] <- "CAFI"
letter_d_index[[9]] <- "I"
d_index[[9]] <- d_index_calculation(graph = river_graph,
                                    barriers_metadata = barriers_metadata, 
                                    weight = "UPLAND_SKM", 
                                    dir_distance_type  = "symmetric",
                                    B_ij_flag = FALSE,
                                    parallel = FALSE)

A quick check to the variation range of CCI for the baseline.

data.frame("Baseline CCI" = do.call(rbind, lab_d_index),
           "min" = t(sapply(d_index, sapply, min))[,4], # baseline in is column 4
           "mean" = t(sapply(d_index, sapply, mean))[,4],
           "max" = t(sapply(d_index, sapply, max))[,4] )
##               Baseline.CCI                 min        mean                 max
## 1            Symmetric DCI  0.0909335515845068 0.092566524   0.126107269330024
## 2           Asymmetric DCI    0.25182200888365 0.254232665   0.290384217177615
## 3            Symmetric IIC 0.00810297139863181 0.008114263 0.00817098728639399
## 4 IIC with uniform weights 0.00892108672458467 0.008945654 0.00905445054107606
## 5         PCI lowland fish  0.0173825827122186 0.017454205  0.0183293467389933
## 6          PCI upland fish  0.0149304263813881 0.014969916  0.0155246121675143
## 7      PCI lowland passive 0.00494555720870629 0.004950211 0.00497254428181167
## 8       PCI upland passive 0.00811555022541386 0.008120642 0.00819221692543936
## 9                     CAFI   0.203495022267748 0.206691818   0.299794408043562

Resulting maps can be plotted to visually compare the results. For comparisons between different indices, barriers ranks are plotted. The lowest rank is the total number of barriers (ca. 80). Barriers with higher ranks should be prioritized for passability improvement.

# Initialize empty list
plot_list <- list()

# iterate through list length
for (i in 1:length(d_index)) {
  
  # Join d_i information with dams shapefile
  network_links_plot <- network_links %>%
    filter(type == "dam") %>% 
    left_join(d_index[[i]]) %>%
    mutate(d_rank = rank(desc(d_index)))

  # plot
  plot_iter <- ggplot() +
    coord_fixed() +
    ggspatial::layer_spatial(river_net_simplified, color = "gray70")+
    ggspatial::layer_spatial(network_links_plot, 
                             aes(color = d_rank, size = 1/d_rank), alpha = 0.8)+
    scale_color_viridis(name = "Ranking", direction = -1)+
    theme_void()+
    guides(size = FALSE) +
    theme(legend.position = "bottom")+
    ggtitle(paste0(letter_d_index[[i]], ") ", lab_d_index[[i]]))
    
  # legend <- cowplot::get_legend(plot_iter)
  
  plot_list[[i]] <- plot_iter #+ theme(legend.position = "none")

}

# plot_list[[length(d_index)+1]] <- legend

plot_final <- ggpubr::ggarrange( plotlist = plot_list, 
                                 common.legend = TRUE, legend = "bottom")
plot_final

ggsave(plot = plot_final, "Figures/barriers_prioritization.jpeg",  
       width = 18, height = 12, units = "cm")

A corrleation plot hepls identifying the correlated indices

d_index_chart <- data.frame("d_i_1" = d_index[[1]]$d_index, 
                            "d_i_2" = d_index[[2]]$d_index,
                            "d_i_3" = d_index[[3]]$d_index,
                            "d_i_4" = d_index[[4]]$d_index,
                            "d_i_5" = d_index[[5]]$d_index,
                            "d_i_6" = d_index[[6]]$d_index,
                            "d_i_7" = d_index[[7]]$d_index,
                            "d_i_8" = d_index[[8]]$d_index,
                            "d_i_9" = d_index[[9]]$d_index)

colnames(d_index_chart) <- c(letter_d_index[[1]], letter_d_index[[2]], letter_d_index[[3]], 
                             letter_d_index[[4]], letter_d_index[[5]], letter_d_index[[6]],
                             letter_d_index[[7]], letter_d_index[[8]], letter_d_index[[9]])

# Spearman corrlations between rankings calculated with different priorities
cor(d_index_chart, method = "spearman")
##            A          B           C          D          E         F          G
## A  1.0000000  0.8636478 0.637634905 0.26204708  0.9408875 0.3693169 0.64327608
## B  0.8636478  1.0000000 0.648567309 0.29080852  0.8753639 0.2381299 0.67678719
## C  0.6376349  0.6485673 1.000000000 0.54545734  0.7551617 0.2080936 0.97127579
## D  0.2620471  0.2908085 0.545457344 1.00000000  0.3773260 0.2458078 0.47141151
## E  0.9408875  0.8753639 0.755161746 0.37732603  1.0000000 0.2660974 0.73870015
## F  0.3693169  0.2381299 0.208093618 0.24580782  0.2660974 1.0000000 0.23057434
## G  0.6432761  0.6767872 0.971275792 0.47141151  0.7387001 0.2305743 1.00000000
## H -0.2640295 -0.3216590 0.009140793 0.03524209 -0.3429779 0.5990565 0.05841429
## I  0.8891517  0.9147256 0.640364507 0.30417163  0.8815929 0.1295633 0.65949971
##              H          I
## A -0.264029452  0.8891517
## B -0.321659038  0.9147256
## C  0.009140793  0.6403645
## D  0.035242088  0.3041716
## E -0.342977925  0.8815929
## F  0.599056538  0.1295633
## G  0.058414289  0.6594997
## H  1.000000000 -0.3937625
## I -0.393762511  1.0000000
corrmorant::corrmorant(d_index_chart, corr_method = "spearman", style = "binned")+
  theme(legend.position = "bottom") +
  scale_color_viridis(name = "Correlation", option = "E", direction = -1, limits = c(-1, 1)) +
  theme(axis.text = element_text(size=5),
        axis.text.x = element_text(angle = 45))

ggsave("Figures/ranking comparison.jpeg",  width = 15, height = 15, units = "cm")