# Overview

This page provides guidance on setting conservation targets and some examples of post-hoc evaluation of network design

The plots shown in this tutorial were created in R using output generated by Marxan Connect. In this example, we automate the iterative runs of Marxan using R to facilitate the process.

require(tidyverse)
require(sf)

# Creating an example connectivity Matrix

Here we create a connectivity matrix which contains flow data. This step is not necessary if you have your own data.

# load the reefs shapefile and calculate area
mutate(area = st_area(.))
## Reading layer reefs' from data source C:\Users\daigl\Documents\GitHub\MarxanConnect\docs\tutorial\targets\reefs.shp' using driver ESRI Shapefile'
## Simple feature collection with 321 features and 15 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -548876.9 ymin: -2104102 xmax: -215060 ymax: -1547444
## projected CRS:  Equidistant_Conic
# generate a distance matrix
distanceMatrix <- round(st_distance(reefs))

# calculate isolation by distance
isolationMatrix <- 1/distanceMatrix^2 %>% matrix(nrow(reefs),nrow(reefs))

# assume dispersal limit is 200km, remove links >200km
isolationMatrix[isolationMatrix<1/200000^2]=0

# create rudimentary probability matrix
probabilityMatrix <- isolationMatrix

probabilityMatrix[is.infinite(probabilityMatrix)] <- max(probabilityMatrix[is.finite(probabilityMatrix)],na.rm = TRUE)/10^5
diag(probabilityMatrix) <- diag(probabilityMatrix)*sqrt(reefs$area)/mean(sqrt(reefs$area))

# row normalize to make it a probability matrix
probabilityMatrix <- probabilityMatrix/rowSums(probabilityMatrix)

