Connectivity tutorial


Connectivity is a key consideration in systematic conservation planning (Margules & Pressey 2000; Briers 2002). This is because isolated and fragmented populations are often more vulnerable to extinction (Dixo et al. 2009; Olds et al. 2012; Hodgson et al. 2009). To promote connectivity in prioritizations, a range of different approaches are available (reviewed in Balbar & Metaxas 2019). These approaches can solely on the spatial configuration of a prioritization to enhance structural connectivity (e.g, reducing the spatial fragmentation of a prioritization; Watts et al. 2009). They can also leverage data – such as environmental, river flow, and telemetry data – to generate prioritizations that promote functional connectivity (e.g., Hermoso et al. 2012; Dwyer et al. 2019; Leonard et al. 2017).

The aim of this tutorial is to show how connectivity can be incorporated into prioritizations using the prioritizr R package. Here we will explore various approaches for incorporating connectivity, and see how they alter the spatial configuration of prioritizations. As you will discover, many of these approaches involve setting threshold or penalty values to specify the relative importance of connectivity compared to other criteria (e.g., overall cost). For more information on calibrating these values, please see the Calibrating trade-offs tutorial.


The dataset used in this tutorial was created for the Coastal Douglas-fir Conservation Partnership (CDFCP; Morrell et al. 2017). Although the original dataset covers a much larger area; for brevity, here we focus only on Salt Spring Island, British Columbia. Briefly, Salt Spring Island supports a diverse and globally unique mix of dry forest and savanna habitats. Today, these habitats are critically threatened due to land conversion, invasive species, and altered disturbance regimes. For more information on the data, please refer to the Marxan tool portal and the tool tutorial.

Extent of Coastal Douglas-fir Conservation Partnership Tool area and location of Salt Spring Island

Extent of Coastal Douglas-fir Conservation Partnership Tool area and location of Salt Spring Island

Let’s begin by loading the packages and data for this tutorial. Since this tutorial requires the prioritizrdata R package, please ensure that it is installed. Specifically, two objects underpin the data for this tutorial. The salt_pu object specifies the planning unit data as a raster layer (i.e., RasterLayer object), and the salt_features object contains biodiversity data represented as a multi-band raster stack (i.e., a RasterStack object).

# load packages

# load planning unit data

# load ecological data

Now we will conduct some preliminary processing. Specifically, we will aggregate from the 100 m resolution to the 300 m resolution. This is to reduce the time needed to generate prioritizations in this tutorial. In practice, we generally recommend consider other criteria too – such as the spatial scale that is relevant for decision making and the resolution of available datasets – when deciding on an appropriate scale for planning units.

# aggregate data to coarser resolution
salt_pu <- aggregate(salt_pu, fact = 3)
salt_features <- aggregate(salt_features, fact = 3)

Next, let’s have a look at the salt_pu object. Here each grid cell represents a planning unit, and the grid cell values denote acquisition costs (BC Assessment 2015). To aid with visualization, we will log-transform the values when plotting them on a map.

# print planning unit data
## class      : RasterLayer 
## dimensions : 94, 67, 6298  (nrow, ncol, ncell)
## resolution : 300, 300  (x, y)
## extent     : 454589.9, 474689.9, 5394414, 5422614  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : salt_pu 
## values     : 0.02552, 8290.383  (min, max)
# plot map showing the planning units costs on a log-scale
plot(log(salt_pu), main = "Planning unit costs (log)", axes = FALSE)

Let’s also look at the salt_features object. This object is a stack of raster layers, with each layer corresponding to a different variable that describes a particular aspect of biodiversity. The first four layers correspond to different ecological communities (i.e., Old Forest, Savannah, Wetland, and Shrub communities), and their cell values indicate the probability of encountering a bird species associated a given community. The fifth layer describes the inverse probability of occurrence of human commensal species. In this tutorial, we will use the first four layers as biodiversity features, and the fifth layer to help parameterize connectivity (wherein higher values denote greater connectivity). So, let’s extract the data and visualize them.

