raptr: Representative and Adequate Prioritisation Toolkit in R

Overview

This vignette illustrates the basic usage of the raptr R package.

To load the raptr R package and learn more about the package, type the following code into R.

# load raptr R package
library(raptr)

# show package overview
?raptr

The raptr R package uses a range of S4 classes to store conservation planning data, parameters, and prioritisations (Table 1).

Table 1: Main classes in the raptr R package

Class	Description
ManualOpts	place-holder class for manually specified solutions
GurobiOpts	stores parameters for solving optimisation problems using Gurobi
RapUnreliableOpts	stores control variables parameters for the unreliable problem formulation
RapReliableOpts	stores control variables for the reliable problem formulation
DemandPoints	stores the coordinates and weights for a given species and attribute space
PlanningUnitPoints	stores the coordinates and ids for planning units that can be used to preserve a given species
AttributeSpace	stores the coordinates for planning units and the demand points for each species
RapData	stores the all the planning unit, species, and attribute space data
RapUnsolved	stores all the data, control variables, and parameters needed to generate prioritisations
RapResults	stores the prioritisations and summary statistics generated after solving a problem
RapSolved	stores the input data and output results

Getting started

This tutorial is designed to provide users with an understanding of how to use the raptr R package to generate and compare solutions. This tutorial uses several additional packages, so first we will run the following code to load them.

# load packages for tutorial
library(parallel)
library(plyr)
library(dplyr)
library(ggplot2)
library(RandomFields)
library(rgeos)

# set seed for reproducibility
set.seed(500)

# set number of threads for computation
threads <- as.integer(max(1, detectCores()-2))

Simulated examples

Data

To investigate the behaviour of the problem, we will generate prioritisations for three simulated species. We will use the unreliable formulation of the problem to understand the basics, and later move onto the reliable formulation. The first species (termed ‘uniform’) will represent a hyper-generalist. This species will inhabit all areas with equal probability. The second species (termed ‘normal’) will represent a species with a single range core. The third species (termed ‘bimodal’) will represent a species with two distinct ecotypes, each with their own range core. To reduce computational time for this example, we will use a 10 \(\times\) 10 grid of square planning units.

# make planning units
sim_pus <- sim.pus(100L)

# simulate species distributions
sim_spp <- lapply(
    c('uniform', 'normal', 'bimodal'),
    sim.species,
    n=1,
    x=sim_pus,
    res=1
)

Let’s see what these species’ distributions look like.

# plot species
plot(
    stack(sim_spp),
    main=c('Uniform species','Normal species','Bimodal species'),
    addfun=function(){lines(sim_pus)},
    nc=3
)

Distribution of three simulated species. Each square represents a planning unit. The colour of each square denotes the probability that individuals from each species occupy it.

Next, we will generate a set of demand points. To understand the effects of probabilities and weights on the demand points, we will generate the demand points in geographic space. These demand points will be the centroids of the planning units. Additionally, we will use the same set of demand points for each species and only vary the weights of the demand points between species. Note that we are only using the same distribution of demand points for different species for teaching purposes. It is strongly recommended to use different demand points for different species in real-world conservation planning exercises. See the case-study section of this tutorial for examples on how to generate suitable demand points.

# generate coordinates for pus/demand points
pu_coords <- gCentroid(sim_pus, byid=TRUE)

# calculate weights
sim_dps <- lapply(
    sim_spp,
    function(x) {
        return(extract(x, pu_coords))
    }
)

# create demand point objects
sim_dps <- lapply(
    sim_dps,
    function(x) {
        return(
            DemandPoints(
                pu_coords@coords,
                weights=c(x)
            )
        )
    }
)

Now, we will construct a RapUnsolved object to store our input data and parameters. This contains all the information to generate prioritisations.

## create RapUnreliableOpts object
# this stores parameters for the unreliable formulation problem (ie. BLM)
sim_ro <- RapUnreliableOpts()

## create RapData object
# create data.frame with species info
species <- data.frame(
    name=c('uniform', 'normal', 'bimodal')
)

## create data.frame with species and space targets
# amount targets at 20% (denoted with target=0)
# space targets at 20% (denoted with target=1)
targets <- expand.grid(
    species=1:3,
    target=0:1,
    proportion=0.2
)

# calculate probability of each species in each pu
pu_probabilities <- calcSpeciesAverageInPus(sim_pus, stack(sim_spp))

## create AttributeSpace object
# this stores the coordinates of the planning units in an attribute space
# and the coordinates and weights of demand points in the space

pu_points <- PlanningUnitPoints(
    coords=pu_coords@coords,
    ids=seq_len(nrow(sim_pus@data))
)

attr_spaces <- AttributeSpaces(
    list(
        AttributeSpace(
            planning.unit.points=pu_points,
            demand.points=sim_dps[[1]],
            species=1L
        ),
        AttributeSpace(
            planning.unit.points=pu_points,
            demand.points=sim_dps[[2]],
            species=2L
        ),
        AttributeSpace(
            planning.unit.points=pu_points,
            demand.points=sim_dps[[3]],
            species=3L
        )
    ),
    name='geographic'
)

# generate boundary data information
boundary <- calcBoundaryData(sim_pus)

## create RapData object
# this stores all the input data for the prioritisation
sim_rd <- RapData(
    sim_pus@data,
    species,
    targets,
    pu_probabilities,
    list(attr_spaces),
    boundary,
    SpatialPolygons2PolySet(sim_pus)
)

## create RapUnsolved object
# this stores all the input data and parameters needed to generate prioritisations
sim_ru <- RapUnsolved(sim_ro, sim_rd)

Single-species prioritisations

Amount-based targets

To investigate the effects of space-based targets, we will generate a prioritisation for each species using only amount-based targets and compare them to prioritisations generated using amount- and space-based targets. To start off, we will generate a prioritisation for the uniform species using amount-based targets. To do this, we will generate a new sim_ru object by subsetting out the data for the uniform species from the sim_ru object. Then, we will update the targets in the new object. Finally, we will solve the object to generate a prioritisation that fulfills the targets for minimal cost.

# create new object with just the uniform species
sim_ru_s1 <- spp.subset(sim_ru, 'uniform')

# update amount targets to 20% and space targets to 0%
sim_ru_s1 <- update(sim_ru_s1, amount.target=0.2, space.target=NA, solve=FALSE)

# solve problem to identify prioritisation
sim_rs_s1_amount <- solve(sim_ru_s1, Threads=threads)

## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Coefficient statistics:
##   Matrix range    [5e-01, 5e-01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+01, 1e+01]
## Found heuristic solution: objective 20
## Presolve removed 1 rows and 100 columns
## Presolve time: 0.00s
## Presolve: All rows and columns removed
## 
## Explored 0 nodes (0 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

## Warning in validityMethod(object): object@space.held contains values less
## than 0, some species are really poorly represented

## show summary
# note the format for this is similar to that used by Marxan
# see ?raptr::summary for details on this table
summary(sim_rs_s1_amount)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    20   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              42               168               10
##   Connectivity_In_Fraction
## 1                0.1909091

# show amount held
amount.held(sim_rs_s1_amount)

##   uniform
## 1     0.2

# show space held
space.held(sim_rs_s1_amount)

##   uniform (Space 1)
## 1        -0.2363636

Now that we have generated a prioritisation, let’s see what it looks like. We can use the spp.plot method to see how the prioritisation overlaps with the uniform species’ distribution. Note that since all planning units have equal probabilities for this species, all planning units have the same fill colour.

# plot the prioritisation and the uniform species' distribution
spp.plot(sim_rs_s1_amount, 1, main='Uniform species')

A prioritisation for the uniformly distributed species generated using amount-based targets (20%). Sqaures represent planning units. Planning units with a green border are selected for prioritisation, and their colour denotes the probability they are inhabited by the species.

The prioritisation for the uniform species appears to be just a random selection of planning units. This behavior is due to the fact that any prioritisation with 20 planning units is optimal. By relying on just amount targets, this solution may preserve a section of the species’ range core, or just focus on the range margin, or some random part of its range–no emphasis is directed towards preserving different parts of the species’ range. This behavior highlights a fundamental limitation of just using amount-based targets. In the absence of additional criteria, conventional reserve selection problems do not contain any additional information to identify the most effective prioritisation.

Now, we will generate a prioritisation for the normally distributed species using amount-based targets. We will use a similar process to what we used for the uniformly distributed species, but for brevity, we will use code to generate solutions immediately after updating the object.

# create new object with just the normal species
sim_ru_s2 <- spp.subset(sim_ru, 'normal')

# update amount targets to 20% and space targets to 0% and solve it
sim_rs_s2_amount <- update(sim_ru_s2, amount.target=0.2, space.target=NA, solve=TRUE, Threads=threads)

## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-02, 7e-01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [7e+00, 7e+00]
## Found heuristic solution: objective 27
## Presolve removed 0 rows and 86 columns
## Presolve time: 0.00s
## Presolved: 1 rows, 14 columns, 14 nonzeros
## Variable types: 0 continuous, 14 integer (0 binary)
## Presolved: 1 rows, 14 columns, 14 nonzeros
## 
## 
## Root relaxation: objective 9.864476e+00, 6 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0    9.86448    0    1   27.00000    9.86448  63.5%     -    0s
## H    0     0                      10.0000000    9.86448  1.36%     -    0s
## 
## Explored 0 nodes (6 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0%

# show summary
summary(sim_rs_s2_amount)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    10   10             10                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              12               192               16
##   Connectivity_In_Fraction
## 1               0.05454545

# show amount held
amount.held(sim_rs_s2_amount)

##      normal
## 1 0.2026153

# show space held
space.held(sim_rs_s2_amount)

##   normal (Space 1)
## 1        0.5909494

Now let’s visualise the prioritisation we made for the normal species.

# plot the prioritisation and the normal species' distribution
spp.plot(sim_rs_s2_amount, 1, main='Normal species')

