1.1 Packages needed

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

• sf,lwgeom,rgdal,raster,elevatr,RANN packagesfor obtaining and processing spatial data;
• igraph package forperforming network analyses;
• tidyversepackages for efficient and easy-to-read data pipeline (seehere andhere for the most frequentlyused functions);
• ggnetworkfor producing easy, ggplot-like, plot of classigraphobjects andggspatialfor easy, ggplot-like plot, of spatial data;
• corrmorantto plot nice correlation plots;
  
• riverconn to calculate the fragmentation indices forthe 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 onsf objects. Thanks to their data.frame-like structure,spatial operations (e.g., selecting and updating columns and rows, joinwith other data frames) are “easier” on such objects. Most spatialobjects are converted to this class as soon as possible in thescript.

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 theigraph object requires at leastsome information on the spatial geometry of the river system and theposition of the barriers.

The river network shapefile can be retrieved fromHydroSHEDS, a globalhydrographic atlas with geo-referenced datasets derived from spacebornelevation data. HydroSHEDS is widely used for global to regionalanalyses (Lehner et al., 2008). River network data (HydroRIVERS) can bedownloaded fromhere.Basins and catchments shapefiles (HydroBASINS) are also available with12 detail levels (Lehner et al., 2013). The Ebro catchment border can beselected from the polygons available in the HydroBASINS layer (use level5 basins shapefile for the identification) and the corresponding networkcan be extracted from the HydroSHEDS layer.

Other options are available for catchments and rivers.TheEuropean Catchments and Rivers network system (ECRINS) datasetprovides good quality data covering most of Europe. The ECRINS data forrivers are structured with segments and nodes, making the creation of anetwork quite straight-forward. ECRINS has a higher resolution thanHydroRIVERS. The recentHydrography90mdataset can also be used. Local sources (e.g. hydrography layers fromlocal authorities) are also generally available, but particular careshould be taken in preprocessing to correct the presence of disconnectedsegments and loops in the shapefile. Different HydroSHEDS andHydroRIVERS versions are available. This tutorial uses HydroRIVERS v1.1:note that the previous versions were slighthly different in terms ofattributes (for further details see the HydroRIVERS technicaldocumentation).

The shapefile with the geographic position of the dams was obtainedfrom the Ebro river public Geographic Information SystemCHE website.

TheAMBERdataset described in Belletti et al. (2020) contains barriers datacompiled from different national and international sources. Metadata onbarriers type (weir, dam, culvert, etc) and height are available butmost barriers are unclassified.

The corresponding shapefiles are available here and can be importedin R (note that the river shapefile “Ebre_riv_modif.shp” has beensligtly modified from the original hydrosheds shapefile -see nextsections on the checks carried out on the shapefiles. The originalshapefile is saved in the Ebro_shape_original/ folder forcomparison).

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 84shape_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 84shape_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 dataprocessing. The source code and the documentation for these functions isfound in Appendix B. Before starting, make sure those functions areloaded in the global environment:get_outlet,create_network,snap_to_river,headwaters_dam,multiple_confluences,gaps_detection,check_components. Detailedinformation on the use of these functions is presented in the nextsections.

2.1 Shapefiles processing

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

Note that different HydroRIVERS versions are available. This tutorialuses HydroRIVERS v1.1, wuere different fields are available for pruningthe river network based on the size (e.g. ORD_STRA, ORD_CLAS, ORD_FLOW).Previous versions have different attributes. For instance, inHydroRIVERS v1.0 the field UP_CELLS (number of upstream cells) can beused to prune the river network. Check the datasetdocumentationfor further details.

# Set a threshold of 80 square kilometersthreshold = 80# Prune HydroRIVERS network based on upstream areashape_river_small <- shape_river[as.numeric(shape_river$UPLAND_SKM) > threshold,]

Next, the dams shapefile can be converted to a shapefile and a uniqueid can be given to each point. As the dams dataset is a polygonshapefile the centroid of each poligon has to be calculated to obtain a“point” layer. The reference system of the dams layer is changed toWGS84, 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 thegeolocalization 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 tworeaches merge into one. In the river graph defined above, confluencesare edges (links) between nodes. To define confluences the pruned rivershapefile is first simplified and its geometry is casted to “POINT”.Note that for large shapefiles the functionrgeos::gLineMerge might be more efficient thansf::st_union.

# Simplify river shapefileshape_river_simple <- shape_river_small %>%  st_as_sf %>%  st_union()# Convert shapefile to point objectriver_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 theoriginal shapefile, a confluence is delineated where three (or more) ofsuch nodes are overlapping.

