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 theZonation
software can be found in Di Minin et al., 2014. I expect the reader to be familiar withZonation
and have read the paper in advance.This html file was written in
RMarkdown
. TheRMarkdown
and sample data are available here.
1.1. Abbreviations used
Species weight
WTLoc
: weighted by Egyptian national red-list statusWTWO
: without red-list weightingUncert
: weighted by predictive uncertaintyNoUncert
: without predictive uncertainty weight
Connectivity
BQP
: boundary quality penaltyBQP-WO
: no connectivity analysesBQP-Low
/BQP-Med
/BQP-Strng
: low/medium/strong BQP curves
Modelling algorithms
Mxnt
: MaxentEN
: elastic net
Sampling bias
BiasWO
: environment-only models (no bias correction)Bias0
: bias-free predictions, implementing model-based bias correction
Surrogate groups
MMls
: MammalsRep
: ReptilesButr
: ButterfliesAllSP
: All three groups together
1.2. Input files
Species weights
- Four
*.csv
files for different weighting options are available in the folderZigInputFiles/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 weightingWTLoc-NoUncert
: weighted by the Egyptian national red-list status onlyWTWO-Uncert
: weighted by predictive uncertainty onlyWTLoc-Uncert
: weighting by both
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")))))
}
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 ofspp
anddat
file will be prepared. - the argument
--use-threads
in the next code chunk represents the number of parallel computaions ofZonation
.
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")