# load Pkg library
using Pkg
# install ConScape
Pkg.add("ConScape")Getting started
Here you will learn
The basics of ConScape: importing data, running basic analyses, and exporting resulting maps.
This tutorial is very similar to the notebook in Appendix A from van Moorter et al. (2022). For a broad overview of the ConScape library, please refer to van Moorter et al. (2022).
In this first notebook we demonstrate the basic workflow as presented in van Moorter et al. (2022) to compute the amount of connected habitat and the movement flow in four steps:
- data import and Grid creation;
- computation of the GridRSP;
- computation of the amount of connected habitat;
- movement flow in two variants (weighted by quality or by proximity).
Setup the environment
Warning
This notebook assumes a running installation of Julia. If you didn’t already install Julia go here.
Install ConScape
In the first time we use ConsCape, we need to install the library. This step can be ignored in the afterwards, unless the user wants to reinstall or update the ConScape library to a new version.
Within the Julia environment, installing ConScape is as simple as:
Load libraries
We continue, and usually would start, by loading the required libraries.
# load libraries
using Pkg
using ConScape
using PlotsThis step is similar to using the library() function in R or the import command in Python.
When setting up the environment, it is also useful to setup the path to the folders where the input data are located and where we want to write the results of our analysis. Here we set the datadir data folder to the folder where the internal ConScape example datasets are saved, after the library is installed.
# path to files
# Pkg.activate(joinpath(ENV["HOMEPATH"], ".julia", "packages", "ConScape", "spkWs", "data"))
# set folders
datadir = joinpath(ENV["HOMEPATH"], ".julia", "packages", "ConScape", "spkWs", "data")
outdir = joinpath(ENV["TMP"], "figures")
# created the output folder, if it does not exist
if !isdir(outdir)
mkdir(outdir)
endStep 1: Data import and grid creation
Import data
We start by importing and checking the input data to be used in ConScape. The first ConScape function is a helper to read maps in ASCII format, the function readasc():
# read habitat quality raster
hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"))
# read movemement probability raster
mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))([NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN], Dict{Any, Any}("cellsize" => 1000.0, "nrows" => 87, "nodata_value" => -9999, "ncols" => 117, "xllcorner" => 110650.0, "yllcorner" => 6.89575e6))The function reads the map as a matrix and the meta data from the ASCII grid as a dictionary. ConScape natively reads ASC files, however, Julia allows easy reading of maps in other file formats through other libraries.
The meta data contain information about the maps: cell size/resolution, number of rows and columns, xy-coordinates of the lower left corner, and no data value:
keys(meta_p)KeySet for a Dict{Any, Any} with 6 entries. Keys:
"cellsize"
"nrows"
"nodata_value"
"ncols"
"xllcorner"
"yllcorner"These meta data can be used to verify that the maps are representing the same geographic domain (i.e. cell size/resolution, number of rows and columns, xy-coordinates of the lower left corner):
collect(values(meta_p))[1:end .!= 3]
collect(values(meta_p))[1:end .!= 3] == collect(values(meta_q))[1:end .!= 3]trueFinally, it is important to ensure that the cells with values match, we conduct the following check and remove non-matching cells:
non_matches = findall(xor.(isnan.(mov_prob), isnan.(hab_qual)))
mov_prob[non_matches] .= NaN
hab_qual[non_matches] .= NaN;
Warning
It is important that all cells or pixels that can be reached in the network (i.e. those with values for the permeability maps) also have matching habitat quality. Otherwise, the ‘NaN’ values may propagate in computations further down. This code ensures this match.
Create a Grid object
Define a ConScape Grid:
adjacency_matrix = ConScape.graph_matrix_from_raster(mov_prob)
g = ConScape.Grid(size(mov_prob)...,
affinities = adjacency_matrix,
source_qualities = hab_qual,
target_qualities = ConScape.sparse(hab_qual),
costs = ConScape.mapnz(x -> -log(x), adjacency_matrix))┌ Info: cost graph contains 4927 strongly connected subgraphs
└ @ ConScape C:\Users\bram.van.moorter\.julia\packages\ConScape\spkWs\src\grid.jl:215
┌ Info: removing 4926 nodes from affinity and cost graphs
└ @ ConScape C:\Users\bram.van.moorter\.julia\packages\ConScape\spkWs\src\grid.jl:225ConScape.Grid of size 87x117
| Affinities |
| Source qualities | Target qualities |
A ConScape Grid describes a graph from a grid of adjacent cells or pixels. It requires four main inputs: the quality of each pixel both as a source and as a target, the affinity between i and j (i.e. probabilities of moving between adjacent pixels i and j), and the cost of moving between between i and j. However, these four inputs can be reduced, for instance, by considering the quality of a pixel identical as a source and target, or by defining the cost as a function of the affinities (e.g. a logarithmic relationship). For our illustration, we introduced those two simplifications and only provided two independent data: the quality of a pixel (identical as source and as target) and the likelihood of moving between adjacent pixels. The likelihood of moving between adjacent pixels i and j was derived from a ‘permeability map’, which describes the permeability of a pixel i (and is similar to the conductivity in circuit theory). The function graph matrix from raster computes the values for an i − j pair from the map either by the average permeability of i and j (AverageWeight) or by the permeability of the target pixel j (TargetWeight; the default); the neighbors of a pixel can be defined either as rook (N4) or as queen (N8; the default).
From the Grid, we can plot the qualities of the pixels:
ConScape.heatmap(g.source_qualities, yflip = true,
title = "Map of habitat uality",
color = cgrad([:white, :green]))
# savefig("figure_grid_outdeg.png")
And the permeability:
ConScape.plot_outdegrees(g, title = "Map of permeability to movement", color = cgrad(:acton))
The information on affinity, cost and quality is stored in the Grid struct as matrices. This is a ‘dense’ matrix (all elements are explicitly stored) in the case of the source qualities:
typeof(g.source_qualities)Matrix{Float64} (alias for Array{Float64, 2})But, a sparse matrix (only the non-zero elements are explicitly stored) for the affinities:
typeof(g.affinities)SparseArrays.SparseMatrixCSC{Float64, Int64}costs:
typeof(g.costmatrix)SparseArrays.SparseMatrixCSC{Float64, Int64}and target qualities:
typeof(g.target_qualities)SparseArrays.SparseMatrixCSC{Float64, Int64}Note the target qualities can be provided either as a sparse (as demonstrated), but also as a dense matrix. This sparse matrix representation forms the basis for the landmark implementation to improve computational performance, as demonstrated in Notebook.
The size of the landscape is given by:
(g.nrows, g.ncols, g.nrows*g.ncols)(87, 117, 10179)Step 2: GridRSP creation
After defining the ConScape Grid of the landscape between adjacent pixels i and j, we need to compute the paths between non-adjacent source and target pixels s and t using the randomized shortest paths framework. This step is computationally intensive and to avoid recomputing central metrics within the framework (e.g. the Z-matrix), we use the ConScape GridRSP struct to store these metrics together with the Grid. In addition to being demanding in terms of processing, this step is also demanding in terms of memory, because sparse affinity and cost matrices are converted into dense matrices. We implemented strategies to reduce memory footprint by using the sparse representation for the target qualities as explained in Van Moorter et al. (2022) and Notebook on performance.
The GridRSP is computed from a Grid with the \(\theta\) parameter to control the amount of randomness in the paths (\(\theta \rightarrow 0\) is random movement, whereas \(\theta \rightarrow \infty\) is optimal movement).
@time h = ConScape.GridRSP(g, θ = 1.0) 10.218829 seconds (1.57 M allocations: 1.952 GiB, 3.51% gc time, 23.50% compilation time)ConScape.GridRSP of size 87x117
As discussed below and in the main text, large values of \(\theta\) may result in numerical instabilities resulting in warning messages.
From this GridRSP we can now easily compute a number of derived metrics.
As discussed in the Notebook on numerical issues and in Van Moorter et al. (2022), large values of \(\theta\) may result in numerical instabilities resulting in warning messages.
From this GridRSP we can now easily compute a number of derived metrics. For instance, we can compute the distance from all pixels in the landscape to a given target pixel.
Let’s take, for illustrative purposes, pixel 4300 as our target:
tmp = zeros(5345)
tmp[4300] = 1
ConScape.plot_values(g, tmp, title = "One target pixel t")The ecological distances from all s to t are:
dists = ConScape.expected_cost(h)
ConScape.plot_values(g, dists[:,4300], title = "Ecological distances to target pixel t")
In many ecological applications we are not simply interested in the ecological distances, but would rather know the proximity between s and t given a movement capacity of the species (e.g. \(\alpha = 75\)):
ConScape.plot_values(g, map!(x -> exp(-x/75), dists[:,4300], dists[:,4300]),
title = "Proximity to target pixel t")
Computation of habitat functionality
One of the main goals of the ConScape library is the computation of the amount of connected habitat in a landscape. The suitability of habitat is typically evaluated based on the local environmental conditions (both biotic and abiotic). We have coined the term ‘functional habitat’ to evaluate not only the suitability, but also the functional connectivity of habitat. In other words, functional habitat is habitat that is simultaneously suitable and functionally connected. The ConScape library was developed primarily to facilitate this task. The functionality of a pixel is computed with the function connected_habitat:
func = ConScape.connected_habitat(h,
connectivity_function = ConScape.expected_cost,
distance_transformation=x -> exp(-x/75));
# func = ConScape.connected_habitat(h, distance_transformation=x -> exp(-x/75));See (Moorter et al. 2021) for a discussion of the expected cost distance and the survival probability. The expected cost distance needs to be transformed to obtain a proximity (\(k_{st} \in [0,1]\)) necessary for the connected habitat computation, an exponential decay is the most common transformation, with a scaling parameter (here: \(\alpha = 75\)).
We plot the results:
ConScape.heatmap(Array(func), yflip = true, title = "Functional habitat")
We can now sum the pixel functionalities to obtain a landscape-level functionality (similar to: (Saura and Pascual-Hortal 2007)):
sum(func)
sum(filter(!isnan, func))590941.1745994696By taking the square root of this number, we can compute the amount of connected habitat similar to (Saura et al. 2011):
sqrt(sum(filter(x -> !isnan(x), func)))768.7269831347601When we compare this value to the amount of ‘unconnected’ habitat, we see that there is some loss of habitat due to movement constraints, i.e.:
100*(1-sqrt(sum(filter(x -> !isnan(x), func)))/
sum(filter(x -> !isnan(x), g.source_qualities)))18.886626540045892percent of the total habitat available.
Computation of movement flow
Finally, in addition to computing habitat functionality or amount of connected habitat, the ConScape library also computes the number of paths going through each pixel weighted by the likelihood of the path, which is a node’s ‘betweenness’ in network sciences. However, ConScape also computes the betweenness weighted by the qualities of the source and target node of a path (i.e. quality weighted) and weighted by the qualities of source and target and their proximity (i.e. proximity weighted). This last version is closely related to what Bodin2010ranking called the ‘generalized betweenness centrality’ using the least-cost paths. A bit tongue in cheek, we could call ConScape’s proximity-weighted betweenness the ‘general generalized betweenness centrality’.
The quality-weighted betweenness (betweenness_qweighted) highlights the pixels that are in between high quality pixels:
ConScape.heatmap(ConScape.betweenness_qweighted(h), yflip = true, title = "Betweenness")
We can also use this function to illustrate the effect of the randomness parameter \(\theta\). By setting the quality of two pixels to one while all others are zero we can visualize the distribution over paths between these two pixels, for instance, two pixels:
tmp = zeros(g.nrows, g.ncols)
tmp[60,70] = 1
tmp[50, 105] = 1
g_tmp = ConScape.Grid(size(mov_prob)...,
affinities=ConScape.graph_matrix_from_raster(mov_prob),
qualities=tmp,
costs=ConScape.MinusLog());
ConScape.heatmap(g_tmp.source_qualities, yflip=true, title = "Source and target pixel")┌ Info: cost graph contains 4927 strongly connected subgraphs
└ @ ConScape C:\Users\bram.van.moorter\.julia\packages\ConScape\spkWs\src\grid.jl:215┌ Info: removing 4926 nodes from affinity and cost graphs
└ @ ConScape C:\Users\bram.van.moorter\.julia\packages\ConScape\spkWs\src\grid.jl:225
We can now over a range of \(\theta\) values show the pixels that are most ‘in-between’ the two focal pixels:
θs = (2.5, 1.0, 0.5, 0.1, 0.01, 0.001)
betw = [copy(mov_prob), copy(mov_prob), copy(mov_prob),
copy(mov_prob), copy(mov_prob), copy(mov_prob)]
for i in 1:length(θs)
h_tmp = ConScape.GridRSP(g_tmp, θ=θs[i]);
betw[i] = ConScape.betweenness_qweighted(h_tmp)
end
We see that as \(\theta\) becomes larger the path distribution becomes narrowly focused on a single path, i.e. the least-cost path, whereas as \(\theta\) becomes smaller the path distribution becomes more diffuse, i.e. a random walk.
Finally, the proximity weighted betweenness (betweenness_kweighted) computes the betweenness where paths are weighted by the quality of the source and target, but also by their proximity (\(k_{st}\)). Here as with the computation of the amount of connected habitat, we need to choose a metric to quantify the proximity (e.g. expected cost or survival probability), and in the case of a distance metric we also need to provide the transformation to proximity (e.g. exponential decay):
kbetw = ConScape.betweenness_kweighted(h,
distance_transformation=x -> exp(-x/75))
ConScape.heatmap(kbetw, yflip=true, title="Betweenness")
Thus, the proximity weighted betweenness highlights the pixels that are most in-between high-quality pixels that are functionally connected together. Hence, this metric is expected to capture best the actual movement flows on the landscape (bodin2010ranking).
For further inspection of the maps produced in ConScape the library also provides a helper export function to export to an ASCII grid, where the meta data can simply be provided as the dictionary from the input maps:
ConScape.writeasc(joinpath(outdir, "Pweighted_flow.asc"), kbetw, meta_p)Summary
We showed a basic workflow to demonstrate the main functionalities in the ConScape library. Step 1 is data import and the representation of the landscape as a Grid with connections between adjacent pixels i and j. Step 2 is the computation of the fundamental matrix in the GridRSP for the RSP distribution over all paths between all source \(s\) and target \(t\) pixels given the randomness of these paths (\(\theta\) parameter). In steps 3 and 4 the main outputs are computed from the : the amount of connected habitat to each pixel and the movement flow through each pixel.
The following notebooks go into some more details of the different components of this basic workflow: Notebook~\(\ref{nbk:landmarks}\) expands on improving computational performance by using a landmark approach, Notebook~\(\ref{nbk:distance}\) discusses in more depth different distance metrics and their transformation to proximities, Notebook~\(\ref{nbk:cost}\) illustrates the use of independent likelihood and cost of movement with mortality risks, and Notebook~\(\ref{nbk:functionality}\) shows the different ways the amount of connected habitat to a pixel can be quantified using .
References
Moorter, Bram van, Ilkka Kivimäki, Manuela Panzacchi, and Marco Saerens. 2021. “Defining and Quantifying Effective Connectivity of Landscapes for Species’ Movements.” Ecography 44 (6): 870–84.
Saura, Santiago, Christine Estreguil, Coralie Mouton, and Mónica Rodrı́guez-Freire. 2011. “Network Analysis to Assess Landscape Connectivity Trends: Application to European Forests (1990–2000).” Ecological Indicators 11 (2): 407–16.
Saura, Santiago, and Lucı́a Pascual-Hortal. 2007. “A New Habitat Availability Index to Integrate Connectivity in Landscape Conservation Planning: Comparison with Existing Indices and Application to a Case Study.” Landscape and Urban Planning 83 (2): 91–103.