Prerequisites

To follow along, you will need several R packages including terra for loading and manipulating spatial data, and prioritizr for the prioritization process. We will also need a solver that prioritizr uses to solve the prioritization problems. For simplicity, we can install the open source solver lpsymphony. Finally we will use the oceandatr to get the data for the prioritization. Use the following code to do all the install work.

install.packages(c("prioritizr", "terra", "slam"))

if (!require(remotes)) install.packages("remotes") #you might already have this: for installing from Github and Bioconductor

remotes::install_bioc("lpsymphony")
options(timeout = 6000) #oceandatr is ~500MB, so need to increase the timeout
remotes::install_github("emlab-ucsb/oceandatr", force = TRUE) #force=TRUE ensures we get the most recent version on Github
remotes::install_github("bio-oracle/biooracler") #for biophysical data downloaded from Bio-Oracle by oceandatr

New we can load the libraries we will need.

library(oceandatr)
library(prioritizr)

Spatial prioritization primer

How do we decide which areas to protect? Or which areas should be for aquaculture, renewable energy development, or any other use? We need to think about our spatial priorities!

We normally have a set of spatial objectives we want to meet, e.g. protect 50% of all fish spawning grounds, protect 30% of all biodiversity, ensure that 50% of an area remains open to fishing. Lets look at the cartoon example below where the area of the ocean is divided into 9 grid squares (our planning grid), and there are 3 species present. Our objective is to protect 1 representative of each species; simple enough. But we also want to do this using the smallest possible area (minimum grid cells). And how about connectivity; we might want to protect grid cells next to each other. This is simple enough to work through for 9 grid cells and 3 species, but imagine thousands of grid cells and hundreds of species! This is where computer algorithms come in.

We will do a simple spatial prioritization using the prioritizr R package. This is widely used for conservation planning and emlab has used it for MPA planning in Bermuda, Montserrat, FSM and Maldives.

Spatial prioritization recipe

We need 3 things to conduct a basic spatial prioritization using prioritizr:

Planning grid: a grid of cells covering the area where we are doing planning, e.g. an EEZ
Features: spatial data for things we want to protect, e.g. habitats
Targets: how much of each feature do we want to protect, e.g. 30% of each habitat
Cost: spatial data for what we want to minimize impact to. In the ocean, we often use some measure of fishing (e.g. fishing effort or catch value), because we are often interested in minimizing the impact of protected areas to fishing communities.

We can get all the spatial data we need using oceandatr.

Bermuda example

As an example, we will use Bermuda, but the same process can be used for any EEZ or area of the ocean.

Planning grid

First we need to get the Bermuda EEZ and create a planning grid for it. We are going to create a raster planning grid, but oceandatr can also create sf grids. When creating the planning grid using the oceandatr function get_grid, we need to enter a coordinate reference system (crs). An equal area crs is normally best for spatial planning and suitable options can be found using projection wizard.

#get the Bermuda EEZ from marineregions.org
bermuda_eez <- get_boundary("Bermuda")

#now create a planning grid with 10 km x 10km grid cells, using an equal area crs
planning_grid <- get_grid(boundary = bermuda_eez, crs = "+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs", resolution = 10000)

#lets have a look
terra::plot(planning_grid)
terra::lines(bermuda_eez |> sf::st_transform(sf::st_crs(planning_grid))) #for plotting the outline of the EEZ

Features

Now we can get some features. Remember these are just spatial data for things that we want to protect. There are several sources of data available via oceandatr and there is a one stop function that can pull down a standard set of features, representing oceanic habitat types. These are what what emlab have been using for offshore prioritizations.

features <- get_features(spatial_grid = planning_grid) |>
  remove_empty_layers()#remove all rasters that are only zero (i.e. data is not present)
## Getting depth zones...
## Getting seamount data...
## Spherical geometry (s2) switched off
## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar
## Spherical geometry (s2) switched on
## Getting knoll data...
## Spherical geometry (s2) switched off
## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar
## Spherical geometry (s2) switched on
## Getting geomorphology data...
## Getting coral data...
## Getting environmental zones data... This could take several minutes

terra::plot(features)

Cost

The cost layer is what we want to minimize impact to. We often use some measure of fishing. Here we will use distance from shore which is a simple proxy of fishing effort for coastal fisheries. We get the inverse of the distance to shore because we want areas nearshore to be the most “costly” since that is where most fishing occurs in coastal fisheries.

#function is from the offshoredatr package
cost <- get_dist(spatial_grid = planning_grid, dist_to = "shore", inverse = TRUE)

terra::plot(cost)

Targets

We need to tell the prioritization how much of each feature to include in the solution, which is potential areas for protection in our case. We will set a target of 30% for all our features to ensure adequate representation of all habitats.

Prioritization

Now we can put together our prioritization “problem”:

prob <- problem(cost, features = features) %>% #setup the problem with cost and features
  add_min_set_objective() %>% #minimum set objective: achieve our objectives at minimum cost
  add_relative_targets(0.3) %>% #30% of each feature must be included in the solution
  add_binary_decisions() %>% #each grid cell is either in or out of the solution
  add_lpsymphony_solver() #the solver to use: can also set add_default_solver() to use the best available on your computer

#solve the problem!
sol <- solve(prob)

Solution

Areas which meet our targets, 30% protection of each feature, at the lowest possible cost, are selected.

terra::plot(sol)
terra::lines(terra::as.polygons(features$seamounts), col = "tan") #add outline of seamounts
terra::lines(bermuda_eez |> sf::st_transform(sf::st_crs(planning_grid))) #and the EEZ

We can check what percentage of the total area of each feature is represented in the solution: the column on the far right has this information.

eval_feature_representation_summary(prob, sol)
## # A tibble: 18 × 5
##    summary feature            total_amount absolute_held relative_held
##    <chr>   <chr>                     <dbl>         <dbl>         <dbl>
##  1 overall continental shelf             6             2         0.333
##  2 overall upper bathyal                 4             2         0.5  
##  3 overall lower bathyal                62            20         0.323
##  4 overall abyssal                    4582          1378         0.301
##  5 overall seamounts                   330            99         0.3  
##  6 overall knolls                     1034           311         0.301
##  7 overall Abyssal_Hills              2058           618         0.300
##  8 overall Abyssal_Plains             2280           707         0.310
##  9 overall Escarpments                 127            39         0.307
## 10 overall Major_ocean_basins         4642          1393         0.300
## 11 overall Ridges                       58            21         0.362
## 12 overall Terraces                      5             2         0.4  
## 13 overall antipatharia                 77            27         0.351
## 14 overall cold_corals                  19             6         0.316
## 15 overall octocorals                   82            25         0.305
## 16 overall enviro_zone_1              1006           307         0.305
## 17 overall enviro_zone_2              1376           413         0.300
## 18 overall enviro_zone_3              2272           682         0.300

Can also see how much of the whole planning area (Bermuda’s EEZ) is included in the solution

#eval_n_summary gives you how many planning units are in the solution. We can then divide this by the total number of planning units
as.numeric(eval_n_summary(prob, sol)[2])/as.numeric(terra::global(planning_grid, "sum", na.rm = TRUE)[1])
## [1] 0.3012462

Rapid intro to spatial prioritization

Jason Flower

2026-03-12