# create flow matrix that includes fecundity and survival
fecundity <- as.numeric(reefs$area) survival <- as.numeric(reefs$area/max(reefs$area))*0.0000001 flowMatrix <- t(t(probabilityMatrix)*survival)*fecundity # write to file write.csv(flowMatrix,"tutorial/targets/reefFlow.csv") # Using Marxan Connect In this example we used the above flow matrix to generate conservation features in Marxan Connect. We appended those conservation features to the existing Marxan files. To replicate this process either follow the “Demographic Data using CF” tutorial or load and modify the targettradeoff.MarCon project file. # Iterating Marxan from R Using R to run Marxan allows us to run Marxan including 1 connectivity related conservation feature at a time, and with varying values of the Connectivity target multiplier (C). C is a tunable ‘contraint’ which scales the targets for connectivity based conservaton features relative to that of normal conservation features. First let’s prepare all the parameters. # set targets for non-connectivity related targets regulartarget <- 0.1 # number of replicates in a Marxan run reps <- 20 # metrics over which to iterate metrics <- c("in_degree_demo_pu", "out_degree_demo_pu", "between_cent_demo_pu", "eig_vect_cent_demo_pu", "google_demo_pu", "self_recruit_demo_pu", "local_retention_demo_pu", "inflow_demo_pu", "outflow_demo_pu", "fa_recipients_demo_pu", "fa_donors_demo_pu", "aa_donors_demo_pu", "aa_recipients_demo_pu") # values of the connectivity multiplier over which to iterate steps <- 1 # connectivity multiplier step size cmax <- 1/regulartarget # maximum connectivity multiplier # get 100% target values # load puvspr puvspr <- read.csv("tutorial/targets/input/puvspr2_appended.dat",stringsAsFactors=FALSE) %>% group_by(species) %>% summarize(target=sum(amount)) %>% ungroup() %>% data.frame() ## summarise() ungrouping output (override with .groups argument) Now we can iterate Marxan. If you’re using the the RMarkdown .Rmd version of this tutorial, you may want to remove the eval=FALSE in the R chunk options since this chunk of code won’t run if it’s still there. It’s there to prevent the lengthy iterations of Marxan when building the website! # delete previous ouputs unlink(x=list.files(path="tutorial/targets/output/",pattern="*",full.names = TRUE)) #iterate over metrics and connectivity multiplier for(m in metrics){ for (c in seq(0,cmax,steps)) { # load specdat specdat <- read.csv("tutorial/targets/input/spec_appended.dat",stringsAsFactors=FALSE) # load old specdat to get names of non-connectivity features oldspecdat <- read.csv("tutorial/targets/input/spec.dat",stringsAsFactors=FALSE) # target for all features # replace 100% targets if(all(specdat$id==puvspr$species)){ specdat$target <- puvspr$target } else { print("things are not lined up properly!") } specdat$target <- specdat$target*regulartarget specdat$target[specdat$name==m] <- specdat$target[specdat$name==m]*c specdat$target[!(specdat$name %in% oldspecdat$name) & specdat$name!=m] <- 0 # save spec file write.csv(specdat,paste0("tutorial/targets/input/c_testingTargetScenario_",m,"_",c,".dat"),quote=FALSE,row.names=FALSE) # load input.dat input.dat <- readLines("tutorial/targets/input.dat") # set parameter iBLMParam <- which(regexpr("BLM",input.dat)==1) iNUMREPSParam <- which(regexpr("NUMREPS",input.dat)==1) iSPECNAMEParam <- which(regexpr("SPECNAME",input.dat)==1) iPUNAMEParam <- which(regexpr("PUNAME",input.dat)==1) iSCENNAMEParam <- which(regexpr("SCENNAME",input.dat)==1) iOUTPUTDIRParam <- which(regexpr("OUTPUTDIR",input.dat)==1) iINPUTDIRParam <- which(regexpr("INPUTDIR",input.dat)==1) input.dat[iBLMParam] <- paste0("BLM ",1) input.dat[iNUMREPSParam] <- paste0("NUMREPS ",reps) input.dat[iSPECNAMEParam] <- paste0("SPECNAME c_testingTargetScenario_",m,"_",c,".dat") input.dat[iPUNAMEParam] <- paste0("PUNAME pu.dat") input.dat[iSCENNAMEParam] <- paste0("SCENNAME c_testingTargetScenario_",m,"_",c) input.dat[iOUTPUTDIRParam] <- gsub("/","\\\\",paste0("OUTPUTDIR ",getwd(),"/tutorial/targets/output")) input.dat[iINPUTDIRParam] <- gsub("/","\\\\",paste0("INPUTDIR ",getwd(),"/tutorial/targets/input")) input.dat[length(input.dat)+1] <- "" # save input.dat writeLines(input.dat,con=file(paste0("tutorial/targets/input",c,".dat")),sep="\n") # run Marxan # TODO: you will need to change to the appropriate Marxan filepath here system(paste0("../Marxan243/Marxan_x64.exe input",c,".dat")) } } # Processing output and plotting The hundreds of Marxan output files must now be loaded, processed and plotted. Here we calculate the ecologically relevant # create table bestfiles <- list.files(path="tutorial/targets/output",pattern="*mvbest*") sumfiles <- list.files(path="tutorial/targets/output",pattern="*sum*") #We want all the cost of every best solution for each scenario #We want the metric value for every best solution for each scenario # create a table to fill in below mastertable <- expand.grid(metric=metrics,c=seq(0,cmax,steps),replicate=sprintf("%05d", 1:reps),Cost=0,Metric=0,"Metapopulation Growth"=0) # load specdat specdat <- read.csv("tutorial/targets/input/spec_appended.dat",stringsAsFactors=FALSE) # iterate over rows in mastertable for(i in 1:nrow(mastertable)){ # extract the amount of the metric held in the network mastertable$Metric[i] <- read.csv(file = paste0("tutorial/targets/output/c_testingTargetScenario_",mastertable$metric[i],"_",mastertable$c[i],"_mv",mastertable$replicate[i],".txt"), header = T) %>% filter(Conservation.Feature==specdat$id[specdat$name==mastertable$metric[i]]) %>%
select(Amount.Held) %>%
unlist() %>%
as.numeric()

# extract the cost of the network
mastertable$Cost[i] <- read.csv(file = paste0("tutorial/targets/output/c_testingTargetScenario_",mastertable$metric[i],"_",mastertable$c[i],"_sum.txt"), header = T) %>% filter(Run_Number==as.numeric(mastertable$replicate[i])) %>%
select(Cost) %>%
unlist()

# Calculate an ecologically relevant conservation objective
# we calculate the metapopulation growth (first real eigenvalue of a connectivity matrix that includes survival and mortality)
# more details on this step in the text below
inreserve <- read.csv(file = paste0("tutorial/targets/output/c_testingTargetScenario_",mastertable$metric[i],"_",mastertable$c[i],"_r",mastertable$replicate[i],".txt"), header = T)$solution==1
networkmatrix <- (flowMatrix)[inreserve,inreserve]
mastertable$"Metapopulation Growth"[i] <- Re(eigen(networkmatrix)$values)

}

