Comparing prioritization methods

1. Setup and background

Conservation planning problems have mostly been dividied into two categories. The minimum-set approach aims to minimize the amount of resources needed to achieve all targets, regardless of how the targets are defined (Underhill 1994). The maximum-coverage approach on the other hand aims at maximizing protection for a given budget (Camm et al. 1996). The quantity maximized can be anything from species richness, to productivity, to accessibility by the public, or any combination of factors i.e., whatever is deemed relevant by the analyst (Arponen et al. 2009). Both approaches are related and are useful for different planning circumstances. Furthermore, there are other types of related approaches, such as maximum utility approches (e.g. Moilanen and Arponen 2011 and Laitila and Moilanen 2012).

Many different heuristic and optimization have been developed over the year to solve minimum-set, maximum-coverage, and maximum utility type of problems. While comparing different methods is usually like comparing apples to oranges, especially from practical conservation planning point of view such comparisons can have use. Here, I compare three methods that fall into maximum-coverage category.

1.1 Zonation tutorial data

For testing purposes, we’re going to use the Zonation tutorial data. The dataset is fairly small and useful for the purpose of testing.

Let’s start by reading in the data and inspecting the individual species features (N=7). The values in each raster are in range [0, 1] and indicate the propabilty of occurrence.

# Set up directory path for the tutorial data
z_tutorial_dir <- "~/dev/git-data/zonation-tutorial"
z_tutorial_data_dir <- file.path(z_tutorial_dir, "data")
tutorial_files <- list.files(path = z_tutorial_data_dir, pattern = "species.+\\.tif", 
                             full.names = TRUE)
sp_rasters <- raster::stack(tutorial_files)
# For some mysterious reason, minmax isn't set automatically for sp_rasters[[1]]
# (species1)
sp_rasters[[1]] <- setMinMax(sp_rasters[[1]])
# Generate template for other spatial data later on
nrows <- nrow(sp_rasters)
ncols <- ncol(sp_rasters)
template <- extent(sp_rasters) %>% 
              raster(nrows = nrows, ncols = ncols, vals = 1)
template[is.na(sp_rasters[[1]])] <- NA
levelplot(sp_rasters, margin = FALSE, scales = list(draw = FALSE),
          col.regions = viridis, layout = c(3, 3))

1.2 Generate cost data

While we will not be using cost data as in actually considering cost data, but we need for as a constraint data for the Gurobi ILP optimization. The cost is 1.0 in each cell, basically meaning that it translates into number of cells. This way it is possible to link an aerial target (such as the top 10% of the landscape) into a budgetary constraint.

cost <- extent(sp_rasters) %>% 
  raster(nrows = nrows, ncols = ncols, vals = 1)
# Alternatively, cost data can be simulated
# cost <- gaussian_field(r, 20, mean = 1000, variance = 500) #%>% 
#  setNames("cost")
# Fill areas normalle NA (NoData) with cost of 0
cost[is.na(sp_rasters[[1]])] <- 0
levelplot(cost, main = "Cost", margin = FALSE, scales = list(draw = FALSE),
          col.regions = viridis)

2. Prioritization methods

Here, we’re testing three different prioritization methods against each other at increasing level of complexity:

1. Rarity-wighted richness
2. Zonation (two methods: core-area Zonation and additive benefit function)
3. Integer programming (ILP) using Gurobi solver

2.1 Rarity-Weighted Richness (RWR)

From Albuquerque and Beier (2015):

“Williams et al. (1996) proposed that the rarity value of a species can be characterized by the inverse of the number of sites or planning units in which it occurs. Thus if a species is found in only 1 site, the species would have the maximum rarity score of 1/1 = 1, and a species that occurs in 20 sites would have a rarity score of 1/20 = 0.05. Williams et al. also proposed that the rarity scores of all species in the site can be summed to yield a single RWR value for the site:

\[ \sum_{n}^{1} (1 / c) \]

where \(c_{i}\) is the number of sites occupied by species \(i\), and the values are summed for the \(n\) species that occur in that site."

While the formulation by Albuquerque and Beier (2015) applies to binary presence/absence data, it can be extended to data where the values in each cell is a quantitative measure, e.g. the probability of occurrence. Hence, the formula becomes:

\[ \sum_{n=1}^{n} \frac{o_{ij}}{\sum_{i} o_{i}} \]

where \(o_{ij}\) is the occurence level of feature \(i\) in cell \(j\) and \(\sum_{i} o_{i}\) is the sum of the occurrence levels of all cells of feature \(i\). Below, after ranking the values are rescaled to range [0, 1]. Values close to 0 are the lowest priority locations, values close to 1 the highest.

# Define a function for occurrence level normalization
ol_normalize <- function(x) {
  min <- raster::minValue(x)
  return((x - min) / raster::cellStats(x - min, "sum"))
}
# Occurrence level normalize all species rasters
sp_rasters_normalized <- ol_normalize(sp_rasters)
# Sum up all normalized occurence levels -> rarity-weighted richness
rwr <- sum(sp_rasters_normalized, na.rm = TRUE)
# Replace 0 with NA
rwr[rwr ==  0] <- NA
# Rank all cell according the RWR
rwr_rank <- rank(raster::getValues(rwr), ties.method = "random", na.last = "keep")
# Scale ranks to range [0, 1]
rwr_rank <- rwr_rank / max(rwr_rank, na.rm = TRUE)
e <- raster::extent(sp_rasters)
# Transform ranks back to a raster
rwr_rank <- raster(e, nrows = nrow(rwr), ncols = ncol(rwr), vals = rwr_rank)

levelplot(rwr_rank, margin = FALSE, main = "RWR", scales = list(draw = FALSE),
                 par.settings = rasterTheme(region = z_colors_RdYlBu$colors), 
                 at = z_colors_RdYlBu$values)

2.2 Zonation (CAZ + ABF)

The RWR priority rank from the prevoius section is compared against two cell-removal rules found in Zonation software: additive benefit function (ABF) and core-area Zonation (CAZ). Whereas ABF somewhat emphasizes the richness of features in any given cell, CAZ somewhat emphasizes rarity. Previously e.g. Albuquerque and Beier (2009) have compared RWR against CAZ variant of Zonation. If step-functions are used as the benefit function in ABF Zonation (which is typically not the case), it equals maximum coverage (Moilanen and Arponen 2011).

We’re not going to run the actual Zonation analysis, but assume that basic runs 01_core_area_zonation and 02_additive_benefit_function have been executed locally. Then, we just do load the results.

For additive benefit function, the \(z\) parameter has been set to 0.25 for all species.

caz_rank <- raster(file.path(z_tutorial_dir, "basic", "basic_output", 
                             "01_core_area_zonation", "01_core_area_zonation.rank.asc"))
levelplot(caz_rank, margin = FALSE, 
          par.settings = rasterTheme(region = z_colors_RdYlBu$colors), 
          scales = list(draw = FALSE), at = z_colors_RdYlBu$values, 
          main = "Core-area Zonation")

result_dir <- "/home/jlehtoma/VirtualBox VMs/trusty64/zonation/zonation-tutorial"
abf_rank <- raster(file.path(result_dir, "basic", "basic_output", 
                             "02_additive_benefit_function", 
                             "output_02_additive_benefit_function.rank.compressed.tif"))
levelplot(abf_rank, margin = FALSE,
          par.settings = rasterTheme(region = z_colors_RdYlBu$colors), 
          scales = list(draw = FALSE),
          at = z_colors_RdYlBu$values, main = "Additive benefit function")

2.3 Prioritization based on integer linear programming (ILP) and Gurobi

For a thorough description of using the Gurobi solver for solving conservation planning problem, see Beyer et al. (2016). Most of the approach and code below are based on the excellent tutorials Integer Programming with Gurobi for Reserve Design and Field Guide to ILP Solvers in R for Conservation Prioritization by Matt Strimas-Mackey, and on the prioritizr package by the same person.

Below, the prioritizaton is done by solving the conservation planning problem in a form a maximum coverage problem (MCP), or in other words, find the set of planning units that maximizes the overall level of representation (as measured by the normalized occurrence levels) across a suite of conservation features, while keeping cost within a fixed budget. It is worth mentioning that the Zonation algorithm is implicitly aiming at a maximum coverage type solution; by minimizing marginal loss one maximizes conservation value remaining at any specific level of cell removal (Moilanen 2007).