A prioritisation for the normally distributed species generated using amount-based targets (20%). See Figure 3 caption for conventions.

The amount-based prioritisation for the normal species focuses only on the species’ range core. This prioritisation fails to secure any peripheral parts of the species’ distribution. As a consequence, it may miss out on populations with novel adaptations to environmental conditions along the species’ range margin.

Now, let’s generate an amount-based target for the bimodally distributed species view it.

# create new object with just the bimodal species
sim_ru_s3 <- spp.subset(sim_ru, 'bimodal')

# update amount targets to 20% and space targets to 0% and solve it
sim_rs_s3_amount <- update(sim_ru_s3, amount.target=0.2, space.target=NA, Threads=threads)

## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e-01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [7e+00, 7e+00]
## Found heuristic solution: objective 21
## Presolve removed 0 rows and 75 columns
## Presolve time: 0.00s
## Presolved: 1 rows, 25 columns, 25 nonzeros
## Variable types: 0 continuous, 25 integer (0 binary)
## Presolved: 1 rows, 25 columns, 25 nonzeros
## 
## 
## Root relaxation: objective 7.919039e+00, 18 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0    7.91904    0    1   21.00000    7.91904  62.3%     -    0s
## H    0     0                       8.0000000    7.91904  1.01%     -    0s
## 
## Explored 0 nodes (18 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0%

# plot the prioritisation and the bimodal species' distribution
spp.plot(sim_rs_s3_amount, 1, main='Bimodal species')

A prioritisation for the bimodally distributed species generated using amount-based targets (20%). See Figure 3 caption for conventions.

# show summary
summary(sim_rs_s3_amount)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL     8    8              8                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1               9               197               14
##   Connectivity_In_Fraction
## 1               0.04090909

# show amount held
amount.held(sim_rs_s3_amount)

##     bimodal
## 1 0.2018391

# show space held
space.held(sim_rs_s3_amount)

##   bimodal (Space 1)
## 1         0.1200278

The amount-based prioritisation for the bimodally distributed species only selects planning units in the bottom left corner of the study area. This prioritisation only preserves individuals belonging to one of the two ecotypes. As a consequence, this prioritisation may fail to preserve a representative sample of the genetic variation found inside this species.

Amount-based and space-based targets

Now that we have generated a prioritisation for each species using only amount-based targets, we will generate a prioritisations using both amount-based and space-targets. To do this we will update the space targets in our amount-based prioritisations to 85%, and store the new prioritisations in new objects.

First, let’s do this for the uniform species.

# make new prioritisation
sim_rs_s1_space <- update(sim_rs_s1_amount, amount.target=0.2, space.target=0.85, Threads=threads)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Coefficient statistics:
##   Matrix range    [5e-01, 8e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 1e+02]
## Found heuristic solution: objective 88
## Presolve removed 36 rows and 0 columns
## Presolve time: 3.75s
## Presolved: 10066 rows, 10100 columns, 40104 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Presolved: 10066 rows, 10100 columns, 40104 nonzeros
## 
## Presolve removed 10066 rows and 10100 columns
## 
## Root relaxation: objective 2.000000e+01, 621 iterations, 0.24 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0      20.0000000   20.00000  0.00%     -    4s
## 
## Explored 0 nodes (1086 simplex iterations) in 4.14 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

# show summary
summary(sim_rs_s1_space)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    20   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1               8               143               69
##   Connectivity_In_Fraction
## 1               0.03636364

# show amount held
amount.held(sim_rs_s1_space)

##   uniform
## 1     0.2

# show space held
space.held(sim_rs_s1_space)

##   uniform (Space 1)
## 1         0.9309091

Let’s take a look at the prioritisation for the uniform species with amount-based and space-based targets. Then, let’s compare the solutions for the amount-based prioritisation with the new prioritisation using both amount and space targets.

# plot the prioritisation and the uniform species' distribution
spp.plot(sim_rs_s1_space, 'uniform', main='Uniform species')

A prioritisation for the uniformly distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.

# plot the difference between old and new prioritisations
plot(sim_rs_s1_amount, sim_rs_s1_space, 1, 1, main='Difference between solutions')

Difference between two prioritisations for the uniformly distributed species. Prioritisation \(X\) was generated using just amount-based targets (20%), and prioritisation \(Y\) was generated using an additional space-based target (85%).

Here we can see that by including a space-target, the prioritisation is spread out evenly across the species’ distribution. Unlike the amount-based prioritisation, this prioritisation samples all the different parts of the species’ distribution.

Now, let’s generate a prioritisation for the normally distributed species that considers amount-based and space-based targets. Then, let’s visualise the new prioritisation and compare it to the old amount-based prioritisation.

# make new prioritisation
sim_rs_s2_space <- update(sim_rs_s2_amount, amount.target=0.2, space.target=0.85, Threads=threads)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-02, 4e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 6e+01]
## Found heuristic solution: objective 91
## Presolve removed 220 rows and 0 columns
## Presolve time: 2.76s
## Presolved: 9882 rows, 10100 columns, 41664 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Presolved: 9882 rows, 10100 columns, 41664 nonzeros
## 
## Presolve removed 9882 rows and 10100 columns
## 
## Root relaxation: objective 1.232570e+01, 4135 iterations, 0.81 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0   12.32570    0  290   91.00000   12.32570  86.5%     -    3s
## H    0     0                      13.0000000   12.32570  5.19%     -    3s
## 
## Explored 0 nodes (4369 simplex iterations) in 3.94 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.300000000000e+01, best bound 1.300000000000e+01, gap 0.0%

# show summary
summary(sim_rs_s2_space)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    13   13             13                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1               5               173               42
##   Connectivity_In_Fraction
## 1               0.02272727

# show amount held
amount.held(sim_rs_s2_space)

##      normal
## 1 0.2014166

# show space held
space.held(sim_rs_s2_space)

##   normal (Space 1)
## 1          0.85809

# plot the prioritisation and the normal species' distribution
spp.plot(sim_rs_s2_space, 'normal', main='Normal species')

A prioritisation for the normally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.

# plot the difference between old and new prioritisations
plot(sim_rs_s2_amount, sim_rs_s2_space, 1, 1, main='Difference between solutions')

Difference between two prioritisations for the normally distributed species. See Figure 7 caption for conventions.

We can see by using both amount-based and space-based targets we can obtain a prioritisation that secures both the species’ range core and parts of its range margin. As a consequence, it may capture any novel adaptations found along the species’ range margin–unlike the amount-based prioritisation.

Finally, let’s generate a prioritisation for the bimodal species using amount-based and space-based targets.

# make new prioritisation
sim_rs_s3_space <- update(sim_rs_s3_amount, amount.target=0.2, space.target=0.85, Threads=threads)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 7e+01]
## Found heuristic solution: objective 86
## Presolve removed 589 rows and 15 columns (presolve time = 5s) ...
## Presolve removed 589 rows and 15 columns
## Presolve time: 6.01s
## Presolved: 9513 rows, 10085 columns, 51461 nonzeros
## Variable types: 0 continuous, 10085 integer (10085 binary)
## Presolved: 9513 rows, 10085 columns, 51461 nonzeros
## 
## Presolve removed 9513 rows and 10085 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    1.0000000e+02   0.000000e+00   1.000000e+02      6s
##     6437    8.9807917e+00   0.000000e+00   0.000000e+00      7s
##     6437    8.9807917e+00   0.000000e+00   0.000000e+00      7s
## 
## Root relaxation: objective 8.980792e+00, 6437 iterations, 1.20 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0    8.98079    0   75   86.00000    8.98079  89.6%     -    7s
## H    0     0                      10.0000000    8.98079  10.2%     -    7s
##      0     0    8.99662    0  125   10.00000    8.99662  10.0%     -    8s
## 
## Cutting planes:
##   Gomory: 1
##   Zero half: 1
## 
## Explored 0 nodes (6761 simplex iterations) in 8.19 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0%

# show summary
summary(sim_rs_s3_space)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    10   10             10                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1               6               186               28
##   Connectivity_In_Fraction
## 1               0.02727273

# show amount held
amount.held(sim_rs_s3_space)

##     bimodal
## 1 0.2127869

# show space held
space.held(sim_rs_s3_space)

##   bimodal (Space 1)
## 1         0.8650189

# plot the prioritisation and the bimodal species' distribution
spp.plot(sim_rs_s3_space, 'bimodal', main='Bimodal species')

A prioritisation for the normally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.

# plot the difference between old and new prioritisations
plot(sim_rs_s3_amount, sim_rs_s3_space, 1, 1, main='Difference between solutions')

Difference between two prioritisations for the bimodally distributed species. See Figure 7 caption for conventions.

Earlier we found that the amount-based prioritisation only preserved individuals from a single ecotype, and would have failed to adequately preserve the intra-specific variation for this species. However, here we can see that by including space-based targets, we can develop prioritisations that secure individuals belonging to both ecotypes. This new prioritisation is much more effective at sampling the intra-specific variation for this species.

Overall, these results demonstrate that under the simplest of conditions, the use of space-based targets can improve prioritisations. However, these prioritisations were each generated for a single species. Prioritisations generated using multiple species may do a better job at preserving the intra-specific variation for individuals species by preserving them in different parts of their range. We will investigate this in the next section.

Multi-species prioritisations

Effects of including space-based targets

So far we have generated prioritisations using only a single species at a time. However, real world prioritisations are often generated using multiple species to ensure that they preserve a comprehensive set of biodiversity. Here, we will generate multi-species prioritisations that preserve all three of the simulated species. First, we will generate a prioritisation using amount-based targets that only aims to preserve 20% of the area they occupy. Then, we will generate a prioritisation that also incoperate space-based targets to also preserve 85% of their geographic distribution. We will then compare the two prioritisations.