That ecologically relevant conservation objective step is very important. Here we calculated the metapopulation growth of a flow matrix that includes survival and mortality using the first eigenvalue of the matrix. Our conservation objective in this case is to have a self-sustaining metapopulation inside the reserve network (i.e. in-reserve network metapopulation growth rate > 1). Please note that this is not the only possible conservation objective, or even the only way to calculate metapopulation growth. This example is for demonstration purposes only. Please see Shima, Noonburg, and Phillips (2010),Drechsler (2009),Kininmonth et al. (2010),Day and Possingham (1995),Figueira and Crowder (2006),Hanski and Ovaskainen (2000),Etienne (2004) for more information on ecologically relevant conservation objectives.

# make it nice! Give the metrics "real" names for plotting
mastertable$metric <- recode(mastertable$metric,
"in_degree_demo_pu" = "In Degree",
"out_degree_demo_pu" = "Out Degree",
"between_cent_demo_pu" = "Betweenness Centrality",
"eig_vect_cent_demo_pu" = "Eigenvector Centrality",
"self_recruit_demo_pu" = "Self Recruitment",
"local_retention_demo_pu" = "Local Rentention",
"inflow_demo_pu" = "In-flow",
"outflow_demo_pu" = "Out-flow",
"fa_recipients_demo_pu" = "Focus Area Recipient",
"fa_donors_demo_pu" = "Focus Area Donor",
"aa_recipients_demo_pu" = "Avoidance Area Recipient",
"aa_donors_demo_pu" = "Avoidance Area Donor")

# calculate the mean and SE, and convert all the values to a relative scale
# (100 = maximum value, 0 = value for network without connectivity based conservation feature)
mastertableMeans <- mastertable %>%
gather(key = "Type", value = "value",-metric,-c,-replicate)%>%
group_by(metric,c,Type) %>%
summarize(means=mean(value),
se=sd(value)/sqrt(n())) %>%
ungroup() %>%
group_by(metric,Type) %>%
mutate(relativeSE = se/(max(means)-means[c==0]*100),
relativeMean = (means-means[c==0])/(max(means)-means[c==0])*100) %>%
ungroup() %>%
dplyr::filter(metric == "In Degree" |
metric == "Betweenness Centrality"|
metric == "Avoidance Area Recipient"
)
## summarise() regrouping output by 'metric', 'c' (override with .groups argument)
# back-calculate what a metapopulation growth of 1 would be on the plot's relative scale
lamdbaM1 <- mean((1-mastertableMeans$means[mastertableMeans$c==0&mastertableMeans$Type=="Metapopulation Growth"])/ (max(mastertableMeans$means[mastertableMeans$Type=="Metapopulation Growth"])-mastertableMeans$means[mastertableMeans$c==0&mastertableMeans$Type=="Metapopulation Growth"]))*100

And finally, let’s see the results!

p <- ggplot(mastertableMeans[mastertableMeans\$c!=0,],aes(colour=Type,fill=Type,shape=Type,linetype=Type,x=c,y=relativeMean))+
geom_point(size=2)+
geom_line(size=1)+
geom_ribbon(aes(ymin=relativeMean-relativeSE,ymax=relativeMean+relativeSE),alpha=0.3,colour="transparent")+
geom_hline(yintercept=0,linetype="dashed",size=0.75) +
geom_hline(yintercept=lamdbaM1,linetype="dashed",size=0.75,colour="#d95f02") +
scale_linetype_manual(values=c("solid","solid","solid"))+
scale_color_manual(values=c("#1b9e77","#d95f02","#7570b3"))+
scale_fill_manual(values=c("#1b9e77","#d95f02","#7570b3"))+
labs(x="Connectivity Target Multiplier (C)", y = "Relative Value (%)")+
scale_x_continuous(breaks = seq(1,cmax,2))+
facet_wrap(~metric) +
theme_classic() +
theme(strip.background = element_blank(), axis.line = element_line(),legend.position="top")
p Optional figure saving:

ggsave("tutorial/targets/C_Multiplier.eps",width = 180, height = 180, units = c("mm"))
ggsave("tutorial/targets/C_Multiplier.png",width = 180, height = 180, units = c("mm"))`