Here, the cost - and hence the budget - is defined to to be constant (1.0) in all cells. Thus, the budget corresponds to a particular number of planning units (cells) and we can maximize, for example, feature representation for the top 10% of the landscape. By repeating the maximum coverage optimization for a range of budgets (5%, 10%, 15%, …, 100% of the landscape), it is possible to create a nested hierarchy of solutions. This hierarchical solutions resembles the rankings resulting from the RWR and Zonation above.

DISCLAIMER: at this point, it is unclear to me whether the approach described above and implemented below (“hierarchical optimization”) is legit. For one thing, strictly speaking the solutions probably are not guaranteed to be nested, especially for more complex problems.

# NAs (NoData) must be raplaced wiht 0s for GUROBI
sp_rasters_filled <- sp_rasters
sp_rasters_filled[is.na(sp_rasters_filled)] <- 0
sp_rasters_filled <- raster::stack(sp_rasters_filled)

# Solve the maximum coverage problem for a range of target budgets
budgets <- seq(0.05, 1.00, 0.05)
results_mc <- list()
for (b in budgets) {
  b_cells <- b * raster::cellStats(cost, "sum")
  mc_model <- prioritizr::maxcover_model(x = cost, features = sp_rasters_filled, budget = b_cells)
  mc_results <- prioritize(mc_model, gap = 0.001)
  ##results <- protectr::gurobi_maxcoverage(cost, sp_rasters_filled, budget = b_cells)
  results_mc[[as.character(b)]] <- mc_results
}

# Function to translate Gurobi solutions to raster objects based on a template 
# raster
process_gurobi_results <- function(results, template_raster) {
  result_raster <- template_raster
  result_raster[] <- results$x
  # Fill in the nodata
  result_raster[is.na(template_raster)] <- NA
  # Make this a categorical raster
  result_raster <- ratify(result_raster)
  rat <- levels(result_raster)[[1]]
  # In some cases (most notably when budget is 100%) all cells might be 
  # selected. Check for this.
  if (nrow(rat) == 2) {
    rat$status <- c("Not Selected", "Selected")
  } else {
    rat$status <- c("Selected")
  }
  levels(result_raster) <- rat
  return(result_raster)
}
# Process the whole list of maximum coverage solutions
mcp_ilp <- raster::stack(lapply(results_mc, process_gurobi_results, template))
names(mcp_ilp) <- paste0("percent", budgets*100)

# Plot a selection of solutions (10%, 20%, 50%)
levelplot(subset(mcp_ilp, c(2, 4, 10)), main = "Maximum coverage solution\n(gap = 0.1%)",
          scales = list(draw = FALSE),
          col.regions = c("grey70", "#d7191c"),
          colorkey = list(space = "bottom", height = 1))

# Define a function for minmax-normalization
normalize <- function(x) {
  min <- raster::minValue(x)
  max <- raster::maxValue(x)
  return((x - min) / (max - min))
}
# Sum up all the layers in the stack -> result is a selection frequency
mcp_ilp_hier <- sum(mcp_ilp, na.rm = TRUE)
# Replace 0s with NAs
mcp_ilp_hier[mcp_ilp_hier == 0] <- NA
# Normalize value into scale [0, 1]
mcp_ilp_hier <- normalize(mcp_ilp_hier)
levelplot(mcp_ilp_hier, margin = FALSE, scales = list(draw = FALSE),
          par.settings = rasterTheme(region = z_colors_RdYlBu$colors), 
          at = z_colors_RdYlBu$values, main = "Gurobi MCP")

NOTE: the color scale above is not completely accurate, this needs to be fixed.

3. Comparison of the prioritizations

3.1 Spatial patterns

Let’s visually compare the top 10% solution patterns in RWR, Zonation and Gurobi solutions.