# make prioritisations
sim_mrs_amount <- update(
    sim_ru,
    amount.target=c(0.2,0.2,0.2),
    space.target=c(0,0,0),
    Threads=threads
)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 30306 rows, 30100 columns and 120000 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 8e+02]
## Found heuristic solution: objective 98
## Presolve removed 0 rows and 0 columns (presolve time = 5s) ...
## Presolve time: 6.67s
## Presolved: 30306 rows, 30100 columns, 120000 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Presolved: 30306 rows, 30100 columns, 120000 nonzeros
## 
## Presolve removed 30306 rows and 30100 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    1.0000000e+02   0.000000e+00   1.000000e+02      8s
##     1545    2.3592178e+01   0.000000e+00   4.512104e+03     10s
##     2021    2.0000000e+01   0.000000e+00   0.000000e+00     11s
##     2021    2.0000000e+01   0.000000e+00   0.000000e+00     11s
## 
## Root relaxation: objective 2.000000e+01, 2021 iterations, 4.43 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0      20.0000000   20.00000  0.00%     -   13s
## 
## Explored 0 nodes (4785 simplex iterations) in 13.04 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

sim_mrs_space <- update(
    sim_ru,
    amount.target=c(0.2,0.2,0.2),
    space.target=c(0.85, 0.85, 0.85),
    Threads=threads
)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 30306 rows, 30100 columns and 120000 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 1e+02]
## Found heuristic solution: objective 100
## Presolve removed 265 rows and 15 columns (presolve time = 5s) ...
## Presolve removed 329 rows and 15 columns (presolve time = 10s) ...
## Presolve removed 360 rows and 15 columns (presolve time = 15s) ...
## Presolve removed 360 rows and 15 columns
## Presolve time: 15.69s
## Presolved: 29946 rows, 30085 columns, 131504 nonzeros
## Variable types: 0 continuous, 30085 integer (30085 binary)
## Presolve removed 67 rows and 1 columns
## Presolved: 29879 rows, 30084 columns, 127194 nonzeros
## 
## Presolve removed 29879 rows and 30084 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    1.0000000e+02   0.000000e+00   1.000000e+02     18s
##     3691    2.0000000e+01   0.000000e+00   0.000000e+00     20s
##     3691    2.0000000e+01   0.000000e+00   0.000000e+00     20s
## 
## Root relaxation: objective 2.000000e+01, 3691 iterations, 4.03 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0   20.00000    0    4  100.00000   20.00000  80.0%     -   20s
## H    0     0                      20.0000000   20.00000  0.00%     -   20s
## 
## Explored 0 nodes (4876 simplex iterations) in 20.20 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

# show summaries
summary(sim_mrs_amount)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    20   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              10               145               65
##   Connectivity_In_Fraction
## 1               0.04545455

summary(sim_mrs_space)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    20   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1               8               142               70
##   Connectivity_In_Fraction
## 1               0.03636364

# show amount held for each prioritisation
amount.held(sim_mrs_amount)

##   uniform    normal   bimodal
## 1     0.2 0.2214677 0.2502513

amount.held(sim_mrs_space)

##   uniform    normal   bimodal
## 1     0.2 0.2262132 0.2213291

# show space held for each prioritisation
space.held(sim_mrs_amount)

##   uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1         0.9145455        0.8898895         0.9234012

space.held(sim_mrs_space)

##   uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1         0.9254545        0.9057832          0.927655

# plot multi-species prioritisation with amount-based targets
plot(sim_mrs_amount, 1, main='Amount-based targets')

A multi-species prioritisation for the uniformly, normally, and bimodally distributed species generated using just amount-based targets (20%). Squares represent planning units. Dark green planning units are selected for preservation.

# plot multi-species prioritisation with amount- and space-based targets
plot(sim_mrs_space, 1, main='Amount and space-based targets')

A multi-species prioritisation for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 12 caption for conventions.

# difference between the two prioritisations
plot(sim_mrs_amount, sim_mrs_space, 1, 1, main='Difference between solutions')

_Difference between two multi-species prioritisations. See Figure 7 caption for conventions.)

Here we can see that the inclusion of space-based targets changes which planning units are selected for prioritisation, but also the number of planning units that are selected. The amount-based prioritisation is comprised of 20 units, and the space-based prioritisation is comprised of 20 units. This result suggests that an adequate and representative prioritisation can be achieved for only a minor increase in cost.

Uncertainty in species’ distributions

The unreliable formulation does not consider the probability that the planning units are occupied by features when calculating how well a given solution secures a representative sample of an attribute space. Thus solutions identified using the unreliable formulation may select regions of an attribute space for a species using planning units that only have a small chance of being inhabited. As a consequence, if the prioritisation is implemented, it may fail to secure regions of an attribute space if individuals do not inhabit these planning units, and ultimately fail to fulfil the space-based targets.

A simple solution to this issue would be to ensure that planning units cannot be assigned to demand points if they have a low probability of occupancy. This can be achieved by setting a probability threshold for planning units, such that planning units with a probability of occupancy below the threshold are effectively set to zero.

# make new prioritisation with probability threshold of 0.1 for each species
sim_mrs_space2 <- solve(
    prob.subset(
        sim_mrs_space,
        species=1:3,
        threshold=c(0.1,0.1,0.1)
    ),
    Threads=threads
)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 30306 rows, 30100 columns and 119974 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e+01]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 1e+02]
## Found heuristic solution: objective 100
## Presolve removed 264 rows and 15 columns (presolve time = 5s) ...
## Presolve removed 270 rows and 15 columns (presolve time = 10s) ...
## Presolve removed 360 rows and 15 columns (presolve time = 16s) ...
## Presolve removed 360 rows and 15 columns
## Presolve time: 19.67s
## Presolved: 29946 rows, 30085 columns, 131478 nonzeros
## Variable types: 0 continuous, 30085 integer (30085 binary)
## Presolve removed 67 rows and 1 columns
## Presolved: 29879 rows, 30084 columns, 127168 nonzeros
## 
## Presolve removed 29879 rows and 30084 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    1.0000000e+02   0.000000e+00   1.000000e+02     22s
##     4242    2.0000000e+01   0.000000e+00   0.000000e+00     25s
##     4242    2.0000000e+01   0.000000e+00   0.000000e+00     25s
## 
## Root relaxation: objective 2.000000e+01, 4242 iterations, 4.87 seconds
## Total elapsed time = 25.08s
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0   20.00000    0  248  100.00000   20.00000  80.0%     -   26s
## H    0     0                      23.0000000   20.00000  13.0%     -   26s
## H    0     0                      21.0000000   20.00000  4.76%     -   27s
## 
## Explored 0 nodes (7683 simplex iterations) in 27.76 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.100000000000e+01, best bound 2.000000000000e+01, gap 4.7619%

# show summary
summary(sim_mrs_space2)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    21   21             21                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              10               139               71
##   Connectivity_In_Fraction
## 1               0.04545455

# plot prioritisation
plot(sim_mrs_space2, 1)

A multi-species prioritisation for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). This priorititisation was generated to be robust against low occupancy probabilities, by preventing planning units with low probabilities from being used to represent demand points. See Figure 12 caption for conventions.

# difference between prioritisations that use and do not use thresholds
plot(sim_mrs_space2, sim_mrs_space, 1, 1, main='Difference between solutions')

Difference between two multi-species prioritisations. See Figure 7 caption for conventions.

But this method requires setting somewhat arbitrary thresholds. A more robust solution to this issue is to actually use the probability that species occupy planning units to generate the prioritisations. This is what the reliable formulation does. First we will try and generate a solution using existing targets and the reliable formulation. To reduce computational time, we will set the maximum backup \(R\)-level to 1.

# make new prioritisation using reliable formulation
sim_mrs_space3 <- try(update(sim_mrs_space, formulation='reliable', max.r.level=1L, Threads=threads))

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 364206 rows, 181900 columns and 3847200 nonzeros
## Coefficient statistics:
##   Matrix range    [8e-04, 1e+02]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [7e-03, 1e+02]
## Presolve removed 338700 rows and 90800 columns (presolve time = 5s) ...
## Presolve removed 338706 rows and 90806 columns
## Presolve time: 7.23s
## 
## Explored 0 nodes (0 simplex iterations) in 10.08 seconds
## Thread count was 1 (of 2 available processors)
## 
## Model is infeasible
## Best objective -, best bound -, gap -

However, this fails. The reason why we cannot generate a prioritisation that fulfills these targets is because even the solution that contains all the planning units is still insufficient when we consider probabilities. The negative maximum targets printed in the error message indicate that planning units have low probabilities of occupancy. To fulfill the targets, we must obtain more planning units with higher probabilities of occupancy. We also could attempt resolving the problem using a higher \(R\)-level. Instead, we will set lower targets and generate solution.

# make new prioritisation using reliable formulation and reduced targets
sim_mrs_space3 <- update(
    sim_mrs_space,
    formulation='reliable',
    max.r.level=1L,
    space.target=-1000,
    Threads=threads
)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 364206 rows, 181900 columns and 3847200 nonzeros
## Coefficient statistics:
##   Matrix range    [8e-04, 1e+02]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [7e-03, 8e+05]
## Presolve removed 333503 rows and 90800 columns (presolve time = 5s) ...
## Presolve removed 333703 rows and 151800 columns (presolve time = 10s) ...
## Presolve removed 333903 rows and 151800 columns (presolve time = 15s) ...
## Presolve removed 333903 rows and 151800 columns (presolve time = 20s) ...
## Presolve removed 333903 rows and 151800 columns
## Presolve time: 20.29s
## Presolved: 30303 rows, 30100 columns, 90300 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Found heuristic solution: objective 27.0000000
## Presolved: 30303 rows, 30100 columns, 90300 nonzeros
## 
## Presolve removed 30303 rows and 30100 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    1.0000000e+02   0.000000e+00   1.000000e+02     25s
##      129    2.0000000e+01   0.000000e+00   0.000000e+00     25s
##      129    2.0000000e+01   0.000000e+00   0.000000e+00     25s
## 
## Root relaxation: objective 2.000000e+01, 129 iterations, 1.37 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0      20.0000000   20.00000  0.00%     -   24s
## 
## Explored 0 nodes (176 simplex iterations) in 24.87 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