# print original data
## class      : RasterBrick 
## dimensions : 94, 67, 6298, 5  (nrow, ncol, ncell, nlayers)
## resolution : 300, 300  (x, y)
## extent     : 454589.9, 474689.9, 5394414, 5422614  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : salt_features.1, salt_features.2, salt_features.3, salt_features.4, salt_features.5 
## min values :       0.3596232,       0.3340984,       0.1420134,       0.4415211,       0.4060181 
## max values :       0.8951858,       0.6440319,       0.5926431,       0.8065580,       0.8969550
# extract connectivity data
salt_con <- salt_features[[nlayers(salt_features)]]

# print connectivity data
## class      : RasterLayer 
## dimensions : 94, 67, 6298  (nrow, ncol, ncell)
## resolution : 300, 300  (x, y)
## extent     : 454589.9, 474689.9, 5394414, 5422614  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : salt_features.5 
## values     : 0.4060181, 0.896955  (min, max)
# plot map showing the connectivity data
plot(salt_con, main = "Connectivity data", axes = FALSE)

# extract ecological communities and use these as features
salt_features <- salt_features[[seq_len(4)]]

# assign names to features
names(salt_features) <- c("Old_Forest", "Savannah", "Wetland", "Shrub")

# print features
## class      : RasterBrick 
## dimensions : 94, 67, 6298, 4  (nrow, ncol, ncell, nlayers)
## resolution : 300, 300  (x, y)
## extent     : 454589.9, 474689.9, 5394414, 5422614  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : Old_Forest,  Savannah,   Wetland,     Shrub 
## min values :  0.3596232, 0.3340984, 0.1420134, 0.4415211 
## max values :  0.8951858, 0.6440319, 0.5926431, 0.8065580
# plot map showing the distribution of the features
plot(salt_features, main = names(salt_features), axes = FALSE)

Baseline problem

In this tutorial, we will explore a few different ways of incorporating connectivity into prioritizations. To enable comparisons among prioritizations based on different approaches, we will first create a baseline problem formulation that we will subsequently customize to incorporate connectivity. Specifically, we will formulate the baseline problem using the minimum set objective. We will use representation targets of 17% – based on Aichi Biodiversity Target 11 – to provide adequate coverage of each ecological community. Additionally, because land properties on Salt Spring Island can either be acquired in their entirety or not at all, we will use binary decision types. This means that planning units are either selected in the solution or not selected in the solution—planning units cannot be partially acquired. Given all these details, let’s formulate the baseline problem.

# create problem
p0 <- problem(salt_pu, salt_features) %>%
      add_min_set_objective() %>%
      add_relative_targets(0.17) %>%
      add_binary_decisions() %>%

# print problem
## Conservation Problem
##   planning units: RasterLayer (2389 units)
##   cost:           min: 0.02552, max: 8290.38315
##   features:       Old_Forest, Savannah, Wetland, Shrub (4 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.17, max: 0.17)]
##   decisions:      Binary decision 
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), node_file_start (-1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (1)]

After formulating the baseline problem, we can solve it to generate a prioritization.

# solve problem
s0 <- solve(p0)
# print solution
## class      : RasterLayer 
## dimensions : 94, 67, 6298  (nrow, ncol, ncell)
## resolution : 300, 300  (x, y)
## extent     : 454589.9, 474689.9, 5394414, 5422614  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : salt_pu 
## values     : 0, 1  (min, max)
# plot solution
  s0, main = "Baseline prioritization", axes = FALSE,
  breaks = c(0, 0.5, 1), col = c("grey70", "darkgreen")

Next, let’s explore some options for incorporating connectivity.

Adding constraints

Let’s explore approaches for promoting connectivity in prioritizations by adding constraints to the baseline problem formulation. These approaches ensure that prioritizations exhibit certain characteristics (e.g., ensure prioritizations form a contiguous reserve; Önal & Briers 2006). This means that, regardless of the optimality gap used to generate a prioritization, the prioritization will always exhibit these characteristics.

Neighbor constraints

Neighbor constraints can be added to ensure that each selected planning unit has a certain number of neighbors surrounding it (using the add_neighbor_constraints() function) (based on Billionnet 2013). The k parameter can be used to specify the required number of neighbors for each selected planning unit. Let’s generate a prioritization by specifying that each planning unit requires at least two neighbors.

Contiguity constraints

