This research was conducted at the Department of Marine Sciences, University of the Aegean, Greece, supported by the European Union’s Horizon 2020 research and innovation programme HORIZON-CL6–2021-BIODIV-01–12, under grant agreement No 101059407, “MarinePlan – Improved transdisciplinary science for effective ecosystem-based maritime spatial planning and conservation in European Seas”.

1 Introduction to the prior3D Package (tutorial)

The prior3D package offers a comprehensive toolset for 3D systematic conservation planning, conducting nested prioritization analyses across multiple depth levels and ensuring efficient resource allocation throughout the water column (Doxa et al. 2022). It provides a structured workflow designed to address biodiversity conservation and management challenges in the 3 dimensions, such as the incorporation of multiple costs at different depth levels, while facilitating users’ choices and parameterization. The process initiates from the deepest level and progressively moves toward the surface, by conducting a step-by-step prioritization analysis at each depth Figure 1. The optimization result at each depth level is considered as a cost layer for the layer above. This approach gives priority to areas chosen in deeper waters when selecting areas at the subsequent upper level, thus creating a cascading prioritization structure. The prior3D package is built upon the prioritizr package (Hanson et al. 2024), using commercial and open-source exact algorithm solvers that ensure optimal solutions to prioritization problems.

Figure 1: Flow chart of the 3D prioritization analysis for the four depth zones considered in the @doxa20224d analysis

Figure 1: Flow chart of the 3D prioritization analysis for the four depth zones considered in the Doxa et al. (2022) analysis

This tutorial will guide you through the key functions of the package, from data preparation to generating informative outputs to address conservation challenges in diverse marine (and terrestrial) ecosystems and enable informed decision-making in biodiversity conservation, restoration and management.

2 Workflow: Running the analysis

The package provides two options for conducting analyses:

  1. Running a Step-by-Step 3D SCP analysis

  2. Running a Comparative Analysis of a 2D and a 3D SCP approach

When opting for the step-by-step analysis (first option), the workflow proceeds as follows:

  1. Use the split_rast() function to convert 2D distribution rasters of biodiversity features into a 3D format.

  2. Use the prioritize_3D() function to set the optimization problem and define its parameters. This function also solves the problem and provides the solution in the form of a map.

  3. Use the evaluate_3D() function to obtain detailed results in a tabular format.

  4. Use the plot_3D() function to generate graphs based on the solution results.

When opting for a comparative analysis of a 2D and a 3D SCP approach (second option), users can use the Compare_2D_3D() function. This function incorporates the aforementioned detailed workflow and applies it to both 2D and 3D approaches, streamlining and simplifying the analysis process for users. By using this function, users provide the input data, define the optimization problem and its parameterization, run the analysis and finally obtain the results in the form of maps, graphs and tables.

The spatial coherence of the solution maps can be evaluated through a post-processing analysis, which can be conducted after either the step-by-step or the comparative analysis. The necessary functions for this assessment are also provided within the package.

3 Illustration example

Let us consider the following dataset as an illustrative example. It represents a subset of the species analyzed in Doxa et al. (2022). For simplicity reasons, we have included only 20 species for demonstration purposes.

Biodiversity features

Two types of input data are needed for the biodiversity features.

  1. Species information tables in tabular form (data.frame). The first data.frame contains information about the features. If the biodiversity features concern species then this data.frame must indicate at least the species name and species classification as pelagic or benthic (mandatory). Additional optional data may include species assignment to prioritization groups and, if available, the species’ bathymetric range (min and max depth at which the species occurs). The second data.frame is a prioritization weights table, where users can assign specific weights to different prioritization groups. These groups can represent any meaningful categorization for the prioritization process, like taxonomical, functional, or conservation status categories, such as those defined by the IUCN.

  2. Biodiversity distribution data in 2D raster form. These rasters contain the information on the spatial distribution of the features across the study area. Biodiversity distribution information can represent either presence-absence data (binary) or any continuous information, such as biomass/abundance, probability of occurrences.

# Import prior3D R package
library(prior3D)

# Species information table
data(biodiv_df)
head(biodiv_df)
##               species_name pelagic min_z max_z
## 1   acanthocybium_solandri       1   -20     0
## 2    acantholabrus_palloni       0  -500   -30
## 3 acanthomysis_longicornis       0  -100    -2
## 4      abraliopsis_morisii       0 -3660     0
## 5          abralia_veranyi       0  -900    -1
## 6     abraliopsis_pfefferi       0  -750    -1
# Biodiversity distribution data in 2D raster form
biodiv_raster <- get_biodiv_raster()
biodiv_raster
## class       : SpatRaster 
## dimensions  : 31, 83, 20  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -5.5, 36, 30.5, 46  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : biodiv_raster.tif 
## names       : aapto~aptos, abiet~etina, abra_alba, abral~ranyi, abral~risii, abral~fferi, ... 
## min values  :        0.01,        0.01,      0.01,        0.01,        0.01,        0.01, ... 
## max values  :        1.00,        0.63,      1.00,        1.00,        1.00,        1.00, ...