# Get the top 10% solutions for RWR and Zonation
rwr_top10 <- rwr_rank >= 0.9
# make this a categorical raster
rwr_top10 <- ratify(rwr_top10)
rat <- levels(rwr_top10)[[1]]
rat$status <- c("Not Selected", "Selected")
levels(rwr_top10) <- rat
abf_top10 <- abf_rank >= 0.9
# make this a categorical raster
abf_top10 <- ratify(abf_top10)
rat <- levels(abf_top10)[[1]]
rat$status <- c("Not Selected", "Selected")
levels(abf_top10) <- rat
p4 <- levelplot(rwr_top10, main = "RWR\ntop 10%",
          scales = list(draw = FALSE),
          col.regions = c("grey70", "#d7191c"),
          colorkey = list(space = "bottom", height = 1))
p5 <- levelplot(abf_top10, main = "Zonation ABF\ntop 10%",
          scales = list(draw = FALSE),
          col.regions = c("grey70", "#d7191c"),
          colorkey = list(space = "bottom", height = 1))
p6 <- levelplot(mcp_ilp[[1]], main = "Gurobi MCP\n10% budget",
          scales = list(draw = FALSE),
          col.regions = c("grey70", "#d7191c"),
          colorkey = list(space = "bottom", height = 1))
grid.arrange(p4, p5, p6, nrow = 1, ncol = 3)

Let’s also check that the solutions above include the same amount of cells included in the top 10%.

cell_count <- lapply(c(rwr_top10, abf_top10, mcp_ilp[[2]]), function(x) table(raster::getValues(x))) 
cell_count <- as.data.frame(do.call("rbind", cell_count))
names(cell_count) <- c("Low", "Top")
cell_count$Method <- c("RWR", "ABF", "MILP MCP")
cell_count <- dplyr::select(cell_count, Method, Low, Top)
knitr::kable(cell_count, align = c('r', 'c', 'c'), row.names = FALSE)

Method	Low	Top
RWR	99253	11029
ABF	99253	11029
MILP MCP	99254	11028

3.2 Quantitative performance

# Reclassify rasters
reclass_raster <- function(x, rcl) {
  
   assert_that(inherits(x, "RasterLayer"),
               inherits(rcl, "data.frame"),
               all(names(rcl) == c("from", "to", "value")))
  
  raster_rcl <- raster::reclassify(x, rcl, include.lowest = TRUE)
  raster_rcl <- ratify(raster_rcl)
  rat_rcl <- levels(raster_rcl)[[1]]
  rat_rcl$status <- rcl$value
  levels(raster_rcl) <- rat_rcl
  return(raster_rcl)
}
zonal_stats <- function(zonal_raster, values_stack, rcl, method) {
  rcl_raster <- reclass_raster(zonal_raster, rcl)
  rcl_zonal <- as.data.frame(raster::zonal(values_stack, rcl_raster, fun = 'sum'))
  # Manually add zone 0, in which nothing has been gained yet
  data_row <- rep(0, raster::nlayers(values_stack) + 1)
  names(data_row) <- c("zone", names(values_stack))
  rcl_zonal <- rbind(rcl_zonal, as.data.frame(t(data_row)))
  rcl_zonal <- rcl_zonal %>% 
    # Sort by zone, Zone 1 = 100-90% ... Zone 10 = 10-0% of the landscape
    arrange(zone) %>% 
    # Gather into long format keeping zone intact
    gather(species, value, -zone) %>%
    group_by(species) %>% 
    # Count cumulative value for distribution sum, that is zone1 + zone2 + ... + zone10
    # Then take the inverse of this, which is 1 -  dsum to show distribution remaining
    mutate(dsum = cumsum(value),
           drem = 1 - dsum) %>% 
    ungroup() %>% 
    # Finally, calculate mean distribution remaining per zone 
    group_by(zone) %>% 
    mutate(min_drem = min(drem),
           mean_drem = mean(drem),
           max_drem = max(drem),
           method = method)
  return(rcl_zonal)
}
# Reclassification table
rcl <- data.frame(from = seq(0, .9, 0.1), 
                  to = seq(0.1, 1.0, 0.1), 
                  value = 1:10)