## Warning in validityMethod(object): object@space.held contains values less
## than 0, some species are really poorly represented

# show summary
summary(sim_mrs_space3)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL    20   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              21               157               42
##   Connectivity_In_Fraction
## 1               0.09545455

# plot prioritisation
plot(sim_mrs_space3, 1)

A multi-species prioritisation for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). This priorititisation was generated to be robust against low occupancy probabilities, by explicitly using the probability of occupancy data when deriving a solution. See Figure 12 caption for conventions.

# difference between prioritisations based on unreliable and reliable formulation
plot(sim_mrs_space3, sim_mrs_space, 1, 1, main='Difference between solutions')

Difference between two multi-species prioritisations. See Figure 7 caption for conventions.

An additional planning unit was selected using the reliable formulation. The prioritisation based on the unreliable formulation had 20 planning units, but the prioritisation based on the reliable formlation has 20 planning units. This difference occurs because the reliable formulation needs to ensure that all selected planning units with a low chance of being occupied have a suitable backup planning unit. While the reliable formulation can deliver more robust prioritisations, it takes much longer to solve conservation planning problems expressed using this formulation than the unreliable formulation. As a consequence, the reliable formulation is only feasible for particularly small problems, such as those involving few features and less than several hundred planning units.

Fragmentation

Fragmentation is an important consideration in real-world planning situations. Up until now, we haven’t considered the effects of fragmentation on the viability of the prioritisation. As a consequence, our prioritisations have tended to contain planning units without any neighbours. We can use the BLM parameter to penalise fragmented solutions.

Let’s generate a new prioritisation that heavily penalises fragmentation. Here, we will update the sim_mrs_amount object with BLM of 100.

# update prioritisation
sim_mrs_amount_blm <- update(sim_mrs_amount, BLM=100, Threads=threads)

## Warning for adding variables: zero or small (< 1e-13) coefficients, ignored
## Optimize a model with 30666 rows, 30280 columns and 120720 nonzeros
## Coefficient statistics:
##   Matrix range    [7e-03, 9e+01]
##   Objective range [1e+02, 4e+02]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 8e+02]
## Found heuristic solution: objective 4498
## Presolve removed 0 rows and 0 columns (presolve time = 5s) ...
## Presolve time: 5.97s
## Presolved: 30666 rows, 30280 columns, 120720 nonzeros
## Variable types: 0 continuous, 30280 integer (30280 binary)
## Presolved: 30666 rows, 30280 columns, 120720 nonzeros
## 
## Presolve removed 30666 rows and 30280 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    2.2100000e+04   0.000000e+00   4.010000e+04      8s
##     1253    1.8797519e+03   0.000000e+00   1.928851e+05     10s
##     3563    1.0786850e+03   0.000000e+00   1.782194e+05     15s
##     5663    6.0999666e+02   0.000000e+00   1.699796e+05     20s
##     7593    4.6444444e+02   0.000000e+00   0.000000e+00     24s
##     7593    4.6444444e+02   0.000000e+00   0.000000e+00     24s
## 
## Root relaxation: objective 4.644444e+02, 7593 iterations, 18.19 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  464.44444    0 1725 4498.00000  464.44444  89.7%     -   24s
## H    0     0                    1531.0000000  464.44444  69.7%     -   63s
##      0     0  519.36948    0 1447 1531.00000  519.36948  66.1%     -   70s
##      0     0  519.36948    0 1451 1531.00000  519.36948  66.1%     -   76s
##      0     0  556.17021    0 1520 1531.00000  556.17021  63.7%     -   80s
##      0     0  556.17021    0 1527 1531.00000  556.17021  63.7%     -   80s
##      0     0  556.17021    0 1540 1531.00000  556.17021  63.7%     -   82s
##      0     0  556.17021    0 1542 1531.00000  556.17021  63.7%     -   83s
##      0     0  556.17021    0 1484 1531.00000  556.17021  63.7%     -   89s
##      0     0  556.17021    0 1484 1531.00000  556.17021  63.7%     -  102s
##      0     0  556.17021    0 1462 1531.00000  556.17021  63.7%     -  105s
##      0     0  556.17021    0 1474 1531.00000  556.17021  63.7%     -  113s
##      0     0  556.17021    0 1482 1531.00000  556.17021  63.7%     -  123s
##      0     0  556.17021    0 1461 1531.00000  556.17021  63.7%     -  124s
##      0     2  559.43897    0 1450 1531.00000  559.43897  63.5%     -  146s
##      8    10  626.75175    8 1435 1531.00000  561.88348  63.3%  1152  150s
##     13    15  682.53542   13 1038 1531.00000  561.88348  63.3%  1218  155s
## H   14    14                    1339.0000000  561.88348  58.0%  1131  156s
## H   17    17                    1338.0000000  561.88348  58.0%  1154  158s
##     20    22  750.43478   17  728 1338.00000  561.88348  58.0%  1107  160s
## H   22    22                    1233.0000000  561.88348  54.4%  1072  161s
##     24    26  960.00000   21  324 1233.00000  561.88348  54.4%  1235  167s
## H   26    26                    1022.0000000  561.88348  45.0%  1141  169s
##     28    30  977.14286   24  455 1022.00000  561.88348  45.0%  1206  170s
## H   47    31                    1020.0000000  561.88348  44.9%   743  171s
##     56    37  663.33333    8 1222 1020.00000  579.68642  43.2%   702  175s
##     66    47  734.58333   18  926 1020.00000  579.68642  43.2%   664  180s
##     83    57 1008.46154   30  507 1020.00000  579.68642  43.2%   602  185s
##    102    74  961.46341   19  792 1020.00000  586.55462  42.5%   546  190s
##    109    74  709.80211    3 1555 1020.00000  589.58255  42.2%   568  195s
##    112    77  807.40280    6 1680 1020.00000  589.58255  42.2%   602  200s
##    115    80  864.16667    9  934 1020.00000  589.58255  42.2%   640  207s
##    122    82     cutoff   14      1020.00000  589.58255  42.2%   642  211s
##    126    82     cutoff   16      1020.00000  589.58255  42.2%   656  216s
##    133    83  699.03030    5 1384 1020.00000  602.27238  41.0%   660  220s
##    162   102  696.75277    6 1360 1020.00000  602.76644  40.9%   591  225s
##    168   108  804.61538   12  988 1020.00000  602.76644  40.9%   608  230s
##    176   116  967.66701   20  812 1020.00000  602.76644  40.9%   624  235s
##    195   124  921.61290   11 1137 1020.00000  612.51839  39.9%   610  240s
##    227   138  874.09836   16 1251 1020.00000  615.17625  39.7%   553  246s
##    246   142  728.94695    7 1196 1020.00000  620.00000  39.2%   544  252s
##    249   143     cutoff   10      1020.00000  620.00000  39.2%   553  256s
##    254   144  907.71930   10 1239 1020.00000  620.00000  39.2%   561  260s
##    262   144  756.00000    3 1383 1020.00000  623.69942  38.9%   576  265s
##    283   155  730.91983    8 1286 1020.00000  625.26316  38.7%   555  270s
##    287   159  846.07230   12 1025 1020.00000  625.26316  38.7%   571  275s
##    298   164  911.38889   17  654 1020.00000  625.26316  38.7%   576  280s
##    305   163  754.28571    4 1251 1020.00000  628.64971  38.4%   582  285s
##    314   172  888.00000   13  913 1020.00000  628.64971  38.4%   582  290s
##    332   183  883.72927   12 1207 1020.00000  637.32899  37.5%   572  295s
##    339   187  991.42857   15 1165 1020.00000  637.32899  37.5%   574  300s
##    351   192  774.16667   15  976 1020.00000  638.18182  37.4%   573  305s
##    369   197  759.13043    5 1240 1020.00000  640.56333  37.2%   563  542s
##    371   199  796.11940    7 1231 1020.00000  640.56333  37.2%   563  545s
##    387   209  945.00000   18  873 1020.00000  640.56333  37.2%   555  550s
##    399   213  841.35356   16  985 1020.00000  648.57143  36.4%   561  556s
##    404   217 1002.35294   21  628 1020.00000  648.57143  36.4%   569  560s
##    410   216     cutoff   23      1020.00000  648.57143  36.4%   579  566s
##    421   223  917.87234   13  897 1020.00000  650.90909  36.2%   574  571s
##    429   223  769.25373    7 1196 1020.00000  653.33333  35.9%   567  576s
##    437   231  934.28571   15  864 1020.00000  653.33333  35.9%   568  580s
##    449   236  843.52941   12  910 1020.00000  659.09091  35.4%   559  585s
##    470   243  738.72752   12 1068 1020.00000  663.33333  35.0%   550  590s
##    476   249  880.86957   18  989 1020.00000  663.33333  35.0%   559  595s
##    483   254     cutoff   25      1020.00000  663.33333  35.0%   567  600s
##    493   258  885.51724   14 1200 1020.00000  680.46099  33.3%   575  606s
##    508   260  762.85714   13  967 1020.00000  681.11111  33.2%   569  610s
##    511   262  875.00000   11 1461 1020.00000  681.11111  33.2%   568  627s
##    513   263  940.68966   19 1725 1020.00000  681.11111  33.2%   565  682s
##    514   264  934.28571   16 1422 1020.00000  681.11111  33.2%   564  686s
##    516   265  800.00000   10 1722 1020.00000  681.11111  33.2%   562  716s
##    517   266  853.33333   17 1427 1020.00000  681.11111  33.2%   561  720s
##    518   267  864.16667   10 1522 1020.00000  681.11111  33.2%   560  733s
##    519   267  905.71429   12 1499 1020.00000  681.11111  33.2%   559  736s
##    520   268  870.00000    9 1431 1020.00000  681.11111  33.2%   558  772s
##    521   269  720.00000   14 1619 1020.00000  681.11111  33.2%   557  775s
##    522   269  720.00000   15 1521 1020.00000  681.11111  33.2%   556  780s
##    523   270  943.07692   16 1474 1020.00000  681.11111  33.2%   555  785s
##    524   271  886.66667   18 1474 1020.00000  681.11111  33.2%   554  807s
##    525   273  681.11111   15 1397 1020.00000  681.11111  33.2%   659  813s
##    526   274  681.11111   15 1646 1020.00000  681.11111  33.2%   665  815s
##    529   276  827.04791   17 1313 1020.00000  681.11111  33.2%   676  825s
##    533   279  772.37456   19 1106 1020.00000  681.11111  33.2%   681  830s
##    535   280  860.58860   20 1518 1020.00000  681.11111  33.2%   685  836s
##    537   281  899.84599   21 1094 1020.00000  681.11111  33.2%   693  843s
##    538   282  944.04672   21  318 1020.00000  681.11111  33.2%   699  848s
##    539   283  952.55367   22  222 1020.00000  681.11111  33.2%   700  850s
##    542   285  989.89481   23  508 1020.00000  681.11111  33.2%   713  860s
##    548   288  973.96573   26  126 1020.00000  681.11111  33.2%   714  866s
##    567   282  720.38876   17 1347 1020.00000  681.11111  33.2%   699  873s
##    569   284  841.09950   18 1401 1020.00000  681.11111  33.2%   705  878s
##    570   284  795.75934   19 1359 1020.00000  681.11111  33.2%   715  883s
##    571   285  863.47826   19 1227 1020.00000  681.11111  33.2%   717  888s
##    572   286  819.09369   20 1036 1020.00000  681.11111  33.2%   720  894s
##    573   286  922.94118   20 1282 1020.00000  681.11111  33.2%   723  901s
##    575   288  922.32558   21  865 1020.00000  681.11111  33.2%   733  906s
##    576   288  865.05666   22  745 1020.00000  681.11111  33.2%   741  917s
##    578   290  872.66549   23  631 1020.00000  681.11111  33.2%   754  931s
##    579   290  918.46879   23  736 1020.00000  681.11111  33.2%   755  935s
##    582   292  920.00000   25    6 1020.00000  681.11111  33.2%   765  943s
##    583   293  964.44444   25 2112 1020.00000  681.11111  33.2%   771  946s
## *  587   259              27     920.0000000  681.11111  26.0%   768  948s
##    592   252  681.11111   17 1128  920.00000  681.11111  26.0%   763  950s
##    594   253  747.49450   18 1275  920.00000  681.11111  26.0%   778  957s
##    596   255  810.84931   19 1305  920.00000  681.11111  26.0%   782  961s
##    599   257  827.77551   21 1396  920.00000  681.11111  26.0%   791  965s
##    616   250  917.72534   35 1224  920.00000  681.11111  26.0%   783  970s
##    642   243  850.37991   17 1309  920.00000  682.60159  25.8%   763  975s
##    656   240  753.22770   17 1487  920.00000  713.43062  22.5%   758  980s
##    659   242  803.36428   19 1045  920.00000  748.00036  18.7%   762  986s
##    660   242  812.09444   18 1516  920.00000  748.00036  18.7%   771 1268s
##    661   243  846.48572   19  710  920.00000  748.00036  18.7%   775 1275s
##    665   246  896.92308   21  525  920.00000  748.00036  18.7%   789 1282s
##    669   247  873.84615   21  641  920.00000  748.00036  18.7%   795 1287s
##    670   246     cutoff   22       920.00000  748.00036  18.7%   804 1290s
##    686   242  803.13722   20 1047  920.00000  755.60659  17.9%   794 1295s
##    689   242  893.09539   21  980  920.00000  755.60659  17.9%   802 1301s
##    690   241  913.60169   22  669  920.00000  755.60659  17.9%   810 1307s
##    696   241  888.29268   21  849  920.00000  776.25048  15.6%   809 1310s
##    704   239  810.87749   19 1388  920.00000  776.82907  15.6%   812 1321s
##    710   236     cutoff   21       920.00000  776.82907  15.6%   814 1326s
##    714   234     cutoff   22       920.00000  787.04312  14.5%   827 1340s
##    724   230  870.00000   19  623  920.00000  808.50758  12.1%   830 1347s
##    725   231  901.81818   20  425  920.00000  808.50758  12.1%   833 1351s
##    734   226  846.66667   20  521  920.00000  809.78912  12.0%   827 1356s
##    735   225  915.65217   21  422  920.00000  809.78912  12.0%   829 1360s
##    737   225     cutoff   22       920.00000  809.78912  12.0%   835 1368s
##    741   222  839.87886   18  980  920.00000  810.86241  11.9%   834 1370s
##    749   221  913.65079   21 1306  920.00000  812.00000  11.7%   833 1375s
##    762   217  910.24089   21  545  920.00000  814.52202  11.5%   832 1381s
##    779   212  880.47423   19  981  920.00000  819.84188  10.9%   822 1386s
##    786   209  857.93103   20  606  920.00000  827.71537  10.0%   823 1392s
##    787   210  901.81818   21 1285  920.00000  827.71537  10.0%   825 1402s
##    794   208  841.59881   18 1529  920.00000  828.93205  9.90%   825 1405s
## 
## Cutting planes:
##   Gomory: 9
##   Zero half: 53
## 
## Explored 795 nodes (704775 simplex iterations) in 1405.60 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 9.200000000000e+02, best bound 8.290000000000e+02, gap 9.8913%

