library(knitr) opts_chunk$set(fig.width=7, fig.height=7,cache=TRUE)
Thank you for your interest in the R package for the bioeconomic evaluation of MPA network design! If you haven't already done so, please refer to the manuscript describing the model results.
If you would like to use this package to your own scenario, please fork this repository and modify as needed, or install the package directly from GitHub.
install.packages("devtools") library(devtools) install_github("remi-daigle/BESTMPA@gh-pages")
We kindly ask that you cite both the original manuscript and this repository in any resulting publication. We recognize that there are some aspects of the model the could be closer to reality and we look forward to seeing how you would improve it. To make changes, we recommend starting with the code in this vignette (or the index.Rmd)
If you would like to replicate our results, download the repository as is and run the simulations by running the code in index.Rmd. This script will contains all the user defined parameters, and the functions are all contained within the R/ directory if you would like to inspect and or modify them individually.
Please note: This package requires the installation of R packages before use (tidyverse, data.table, ggplot2 grid, Grid2Polygons, gridExtra, igraph, maptools, Matrix, raster, readr, rgdal*, rgeos, sp, spdep, tidyr)
*installation of the rgdal package is notoriously difficult, please see here and here for further information should you have any problems.
require(BESTMPA) require(rgdal) require(rgeos) require(raster) require(igraph) require(ggplot2) require(tidyverse) require(data.table) require(gridExtra) require(grid) require(spdep) require(Matrix) require(Grid2Polygons) require(plotrix)
Before getting started, using the housekeeping function will prepare the environment and your file structure for modelling. This function will conveniently create a results folder, clear the environment, figures, delete previous results, and increase the memory limit if desired.
results_folder <- "BESTMPA_results" housekeeping(results_folder=results_folder,env=TRUE,fig=TRUE,delete=FALSE,memorylimit=memory.limit())
We need to define the spatial domain of the model.
We're going to use the Atlantic part of the Canadian Exclusive Economic Zone from the Maritime Boundaries Geodatabase
EEZ <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="eez_iho_union_v2") plot(EEZ, col='blue')
Let's Remove Canadian part of the Davis Strait for EEZ since it's beyond cod's normal habitat.
plot(EEZ, col='blue') EEZ <- EEZ[EEZ$marregion!="Canadian part of the Davis Strait",] plot(EEZ, col='red',add=T)
The BESTMPA package is able to exclude certain habitats, or include them by default or as a priority (i.e. that habitat will be selected before proceeding to others). We're loading cod habitat and breeding grounds generated by georeferencing figure 2 from Lough (2004).
Habitats <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="cod_habitat") plot(Habitats,col='green') Breeding <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="cod_breeding")[2:15,] # remove first Breeding zone, outside EEZ plot(Breeding,col='yellow',add=T)
Define cell_size
in m, minimum size of MPA's in number of cells
, default projection and which areas to include in the grid's dataframe
# cell size in m cell_size <- 20000 # minimum size of MPA cells=1 # default projection proj <- "+proj=utm +zone=20 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0" EEZ <- spTransform(EEZ,proj) # create grid p <- initgrid(EEZ,cell_size,proj,areas=c(Habitats=Habitats,Breeding=Breeding)) plot(p) plot(p[p$Habitats==1,],col='green',add=T) plot(p[p$Breeding==1,],col='yellow',add=T)
Define basic parameters for protection scenarios:
# load existing MPAS oldMPA <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="MPAs_mar") oldMPA <- spTransform(oldMPA,CRS(proj)) #remove those smaller than cell size oldMPA <- oldMPA[gArea(oldMPA,byid = TRUE)>=cell_size^2,] oldMPA <- apply(gCovers(oldMPA,p,byid = TRUE),1,any)|apply(gOverlaps(oldMPA,p,byid = TRUE),1,any) # number of replicates replicates <- 1:100 # target protection level in proportion (e.g. 0.2 is 20% protection) MPA_coverage <- 0.10 # fixed distance for setting MPA distance in km in fixed distance scenario fixdist <- 75000
This is how the four scenarios were created. This portion of code is a little time consuming; therefore, I write the grid to file. If using the index.Rmd and want to generate new scenarios, remove the eval=FALSE
option from the chunk below
#### set status-quo scenario #### p <- addscenario(domain=p,included=oldMPA,replicates=replicates,name="MPA_SQ",cells=cells) #### set maximum distance #### p <- addscenario(domain=p,included=oldMPA,MPA_coverage=MPA_coverage,replicates=replicates,name="MPA_MD",cell_size=cell_size,cells=cells) #### set fixed distance #### p <- addscenario(domain=p,included=oldMPA,MPA_coverage=MPA_coverage,replicates=replicates,dist=fixdist,name="MPA_FX",cell_size=cell_size,cells=cells) #### set targeted scenario #### p <- addscenario(domain=p,included=oldMPA,priority=p$Habitats,excluded=p$Breeding,MPA_coverage=MPA_coverage,replicates=replicates,dist=fixdist,name="MPA_TR",cell_size=cell_size,cells=cells) writeOGR(p,dsn="shapefiles",layer="master_polygon",driver="ESRI Shapefile",overwrite_layer = TRUE)
p <- readOGR(dsn="shapefiles",layer="master_polygon")
Or you can plot all of the scenarios statically
# protect_scen_colour <- c("purple","green","red","blue") protect_scen_colour <- c("#D01C8B","#F1B6DA","#B8E186","#4DAC26") protect_scen_names <- c("Status Quo", "Maximum Distance", "Fixed Distance" , "Targeted") res <- 400 qual <- 100 textsize <- 10 # list scenarios for looping scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_" scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))] if(!file.exists(paste0(results_folder,"/figures"))) dir.create(paste0(results_folder,"/figures")) for(r in replicates){ if(r==1){ jpeg(paste0(results_folder,'/figures/f1_scenarios.jpg'),height=20,width=17,units="cm",res=res,qual=qual) } else { jpeg(paste0(results_folder,'/figures/fA',r-1,'_scenarios.jpg'),height=20,width=17,units="cm",res=res,qual=qual) } layout(matrix(c(1,2,3,4),nrow=2)) par(mar=c(1,1,1,1)) nam <- names(p)[grep(paste0("MPA_.._",r,"$"),names(p))] for(s in nam){ plot(EEZ) plot(p[p[[s]]==1,],col=protect_scen_colour[nam==s],add=T,border='transparent') title(protect_scen_names[which(nam==s)]) } dev.off() }
r=1 layout(matrix(c(1,2,3,4),nrow=2)) par(mar=c(1,1,1,1)) nam <- names(p)[grep(paste0("MPA_.._",r,"$"),names(p))] for(s in nam){ plot(EEZ) plot(p[p[[s]]==1,],col=protect_scen_colour[nam==s],add=T,border='transparent') title(protect_scen_names[which(nam==s)]) }
# total model run time in years (e.g. 2001:2100 would be 100 years) time <- 2021:2071 spinup <- 20 # number of years before "time" the model starts, results from spin-up years are not saved, all scenarios start as status quo tot_time <- (min(time)-spinup):max(time)
# Von Bertalanffy growth model parameters (Knickle and Rose 2013) # Lt = Linf * {1-exp(-k*(t-t0))}, where Lt is length (cm) at age t (year), Linf is the asymptotic length (cm), k is the VB growth coefficient (1/year), and t0 is the x intercept (year). Linf = 112.03 (95% CI = 10.46). k = 0.13 (95% CI = 0.021). t0 = 0.18). Linf_mean <- 112.03 Linf_SD <- 10.46/1.96 k_mean <- 0.13 k_SD <- 0.021/1.96 t0 <- 0.18 # calculate new weight - Length-weight relationship (Knickle and Rose 2013) fish$weight <- l_to_w_int * fish$length^l_to_w_power l_to_w_int <- 0.000011 l_to_w_power <- 2.91 # minimum age at maturity logistic equation age_mat_steepness <- 2.5 age_mat_sigmoid <- 4 # Fecundity (size dependent). 0.5 million eggs per kg of female fecundity <- 0.5*10^6 # Maximum age considered in matrix maxage <- 20 # calculate which breeding area is closest to each grid cell breeding_near <- apply(gDistance(spTransform(Breeding,proj),p,byid = T),1,which.min) # initial abundance initial_abun <- 250*10^6
# natural mortality (Swain & Chouinard 2008) M <- rnorm(10000,0.5938,0.0517) M <- M[M<=1&M>=0] # eliminate any possible values of M >1 or M <0 # larval mortality (Mountain et al. 2008) lM <- rbeta(10000,1000,1.2) #larval mortality of 99.88% (range 98.98-99.99%) #hist(lM);mean(lM);min(lM);max(lM) # Beverton-Holt model for carrying capacity based recruitment mortality, carrying capacity is the mean of North American carrying capacities in Table 3 of Myers et al. 2001 (mean of log CC=-1.202222222 tonnes/km^2 SD=0.9199018667) # Habitat carrying capacity, in kg of virtual fish per cell (4779 is the number of cells in the default grid). This could be substituted with "known" habitat carrying capacity. CCs <- (10^rnorm(10000,-1.202222222,0.9199018667))*(cell_size^2)/1000 CCs <- CCs[CCs>0][1:4779] #enforce no negative CCs # adjust CCs to be 0 outside habitat CCs <- CCs*p$Habitats
# larval dispersal kernels are assumed to be exponential, e_fold_larvae is the e folding scale in km (the distance at which there will be fewer settlers by a factor of e). We assume that scale to be sqrt of 2cm/s*90d (avg current velocity * PLD) because we assume that the current is like a random walk e_fold_larvae <- 2/100000*60*60*24*90*1000 # adult dispersal kernels are also assumed to be exponential, e_fold_adult (in km) was calculated from data in Lawson & Rose 2000 e_fold_adult <- sqrt(74.139*1000) # minimium age for adult migration (minimum size is 50 cm) Lawson & Rose 2000 min_age_migration <- 6 # generate connectivity matrix #adult cma <- initcm(p,e_fold_adult,cell_size) #larvae cml <- initcm(p,e_fold_larvae,cell_size)
Specify spatial distribution of fish_licenses
for shore distance calculation in effort calculation
# can be spatial points or spatial polygons data frame (sp package) fish_communities <- getData('GADM', country="CAN", level=1) # provinces fish_communities <- fish_communities[fish_communities$NAME_1=="New Brunswick"| fish_communities$NAME_1=="Newfoundland and Labrador"| fish_communities$NAME_1=="Nova Scotia"| fish_communities$NAME_1=="Prince Edward Island"| fish_communities$ID_1==11,] # this means Quebec, the accent make bad things happen when formatting # fish_communities <- getBigPolys(gSimplify(fish_communities,tol=0.05,topologyPreserve = T)) fish_communities <-gSimplify(fish_communities,tol=0.05,topologyPreserve = T) plot(fish_communities) # number of licenses per region in fish_communities fish_licenses <- c(866, 4714, 3002, 879, 963) # fisheries mortality at Maximum Sustainable Yield (Mountain et al. 2008) FMSY <- 0.28 # quota set to fraction of FMSY as per precautionary principle FMSY_buffer <- 2/3 # number of years to use in biomass estimate biomass_est_n_years <- 5 # minimium age caught by nets (minimum size is 38 cm) Feekings et al. 2013 min_age_catch <- 4 # calculate distances from shore distance <- gDistance(spTransform(fish_communities,proj),p,byid = T) # list scenarios for looping scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_" scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))]
If using the index.Rmd and want to generate new data, remove the eval=FALSE
option from the chunk below
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish <- growth(fish,y,tot_time,initial_abun,p,maxage,recruits,M) #### reproduction #### recruits <- reproduction(fish,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish <- dispersal(fish,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits <- larvaldispersal(recruits,cml,lM,CCs,fish) ############################ Harvesting ################################################# quota <- estimatequota(fish,maxage,y,tot_time,FMSY,FMSY_buffer) fish <- harvest(fish,quota,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder,y=y,time=time) } } close(pb)
########################################### cost evaluation ############################################## # normal operating cost ($) per fisherman from (Department of Fisheries and Oceans, 2007, Table A.19 Mixed Fishery Fleet) # labour (CAD $46587) and fuel ($9008) vary with distance, the remainder does not fish_operating_cost_ratio <- (46587+9008)/105054 # (labour+fuel)/total # default profitability to calibrate operating cost (Department of Fisheries and Oceans, 2007, Table A.19 and 5.7 Mixed Fishery Fleet) Status_quo_profitability <- 166184/105054 #$ catch value/$ operating expense # landed value for cod CAD/kg (http://www.dfo-mpo.gc.ca/stats/commercial/sea-maritimes-eng.htm) fish_landed_value <- 1.24 # Social discount rates (choose 3) SDR <- c(0.015,0.03,0.06) textsize <- 10
# load fish catches filenames <- list.files(results_folder,pattern="kg_") catch_summary <- rbindlist(lapply(filenames,processkg,results_folder,distance,fish_landed_value)) catch_summary <- calcnetvalue(catch_summary,Status_quo_profitability,fish_operating_cost_ratio) catch_summary <- SDRvalues(catch_summary,t=time,value=catch_summary$netvalue,SDR) catch_summary$scenario <- ordered(catch_summary$scenario,levels=c("SQ","MD","FX","TR")) filenames <- list.files(results_folder,pattern="fish_") fish_summary <- rbindlist(lapply(filenames,processfish,results_folder)) fish_summary$scenario <- ordered(fish_summary$scenario,levels=c("SQ","MD","FX","TR")) levels(catch_summary$scenario) <- protect_scen_names levels(fish_summary$scenario) <- protect_scen_names gathered <- catch_summary %>% dplyr::select(scenario,replicate,year,starts_with("SDRcumsum")) %>% gather(SDR,cumsumvalue,starts_with("SDRcumsum"))
table1 <- catch_summary %>% group_by(scenario) %>% summarize(Catch=paste0(round(mean(kg)/1000)," (\u00b1 ",round(std.error(kg)/1000),")"), Distance=paste0(round(mean(distanceGV)/1000)," (\u00b1 ",round(std.error(distanceGV)/1000),")")) table1 <- cbind(table1,fish_summary %>% group_by(scenario) %>% summarize(Biomass=paste0(round(mean(biomass)/1000)," (\u00b1 ",round(std.error(biomass)/1000),")")) %>% select(Biomass)) table1 <- cbind(table1,catch_summary %>% group_by(scenario,replicate) %>% summarize(Mort=sum(kg<2000000)/n()*100) %>% group_by(scenario) %>% summarize(Moratorium=paste0(round(mean(Mort),2)," (\u00b1 ",round(std.error(Mort)),")")) %>% select(Moratorium)) table1[,c(1,4,2,3,5)]
########## Figure 3 - plot of biomass over time ################################# p1 <- ggplot(fish_summary,aes(x=year,y=biomass/1000,colour=scenario))+ geom_smooth(cex=2)+ theme_classic()+ theme(legend.position="top", text=element_text(size=textsize))+ labs(x="Year",y="Total Stock Biomass (t)")+ scale_colour_manual(values=protect_scen_colour, labels=protect_scen_names,name="") ggsave(paste0(results_folder,'/figures/f3_biomass.jpg'),p1,height=8,width=17,units="cm",dpi=res) p1 ########## Figure 4 - plot of catch biomass over time ################################# p1 <- ggplot(catch_summary,aes(x=year,y=kg/1000,colour=scenario))+ geom_smooth(cex=2)+ theme_classic()+ theme(legend.position="top", text=element_text(size=textsize))+ labs(x="Year",y="Total catch (t)")+ scale_colour_manual(values=protect_scen_colour, labels=protect_scen_names,name="") ggsave(paste0(results_folder,'/figures/f4_catch.jpg'),p1,height=8,width=17,units="cm",dpi=res) p1 ########## Figure 5 - distance from shore over time ################################# p1 <- ggplot(catch_summary,aes(x=year,y=distanceGV/1000,colour=scenario))+ geom_smooth(cex=2)+ theme_classic()+ theme(legend.position="top", text=element_text(size=textsize))+ labs(x="Year",y="Mean distance from shore (km)")+ scale_colour_manual(values=protect_scen_colour, labels=protect_scen_names,name="") ggsave(paste0(results_folder,'/figures/f5_distance.jpg'),p1,height=8,width=17,units="cm",dpi=res) p1 ########## Figure 6 - Social Discount Rate ################################# jpeg(paste0(results_folder,'/figures/f6_SDR.jpg'),height=20,width=17,units="cm",res=res,qual=qual) p1 <- ggplot(gathered,aes(x=year,y=cumsumvalue/10^6,colour=scenario,xlab="test"))+ geom_smooth(cex=1.5)+ theme_classic()+ scale_colour_manual(values=protect_scen_colour,labels = protect_scen_names,name="")+ facet_wrap(~SDR,ncol=1,scale="free_x")+ theme(strip.background=element_blank(), legend.position="none", strip.text=element_text(face="bold"), text=element_text(size=textsize))+ labs(x="Year",y=expression(paste("Cumulative Net Present Value (10"^"6"," CAD)",sep=""))) p2 <- ggplot(gathered[gathered$year==max(gathered$year),],aes(x=scenario,y=cumsumvalue/10^6,fill=scenario))+ # geom_jitter(size=0.5)+ # geom_violin(alpha=0.7)+ geom_boxplot()+ theme_classic()+ scale_fill_manual(values=protect_scen_colour,labels = protect_scen_names,name="")+ facet_wrap(~SDR,ncol=1,scale="free_x")+ theme(strip.background=element_blank(), legend.position="none", strip.text=element_text(face="bold"), axis.text.x=element_text(angle=20,hjust=1), text=element_text(size=textsize))+ labs(x="",y=expression(paste("Net Present Value (10"^"6"," CAD)",sep=""))) g <- ggplotGrob(p1 + theme(legend.position="bottom"))$grobs legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]] lheight <- sum(legend$height) grid.arrange(legend,arrangeGrob(p1,p2,ncol=2),heights=unit.c(lheight, unit(1, "npc") - lheight)) dev.off()
knitr::include_graphics(paste0(results_folder,'/figures/f6_SDR.jpg'), dpi=400)
results_folder_sens <- "BESTMPA_sensitivity" housekeeping(results_folder=results_folder_sens,env=FALSE,fig=TRUE,delete=FALSE,memorylimit=memory.limit()) sens_factor <- 0.1 # list scenarios for looping scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_" #select only SQ and TR scenarios scenarios <- scenarios[c(starts_with("MPA_SQ",vars=scenarios),starts_with("MPA_TR",vars=scenarios))] scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))] # re-order numerically # select only first 1 scenarios of each scenarios <- scenarios[1:50]
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity*(1-sens_factor),age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity*(1+sens_factor),age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="fecundity_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="fecundity_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="initial_abun_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="initial_abun_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid*(1-sens_factor),l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid*(1+sens_factor),l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="age_mat_sigmoid_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="age_mat_sigmoid_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean*(1-sens_factor),k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean*(1+sens_factor),k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="k_mean_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="k_mean_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=(min_age_migration-1):maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=(min_age_migration+1):maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_migration_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_migration_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=(min_age_catch-1):maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_catch_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=(min_age_catch+1):maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_catch_plus") } } close(pb)
# generate connectivity matrix #adult cma_minus <- initcm(p,e_fold_adult*(1-sens_factor),cell_size) cma_plus <- initcm(p,e_fold_adult*(1+sens_factor),cell_size) pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma_minus,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma_plus,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_adult_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_adult_plus") } } close(pb)
# generate connectivity matrix #larvae cml_minus <- initcm(p,e_fold_larvae*(1-sens_factor),cell_size) cml_plus <- initcm(p,e_fold_larvae*(1+sens_factor),cell_size) pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml_minus,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml_plus,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_larvae_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_larvae_plus") } } close(pb)
pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0) for(s in scenarios){ for(y in tot_time){ gc() # keeps memory from getting clogged setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y)) ############################ Growth and reproduction ################################################# fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M) fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M) #### reproduction #### recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p) ############################ Dispersal ################################################# #### adult dispersal #### fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage) fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage) #### larval dispersal #### recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus) recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus) ############################ Harvesting ################################################# quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY*(1-sens_factor),FMSY_buffer) quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY*(1+sens_factor),FMSY_buffer) fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="FMSY_minus") fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="FMSY_plus") } } close(pb)
# load fish catches filenames <- list.files(results_folder_sens,pattern="kg_") catch_summary_sens <- rbindlist(lapply(filenames,processkg,results_folder_sens,distance,fish_landed_value)) catch_summary_sens <- calcnetvalue(catch_summary_sens,Status_quo_profitability,fish_operating_cost_ratio) catch_summary_sens <- SDRvalues(catch_summary_sens,t=time,value=catch_summary_sens$netvalue,SDR) filenames <- list.files(results_folder_sens,pattern="fish_") fish_summary_sens <- rbindlist(lapply(filenames,processfish,results_folder_sens)) gathered <- catch_summary_sens %>% dplyr::select(scenario,replicate,year,starts_with("SDRcumsum")) %>% gather(SDR,cumsumvalue,starts_with("SDRcumsum")) radialmeans <- catch_summary_sens %>% filter(year==max(tot_time)) %>% group_by(scenario) %>% summarize(meanSDR=mean(`SDRcumsum = 0.015`)) %>% mutate(variable=gsub("^[^_]+_|_[^_]+$", "", scenario), scenario=gsub("_.*_", "_", scenario)) %>% spread(key=scenario,value=meanSDR) %>% ungroup() %>% as.data.frame() catch_summary_means <- catch_summary %>% filter(year==max(tot_time)) %>% group_by(scenario) %>% summarize(meancumsum = mean(`SDRcumsum = 0.015`)) %>% spread(scenario,meancumsum) catch_summary_means <- catch_summary_means[rep(1,nrow(radialmeans)),]
jpeg(paste0(results_folder,'/figures/f7_sensitivity.jpg'),height=17,width=17,units="cm",res=res,qual=qual) radial.plot(t(cbind(catch_summary_means$`Status Quo`,catch_summary_means$Targeted)), labels=as.character(radialmeans[,1]), rp.type="p", radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2), poly.col="transparent", line.col=protect_scen_colour[c(1,4)], show.grid.labels=1, lwd=2, lty=1) radial.plot(t(radialmeans[,2:5]), labels=as.character(radialmeans[,1]), rp.type="s", radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2), point.symbols=c("-","+","-","+"), point.col=protect_scen_colour[c(1,1,4,4)], show.grid.labels=1, lwd=2, cex=2.5, add=T) legend(max(radialmeans[,2:5])*0.7, max(radialmeans[,2:5])*2.2, c("Status Quo","Status Quo 90%","Status Quo 110%","Targeted","Targeted 90%","Targeted 110%"), col=protect_scen_colour[c(1,1,1,4,4,4)], lty=c(1,NA,NA,1,NA,NA), pch=c(NA,"-","+",NA,"-","+"), lwd=2, pt.cex=2.5, bty="n") dev.off()
########## Figure 6 - Social Discount Rate ################################# radial.plot(t(cbind(catch_summary_means$`Status Quo`,catch_summary_means$Targeted)), labels=as.character(radialmeans[,1]), rp.type="p", radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2), poly.col="transparent", line.col=protect_scen_colour[c(1,4)], show.grid.labels=1, lwd=2, lty=1) radial.plot(t(radialmeans[,2:5]), labels=as.character(radialmeans[,1]), rp.type="s", radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2), point.symbols=c("-","+","-","+"), point.col=protect_scen_colour[c(1,1,4,4)], show.grid.labels=1, lwd=2, cex=2.5, add=T) legend(max(radialmeans[,2:5])*0.7, max(radialmeans[,2:5])*2.2, c("Status Quo","Status Quo 90%","Status Quo 110%","Targeted","Targeted 90%","Targeted 110%"), col=protect_scen_colour[c(1,1,1,4,4,4)], lty=c(1,NA,NA,1,NA,NA), pch=c(NA,"-","+",NA,"-","+"), lwd=2, pt.cex=2.5, bty="n")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.