rwr_zonal <- zonal_stats(rwr_rank, sp_rasters_normalized, rcl, "RWR")
abf_zonal <- zonal_stats(abf_rank, sp_rasters_normalized, rcl, "ABF")
caz_zonal <- zonal_stats(caz_rank, sp_rasters_normalized, rcl, "CAZ")
mcp_ilp_zonal <- zonal_stats(mcp_ilp_hier, sp_rasters_normalized, rcl, "ILP")
# Bind all zonal stats together
all_zonal <- dplyr::bind_rows(rwr_zonal, abf_zonal, caz_zonal, mcp_ilp_zonal) %>% 
  dplyr::select(zone, mean_drem, method)
ggplot(rwr_zonal, aes(x = zone, y = drem, color = species)) + geom_line() + ggtitle("RWR") +
  ylab("Distribution remaining\n") + xlab("\n Landscape selected") + 
  scale_x_continuous(breaks = 0:10, labels = paste0(seq(100, 0, -10), "%"))
ggplot(abf_zonal, aes(x = zone, y = drem, color = species)) + geom_line() + ggtitle("ABF") +
  ylab("Distribution remaining\n") + xlab("\n Landscape selected") + 
  scale_x_continuous(breaks = 0:10, labels = paste0(seq(100, 0, -10), "%"))
ggplot(caz_zonal, aes(x = zone, y = drem, color = species)) + geom_line() + ggtitle("CAZ") +
  ylab("Distribution remaining\n") + xlab("\n Landscape selected") + 
  scale_x_continuous(breaks = 0:10, labels = paste0(seq(100, 0, -10), "%"))
ggplot(mcp_ilp_zonal, aes(x = zone, y = drem, color = species)) + geom_line() + 
  ggtitle("Gurobi MCP") + ylab("Distribution remaining\n") + xlab("\n Landscape selected") + 
  scale_x_continuous(breaks = 0:10, labels = paste0(seq(100, 0, -10), "%"))
ggplot(all_zonal, aes(x = zone, y = mean_drem, color = method)) + geom_line() +
  ylab("Mean distribution remaining\n") + xlab("\n Landscape selected") + 
  ggtitle("Mean performance among the methods") +
  scale_x_continuous(breaks = 0:10, labels = paste0(seq(100, 0, -10), "%"))

What is striking here, is that there is very little difference between the methods. CAZ - as expected - has on average a little lower performance, others are very similar.

4. Literature cited

• Albuquerque, F. S., & Beier, P. (2015). Rarity-Weighted Richness: A Simple and Reliable Alternative to Integer Programming and Heuristic Algorithms for Minimum Set and Maximum Coverage Problems in Conservation Planning. Plos One, 10(3), e0119905. http://doi.org/10.1371/journal.pone.0119905
• Arponen, A., Cabeza, M., Eklund, J., Kujala, H., & Lehtomäki, J. (2010). Costs of integrating economics and conservation planning. Conservation Biology, 24(5), 1198–204. http://doi.org/10.1111/j.1523-1739.2010.01539.x
• Beyer, H. L., Dujardin, Y., Watts, M. E., & Possingham, H. P. (2016). Solving conservation planning problems with integer linear programming. Ecological Modelling, 328, 14–22. http://dx.doi.org/10.1016/j.ecolmodel.2016.02.005
• Camm, J. D., Polasky, S., Solow, A., & Csuti, B. (1996). A note on optimal algorithms for reserve site selection. Biological Conservation, 78(3), 353–355. http://doi.org/10.1016/0006-3207(95)00132-8
• Laitila, J., & Moilanen, A. (2012). Use of many low-level conservation targets reduces high-level conservation performance. Ecological Modelling, 247, 40–47. http://doi.org/10.1016/j.ecolmodel.2012.08.010
• Moilanen, A. (2007). Landscape Zonation, benefit functions and target-based planning: unifying reserve selection strategies. Biological Conservation, 134, 571–579. http://dx.doi.org/10.1016/j.biocon.2006.09.008
• Moilanen, A., & Arponen, A. (2011). Setting conservation targets under budgetary constraints. Biological Conservation, 144(1), 650–653. http://doi.org/10.1016/j.biocon.2010.09.006
• Underhill, L. G. 1994. Optimal and suboptimal reserve selection algorithms. Biological Conservation 70:85–87. http://dx.doi.org/10.1016/0006-3207(94)90302-6