Contiguity constraints can be added to ensure that all planning units form a single contiguous reserve (using the add_contiguity_constraints() function) (similar to Önal & Briers 2006). These constraints are extremely complex. As such, they can only be applied to small conservation planning problems and the Gurobi solver is required to solve them in a feasible period of time. Since it would take a long time to generate a near-optimal prioritization for this dataset with contiguity constraints, we will also tell the solver to simply return the first solution that it finds which meets the representation targets and the contiguity constraints.

## Warning in res(x, ...): overwriting previously defined solver

There is also an even more complex version of the contiguity constraints that is available. These constraints – termed feature contiguity constraints (similar to Cerdeira et al. 2010) – can be added to ensure that all of the selected planning units used to the reach representation targets within a prioritization form a contiguous network for each feature (using the add_feature_contiguity_constraints() function). In other words, they ensure that each feature can disperse through the prioritization to access a target threshold amount of habitat. However, these constraints are extraordinarily complex, only feasible for small problems, and require preprocessing routines to identify initial solutions. As such, we will not consider them in this tutorial.

Linear constraints

Linear constraints can be used to specify that the prioritizations must meet an arbitrary set of criteria. As such, they can be used to ensure that prioritizations provide adequate coverage of planning units that have facilitate a high level of connectivity. Recall that the salt_con data are used to describe connectivity across the study area. Since higher values denote planning units with greater connectivity, we could use linear constraints to ensure that the total sum of connectivity values – based on this dataset – meets a particular threshold (e.g. cover at least 30% of the total amount). This would effectively be treating connectivity as an additional feature (similar to Daigle et al. 2020).

## [1] 520.7628

Although using continuous values has the advantage that the prioritization process can explicitly account for differences in the relative amount of connectivity facilitated by different planning units, the disadvantage is that the prioritization could potentially focus on selecting lots of planning units with low connectivity values. To avoid this result, one strategy is to convert the continuous values into binary values using a threshold limit (similar to Carroll 2021). By applying such a threshold limit, linear constraints can then be used to ensure that the prioritization selects a minimum amount of planning units with high connectivity values (i.e., those with connectivity values that are equal to or greater than the threshold limit).

Another strategy is to clamp the continuous values below a threshold limit are assigned a value of zero (similar to Hanson et al. 2020). This strategy has the advantage that (i) the prioritization won’t focus on selecting lots of planning units with low connectivity values to meet the constraint threshold, and (ii) the optimization process can use semi-continuous values to distinguish between places that can facilitate a moderate amount and a high amount of connectivity.

If we were concerned that the prioritization did not facilitate a high enough level of connectivity, we could increase the threshold value or the threshold_limit value. For example, let’s increase the threshold_limit value used to clamp the continuous connectivity values.

Despite the advantages of clamping the connectivity values, we can see that the prioritization has a relatively high level of spatial fragmentation. In fact, all prioritizations generated using the linear constraints can potentially have this issue. This is because linear constraints do not explicitly account for the spatial arrangement of the planning units. As such, we recommend combining the linear constraints approach with another approach (e.g., the boundary penalties approach discussed below; Carroll 2021).

Adding penalties

Now let’s explore approaches for promoting connectivity in prioritizations by adding penalties to the baseline problem formulation. These approaches involve penalizing solutions according to exhibit certain (e.g., penalize spatial fragmentation of prioritizations; Watts et al. 2009). Unlike constraint-based methods for incorporating connectivity – if the optimality gap used to generate a prioritization is too high – they may not necessarily produce prioritizations that exhibit desirable characteristics.

Boundary penalties

Boundary penalties can be added to used to reduce the spatial fragmentation of prioritizations (using the add_boundary_penalties() function). Specifically, these penalties update the problem formulation to penalize solutions that have a high total amount of exposed boundary length (Ball et al. 2009). Since boundary data often have large values which can degrade solver performance and result in excessive run times (see the Calibrating trade-offs tutorial for details), we will first precompute rescale the boundary data.

Next, let’s generate a prioritization using boundary penalties. To specify the relative importance reducing spatial fragmentation – compared with the primary objective of a problem (e.g. minimizing cost) – we need to a value for the penalty parameter is used. Setting a higher value for penalty indicates that it is more important to avoid highly fragmented solutions. Let’s generate a prioritization with a penalty value of 0.001.