# show summary of prioritisation
summary(sim_mrs_amount_blm)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL   920   20             20                220
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1              36               171               13
##   Connectivity_In_Fraction
## 1                0.1636364

# show amount held for each prioritisation
amount.held(sim_mrs_amount_blm)

##   uniform    normal   bimodal
## 1     0.2 0.2122704 0.2035081

# show space held for each prioritisation
space.held(sim_mrs_amount_blm)

##   uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1         0.2606061        0.1936169        0.06415046

# plot prioritisation
plot(sim_mrs_amount_blm, 1)

A multi-species prioritisation for the uniformly, normally, and bimodally distributed species generated using only amount-based targets (20%). Additionally, this priorititisation was specified to have high connectivity, by using a high \(BLM\) parameter. See Figure 12 caption for conventions.

# difference between the two prioritisations
plot(sim_mrs_amount_blm, sim_mrs_amount, 1, 1, main='Difference between solutions')

Difference between two multi-species prioritisations. See Figure 7 caption for conventions.

Here we can see that the prioritisation generated using a BLM parameter of 100 is much more clustered than the prioritisation generated using a BLM of 0. In practice, conservation planners will need to try a variety of BLM parameters to find a suitable prioritisation.

Case-study examples

Introduction

Here we will investigate how space-based targets can affect prioritisations using a more realistic dataset. We will generate prioritisations for the four bird species– blue-winged kookaburra, brown-backed honeyeater, brown falcon, pale-headed rosella–in Queensland, Australia. This region contains a broad range of different habitats–such as rainforests, woodlands, and deserts–making it ideal for this tutorial. First, we will generate a typical amount-based prioritisation that aims to capture 20% of the species’ distributions. Then we will generate a prioritisation that also aims to secure populations in representative parts of the species’ distributions in terms of their geographic location and their environmental heterogeneity. To do this we will generate a prioritisation using 20% amount-based targets and 85% space-based targets. Finally, we will compare these prioritisations to Australia’s existing protected network.

Data

Survey data for the species were obtained from BirdLife Australia. The survey data was rarefied using a 100 km\(^2\) grid, wherein the survey with the greatest number of repeat visits in each grid cell was chosen. To model the distribution of each species, environmental data were obtained at survey location (site). Specifically, climatic data (bio1, bio4, bio15, bio16, bio17) and classifications of the vegetation at the site were used. Occupancy-detection models were fit using Stan using manually tuned parameters (adapt deta=0.9, maximum treedepth=20, chains=4 , warmup iterations=1000, total iterations=1500) with five-fold cross-validation. In each replicate, data were partitioned into training and test sites. A full model was fit using quadratic terms for environmental variables in the site-component, and an intercept in the detection-component of the model. The full model was then subject to a step-wise backwards term deletion routine. Terms were retained when their inclusion resulted in a model with a greater area under the curve (AUC) value as calculated using the test data. Maps were generated for each species as an average of the predictions from the best model in each best replicate. To further improve the accuracy of these maps, areas well outside of the species’ known distributions were set to 0. For each species, this was achieved by masking out biogeographic regions where the species was not observed, and regions that did not have a neighbouring planning unit where the species was observed. The maps were then resampled (10 km\(^2\) resolution) and cropped to the study area. The resulting maps are stored in the cs_spp object.

# load data
data(cs_spp)

# remove areas with low occupancy (prob < 0.6)
cs_spp[cs_spp<0.6] <- NA

# plot species' distributions
plot(cs_spp, main=c(
    "Blue-winged kookaburra", "Brown-backed honeyeater", 
    "Brown falcon", "Pale-headed rosella"
))

Distribution map for four Australian bird species. Pixel colours denote probability of occupancy.

Planning units (50km\(^2\) resolution) were generated across Australia, and then clipped to the Queensland state borders and coastline. Note that we are using excessively coarse planning units so that our examples will complete relatively quickly. In a real-world planning exercise, we would use much finer planning units. To compare our prioritisations to Queensland’s existing protected area network, this network was intersected with the planning units. Planning units with more than 50% of their area inside a protected area had their status set to 2 (following conventions in Marxan). Since we do not have cost data, the prioritisations will aim to select the minimum number of planning units required to meet the targets. The planning units are stored in the cs_pu object.

# load data
data(cs_pus)

## plot planning units
# convert SpatialPolygons to PolySet for quick plotting
cs_pus2 <- SpatialPolygons2PolySet(cs_pus)

# create vector of colours for planning units
# + light green: units not already inside reserve
# + yellow: units already inside reserve
cs_pus_cols <- rep('#c7e9c0', nrow(cs_pus@data))
cs_pus_cols[which(cs_pus$status==2)] <- 'yellow'

