Repository files navigation

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

Theprior3D package offers a comprehensive toolset for 3D systematicconservation planning, conducting nested prioritization analyses acrossmultiple depth levels and ensuring efficient resource allocationthroughout the water column (Doxa et al. 2022). Itprovides a structured workflow designed to address biodiversityconservation and management challenges in the 3 dimensions, such as theincorporation of multiple costs at different depth levels, whilefacilitating users’ choices and parameterization. The process initiatesfrom the deepest level and progressively moves toward the surface, byconducting a step-by-step prioritization analysis at each depth Figure1. The optimization result at each depth level is considered asa cost layer for the layer above. This approach gives priority to areas chosenin deeper waters when selecting areas at the subsequent upper level, thuscreating a cascading prioritization structure. Theprior3D package is builtupon theprioritizr package (Hanson et al. 2024), usingcommercial and open-source exact algorithm solvers that ensure optimal solutionsto prioritization problems.

Figure 1: Flow chart of the 3Dprioritization analysis for the four depth zones considered in the Doxaet al. (2022) analysis

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

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 workflowproceeds as follows:

1. Use thesplit_rast() function to convert 2D distribution rastersof biodiversity features into a 3D format.

2. Use theprioritize_3D() function to set the optimization problemand define its parameters. This function also solves the problem andprovides the solution in the form of a map.

3. Use theevaluate_3D() function to obtain detailed results in atabular format.

4. Use theplot_3D() function to generate graphs based on thesolution results.

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

Thespatial coherence of the solution maps can be evaluated througha post-processing analysis, which can be conducted after either thestep-by-step or the comparative analysis. The necessary functions forthis assessment are also provided within the package.

All the functions of the package prior3D can be installed in R via:

install.packages("prior3D")

Alternatively, development version of the package can be installed using:

if (!require(remotes)) install.packages("remotes")remotes::install_github("cadam00/prior3D")

Doxa A, Adam C, Nagkoulis N, Mazaris AD, Katsanevakis S, 2025. prior3D: An Rpackage for three-dimensional conservation prioritization. Ecological Modelling499: 110919.https://doi.org/10.1016/j.ecolmodel.2024.110919

Let us consider the following dataset as an illustrative example. Itrepresents a subset of the species analyzed in Doxa et al.(2022). For simplicity reasons, we have included only20 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 firstdata.frame contains information about the features. If thebiodiversity features concern species then this data.frame mustindicate at least the species name and species classification aspelagic or benthic (mandatory). Additional optional data may includespecies assignment to prioritization groups and, if available, thespecies’ bathymetric range (min and max depth at which the speciesoccurs). The seconddata.frame is a prioritization weights table,where users can assign specific weights to different prioritizationgroups. These groups can represent any meaningful categorization forthe prioritization process, like taxonomical, functional, orconservation status categories, such as those defined by the IUCN.

2. Biodiversity distribution data in 2D raster form. These rasterscontain the information on the spatial distribution of the featuresacross the study area. Biodiversity distribution information canrepresent either presence-absence data (binary) or any continuousinformation, such as biomass/abundance, probability of occurrences.

# Import prior3D R packagelibrary(prior3D)# Species information tabledata(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 formbiodiv_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 theMediterranean 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 tothe 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 plainsand abyssal zone) (Figure2).

Figure 2: The study area and theconsidered depth zones

To conduct the analysis, a SpatRaster object containing bathymetric datafor the planning site is needed. This raster should represent depthswith negative values and match the extent and resolution of thebiodiversity rasters. Alternatively, if bathymetry maps of greaterresolution and broader extent are available, they can also be used, asthe prior3D functions internally conduct cropping and resampling tomatch the biodiversity data. Producing the final depth raster thatdelineates the desired depth zones is also produced by the prior3Dfunctions.

# Biodiversity distribution data in 2D raster formdepth_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

Transforming Biodiversity Distributions into Multilevel (3D) Data

Thesplit_rast() function is used to convert 2D distributions ofbiodiversity features (rasters) into a 3D format.

# Splitting features' 2D distributions into 3D onessplit_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 eachbathymetric layer, necessary for the analysis next steps.

The 3D prioritization algorithm is implemented using theprioritize_3D() function, the core function of theprior3Dpackage. This function uses the list generated from thesplit_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",#"shuffle"portfolio_opts=list(number_solutions=10))

## Budget: 0.3

Notes:

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

budget_weights: Theprioritize_3D() function allows users to specifyhow the budget is distributed among depth levels. Three allocationmethods are available:

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

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

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

Prioritization Maps

Theprioritize_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 protect10%, 30% and 50% of the area) result in cumulative maps, illustratingareas selected by various budget levels. Although this output follows adifferent approach, it resembles to typical Zonation output maps.

Figure 3: Prioritization maps for single andmultiple budget percentages

# Create plot of outputs for a single budget percentageplot_3D(single_3D,to_plot="all",add_lines=FALSE)