Planning site, planning units and depth levels

In the illustration example, we consider as our planning site the Mediterranean Sea, with 0.5°x0.5° cells as our planning unites (PUs). We consider four depth levels: (i) 0 to 40 m (infralittoral zone, extending to the lower limit of photophilic algae and seagrasses), (ii) 40 to 200 m (circalittoral zone, continental shelf, animal-dominated), (iii) 200 to 2000 m (~continental slope), and (iv) exceeding 2000 m in depth (lower bathyal plains and abyssal zone) (Figure 2).

Figure 2: The study area and the considered depth zones

Figure 2: The study area and the considered depth zones

To conduct the analysis, a SpatRaster object containing bathymetric data for the planning site is needed. This raster should represent depths with negative values and match the extent and resolution of the biodiversity rasters. Alternatively, if bathymetry maps of greater resolution and broader extent are available, they can also be used, as the prior3D functions internally conduct cropping and resampling to match the biodiversity data. Producing the final depth raster that delineates the desired depth zones is also produced by the prior3D functions.

# Biodiversity distribution data in 2D raster form
depth_raster <- get_depth_raster()
depth_raster
## class       : SpatRaster 
## dimensions  : 31, 83, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -5.5, 36, 30.5, 46  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : depth_raster.tif 
## name        : depth_raster 
## min value   :  -4082.70312 
## max value   :     -6.60191

4 Step-by-Step 3D SCP analysis

4.1 Step 1: Data Preparation

Transforming Biodiversity Distributions into Multilevel (3D) Data

The split_rast() function is used to convert 2D distributions of biodiversity features (rasters) into a 3D format.

# Splitting features' 2D distributions into 3D ones
split_features <- split_rast(biodiv_raster,
                             depth_raster,
                             breaks = c(0, -40, -200, -2000, -Inf),
                             biodiv_df,
                             val_depth_range=TRUE)

The output is a list containing species distributions for each bathymetric layer, necessary for the analysis next steps.

4.2 Step 2: 3D Prioritization Algorithm

The 3D prioritization algorithm is implemented using the prioritize_3D() function, the core function of the prior3D package. This function uses the list generated from the split_rast() function and other necessary inputs.

single_3D <- prioritize_3D(split_features  = split_features,
                           depth_raster    = depth_raster,
                           breaks          = c(0, -40, -200, -2000, -Inf),
                           biodiv_df       = biodiv_df,
                           budget_percents = 0.3,
                           budget_weights  = "richness",
                           threads         = parallel::detectCores(),
                           portfolio       = "gap",
                           portfolio_opts  = list(number_solutions = 10))
## Budget: 0.3

Notes:

budget_percent: Contrarily to its strict economic definition, budget reflects the desired level of protection to be modeled. It ranges from 0 to 1, with 0 indicating no resources available for protection, while 1 signifies resources sufficient to protect the entire study area. Typically, setting a budget of 0.3 corresponds to the 30% conservation target (i.e. 30% of the total area set aside for conservation). Users also have the flexibility to define multiple budget levels using a vector, allowing for the exploration of various protection scenarios. For instance, a vector like c(0.1, 0.3, 0.5) represents three scenarios where 10%, 30%, and 50% of the study area are designated for protection.

budget_weights: The prioritize_3D() function allows users to specify how the budget is distributed among depth levels. Three allocation methods are available:

  1. Equal Distribution: Allocates an equal share of the budget to each depth level (budget_weights ="equal").

  2. Proportional to Area: Allocates budget based on the spatial extent of each depth level (budget_weights ="area").

  3. Proportional to Species Richness: Prioritizes budget allocation to depth levels with higher species diversity (number of species). (budget_weights = "richness")

4.3 Step 3: Generating Outputs

Prioritization Maps

The prioritize_3D() function is used to generate prioritization maps. Single budget settings (ex. total_budget=0.3) produce standard maps, as typical Marxan outputs. Multiple budgets, by using a vector (ex. c(0.1,0.3,0.5), indicating available resources sufficient to protect 10%, 30% and 50% of the area) result in cumulative maps, illustrating areas selected by various budget levels. Although this output follows a different approach, it resembles to typical Zonation output maps.