# set plotting window
par(mar=c(0.1, 0.1, 4.1, 0.1))

# plot polygons
PBSmapping::plotPolys(
    cs_pus2, col=cs_pus_cols, border='gray30', 
    xlab='', ylab='', axes=FALSE, 
    main='Case-study planning units',
    cex=1.8
)

_Planning units for the case-study examples. Polygons denote planning units. Yellow units have more than 50% of their area already in a reserve.)

# reset plotting window
par(mar=c(5.1, 4.1, 4.1, 2.1))

To map the distribution of environmental conditions across the species’ range, 21 bioclimatic layers were obtained. These layers were cropped to Australia and subject to detrended correspondence analysis to produce two new variables. These layers are stored in the cs_space object.

# load data
data(cs_space)

# plot variables
plot(cs_space, main=c('DC1', 'DC2'), legend=FALSE, axes=FALSE)

Broad-scale environmental variation across Australia. The variable DC1 describes the transition from wet and cool to dry and hot conditions. The variable DC2 describes the transition from wet and hot to dry and cool conditions.

Analysis

To simplify the process of formatting data and generating prioritisations, we can use the rap function. First, we will generate an amount-based prioritisation that aims to capture 10% of the rosella’s range. We will use 50 demand points to map the geographic and environmental spaces. Be warned, the examples hereafter can take 5-10 minutes to run.

# make amount-based prioritisation
# and ignore existing protected areas by discarding values in the 
# status (third) column of the attribute table
cs_rs_amount <- rap(
    cs_pus[,-2], cs_spp, cs_space,
    amount.target=0.2, space.target=NA, n.demand.points=50L,
    include.geographic.space=TRUE, formulation='unreliable',
    kernel.method='hypervolume', quantile=0.7, solve=FALSE
)

## Warning in (function (pus, species, spaces = NULL, amount.target = 0.2, :
## argument to pus does not have a 'status' column, creating default with all
## status=0L

## Evaluating probability density...
## Building tree... done.
## Querying tree... 1.0101e-06  0.010102  0.020203  0.030304  0.0404051  0.0505061  0.0606071  0.0707081  0.0808091  0.0909101  0.101011  0.111112  0.121213  0.131314  0.141415  0.151516  0.161617  0.171718  0.181819  0.19192  0.202021  0.212122  0.222223  0.232324  0.242425  0.252526  0.262627  0.272728  0.282829  0.29293  0.303031  0.313132  0.323233  0.333334  0.343435  0.353536  0.363637  0.373738  0.383839  0.39394  0.404041  0.414142  0.424243  0.434344  0.444445  0.454546  0.464647  0.474748  0.484849  0.494951  0.505052  0.515153  0.525254  0.535355  0.545456  0.555557  0.565658  0.575759  0.58586  0.595961  0.606062  0.616163  0.626264  0.636365  0.646466  0.656567  0.666668  0.676769  0.68687  0.696971  0.707072  0.717173  0.727274  0.737375  0.747476  0.757577  0.767678  0.777779  0.78788  0.797981  0.808082  0.818183  0.828284  0.838385  0.848486  0.858587  0.868688  0.878789  0.88889  0.898991  0.909092  0.919193  0.929294  0.939395  0.949496  0.959597  0.969698  0.979799  0.9899  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70
## Evaluating probability density...
## Building tree... done.
## Querying tree... 9.90099e-06  0.0990198  0.19803  0.29704  0.39605  0.495059  0.594069  0.693079  0.792089  0.891099  0.990109  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.68
## Evaluating probability density...
## Building tree... done.
## Querying tree... 2.60213e-07  0.00260239  0.00520453  0.00780666  0.0104088  0.0130109  0.0156131  0.0182152  0.0208173  0.0234195  0.0260216  0.0286237  0.0312259  0.033828  0.0364301  0.0390323  0.0416344  0.0442365  0.0468387  0.0494408  0.0520429  0.0546451  0.0572472  0.0598493  0.0624515  0.0650536  0.0676557  0.0702579  0.07286  0.0754621  0.0780643  0.0806664  0.0832685  0.0858707  0.0884728  0.0910749  0.0936771  0.0962792  0.0988813  0.101483  0.104086  0.106688  0.10929  0.111892  0.114494  0.117096  0.119698  0.122301  0.124903  0.127505  0.130107  0.132709  0.135311  0.137913  0.140515  0.143118  0.14572  0.148322  0.150924  0.153526  0.156128  0.15873  0.161333  0.163935  0.166537  0.169139  0.171741  0.174343  0.176945  0.179547  0.18215  0.184752  0.187354  0.189956  0.192558  0.19516  0.197762  0.200365  0.202967  0.205569  0.208171  0.210773  0.213375  0.215977  0.218579  0.221182  0.223784  0.226386  0.228988  0.23159  0.234192  0.236794  0.239397  0.241999  0.244601  0.247203  0.249805  0.252407  0.255009  0.257612  0.260214  0.262816  0.265418  0.26802  0.270622  0.273224  0.275826  0.278429  0.281031  0.283633  0.286235  0.288837  0.291439  0.294041  0.296644  0.299246  0.301848  0.30445  0.307052  0.309654  0.312256  0.314858  0.317461  0.320063  0.322665  0.325267  0.327869  0.330471  0.333073  0.335676  0.338278  0.34088  0.343482  0.346084  0.348686  0.351288  0.35389  0.356493  0.359095  0.361697  0.364299  0.366901  0.369503  0.372105  0.374708  0.37731  0.379912  0.382514  0.385116  0.387718  0.39032  0.392922  0.395525  0.398127  0.400729  0.403331  0.405933  0.408535  0.411137  0.41374  0.416342  0.418944  0.421546  0.424148  0.42675  0.429352  0.431954  0.434557  0.437159  0.439761  0.442363  0.444965  0.447567  0.450169  0.452772  0.455374  0.457976  0.460578  0.46318  0.465782  0.468384  0.470986  0.473589  0.476191  0.478793  0.481395  0.483997  0.486599  0.489201  0.491804  0.494406  0.497008  0.49961  0.502212  0.504814  0.507416  0.510018  0.512621  0.515223  0.517825  0.520427  0.523029  0.525631  0.528233  0.530836  0.533438  0.53604  0.538642  0.541244  0.543846  0.546448  0.54905  0.551653  0.554255  0.556857  0.559459  0.562061  0.564663  0.567265  0.569868  0.57247  0.575072  0.577674  0.580276  0.582878  0.58548  0.588082  0.590685  0.593287  0.595889  0.598491  0.601093  0.603695  0.606297  0.6089  0.611502  0.614104  0.616706  0.619308  0.62191  0.624512  0.627114  0.629717  0.632319  0.634921  0.637523  0.640125  0.642727  0.645329  0.647932  0.650534  0.653136  0.655738  0.65834  0.660942  0.663544  0.666147  0.668749  0.671351  0.673953  0.676555  0.679157  0.681759  0.684361  0.686964  0.689566  0.692168  0.69477  0.697372  0.699974  0.702576  0.705179  0.707781  0.710383  0.712985  0.715587  0.718189  0.720791  0.723393  0.725996  0.728598  0.7312  0.733802  0.736404  0.739006  0.741608  0.744211  0.746813  0.749415  0.752017  0.754619  0.757221  0.759823  0.762425  0.765028  0.76763  0.770232  0.772834  0.775436  0.778038  0.78064  0.783243  0.785845  0.788447  0.791049  0.793651  0.796253  0.798855  0.801457  0.80406  0.806662  0.809264  0.811866  0.814468  0.81707  0.819672  0.822275  0.824877  0.827479  0.830081  0.832683  0.835285  0.837887  0.840489  0.843092  0.845694  0.848296  0.850898  0.8535  0.856102  0.858704  0.861307  0.863909  0.866511  0.869113  0.871715  0.874317  0.876919  0.879521  0.882124  0.884726  0.887328  0.88993  0.892532  0.895134  0.897736  0.900339  0.902941  0.905543  0.908145  0.910747  0.913349  0.915951  0.918553  0.921156  0.923758  0.92636  0.928962  0.931564  0.934166  0.936768  0.939371  0.941973  0.944575  0.947177  0.949779  0.952381  0.954983  0.957585  0.960188  0.96279  0.965392  0.967994  0.970596  0.973198  0.9758  0.978403  0.981005  0.983607  0.986209  0.988811  0.991413  0.994015  0.996617  0.99922  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70
## Evaluating probability density...
## Building tree... done.
## Querying tree... 1.64204e-06  0.016422  0.0328424  0.0492627  0.0656831  0.0821034  0.0985238  0.114944  0.131365  0.147785  0.164205  0.180626  0.197046  0.213466  0.229887  0.246307  0.262727  0.279148  0.295568  0.311989  0.328409  0.344829  0.36125  0.37767  0.39409  0.410511  0.426931  0.443351  0.459772  0.476192  0.492612  0.509033  0.525453  0.541874  0.558294  0.574714  0.591135  0.607555  0.623975  0.640396  0.656816  0.673236  0.689657  0.706077  0.722498  0.738918  0.755338  0.771759  0.788179  0.804599  0.82102  0.83744  0.85386  0.870281  0.886701  0.903122  0.919542  0.935962  0.952383  0.968803  0.985223  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70
## Evaluating probability density...
## Building tree... done.
## Querying tree... 1.0101e-06  0.010102  0.020203  0.030304  0.0404051  0.0505061  0.0606071  0.0707081  0.0808091  0.0909101  0.101011  0.111112  0.121213  0.131314  0.141415  0.151516  0.161617  0.171718  0.181819  0.19192  0.202021  0.212122  0.222223  0.232324  0.242425  0.252526  0.262627  0.272728  0.282829  0.29293  0.303031  0.313132  0.323233  0.333334  0.343435  0.353536  0.363637  0.373738  0.383839  0.39394  0.404041  0.414142  0.424243  0.434344  0.444445  0.454546  0.464647  0.474748  0.484849  0.494951  0.505052  0.515153  0.525254  0.535355  0.545456  0.555557  0.565658  0.575759  0.58586  0.595961  0.606062  0.616163  0.626264  0.636365  0.646466  0.656567  0.666668  0.676769  0.68687  0.696971  0.707072  0.717173  0.727274  0.737375  0.747476  0.757577  0.767678  0.777779  0.78788  0.797981  0.808082  0.818183  0.828284  0.838385  0.848486  0.858587  0.868688  0.878789  0.88889  0.898991  0.909092  0.919193  0.929294  0.939395  0.949496  0.959597  0.969698  0.979799  0.9899  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70