We can see that the resulting prioritization is still relatively fragmented, so let’s try generating another prioritization with a higher penalty value.

Although the prioritization is now less fragmented, it has also selected a greater number of planning units. Let’s calculate the cost of the prioritizations to see how they vary in overall cost.

## # A tibble: 1 × 2
##   summary  cost
##   <chr>   <dbl>
## 1 overall  46.0
## # A tibble: 1 × 2
##   summary  cost
##   <chr>   <dbl>
## 1 overall  46.0
## # A tibble: 1 × 2
##   summary  cost
##   <chr>   <dbl>
## 1 overall  52.3

We can see that the cost of the prioritizations increase with when we use higher penalty values. This is because there is a trade-off between the cost of a prioritization and the level of spatial fragmentation. Although it can be challenging to find the best balance, there are qualitative and quantitative methods available to help navigate such trade-offs. Please see the Calibrating trade-offs tutorial for a details on these methods.

Connectivity penalties

Connectivity penalties can be used to promote connectivity in prioritizations (using the add_connectivity_penalties() function). These penalties use connectivity scores to parametrize the strength of connectivity between pairs of planning units (Beger et al. 2010). Thus higher scores denote a greater level of connectivity between different planning units. For example, previous studies have parametrized connectivity scores using habitat quality, environmental, and river flow data (e.g. Leonard et al. 2017; Alagador et al. 2012; Hermoso et al. 2012). Although there are many approaches to calculate connectivity scores, one approach involves using conductance data – data that describe how much each planning unit facilitates movement (opposite of landscape resistance data) – and calculating scores for each pair of planning units by averaging their conductance values (implemented using the connectivity_matrix() function).

Let’s compute connectivity scores by treating the salt_con object as conductance data. This means that we assume that neighboring planning units with higher values in the salt_con object are capable of facilitating a greater amount of connectivity. Note that the data used to compute connectivity scores must conform to the same spatial properties as the planning unit data (e.g., resolution, spatial extent, coordinate reference system). Also, although we are using raster data here, these scores can also be computed for vector data too (e.g., sf::st_sf() objects). Similar to the boundary data, we will also rescale the connectivity scores to avoid numerical issues during optimization.

After computing the connectivity scores, we can use them to generate prioritizations using connectivity penalties. Similar to the boundary penalties, we use the penalty parameter to specify the relative importance of promoting connectivity relative to the primary objective of a problem (i.e., minimizing overall cost). Let’s generate a prioritization with a penalty value of 0.001.

Now let’s try generating another prioritization with a higher penalty value.

We can see that increasing the penalty parameter causes the prioritizations to select planning units in regions with greater connectivity values (i.e., per the salt_con object). As discussed with the boundary penalties, increasing the penalty value tells the optimization process to focus more on promoting connectivity—meaning that it won’t focus as much on the primary objective (i.e., because the primary objective is to minimize overall costs). For details on calibrating these trade-offs please see the Calibrating trade-offs tutorial. Note that you will need to the use eval_connectivity_summary() function – instead of the eval_boundary_summary() function – when adapting the tutorial code for connectivity penalties.


Hopefully, this tutorial has provided a helpful introduction for incorporating connectivity into prioritizations. Broadly speaking, we recommend using the boundary penalties or the connectivity penalties to ensure that prioritizations explicitly account for the spatial configuration of selected planning units. Additionally, though not fully explored here, the connectivity penalties are a very flexible approach for promoting connectivity. For instance, in addition to parametrizing pair-wise connectivity scores for neighboring planning units, they can also be used to parametrize pair-wise connectivity scores between more distant planning units. Thus connectivity penalties could be used to parametrize connectivity across both small scales and large spatial scales (e.g., using a scaling procedure wherein connectivity scores between pairs of planning units decline with the distance between them).


Alagador, D., Trivino, M., Cerdeira, J.O., Brás, R., Cabeza, M. & Araújo, M.B. (2012). Linking like with like: Optimising connectivity between environmentally-similar habitats. Landscape Ecology, 27, 291–301.

Balbar, A.C. & Metaxas, A. (2019). The current application of ecological connectivity in the design of marine protected areas. Global Ecology and Conservation, 17, e00569.

