-
Notifications
You must be signed in to change notification settings - Fork 25
/
README.Rmd
361 lines (275 loc) · 19.3 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
---
output:
rmarkdown::github_document:
html_preview: no
---
<!--- README.md is generated from README.Rmd. Please edit that file -->
# prioritizr <img src="man/figures/logo.png" align="right" style="width:140px"/>
## Systematic Conservation Prioritization in R
<!-- badges: start -->
[![lifecycle](https://img.shields.io/badge/Lifecycle-stable-brightgreen.svg)](https://lifecycle.r-lib.org/articles/stages.html)
[![R-CMD-check-Ubuntu](https://img.shields.io/github/actions/workflow/status/prioritizr/prioritizr/R-CMD-check-ubuntu.yaml?branch=main&label=Ubuntu)](https://github.com/prioritizr/prioritizr/actions)
[![R-CMD-check-Windows](https://img.shields.io/github/actions/workflow/status/prioritizr/prioritizr/R-CMD-check-windows.yaml?branch=main&label=Windows)](https://github.com/prioritizr/prioritizr/actions)
[![R-CMD-check-macOS](https://img.shields.io/github/actions/workflow/status/prioritizr/prioritizr/R-CMD-check-macos.yaml?branch=main&label=macOS)](https://github.com/prioritizr/prioritizr/actions)
[![Documentation](https://img.shields.io/github/actions/workflow/status/prioritizr/prioritizr/documentation.yaml?branch=main&label=Documentation)](https://github.com/prioritizr/prioritizr/actions)
[![Coverage Status](https://img.shields.io/codecov/c/github/prioritizr/prioritizr?label=Coverage)](https://app.codecov.io/gh/prioritizr/prioritizr/branch/main)
[![CRAN-Status-Badge](http://www.r-pkg.org/badges/version/prioritizr)](https://CRAN.R-project.org/package=prioritizr)
<!-- badges: end -->
```{r, include = FALSE}
devtools::load_all()
library(prioritizrdata)
plot <- function(...) terra::plot(..., cex = 1.5)
```
```{r, include = FALSE}
# set parameters for figure sizes
h <- 4.5
w <- 6.5
ow <- "500"
knitr::opts_chunk$set(
fig.height = h,
fig.width = w,
out.width = ow,
fig.path = "man/figures/README-", fig.align = "center"
)
```
```{r, include = FALSE}
# set up print method
print <- function(x, ...) {
if (inherits(x, "ConservationProblem")) {
prioritizr::knit_print.ConservationProblem(x)
} else if (inherits(x, "OptimizationProblem")) {
prioritizr::knit_print.OptimizationProblem(x)
} else {
base::print(x)
}
}
```
The _prioritizr R_ package uses mixed integer linear programming (MILP) techniques to provide a flexible interface for building and solving conservation planning problems. It supports a broad range of objectives, constraints, and penalties that can be used to custom-tailor conservation planning problems to the specific needs of a conservation planning exercise. Once built, conservation planning problems can be solved using a variety of commercial and open-source exact algorithm solvers. In contrast to the algorithms conventionally used to solve conservation problems, such as heuristics or simulated annealing, the exact algorithms used here are guaranteed to find optimal solutions. Furthermore, conservation problems can be constructed to optimize the spatial allocation of different management actions or zones, meaning that conservation practitioners can identify solutions that benefit multiple stakeholders. Finally, this package has the functionality to read input data formatted for the _Marxan_ conservation planning program, and find much cheaper solutions in a much shorter period of time than _Marxan_.
## Installation
#### Official version
The latest official version of the _prioritizr R_ package can be installed from the [Comprehensive R Archive Network (CRAN)](https://cran.r-project.org/) using the following _R_ code.
```{r, eval = FALSE}
install.packages("prioritizr", repos = "https://cran.rstudio.com/")
```
#### Developmental version
The latest development version can be installed to gain access to new functionality that is not yet present in the latest official version. Please note that the developmental version is more likely to contain coding errors than the official version. To install the developmental version, you can install it directly from the [GitHub online code repository](https://github.com/prioritizr/prioritizr) or from the [R Universe](https://prioritizr.r-universe.dev/prioritizr). In general, we recommend installing the developmental version from the [R Universe](https://prioritizr.r-universe.dev/prioritizr). This is because installation via [R Universe](https://prioritizr.r-universe.dev/prioritizr) does not require any additional software (e.g., [RTools](https://cran.r-project.org/bin/windows/Rtools/) for Windows systems, or [Xcode and gfortran](https://mac.r-project.org/tools/) for macOS systems).
* To install the latest development version from [R Universe](https://prioritizr.r-universe.dev/prioritizr), use the following _R_ code.
```{r, eval = FALSE}
install.packages(
"prioritizr",
repos = c(
"https://prioritizr.r-universe.dev",
"https://cloud.r-project.org"
)
)
```
* To install the latest development version from [GitHub](https://github.com/prioritizr/prioritizr), use the following _R_ code.
```{r, eval = FALSE}
if (!require(remotes)) install.packages("remotes")
remotes::install_github("prioritizr/prioritizr")
```
## Citation
Please cite the _prioritizr R_ package when using it in publications. To cite the latest official version, please use:
```{r, include = FALSE}
# define helpers
is_rel_version <- function(x) {
length(strsplit(as.character(x), ".", fixed = TRUE)[[1]]) == 3
}
package_version_year <- function(x) {
if (x %in% rownames(available.packages())) {
return(live_package_version_year(x))
} else {
return(archive_package_version_year(x))
}
}
live_package_version_year <- function(x) {
# fetch data
d <- paste0("http://cran.r-project.org/web/packages/", x, "/DESCRIPTION")
d <- read.dcf(url(d))
c(
as.character(d[, "Version"]),
as.character(format(as.POSIXct(d[, "Date/Publication"]), "%Y"))
)
}
archive_package_version_year <- function(x) {
# fetch data (copied from remotes R package)
tf <- tempfile(fileext = ".gz")
on.exit(unlink(tf), add = TRUE)
download.file(tf, url = sprintf("%s/src/contrib/Meta/archive.rds", repo))
con <- gzfile(tf, "rb")
on.exit(close(con), add = TRUE)
d <- readRDS(con)[[x]]
# extract version
rel_version <- basename(rownames(d))
rel_version <- sapply(strsplit(rel_version, "_", fixed = TRUE), `[[`, 2)
rel_version <- sapply(strsplit(rel_version, ".tar", fixed = TRUE), `[[`, 1)
rel_version <- max(rel_version)
# extract year
rel_year <- as.character(max(format(d[, "mtime"], "%Y")))
# return year version
c(rel_version, rel_year)
}
# get version numbers and years
dev_version_year <- c(
as.character(packageDescription("prioritizr")$Version),
format(Sys.time(), "%Y")
)
if (is_rel_version(dev_version_year[[1]])) {
rel_version_year <- dev_version_year
} else {
rel_version_year <- package_version_year("prioritizr")
}
```
```{r, echo = FALSE, results = "asis", comment = ""}
cat(paste0("> Hanson JO, Schuster R, Morrell N, Strimas-Mackey M, Edwards BPM, Watts ME, Arcese P, Bennett J, Possingham HP (",rel_version_year[2],"). prioritizr: Systematic Conservation Prioritization in R. R package version ", rel_version_year[1],". Available at https://CRAN.R-project.org/package=prioritizr.\n"))
```
Alternatively, to cite the latest development version, please use:
```{r, echo = FALSE, results = "asis", comment = ""}
cat(paste0("> Hanson JO, Schuster R, Morrell N, Strimas-Mackey M, Edwards BPM, Watts ME, Arcese P, Bennett J, Possingham HP (", dev_version_year[2],"). prioritizr: Systematic Conservation Prioritization in R. R package version ", dev_version_year[1], ". Available at https://github.com/prioritizr/prioritizr.\n"))
```
Additionally, we keep a [record of publications](https://prioritizr.net/articles/publication_record.html) that use the _prioritizr R_ package. If you use this package in any reports or publications, please [file an issue on GitHub](https://github.com/prioritizr/prioritizr/issues/new) so we can add it to the record.
## Usage
Here we provide a short example showing how the _prioritizr R_ package can be used to build and solve conservation problems. Specifically, we will use an example dataset available through the _prioritizrdata R_ package. Additionally, we will use the _terra R_ package to perform raster calculations. To begin with, we will load the packages.
```{r, eval = FALSE}
# load packages
library(prioritizr)
library(prioritizrdata)
library(terra)
```
We will use the Washington dataset in this example. To import the planning unit data, we will use the `get_wa_pu()` function. Although the _prioritizr R_ package can support many different types of planning unit data, here our planning units are represented as a single-layer raster (i.e., `terra::rast()` object). Each cell represents a different planning unit, and cell values denote land acquisition costs. Specifically, there are `r withr::with_options(list(scipen = 1000), as.character(terra::global(!is.na(get_wa_pu()), "sum", na.rm = TRUE)[[1]]))` planning units in total (i.e., cells with non-missing values).
```{r "planning_units"}
# import planning unit data
wa_pu <- get_wa_pu()
# preview data
print(wa_pu)
# plot data
plot(wa_pu, main = "Costs", axes = FALSE)
```
Next, we will use the `get_wa_features()` function to import the conservation feature data. Although the _prioritizr R_ package can support many different types of feature data, here our feature data are represented as a multi-layer raster (i.e., `terra::rast()` object). Each layer describes the spatial distribution of a feature. Here, our feature data correspond to different bird species. To account for migratory patterns, the breeding and non-breeding distributions of species are represented as different features. Specifically, the cell values denote the relative abundance of individuals, with higher values indicating greater abundance.
```{r "features", fig.height = w * 1.2, fig.width = w * 1.5, out.width = "800"}
# import feature data
wa_features <- get_wa_features()
# preview data
print(wa_features)
# plot the first nine features
plot(wa_features[[1:9]], nr = 3, axes = FALSE)
```
Let's make sure that you have a solver installed on your computer. This is important so that you can use optimization algorithms to generate spatial prioritizations. If this is your first time using the _prioritizr R_ package, please install the HiGHS solver using the following _R_ code. Although the HiGHS solver is relatively fast and easy to install, please note that you'll need to install the [Gurobi software suite and the _gurobi_ _R_ package](https://www.gurobi.com/) for best performance (see the [Gurobi Installation Guide](https://prioritizr.net/articles/gurobi_installation_guide.html) for details).
```{r, eval = FALSE}
# if needed, install HiGHS solver
install.packages("highs", repos = "https://cran.rstudio.com/")
```
Now, let's generate a spatial prioritization. To ensure feasibility, we will set a budget. Specifically, the total cost of the prioritization will represent a 5% of the total land value in the study area. Given this budget, we want the prioritization to increase feature representation, as much as possible, so that each feature would, ideally, have 20% of its distribution covered by the prioritization. In this scenario, we can either purchase all of the land inside a given planning unit, or none of the land inside a given planning unit. Thus we will create a new `problem()` that will use a minimum shortfall objective (via `add_min_shortfall_objective()`), with relative targets of 20% (via `add_relative_targets()`), binary decisions (via `add_binary_decisions()`), and specify that we want near-optimal solutions (i.e., 10% from optimality) using the best solver installed on our computer (via `add_default_solver()`).
```{r}
# calculate budget
budget <- terra::global(wa_pu, "sum", na.rm = TRUE)[[1]] * 0.05
# create problem
p1 <-
problem(wa_pu, features = wa_features) %>%
add_min_shortfall_objective(budget) %>%
add_relative_targets(0.2) %>%
add_binary_decisions() %>%
add_default_solver(gap = 0.1, verbose = FALSE)
# print problem
print(p1)
```
After we have built a `problem()`, we can solve it to obtain a solution.
```{r "minimal_solution"}
# solve the problem
s1 <- solve(p1)
# extract the objective
print(attr(s1, "objective"))
# extract time spent solving the problem
print(attr(s1, "runtime"))
# extract state message from the solver
print(attr(s1, "status"))
# plot the solution
plot(s1, main = "Solution", axes = FALSE)
```
After generating a solution, it is important to evaluate it. Here, we will calculate the number of planning units selected by the solution, and the total cost of the solution. We can also check how many representation targets are met by the solution.
```{r "eval_performance"}
# calculate number of selected planning units by solution
eval_n_summary(p1, s1)
# calculate total cost of solution
eval_cost_summary(p1, s1)
# calculate target coverage for the solution
p1_target_coverage <- eval_target_coverage_summary(p1, s1)
print(p1_target_coverage)
# check percentage of the features that have their target met given the solution
print(mean(p1_target_coverage$met) * 100)
```
Although this solution helps meet the representation targets, it does not account for existing protected areas inside the study area. As such, it does not account for the possibility that some features could be partially -- or even fully -- represented by existing protected areas and, in turn, might fail to identify meaningful priorities for new protected areas. To address this issue, we will use the `get_wa_locked_in()` function to import spatial data for protected areas in the study area. We will then add constraints to the `problem()` to ensure they are selected by the solution (via `add_locked_in_constraints()`).
```{r "locked_in_constraints"}
# import locked in data
wa_locked_in <- get_wa_locked_in()
# print data
print(wa_locked_in)
# plot data
plot(wa_locked_in, main = "Existing protected areas", axes = FALSE)
# create new problem with locked in constraints added to it
p2 <-
p1 %>%
add_locked_in_constraints(wa_locked_in)
# solve the problem
s2 <- solve(p2)
# plot the solution
plot(s2, main = "Solution", axes = FALSE)
```
This solution is an improvement over the previous solution. However, there are some places in the study area that are not available for protected area establishment (e.g., due to land tenure). As a consequence, the solution might not be practical for implementation, because it might select some places that are not available for protection. To address this issue, we will use the `get_wa_locked_out()` function to import spatial data describing which planning units are not available for protection. We will then add constraints to the `problem()` to ensure they are not selected by the solution (via `add_locked_out_constraints()`).
```{r "locked_out_constraints"}
# import locked out data
wa_locked_out <- get_wa_locked_out()
# print data
print(wa_locked_out)
# plot data
plot(wa_locked_out, main = "Areas not available for protection", axes = FALSE)
# create new problem with locked out constraints added to it
p3 <-
p2 %>%
add_locked_out_constraints(wa_locked_out)
# solve the problem
s3 <- solve(p3)
# plot the solution
plot(s3, main = "Solution", axes = FALSE)
```
This solution is even better then the previous solution. However, we are not finished yet. The planning units selected by the solution are fairly fragmented. This can cause issues because fragmentation increases management costs and reduces conservation benefits through edge effects. To address this issue, we can further modify the problem by adding penalties that punish overly fragmented solutions (via `add_boundary_penalties()`). Here we will use a penalty factor (i.e., boundary length modifier) of 0.003, and an edge factor of 50% so that planning units that occur on the outer edge of the study area are not overly penalized.
```{r "boundary_penalties"}
# create new problem with boundary penalties added to it
p4 <-
p3 %>%
add_boundary_penalties(penalty = 0.003, edge_factor = 0.5)
# solve the problem
s4 <- solve(p4)
# plot the solution
plot(s4, main = "Solution", axes = FALSE)
```
Now, let's explore which planning units selected by the solution are most important for cost-effectively meeting the targets. To achieve this, we will calculate importance (irreplaceability) scores using the Ferrier method. Although this method produces scores for each feature separately, we will examine the total scores that summarize overall importance across all features.
```{r "importance"}
# calculate importance scores
rc <-
p4 %>%
eval_ferrier_importance(s4)
# print scores
print(rc)
# plot the total importance scores
## note that gray cells are not selected by the prioritization
plot(
rc[["total"]], main = "Importance scores", axes = FALSE,
breaks = c(0, 1e-10, 0.005, 0.01, 0.025),
col = c("#e5e5e5", "#fff7ec", "#fc8d59", "#7f0000")
)
```
This short example demonstrates how the _prioritizr R_ package can be used to build and customize conservation problems, and then solve them to generate solutions. Although we explored just a few different functions for modifying a conservation problem, the package provides many functions for specifying objectives, constraints, penalties, and decision variables, so that you can build and custom-tailor conservation planning problems to suit your planning scenario.
## Learning resources
The [package website](https://prioritizr.net/index.html) contains information on the _prioritizr R_ package. Here you can find [documentation for every function and built-in dataset](https://prioritizr.net/reference/index.html), and [news describing the updates in each package version](https://prioritizr.net/news/index.html). It also contains the following articles and tutorials.
* [**Getting started**](https://prioritizr.net/articles/prioritizr.html): Short tutorial on using the package.
* [**Package overview**](https://prioritizr.net/articles/package_overview.html): Introduction to systematic conservation planning and a comprehensive overview of the package.
* [**Connectivity tutorial**](https://prioritizr.net/articles/connectivity_tutorial.html): Tutorial on incorporating connectivity into prioritizations.
* [**Calibrating trade-offs tutorial**](https://prioritizr.net/articles/calibrating_trade-offs_tutorial.html): Tutorial on running calibration analyses to satisfy multiple criteria.
* [**Management zones tutorial**](https://prioritizr.net/articles/management_zones_tutorial.html): Tutorial on incorporating multiple management zones and actions into prioritizations.
* [**Gurobi installation guide**](https://prioritizr.net/articles/gurobi_installation_guide.html): Instructions for installing the _Gurobi_ optimization suite for generating prioritizations.
* [**Solver benchmarks**](https://prioritizr.net/articles/solver_benchmarks.html): Performance comparison of optimization solvers for generating prioritizations.
* [**Publication record**](https://prioritizr.net/articles/publication_record.html): List of publications that have cited the package.
Additional resources can also be found in [online repositories under the _prioritizr_ organization](https://github.com/prioritizr). These resources include [slides for talks and seminars about the package](https://github.com/prioritizr/teaching). Additionally, workshop materials are available too (e.g., the [Carleton 2023 workshop](https://prioritizr.github.io/workshop/)).
## Getting help
If you have any questions about the _prioritizr R_ package or suggestions for improving it, please [post an issue on the code repository](https://github.com/prioritizr/prioritizr/issues/new).