1. Introduction

  • This code aimed at automatizing the creation of input files to run the sensitivity analyses of Zonation, similar to those used in El-Gabbas et al. paper:

    El-Gabbas, Ahmed; Gilbert, Francis; and Dormann, Carsten F. (2020) Spatial conservation prioritisation in data-poor countries: a quantitative sensitivity analysis using multiple taxa. BMC Ecology (under review).

  • Please note that this appendix is not intended to provide an introduction to Zonation software. An introduction to the Zonation software can be found in Di Minin et al., 2014. I expect the reader to be familiar with Zonation and have read the paper in advance.

  • This html file was written in RMarkdown. The RMarkdown and sample data are available here.


1.1. Abbreviations used


Species weight

  • WTLoc: weighted by Egyptian national red-list status
  • WTWO: without red-list weighting
  • Uncert: weighted by predictive uncertainty
  • NoUncert: without predictive uncertainty weight


Connectivity

  • BQP: boundary quality penalty
  • BQP-WO: no connectivity analyses
  • BQP-Low /BQP-Med /BQP-Strng: low/medium/strong BQP curves


Modelling algorithms


Sampling bias


Surrogate groups

  • MMls: Mammals
  • Rep: Reptiles
  • Butr: Butterflies
  • AllSP: All three groups together


1.2. Input files

# loading required packages
require(raster)
require(tidyverse)
require(readr)
require(DT)


Species weights

  • Four *.csv files for different weighting options are available in the folder ZigInputFiles/Data/. Below is an example on one of these files: ZigInputFiles/Data/WeightsData_EN_Bias0.csv:

Setting files

  • Four *.dat setting files are available in the folder: ZigInputFiles/Dat/*.dat. Below is an example on one of these files (ZigInputFiles/Dat/ABF_Mask.dat):
[Settings]
removal rule = 2
warp factor = 10
edge removal = 0
add edge points = 50
use SSI = 0
SSI file name = SSI_list.txt
use planning unit layer = 0
planning unit layer file = PLU_file.asc
initial removal percent = 0.0
use cost = 1
cost file = Maps\CostLayer.tif
use mask = 1
mask file = Maps\Mask_PAs.tif
use boundary quality penalty = 1
BQP profiles file = Maps\BQPcurves.txt
BQP mode = 1
BLP = 0
use tree connectivity = 0
tree connectivity file = tree.txt
use interactions = 0
interaction file = interact.spp
annotate name = 0
logit space = 0
treat zero-areas as missing data = 0
z = 0.25
resample species = 0
[Info-gap settings]
Info-gap proportional = 0
use info-gap weights  = 0
Info-gap weights file = UCweights.spp
[Outputs]
output proportional loss ranking = 1


Predicted distribution maps

  • In this reproducible code, I used 5 species for each surrogate group (butterflies, mammals, reptiles).
  • The folder ZigInputFiles/Maps/ contains predicted distribution maps *.tif for each combination of species distribution model algorithm and sampling bias correction.

Example map:

Cost layer

  • The file ZigInputFiles/Maps/CostLayer.tif identifies urban areas which were given low priority values.

Current Protected Areas

  • The file ZigInputFiles/Maps/Mask_PAs.tif discriminates the Egyptian Protected Areas from other areas.

Response curve

  • The file ZigInputFiles/Maps/BQPcurves.txt contains data for different BQP curves used.

Response curves used in the boundary-quality penalty connectivity analyses. These curves describe a range of species sensitivity to habitat loss in neighbour cells, ranging from no response (1) to strong response (4). The x-axis shows the percentage of neighbour cells (specified by three values of radii) remaining: 100 represents no habitat loss; while 0 represents total habitat loss of all neighbour cells. The y-axis shows the percentage reduction in local cell value in response to habitat loss.


2. Preparing feature list files (.spp)

“The feature list file (.spp) includes a list of all biodiversity features (species, ecosystems, etc.) included in the analysis. This file is a text file that can be created in Notepad. [The file name extension .spp is not compulsory but we often use it to indicate that the file is a Zonation "species list file".] Each row in the file corresponds to a biodiversity feature. For each biodiversity feature, there are six columns of information which require values to be entered. Frequently, dummy values are used in columns that are not relevant for the particular case.” Di Minin et al., 2014

An example .spp file is shown below. For more information, see Di Minin et al., 2014

1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp1.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp2.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp3.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp4.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp5.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp6.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp7.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp8.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp9.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp10.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp11.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp12.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp13.tif
1   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp14.tif
3   1   2   1   1   ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp15.tif

The following code creates 640 .spp files:

  • 4 weights × 10 BQP combinations × 2 modelling algorithms × 2 sampling bias scenarios × 4 surrogate options


2.1. Preparing file names

The four weighting options are:

  • WTWO-NoUncert: no weighting
  • WTLoc-NoUncert: weighted by the Egyptian national red-list status only
  • WTWO-Uncert: weighted by predictive uncertainty only
  • WTLoc-Uncert: weighting by both
Weights <- c("WTWO-NoUncert", "WTLoc-NoUncert", "WTWO-Uncert", "WTLoc-Uncert")

Four BQP curves and 3 radii (10 options in total)

BQP_Curve <- expand.grid(c("BQP-WO", "BQP-Low", "BQP-Med", "BQP-Strng"), 1:3)
BQP_Curve <- paste0(BQP_Curve[,1], "-", BQP_Curve[,2])
BQP_Curve <- BQP_Curve[!BQP_Curve %in% c("BQP-WO-2", "BQP-WO-3")]

Two modelling algorithms × 2 sampling bias scenarios × 4 surrogate options

SpeciesMaps <- expand.grid(
  c("Mxnt", "EN"),
  c("BiasWO", "Bias0"),
  c("MMls", "Rep", "Butr", "AllSP")) 

SpeciesMaps <- paste0(
  SpeciesMaps[,1], "-",
  SpeciesMaps[,2], "-",
  SpeciesMaps[,3])

# preparing spp file names

AllSppFiles <- expand.grid(
  Weights = Weights,
  BQP_Curve = BQP_Curve,
  SpeciesMaps = SpeciesMaps)

AllSppFiles <- paste0(
  AllSppFiles[,1], "__",
  AllSppFiles[,2], "__",
  AllSppFiles[,3], ".spp")

The names of the 640 spp files:

Each .spp file has different number of rows based on the number of species: 5 species for analyses of each of butterflies, reptiles, and mammals; 15 species for all species analysis.

AllSppData <- vector(mode = "list", length = length(AllSppFiles))
names(AllSppData) <- AllSppFiles

Index_AllSP <- grep(pattern = "AllSP", x = names(AllSppData))
Index_MMls <- grep(pattern = "MMls", x = names(AllSppData))
Index_Rep <- grep(pattern = "Rep", x = names(AllSppData))
Index_Butr <- grep(pattern = "Butr", x = names(AllSppData))

for(i in Index_AllSP){
  AllSppData[[i]] <- data.frame(matrix(ncol = 6, nrow = 15))
}

for(i in c(Index_MMls, Index_Rep, Index_Butr)){
  AllSppData[[i]] <- data.frame(matrix(ncol = 6, nrow = 5))
}

# renaming columns
AllSppData <- lapply(AllSppData, FUN = function(x){
  names(x) <- c("Weight", "Alpha", "BQP_CurveNum",
                "BQP_Buff", "ABF_Rate", "SpeciesMaps")
  x$Alpha <- 1
  x$ABF_Rate <- 1
  x
})

Next, we change the options within each of the 640 .spp files

  • BQP radius
  • BQP response curves
  • Species
  • Weights


2.2. BQP radius options

Index_BQP_Buff1 <- grep(
  pattern = "__BQP-.*-1__", x = names(AllSppData))
for(i in Index_BQP_Buff1){
  AllSppData[[i]]$BQP_Buff <- 1
}

Index_BQP_Buff2 <- grep(
  pattern = "__BQP-.*-2__", x = names(AllSppData))
for(i in Index_BQP_Buff2){
  AllSppData[[i]]$BQP_Buff <- 2
}

Index_BQP_Buff3 <- grep(
  pattern = "__BQP-.*-3__", x = names(AllSppData))
for(i in Index_BQP_Buff3){
  AllSppData[[i]]$BQP_Buff <- 3
}


2.3. BQP response curves options

BQP_CurveNum1 <- grep(
  pattern = "BQP-WO", x = names(AllSppData))
for(i in BQP_CurveNum1){
  AllSppData[[i]]$BQP_CurveNum <- 1
}

BQP_CurveNum2 <- grep(
  pattern = "BQP-Low", x = names(AllSppData))
for(i in BQP_CurveNum2){
  AllSppData[[i]]$BQP_CurveNum <- 2
}

BQP_CurveNum3 <- grep(
  pattern = "BQP-Med", x = names(AllSppData))
for(i in BQP_CurveNum3){
  AllSppData[[i]]$BQP_CurveNum <- 3
}

BQP_CurveNum4 <- grep(
  pattern = "BQP-Strng", x = names(AllSppData))
for(i in BQP_CurveNum4){
  AllSppData[[i]]$BQP_CurveNum <- 4
}


2.4. Species options

In this example, there are 5 species for each group. Species 1:5 for reptiles; Species 6:10 for butterflies; 11:15 for mammals


2.4.1. All species

SpNums <- 1:15

#-------------------------

FileNames_AllSp_Mxnt_NoBias <- paste0(
  "ZigInputFiles/Maps/Mxnt_NoBias/Mxnt_NoBias_Sp",
  SpNums, ".tif")
Index_AllSP_Mxnt_NoBias <- grep(
  pattern = "Mxnt-BiasWO-AllSP",
  x = names(AllSppData))
for (i in Index_AllSP_Mxnt_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_AllSp_Mxnt_NoBias
}

#-------------------------

FileNames_AllSp_Mxnt_Bias0 <- paste0(
  "ZigInputFiles/Maps/Mxnt_Bias0/Mxnt_Bias0_Sp",
  SpNums, ".tif")
Index_AllSP_Mxnt_Bias0 <- grep(
  pattern = "Mxnt-Bias0-AllSP",
  x = names(AllSppData))
for (i in Index_AllSP_Mxnt_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_AllSp_Mxnt_Bias0
}

#-------------------------

FileNames_AllSp_En_NoBias <- paste0(
  "ZigInputFiles/Maps/En_NoBias/En_NoBias_Sp",
  SpNums, ".tif")
Index_AllSP_En_NoBias <- grep(
  pattern = "EN-BiasWO-AllSP",
  x = names(AllSppData))
for (i in Index_AllSP_En_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_AllSp_En_NoBias
}

#-------------------------

FileNames_AllSp_En_Bias0 <- paste0(
  "ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp",
  SpNums, ".tif")
Index_AllSP_En_Bias0 <- grep(
  pattern = "EN-Bias0-AllSP",
  x = names(AllSppData))
for (i in Index_AllSP_En_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_AllSp_En_Bias0
}


2.4.2. Reptiles

SpNums <- 1:5

#-------------------------

FileNames_REPT_Mxnt_NoBias <- paste0(
  "ZigInputFiles/Maps/Mxnt_NoBias/Mxnt_NoBias_Sp",
  SpNums, ".tif")
Index_REPT_Mxnt_NoBias <- grep(
  pattern = "Mxnt-BiasWO-Rep",
  x = names(AllSppData))
for (i in Index_REPT_Mxnt_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_REPT_Mxnt_NoBias
}

#-------------------------

FileNames_REPT_Mxnt_Bias0 <- paste0(
  "ZigInputFiles/Maps/Mxnt_Bias0/Mxnt_Bias0_Sp",
  SpNums, ".tif")
Index_REPT_Mxnt_Bias0 <- grep(
  pattern = "Mxnt-Bias0-Rep",
  x = names(AllSppData))
for (i in Index_REPT_Mxnt_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_REPT_Mxnt_Bias0
}

#-------------------------

FileNames_REPT_En_NoBias <- paste0(
  "ZigInputFiles/Maps/En_NoBias/En_NoBias_Sp",
  SpNums, ".tif")
Index_REPT_En_NoBias <- grep(
  pattern = "EN-BiasWO-Rep",
  x = names(AllSppData))
for (i in Index_REPT_En_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_REPT_En_NoBias
}

#-------------------------

FileNames_REPT_En_Bias0 <- paste0(
  "ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp",
  SpNums, ".tif")
Index_REPT_En_Bias0 <- grep(
  pattern = "EN-Bias0-Rep",
  x = names(AllSppData))
for (i in Index_REPT_En_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_REPT_En_Bias0
}


2.4.3. Butterflies

SpNums <- 6:10

#-------------------------

FileNames_BTTR_Mxnt_NoBias <- paste0(
  "ZigInputFiles/Maps/Mxnt_NoBias/Mxnt_NoBias_Sp",
  SpNums, ".tif")
Index_BTTR_Mxnt_NoBias <- grep(
  pattern = "Mxnt-BiasWO-Butr",
  x = names(AllSppData))
for (i in Index_BTTR_Mxnt_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_BTTR_Mxnt_NoBias
}

#-------------------------

FileNames_BTTR_Mxnt_Bias0 <- paste0(
  "ZigInputFiles/Maps/Mxnt_Bias0/Mxnt_Bias0_Sp",
  SpNums, ".tif")
Index_BTTR_Mxnt_Bias0 <- grep(
  pattern = "Mxnt-Bias0-Butr",
  x = names(AllSppData))
for (i in Index_BTTR_Mxnt_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_BTTR_Mxnt_Bias0
}

#-------------------------

FileNames_BTTR_En_NoBias <- paste0(
  "ZigInputFiles/Maps/En_NoBias/En_NoBias_Sp",
  SpNums, ".tif")
Index_BTTR_En_NoBias <- grep(
  pattern = "EN-BiasWO-Butr",
  x = names(AllSppData))
for (i in Index_BTTR_En_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_BTTR_En_NoBias
}

#-------------------------

FileNames_BTTR_En_Bias0 <- paste0(
  "ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp",
  SpNums, ".tif")
Index_BTTR_En_Bias0 <- grep(
  pattern = "EN-Bias0-Butr",
  x = names(AllSppData))
for (i in Index_BTTR_En_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_BTTR_En_Bias0
}


2.4.4. Mammals

SpNums <- 11:15

#-------------------------

FileNames_MMLS_Mxnt_NoBias <- paste0(
  "ZigInputFiles/Maps/Mxnt_NoBias/Mxnt_NoBias_Sp",
  SpNums, ".tif")
Index_MMLS_Mxnt_NoBias <- grep(
  pattern = "Mxnt-BiasWO-MMls",
  x = names(AllSppData))
for (i in Index_MMLS_Mxnt_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_MMLS_Mxnt_NoBias
}

#-------------------------

FileNames_MMLS_Mxnt_Bias0 <- paste0(
  "ZigInputFiles/Maps/Mxnt_Bias0/Mxnt_Bias0_Sp",
  SpNums, ".tif")
Index_MMLS_Mxnt_Bias0 <- grep(
  pattern = "Mxnt-Bias0-MMls",
  x = names(AllSppData))
for (i in Index_MMLS_Mxnt_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_MMLS_Mxnt_Bias0
}

#-------------------------

FileNames_MMLS_En_NoBias <- paste0(
  "ZigInputFiles/Maps/En_NoBias/En_NoBias_Sp",
  SpNums, ".tif")
Index_MMLS_En_NoBias <- grep(
  pattern = "EN-BiasWO-MMls",
  x = names(AllSppData))
for (i in Index_MMLS_En_NoBias){
  AllSppData[[i]]$SpeciesMaps <- FileNames_MMLS_En_NoBias
}

#-------------------------

FileNames_MMLS_En_Bias0 <- paste0(
  "ZigInputFiles/Maps/En_Bias0/En_Bias0_Sp",
  SpNums, ".tif")
Index_MMLS_En_Bias0 <- grep(
  pattern = "EN-Bias0-MMls",
  x = names(AllSppData))
for (i in Index_MMLS_En_Bias0){
  AllSppData[[i]]$SpeciesMaps <- FileNames_MMLS_En_Bias0
}


2.5. Weight options

These files contain the different values of per-species weight

WeightsData_EN_Bias0 <- read.csv(
  "ZigInputFiles/Data/WeightsData_EN_Bias0.csv",
  sep = ",", header = T)

WeightsData_EN_WOBias <- read.csv(
  "ZigInputFiles/Data/WeightsData_EN_WOBias.csv",
  sep = ",", header = T)

WeightsData_Mxnt_Bias0 <- read.csv(
  "ZigInputFiles/Data/WeightsData_Mxnt_Bias0.csv",
  sep = ",", header = T)

WeightsData_Mxnt_WOBias <- read.csv(
  "ZigInputFiles/Data/WeightsData_Mxnt_WOBias.csv",
  sep = ",", header = T)


2.5.1. No predictive uncertainty & No Red-List weight

Index_WTWO_NoUncert <- grep(
  pattern = "WTWO-NoUncert",
  x = names(AllSppData))

for (i in Index_WTWO_NoUncert){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_WO,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <-
        WeightsData_Mxnt_WOBias$WT_WO[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <-
          WeightsData_Mxnt_WOBias$WT_WO[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <-
            WeightsData_Mxnt_WOBias$WT_WO[11:15],
          print("ERROR")))))
}


2.5.2. No predictive uncertainty & with Red-List weight 

Index_WTLoc_NoUncert <- grep(
  pattern = "WTLoc-NoUncert",
  x = names(AllSppData))

for (i in Index_WTLoc_NoUncert){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <-
      WeightsData_Mxnt_WOBias$WT_Local,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <-
        WeightsData_Mxnt_WOBias$WT_Local[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <-
          WeightsData_Mxnt_WOBias$WT_Local[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <-
            WeightsData_Mxnt_WOBias$WT_Local[11:15],
          print("ERROR")))))
}


2.5.3. With predictive uncertainty & No Red-List weight

2.5.3.1. With no sampling bias | Maxent

Index_WTWO_Uncert_Mxnt_BiasWO <- grep(
  pattern = "WTWO-Uncert.*Mxnt-BiasWO",
  x = names(AllSppData))

for (i in Index_WTWO_Uncert_Mxnt_BiasWO){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_WO_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_WO_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_WO_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_WO_Uncert[11:15],
          print("ERROR")))))
}

2.5.3.2. With no sampling bias | Elastic net

Index_WTWO_Uncert_EN_BiasWO <- grep(
  pattern = "WTWO-Uncert.*EN-BiasWO",
  x = names(AllSppData))

for (i in Index_WTWO_Uncert_EN_BiasWO){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_WO_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_WO_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_WO_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_WO_Uncert[11:15],
          print("ERROR")))))
}

2.5.3.3. With correction for sampling bias | Maxent

Index_WTWO_Uncert_Mxnt_Bias0 <- grep(
  pattern = "WTWO-Uncert.*Mxnt-Bias0",
  x = names(AllSppData))

for (i in Index_WTWO_Uncert_Mxnt_Bias0){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_WO_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_WO_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_WO_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_WO_Uncert[11:15],
          print("ERROR")))))
}

2.5.3.4. With correction for sampling bias | Elastic net

Index_WTWO_Uncert_EN_Bias0 <- grep(
  pattern = "WTWO-Uncert.*EN-Bias0",
  x = names(AllSppData))

for (i in Index_WTWO_Uncert_EN_Bias0){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_WO_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_WO_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_WO_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_WO_Uncert[11:15],
          print("ERROR")))))
}


2.5.4. With predictive uncertainty & Red-List weight

2.5.4.1. With no sampling bias | Maxent

Index_WTLoc_Uncert_Mxnt_BiasWO <- grep(
  pattern = "WTLoc-Uncert.*Mxnt-BiasWO",
  x = names(AllSppData))

for (i in Index_WTLoc_Uncert_Mxnt_BiasWO){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_Local_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_Local_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_Local_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_Mxnt_WOBias$WT_Local_Uncert[11:15],
          print("ERROR")))))
}

2.5.4.2. With no sampling bias | Elastic net

Index_WTLoc_Uncert_EN_BiasWO <- grep(
  pattern = "WTLoc-Uncert.*EN-BiasWO",
  x = names(AllSppData))

for (i in Index_WTLoc_Uncert_EN_BiasWO){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_Local_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_Local_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_Local_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_EN_WOBias$WT_Local_Uncert[11:15],
          print("ERROR")))))
}

2.5.4.3. With correction for sampling bias | Maxent

Index_WTLoc_Uncert_Mxnt_Bias0 <- grep(
  pattern = "WTLoc-Uncert.*Mxnt-Bias0",
  x = names(AllSppData))

for (i in Index_WTLoc_Uncert_Mxnt_Bias0){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_Local_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_Local_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_Local_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_Mxnt_Bias0$WT_Local_Uncert[11:15],
          print("ERROR")))))
}

2.5.4.4. With correction for sampling bias | Elastic net

Index_WTLoc_Uncert_EN_Bias0 <- grep(
  pattern = "WTLoc-Uncert.*EN-Bias0",
  x = names(AllSppData))

for (i in Index_WTLoc_Uncert_EN_Bias0){
  ifelse(
    grepl(pattern = "AllSP", x = names(AllSppData[i])),
    AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_Local_Uncert,
    
    ifelse(
      grepl(pattern = "Rep", x = names(AllSppData[i])),
      AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_Local_Uncert[1:5],
      
      ifelse(
        grepl(pattern = "Butr", x = names(AllSppData[i])),
        AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_Local_Uncert[6:10],
        
        ifelse(
          grepl(pattern = "MMls", x = names(AllSppData[i])),
          AllSppData[[i]]$Weight <- WeightsData_EN_Bias0$WT_Local_Uncert[11:15],
          print("ERROR")))))
}


2.6. Exporting .spp files

for (i in 1:length(AllSppData)){
  write.table(
    x = AllSppData[[i]],
    file = paste0("output_spp/", names(AllSppData[i])),
    sep = "\t", row.names = FALSE, col.names = FALSE)
}


3. The run settings file (.dat)

8 files are available in ‘Maps’ folder (they were written outside R). An example file is shown above.


4. Preparing project files [.sh or bat files]

The project file specifies two main input files and the output files.


4.1. Centos Linux cluster

Example .sh file:

#!/bin/bash
zig4 -r "ABF_Mask.dat" "WTLoc-NoUncert__BQP-Low-1__EN-Bias0-AllSP.spp" "output/ABF_Mask__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-AllSP.txt" 0 1 1 1 --use-threads=16


  • These .sh files were prepared to work in a Centos Linux cluster servers; using Slurm Workload Manager.
  • This code produce one .sh file for each combination of spp and dat file will be prepared.
  • the argument --use-threads in the next code chunk represents the number of parallel computaions of Zonation.
SppFiles <- list.files(
  path = "output_spp\\",
  pattern = ".*.spp", full.names = FALSE)

SppFilesInit <- gsub(
  pattern = ".spp", replacement = "", x = SppFiles)

datFiles <- list.files(
  path = "ZigInputFiles\\Dat\\",
  pattern = ".*.dat", full.names=FALSE)

datFilesInit <- gsub(
  pattern = ".dat", replacement = "", x = datFiles)

AllShFiles <- expand.grid(
  dat = datFilesInit, spp = SppFilesInit,
  stringsAsFactors = FALSE)

AllShFiles$sh <- paste0(
  AllShFiles$dat, "__", AllShFiles$spp, ".sh")

AllShFiles$OutFile <- paste0(
  "output/", AllShFiles$dat, "__",
  AllShFiles$spp, ".txt")

AllShFiles$dat <- paste0(AllShFiles$dat, ".dat")
AllShFiles$spp <- paste0(AllShFiles$spp, ".spp")

for (i in 1:nrow(AllShFiles)){
  Currdat <- AllShFiles[i,"dat"]
  CurrSpp <- AllShFiles[i,"spp"]
  CurrOutSh <- AllShFiles[i,"sh"]
  CurrOutTxt <- AllShFiles[i,"OutFile"]
  
  cat(
    paste0(
      '#!/bin/bash\nzig4 -r "',
      Currdat, '" "', CurrSpp, '" "',
      CurrOutTxt, '" ', '0 1 1 1 --use-threads=16'),
    file = paste0("output_sh\\", CurrOutSh))
  rm(Currdat, CurrSpp, CurrOutSh, CurrOutTxt)
}

Now, there should be 2,560 .sh files:

  • 4 surrogate taxa (butterflies vs reptiles vs mammals vs the three groups together) ×
  • 2 modelling algorithms (Maxent vs elastic net) ×
  • 2 biases (without vs with sampling-bias correction) ×
  • 2 cell-removal rules (ABF vs CAZ) ×
  • 4 weighting schemes (with vs without Red-List or predictive-uncertainty weights) ×
  • 10 connectivity options ×
  • 2 PA-integration options (without vs with PA masking)


4.2. Windows PC

Example .bat file:

call zig4.exe -r "ABF_Mask.dat" "WTLoc-NoUncert__BQP-Low-1__EN-Bias0-AllSP.spp" "output/ABF_Mask__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-AllSP.txt" 0 1 1 1 --use-threads=16


This is a similar code to produce .bat files to work under windows environment

AllShFiles$bat <- paste0(
  gsub(pattern = ".dat", replacement = "",
       x = AllShFiles$dat),
  "__",
  gsub(pattern = ".spp", replacement = "",
       x = AllShFiles$spp),
  ".bat")

for (i in 1:nrow(AllShFiles)){
  Currdat <- AllShFiles[i,"dat"]
  CurrSpp <- AllShFiles[i,"spp"]
  CurrOutBat <- AllShFiles[i,"bat"]
  CurrOutTxt <- AllShFiles[i,"OutFile"]
  
  cat(
    paste0(
      'call zig4.exe -r "', Currdat, '" "', CurrSpp,
      '" "', CurrOutTxt, '" ', '0 1 1 1 --use-threads=16'),
    file = paste0("output_bat\\", CurrOutBat))
  
  rm(Currdat, CurrSpp, CurrOutBat, CurrOutTxt)
}


5. Running analyses on parallel

This code creates .moab files to facilitate the parallel run of 2,560 Zonation runs on a Centos Linux cluster

An example .moab file:

#!/bin/sh
########### Begin MOAB/Slurm header ##########
#
# Give job a reasonable name
#MOAB -N Task_1
#
# Request number of nodes and CPU cores per node for job
#MOAB -l nodes=1:ppn=16, pmem=48gb
#
# Estimated wallclock time for job
#MOAB -l walltime=00:12:00:00
#
# Write standard output and errors in same file
#MOAB -j oe
#
# Send mail when job begins, aborts and ends
#MOAB -m bae
#
########### End MOAB header ##########

echo "Working Directory:                    $PWD"
echo "Running on host                       $HOSTNAME"
echo "Job id:                               $MOAB_JOBID"
echo "Job name:                             $MOAB_JOBNAME"
echo "Number of nodes allocated to job:     $MOAB_NODECOUNT"
echo "Number of cores allocated to job:     $MOAB_PROCCOUNT"

module load geo/zonation/4.0.0
cd ZonationEgypt/

bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-AllSP.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-Butr.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-MMls.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-Bias0-Rep.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-BiasWO-AllSP.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-BiasWO-Butr.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-BiasWO-MMls.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__EN-BiasWO-Rep.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-Bias0-AllSP.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-Bias0-Butr.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-Bias0-MMls.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-Bias0-Rep.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-BiasWO-AllSP.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-BiasWO-Butr.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-BiasWO-MMls.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-1__Mxnt-BiasWO-Rep.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-2__EN-Bias0-AllSP.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-2__EN-Bias0-Butr.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-2__EN-Bias0-MMls.sh" &
bash "ABF_MaskNo__WTLoc-NoUncert__BQP-Low-2__EN-Bias0-Rep.sh" &
wait
echo "All are complete"


5.1. Analyses without PAs masking

This code split 1280 Zonation runs into 64 files, each run 20 Zonation run in parallel

FirstRuns <- list.files(
  path = "output_sh\\", pattern = ".*.sh"
)

FirstRuns <- FirstRuns[grep(pattern = ".*_MaskNo__",
                            x = FirstRuns)]

FirstRuns1 <- split(
  FirstRuns, ceiling(seq_along(1:length(FirstRuns))/20))

names(FirstRuns1) <- paste0("Task_", 1:64, ".moab")

for (ii in 1:length(FirstRuns1)){
  CurrFile <- paste0(
    "output_moab\\", names(FirstRuns1)[ii])
  
  cat("#!/bin/sh\n", file = CurrFile, append = FALSE)
  
  cat(paste0(
    "########### Begin MOAB/Slurm header ##########\n#\n# Give job a reasonable name\n#MOAB -N Task_", 0+ii),
    file = CurrFile, append = TRUE)
  
  cat(paste0(
    '\n#\n# Request number of nodes and CPU cores per node for job\n#MOAB -l nodes=1:ppn=16, pmem=48gb\n#\n# Estimated wallclock time for job\n#MOAB -l walltime=00:12:00:00\n#\n# Write standard output and errors in same file\n#MOAB -j oe\n#\n# Send mail when job begins, aborts and ends\n#MOAB -m bae\n#\n########### End MOAB header ##########\n\necho "Working Directory:                    $PWD"\necho "Running on host                       $HOSTNAME"\necho "Job id:                               $MOAB_JOBID"\necho "Job name:                             $MOAB_JOBNAME"\necho "Number of nodes allocated to job:     $MOAB_NODECOUNT"\necho "Number of cores allocated to job:     $MOAB_PROCCOUNT"\n\nmodule load geo/zonation/4.0.0\ncd ZonationEgypt/\n\n'),
    file = CurrFile, append = TRUE)
  
  for(i in 1:length(FirstRuns1[[ii]])){
    cat(paste0(
      'bash "', FirstRuns1[[ii]][i], '" &\n'),
      file = CurrFile, append = TRUE)
  }
  
  cat(paste0('wait\necho "All are complete"'),
      file = CurrFile, append = TRUE)
}

Tasks <- paste0("msub Task_", 1:128, ".moab")
cat(Tasks, file = "output_tasks\\Tasks1_128.txt", sep = "\n")


5.2. Analyses with PAs masking

FirstRuns <- list.files(
  path = "output_sh\\", pattern = ".*.sh")

FirstRuns <- FirstRuns[grep(pattern = ".*_Mask__",
                            x = FirstRuns)]

FirstRuns1 <- split(
  FirstRuns, ceiling(seq_along(1:length(FirstRuns))/20))

names(FirstRuns1) <- paste0("Task_", 65:128, ".moab")

for (ii in 1:length(FirstRuns1)){
  CurrFile <- paste0(
    "output_moab\\", names(FirstRuns1)[ii])
  
  cat("#!/bin/sh\n", file = CurrFile, append = FALSE)
  
  cat(paste0(
    "########### Begin MOAB/Slurm header ##########\n#\n# Give job a reasonable name\n#MOAB -N Task_", 64+ii),
    file = CurrFile, append = TRUE)
  
  cat(paste0(
    '\n#\n# Request number of nodes and CPU cores per node for job\n#MOAB -l nodes=1:ppn=16, pmem=48gb\n#\n# Estimated wallclock time for job\n#MOAB -l walltime=00:12:00:00\n#\n# Write standard output and errors in same file\n#MOAB -j oe\n#\n# Send mail when job begins, aborts and ends\n#MOAB -m bae\n#\n########### End MOAB header ##########\n\necho "Working Directory:                    $PWD"\necho "Running on host                       $HOSTNAME"\necho "Job id:                               $MOAB_JOBID"\necho "Job name:                             $MOAB_JOBNAME"\necho "Number of nodes allocated to job:     $MOAB_NODECOUNT"\necho "Number of cores allocated to job:     $MOAB_PROCCOUNT"\n\nmodule load geo/zonation/4.0.0\ncd ZonationEgypt/\n\n'),
    file = CurrFile, append = TRUE)
  
  for(i in 1:length(FirstRuns1[[ii]])){
    cat(paste0(
      'bash "', FirstRuns1[[ii]][i], '" &\n'),
      file = CurrFile, append = TRUE)
  }
  
  cat(paste0('wait\necho "All are complete"'),
      file = CurrFile, append = TRUE)
}

Tasks <- paste0("msub Task_", 65:128, ".moab")
cat(Tasks, file = "output_tasks\\Tasks65_128.txt", sep = "\n")