dci User Guide
library(dci)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
The dci package provides an R interface for the measurement of
connectivity in river networks with the dendritic connectivity index
(Cote et al., 2009). The main class of the package is the river_net
class which links river lines to barriers and the outlet of the river
network. It inherits the structure of the sfnetwork class from the sfnetworks package
which itself inherits the tbl_graph class from the tidygraph
package further inherited by the igraph class from the igraph package. These dependencies on
other packages provide a rich set of tools and algorithms which remain
applicable to all instances of the river_net class. Therefore, further
and more complex analyses can be carried out with functions from the
igraph, tidygraph, and sfnetworks
packages. Further, the tbl_graph class from tidygraph is
created with the tidyverse in
mind and, thus, allows for the development of workflows using
tidyverse packages.
This vignette illustrates the dci package using river
and barrier data from the Yamaska watershed in Southern Quebec.
Required data
DCI analysis with the dci package requires three pieces
of spatial data:
- River lines: Data representing river lines which will make
up the edges of the spatial network. The
import_rivers
function of the package will extract the largest fully connected
component of the provided data and other isolated rivers will be
discarded.
- Barrier points: Points representing barriers along the
river network. Points are not required to be perfectly snapped to the
river lines, the tolerance parameter of the
river_net
function can perform this snapping.
- Outlet point: A single point representing the downstream
outlet of the river network. This point should be accurate for the
calculation of the diadromous DCI.
Importing
A river_net object is constructed using three pieces of
required data: river lines, barrier points, and the outlet location.
These input spatial data must first be imported with
import_rivers and import_points functions.
Rivers can be imported from either the path to a shapefile or an
object of the sf class. This function performs some basic
error checking, extracts the largest fully connected component in the
supplied data, and calculates river lengths. A plot with changes
highlighted in red is printed alongside the function unless the
quiet parameter is set to TRUE. River
importing is illustrated below with the path to a shapefile representing
the Yamaska river lines.
# Load Yamaska river data
data("yamaska_rivers", package = "dci")
# Import rivers
rivers_in <- import_rivers(rivers = yamaska_rivers)
Similarly, all point data can be imported with the
import_points function along with the type
parameter specifying the type of point. Possible values are: “bars” for
barriers and “out” for the outlet.
# Load Yamaska barrier and outlet data
data("yamaska_barriers", package = "dci")
data("yamaska_outlet", package = "dci")
# Import barriers
barriers_in <- import_points(pts = yamaska_barriers, type = "bars")
outlet_in <- import_points(pts = yamaska_outlet, type = "out")
barriers_in
#> Simple feature collection with 14 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -353726.8 ymin: 233223.5 xmax: -337065.2 ymax: 241340.8
#> Projected CRS: NAD83 / Quebec Lambert
#> # A tibble: 14 × 4
#> pass_1 pass_2 geometry type
#> * <dbl> <dbl> <POINT [m]> <chr>
#> 1 0.3 0 (-353726.8 237532.9) barrier
#> 2 0.5 0 (-352055.2 236475.2) barrier
#> 3 0.1 0 (-351647.6 233223.5) barrier
#> 4 0.6 0 (-349660.2 235578.7) barrier
#> 5 0.8 0 (-349237.3 237746.8) barrier
#> 6 0.9 0 (-349254.4 239658.6) barrier
#> 7 0.8 0 (-344910.5 241340.8) barrier
#> 8 0.7 0 (-344990.2 239352.8) barrier
#> 9 0.4 0 (-343315.7 236427.2) barrier
#> 10 0.3 0 (-339418.4 235113.6) barrier
#> 11 0.9 0 (-337065.2 235843) barrier
#> 12 0.5 0 (-345792.5 233877.2) barrier
#> 13 0.6 0 (-346344.7 234715.3) barrier
#> 14 0.6 0 (-347463.1 234803.6) barrier
Network construction
Once all required data have been pre-processed with the import
functions they can be joined together in the river_net
function to generate an object of the river_net class.
Messages are printed indicating the topological corrections that are
being performed on the river network. These are explained in more detail
in the next section. Additionally, this function will split rivers at
new barrier node locations.
# Combine rivers, barriers, and outlet
# Warnings about river splitting suppressed since rivers are already split at barrier node locations
net <- suppressWarnings(
river_net(rivers = rivers_in, barriers = barriers_in, outlet = outlet_in)
)
#> Iteration 1:
#> 36 divergences corrected.
#> 4 complex confluences found.
#> Iteration 2:
#> No divergences detected.
#> No complex confluences found.
plot(net)
Dendritic topology
Since the river_net class inherits from the
sfnetworks class it shares the enforcement of certain
network rules from the sfnetworks package: nodes should
have POINT geometries, edges should be spatially explicit
and have “LINESTRING” geometries, and both nodes and edges should have
the same CRS. In addition to these, the dci package
enforces additional rules pertaining to the dendritic topology of river
networks. A dendritic topology is defined as a directional network
topology in which at most two upstream edges combine at a node to form a
single outgoing downstream edge. In the case of a barrier there would be
one incoming edge while in the case of natural confluences two incoming
edges combine. This topological rule is enforced because it allows for
the design of custom river network path finding algorithms which run
quicker than traditional network algorithms. This is explained in more
detail in the labeling section.
The river_net provides functionality for correcting some
common topological errors. These errors are illustrated below:
- Divergent rivers - Divergent rivers occur where a single
upstream river splits into two downstream rivers. This often occurs
because of large obstacles in the waterway of large rivers causing a
single river to split and recombine further downstream. If the
river_net function is run with the clean
parameter set to TRUE divergences are corrected by deleting
the shortest of the two downstream rivers. If more fine-tuned
corrections are preferred the enforce_dendritic function
can be called directly on the rivers imported by
import_rivers. This functionality is illustrated
below.
- Complex confluences - Complex confluences occur when a
confluence node has over 3 incident edges. These often occur as a
limitation of data resolution where two very close confluences merge
into one. If the
river_net function is run with the
clean parameter set to TRUE complex
confluences are corrected by moving one of the incident edges further
downstream on the main channel. When a complex confluence is present
with over 4 incident edges the function will quit and a manual
correction will be required.
By default the river_net constructor performs these
topological corrections automatically. However, if manual corrections
are desired the enforce_dendritic function can be called
directly with the river lines data and the correct
parameter set to FALSE. Then, when constructing the
river_net object, the check parameter can be
set to FALSE since the topology has already been
corrected.
# Identify topological errors in rivers
river_errors <- enforce_dendritic(rivers = rivers_in, correct = FALSE)
#> 27 divergences have been found.
#> 7 complex confluences found.
river_errors
#> Simple feature collection with 620 features and 4 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: -357300.6 ymin: 230498.2 xmax: -335432.5 ymax: 245955.3
#> Projected CRS: NAD83 / Quebec Lambert
#> # A tibble: 620 × 5
#> qual riv_length geometry complexID divergent
#> <dbl> <dbl> <LINESTRING [m]> <int> <int>
#> 1 3 262. (-347139.6 241142.8, -347058.6 240893.6) NA NA
#> 2 5 51.8 (-346605.8 241063.7, -346602.3 241061.8… 1 3
#> 3 1 501. (-347058.6 240893.6, -347052.4 240896.2… 1 NA
#> 4 8 1345. (-348305.1 240996.5, -348304.5 240994.9… NA NA
#> 5 7 75.4 (-344551 237776.1, -344566.4 237776, -3… NA NA
#> 6 7 599. (-344048.9 237965.3, -344620.7 237788.2) NA NA
#> 7 6 204. (-342723.6 236820.6, -342727.1 236872, … NA NA
#> 8 3 59.3 (-354145.2 237109.7, -354140.1 237111.4… NA NA
#> 9 0 475. (-354088.8 237127.9, -354066.1 237137.3… NA NA
#> 10 5 323. (-354207.6 237379.7, -354219.2 237343.4… NA NA
#> # ℹ 610 more rows
In the above example a manual dendritic correction is requested which
returns the input river lines augmented with a divergent
variable which links divergent rivers from the same pair and a
complex variable which indicates which rivers take part in
complex confluences. These variables can then be used in a dedicated GIS
software or R and the corrected rivers can be used to construct a
river_net object. This manual correction might be
beneficial in situations where an important barrier or point of interest
is lost during automatic correction.
# Visualize divergence locations
plot(sf::st_geometry(rivers_in), col = "grey75")
plot(river_errors["divergent"], add = TRUE)
# Visualize complex confluence locations
plot(sf::st_geometry(rivers_in), col = "grey75")
plot(river_errors["complexID"], add = TRUE)
Labeling
There are 2 types of labeling performed on a river_net
object.
Node labeling - Node labeling makes use of a binary code to,
essentially, give an address to each node in the network relative to the
outlet. Labeling begins at the outlet with a label of 0 and proceeds
upstream. At the first split, the two children of the outlet node are
given labels 00 and 01. These labels allow for rapid path tracing in
river networks to gather the watercourse distance between two nodes or
the passability of barriers which lie on that path. The task of finding
the least-cost path between points is a computationally intensive
network theory problem and it’s an important one when considering
connectivity in terrestrial landscapes. However, river networks which
follow a dendritic topology exhibit an important property - there is
only one path between any two points in the network. Common path-finding
algorithms like Djikstra’s algorithm are made to be general and don`t
make these types of assumptions about network structure. This labeling
allows for the design of custom, rapid algorithms to gather the path
between nodes in the network.
Membership labeling - Membership labeling is much more simple
than the node labeling. Membership labels simply classify nodes into
specific isolated river segments bounded on all sides by either barriers
or terminal nodes of the network (sources and the outlet). These labels
are simply integers from 0 to the number of segments.
Calculating the DCI
Once a valid river_net object is constructed all the
elements are present to calculate the DCI. The DCI comes in a few
flavours and they can all be derived from the single
calculate_dci function. The parameters of the function will
dictate what type of connectivity is being estimated:
- the
form parameter - The form of the DCI relates to the
dispersal style of the target organism. The accepted forms are
“potamodromous” and “diadromous”. This will dictate the form of the DCI
equation that is calculated, more details can be found in (Cote et al.,
2009). This is the only parameter that is required to calculate the
DCI.
- the
pass parameter - The pass parameter
allows for the selection of a relevant passability variable in the
barrier data. Passability probabilities for barriers must range from 0
to 1. If they are not within that range the data will automatically be
normalized. If the parameter is left blank all barriers will receive a
passability of 0.
- the
weight parameter - Similar to the pass
parameter, the weight parameter allows for the selection of
a relevant weighting variable in the river data. Again, these values are
required to be between 0 and 1 since they are applied to the river
lengths. Weights of 0 will fully exclude rivers while weights of 1 will
fully include them. A possible application is to use habitat suitability
estimates as river weighting to achieve a more specific measure of
landscape connectivity for a species.
- the
threshold parameter - The threshold
parameter gives the option of considering a dispersal limit when
calculating the DCI. The thresholded DCI prunes any connections between
segments that are further away than the dispersal limit. The
non-thresholded DCI can be thought of the DCI calculated with a
dispersal limit of infinity.
In the example below the DCI is calculated for a potamodromous
species. No weighting is used in this case. The global DCI is printed
and a dataframe is returned containing each segment’s contribution to
that global DCI - these are the segmental DCI scores.
yamaska_res <- calculate_dci(net = net,
form = "pot",
pass = "pass_1",
weight = NULL,
threshold = NULL)
#> potamodromous DCI: 55.9457980804545
yamaska_res
#> segment DCI DCI_rel
#> 1 0 12.56427348 22.45793950
#> 2 1 0.67342410 1.20370809
#> 3 2 0.19304043 0.34504902
#> 4 3 2.66030823 4.75515288
#> 5 4 3.03190313 5.41935808
#> 6 5 23.97697527 42.85750868
#> 7 6 0.39179792 0.70031697
#> 8 7 0.12409826 0.22181873
#> 9 8 2.63413959 4.70837789
#> 10 9 1.14102462 2.03951799
#> 11 10 2.78996621 4.98690930
#> 12 11 0.02468512 0.04412328
#> 13 12 2.57142913 4.59628643
#> 14 13 1.79670689 3.21151355
#> 15 14 1.37202572 2.45241960
The calculate_dci function returns a table where each
river segment is assigned a raw DCI score along with a standardized
relative DCI score. The raw scores are useful when comparisons are
desired between different watersheds. On the other hand, the relative
scores allow for a clearer picture of differences within a watershed. An
additional function, export_dci, is provided to rejoin
these results to the original line or point data. This is explained in
the next section.
Distance thresholds
The threshold parameter of the
calculate_dci function allows for the inclusion of a
distance threshold above which river fragments aren’t considered
connected. In terms of the DCI equation, this would mean that not all
pairs of segments will participate in the summation. In addition, the
threshold distance is used to generate an additional weighting term
based on the amount of habitat in the source fragment that can be
considered connected to the exit node. This is done by generating an
isodistance map around the source fragment`s exit node and gathering the
total distance of rivers from the source fragment within that
isodistance map. This length is then divided by the total length of the
river network.
The use of a distance threshold can be useful to derive a clearer
picture of connectivity when a species` dispersal limit is known. Below
is an example done with the same Yamaska watershed using a distance
threshold of 1.5km. This distance is written in the default map units
which are meters.
yamaska_res_thr <- calculate_dci(net = net,
form = "pot",
pass = "pass_1",
weight = NULL,
threshold = 1500)
#> potamodromous DCI with distance limit of 1500: 20.8405089609362
yamaska_res_thr
#> segment DCI DCI_rel
#> 1 0 3.72060622 17.85276083
#> 2 1 0.12702332 0.60950202
#> 3 2 0.06085524 0.29200457
#> 4 3 0.52654957 2.52656770
#> 5 4 1.12481798 5.39726732
#> 6 5 11.61246233 55.72062734
#> 7 6 0.19000695 0.91171932
#> 8 7 0.02638980 0.12662744
#> 9 8 0.66217202 3.17733132
#> 10 9 0.45759080 2.19567957
#> 11 10 0.70825539 3.39845535
#> 12 11 0.00516824 0.02479901
#> 13 12 0.73075801 3.50643072
#> 14 13 0.52481637 2.51825124
#> 15 14 0.36303672 1.74197624
River weights
An additional option with the DCI calculation is to include a
weighting term with the river data. By default the DCI only considers
the length of rivers as a measure of amount of habitat. However, certain
cases might benefit from a more specific definition of amount of
habitat. The weight parameter of the
calculate_dci function can use one of the variables
included with the river data as a weighting term applied to the river
lengths. These values are rescaled between 0 and 1 before applying the
weighting. An example using a habitat suitability weighting is
demonstrated below.
yamaska_res_w <- calculate_dci(net = net,
form = "pot",
pass = "pass_1",
weight = "qual",
threshold = NULL)
#> potamodromous DCI: 55.8242406562483
yamaska_res_w
#> segment DCI DCI_rel
#> 1 0 13.03629720 23.35239503
#> 2 1 1.34188653 2.40377033
#> 3 2 0.22867504 0.40963395
#> 4 3 3.23684297 5.79827497
#> 5 4 3.97779632 7.12557175
#> 6 5 22.64545088 40.56562276
#> 7 6 0.51185358 0.91690199
#> 8 7 0.13546608 0.24266533
#> 9 8 3.13011120 5.60708245
#> 10 9 0.98917111 1.77193832
#> 11 10 2.55680391 4.58009617
#> 12 11 0.01010649 0.01810412
#> 13 12 1.93304252 3.46272962
#> 14 13 1.20238210 2.15387094
#> 15 14 0.88835474 1.59134227
Visualizing results
More sophisticated visualization of river networks and the results of
the DCI analysis are beyond the scope of this package. However, some
code to visualize these networks and DCI results with the ggplot2
package is provided below.
# Visualize river network with nodes coloured according to their type
ggplot() +
geom_sf(data = net %>% tidygraph::activate(edges) %>% sf::st_as_sf()) +
geom_sf(data = net %>% tidygraph::activate(nodes) %>% sf::st_as_sf(), aes(col = type))
# Visualize the DCI color coded river network with the barriers
ggplot() +
geom_sf(data = res_riv, aes(col = DCI)) +
geom_sf(data = net %>% tidygraph::activate(nodes) %>% sf::st_as_sf() %>% dplyr::filter(type == "barrier"))
References
Cote, D., Kehler, D. G., Bourne, C., & Wiersma, Y. F. (2009). A
new measure of longitudinal connectivity for stream networks. Landscape
Ecology, 24(1), 101-113.
