Tutorial Prioritzr
R Markdown
#Load the packages and data used in this tutorial
# load packages
library(prioritizrdata)## Warning: package 'prioritizrdata' was built under R version 4.1.2## Loading required package: raster## Loading required package: splibrary(prioritizr)## Warning: package 'prioritizr' was built under R version 4.1.2## Loading required package: sf## Warning: package 'sf' was built under R version 4.1.2## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1## Loading required package: proto## Warning: package 'proto' was built under R version 4.1.2## Warning: multiple methods tables found for 'ncell'library(vegan)## Warning: package 'vegan' was built under R version 4.1.2## Loading required package: permute## Warning: package 'permute' was built under R version 4.1.2## Loading required package: lattice## This is vegan 2.5-7library(cluster)## Warning: package 'cluster' was built under R version 4.1.2# load planning unit data
data(tas_pu)
# load feature data
data(tas_features)# print planning unit data
print(tas_pu)## class : SpatialPolygonsDataFrame
## features : 1130
## extent : 298809.6, 613818.8, 5167775, 5502544 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 5
## names : id, cost, status, locked_in, locked_out
## min values : 1, 0.192488262910798, 0, 0, 0
## max values : 1130, 61.9272727272727, 2, 1, 0# plot map of planning unit costs
plot(st_as_sf(tas_pu[, "cost"]), main = "Planning unit costs")
# plot map of planning unit coverage by protected areas
plot(st_as_sf(tas_pu[, "locked_in"]), main = "Protected area coverage")
#the feature data are expressed as a stack of 62 rasters (i.e., a RasterStack object). Each layer in the stack corresponds to one of 62 different vegetation communities.
#print planning unit data
print(tas_features)## class : RasterStack
## dimensions : 398, 359, 142882, 62 (nrow, ncol, ncell, nlayers)
## resolution : 1000, 1000 (x, y)
## extent : 288801.7, 647801.7, 5142976, 5540976 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : tas_features.1, tas_features.2, tas_features.3, tas_features.4, tas_features.5, tas_features.6, tas_features.7, tas_features.8, tas_features.9, tas_features.10, tas_features.11, tas_features.12, tas_features.13, tas_features.14, tas_features.15, ...
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## max values : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...# plot map of the first four vegetation classes
plot(tas_features[[1:4]], main = paste("Feature", 1:4))
#Problem formulation
# build problem
p1 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
add_min_set_objective() %>%
add_boundary_penalties(penalty = 0.005) %>%
add_relative_targets(0.17) %>%
add_locked_in_constraints("locked_in") %>%
add_binary_decisions()
# print the problem
print(p1)## Conservation Problem
## planning units: SpatialPolygonsDataFrame (1130 units)
## cost: min: 0.19249, max: 61.92727
## features: tas_features.1, tas_features.2, tas_features.3, ... (62 features)
## objective: Minimum set objective
## targets: Relative targets [targets (min: 0.17, max: 0.17)]
## decisions: Binary decision
## constraints: <Locked in planning units [257 locked units]>
## penalties: <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (0.005), zones]>
## portfolio: default
## solver: default#Prioritization (The prioritizr R package supports a range of different exact algorithm solvers, including Gurobi, IBM CPLEX, CBC, Rsymphony, and lpsymphony. Although there are benefits and limitations associated with each of these different solvers, they should return similar results. Note that you will need at least one solver installed on your system to generate prioritizations)
#Install the selected solver# solve problem
s1 <- solve(p1)# plot map of prioritization
plot(st_as_sf(s1[, "solution_1"]), main = "Prioritization",
pal = c("grey90", "darkgreen"))
# create column with existing protected areas
tas_pu$pa <- round(tas_pu$locked_in)
# calculate feature representation statistics based on existing protected areas
tc_pa <- eval_target_coverage_summary(p1, tas_pu[, "pa"])
print(tc_pa)## # A tibble: 62 x 9
## feature met total_amount absolute_target absolute_held absolute_shortf~
## <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 tas_features.1 FALSE 33.9 5.77 0.556 5.21
## 2 tas_features.2 FALSE 170. 28.9 13.5 15.4
## 3 tas_features.3 FALSE 24.0 4.08 2.00 2.08
## 4 tas_features.4 FALSE 32.8 5.57 1.37 4.19
## 5 tas_features.5 FALSE 24.8 4.21 0 4.21
## 6 tas_features.6 FALSE 22.0 3.74 0 3.74
## 7 tas_features.7 FALSE 16.4 2.78 0 2.78
## 8 tas_features.8 FALSE 43.0 7.31 5.12 2.19
## 9 tas_features.9 FALSE 388. 66.0 22.4 43.5
## 10 tas_features.10 FALSE 14.5 2.47 0 2.47
## # ... with 52 more rows, and 3 more variables: relative_target <dbl>,
## # relative_held <dbl>, relative_shortfall <dbl># calculate feature representation statistics based on the prioritization
tc_s1 <- eval_target_coverage_summary(p1, s1[, "solution_1"])
print(tc_s1)## # A tibble: 62 x 9
## feature met total_amount absolute_target absolute_held absolute_shortf~
## <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 tas_features.1 TRUE 33.9 5.77 9.97 0
## 2 tas_features.2 TRUE 170. 28.9 45.7 0
## 3 tas_features.3 TRUE 24.0 4.08 8.00 0
## 4 tas_features.4 TRUE 32.8 5.57 11.4 0
## 5 tas_features.5 TRUE 24.8 4.21 4.98 0
## 6 tas_features.6 TRUE 22.0 3.74 4.00 0
## 7 tas_features.7 TRUE 16.4 2.78 4.00 0
## 8 tas_features.8 TRUE 43.0 7.31 8.12 0
## 9 tas_features.9 TRUE 388. 66.0 120. 0
## 10 tas_features.10 TRUE 14.5 2.47 6.01 0
## # ... with 52 more rows, and 3 more variables: relative_target <dbl>,
## # relative_held <dbl>, relative_shortfall <dbl># explore representation by existing protected areas
## calculate number of features adequately represented by existing protected
## areas
sum(tc_pa$met)## [1] 16## summarize representation (values show percent coverage)
summary(tc_pa$relative_held * 100)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 4.845 13.909 16.995 58.488## visualize representation (values show percent coverage)
hist(tc_pa$relative_held * 100,
main = "Feature representation by existing protected areas",
xlim = c(0, 100),
xlab = "Percent coverage of features (%)")
# explore representation by prioritization
## summarize representation (values show percent coverage)
summary(tc_s1$relative_held * 100)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 17.34 23.83 32.49 37.62 47.17 100.00## calculate number of features adequately represented by the prioritization
sum(tc_s1$met)## [1] 62## visualize representation (values show percent coverage)
hist(tc_s1$relative_held * 100,
main = "Feature representation by prioritization",
xlim = c(0, 100),
xlab = "Percent coverage of features (%)")
# calculate irreplaceability
irrep_s1 <- eval_ferrier_importance(p1, s1["solution_1"])
print(irrep_s1)## class : SpatialPolygonsDataFrame
## features : 1130
## extent : 298809.6, 613818.8, 5167775, 5502544 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 63
## names : tas_features.1, tas_features.2, tas_features.3, tas_features.4, tas_features.5, tas_features.6, tas_features.7, tas_features.8, tas_features.9, tas_features.10, tas_features.11, tas_features.12, tas_features.13, tas_features.14, tas_features.15, ...
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## max values : 0.00746295453528512, 0.00576627072867986, 0.00562272142157942, 0.00632967078700159, 0.00234251072270844, 0.00239546652023844, 0.00375553966968299, 0.0038621328449519, 0.00448956454453004, 0.00737063488756787, 0.249023083885638, 0.0049642250918506, 0.00574350763222759, 0.00622666526653614, 0.00648766470610201, ...# manually coerce values for planning units not selected in prioritization
# to NA, so that they are shown in white
irrep_s1$plot_total <- irrep_s1$total
irrep_s1$plot_total[s1$solution_1 < 0.5] <- NA_real_
# plot map of overall importance scores
plot(st_as_sf(irrep_s1[, "plot_total"]), main = "Overall importance")
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.