Exploring Zonation 5 Settings and Outputs
Thiago Cavalcante
2023-12-12
1 Setting the stage
Spatial conservation prioritization is a crucial process in modern conservation science that involves the use of quantitative techniques for systematic assessment of conservation priorities. Considering a range of ecological and spatial criteria, limited resources and efforts can be efficiently allocated through the use of spatial prioritization softwares, such as Zonation, a widely employed tool in conservation planning. This maximizes the impact of conservation actions, ensuring a more effective approach to safeguarding biodiversity and ecological integrity. Zonation applies spatial analysis techniques to prioritize areas for conservation based on their ecological importance (i.e., considering the ecological features they contain). This is accomplished by systematically ranking sites according to their ability to optimize the representation of diverse species and/or habitats, facilitating the identification of critical conservation areas (Moilanen, Lehtinen, et al. 2022).
This illustrative workflow aims to serve as a general introduction for Zonation analysis, emphasizing the utilization of species distribution maps. We started by employing basic Zonation runs to establish fundamental conservation priorities, using species distribution models as our primary input. Subsequently, we explore more complex configurations, considering factors such as weighting schemes, dispersal and projected future distributions. While offering basic Zonation configurations, this tutorial is primarily focused on integrating considerations for connectivity and climate change into prioritization exercises. Most of this analytical approach was inspired by and built upon the ideas from Kujala et al. (2013). Additionally, at least three important references (among many others, as one might expect) are highly valuable for following up on and facilitating further exploration (Moilanen, Lehtinen, et al. 2022; Lehtomäki et al. 2016; Moilanen, Kohonen, et al. 2022). I hope it can provide a good starting point for those interested in applying Zonation in their own projects.
1.1 Species distribution models
Species distribution models (SDMs) used in this tutorial were developed within the framework of the NaturaConnect project (https://naturaconnect.eu/). SDMs were built for nine mammal species across Europe (Figure 1.1).
Figure 1.1: Modelling domains for SDMs highlighted in translucent blue. Map projection WGS84.
The species were Bison bonasus, Canis lupus, Capra ibex, Cervus elaphus, Crocidura sicula, Crocidura zimmermanni, Gulo gulo, Lepus europaeus, and Spalax antiquus (sp1-sp9, respectively: Figure 1.2).
Figure 1.2: Current habitat suitability modelled for the nine species across Europe.
Occurrence data were obtained from GBIF and filtered based on IUCN ranges for Canis lupus, Cervus elaphus, Lepus europaeus, Capra ibex, while for species with presence in fewer than 50 1km² resolution cells (Crocidura sicula, Gulo gulo, Bison bonasus, Crocidura zimmermanni, Spalax antiquus), only IUCN range was used. For modelling, occurrence data was subsampled aiming for an ideal 3:1 ratio of presences to absences. The chosen algorithm was XGBOOST. For cross-validation, a random split (80% of the dataset) with three replicates was used. The three cross-validation replicates were ensembled using the mean. Future climate projections (2041-2070) are under the SSP370 scenario and the GFDL climate model. SDMs were upscaled to 10km² resolution. Please note that these SDMs were built exclusively for testing purposes.
2 Replicability
All files, setups and R functions used in this tutorial are freely available on GitHub (https://github.com/thiago-cav/zonation_workflow). This framework uses a single root folder for Zonation and R softwares, allowing easy replication of all setups and R functions. The project folders structure (important for this workflow both from the Zonation and the R software perspective) was inspired by Lehtomäki et al. (2016) recommendations. To replicate all the analysis variants with your own data, simply replace the SDM files in the ‘data’ folder. If you use the same file naming convention and have the same number of species as in this workflow (see feature_list.txt files in the ´setup´ folders), you don’t need to do anything else besides running the Zonation for each setup. In Windows, the Zonation call can be launched using the batch files (.cmd) available in the ´setup´ folders (see Moilanen, Kohonen, et al. 2022 for more details). If you have more species or other features, you will need to edit all the files contained in each setup folder (see editing steps and instructions below).
For now, there are two R functions available (used to plot the figures below). The first function ‘plot_rankmaps’ can be used to access the priority rank maps from Zonation. The function has two arguments: variants and cols. You can choose the variants you want to plot and specify the structure of the plot in terms of number of columns.
# Usage example:
plot_rankmaps(variants = 1:6, cols = 3)# plot rankmaps from setup folders 1 to 6 in three columns (default is two)The other function ‘compare_perf_curves’ was created to compare mean and miminum performance curves between different Zonation variants. It follows the same structure as the previous one, but the curves are plotted in the same graph. In this case, the only argument needed is the Zonation variants.
# Usage example:
compare_perf_curves(variants = 3:6)# plot mean performance curves from setup folders 3 to 6The two functions can be found in the ‘r_functions’ folder within the repository, with reasonably informative documentation. Please remember to load the functions and install/load the external packages: cptcity, ggplot2, and terra. Also make sure you are in the root folder before using the functions.
3 Minimal setups
3.1 Setup 1
A basic analysis using the species distribution layers for the present time (hereafter “baseline”) and the ABF (additive benefit function) marginal loss rule. ABF emphasize high mean coverage across features, or areas with high feature (species) richness. The feature weights option was activated, but with each feature being assigned the same weight. The Zonation command in the batch file was:
z5 -w --mode=ABF --gui settings.z5 output
The priority rank map derived from a Zonation analysis, as shown in Figure 3.1, reveals the relative importance of different areas for nature conservation, taking into account the input features and analysis settings used in the Zonation process. Every grid cell in the map (10km² resolution) has a rank value, which ranges from 0 (lowest priority, depicted in blue) to 1 (highest priority, in red).
Figure 3.1: Priority rank map from Setup 1.
Since we are focusing on just nine species in this Zonation workflow, emphasizing single-feature areas could be an important goal. Therefore, our first adjustment involved the choice of a different marginal loss rule variant. For Setup 2, we opted for the Core Area Zonation CAZMAX, which prioritizes high-occurrence areas for all features, including exclusive single-feature zones, even at the cost of losing some average coverage (see Section 3.4 at Moilanen, Kohonen, et al. 2022).
3.2 Setup 2
A basic analysis using the baseline distribution layers and the CAZMAX (core area zonation) marginal loss rule. The feature weights option was activated, but with each feature being assigned the same weight.
The command was:
z5 -w --mode=CAZMAX --gui settings.z5 output
Figure 3.2: Priority rank map from Setup 1 and Setup 2.
In Figure 3.2, it is evident that changing the marginal loss rule leads to distinct variations in the distribution of priority areas, particularly noticeable in the western and northeastern regions of Europe. For example, we can observe a large high-priority area in northeastern Europe that is exclusive to a few species (e.g., sp7; Figure 1.2).
3.3 Setup 3
In Setup 3, we fine-tuned the analysis by assigning weights to individual features. These weights reflect perceived priority in the Zonation analysis. Features with higher assigned weights hold greater significance in the prioritization process. We designated weights for each taxon to reflect increasing levels of vulnerability based on IUCN categories. We followed an ascending order, increasing weights on a logarithmic scale: 1.0 for Least Concern, 2.0 for Near Threatened, and 4.0 for Endangered species (Table 3.1).
| weight | filename |
|---|---|
| 2 | ../data/sp1.tif |
| 1 | ../data/sp2.tif |
| 1 | ../data/sp3.tif |
| 1 | ../data/sp4.tif |
| 1 | ../data/sp5.tif |
| 4 | ../data/sp6.tif |
Please note that all tables in this tutorial were plotted for illustration purposes only. For precise configuration details and formatting styles, refer to the original files in the setup folders within the repository.
The command for this setup was:
z5 -w --mode=CAZMAX --gui settings.z5 output
Figure 3.3: Priority rank map from Setup 2 and Setup 3.
When assigning weights to individual features, we can observe conspicuous variations in the distribution of priority areas (Figure 3.3). This occurs because the marginal loss rule is now counterbalanced by species weights (see Section 3.4.3 at Moilanen, Kohonen, et al. 2022 for more details).
3.4 Setup 4
Using the same configuration from Setup 3, we can also identify priority areas for the future by just replacing baseline distributions with their projected future distributions (Table 3.2).
| weight | filename |
|---|---|
| 2 | ../data/sp1_futureM.tif |
| 1 | ../data/sp2_futureM.tif |
| 1 | ../data/sp3_futureM.tif |
| 1 | ../data/sp4_futureM.tif |
| 1 | ../data/sp5_futureM.tif |
| 4 | ../data/sp6_futureM.tif |
Figure 3.4: Baseline and future conservation priorities for nine mammal species (Canis lupus, Cervus elaphus, Crocidura sicula, Gulo gulo, Lepus europaeus, Capra ibex, Bison bonasus, Crocidura zimmermanni, and Spalax antiquus) across Europe.
In addition to the priority rank map, another crucial output derived from a Zonation analysis is the performance curve. Performance curves describe how much of the distribution of each feature is covered at different priority rank levels. This can be evaluated for each feature or across features. For example, we can compare mean and minimum performance curves for different Zonation 5 variants (Figure 3.5).
Figure 3.5: Mean and minimum performance curves obtained from Setups 1 and 3. The dashed black and red lines indicate the top 30% and 10% priority areas, respectevely (i.e., rank range = 0.7-1.0 and 0.9-1.0).
We can see that Setup 1 (ABF marginal loss rule without assigning weights) exhibits the highest mean coverage throughout the priority ranking in comparison to Setup 3 (CAZMAX marginal loss rule with different weights to individual features). This is expected because in Setup 1, we are not emphasizing any particular species using the weights, and ABF tends to promote high average coverage. However, the trade-off between average and minimum coverage is being affected by the feature weighting scheme (i.e., ABF is showing highest minimum performance). Moreover, the performance curves closely resemble a diagonal pattern, indicating a suboptimal performance likely associated with the wide and flat distribution of certain species across Europe (Figure 1.2).
4 Connectivity setups
4.1 Overview
Zonation 5 incorporates a range of connectivity methods that are based on dispersal kernels, which model how the dispersal of individuals declines with increasing distance from source locations. Different kernel shapes are available in Zonation 5, and they can be used to calculate connectivity for a single feature, for pairs of features, or for a group of interacting features (see Section 3.7 at Moilanen, Kohonen, et al. 2022 for a general description of the connectivity measures).
4.2 Dispersal kernels
There are four types of dispersal kernels currently available in Zonation 5: uniform, triangular, negative exponential, and Gaussian. In all kernels, we need to specify the \(\alpha\) parameter, which denotes a value greater than zero that defines the distance decay (i.e., the bandwidth) of the kernel. Additionally, the cut parameter in the 2D Gaussian and negative exponential kernels specifies the number of decay distances at which the kernel is truncated to zero, effectively limiting the influence of grid cells beyond this distance. The choice of the cut parameter is often influenced by the spatial extent and the scale at which you want to model the dispersal. If your data covers a large area, setting the cut parameter to a value that corresponds to a few times the dispersal distance decay could be reasonable. This helps in reducing computation and trimming the thin, long tails of dispersal kernels that hold high uncertainty.
In Zonation 5, we need to specify the \(\alpha\) parameter in grid cell units (instead of directly in e.g. meters). For instance, let’s assume a cell size of 1km in the feature distribution file and a feature-specific dispersal capability of 3km. This sets the \(\alpha\) value to 3. However, if we have a species with a total dispersal capability (\(D\)) of 100 km, but the feature raster cell size is in degrees (e.g., 0.5 degrees), we need to convert the feature raster cell size from degrees to kilometers. This can be achieved by using an approximate conversion factor at the equator, where 1 degree is approximately equal to 111.12 km. Thus, we calculate the cell size in kilometers (\(d\)) by multiplying the feature raster cell size in degrees with the conversion factor:
\[ d = 0.5 \times 111.12 \approx 55.56 \text{ km} \]
With this information, we can move on to determining the \(\alpha\) parameter:
\[ \alpha = \frac{D}{d} = \frac{100}{55.56} \approx 1.8 \]
Please note that the variation in distance per degree of latitude is greatest near the poles, leading to less precise values for very high latitudes.
4.3 Setup 5
Baseline Connectivity
In this setup, we implemented the weight groups option and the single-feature connectivity technique of Zonation 5. We need some extra files and supplementary information for this purpose. Single-feature connectivity allows to perform various connectivity-related operations on individual feature layers and is commonly used to model metapopulation-like connectivity (see Section 5.6.2 in Moilanen, Kohonen, et al. 2022). It is a way to understand how well different parts of the habitat are connected to each other by using the dispersal kernels described above. The rationale behind this setup is that highest prioritization values are given to areas within the baseline distributions that are highly suitable and geographically close to each other given species dispersal limitations.
For this exercise, we retrieved estimated dispersal from the COMBINE dataset (Soria et al. 2021) to incorporate detailed species-specific information. Dispersal in this study was defined as the distance covered by a species from its place of birth to its breeding location. We used this distance to parameterize the decay distance (\(\alpha\)) in the dispersal kernels, thus constraining movement to a maximum number of cells as defined by species dispersal abilities. We set the cut parameter to three times the dispersal distance decay.
In the feature list, we need to include the columns “sf_conn” and “wgrp” (Table 4.1). The “sf_conn” is used to specify the type of interaction connectivity transform associated with each feature, aligning with the information provided in a connectivity link file (see details below). A value of -1 in the “sf_conn” column indicate that single feature connectivity is not applied to that feature.
| sf_conn | wgrp | weight | filename |
|---|---|---|---|
| -1 | 1 | 2 | ../data/sp1.tif |
| -1 | 1 | 1 | ../data/sp2.tif |
| -1 | 1 | 1 | ../data/sp3.tif |
| -1 | 1 | 1 | ../data/sp4.tif |
| -1 | 1 | 1 | ../data/sp5.tif |
| -1 | 1 | 4 | ../data/sp6.tif |
| -1 | 1 | 1 | ../data/sp7.tif |
| -1 | 1 | 1 | ../data/sp8.tif |
| -1 | 1 | 4 | ../data/sp9.tif |
| 1 | 2 | 2 | ../data/sp1.tif |
| 2 | 2 | 1 | ../data/sp2.tif |
| 3 | 2 | 1 | ../data/sp3.tif |
The “wgrp” column in this feature list specifies the weight group each feature is connected to. This detail is outlined in a separated text file (weight_groups.txt) with two columns: (i) group ID numbers and (ii) weight assigned to each group. To calibrate the emphasis between local habitat quality and connectivity (which involves significant uncertainty in terms of species-specific dispersal), we assigned a weight of 1.5 to the untransformed baseline layers and a weight of 1.0 to the connectivity layers:
| 1 | 1.5 |
| 2 | 1.0 |
We also need a connectivity file (connectivity_list.txt), which is a plain text file with three columns (Table 4.2).
| 1 | 1 | negexp(4, 12) |
| 2 | 1 | negexp(4, 12) |
| 3 | 1 | negexp(2, 6) |
| 4 | 1 | negexp(2, 6) |
| 5 | 1 | negexp(1, 3) |
| 6 | 1 | negexp(1, 3) |
Column Explanations:
ID number that links with the ‘sf_conn’ column of the feature list file.
Whether the connectivity-transformed feature is multiplied by local habitat quality or not 1/0.
The dispersal kernel used.
For this exercise, we employed the negative exponential kernel. We set a value of 1 in the second column of the connectivity list, which means that the kernel-smoothed distribution was multiplied (pixel-wise) by the original input layers highlighting areas that have high local habitat quality and are well connected.
To sum up, we are calibrating the emphasis on habitat quality (from the SDMs) by (i) multiplying the SDMs with the connectivity layers (second column in the connectivity list), (ii) incorporating all models in the prioritization (SDMs + connectivity-transformed layers), and (iii) assigning a high weight to non-transformed SDMs. My assumption in this exercise was that uncertainties are greater for dispersal (i.e., connectivity parameters) compared to habitat quality.
The Zonation 5 call was:
z5 -dwW --mode=CAZMAX --gui settings.z5 output
Figure 4.1: Baseline conservation priorities (with connectivity) for nine mammal species (Canis lupus, Cervus elaphus, Crocidura sicula, Gulo gulo, Lepus europaeus, Capra ibex, Bison bonasus, Crocidura zimmermanni, and Spalax antiquus) across Europe.
While it seems that connectivity had a relatively modest impact on the prioritization solution, the distinctions are visually discernible in the maps (Figure 4.1). This can be attributed to the large cell size used for this workflow (10 x 10 km cells). With cells of this magnitude, it is reasonable to assume that they are large enough to encompass population dynamics for most species, making connectivity considerations less crucial. However, connectivity is typically considered desirable in spatial conservation planning, as it plays a significant role in maintaining biodiversity and ecosystem health.
Figure 4.2: Comparison of baseline conservation priorities for nine mammal species (Canis lupus, Cervus elaphus, Crocidura sicula, Gulo gulo, Lepus europaeus, Capra ibex, Bison bonasus, Crocidura zimmermanni, and Spalax antiquus) across Europe, illustrating scenarios with and without connectivity (left and right, respectively).
4.4 Setup 6
Climate resilient priority areas
Now let’s dive into the more complex Zonation configuration, where things get a bit wild and undeniably interesting. In this setup, we employed the functionalities of individual weights, weight groups, output groups and output transformed layers. We also used the Interaction Connectivity technique in Zonation 5, which allows calculation of connectivity between two distributions (see Section 5.6.3 in Moilanen, Kohonen, et al. 2022). Interaction connectivity covers a broad range of applications, such as biotic interactions, connectivity from disturbance sources (e.g., roads), and connectivity to protected areas. In this setup, we employed to model connectivity between the baseline distributions and projected future distributions.
We cover four theoretical scenarios. The first two, the baseline and future core areas, have already been introduced in Setups 3 and 4. The third scenario involves implementing baseline to future connectivity. The rationale behind this approach is that highest prioritization values are given to areas within baseline distributions that are highly suitable and geographically close to the expected future distribution given species dispersal limitations. These areas are expected to function as source areas from where dispersal to future distribution core areas might occur (Kujala et al., 2013). And the last one is the future to baseline connectivity. In this scenario, the highest values are given to highly suitable future areas that are well connected to baseline areas during the prioritization. These areas are expected to function as stepping stones assisting species in reaching their future core areas (Kujala et al., 2013). Considering the time step interval for future projections, we scaled up the dispersal estimates from COMBINE dataset by a factor of 50 to provide a rough estimate of the projected maximum dispersal distance for each individual species up until the mid-21st century (2041-2070). To wrap it up, for each species, we have four distribution layers to be covered in the priority setting: baseline, future, source areas (connectivity from the baseline to the future), and stepping stones (connectivity from the future to the baseline).
We employed the grouping weight feature in Zonation 5 to implement a complementary weighting strategy. The idea here is to acknowledge that the future is inherently more uncertain than the present (Kujala et al. 2013). The same is valid for habitat quality in comparison to connectivity.
The weighting scheme was defined as follows:
\[w(B_j) > w(C_{BF_j}) \geq w(F_j) > w(C_{FB_j})\]
, where j acts as the index for species, \(w(B_{j})\) represents the weight assigned to the baseline distribution, \(w(C_{BF_j})\) stands for the weight assigned to baseline to future connectivity, \(w(F_{j})\) denotes the weight assigned to the future distribution, and \(w(C_{FB_j})\) for the weight assigned to future to baseline connectivity. In this weighting scheme, I am assuming that the connectivity from the future to the present is marked by a “dual uncertainty” because it relies on uncertainties from the projected future distributions and species-specific dispersal estimates. We followed a descending logarithmic scale (4.0, 2.0, 1.0).
The first difference from Setup 5 is that we use a distinct column name in the feature list (Table 4.3). Instead of “sf_conn”, we use “ia_conn1”. The function of this column remains unchanged from the previous configuration – it specifies the interaction connectivity transform associated with each feature, aligning with the information provided in the connectivity link file.
| ia_conn1 | wgrp | weight | group | tr_out | filename |
|---|---|---|---|---|---|
| -1 | 1 | 2 | 1 | 0 | ../data/sp1.tif |
| -1 | 1 | 1 | 1 | 0 | ../data/sp2.tif |
| -1 | 1 | 1 | 1 | 0 | ../data/sp3.tif |
| -1 | 1 | 1 | 1 | 0 | ../data/sp4.tif |
| -1 | 1 | 1 | 1 | 0 | ../data/sp5.tif |
| -1 | 1 | 4 | 1 | 0 | ../data/sp6.tif |
You can notice that we have two new columns in the feature list: “group” and “tr_out”. We used the output groups option (column “group”) to generate separated summary curves for the four different scenarios of this setup. These group summary curves can be used, for example, to evaluate the trade-offs between baseline and future conservation weighting schemes (see Kujala et al., 2013 for some ideas). The output transformed layers option (column “tr_out”) was also employed to access the pre-processed connectivity layers as individual files. These layers can be used for error checking (verifying transform performance) and in various other applications beyond Zonation.
The connectivity link file for this setup has six columns instead of three (Table 4.4).
| 1 | 1 | 1 | 1 | negexp(200, 600) | ../data/sp1_futureM.tif |
| 2 | 1 | 1 | 1 | negexp(200, 600) | ../data/sp2_futureM.tif |
| 3 | 1 | 1 | 1 | negexp(100, 300) | ../data/sp3_futureM.tif |
| 4 | 1 | 1 | 1 | negexp(100, 300) | ../data/sp4_futureM.tif |
| 5 | 1 | 1 | 1 | negexp(50, 150) | ../data/sp5_futureM.tif |
| 6 | 1 | 1 | 1 | negexp(50, 150) | ../data/sp6_futureM.tif |
Column Explanations:
ID number that links with the ‘ia_conn1’ column of the feature list file.
Whether the interaction is positive or negative, indicated by 1 or 0, respectively.
The scaling factor β for the strength of the negative interaction; use 1.0 by default.
Whether to renormalize the negative interaction or not, indicated by 1/0.
The dispersal kernel used for the interaction.
The transforming file.
In this exercise, we assigned a value of 1 to columns 2 through 4 to denote positive interactions. However, this type of Zonation connectivity covers a broader range of applications as previously disclosed, including accounting for negative interactions. Please be also aware that for the Interaction Connectivity technique, we require additional raster files, as specified in column six. Since we are modeling connectivity between baseline and future distributions, the interaction connectivity transform will be unique for each feature. However, in many cases, the same type of connectivity interaction will impact multiple species or habitats (e.g., protected areas network or negative connectivity to roads). In these cases, the same link layer is applied to multiple input features.
The Zonation 5 call was:
z5 -iwWgf --mode=CAZMAX --gui settings.z5 output
Figure 4.3: Climate resilient conservation priorities for nine mammal species (Canis lupus, Cervus elaphus, Crocidura sicula, Gulo gulo, Lepus europaeus, Capra ibex, Bison bonasus, Crocidura zimmermanni, and Spalax antiquus) across Europe.
6 Acknowledgments
A mega thanks to Heini Kujala and Atte Moilanen for their insightful comments and feedback, which greatly enriched this tutorial. And a big special appreciation to Wilfried Thuiller and Rémi Patin from LECA-CNRS (Grenoble-FRA), who kindly shared the test SDMs. Many thanks also to the CBIG team for their valuable discussions exploring the intricate details of Zonation’s configuration.