## Warning in (function (data, repsperpoint = NULL, bandwidth, quantile = 0, : 
##  Dimensions x and y are highly correlated (r=-0.91)
## Consider removing some axes.

## Evaluating probability density...
## Building tree... done.
## Querying tree... 9.90099e-06  0.0990198  0.19803  0.29704  0.39605  0.495059  0.594069  0.693079  0.792089  0.891099  0.990109  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.69
## Evaluating probability density...
## Building tree... done.
## Querying tree... 2.60213e-07  0.00260239  0.00520453  0.00780666  0.0104088  0.0130109  0.0156131  0.0182152  0.0208173  0.0234195  0.0260216  0.0286237  0.0312259  0.033828  0.0364301  0.0390323  0.0416344  0.0442365  0.0468387  0.0494408  0.0520429  0.0546451  0.0572472  0.0598493  0.0624515  0.0650536  0.0676557  0.0702579  0.07286  0.0754621  0.0780643  0.0806664  0.0832685  0.0858707  0.0884728  0.0910749  0.0936771  0.0962792  0.0988813  0.101483  0.104086  0.106688  0.10929  0.111892  0.114494  0.117096  0.119698  0.122301  0.124903  0.127505  0.130107  0.132709  0.135311  0.137913  0.140515  0.143118  0.14572  0.148322  0.150924  0.153526  0.156128  0.15873  0.161333  0.163935  0.166537  0.169139  0.171741  0.174343  0.176945  0.179547  0.18215  0.184752  0.187354  0.189956  0.192558  0.19516  0.197762  0.200365  0.202967  0.205569  0.208171  0.210773  0.213375  0.215977  0.218579  0.221182  0.223784  0.226386  0.228988  0.23159  0.234192  0.236794  0.239397  0.241999  0.244601  0.247203  0.249805  0.252407  0.255009  0.257612  0.260214  0.262816  0.265418  0.26802  0.270622  0.273224  0.275826  0.278429  0.281031  0.283633  0.286235  0.288837  0.291439  0.294041  0.296644  0.299246  0.301848  0.30445  0.307052  0.309654  0.312256  0.314858  0.317461  0.320063  0.322665  0.325267  0.327869  0.330471  0.333073  0.335676  0.338278  0.34088  0.343482  0.346084  0.348686  0.351288  0.35389  0.356493  0.359095  0.361697  0.364299  0.366901  0.369503  0.372105  0.374708  0.37731  0.379912  0.382514  0.385116  0.387718  0.39032  0.392922  0.395525  0.398127  0.400729  0.403331  0.405933  0.408535  0.411137  0.41374  0.416342  0.418944  0.421546  0.424148  0.42675  0.429352  0.431954  0.434557  0.437159  0.439761  0.442363  0.444965  0.447567  0.450169  0.452772  0.455374  0.457976  0.460578  0.46318  0.465782  0.468384  0.470986  0.473589  0.476191  0.478793  0.481395  0.483997  0.486599  0.489201  0.491804  0.494406  0.497008  0.49961  0.502212  0.504814  0.507416  0.510018  0.512621  0.515223  0.517825  0.520427  0.523029  0.525631  0.528233  0.530836  0.533438  0.53604  0.538642  0.541244  0.543846  0.546448  0.54905  0.551653  0.554255  0.556857  0.559459  0.562061  0.564663  0.567265  0.569868  0.57247  0.575072  0.577674  0.580276  0.582878  0.58548  0.588082  0.590685  0.593287  0.595889  0.598491  0.601093  0.603695  0.606297  0.6089  0.611502  0.614104  0.616706  0.619308  0.62191  0.624512  0.627114  0.629717  0.632319  0.634921  0.637523  0.640125  0.642727  0.645329  0.647932  0.650534  0.653136  0.655738  0.65834  0.660942  0.663544  0.666147  0.668749  0.671351  0.673953  0.676555  0.679157  0.681759  0.684361  0.686964  0.689566  0.692168  0.69477  0.697372  0.699974  0.702576  0.705179  0.707781  0.710383  0.712985  0.715587  0.718189  0.720791  0.723393  0.725996  0.728598  0.7312  0.733802  0.736404  0.739006  0.741608  0.744211  0.746813  0.749415  0.752017  0.754619  0.757221  0.759823  0.762425  0.765028  0.76763  0.770232  0.772834  0.775436  0.778038  0.78064  0.783243  0.785845  0.788447  0.791049  0.793651  0.796253  0.798855  0.801457  0.80406  0.806662  0.809264  0.811866  0.814468  0.81707  0.819672  0.822275  0.824877  0.827479  0.830081  0.832683  0.835285  0.837887  0.840489  0.843092  0.845694  0.848296  0.850898  0.8535  0.856102  0.858704  0.861307  0.863909  0.866511  0.869113  0.871715  0.874317  0.876919  0.879521  0.882124  0.884726  0.887328  0.88993  0.892532  0.895134  0.897736  0.900339  0.902941  0.905543  0.908145  0.910747  0.913349  0.915951  0.918553  0.921156  0.923758  0.92636  0.928962  0.931564  0.934166  0.936768  0.939371  0.941973  0.944575  0.947177  0.949779  0.952381  0.954983  0.957585  0.960188  0.96279  0.965392  0.967994  0.970596  0.973198  0.9758  0.978403  0.981005  0.983607  0.986209  0.988811  0.991413  0.994015  0.996617  0.99922  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70

## Warning in (function (data, repsperpoint = NULL, bandwidth, quantile = 0, : 
##  Dimensions x and y are highly correlated (r=-0.86)
## Consider removing some axes.

## Evaluating probability density...
## Building tree... done.
## Querying tree... 1.64204e-06  0.016422  0.0328424  0.0492627  0.0656831  0.0821034  0.0985238  0.114944  0.131365  0.147785  0.164205  0.180626  0.197046  0.213466  0.229887  0.246307  0.262727  0.279148  0.295568  0.311989  0.328409  0.344829  0.36125  0.37767  0.39409  0.410511  0.426931  0.443351  0.459772  0.476192  0.492612  0.509033  0.525453  0.541874  0.558294  0.574714  0.591135  0.607555  0.623975  0.640396  0.656816  0.673236  0.689657  0.706077  0.722498  0.738918  0.755338  0.771759  0.788179  0.804599  0.82102  0.83744  0.85386  0.870281  0.886701  0.903122  0.919542  0.935962  0.952383  0.968803  0.985223  done.
## Finished evaluating probability density.
## Beginning volume calculation... done. 
## Quantile requested: 0.70   obtained: 0.70

# generate prioritisation
cs_rs_amount <- solve(cs_rs_amount, MIPGap=0.1, Threads=threads)

## Optimize a model with 4 rows, 762 columns and 1336 nonzeros
## Coefficient statistics:
##   Matrix range    [4e+02, 2e+04]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+05, 3e+06]
## Found heuristic solution: objective 159
## Presolve time: 0.01s
## Presolved: 4 rows, 762 columns, 1336 nonzeros
## Variable types: 0 continuous, 762 integer (762 binary)
## Presolved: 4 rows, 762 columns, 1336 nonzeros
## 
## 
## Root relaxation: objective 1.354267e+02, 597 iterations, 0.01 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  135.42669    0    4  159.00000  135.42669  14.8%     -    0s
## H    0     0                     136.0000000  135.42669  0.42%     -    0s
## 
## Explored 0 nodes (597 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.360000000000e+02, best bound 1.360000000000e+02, gap 0.0%

# show summary
summary(cs_rs_amount)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL   136  136            136           98882414
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1         9786021          81270163          7826230
##   Connectivity_In_Fraction
## 1               0.09896624

# plot prioritisation
plot(cs_rs_amount, 1)

Multi-species prioritisation generated for four bird species using amount-based targets (20%). See Figure 12 captions for conventions.

We can also see how well the prioritisation secures the species’ distributions in the geographic and environmental attribute spaces.

# plot prioritisation in geographic attribute space
p1 <- space.plot(cs_rs_amount, 1, 2, main='Blue-winged\nkookaburra')
p2 <- space.plot(cs_rs_amount, 2, 2, main='Brown-backed\nhoneyeater')
p3 <- space.plot(cs_rs_amount, 3, 2, main='Brown falcon')
p4 <- space.plot(cs_rs_amount, 4, 2, main='Pale-headed\nrosella')
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2)

Distribution of amount-based prioritisation in the geographic attribute space. Points denote combinations of environmental conditions. Green and grey points represent planning unit selected for and not selected for prioritisation (respectively). Blue points denote demand points, and their size indicates their weighting.

# plot prioritisation in environmental attribute space
p1 <- space.plot(cs_rs_amount, 1, 1, main='Blue-winged\nkookaburra')
p2 <- space.plot(cs_rs_amount, 2, 1, main='Brown-backed\nhoneyeater')
p3 <- space.plot(cs_rs_amount, 3, 1, main='Brown falcon')
p4 <- space.plot(cs_rs_amount, 4, 1, main='Pale-headed\nrosella')
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2)

Distribution of amount-based prioritisation in the environmental attribute space. See Figure 28 caption for conventions.

Next, let’s generate a prioritisation using amount- and space-based targets. This prioritisation will secure 70% of the species distribution in geographic and environmental space.