Figure 4: Output for 30% budget percentage

And for multiple budgets

# Create plot of outputs for multiple budget percentagesmultiple_3D<- prioritize_3D(split_features=split_features,depth_raster=depth_raster,breaks= c(0,-40,-200,-2000,-Inf),biodiv_df=biodiv_df,budget_percents= seq(0,1,0.1),budget_weights="richness",threads=parallel::detectCores(),portfolio="gap",portfolio_opts=list(number_solutions=10))

## Budget: 0## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Budget: 0.1## Budget: 0.2## Budget: 0.3## Budget: 0.4## Budget: 0.5## Budget: 0.6## Budget: 0.7## Budget: 0.8## Budget: 0.9## Budget: 1## Warning: Problem failed presolve checks.#### These checks indicate that solutions might not identify meaningful priority## areas:#### ✖ Budget is greater than the total cost of selecting all planning units.## → Maybe you made a mistake when setting the `budget` in the objective function?#### ℹ For more information, see `presolve_check()`.

plot_3D(multiple_3D)

Figure 5: Output for multiple budgetpercentages

To facilitate comparisons between 3D and 2D approaches, thecompare_2D_3D() function is provided in the package. This functionenables users to conduct all the above mentioned steps of analysis (datageneration, setting and solving the optimization problem and producingoutputs), by executing both 2D and 3D approaches, with similar settings,that facilitate comparisons. The functionplot_Compare_2D_3D()generates corresponding maps and graphs for both approaches.