# Check points that overlapjoins_selection <- river_to_points %>%  st_equals() %>%  Filter(function(x){length(x) > 2}, .) %>%  lapply(., FUN = min) %>%  unlist() %>%  unique()# Filter original point shapefile to retain only confluencesriver_joins <- river_to_points %>%   filter(id %in% joins_selection)

The last step is the splitting of the simplified shapefile where itintersects a confluence. This allows to identify unequivocally the riversections that are comprised between two confluences. Such elements willthen became vertices (nodes) of the graph. The functionlwgeom::st_split is used. This function splits an object ofclass ‘sf’ with ‘MULTILINESTRING’ geometry into and object of class ‘sf’containing multiple ‘LINESTRING’ geometries records based on theposition of the points. Note that for this function to work properly,the point object must perfectly overlap with the polyline (this can bechecked with the functionsf::st_distance applied to thepoints and the polyline, which must return a vector of zeroes).

# Split polylineshape_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 ideais to obtain a shapefile that is sliced in points where dams orconfluences are located. However, dams shapefiles require special carein the processing of the data. First, dams need to be snapped to theriver network. Second, dams that are too far away from the river networkneed to be removed from the analysis because they are probably locatedin smaller streams that were removed beforehand.

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

# Snap damsdams_snapped <- snap_to_river(dams_to_points,                              shape_river_simple %>% st_sf(),                              max_dist = 1000)# Retain dams that were snappeddams_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 samegeographic point, and to assign passablity values to the obtainedsnapped points. Here, uniform (not dam-dependent) passabilities wereassigned: 0.8 for downstream passability, and 0.1 for upstreampassability. An id is assigned to each of the geographic points.Finally, each dam is assigned the value of the upstream area with aspatial join with the original river shapefile (HydroSheds expresses theupstream 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 functionheadwaters_dam (see Appendix B) canidentify dams that are located in the headwaters, meaning that they donot have any river network segment upstream. When no upstream segment ispresent, the network can not be created (because the dam would be anedge without an upstream vertex). When calling the flag_headwater youshould 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 theposition of dams and confluences. First a combined ‘junctions’ datasetis created (the network_links object), then the functionlwgeom::st_split is applied again. note that theconfluences are assigned NA values to the fields present in the damsshapefile. Note also that in the river network shapefile, the fieldlength is calculated. Finally, a spatial join is performed with theoriginal river shapefile to assign to each reach the correspondingupstream area.

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

# Create junction point shapefilenetwork_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 networkriver_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 confluenceswhere four or more rivers join. Such type of confluences can causeproblems in the next steps. The functionmultiple_confluences (Appendix B) identifies suchproblematic confluences so that they can be corrected in a GISenvironment. The correction consists in editing the river shapefile tomove re-locate one of the reaches (either upstream or downstream) on theintersections that have more than three reaches (this can be easilyachieved once these intersections are located by the followingfunction).

Once edited, re-run the tutorial and check again that there are nomore 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 thepresence of gaps between segments. This issue can cause problems in thenetwork_creation function, which expects a shapefile withno gaps. The functionsgaps_detection andcheck_components (Appendix B) can be used to detect suchproblematic points in the shapefile. In particular, the functioncheck_components returns an sf file that is a copy of theoriginal river_net_simplified with the additional field component thatcontains a progressive number identifying multiple components. A correctshapefile should contain only ‘1’ values in this field. A plot can helpidentifying the location of those components and a manual correction ina gis environment can fix the problem.

Once edited, re-run the tutorial and check again that the shapefilehas 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 addedto the river shapefile. Elevation is an example. Here we use thefunctionelevatr::get_elev_raster to get the elevationraster. Then, the raster is converted to data frame, and elevation isjoined with the shapefile based on the nearest point. This procedure isapproximated, but is fine for most of the large-scale analyses.

# get DEM and transform to data frame with coordinateselevation <- get_elev_raster(shape_basin, z = 8)catchment_DEM <- raster::as.data.frame(elevation, xy = TRUE)# Get coordinates of the river network segmentsriver_net_simplified_centroids <- river_net_simplified %>%  st_as_sf() %>%  st_centroid()# Get coordinates of both elements for joiningCoord_Edges <- st_coordinates(river_net_simplified_centroids) #Coordinates of the joinsCoord_DEM <- catchment_DEM[,1:2] #Coordinates of the dams# Matching each centroid with its closer altittude point to later obtain the altitudesmatching_altitudes <- RANN::nn2(data=Coord_DEM, query = Coord_Edges, k=1, radius = 1)[[1]] # Get values and add to the river shapefilecatchment_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 columnsriver_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 theigraph object

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

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