# make amount- and space-based prioritisation
cs_rs_space <- update(cs_rs_amount, space.target=0.7, Threads=threads)

## Optimize a model with 134012 rows, 134362 columns and 535736 nonzeros
## Coefficient statistics:
##   Matrix range    [1e-05, 2e+04]
##   Objective range [1e+00, 1e+00]
##   Bounds range    [1e+00, 1e+00]
##   RHS range       [1e+00, 3e+06]
## Presolve removed 1748 rows and 1568 columns (presolve time = 5s) ...
## Presolve removed 2342 rows and 1568 columns (presolve time = 10s) ...
## Presolve removed 2564 rows and 1568 columns (presolve time = 15s) ...
## Presolve removed 2614 rows and 1568 columns (presolve time = 20s) ...
## Presolve removed 2861 rows and 1568 columns (presolve time = 25s) ...
## Presolve removed 2937 rows and 1568 columns (presolve time = 30s) ...
## Presolve removed 3135 rows and 1568 columns (presolve time = 35s) ...
## Presolve removed 3176 rows and 1568 columns (presolve time = 40s) ...
## Presolve removed 3215 rows and 1568 columns (presolve time = 45s) ...
## Presolve removed 3249 rows and 1568 columns (presolve time = 50s) ...
## Presolve removed 3297 rows and 1568 columns (presolve time = 55s) ...
## Presolve removed 3314 rows and 1568 columns (presolve time = 60s) ...
## Presolve removed 3357 rows and 1568 columns (presolve time = 65s) ...
## Presolve removed 4366 rows and 1568 columns (presolve time = 70s) ...
## Presolve removed 5220 rows and 1568 columns (presolve time = 75s) ...
## Presolve removed 5583 rows and 1568 columns (presolve time = 80s) ...
## Presolve removed 5773 rows and 1568 columns (presolve time = 85s) ...
## Presolve removed 6009 rows and 1568 columns (presolve time = 90s) ...
## Presolve removed 6009 rows and 1568 columns (presolve time = 95s) ...
## Presolve removed 6044 rows and 1568 columns (presolve time = 100s) ...
## Presolve removed 6099 rows and 1568 columns (presolve time = 105s) ...
## Presolve removed 6113 rows and 1568 columns (presolve time = 110s) ...
## Presolve removed 6123 rows and 1568 columns (presolve time = 115s) ...
## Presolve removed 6123 rows and 1568 columns (presolve time = 120s) ...
## Presolve removed 6123 rows and 1568 columns (presolve time = 125s) ...
## Presolve removed 6123 rows and 1568 columns (presolve time = 130s) ...
## Presolve removed 6157 rows and 1568 columns (presolve time = 135s) ...
## Presolve removed 6161 rows and 1568 columns (presolve time = 140s) ...
## Presolve removed 6169 rows and 1568 columns (presolve time = 145s) ...
## Presolve removed 6805 rows and 1568 columns (presolve time = 150s) ...
## Presolve removed 7161 rows and 1568 columns (presolve time = 155s) ...
## Presolve removed 7173 rows and 1568 columns (presolve time = 160s) ...
## Presolve removed 7173 rows and 1568 columns (presolve time = 165s) ...
## Presolve removed 7178 rows and 1568 columns (presolve time = 170s) ...
## Presolve removed 7208 rows and 1568 columns (presolve time = 175s) ...
## Presolve removed 7208 rows and 1568 columns (presolve time = 180s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 185s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 190s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 198s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 200s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 205s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 210s) ...
## Presolve removed 7262 rows and 1568 columns (presolve time = 215s) ...
## Presolve removed 7262 rows and 1570 columns (presolve time = 220s) ...
## Presolve removed 7262 rows and 1570 columns (presolve time = 225s) ...
## Presolve removed 7262 rows and 1570 columns (presolve time = 230s) ...
## Presolve removed 7262 rows and 1570 columns
## Presolve time: 234.47s
## Presolved: 126750 rows, 132792 columns, 568881 nonzeros
## Variable types: 0 continuous, 132792 integer (132792 binary)
## Found heuristic solution: objective 246.0000000
## Presolve removed 7 rows and 0 columns (presolve time = 5s) ...
## Presolve removed 7 rows and 0 columns
## Presolved: 126743 rows, 132792 columns, 568837 nonzeros
## 
## Presolve removed 126389 rows and 124881 columns
## 
## Root simplex log...
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##        0    0.0000000e+00   6.211282e+00   9.741038e+08    246s
##     7098    1.3780555e+02   0.000000e+00   1.564206e+03    250s
##    13830    1.3677600e+02   0.000000e+00   4.705039e+03    255s
##    18582    1.3580963e+02   0.000000e+00   1.983081e+02    260s
##    23928    1.3566193e+02   0.000000e+00   7.898848e+01   2024s
##    24324    1.3564946e+02   0.000000e+00   1.515519e+02   2026s
##    26502    1.3561485e+02   0.000000e+00   4.040719e+01   2030s
##    30660    1.3555551e+02   0.000000e+00   1.946716e+02   2035s
##    33603    1.3552511e+02   2.271999e-02   0.000000e+00   2041s
##    33937    1.3552511e+02   0.000000e+00   0.000000e+00   2041s
##    33937    1.3552511e+02   0.000000e+00   0.000000e+00   2041s
## 
## Root relaxation: objective 1.355251e+02, 33937 iterations, 1804.15 seconds
## Total elapsed time = 2045.85s
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  135.52511    0  199  246.00000  135.52511  44.9%     - 2046s
## H    0     0                     136.0000000  135.52511  0.35%     - 2048s
## 
## Explored 0 nodes (48208 simplex iterations) in 2048.82 seconds
## Thread count was 1 (of 2 available processors)
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.360000000000e+02, best bound 1.360000000000e+02, gap 0.0%

# show summary
summary(cs_rs_space)

##   Run_Number  Status Score Cost Planning_Units Connectivity_Total
## 1          1 OPTIMAL   136  136            136           98882414
##   Connectivity_In Connectivity_Edge Connectivity_Out
## 1         9738090          81118598          8025726
##   Connectivity_In_Fraction
## 1               0.09848151

# plot prioritisation
plot(cs_rs_space,1)

# plot prioritisation in geographic attribute space
p1 <- space.plot(cs_rs_space, 1, 2, main='Blue-winged\nkookaburra')
p2 <- space.plot(cs_rs_space, 2, 2, main='Brown-backed\nhoneyeater')
p3 <- space.plot(cs_rs_space, 3, 2, main='Brown falcon')
p4 <- space.plot(cs_rs_space, 4, 2, main='Pale-headed\nrosella')
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2)

Distribution of the amount- and space-based prioritisation in the geographic attribute space. See Figure 28 caption for conventions.

# plot prioritisation in environmental attribute space
p1 <- space.plot(cs_rs_space, 1, 1, main='Blue-winged\nkookaburra')
p2 <- space.plot(cs_rs_space, 2, 1, main='Brown-backed\nhoneyeater')
p3 <- space.plot(cs_rs_space, 3, 1, main='Brown falcon')
p4 <- space.plot(cs_rs_space, 4, 1, main='Pale-headed\nrosella')
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=2)

Distribution of the amount- and space-based prioritisation in the environmental attribute space. See Figure 28 caption for conventions.

Let’s compare these prioritisations with Queensland’s existing protected areas system. To do this, we can create update the cs_rs_space with manually specified solutions to create a RapSolved object to represent the Queensland’s reserve network.

# generate vector with Australia's selections
aus_selections <- which(cs_pus$status>0)

# create new object with Australia's network
cs_rs_aus <- update(cs_rs_amount, b=aus_selections)

Now, let’s plot the performance metrics for these prioritisations.

# define standard error function
se=function(x){sd(x,na.rm=TRUE)/sqrt(sum(!is.na(x)))}

# create a table to store the values for the 3 prioritisations
cs_results <- data.frame(
    name=factor(rep(rep(c('Amount-based\nprioritisation', 
        'Amount & space-based\nprioritisation', 'Queensland reserve\nnetwork'),
        each=4),3), levels=c(c('Amount-based\nprioritisation', 
        'Amount & space-based\nprioritisation', 'Queensland reserve\nnetwork'))),
    variable=rep(c('Amount', 'Geographic space', 'Environmental space'), each=12),
    species=colnames(amount.held(cs_rs_amount)),
    value=c(
        amount.held(cs_rs_amount)[1,], amount.held(cs_rs_space)[1,],
            amount.held(cs_rs_aus)[1,],
        space.held(cs_rs_amount, space=2)[1,], space.held(cs_rs_space, space=2)[1,], 
            space.held(cs_rs_aus, space=2)[1,],
        space.held(cs_rs_amount, space=1)[1,], space.held(cs_rs_space, space=1)[1,],
            space.held(cs_rs_aus, space=1)[1,]
    )
) %>% group_by(name,variable) %>%
    summarise(mean=mean(value),se=se(value))

# plot the performance metrics
ggplot(aes(x=variable, y=mean, fill=name), data=cs_results) +
    geom_bar(position=position_dodge(0.9), stat='identity') +
    geom_errorbar(
        aes(ymin=mean-se, ymax=mean+se), position=position_dodge(0.9),
        width=0.2) +
    xlab('Property of species') +
    ylab('Proportion held in\nselected planning units (%)') +
    scale_fill_discrete(name='') +
    theme_classic() +
    theme(legend.position='bottom',legend.direction='horizontal',
        axis.line.x=element_line(),axis.line.y=element_line())

Prioritisations were generated using amount-based targets (20%), and with additional space-based targets (85%). These are compared to the Queensland reserve network. Data represent means and standard errors for the four species in each prioritisation.

We can see that a greater number of planning units is needed to satisfy the space-based targets. The prioritisation generated using just amount-based targets has 136 planning units, and the prioritisations using amount-based and space-based targets has 136 planning units. These results suggest only a few additional planning units can result in significant improvements in the effectivenss of a prioritisation.