Ball, I., Possingham, H. & Watts, M.E. (2009). Marxan and relatives: Software for spatial conservation prioritisation. Spatial Conservation Prioritisation: Quantitative Methods & Computational Tools (eds A. Moilanen, K.A. Wilson & H. Possingham), pp. 185–189. Oxford University Press, Oxford, UK.

BC Assessment. (2015). Property Information Services. URL [accessed 13 June 2016].

Beger, M., Linke, S., Watts, M., Game, E., Treml, E., Ball, I. & Possingham, H.P. (2010). Incorporating asymmetric connectivity into spatial decision making for conservation. Conservation Letters, 3, 359–368.

Billionnet, A. (2013). Mathematical optimization ideas for biodiversity conservation. European Journal of Operational Research, 231, 514–534.

Briers, R.A. (2002). Incorporating connectivity into reserve selection procedures. Biological Conservation, 103, 77–83.

Carroll, K.A. (2021). Systematic prioritization protocol applied to wolverine habitat connectivity. STAR Protocols, 2, 100882.

Cerdeira, J.O., Pinto, L.S., Cabeza, M. & Gaston, K.J. (2010). Species specific connectivity in reserve-network design using graphs. Biological Conservation, 143, 408–415.

Daigle, R.M., Metaxas, A., Balbar, A.C., McGowan, J., Treml, E.A., Kuempel, C.D., Possingham, H.P. & Beger, M. (2020). Operationalizing ecological connectivity in spatial conservation planning with Marxan Connect. Methods in Ecology and Evolution, 11, 570–579.

Dixo, M., Metzger, J.P., Morgante, J.S. & Zamudio, K.R. (2009). Habitat fragmentation reduces genetic diversity and connectivity among toad populations in the brazilian atlantic coastal forest. Biological Conservation, 142, 1560–1569.

Dwyer, R.G., Campbell, H.A., Pillans, R.D., Watts, M.E., Lyon, B.J., Guru, S.M., Dinh, M.N., Possingham, H.P. & Franklin, C.E. (2019). Using individual-based movement information to identify spatial conservation priorities for mobile species. Conservation Biology, 33, 1426–1437.

Hanson, J.O., Marques, A., Verı́ssimo, A., Camacho-Sanchez, M., Velo-Antón, G., Martı́nez-Solano & Carvalho, S.B. (2020). Conservation planning for adaptive and neutral evolutionary processes. Journal of Applied Ecology, 57, 2159–2169.

Hermoso, V., Kennard, M.J. & Linke, S. (2012). Integrating multidirectional connectivity requirements in systematic conservation planning for freshwater systems. Diversity and Distributions, 18, 448–458.

Hodgson, J.A., Thomas, C.D., Wintle, B.A. & Moilanen, A. (2009). Climate change, connectivity and conservation decision making: Back to basics. Journal of Applied Ecology, 46, 964–969.

Leonard, P.B., Baldwin, R.F. & Hanks, R.D. (2017). Landscape-scale conservation design across biotic realms: Sequential integration of aquatic and terrestrial landscapes. Scientific Reports, 7, 1–12.

Margules, C.R. & Pressey, R.L. (2000). Systematic conservation planning. Nature, 405, 243–253.

Morrell, N., Schuster, R., Crombie, M. & Arcese, P. (2017). A Prioritization Tool for the Conservation of Coastal Douglas-fir Forest and Savannah Habitats of the Georgia Basin. The Nature Trust of British Colombia, Coastal Douglas Fir Conservation Partnership, and the Department of Forest and Conservation Sciences, University of British Colombia, URL [accessed 9 October 2017].

Olds, A.D., Connolly, R.M., Pitt, K.A. & Maxwell, P.S. (2012). Habitat connectivity improves reserve performance. Conservation Letters, 5, 56–63.

Önal, H. & Briers, R.A. (2006). Optimal selection of a connected reserve network. Operations Research, 54, 379–388.

Watts, M.E., Ball, I.R., Stewart, R.S., Klein, C.J., Wilson, K., Steinback, C., Lourival, R., Kircher, L. & Possingham, H.P. (2009). Marxan with Zones: Software for optimal conservation based land- and sea-use zoning. Environmental Modelling & Software, 24, 1513–1521.