out_2D_3D<- Compare_2D_3D(biodiv_raster=biodiv_raster,depth_raster=depth_raster,breaks= c(0,-40,-200,-2000,-Inf),biodiv_df=biodiv_df,budget_percents= seq(0,1,0.1),budget_weights="richness",threads=parallel::detectCores(),portfolio="gap",#"shuffle"portfolio_opts=list(number_solutions=10))

## Budget: 0## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Warning: Portfolio could only find 1 out of 10 solutions.## Budget: 0.1## Budget: 0.2## Budget: 0.3## Budget: 0.4## Budget: 0.5## Budget: 0.6## Budget: 0.7## Budget: 0.8## Budget: 0.9## Budget: 1## Warning: Problem failed presolve checks.#### These checks indicate that solutions might not identify meaningful priority## areas:#### ✖ Budget is greater than the total cost of selecting all planning units.## → Maybe you made a mistake when setting the `budget` in the objective function?#### ℹ For more information, see `presolve_check()`.## Warning: Problem failed presolve checks.#### These checks indicate that solutions might not identify meaningful priority## areas:#### ✖ Budget is greater than the total cost of selecting all planning units.## → Maybe you made a mistake when setting the `budget` in the objective function?#### ℹ For more information, see `presolve_check()`.

plot_Compare_2D_3D(out_2D_3D,to_plot="all",add_lines=TRUE)

Figure 6: Comparison of 2D vs 3D approachfor multiple budget percentages

The spatial coherence of the prioritization output (optimizationsolution) maps is assessed using three metrics: average surfaceroughness (SA), surface kurtosis (SKU), and the RAO index. These can beused for comparison among solution 2D and 3D solutions.

High SA values signify that there is a high spatial heterogeneity,indicating lower spatial coherence. High SKU indicate high spatialcoherence. Both SA and SKU are calculated using R packagegeodiv(Smith et al. 2023), applyinggeodiv::focal_metricsfunctions “sa” and “sku” to optimization solution rasters.

High RAO values suggest increased spatial heterogeneity, thus lowspatial coherence. To compute the RAO metric, a moving window approachis employed on optimization solution maps, using functionrasterdiv::paRao from R packagerasterdiv (Rocchini, Thouverai,et al. 2021;Rocchini, Marcantonio, et al.2021). The dimensions of the window chosen is 3×3.The new raster, which is a result of the application of the algorithm,is used to get an average RAO value for the whole raster.

coherence(out_2D_3D,w=3)

## Progress metrics:  1 / 1## Progress metrics:  1 / 1##  sa2D  sa3D sa2Dw sa3Dw## 3.806 3.503 2.413 2.118

Figure 7: SA

coherence(out_2D_3D,w=3,metric="sku")

## Progress metrics:  1 / 1## Progress metrics:  1 / 1## sku2D  sku3D sku2Dw sku3Dw## 0.347  1.138 -0.536 -0.374

Figure 8: SKU

coherence(out_2D_3D,w=3,metric="rao")

## 2D RAO###### Processing alpha: 1 Moving Window: 3###### Processing alpha: 1 Moving Window: 3#### [============================>----------------]  64% in  0s## [============================>----------------]  65% in  0s## [=============================>---------------]  66% in  0s## [=============================>---------------]  67% in  0s## [==============================>--------------]  69% in  0s## [==============================>--------------]  70% in  0s## [===============================>-------------]  71% in  0s## [================================>------------]  72% in  0s## [================================>------------]  73% in  0s## [=================================>-----------]  75% in  0s## [=================================>-----------]  76% in  0s## [==================================>----------]  77% in  0s## [==================================>----------]  78% in  0s## [===================================>---------]  80% in  0s## [===================================>---------]  81% in  0s## [====================================>--------]  82% in  0s## [====================================>--------]  83% in  0s## [=====================================>-------]  84% in  0s## [=====================================>-------]  86% in  0s## [======================================>------]  87% in  0s## [=======================================>-----]  88% in  0s## [=======================================>-----]  89% in  0s## [========================================>----]  90% in  0s## [========================================>----]  92% in  0s## [=========================================>---]  93% in  0s## [=========================================>---]  94% in  0s## [==========================================>--]  95% in  0s## [==========================================>--]  96% in  0s## [===========================================>-]  98% in  0s## [===========================================>-]  99% in  0s## [=============================================] 100% in  0s#### 3D RAO###### Processing alpha: 1 Moving Window: 3###### Processing alpha: 1 Moving Window: 3#### [=============================>---------------]  67% in  0s## [==============================>--------------]  69% in  0s## [==============================>--------------]  70% in  0s## [===============================>-------------]  71% in  0s## [================================>------------]  72% in  0s## [================================>------------]  73% in  0s## [=================================>-----------]  75% in  0s## [=================================>-----------]  76% in  0s## [==================================>----------]  77% in  0s## [==================================>----------]  78% in  0s## [===================================>---------]  80% in  0s## [===================================>---------]  81% in  0s## [====================================>--------]  82% in  0s## [====================================>--------]  83% in  0s## [=====================================>-------]  84% in  0s## [=====================================>-------]  86% in  0s## [======================================>------]  87% in  0s## [=======================================>-----]  88% in  0s## [=======================================>-----]  89% in  0s## [========================================>----]  90% in  0s## [========================================>----]  92% in  0s## [=========================================>---]  93% in  0s## [=========================================>---]  94% in  0s## [==========================================>--]  95% in  0s## [==========================================>--]  96% in  0s## [===========================================>-]  98% in  0s## [===========================================>-]  99% in  0s## [=============================================] 100% in  0s#### $rao2D.mean## [1] 1.872#### $rao3D.mean## [1] 1.709#### Warning messages:## 1: In rasterdiv::paRao(x = raster_data_2D, window = w, na.tolerance = 1,  :##   Simplify=0. Rounding data to 0 decimal places.## 2: In rasterdiv::paRao(x = raster_data_2D, window = w, na.tolerance = 1,  :##   Input data are float numbers. Converting data to integer matrices.## 3: In rasterdiv::paRao(x = raster_data_3D, window = w, na.tolerance = 1,  :##   Simplify=0. Rounding data to 0 decimal places.## 4: In rasterdiv::paRao(x = raster_data_3D, window = w, na.tolerance = 1,  :##   Input data are float numbers. Converting data to integer matrices.

Figure 9: RAO

Doxa, Aggeliki, Vasiliki Almpanidou, Stelios Katsanevakis, Ana MQueirós, Kristin Kaschner, Cristina Garilao, Kathleen Kesner-Reyes, andAntonios D Mazaris. 2022. “4D marineconservation networks: Combining 3D prioritization of present and futurebiodiversity with climatic refugia.”Global Change Biology 28 (15): 4577–88.https://doi.org/10.1111/gcb.16268

Hanson, Jeffrey O, Richard Schuster, Nina Morrell, MatthewStrimas-Mackey, Brandon P M Edwards, Matthew E Watts, Peter Arcese,Joseph Bennett, and Hugh P Possingham. 2024.prioritizr: Systematic ConservationPrioritization in R.https://prioritizr.net

Rocchini, Duccio, Matteo Marcantonio, Daniele Da Re, Giovanni Bacaro,Enrico Feoli, Giles Foody, Reinhard Furrer, et al. 2021.“From zero to infinity: Minimumto maximum diversity of the planet by spatio-parametric Rao’s quadraticentropy.”Global Ecology and Biogeography 30 (5): 2315.https://doi.org/10.1111/geb.13270

Rocchini, Duccio, Elisa Thouverai, Matteo Marcantonio, MartinaIannacito, Daniele Da Re, Michele Torresani, Giovanni Bacaro, et al.2021. “rasterdiv - An Information Theory tailored R package for measuring ecosystemheterogeneity from space: To the origin and back.”Methods in Ecologyand Evolution 12 (6): 2195.https://doi.org/10.1111/2041-210X.13583

Smith, Annie C., Phoebe Zarnetske, Kyla Dahlin, Adam Wilson, and AndrewLatimer. 2023.Geodiv: Methodsfor Calculating Gradient Surface Metrics.https://doi.org/10.32614/CRAN.package.geodiv