Next, the name attribute is converted from numeric to charachter toavoid issues with theindex_calculation function.

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

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

# Function for organism that prefers high elevationsuit_fun_high <- function(x){dnorm(x, mean = 1500, sd = 500)*410*sqrt(3*pi)}# Function for organism that prefers low elevationsuit_fun_low <- function(x){exp(- 0.001*x)}# Calculate HSI for the igraph nodesV(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 nodesV(river_graph)$WUL_low <- V(river_graph)$HSI_low * V(river_graph)$lengthV(river_graph)$WUL_high <- V(river_graph)$HSI_high * V(river_graph)$length# Plot the two response functions1: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 theigraph object

Theggnetwork package can be used to fortifyigraph objects for ggplot-friendly plots. For a nicenetwork layout that resembles the geometry of the starting network, thecentroids of the reaches are extracted. The habitat suitabilities areplotted.

# Extract reaches centroidsriver_net_simplified_centroids <- river_net_simplified %>%  st_as_sf() %>%  st_centroid() # get the centroids coordinatescoordinates <- 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 objectgg0 <- 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 funtionindex_calculation from theriverconnpackage is used to calculate reach-scale connectivity indices (RCI).Refer to the package vignette for details on the calculation of theindices. 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 ondispersal distance is used); (D) Integral Index of connectivity withuniform weights; (E) lowland-preferring organisms with active aquaticdispersal (e.g. fish), (F) upland organisms with active aquaticdispersal (e.g. fish), (G) lowland-preferring organism with passiveaquatic dispersal (e.g., invertebrate larval stage), (H) upland organismwith passive aquatic dispersal (e.g., invertebrate larval stage,bivalve). Note that for all those indices, the optionindex_mode = "from" was used, meaning paths exiting fromthe node are considered for calculating the coincidence probabilities.Thus, such indices can be interpreted as the potential for each node tosupport dispersal of organisms to the whole river network. Analternative could be to specifyindex_mode = "to". In thiscase, only incoming links are considered for the calculation of theindex, yielding a prioritization based on the node potential forrecolonization. Note that when symmetric dispersal and symmetricbarriers passability are used, the two indices should overlap.

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

# Initialize list where to store all the index calculation outputsindex <- 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 definedE(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 Connectivitylab_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 IICV(river_graph)$unif_w <- 1# 4: Integral Index of Connectivity with uniform weightslab_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 fishlab_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 fishlab_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. Forcomparisons between different indices, reach ranks are plotted. Thelowest rank is the total number of reaches (ca. 650). Reaches withhigher rankings offer higher potential for being dispersal sources.

# Initialize empty listplot_list <- list()# iterate through list lengthfor (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]] <- legendplot_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 prioritiescor(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.0000000corrmorant::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 functiond_index_calculation from theriverconnpackage is used to identify barriers whose passability improvement canbest benefit the catchment connectivity. Refer to the package vignettefor details on the calculation of the indices. Check the packagevignette for details on the indices used to this end. Here we use thesame setups described above (section 4.1).

First dams metadata are to be defined. The object is saved as adataframe that stores the id of the dams to be individually tested in acolumn, and the relative improvements in upstream and downstreampassabilies (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 maytake a while depending on the computer where it is run. There is theoption to parallelize the process (marking ->parallel=TRUE and the number of cores devoted -> ncores=). It takes 4 minutes withparallel=FALSE for this examplein a computer with this procesor: 11th Gen Intel Core i7, 3.30 GHz and32 GbRAM. On the other hand it takes 1.5 minutes when parallelizing with7 cores (parallel = TRUE,ncores = 7). Toparallelize calculations, thedoParallel package isrequired.

# Initialize list where to store all the index calculation outputsd_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 definedE(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 Connectivitylab_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 IICV(river_graph)$unif_w <- 1# 4: Integral Index of Connectivity with uniform weightslab_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 fishlab_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 fishlab_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. Forcomparisons between different indices, barriers ranks are plotted. Thelowest rank is the total number of barriers (ca. 80). Barriers withhigher ranks should be prioritized for passability improvement.

# Initialize empty listplot_list <- list()# iterate through list lengthfor (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]] <- legendplot_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 prioritiescor(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.0000000corrmorant::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")