#-----------------------------------------# # GRIP 2.0 # #-----------------------------------------# # A program to generate sets of # random input parameters (GRIP) # for sensitivity analysis of # population viability analysis # (PVA) models created with # RAMAS GIS 5.0 and RAMAS GIS 6.0. # GRIP 2.0 applies a # global sensitivity analysis and # varies spatial and non-spatial # parameters and modifies # habitat-based parameters on a # raster reference landscape. #-----------------------------------------# # Project Aims # #-----------------------------------------# # GRIP 2.0 was developed to # stimulate interest in # incorporating options for # sensitivity analysis of spatial # parameters into PVA software # packages, and to provide # interim tools for use in species # at risk recovery planning and # critical habitat assessments. #-----------------------------------------# # Code modified: February 2016 # # # # Authors: # # Ilona Naujokaitis-Lewis and # # Janelle M. R. Curtis # # # #-----------------------------------------# # Contact Information: # Ilona Naujokaitis-lewis # Email: ilona.naujo.lewis@gmail.com # Mail: National Wildlife Research Centre, # Environment and Climate Change Canada, # Carleton University, # 1125 Colonel By Drive, # Ottawa, Ontario, Canada K1A 5B6 # Janelle M. R. Curtis # Email: janelle.curtis@dfo-mpo.gc.ca # Mail: Pacific Biological Station # Fisheries and Oceans Canada # 3190 Hammond Bay Road, Nanaimo # British Columbia, Canada V9T 6N7 #----------------------------------# # Notes to Users # #----------------------------------# # Please report bugs, # and contact us with questions, # comments, or suggestions to # improve the code. # We would like to track use # of GRIP 2.0. # Please let us know if you are # planning to use or # customize this code. # or better yet: please cite! #-----------------------------------------# # Recommended citation: TBA. #-----------------------------------------# #-----------------------------------------# #-----------------------------------------# #-----------------------------------------# # INSTRUCTIONS AND NOTES # #-----------------------------------------# #-----Overview of GRIP 2.0----------------# # Grip 2.0 builds on GRIP 1.0 # (Curtis & Naujokaitis-Lewis 2008) # and is designed to vary parameter # values of PVA models created with # RAMAS GIS v5.0 and v6.0. # For more information about RAMAS # software go to www.ramas.com # GRIP 2.0 was designed to run as an # 'R' script and was tested on # versions R 2.7.0 to R 3.23.0. # Should there be incompatibilities with # different versions of R, we # suggest that users 1) use 1 of # the versions that it was tested on, # or, 2) let the authors know of any # problems and we will do our best # to respond in a timely manner and assist in any way. #-----------------------------------------# #---GRIP 2.0: Major steps-----------------# #-----------------------------------------# # GRIP 2.0: generates random sets of # input parameters in the following major steps: # 1. The program READS IN the following files # (A) RAMAS Spatial 5.0/6.0 PVA model (*.ptc file), # (B) RAMAS Metapop file (*.mp), # (C) a habitat suitability map (*.asc), # (D) and a patch ID map (*.asc). # The first two files ((A): *.ptc and (B) *.mp) # require certain configurations. # See section 'Prior to running GRIP 2.0'. # 2. Spatial parameters are varied on # the habitat suitability (raster) map # and a modified version of RAMAS Spatial file # is created. The two outputs at this stage include: # (1) 'ptc_*.ptc' (generic file name), and # (2) 'HSmap_*.asc' # One file is created for each of the '*' replicate simulations. # Parameters varied during this stage include: # spatial parameters: number of populations #(we use the terms population and patches synonymously), # habitat suitability values, habitat suitability # threshold, neighbourhood distance, # method for calculating distances between patches. # The majority of these parameters are varied # directly on raster-based habitat suitability map. # A new habitat suitability map is created # (based on the varied parameters), which is # used to define the new metapopulation structure. # Each replicate RAMAS Spatial Data file (ptc_*.ptc) # is run in RAMAS Spatial 6.0 using a # batch file created by GRIP 2.0. # Results from runs of each replicate file are # recorded, saved, and collated into two tables, # "PTCvariables.txt" and "spatialvar.txt". # Results include: total amount of habitat, # mean patch size, average patch quality, # number of populations, mean core area, # mean edge to area ratio, mean patch edge, # mean shape index, and percent area habitat # within the landscape. # The GRIP 2.0 code can be customized by the # user to collate additional variables of interest. # 3. A 'pref_metapop' file is created, # which is created from the output # from the previously run ptc_*.ptc file # (called pref_*.mp, one is created for each # of the n replicate simulations). # This is created by GRIP 2.0 because # a copy of the metapop file is not created # automatically when RAMAS Spatial models # are run using batch files. Because GRIP 2.0 # uses output from the RAMAS Spatial Data # program and the user-specified Metapop Link # file, users must ensure that parameters # specified in the original Metapop Link # file (*.mp) are specified the same way in the # RAMAS GIS model (*.ptc). # See more details on ensuring consistency # between the two files below. # 4. A 'prefile' file is created, which # is an interim replicate simulation file # (called interim_*.mp, one is created # for each of the n replicate simulations). # Key steps during the 'prefile' creation # stage of GRIP 2.0 include varying # parameters such as initial abundances and Rmax. # In addition, a population with an arbitrarily # low K is added to simulate dispersal survival # (see Curtis and Naujokaitis-Lewis 2008). # 5. A 'repfile' file is created, which is # a final replicate simulation file to be run # during the RAMAS Metapop simulations # (called rep_*.mp, one is created for # each of the n replicate simulations). # Key steps at this stage of GRIP 2.0 include # varying parameters including: frequency, # intensity and spatial extent of catastrophes, # dispersal rates, number of population pairs # connected by dispersal, correlations in # vital rates, and vital rates. # See Curtis and Naujokaitis-Lewis 2008 for details. # 6. All rep_*.mp files are run in # RAMAS Metapop using a batch file # created by GRIP 2.0. # Results from each of the n repfiles # are recorded, saved, and collated into one table # (data.csv) for further analysis. #-----------------------------------------# # GRIP 2.0 CODE STARTS HERE # #-----------------------------------------# # Install/load libraries # If not already installed, use install.packages("xxx") library(sp) library(rgdal) library(spatial) library(spatstat) library(raster) library(plotrix) library(adehabitatHS) library(ade4) library(SDMTools) library(maptools) library(RandomFields) # Removes all objects from the R workspace rm(list=ls(all=TRUE)) #-----------------------------------------# # The following lines of code # # require USER-specification # #-----------------------------------------# # Name of original RAMAS Spatial file # Model was run, and includes saved results # Needed to reference certain parameters/values Original_Spatial <- "WhPine.ptc" #e.g. "Whitebarkpine.ptc" # Reference habitat suitability map # File exported from results of # RAMAS GIS output model # File could be of any type allowed by RAMAS GIS # If ensemble approach to estimate habitat suitability # was applied, this file contains mean/median/etc predictions Original_HSmap <- "WhPine_HS.ASC" #e.g. "Whitebarkpine_HS.ASC" # Reference patch map file exported # from results of RAMAS GIS output model Original_Patchmap <- "WhPine_PA.ASC" #e.g."Whitebarkpine_PA.ASC" # Reference mask map (ascii) file, if it is in use. # The mask map is used if there are areas within the # landscape where patches cannot be created. # For e.g., for a terrestrial species patches #can only be created on land, so a mask map # might include any water bodies included # within the modeled landscape extent. # If a mask is used, remove the # sign # from the first mask object. # If a mask map is not used, use the line of code # specifying 'mask<-0', i.e. remove '#' sign # from this line of code. # Origmask<-"Speciesname_mask.ASC" #e.g "Whitebarkpine_mask.ASC" Origmask <- 0 # Specify the number of replicate simulations Nreps <- 5 # Specify the number of stochastic runs # within each simulation (default is 1000) Nruns <- 1000 #---------------------------------------------# # Specify method to vary # # habitat suitability values # #---------------------------------------------# # Users can select the method by which to vary HS values # These include: # "SACexp": generate a landscape using randomfield package # Users can specify the model and degree to which # Spatial autocre # "rnorm": each new HS sampled randomly based on gaussian distribution # If rnorm selected,New habitat suitability values are sampled from normal distribution with mean = mean of HS values in original landscape, SD = standard deviaion of HS values. If the value is < the HS threshold, it is replaced with the threshold value. # "ensemble": uses ensemble predictions and the # measure of uncertainty to vary new HS values # the current implementation of GRIP2 resamples # HS values using a random normal variate with # the mean based on the ensemble prediction # with a SD based on the uncertainty estimate from the # ensemble model (i.e., the model-based measure of variation) # HS resampling is done on a cell-by-cell basis # users can specify alternative distributions by # which to sample from, but 'rnorm' is the current implementation # 'none': users can opt to not vary HS # Select which method to vary habitat suitability # values by choosing one of the following: # Ensure that the line with the method # for varying HS does not have a '#' # in front of it (i.e., keep the '#' # in front of all lines of code you # do NOT want to be run in R #HSvary <- "SACexp" #HSvary <- "rnorm" HSvary <- "ensemble" #HSvary <- "none" # If ensemble methods for used for HS sampling/variation # are applied, specify the name of the input raster file # containing the measure of variation if (HSvary == "Yes") sdHSmap <- "ensemble_sd.ASC" # input name of SD map #---END OF USER-SPECIFIED CODE---------------# #--------------------------------------------# #--------------------------------------------# # Create batch file to remove all previous # # files from working directory. Includes: # # prefiles, repfiles, saved results files # # created by GRIP 2.0 # #--------------------------------------------# write("",file="clearspace.bat", append = FALSE) write("del HSHist*", file = "clearspace.bat", append = TRUE) write("del PatchSummary*", file = "clearspace.bat", append = TRUE) write("del LandInd*", file = "clearspace.bat", append = TRUE) write("del ptc*", file = "clearspace.bat", append = TRUE) write("del HSmap*", file = "clearspace.bat", append = TRUE) write("del pref*", file = "clearspace.bat", append = TRUE) write("del Patch.RES", file = "clearspace.bat", append = TRUE) write("del patchBAT.REC", file = "clearspace.bat", append = TRUE) write("del batch_ptc.bat", file = "clearspace.bat", append = TRUE) write("del PTCvariables.txt", file = "clearspace.bat", append = TRUE) write("del spatialvar.txt", file = "clearspace.bat", append = TRUE) write("del LandIndices*", file = "clearspace.bat", append = TRUE) write("del PatchSummary*", file = "clearspace.bat", append = TRUE) write("del Abund*", file = "clearspace.bat", append = TRUE) write("del FinalStageN*", file = "clearspace.bat", append = TRUE) write("del Harvest*", file = "clearspace.bat", append = TRUE) write("del IntExpRisk*", file = "clearspace.bat", append = TRUE) write("del IntExtRisk*", file = "clearspace.bat", append = TRUE) write("del IntPerDec*", file = "clearspace.bat", append = TRUE) write("del LocalOcc*", file = "clearspace.bat", append = TRUE) write("del LocExtDur*", file = "clearspace.bat", append = TRUE) write("del MetapopOcc*", file = "clearspace.bat", append = TRUE) write("del QuasiExp*", file = "clearspace.bat", append = TRUE) write("del QuasiExt*", file = "clearspace.bat", append = TRUE) write("del TerExpRisk*", file = "clearspace.bat", append = TRUE) write("del TerExtRisk*", file = "clearspace.bat", append = TRUE) write("del TerPerDec*", file = "clearspace.bat", append = TRUE) write("del interim*", file = "clearspace.bat", append = TRUE) write("del rep*", file = "clearspace.bat", append = TRUE) write("del metapop.res", file = "clearspace.bat", append = TRUE) write("del metabat.rec", file = "clearspace.bat", append = TRUE) write("del batch_file.bat", file = "clearspace.bat", append = TRUE) write("del vital_rates.txt", file = "clearspace.bat", append = TRUE) write("del predictions.txt", file = "clearspace.bat", append = TRUE) write("del variable*", file = "clearspace.bat", append = TRUE) write("del data.csv", file = "clearspace.bat", append = TRUE) # Runs the clearspace.bat file system("clearspace.bat") # Note: you can also clear the working directory manually by simply double-clicking on clearspace.bat in the working directory after sourcing GRIP. #---END batch clear working directory--------# #--------------------------------------------# #---Track parameters varied------------------# # Writes column names of the varied input # parameters to "PTCvariables.txt" PTCvarName <- list("Npopsamp","HSthreshold","NghbdDistance","HSvary","PatchChange","RadiusAddPatch") write.table(PTCvarName, file="PTCvariables.txt", col.names = FALSE, row.names= FALSE, sep=",", quote=FALSE, append=FALSE) #---Read in batch file-----------------------# # Reads the batch_GRIP2.txt file required for # creating the batch file for running input files # in RAMAS Spatial Module 5.0 Start <- readLines("batch_GRIP2_v6.txt") #--------------------------------------------# # in final version, TO DELETE #---Functions--------------------------------# # Mirror matrix (left-right) mirror.matrix <- function(x) { xx <- as.data.frame(x); xx <- rev(xx); xx <- as.matrix(xx); xx; } # Rotate matrix 90 clockworks rotate90.matrix <- function(x) { t(mirror.matrix(x)) } #---Dilate function-----# # Applied to 'owin' object types # Needed for adding new patches to the landscape # 'dilate2.owin' function: ensures dimensions # of new window same as original landscape # for when a patch is added is of the dilate2.owin <- function (w, r, ..., include.bb=FALSE) { verifyclass(w, "owin") if (r < 0) stop("r must be nonnegative") if (r == 0) return(w) if (w$type != "mask") w <- as.mask(w, ...) epsilon <- sqrt(w$xstep^2 + w$ystep^2) r <- max(r, epsilon) if (!include.bb) bb <- bounding.box(w) else bb <- w newbox <- owin(bb$xrange, bb$yrange) w <- rebound.owin(w, newbox) distant <- distmap(w) dil <- w dil$m <- (distant$v <= r) return(dil) } #-----------------------------------------------# # START: Generate alternative landscapes # #-----------------------------------------------# for (y in 1:Nreps){ ptcfile <- paste("ptc", "_",y,".ptc", sep="") new_HSfunction <- paste("[HSmap","_",y,"]",sep="") newHSmap_name <- paste("HSmap","_",y,sep="") newHSmap_asc <- paste("HSmap","_",y,".ASC", sep="") # Creates a batch file # to run replicate simulation files in RAMAS batch_rep <- paste(Start[1], " ",Start[2], " ", Start[3], " ", Start[4],"ptc_", y, ".ptc", Start[5], sep="") write(batch_rep, file="batch_ptc.bat", append = TRUE) #LandIndices <- paste("rename LandIndices.txt ", "LandIndices_", y, ".txt", sep="") #PatchSummary <- paste("rename PatchSummary.txt ", "PatchSummary_", y, ".txt", sep="") # write(LandIndices, "batch_ptc.bat", append = TRUE) # write(PatchSummary, "batch_ptc.bat", append = TRUE) #---START preliminary set-up---------------# # Reads in the original input Habitat Suitability (HS) map origHSmap <- readGDAL(Original_HSmap) # Reads in the original input Patch map origPatchmap <- readGDAL(Original_Patchmap) # Reads in the mask map, if it exists if (Origmask==0){ Origmask<-0 } else { maskmap<-readGDAL(Origmask) } # Changes names of data in SpatialDataFrames names(origHSmap) = paste("HSvalue", sep="") names(origPatchmap) = paste("ID", sep="") # Retrieve grid parameters and assign object to grid cell size gridpara <- gridparameters(origPatchmap) cellsize<-gridpara[1,2] celldimX<-gridpara[1,3] celldimY<-gridpara[2,3] # Reads in the whole input file (*.ptc) as an object inputprefile <- readLines(Original_Spatial, n=-1) # Note: used for referencing lines with the GREP function and for writing lines that are not changed by GRIP. # Returns line number in the input file with "Results" Results <- grep("Results", inputprefile) # Returns line number in the input file with "Migration" Migration<-grep("Migration", inputprefile) # Returns line number in the input file with "End of file" Endoffile<-grep("End of file", inputprefile) # Reads in the original number of populations # in the PVA model as a numeric object Original_Npops <- scan(Original_Spatial, what = "list", skip = Results, nlines = 1) Original_Npops <- as.numeric(Original_Npops[1]) # Reads in number of decimals for habitat suitability map export HSdecimal <- as.numeric(scan(Original_Spatial,what="numeric",skip=12,nlines=1)) # Reads in original Habitat Suitability Threshold # for patches identification algorithm HSthreshold <- as.numeric(scan(Original_Spatial, what = "list", skip = 9, nlines = 1)) # Sample new habitat suitability threshold # applying the original threshold as a mean and a 10% CV new_HSthreshold <- round(rnorm(1,mean=HSthreshold,sd=0.1*HSthreshold),digits=HSdecimal) # Reads in the results population dataframe # from the '*.ptc' results section pop_data <- read.table(Original_Spatial, skip = (Results+Original_Npops+3), nrow = Original_Npops,sep="") # Mean area of patches of original landscape (# of cells) avg_area_pop<-mean(pop_data$V2) # Standard deviation of patch areas of # original landscape (# of cells) sd_area_pop<-sd(pop_data$V2) # Calculates the minimum (greater than 0) and maximum possible habitat suitability values for the # model based on the habitat suitability map # Rounds Patchmap values to 0 decimal places; these values are always whole numbers between 1:n, where n is the number of patches identified in the original model origPatchmap$ID <- round(origPatchmap$ID,0) # Rounds HSmap decimals places to those specified in the original (*.ptc) model origHSmap$HSvalue <- round(origHSmap$HSvalue,HSdecimal) # Coerces NA values to 0s origHSmap$HSvalue[is.na(origHSmap$HSvalue)]<-0 # Creates duplicate object for backup Patchmap <- origPatchmap HSmap <- origHSmap # Calculates the mean and SD habitat suitability across habitat # Takes the mean of cells that are # => habitat suitability threshold of the original model mean_HS <- mean( origHSmap$HSvalue[ origHSmap$HSvalue >= HSthreshold] ) #same as:avgHS_cell<-Tot_HS/Tot_cells sd_HS <- sd( origHSmap$HSvalue[ origHSmap$HSvalue >= HSthreshold ]) # set SD of HS to 10% CV if it is 0 in the original landscape map if (sd_HS == 0) sd_HS<-0.1 # Reads in original neighbourhood distance for patch identification algorithm (distance is equivalent to number of cells) nghbd_dist<-as.numeric(scan(Original_Spatial, what = "list", skip = 10, nlines = 1)) # Sample new neighbourhood distance # based on the original value and a CV of 10% line11_nghbd_dist<-round(rnorm(1,mean=nghbd_dist,sd=0.1*nghbd_dist),digits=2) # If newly sampled neighbourhood distance # is less than 1, coerces to 1 # RAMAS requires a minimum distance of 1.00 if (line11_nghbd_dist<1) line11_nghbd_dist<-1 # Reads in friction map file (should not be specified) fric_map<-scan(Original_Spatial,what="character",skip=25,nlines=1) # Vary method to measure distance between patches # Creates a vector of the 3 methods RAMAS uses distance<-c("Edge to edge","Center to edge","Default: Center to center") # Randomly samples distance method line27_distance_measure<-sample(distance,1) # Reads in number of input maps # This should include the habitat suitability map and maps used to define how input parameters are a function of habitat suitability nu_maps <- as.numeric(scan(Original_Spatial,what="numeric",skip=27,nlines=1)) # Sample the new number of populations (random normal) # with mean = original number of populations # with a 50% CV Npops <- abs(round(rnorm(1, mean=Original_Npops, sd=0.5*Original_Npops), digits=0)) # Difference between number of populations in the original versus the new landscape Diff <- Original_Npops - Npops # Ensures that Diff is not equal to the number of original populations, as this would translate to a landscape with no populations # if Npops sampled = 0, change so that Npops = 2 # as there must be a minimum of 2 for dispersal matrix # and error-free code functioning if (Diff == Original_Npops){ Diff <- Original_Npops-2 Npops <- 2 } # if Npops = 1, change so that Npops = 2 if (Diff == (Original_Npops-1)){ Diff <- Original_Npops-2 Npops <- 2 } #END preliminary set-up---------------# #-------------------------------------# # Part 1: # # Remove or add populations # # to landscapes # # or change names of objects, # # if patch number stays same # # as original # #-------------------------------------# #-------------------------------------# # REMOVE POPULATIONS # #-------------------------------------# if (Original_Npops > Npops){ # Randomly selects which populations to remove # and sorts them in descending order Pops_removed <- sample(seq(1, Original_Npops), Diff, replace=FALSE) Pops_removed <- sort(Pops_removed, decreasing = TRUE) # Coerce NAs in the spatial grid object to 0s HSmap$HSvalue[is.na(HSmap$HSvalue)] <- 0 # Add HSvalues to raster Patchmap$HSvalue <- HSmap$HSvalue # Create a vector containing the names of # the expressions to be evaluated Npopslength <- length(Pops_removed) pnames <- paste("(Patchmap$ID == Pops_removed[", 1:Npopslength, "])",sep="") # Combine them into a single string that # represents the expression to be evaluated: pstring <- paste(pnames, collapse=" | ") # Convert this into an R language expression: pexpr <- parse(text = pstring) # Execute eval peval <- do.call("eval", list(pexpr)) # ID2: patches that need to be retained Patchmap$ID2 <- ifelse(peval==TRUE, 0, Patchmap$ID) # Replace HS values of patches to be removed with 0 Patchmap$HSvalue2 <- ifelse(Patchmap$ID2==0, 0, Patchmap$HSvalue) # Convert to raster stack test <- stack(Patchmap) # Create new patchmap object new_patchmap <- test$ID2 names(new_patchmap) <- "ID" # Create new hsmap object new_HSmap <- test$HSvalue2 names(new_HSmap) <- "HSvalue" # Assign '0' to NA set0toNA <- function(x) { x[x==0] <- NA; return(x) } new_HSmap <- calc(new_HSmap, set0toNA) new_patchmap <- calc(new_patchmap, set0toNA) } #---END remove populations------------# #-------------------------------------# # ADD POPULATIONS # #-------------------------------------# if (Original_Npops < Npops) { # Ensure that object is gridded gridded(Patchmap) <- TRUE # Coerces 0s to NAs Patchmap$ID[Patchmap$ID==0]<-NA # Retrieve min and max of (x,y) values from # bounding box matrix of spatial data Patchmap_bbox <- bbox(Patchmap) min_x <- Patchmap_bbox[1,1] max_x <- Patchmap_bbox[1,2] min_y <- Patchmap_bbox[2,1] max_y <- Patchmap_bbox[2,2] # Creates a binary mask of dilated patches and matrix, # patches are dilated by 1 cell around all sides # Convert to a matrix mt <- as.matrix(Patchmap) # Rotate image for correct orientation mt_rot <- rotate90.matrix(mt) #dont use: mt_rot <- mirror.matrix(mt) # Convert to an ascii class object PaIDasc <- as.asc(mt_rot) #image(PaIDasc) #image(mt_rot) # Coerce 0s to NAs PaIDasc[PaIDasc==0] <- NA # Create new object PaIDasc_D <- PaIDasc #image(PaIDasc_D) # Ensure NAs are translated to value '0' PaIDasc_D[is.na(PaIDasc_D)] <- 0 #image(PaIDasc_D) # Coerce to logical matrix PaIDasc2_mt <- matrix(PaIDasc_D,nrow=nrow(PaIDasc_D),ncol=ncol(PaIDasc_D)) PaIDasc2mt_logical <- as.matrix(PaIDasc2_mt==0) # Create a mask layer mask2 <- owin(xrange=c(min_x,max_x), yrange=c(min_y,max_y), mask=PaIDasc2mt_logical) # Samples centre cell for x number of patches to add, samples on patch map layer, but only in area defined as non-habitat matrix # Generate a random point pattern, a realisation of the Simple Sequential Inhibition (SSI) process. # r=Inhibition distance; n=# of points to generate # Generic sampling function: rSSI(r, n, win = owin(c(0,1),c(0,1)), giveup = 1000) # If the new point lies closer than r units from an existing point, then it is rejected and another random point is generated # r is set at 5 times the original cell size to insure that new patch centres are not located too close to one another pts2 <- rSSI(5*cellsize, -Diff, as.mask(mask2), giveup = 1000) #plot(mask2) #plot(mask) #plot(pts2,add=TRUE, col = 'yellow') #---End sample new patch centroids if no landscape mask #-----------------------------------------# # Add new patches to the spatial grid # #-----------------------------------------# # Retrieves grid topology parameters ) Patchmap_topo <- getGridTopology(origPatchmap) # Produce a vector containing the unique IDs (i.e 1:length(origPatchmap)) of new patch centres coords_p <- SpatialGrid(Patchmap_topo) pts2.sp <- as(pts2, "SpatialPoints") # Extract row number of each new patch centroid pts2_overlay <- sp::over( pts2.sp, coords_p) #---Add new patches to landscape---------# # Loop adds patches 1 at a time # Loop enables varying patch size for (y in 1:length(pts2_overlay)) { patch_A <- paste("patchadd","_",y, sep="") # Create objects # Creates empty grid with correct dimensions SGDF_empty <- SpatialGridDataFrame(Patchmap_topo, data=data.frame(new_h=rep(NA, length(origPatchmap$ID)))) # Adds each new patch centre point SGDF_empty$new_h[pts2_overlay[y]] <- 1 # Coerces object of class "SpatialGridDataFrame" to class "im" Newpop_im <- as(SGDF_empty, "im") # Coerces "im" object to "owin" object Newpop_owin <- as.owin(Newpop_im,win=owin(c(min_x,max_x),c(min_y,max_y))) # Sample size of the new patch to be added # Current protocol is a function of the average and standard deviation of the area (# of cells) of patches within the original landscape # The new radius is sampled from normal distribution with mean = the square root of the mean area of current patches divided by 2, and standard deviation = to a similar function # adjusts radius for cell size #radius<-abs(round(rnorm(1,((avg_area_pop*20^2)^0.5)/2,((sd_area_pop*20^2)^0.5)/2),0)) # generic version: radius <- abs(round(rnorm(1,((avg_area_pop*cellsize^2)^0.5)/2,((sd_area_pop*cellsize^2)^0.5)/2),0)) # re-sample if radius is = 0 if (radius == 0) radius<-abs(round(rnorm(1,(avg_area_pop^0.5)/2,(sd_area_pop^0.5)/2),0)) # Dilates each patch centre based on a randomly sampled radius Newpop_owin_50 <- dilate2.owin(Newpop_owin, r=radius,win=as.mask(mask),include.bb=TRUE)# # Coerces "owin" to "im" object class pts_extract_I <- as.im(Newpop_owin_50) # For the first patch, create a duplicate of pts_extract_I, but with a different name if (patch_A == "patchadd_1") p_im <- pts_extract_I # Combines images of patches creates using eval.im function p_im <- eval.im(p_im|pts_extract_I) } #---END add patches to grid # Coerces "im" to a spatial grid dataframe object class xx2 <- as(p_im, "SpatialGridDataFrame") # SGDF consisting of new patches only #image(xx2) # Coerce 0s to NAs Patchmap$ID[Patchmap$ID==0] <- NA #replace 0s with NA # Convert to raster newp.ra <- raster(xx2)# new patches to be added to landscape patch.ra <- raster(origPatchmap) # original patch map # Add patch maps together new_patchmap <- sum(patch.ra, newp.ra, na.rm = T) names(new_patchmap) = "ID" # Change HSmap to raster object # New patches added do NOT yet have HS values new_HSmap <- raster(HSmap) #plot(new_patchmap) #plot(new_HSmap) #getValues(new_patchmap) #plot(patch.ra) #click(patch.ra) } # END add new patches to landscape #-----------------------------------# #-------------------------------------# # NO CHANGE IN NO. of POPULATIONS # #-------------------------------------# # START number of patches not changed # Code modifies object names #to ensure consistency if (Original_Npops == Npops) { # Create new patchmap object new_patchmap <- Patchmap names(new_patchmap) <- "ID" # Create new hsmap object new_HSmap <- HSmap names(new_HSmap) <- "HSvalue" } # END number of patches not changed #-------------------------------------# # Part 2: # # Vary map habitat suitability # #-------------------------------------# # If patches are being added #if ( Original_Npops < Npops ) { # Function: set values below 100 to NA. #set0toNA <- function(x) { x[x==0] <- NA; return(x) } # Function: set NA values to 0 #setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} #or use this: NAvalue(new_HSmap) <- 0 #---------------------------------------# # HS variation #1 # # Generate landscape of autocorrelated # # HS values # # using package 'randomFields' # #---------------------------------------# if (HSvary == "SACexp"){ # Define variogram exponential model # the model includes nugget effect and the mean: model <- RMexp(var=1, scale=100) #+ #RMnugget(var=1) + # nugget #RMtrend(mean=0) # and mean ## define the locations: # dimensions of landscape celldimX<-gridpara[1,3] celldimY<-gridpara[2,3] from <- 1 x.seq <- seq(from, celldimX, length=celldimX) y.seq <- seq(from, celldimY, length=celldimY) # RFsimulate produces kriged raster based on variogram # RFoptions(spConform=TRUE) # produces output as sp object, slower than FALSE simu <- RFsimulate(model, x=x.seq, y=y.seq) #plot(simu) # Add topology Patchmap_topo <- getGridTopology(origPatchmap) new <- SpatialGridDataFrame(Patchmap_topo, data=data.frame(simu)) newHS <- raster(new) # to raster names(newHS)<- "newHS" #------------------------------------# # Rescale HS values # #------------------------------------# # Extract max/min values of new raster max_HS <- cellStats(newHS, stat='max', na.rm=TRUE) # Rescale HS values to range from new sampled HS threshold # to max HS value of varied HS landscape; # typically this will be 1 newrange <- c(new_HSthreshold, max_HS) rescaleHS <- plotrix::rescale(newHS@data@values, newrange) # extract coordinates newHS from raster coord <- xyFromCell(newHS,1:ncell(newHS)) # extract values # vals <- extract(newHS,1:ncell(newHS)) # Convert dtf to raster # create spatial points data frame spg <- data.frame(coord, rescaleHS) coordinates(spg) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(spg) <- TRUE # coerce to raster rescale.hs <- raster(spg) # habitat suitability surface # Extract patch locations from habitat suitability surface # Reclass cells in patches to a value of '1' #hold <- new_patchmap[[1]] hold <- calc(new_patchmap, function(x) { x[x >= 1] <- 1; return(x) } ) #plot(hold) #plot(new_patchmap) #plot(new_HSmap) # combine rasters by multiplying newHS.patch <- hold * rescale.hs #plot(newHS.patch) #cellStats(newHS.patch, stat='range', na.rm=TRUE) } # END sample HS using 'RandomFields' exponential model #---------------------------------------# #---------------------------------------# # HS variation #2 # # Vary HS map using # # random normal sampling # #---------------------------------------# if (HSvary == "rnorm"){ # Create binary patch map (matrix:0s and patches:1s) # Reclass cells in patches to a value of '1' hold <- calc(new_patchmap, function(x) { x[x >= 1] <- 1; return(x) } ) #plot(hold) # Resample cells based on random normal distribution # No spatial autocorrelation # ncells: maximum number of cells to be sampled = to number of cells in the raster newHS <- calc(hold, function(x) ifelse(x > 0, rnorm( ncell(x), mean_HS, sd_HS), NA) ) #------------------------------------# # Rescale HS values # #------------------------------------# # Extract max/min values of new raster max_HS <- cellStats(newHS, stat='max', na.rm=TRUE) # Rescale HS values to range from new sampled HS threshold # to max HS value of varied HS landscape; # typically this will be 1 newrange <- c(new_HSthreshold, max_HS) rescaleHS <- plotrix::rescale(newHS@data@values, newrange) # extract coordinates 'newHS' from raster coord <- xyFromCell(newHS, 1:ncell(newHS)) # extract values # vals <- extract(newHS,1:ncell(newHS)) # Convert dtf to raster # create spatial points data frame spg <- data.frame(coord, rescaleHS) coordinates(spg) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(spg) <- TRUE # coerce to raster newHS.patch <- raster(spg) # Plots: to delete #par(mfrow = c(2,2)) #plot(r2) #hist(values(r2) ) #plot(newHS.patch) #hist(values(newHS.patch) ) } # END HS sample by random normal distribution, no SAC #---------------------------------------# #---------------------------------------# # HS variation #3 # # Vary HS map using ensemble model # # predictions and uncertainty raster # #---------------------------------------# if (HSvary == "ensemble"){ # Read in map describing variation in ensemble #sdHSmap <- new_HSmap*0.1 sdHSmap <- raster(origHSmap)*0.1 #sdHSmap <- readGDAL(sdHSmap) # Modify sdHSmap to reflect new (lower or same #) Npops if (Original_Npops >= Npops){ # Remove cells from original sdHSmap that were originally # patches, to '0' # Cells were removed from 'new_HSmap' object # during the remove patches routine sdHSmap <- (new_HSmap > 0) * (sdHSmap) + (new_HSmap <= 0)*new_HSmap setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} new_HSmap <- calc(new_HSmap, setNAto0) sdHSmap <- calc(sdHSmap, setNAto0) # Sample new HS values based on a random normal variate # where the mean is the ensemble prediction # and SD is uncertainty of ensemble prediction # Sampled on a cell by cell basis hh <- rnorm(ncell(new_HSmap), getValues(new_HSmap), getValues(sdHSmap)) newHS.patch <- setValues(new_HSmap, hh) # Assign NA to '0' set0toNA <- function(x) { x[x==0] <- NA; return(x) } newHS.patch <- calc(newHS.patch, set0toNA) #plot(newHS.patch) } else { setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} new_HSmap <- calc(new_HSmap, setNAto0) sdHSmap <- calc(sdHSmap, setNAto0) # Vary HS based on mean ensemble prediction and SD of ensemble models # 'new_HSmap' does NOT have new patches added hh <- rnorm(ncell(new_HSmap), getValues(new_HSmap), getValues(sdHSmap)) newHS.patch <- setValues(new_HSmap, hh) # Cells where patches were added to the landscape # need to be assigned HS values # Each HS values (i.e., new cell where patches are added) # is assigned the mean value across the newly sampled HS values # i.e., base on newHS.patch # Users can change this as required # Mean of newly varied HS values # NOTE: mean calculated on any cell > 0 of modified HS landscape mean_newHS <- mean( newHS.patch[[1]][ newHS.patch[[1]] > 0] ) # Extract areas where new patches added but no previous patches existed new_patches_only <- newp.ra - patch.ra # Replace any negative values with NA new_patches_only[new_patches_only < 1] <- NA #unique(new_patches_only@data@values, na.rm=T) #plot(new_patches_only) # Function: set values above 0 to mean HS value setHS <- function(x) { x[x>0] <- mean_newHS; return(x) } new_patches_only_HS <- calc(new_patches_only, setHS) # Function set NA values to 0 setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} new_patches_only_HS <- calc(new_patches_only_HS, setNAto0) # Add new patches that do not overlap original patch cells # and varied HS raster newHS.patch <- new_patches_only_HS + newHS.patch # Assign NA to '0' set0toNA <- function(x) { x[x==0] <- NA; return(x) } newHS.patch <- calc(newHS.patch, set0toNA) } #------------------------------------# # Rescale HS values # #------------------------------------# # Extract max/min values of new raster max_HS <- cellStats(newHS.patch, stat='max', na.rm=TRUE) # Rescale HS values to range from new sampled HS threshold # to max HS value of varied HS landscape; # typically this will be 1 newrange <- c(new_HSthreshold, max_HS) rescaleHS <- plotrix::rescale(newHS.patch@data@values, newrange) #unique(rescaleHS, na.rm=F)[1:10] # extract coordinates 'newHS.patch' from raster coord <- xyFromCell(newHS.patch, 1:ncell(newHS.patch)) # extract values # vals <- extract(newHS,1:ncell(newHS)) # Convert dtf to raster # create spatial points data frame spg <- data.frame(coord, rescaleHS) coordinates(spg) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(spg) <- TRUE # coerce to raster newHS.patch <- raster(spg) #-----End rescale raster #rm(newHS.patch) #plot(newHS.patch) } #-----End vary ensemble HS-----# #---------------------------------------# #---------------------------------------# # HS variation #4 # # HS not varied # #---------------------------------------# if (HSvary == "none"){ if (Original_Npops >= Npops){ newHS.patch <- new_HSmap } else { # Assign HS values to new patches added to the landscape # 'new_HSmap' does NOT have new patches added # Mean of newly varied HS values # NOTE: mean calculated on any cell > 0 of modified HS landscape mean_newHS <- mean( new_HSmap[[1]][ new_HSmap[[1]] > 0] ) # Extract areas where new patches added but no previous patches existed new_patches_only <- newp.ra - patch.ra # Replace any negative values with NA # Negative values overlap with currently existing patches new_patches_only[new_patches_only < 1] <- NA #unique(new_patches_only@data@values, na.rm=T) #plot(new_patches_only) # Function: set values above 0 to mean HS value setHS <- function(x) { x[x>0] <- mean_newHS; return(x) } new_patches_only_HS <- calc(new_patches_only, setHS) # Function set NA values to 0 setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} new_patches_only_HS <- calc(new_patches_only_HS, setNAto0) # Add new patches that do not overlap original patch cells # and varied HS raster newHS.patch <- new_patches_only_HS + new_HSmap # Assign '0' to NA set0toNA <- function(x) { x[x==0] <- NA; return(x) } newHS.patch <- calc(newHS.patch, set0toNA) #plot(new_patches_only_HS) #plot(newHS.patch) } } #---------------------------------------# #par(mfrow=c(3,1)) #plot(new_HSmap) #plot(sdHSmap) #plot(newHS.patch) # click(newHS.patch) #plot(new_patches_only_HS) # click(new_patches_only_HS) #plot(newHS.patch2) #---------------------------------------# # Write new HS map to ascii files # # format required for # # RAMAS GIS spatial module # #---------------------------------------# # Set NA values to 0 setNAto0 <- function(x) { x[is.na(x)] <- 0; return(x)} newHS.patch <- calc(newHS.patch, setNAto0) r <- writeRaster(newHS.patch, filename = newHSmap_asc, datatype='ascii', bylayer = TRUE, suffix="names", overwrite = TRUE ) #---------------------------------------# # Write newly varied .ptc file # # (ptc_"y".ptc) # #---------------------------------------# # Write lines 1 through 7 to the ptcfile write(inputprefile[1], file=ptcfile, append=FALSE) for (i in 2:7) { write(inputprefile[i], file=ptcfile, append=TRUE) } # Line 8: write new Habitat Suitability Function write(new_HSfunction, file =ptcfile, append = TRUE) # Line 9: write a blank line write(inputprefile[9], file =ptcfile, append = TRUE) # Line 10: write new Habitat Suitability Threshold # The subtraction is required as a work around # to ensure that the proper HS threshold is set. # For some reason, after GRIP 2.0 is run, # RAMAS does not read the HS values in # the new HS map as they appear. Perhaps due to # how computers code real numbers. # RAMAS is set up so that any cell in the # HS map >= HS threshold will be included as a patch. # GRIP 2.0 respects this constraint but because # RAMAS reads the HS values differently, # 0.001 is subtracted from the new HS threshold value. # This should not be an issue if HS maps are created # outside of GRIP 2.0 but we advise everyone to # check their models to ensure the accuracy of this statement. new_HSthreshold <- new_HSthreshold-0.001 write(new_HSthreshold , file =ptcfile, append = TRUE) # Line 11: write new Neighbourhood Distance write(line11_nghbd_dist, file =ptcfile, append = TRUE) # Line 12: write habitat suitability map colour write(inputprefile[12], file =ptcfile, append = TRUE) # Line 13: write decimals for habitat suitability map export write(inputprefile[13], file =ptcfile, append = TRUE) # Write lines 14 to 24 for (i in 14:24) { write(inputprefile[i], file=ptcfile, append=TRUE) } # Write if habitat based distances from a friction map write("No", file =ptcfile, append = TRUE) # Write friction map file (should not be specified) write(fric_map, file =ptcfile, append = TRUE) # Line 27: write new method for measuring distances between patches write(line27_distance_measure, file =ptcfile, append = TRUE) # Write the number of input maps included in model write(nu_maps, file =ptcfile, append = TRUE) # Write Habitat Suitability map name write(paste(newHSmap_name), file = ptcfile, append = TRUE) #Write habitat suitability file name write(paste(newHSmap_asc),file = ptcfile, append = TRUE) # Write habitat suitability file type write("ARC/INFO,ConstantMap", file = ptcfile, append = TRUE) # Write remaining lines: from line 32 to end for (i in 32:(Migration+5)) { write(inputprefile[i], file=ptcfile, append=TRUE) } # Write "End of file" write("-End of file-", file = ptcfile, append = TRUE) #---------------------------------------# # END OF PREFILE CREATION # #---------------------------------------# #---------------------------------------# #---------------------------------------# # Collects results from parameters # # varied in the previous code # #---------------------------------------# # Creates an empty dataframe, where the number of columns is = to the # of variables collected PTCVars<-as.data.frame(matrix(,1,6)) # Records the number of sampled populations PTCVars[1]<-Npops # Records the new habitat suitability threshold PTCVars[2]<-new_HSthreshold+0.001 # Records the new neighbourood distance PTCVars[3]<-line11_nghbd_dist # Records and codes for the new distance measure if ( line27_distance_measure == "Default: Center to center") PTCVars[4]<-1 if ( line27_distance_measure == "Center to edge") PTCVars[4]<-2 if ( line27_distance_measure == "Edge to edge") PTCVars[4]<-3 # records HS sample method PTCVars[5] <- HSvary # Calculates and records the size of the sampled radius if populations are added explicitly to the landscape if (Diff < 0) PTCVars[6] <- radius else PTCVars[6]<- NA # Binds and writes variables to PTCvariables.txt data.frame(do.call("cbind",PTCVars)) write.table(PTCVars,file="PTCvariables.txt",append=TRUE,sep=",",col.names = FALSE, row.names= FALSE) gc() } # END generate new replicate landscapes #---------------------------------------# #---------------------------------------# #---------------------------------------# # Initiates the batch file to run # replicate simulations in RAMAS Spatial. # The output of this batch file is a set of # files required to run RAMAS Spatial. # For each replicate, this includes a '*.ptc' # file and a habitat suitability map ('HS_*.ASC'). # Once run in RAMAS Spatial, the outputs # will be used to create and then # vary parameters in RAMAS Metapop. system("batch_ptc.bat", wait = TRUE) #---------------------------------------# #---------------------------------------# #---------------------------------------# # Collates and writes potential # # explanatory variables from the # # results files created by RAMAS Spatial# #---------------------------------------# # Writes the names of the varied input # parameters and other variables to the # data output file called "spatialvar.txt" SpatialVarName <- list("Npops", "TotalHS", "meanTotalHS", "sdTotalHS", "minTotalHS", "maxTotalHS", "meanHSPop", "sdHSPop", "minHSPop", "maxHSPop","totalPatchAreaC", "meanAreaCPop", "sdAreaCPop", "minAreaCPop", "maxAreaCPop", "meanAreaPop", "sdAreaPop", "minAreaPop", "maxAreaPop", "totalPopArea", "meanCoreArea","sdCoreArea","minCoreArea", "maxCoreArea", "meanEdgeAreaRatio", "totalEdge ", "meanEdge", "totalCoreArea", "meanShapeIndex", "sdShapeIndex", "minShapeIndex", "maxShapeIndex", "meanFractalDim", "areaPercLands") write.table(SpatialVarName, file="spatialvar.txt", col.names = FALSE, row.names= FALSE, sep=",", quote=FALSE, append=FALSE) #-----------------------------------------# # Collect results from each # # replicate simulation file # #-----------------------------------------# for(y in 1:Nreps) { # Ensures ptcfile is linked to the right ptc file ptcfile <- paste("ptc", "_",y,".ptc", sep="") # Reads in the replicate ptc file ptcinput <- readLines(ptcfile, n=-1) # Gives the line number in the input file that has "Results" on it Results <- grep("Results", ptcinput) # Gives the line number in the input file that has "Migration" on it Migration <- grep("Migration", ptcinput) # Gives the line number in the input file that has "End of file" on it Endoffile <- grep("End of file", ptcinput) # Reads in the number of populations in the replicate ptc file as a numeric object Npops <- scan(ptcfile, what = "list", skip = Results, nlines = 1) Npops <- as.numeric(Npops[1]) # Reads in the population dataframe pop <- read.table(ptcfile, skip = (Results+1), nrow = Npops,sep=",") # Reads in the population results dataframe pop_data <- read.table(ptcfile, skip = (Results+Npops+3), nrow = Npops, sep="") #V1 total habitat suitability/patch #V2 area (# of cells) #V3 patch minX #V4 patchminY #V5 patchmaxX #V6 patchmaxY #V7 patch # # Reference line numbers LandsIndices <- grep("Landscape indices", ptcinput) Mapaverages <- grep("Map averages", ptcinput) #gives line number that has "Map averages" # Reads in the "Landscape indices" results dataframe LandsResults <- read.table(ptcfile, skip = (LandsIndices), nrow = Npops, sep="") # Reads in one line of results containing a number of Landscape Indices LandsResults2 <- scan(ptcfile, what = "list", skip = (LandsIndices+Npops), nlines = 1) # Reads in one line of results containing a number of Landscape Indices LandsResults3 <- scan(ptcfile, what = "list", skip = (LandsIndices-Npops-2), nlines = 1) # Retrieves cell length (dimensions in km) cell_length <- as.numeric(LandsResults3[1]) # Calculates area of 1 cell cell_area <- cell_length^2 #################################################################################################### # Creates an empty dataframe, where the number of columns is = to the # of variables collected SpatialVars <- as.data.frame(matrix(,1,34)) # Records the number of populations (= number of habitat patches) identified by RAMAS Spatial SpatialVars[1] <- Npops # Npops # Calculates total habitat suitability summed over across all habitat patches SpatialVars[2] <- sum(pop_data$V1) # TotalHS # Calculates the mean of the total HS from each patch SpatialVars[3] <- mean(pop_data$V1) # meanTotalHS # Calculates SD of the total HS from each patch SpatialVars[4] <- sd(pop_data$V1) # sdTotalHS # Calculates minimum value of the total HS from each patch SpatialVars[5] <- min(pop_data$V1) # minTotalHS # Calculates maximum value of the total HS from each patch SpatialVars[6] <- max(pop_data$V1) # maxTotalHS # Calculates mean of the average habitat suitability within each patch SpatialVars[7] <- mean(pop_data$V1/pop_data$V2) # meanHSPop # Calculates sd of the average habitat suitability within each patch SpatialVars[8] <- sd(pop_data$V1/pop_data$V2) # Calculates min of the average habitat suitability within each patch SpatialVars[9] <- min(pop_data$V1/pop_data$V2) # Calculates max of the average habitat suitability within each patch SpatialVars[10] <- max(pop_data$V1/pop_data$V2) # Calculates the total number of cells across all patches (area, but unit is # of cells) SpatialVars[11] <- sum(pop_data$V2) #totalPatchAreaC # number of cells # Calculates mean number of cells within patches SpatialVars[12] <- mean(pop_data$V2) #meanAreaCPop # Calculates standard deviation of the number of cells within patches SpatialVars[13] <- sd(pop_data$V2) #sdAreaCPop # Calculates minimum number of cells within patches SpatialVars[14] <-min(pop_data$V2) #minAreaCPop # Calculates maximum number of cells within patches SpatialVars[15] <- max(pop_data$V2) # maxAreaCPop # Calculates the mean area of patches within landscape (in km^2) SpatialVars[16] <- mean(pop_data$V2)*cell_area # meanAreaPop # Calculates standard deviation of patch areas within landscape SpatialVars[17] <- sd(pop_data$V2)*cell_area # sdAreaPop # Calculates minimum patch area within landscape SpatialVars[18] <- min(pop_data$V2)*cell_area #minAreaPop # Calculates maximum patch area within landscape SpatialVars[19] <- max(pop_data$V2)*cell_area # maxAreaPop # Calculates total patch area within landscape SpatialVars[20] <- sum(pop_data$V2)*cell_area # totalPopArea # Calculates the mean core area (km^2) across all patches (all habitat except edges) SpatialVars[21] <- mean(LandsResults$V4)*cell_area #meanCoreArea # Calculates the sd core area (km^2) across all patches (all habitat except edges) SpatialVars[22] <- sd(LandsResults$V4)*cell_area # Calculates the min core area (km^2) across all patches (all habitat except edges) SpatialVars[23] <- min(LandsResults$V4)*cell_area # Calculates the mean core area (km^2) across all patches (all habitat except edges) SpatialVars[24] <- max(LandsResults$V4)*cell_area # Calculates mean edge to area ratio (1:km^2) across all patches SpatialVars[25] <- mean(LandsResults$V5) # meanEdgeAreaRatio # Calculates total edge length across all populations SpatialVars[26] <- as.numeric(LandsResults2[1]) #totalEdge # Calculates mean edge length across all populations (km) SpatialVars[27] <- as.numeric(LandsResults2[5]) # meanEdge (km) # Calculates total core area across all populations SpatialVars[28] <- sum(LandsResults$V4)*cell_area #totalCoreArea (km^2) # Calculates mean shape index across all populations SpatialVars[29] <- as.numeric(LandsResults2[3]) #meanShapeIndex # Calculates SD of shape index across all patches SpatialVars[30] <- sd(LandsResults$V2) #sdShapeIndex # Calculates minimum shape index across all patches SpatialVars[31] <- min(LandsResults$V2) #minShapeIndex # Calculates maximum shape index across all patches SpatialVars[32] <- max(LandsResults$V2) #maxShapeIndex # Calculates mean fractal dimension across all populations SpatialVars[33] <- as.numeric(LandsResults2[4]) #meanFractalDim # Records the area of patches as a % of the landscape LandDim <- celldimY*celldimX SpatialVars[34] <- round((sum(pop_data$V2)/LandDim)*100,digits=2) #Binds and writes variables to spatialvar.txt files data.frame(do.call("cbind",SpatialVars)) write.table(SpatialVars,file="spatialvar.txt",append=TRUE,sep=",",col.names = FALSE, row.names= FALSE) } #---------------------------------------# # Code to create the Metapop file from # # RAMAS Spatial results # #---------------------------------------# # Creates a loop to run through each # '*.ptc' replicate simulation file and # the Metapop link file to create a # new Metapop file used to vary other # spatial and non-spatial parameters # This code is required for RAMAS to run # as the Metapop file cannot be automatically # saved and exported in batch format in RAMAS GIS for (y in 1:Nreps){ pref_metapop <- paste("pref", "_",y,".mp", sep="") # Ensures ptcfile is linked to the right replicate ptc file ptcfile <- paste("ptc", "_",y,".ptc", sep="") # Reads in the replicate ptc file ptcinput <- readLines(ptcfile, n=-1) # Creates an object of the metapop file used as the link to the metapop module in RAMAS spatial module metapop_link <- ptcinput[20] # Reads in the whole link to RAMAS Metapop (*.mp) file as the object inputprefile metapop_input <- readLines(metapop_link, n=-1) # Finds the line that Results is on, in the replicate ptc file Results <- grep("Results", ptcinput) # Gives the line number in the input file that has "Migration" on it Migration_linkmp <- grep("Migration", metapop_input) # Reads in the original number of populations in the Metapop Link model Original_Npops <- (Migration_linkmp-41)-4 # Writes the first 44 lines from the .mp used as the link to metapop file in ramas spatial write(metapop_input[1], file=pref_metapop, append=FALSE) for (i in 2:32) { write(metapop_input[i], file=pref_metapop, append=TRUE) } if(metapop_input[34]=="UD") { write("No", file=pref_metapop, append=TRUE) } else { write("Yes", file=pref_metapop, append=TRUE) } for (i in 34:44) { write(metapop_input[i], file=pref_metapop, append=TRUE) } # Read in # of populations defined by Spatial module, in *.ptc file (just run in RAMAS spatial module in the batch mode) Npops<-scan(ptcfile, what = "list", skip = Results, nlines = 1) Npops<-as.numeric(Npops[1]) # Reads in the population dataframe from the replicate .ptc file pop<-read.table(ptcfile,skip=Results+1,sep=",",nrow=Npops) # Writes the population dataframe write.table(pop,file=pref_metapop,append=TRUE,sep=",",row.names=FALSE,col.names=FALSE,na="",quote=FALSE) # Writes the lines corresponding to Migration and Correlation: includes title, distance function # Reads this information from the new .ptc file, and therefore the user MUST specify the distance function in the original .ptc file, otherwise this will not be read in #Alternatively, users could change the code to read from the link to metapop file Migration<-grep("Migration", ptcinput) # Gives the line number in input file that has "Migration" on it Migration<-grep("Migration", ptcinput) #Gives the line number in input file that has "Correlation" on it Correlation<-grep("Correlation", ptcinput) # Reads in the dispersal-distance parameters M_parameters<-scan(ptcfile, skip = (Migration[1]+1), nlines = 1, sep=",") # Gives the line number in Spatial input file (*.ptc) that has "Landscape_indices" on it Landscape_indices<-grep("Landscape indices", ptcinput) # Creates an empty distance matrix Pairwise_Distance <- array(0, c(Npops,Npops)) # Reads in distance matrix; part of spatial .ptc results if ( ptcinput[27] == "Center to edge") { Pairwise_Distance <- matrix(scan(ptcfile, nlines = Npops, skip = (Landscape_indices-Npops-1)),ncol = Npops, byrow = TRUE) } else { # Scans in pairwise distance matrix (only lower half exists) from varied Spatial Data file, which includes results Pairwise_Distance[upper.tri(Pairwise_Distance, diag = TRUE)] <- scan(ptcfile, nlines = Npops, skip = (Landscape_indices-Npops-1)) # Creates a full matrix of pairwise distances among populations, where the lower and upper halves of the matrix are symmetrical Pairwise_Distance <- Pairwise_Distance + t(Pairwise_Distance) } # Creates an empty dispersal matrix Dispersal_matrix <- array(0, c(Npops,Npops)) # Adjusts the dispersal rates according to the distances among population pairs calculated by the Spatial program # Distances could be calculated based on either center-center (default), edge-edge, or center-edge distances for (i in 1:(Npops)) { for (j in 1:(Npops)) { Dispersal_matrix[i,j]<-(M_parameters[1]*exp((-Pairwise_Distance[i,j]^M_parameters[3])/M_parameters[2])) if(M_parameters[4]>0) { if(Pairwise_Distance[i,j]>M_parameters[4]) { Dispersal_matrix[i,j]<-0 } } } } # Also enforces a maximum distance in this code, setting dispersal to 0 for population pairs separated by more than this distance # Ensures diagonal elements on the dispersal matrix are 0 diag(Dispersal_matrix)<-c(rep(0,length(Dispersal_matrix[1,]))) # Writes "Migration" on a line write("Migration", file=pref_metapop,append=TRUE) # Writes the line after "Migration" write("FALSE", file=pref_metapop,append=TRUE) # Writes Dispersal distance function parameters write(ptcinput[Migration+2], file = pref_metapop, append = TRUE) # Writes Dispersal Matrix Extra_NA<-c(rep(NA, length(Dispersal_matrix[1,]))) write.table((cbind(Dispersal_matrix, Extra_NA, deparse.level=0)), file = pref_metapop, sep = ",", append = TRUE, quote=FALSE, row.names=FALSE, col.names=FALSE, na="") # Reads in the correlation matrix, if it's there or creates one with the correct dimensions D_autofill<-scan(ptcfile, what = "list", skip = (Correlation[1]), nlines = 1) D_autofill<-as.logical(D_autofill) if (any(D_autofill)) { Correlation_matrix<-matrix(0,Npops,Npops) } else { Correlation_matrix<-matrix(scan(ptcfile, sep = ",", nlines = Npops, skip =Correlation+2), nrow = Npops, byrow = TRUE) Correlation_matrix<-Correlation_matrix[,-(length(Correlation_matrix[1,]-1))] } # Writes "Correlation" on a line write("Correlation", file=pref_metapop,append=TRUE) # Writes the line after "Correlation" write("FALSE", file=pref_metapop,append=TRUE) # Writes the correlation distance function parameters write(ptcinput[Migration+5], file=pref_metapop, append=TRUE) diag(Correlation_matrix)<-c(rep(1, length(Correlation_matrix[1,]))) Correlation_matrix<-cbind(Correlation_matrix, Extra_NA, deparse.level=0) write.table(Correlation_matrix, file = pref_metapop, sep = ",", append = TRUE, quote=FALSE, row.names=FALSE, col.names=FALSE, na="") # Gives the line number in Metapop link file that has "Stage Matrix" on it Line_N_stage_matrices<-grep("stage matrix", metapop_input) # Gives the line number in Metapop link file that has "Constraints Matrix" on it Line_Constraints<-grep("Constraints", metapop_input) # Writes the stage and standard deviation information from the Metapop file link file (used as an input to RAMAS Spatial) for (i in Line_N_stage_matrices:(Line_Constraints-1)) { write(metapop_input[i], file = pref_metapop, append = TRUE) } # Finds the line that specifies the number of stages modeled N_stages<-scan(metapop_link, what = "list", skip = 9, nlines = 1) N_stages<-as.numeric(N_stages[1]) # Writes in the Constraints Matrix, relative dispersal indices and catastrophe multipliers for (i in Line_Constraints:(Line_Constraints+(N_stages[1]*3+3))) { write(metapop_input[i], file = pref_metapop, append = TRUE) } # Read in the initial abundances for each population from the population dataframe from the replicate .ptc file # This does not take into account the distribution by stages Initial_Abundances<-pop$V4 # Reads in the stage-specific initial abundances from the link metapop file # RAMAS Spatial assumes that the initial abundance distribution by stage is the same for all populations, unless these are specfied manually # However, if > 1 population is specified in the Metapop Link .mp model, the following code will calculate an average among all populations included Initial_Abundances_dist<-matrix(scan(metapop_link, skip=(Line_Constraints+(N_stages[1]*3+3)), nlines=Original_Npops), ncol=N_stages[1], byrow=TRUE) # Calculates the average stage structure (i.e. proportion of initial individuals in each stage at the beginning of the simulation) Mean_stage_abundances<-apply(Initial_Abundances_dist, 2, mean, na.rm = TRUE) Prop_stage_abundances<-Mean_stage_abundances/sum(Mean_stage_abundances) # Pastes population initial abundances into stage-specific abundances lines if only one stage is modeled if(N_stages[1]==1) { for (i in 1:Npops) { Initial_Abundances[i]<-pop$V4[i] } } # Calculates the initial stage-specific abundances, if more than one stage is modeled if(N_stages[1]>1) Initial_Abundances<-matrix(0,Npops,N_stages[1]) { for (i in 1:Npops) { Initial_Abundances[i,]<-Prop_stage_abundances*pop$V4[i] #for dispersal mortality population; Initial_Abundances[Npops+1,]<-c(rep(0, N_stages[1])) } Initial_Abundances<-round(Initial_Abundances, digits=0) # Ensures that the number if individuals in the initial abundances by stage matrix is the same as the initial abundances in the population dataframe after rounding To_add<-vector("numeric", length = Npops) for (i in 1:Npops) { To_add[i]<-pop$V4[i]-(sum(Initial_Abundances[i,])) } for (i in 1:Npops) { Initial_Abundances[i,1]<-Initial_Abundances[i,1]+To_add[i] } } # Write the stage-specific initial abundances matrix write.table(Initial_Abundances, file=pref_metapop, append = TRUE, col.names = FALSE, row.names=FALSE) # Reads and writes the information from the stages menu # Original number of populations and number of stages are from the link metapop file for (i in (Line_Constraints+1+(N_stages[1]*3+3)+Original_Npops):(Line_Constraints+(N_stages[1]*3+3)+Original_Npops+5*N_stages[1])) { write(metapop_input[i], file =pref_metapop, append = TRUE) } # Writes that there is no population management write("0 (pop mgmnt)", file = pref_metapop, append = TRUE) # Some users will want to include population management # Writes the extinction and explosion thresholds and timestep datum # Writes a generic extinction threshold of zero write("0",file = pref_metapop, append = TRUE) # Writes a generic explosion threshold of zero write("0",file = pref_metapop, append = TRUE) # Writes a generic timestep datum equal to 3.000 (as required for 100 timesteps) write("3", file = pref_metapop, append = TRUE) # Writes end-of-file write ("-End of file-", file = pref_metapop, append = TRUE) gc() } # END OF METAPOP FILE CREATION: pref_rep*.mp #-------------------------------------------# #-------------------------------------------# # CODE TO VARY METAPOP FILE PARAMETERS # #-------------------------------------------# # Empty list in which to save the stage names for survival and fecundity rates Vital_rates_names<-list() # Reads in the whole original PVA input file as the object metapop_input metapop_input<-readLines(metapop_link, n=-1) #Note: used for referencing lines with the GREP function and for writing lines that are not changed by GRIP. # Finds the line that specifies the number of stages modeled N_stages<-scan(metapop_link, what = "list", skip = 9, nlines = 1) N_stages<-as.numeric(N_stages[1]) # Writes the list of vital rates as headings for the data output file called "vital_rates.txt" if(N_stages[1]>1) { for (i in 1:N_stages[1]) { Vital_rates_names[[i]]<-paste("Fec",i, sep="") Vital_rates_names[[i+N_stages[1]]]<-paste("Surv",i,sep="") } } else { Vital_rates_names[[1]]<-"Pop_growth" } write.table(Vital_rates_names, file="vital_rates.txt", col.names = FALSE, row.names= FALSE, sep=",", quote=FALSE, append=FALSE) # Writes the names of the varied input parameters and other parameters potentially of interest to the data output file called "variables.txt" VarNames<-list("Npop","mean_K","sd_K","max_K","min_K","mean_disCC","max_disCC","min_disCC","mean_abun","sd_abun","max_abun","min_abun","ext_1","ext_2","cat_corr","prob1","local_mult_new","Cat1Abmult","Cat1VRmult","Prob2","local_mult_new2","Cat2Abmult","Cat2VRmult","Transform_dispersal","Connections","Dispersal_surv","Transform_correlation","MCP", "N_sinks", "Mean_dispersal_rate", "min_dispersal_rate", "max_dispersal_rate", "Mean_Correlations", "Min_Correlations", "Max_Correlations", "Mean_connected_dis", "Min_connected_dis", "Max_connected_dis","Selected_Rmax", "meanDis_DistCC_CE_EE", "maxDis_DistCC_CE_EE", "minDis_DistCC_CE_EE") write.table(VarNames, file="variables.txt", col.names = FALSE, row.names= FALSE, sep=",", quote=FALSE, append=FALSE) #-----------------------------------------# # Reads the batch.txt # #-----------------------------------------# # Reads the batch.txt file required for creating the batch file for running input files in RAMAS Metapop 6.0 Start <- readLines( "batch_v6.txt", n = -1L) #-----------------------------------------# # Varies metapop model parameters # #-----------------------------------------# # Changes the name of each replicate simulation input file. For example, the first replicate simulation is called "rep_1.mp" # Loops through all the code from 1 to Nreps (number of replicate simulations, specified above) for (y in 1:Nreps){ # Interim .mp file prefile <- paste("interim", "_",y,".mp", sep="") # Pre-final .mp file pref_metapop <- paste("pref", "_",y,".mp", sep="") # Final replicate .mp file filename<-paste("rep", "_",y,".mp", sep="") # Links to the correct .ptc file # created in the previous script section ptcfile <- paste("ptc", "_",y,".ptc", sep="") # Object for folder location of reps # Change to forward slashes (gsub fx) file.loc <- getwd() file.loc2 <- gsub("/", "\\", file.loc, fixed = TRUE) # adds '\\' for dos batch file # Writes a batch file to run replicate simulation input files in RAMAS Metapop batch_rep <- paste( Start[1], " \"", Start[2], "\" ", " \"", Start[3], "\" ", " \"", "rep_", y, ".mp", "\" ", Start[4], "\"", file.loc2, "\" ", Start[5], sep="" ) write(batch_rep, file= "batch_file.bat", append = TRUE) # Ensures results are saved in text files # that correspond to appropriate replicate input files Abund <- paste( "rename", " \"", file.loc2, "\\Abund.txt", "\" ", "Abund_", y, ".txt", sep="") FinalStageN <- paste( "rename", " \"", file.loc2, "\\FinalStageN.txt", "\" ", "FinalStageN_", y, ".txt", sep="") Harvest <- paste( "rename", " \"", file.loc2, "\\Harvest.txt", "\" ", "Harvest_", y, ".txt", sep="") HarvestRisk <- paste( "rename", " \"", file.loc2, "\\HarvestRisk.txt", "\" ", "HarvestRisk_", y, ".txt", sep="") IntExpRisk <- paste( "rename", " \"", file.loc2, "\\IntExpRisk.txt", "\" ", "IntExpRisk_", y, ".txt", sep="") IntExtRisk <- paste( "rename", " \"", file.loc2, "\\IntExtRisk.txt", "\" ", "IntExtRisk_", y, ".txt", sep="") IntPerDec <- paste( "rename", " \"", file.loc2, "\\IntPerDec.txt", "\" ", "IntPerDec_", y, ".txt", sep="") LocalOcc <- paste( "rename", " \"", file.loc2, "\\LocalOcc.txt", "\" ", "LocalOcc_", y, ".txt", sep="") LocExtDur <- paste( "rename", " \"", file.loc2, "\\LocExtDur.txt", "\" ", "LocExtDur_", y, ".txt", sep="") MetapopOcc <- paste( "rename", " \"", file.loc2, "\\MetapopOcc.txt", "\" ", "MetapopOcc_", y, ".txt", sep="") QuasiExp <- paste( "rename", " \"", file.loc2, "\\QuasiExp.txt", "\" ", "QuasiExp_", y, ".txt", sep="") QuasiExt <- paste( "rename", " \"", file.loc2, "\\QuasiExt.txt", "\" ", "QuasiExt_", y, ".txt", sep="") TerExpRisk <- paste( "rename", " \"", file.loc2, "\\TerExpRisk.txt", "\" ", "TerExpRisk_", y, ".txt", sep="") TerExtRisk <- paste( "rename", " \"", file.loc2, "\\TerExtRisk.txt", "\" ", "TerExtRisk_", y, ".txt", sep="") TerPerDec <- paste( "rename", " \"", file.loc2, "\\TerPerDec.txt", "\" ", "TerPerDec_", y, ".txt", sep="") write(Abund, "batch_file.bat", append = TRUE) write(FinalStageN, "batch_file.bat", append = TRUE) write(Harvest, "batch_file.bat", append = TRUE) write(HarvestRisk, "batch_file.bat", append = TRUE) write(IntExpRisk, "batch_file.bat", append = TRUE) write(IntExtRisk, "batch_file.bat", append = TRUE) write(IntPerDec, "batch_file.bat", append = TRUE) write(LocalOcc, "batch_file.bat", append = TRUE) write(LocExtDur, "batch_file.bat", append = TRUE) write(MetapopOcc, "batch_file.bat", append = TRUE) write(QuasiExp, "batch_file.bat", append = TRUE) write(QuasiExt, "batch_file.bat", append = TRUE) write(TerExpRisk, "batch_file.bat", append = TRUE) write(TerExtRisk, "batch_file.bat", append = TRUE) write(TerPerDec, "batch_file.bat", append = TRUE) # Reads in the replicate ptc file ptcinput <- readLines(ptcfile, n=-1) # Creates an object of the metapop file used as the link to the metapop module in RAMAS spatial module metapop_link<-ptcinput[20] # Reads in the whole link to RAMAS Metapop (*.mp) file as the object inputprefile metapop_input<-readLines(metapop_link, n=-1) # Finds the line that Results is on, in the replicate ptc file Results<-grep("Results", ptcinput) # Gives the line number in the Metapop Link file that has "Migration" on it Migration_linkmp <- grep("Migration", metapop_input) # Reads in the number of populations in the replicate ptc file as a numeric object Npops_ptc <- scan(ptcfile, what = "list", skip = Results, nlines = 1) Npops_ptc <- as.numeric(Npops_ptc[1]) # Reads in the input file as an object InputFile<-readLines(pref_metapop, n=-1) # Gives the line number in the input file that has "Migration" on it Migration<-grep("Migration", InputFile) # Gives the line number in input file that has "Correlation" on it Correlation<-grep("Correlation", InputFile) # Calculates the number of rows in the dispersal matrix, if it exists this will be > 0 Nrows_dispersal_matrix<-(Correlation-Migration-3) # Gives the line number in input file that has "Constraints Matrix" on it Line_Constraints<-grep("Constraints", InputFile) # Finds the line that specifies the number of stages modeled N_stages<-scan(pref_metapop, what = "list", skip = 9, nlines = 1) N_stages<-as.numeric(N_stages[1]) # Reads in the number of populations in the PVA model; this is the new number of populations resulting from the varied .ptc file # This object was specified earlier, but need to re-initialize it here for every .mp rep file that is read and subsequently varied Npops <- (Migration-41)-4 gc() ####################################################################################################### # The following series of code (to line XXXX) adjusts the population dataframe, dispersal matrix, and correlation matrix for the addition of a dispersal mortality sink population; and initializes some objects that will be varied later in the program #Read in from Spatial data metapop file #Reads in the population dataframe pop <- read.table(pref_metapop,skip=44,sep=",",nrow=Npops) # Calculates means and standard deviations that will be required to parameterize sink population # X-coordinate Mean_x<-mean(pop$V2, na.rm = TRUE) SD_x<-sd(pop$V2, na.rm = TRUE) # Y-coordinate Mean_y<-mean(pop$V3, na.rm = TRUE) SD_y<-sd(pop$V3, na.rm = TRUE) # Rmax Mean_Rmax<-mean(pop$V6, na.rm = TRUE) SD_Rmax<-sd(pop$V6, na.rm = TRUE) # Creates a vector of Allee parameters Allee<-c(rep(0,Npops)) for (i in 1:Npops) { Allee[i]<-pop$V9[i] } # Copies the specified model of density-dependence to all populations, if this model is not population-specific #pop$V5<-as.character(pop$V5) #if(InputFile[34]=="UD") # { # } else # { # if(InputFile[33]=="No") # { # pop$V5<-c(rep(InputFile[34], length(pop$V5))) # } # } # Creates lists of values for any user-defined parameters Extra_parameters<-list() N_extras<-0 if(length(pop[1,])>25) { N_extras<-(length(pop[1,])-25) for(i in 1:N_extras) { Extra_parameters[[i]]<-pop[,24+i] } } # Determines whether or not there are temporal trends files for K Levels_10<-0 if(length(levels(pop$V10))>0) { Levels_10<-length(levels(pop$V10)) pop$V10<-as.character(pop$V10) } # Creates a vector of the trends in K from which to sample values below K_trend<-c(rep(0,Npops)) for (i in 1:Npops) { K_trend[i]<-pop$V10[i] } # Determines whether or not there are temporal trends files for the local population multiplier of catastrophe 1 Levels_12<-0 if(length(levels(pop$V12))>0) { Levels_12<-length(levels(pop$V12)) pop$V12<-as.character(pop$V12) } # Creates a vector of the local population multiplier for catastrophe 1 lpm<-c(rep(0,Npops)) for (i in 1:Npops) { lpm[i]<-pop$V12[i] } # Determines whether or not there are temporal trends files for the probability of catastrophe 1 Levels_13<-0 if(length(levels(pop$V13))>0) { Levels_13<-length(levels(pop$V13)) pop$V13<-as.character(pop$V13) } # Creates a vector of the probabilities of catastrophe 1 lpc<-c(rep(0,Npops)) for (i in 1:Npops) { lpc[i]<-pop$V13[i] } # Determines whether or not there are temporal trends files for the relative fecundities Levels_16<-0 if(length(levels(pop$V16))>0) { Levels_16<-length(levels(pop$V16)) pop$V16<-as.character(pop$V16) } # Creates a vector of the relative fecundities rel_fec<-c(rep(0,Npops)) for (i in 1:Npops) { rel_fec[i]<-pop$V16[i] } # Determines whether or not there are temporal trends files for the relative survivals Levels_17<-0 if(length(levels(pop$V17))>0) { Levels_17<-length(levels(pop$V17)) pop$V17<-as.character(pop$V17) } # Creates a vector of the relative survivals rel_surv<-c(rep(0,Npops)) for (i in 1:Npops) { rel_surv[i]<-pop$V17[i] } # Determines whether or not there are temporal trends files for the local population multiplier of catastrophe 2 Levels_19<-0 if(length(levels(pop$V19))>0) { Levels_19<-length(levels(pop$V19)) pop$V19<-as.character(pop$V19) } # Creates a vector of the local population multipliers for catastrophe 2 lpm2<-c(rep(0,Npops)) for (i in 1:Npops) { lpm2[i]<-pop$V19[i] } # Determines whether or not there are temporal trends files for the probability of catastrophe 2 Levels_20<-0 if(length(levels(pop$V20))>0) { Levels_20<-length(levels(pop$V20)) pop$V20<-as.character(pop$V20) } # Creates a vector of the probabilities of catastrophe 2 lpc2<-c(rep(0,Npops)) for (i in 1:Npops) { lpc2[i]<-pop$V20[i] } # Reads in the dispersal matrix # The code is set up to read in the dispersal matrix if it is there, and if not a new empty matrix is created; # There should always be one there as it is reading from the new varied .ptc file if (Nrows_dispersal_matrix>0) { Dispersal_matrix<-matrix(scan(pref_metapop, sep=",", nlines = (Correlation-Migration-3), skip = Migration+2), ncol = (Correlation-Migration-2), byrow = TRUE) Dispersal_matrix<-Dispersal_matrix[,-(length(Dispersal_matrix[1,]-1))] } # If the dispersal matrix is not there, creates an empty dispersal matrix with the correct dimensions if (Nrows_dispersal_matrix==0) { Dispersal_matrix<-matrix(0,Npops,Npops) } # Reads in the correlation matrix, if it's there or creates one with the correct dimensions D_autofill<-scan(pref_metapop, what = "list", skip = (Correlation[1]), nlines = 1) D_autofill<-as.logical(D_autofill) if (any(D_autofill)) { Correlation_matrix<-matrix(0,Npops,Npops) } else { Correlation_matrix<-matrix(scan(pref_metapop, sep = ",", nlines = Npops, skip =Correlation+2), nrow = Npops, byrow = TRUE) Correlation_matrix<-Correlation_matrix[,-(length(Correlation_matrix[1,]-1))] } # Reads in the stage-specific initial abundances Initial_Abundances <- matrix(scan(pref_metapop, skip=(Line_Constraints+(N_stages[1]*3+3)), nlines=Npops), ncol=N_stages[1], byrow=TRUE) # Reads the stage and standard deviation matrices and their respective descriptions #Gives the line on which the number of stages matrices is listed Line_N_stage_matrices<-grep("stage matrix", InputFile) N_stage_matrices<-scan(pref_metapop, skip=(Line_N_stage_matrices-1), nlines = 1, what="list") #Creates a numeric object that represents the number of stage and standard deviation matrices N_matrices<-as.numeric(N_stage_matrices[1]) # Reads in the 4-line descriptions of the Stage_matrices and saves them in a vector of lists # Reads in the stage matrices and saves them in a list Description_Stage_matrix<-vector("list", N_matrices) Stage_matrices<-list() for (i in 1:N_matrices) { Description_Stage_matrix[[i]]<-scan(pref_metapop, skip = (Line_N_stage_matrices+((i-1)*4)+(i-1)*N_stages[1]), nlines = 4, what = "list", sep= ",") Stage_matrices[[i]]<-matrix(scan(pref_metapop, skip = (Line_N_stage_matrices+4+((i-1)*4)+((i-1)*N_stages[1])), nlines = (N_stages[1])), N_stages[1], N_stages[1], byrow = TRUE) } # Gives the line on which the number of standard deviation matrices is listed Line_N_stdev_matrices<-grep("st.dev. matrix", InputFile) # Gives the information on this line N_stdev_matrices<-scan(pref_metapop, skip=(Line_N_stdev_matrices-1), nlines = 1, what="list") # Reads in the 1-line descriptions of the Stdev_matrices and saves them in a vector of lists # Reads in the standard deviation matrices and saves them in a list Description_stdev_matrix<-vector("list", N_matrices) Stdev_matrices<-list() for (i in 1:N_matrices) { Description_stdev_matrix[[i]]<-scan(pref_metapop, skip = (Line_N_stdev_matrices+((i-1))+(i-1)*N_stages[1]), nlines = 1, what = "list", sep= ",") Stdev_matrices[[i]]<-matrix(scan(pref_metapop, skip = (Line_N_stdev_matrices+1+(i-1)+((i-1)*N_stages[1])), nlines = (N_stages[1])), N_stages[1], N_stages[1], byrow = TRUE) } # Add dispersal mortality (i.e. sink) population # Creates an extra row for simulating dispersal mortality in the population dataframe (i.e. sink population) pop <- rbind(pop, matrix(, 1, length(pop[1,])), deparse.level = 0) # Creates an extra row and column in the dispersal matrix Dispersal_matrix <- rbind(Dispersal_matrix, matrix(0, 1, length(Dispersal_matrix[1,]))) Dispersal_matrix<-cbind(Dispersal_matrix, matrix(0, length(Dispersal_matrix[,1]),1)) # Creates an extra row and column in the correlation matrix Correlation_matrix<-rbind(Correlation_matrix, matrix(0, 1, length(Correlation_matrix[1,]))) Correlation_matrix<-cbind(Correlation_matrix, matrix(0, length(Correlation_matrix[,1]),1)) # Creates an extra population for the stage-specific abundances matrix Initial_Abundances<-rbind(Initial_Abundances, matrix(0, 1, N_stages[1]), deparse.level=0) # Assigns input parameter values in the population dataframe for the dispersal mortality population # Creates population names for all populations, including sink population pop$V1<-as.character(pop$V1) for(i in 1:Npops) { pop$V1[i]<-paste("Population_",i, sep="") } pop$V1[Npops+1]<-"Dispersal_Mortality" # Assigns new spatial coordinates to the new sink population # Coordinates are arbitrary values as the sink population is used to track dispersal mortality pop$V2[Npops+1]<-rnorm(1,Mean_x, 0.1*Mean_x) pop$V3[Npops+1]<-rnorm(1, Mean_y, 0.1*Mean_y) # Rounds coordinates to 3 decimal places pop$V2<-round(pop$V2, digits=3) pop$V3<-round(pop$V3, digits=3) # Varies initial population sizes (assumes normal distribution) and sets constraints # Reads in habitat specific initial abundance function (line 16) in_abund<-scan(ptcfile,what="numeric",skip=15,nlines=1) # Reads in the original number of populations in the Metapop Link model Original_Npops<-(Migration_linkmp-41)-4 # Reads in the population dataframe from the replicate .ptc file pop_link<-read.table(metapop_link,skip=44,sep=",",nrow=Original_Npops) # The initial abundances for each population will be taken from the RAMAS Spatial results (i.e. the replicate .ptc file) and the following code will adjust the initial abundance distribution by stages # If initial abundances are not a function of habitat but instead are constants (either a 0 or 1) then initial abundances will be sampled from the Metapop Link model. The user should have specified more than 1 population in the Metapop Link file in this case. The new abundance will be sampled using the mean and SD of initial abundances across all populations; a 10% CV is applied if there is SD = 0 # If initial abundances are a function of habitat, initial abundances were calculated by RAMAS Spatial and new initial abundances are sampled using the initial abundance (for that population) as the mean and a 10% CV is applied to each # Calculates mean and SD of initial abundances, if initial abundances are not a function of habitat; sets SD to 0.1 if SD = 0 Mean_abundances<-mean(pop_link$V4, na.rm = TRUE) if (length(pop_link$V4)==1) { SD_abundances<-(Mean_abundances*0.1) } else { SD_abundances<-sd(pop_link$V4, na.rm = TRUE) if(SD_abundances==0) { SD_abundances<-(Mean_abundances*0.1) } } # NOTE: If no variability in initial population sizes, a 10% CV is assumed # code to sample initial abundances if it is not a function of habitat (i.e. it was set as a constant in the original ptc model) if (in_abund==0|in_abund==1) { for (i in 1:Npops) { pop$V4[i]<-as.integer(rnorm(1, mean=Mean_abundances, sd=SD_abundances)) if(pop$V4[i]<0) { pop$V4[i]<-0 } # Sets constraints for initial abundances } # Sets dispersal mortality initial abundance to 0 pop$V4[Npops+1]<-0 # Rounds initial abundances to 0 decimal places pop$V4<-round(pop$V4, digits=0) } else { for (i in 1:Npops) { pop$V4[i]<-as.integer(rnorm(1, mean=pop$V4[i], sd=0.1*pop$V4[i])) if(pop$V4[i]<0) { pop$V4[i]<-0 } # Sets constraints for initial abundances } # Sets dispersal mortality initial abundance to 0 pop$V4[Npops+1]<-0 # Rounds initial abundances to 0 decimal places pop$V4<-round(pop$V4, digits=0) } # Pastes population initial abundances into stage-specific abundances lines if only one stage is modeled if(N_stages[1]==1) { for (i in 1:Npops) { Initial_Abundances[i]<-pop$V4[i] } } # Calculates the initial stage-specific abundances, if more than one stage is modeled if(N_stages[1]>1) { for (i in 1:Npops) { Initial_Abundances[i,]<-Prop_stage_abundances*pop$V4[i] Initial_Abundances[Npops+1,]<-c(rep(0, N_stages[1])) } Initial_Abundances<-round(Initial_Abundances, digits=0) # Ensures that the number if individuals in the initial abundances by stage matrix is the same as the initial abundances in the population dataframe after rounding. To_add<-vector("numeric", length = Npops) for (i in 1:Npops) { To_add[i]<-pop$V4[i]-(sum(Initial_Abundances[i,])) } for (i in 1:Npops) { Initial_Abundances[i,1]<-Initial_Abundances[i,1]+To_add[i] } } gc() ######################################################################## ##Samples a new model of density dependence DD_Models<-c("CE", "BH", "LO") #creates a vector of models of density dependence to sample from #Copies the model of density dependence to all populations except sink population pop$V5<-as.character(pop$V5) Selected_DD<-sample(DD_Models,1) for(i in 1:Npops) { pop$V5[i]<-Selected_DD } pop$V5<-as.character(pop$V5) #coerces proper object class # Adjusts model of DD depending on whether the original type was user-specified (DD) or exponential (the latter adjustment is needed to ensure error-free coding) for(i in 1:Npops) { if(InputFile[34]=="UD") { } else { if(pop$V5[i]=="EX") { pop$V5[i]<-"CE" pop$V7[i]<-500000000 } } } # Adjusts dispersal mortality population model of DD if(InputFile[34]=="UD") { } else { pop$V5[Npops+1]<-"CE" } # Assigns a value of Rmax to dispersal mortality population (assumes Rmax >=1, and samples from a normal distribution) pop$V6[Npops+1]<-rnorm(1, Mean_Rmax, SD_Rmax) if(pop$V6[Npops+1]<1) { pop$V6[Npops+1]<-1 } pop$V6<-round(pop$V6, digits=2) #Ensures Rmax values are rounded to 1 decimal places # Arbitrarily sets maximum K of sink population to arbitrarily low number pop$V7[Npops+1]<-2 # Ensures values are rounded to 0 decimal places pop$V7<-round(pop$V7, digits=0) # Arbitrarily sets standard deviation in K for sink population to 0 pop$V8[Npops+1]<-0 # Randomly assigns Allee parameter value to the sink population pop$V9[Npops+1]<-sample(Allee,1) # Randomly assigns a temporal trend in K to the sink population if(Levels_10>0) { pop$V10[Npops+1]<-sample(K_trend,1) } else { pop$V10[Npops+1]<-"" } # Assigns a 0 slope for density-dependent dispersal as a function of source pop size pop$V11[Npops+1]<-0 # Randomly samples a catastrophe 1 local multiplier for the sink population if(Levels_12>0) { pop$V12[Npops+1]<-sample(lpm,1) } else { pop$V12[Npops+1]<-1 } # Randomly samples a local probability of catastrophe 1 for sink population if(Levels_13>0) { pop$V13[Npops+1]<-sample(lpc,1) } else { pop$V13[Npops+1]<-0 } # Ensures all populations are included in the summation of individuals, except for the sink population for (i in 1:Npops) { pop$V14[i]<-"TRUE" } pop$V14[Npops+1]<-"FALSE" # Arbitrarily sets stage matrix for sink population to the first matrix pop$V15[Npops+1]<-1 # Rounds values to 0 decimal places pop$V15<-round(pop$V15, digits=0) # Assigns a fecundity multiplier to sink population pop$V16[Npops+1]<-1 # Assigns a survival multiplier to sink population pop$V17[Npops+1]<-1 # Assigns a local threshold of 0 to sink population pop$V18<-c(rep(0, Npops+1)) # Randomly samples a catastrophe 2 local multiplier for the sink population if(Levels_19>0) { pop$V19[Npops+1]<-sample(lpm2,1) } else { pop$V19[Npops+1]<-1 } if(Levels_20>0) { pop$V20[Npops+1]<-sample(lpc2,1) } else { pop$V20[Npops+1]<-0 } pop$V20[Npops+1]<-0 pop$V22[Npops+1]<-500000000 pop$V23[Npops+1]<-0 pop$V24[Npops+1]<-0 if(length(pop[1,])>25) { for(i in 1:N_extras) { pop[Npops+1,24+i]<-sample(Extra_parameters[[i]],1) } } # Assigns a standard deviation matrix type for each population pop$V21<-pop$V15 pop$V21<-round(pop$V21, digits=0) # Creates an interim replicate simulation file called pref_*.mp that will be used to vary the other input parameters # Writes first 32 lines of the input file (pref_*.mp): write(InputFile[1], file=prefile, append=FALSE) for (i in 2:44) { write(InputFile[i], file=prefile, append=TRUE) } # Writes model of density dependence lines # write("Yes", file=prefile, append=TRUE) # write(Selected_DD, file=prefile, append=TRUE) # Writes lines 35-44 lines of the input file (pref_*.mp): # for (i in 35:44) # { # write(InputFile[i], file=prefile, append=TRUE) # } # Writes in the new population dataframe write.table(pop,file=prefile,append=TRUE,sep=",",row.names=FALSE,col.names=FALSE,na="", quote=FALSE) # Writes "Migration" on a line write("Migration", file = prefile, append = TRUE) # Writes the line after "Migration" write("FALSE", file = prefile, append = TRUE) # Writes dispersal distance function parameters write(InputFile[Migration+2], file = prefile, append = TRUE) # Writes Dispersal Matrix Extra_NA<-c(rep(NA, length(Dispersal_matrix[1,]))) write.table((cbind(Dispersal_matrix, Extra_NA, deparse.level=0)), file = prefile, sep = ",", append = TRUE, quote=FALSE, row.names=FALSE, col.names=FALSE, na="") # Writes "Correlation" write("Correlation", file = prefile, append = TRUE) # Writes the line after "Correlation" write("FALSE", file = prefile, append = TRUE) # Writes Correlation_distance function parameters write(InputFile[Correlation+2], file = prefile, append = TRUE) # Writes Correlation_Matrix diag(Correlation_matrix)<-c(rep(1, length(Correlation_matrix[1,]))) Correlation_matrix<-cbind(Correlation_matrix, Extra_NA, deparse.level=0) write.table(Correlation_matrix, file = prefile, sep = ",", append = TRUE, quote=FALSE, row.names=FALSE, col.names=FALSE, na="") # Writes the stage and standard deviation information for (i in Line_N_stage_matrices:(Line_Constraints-1)) { write(InputFile[i], file = prefile, append = TRUE) } # Writes in the Constraints Matrix, relative dispersal indices and catastrophe multipliers for (i in Line_Constraints:(Line_Constraints+(N_stages[1]*3+3))) { write(InputFile[i], file = prefile, append = TRUE) } # Writes Initial_Abundance matrix write.table(Initial_Abundances, prefile, append = TRUE, col.names = FALSE, row.names=FALSE) # Reads and writes in the information from the "stages menu" for (i in (Line_Constraints+1+(N_stages[1]*3+3)+Npops):(Line_Constraints+(N_stages[1]*3+3)+Npops+5*N_stages[1])) { write(InputFile[i], file = prefile, append = TRUE) } # Ensures there is no population management modeled write("0 (pop mgmnt)", file = prefile, append = TRUE) #Note: some users may be interested in knowing how sensitive model predictions are to different levels of management # Writes a generic extinction threshold of zero write("0",file = prefile, append = TRUE) # Writes a generic explosion threshold of zero write("0",file = prefile, append = TRUE) # Writes the timestep datum as specified in the original input file Mgmnt <- grep("mgmnt", InputFile) N_Mgmnt <- scan(pref_metapop, what = "list", skip = (Mgmnt-1), nlines = 1) N_Mgmnt <- as.numeric(N_Mgmnt)[1] write(InputFile[Mgmnt+N_Mgmnt[1]+1+2], file = prefile, append = TRUE) # Writes end-of-file write ("-End of file-", file = prefile, append = TRUE) gc() ###--------------------- # # Re-initializes the interim replicate simulation input file (i.e. prefile) # ###--------------------- #Assigns the prefile name to the object New_Input New_Input<-prefile #Code to read in the whole prefile as an object InputFile<-readLines(New_Input, n=-1) #Gives the line number in input file that has "Migration" on it Migration<-grep("Migration", InputFile) #Gives the line number in input file that has "Correlation" on it Correlation<-grep("Correlation", InputFile) #Gives the number of populations, including the dispersal mortality population, in the input file Npops<-(Migration-41)-4 #Reads in the population dataframe pop<-read.table(New_Input,skip=44,sep=",",nrow=Npops) #Identifies the line number that the term "Constraints Matrix" is located on within the original input file Line_Constraints<-grep("Constraints Matrix",InputFile) #Reads in the number of stages (N_stages) as a numeric object # this code is slightly different from the line 10 scan specified above as it coerces it into a numeric object line10_numeric <- as.numeric(scan(New_Input,what="numeric",skip=9,nlines=1)) N_stages <- line10_numeric[1] #Reads in the model description Description <- readLines(New_Input, n=6) #Reads in line 12 (regional probability of catastrophe 1) line12 <- as.numeric(scan(New_Input,what="numeric",skip=11,nlines=1,sep=",")) #Reads in line 13 (spatial extent of catastrophe 1) line13 <- scan(New_Input,what="character",skip=12,nlines=1) #Reads in line14 (parameters affected by catastrophe 1) line14<-scan(New_Input,skip=13,nlines=1,sep=",",what="character") #Reads in line 19 (probability of REGIONAL catastrophe 2 occurrence) line19<-as.numeric(scan(New_Input,what="numeric",skip=18,nlines=1,sep=",")) #Reads in line 20 (spatial extent of catastrophe 2) line20 <- scan(New_Input,what="character",skip=19,nlines=1) #Reads in line 21 (parameters affected by catastrophe 2) line21 <- scan(New_Input,skip=20,nlines=1,sep=",",what="character") #Varies probability and intensity of catastrophe 1 #Randomly varies spatial extent of catastrophe 1 ext<-c("Regional","Local") line13_new<-sample(ext,1) ext_1<-line13_new prob1<-0 vect<-c(seq(1, (Npops-1))) #Determines what probability of catastrophe 1 to vary by 10% CV, if Catastrophe 1 occurs if(length(levels(pop$V13))==0) { if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"| line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { if(line13=="Regional") { prob1<-line12 } else { prob1<-mean(lpc) } #Samples probability of catastrophe 1 (assumes a normal distribution) prob1<-rnorm(1,mean=prob1,sd=0.1*prob1) #Ensures prob1 is constrained to between 0.001 and 1 if(prob1>1) { prob1<-1 } if(prob1<0.001) { prob1<-0.001 } } # Assigns the newly sampled probability of occurrence of catastrophe 1 to each population # Note: if spatial extent of catastrophe is local, all populations have the same prob1 if (line13_new=="Regional") { line12<-prob1 line12<-round(line12, digits=3) } if (line13_new=="Local") { # replaces new probability of occurrence values into the appropriate column in the population dataframe table pop$V13<-c(rep(prob1,times=length(pop[,13]))) pop$V13<-round(pop$V13, digits=3) } } # Varies the local multiplier for catastrophe 1 in the population dataframe, if there are no temporal trend files specified #Only multipliers <1 are varied, code is set so that if the original model indicates a local multiplier # or a stage-multiplier of 1 (i.e. no effect of catastrophe) this constraint is observed # If original model indicates a local or stage-specific multiplier >1, this constraint is observed # If both catastrophes 1 and 2 are modeled, the local and stage-specific multipliers (and probability of occurrences) # are sampled separately if(length(levels(pop$V12))==0) { #Calculates the mean local multiplier for catastrophe 1 mult1<-mean(lpm) if(mult1<1) { #Samples local multiplier of catastrophe 1 (assumes a normal distribution) mult1<-rnorm(1,mean=mult1,sd=0.1*mult1) #Code to ensure that the new mult1 is constrained to between 0 and 0.999999. if(mult1>0.999999) { mult1<-0.999999 } if(mult1<0) { mult1<-0 } } #Note: the multiplier is varied only if the average of the local population multipliers is <1. # Code that identifies if a parameter that catastrophe 1 affects has been specified and varies the local # population multiplier if the mean of all lpm is <1 (multiplier is sampled from a normal distribution) if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"| line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { # replaces mult1 into the appropriate column in the population dataframe table pop$V12<-c(rep(mult1,times=length(pop[,12]))) pop$V12<-round(pop$V12, digits=3) } } # If catastrophe 1 affects abundances, varies the Stage-specific multiplier # Scans in Abundances population-specific multipliers for catastrophe 1 Cat1mult_Abun<-scan(New_Input,skip=Line_Constraints+(2*N_stages)+1,nlines=1) Cat1mult_Abun_new<-Cat1mult_Abun if (any(line14=="Abundances"|line14=="abundances")) { # The new multiplier is sampled from a uniform distribution and this value is replaced in each element <1 # Elements >=1 remain the same Cat1multA<-runif(1,min=0,max=0.999999) for(i in 1:length(Cat1mult_Abun)) { if (Cat1mult_Abun[i] <1) { Cat1mult_Abun_new[i]<-c(rep((Cat1multA),times=length(Cat1mult_Abun[i]))) } } } Cat1mult_Abun_new<-round(Cat1mult_Abun_new, digits = 6) # If catastrophe 1 affects vital rates, varies the Stage-specific multiplier matrix #Scans in and reads Vital Rates stage specific multiplier for Catastrophe 1 Cat1mult_VR<-matrix(scan(New_Input, sep="", nlines = (N_stages), skip = Line_Constraints+N_stages+1), ncol = (N_stages), byrow = TRUE) Cat1mult_VR_new<-Cat1mult_VR if (any(line14=="Vital Rates"|line14=="vital rates")) { # The new multiplier is sampled from a uniform distribution and this value is replaced in each element <1 of the matrix # Elements >=1 remain the same Cat1multVR<-runif(1,min=0,max=0.999999) for(i in 1:length(Cat1mult_VR[,1])) { for(j in 1:length(Cat1mult_VR[1,])) { if (Cat1mult_VR[i,j] <1) { Cat1mult_VR_new[i,j] <- c(rep(Cat1multVR,times=length(Cat1mult_VR[i,j]))) } } } } # Ensures values are rounded to 6 decimal places Cat1mult_VR_new<-round(Cat1mult_VR_new, digits = 6) # Varies probability and intensity of catastrophe 2 # Randomly varies spatial extent of catastrophe 2 ext<-c("Regional","Local") line20_new<-sample(ext,1) ext_2<-line20_new prob2<-0 vect<-c(seq(1, (Npops-1))) # Determines what probability of catastrophe 2 to vary by 10% CV, if Catastrophe 2 occurs, if there are no Temporal Trend files associated with probability of cat 2 occurrence if(length(levels(pop$V20))==0) { if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { if(line20=="Regional") { prob2<-line19 } else { prob2<-mean(lpc2) } # Samples probability of catastrophe 2 (assumes a normal distribution) prob2<-rnorm(1,mean=prob2,sd=0.1*prob2) # Code to ensure that prob2 is constrained to between 0.001 and 1. if(prob2>1) { prob2<-1 } if(prob2<0.001) { prob2<-0.001 } } # Assigns the newly sampled probability of occurrence of catastrophe 1 to each population # Note: if spatial extent of catastrophe is local, all populations have the same prob1 if (line20_new=="Regional") { line19<-prob2 line19<-round(line19, digits = 3) } if (line20_new=="Local") { #Replaces new probability of occurrence values into the appropriate column in the population dataframe table pop$V20<-c(rep(prob2,times=length(pop[,20]))) pop$V20<-round(pop$V20, digits = 3) } } # Varies the local multiplier for catastrophe 2 in the population dataframe, if there are no temporal trend files specified #O nly multipliers <1 are varied, code is set so that if the original model indicates a local multiplier or a stage-multiplier of 1 (i.e. no effect of catastrophe) this constraint is observed # If original model indicates a local or stage-specific multiplier >1, this constraint is observed # If both catastrophes 1 and 2 are modeled, the local and stage-specific multipliers (and probability of occurrences) are sampled separately if(length(levels(pop$V19))==0) { # Calculates the mean local multiplier for catastrophe 2 mult2<-mean(lpm2) if(mult2<1) { # Samples local multiplier of catastrophe 2 (assumes a normal distribution) mult2<-rnorm(1,mean=mult2,sd=0.1*mult2) # Ensures that mult2 is constrained to between 0 and 0.999999. if(mult2>0.999999) { mult2<-0.999999 } if(mult2<0) { mult2<-0 } } # Note: mult2 is varied only if average of the local population multipliers is <1. # Code that identifies if a parameter that catastrophe 2 affects has been specified and varies the local population multiplier if the mean of all local population multipliers is <1 (multiplier is sampled from a normal distribution) if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { pop$V19<-c(rep(mult2,times=length(pop[,19]))) pop$V19<-round(pop$V19, digits=3) } } # If catastrophe 2 affects abundances, varies the stage-specific multiplier # scans in Abundances population specific multipliers for catastrophe 2 Cat2mult_Abun<-scan(New_Input,skip=Line_Constraints+(3*N_stages)+2,nlines=1) Cat2mult_Abun_new<-Cat2mult_Abun if (any( line21=="Abundances"|line21=="abundances")) { # The new multiplier is sampled from a uniform distribution and this value is replaced in each element <1 # Elements >=1 remain the same Cat2multA<-runif(1,min=0,max=0.999999) for(i in 1:length(Cat2mult_Abun)) { if (Cat2mult_Abun[i] <1) { Cat2mult_Abun_new[i]<-c(rep((Cat2multA),times=length(Cat2mult_Abun[i]))) } } } # Ensures values are rounded to 6 decimal places Cat2mult_Abun_new<-round(Cat2mult_Abun_new, digits = 6) # If catastrophe 2 affects vital rates, varies the stage-specific multiplier matrix # Reads Vital Rates stage specific multiplier for Catastrophe 2 Cat2mult_VR<-matrix(scan(New_Input, sep="", nlines = (N_stages), skip = Line_Constraints+2*N_stages+2), ncol = (N_stages), byrow = TRUE) Cat2mult_VR_new<-Cat2mult_VR if (any(line21=="Vital Rates"|line21=="vital rates" )) { # The new multiplier is sampled from a uniform distribution and this value is replaced in each element <1 of the matrix # Elements >=1 remain the same Cat2multVR<-runif(1,min=0,max=0.999999) for(i in 1:length(Cat2mult_VR[,1])) { for(j in 1:length(Cat2mult_VR[1,])) { if (Cat2mult_VR[i,j] <1) { Cat2mult_VR_new[i,j] <- c(rep(Cat2multVR,times=length(Cat2mult_VR[i,j]))) } } } } # Ensures values are rounded to 6 decimal places Cat2mult_VR_new<-round(Cat2mult_VR_new, digits = 6) # Randomly varies correlations between catastrophes 1 and 2, if both were modeled # Zero = no correlation # MaxNegative = negatively correlated to maximum extent possible # MaxPositive = positively correlated to maximum extent possible line25<-scan(New_Input,what="character",skip=24,nlines=1,sep=",") temporal_catas<-line25[1] corr_catas<-line25[2] cc<-c("Zero","MaxNegative","MaxPositive") corr_catas<-sample(cc,1) line25<-c(temporal_catas,corr_catas) # Identifies whether vital rates are sampled from normal or lognormal distributions EV_Distn_lognorm<-grep("Lognormal",InputFile) EV_Distn_norm<-grep("Normal", InputFile) if (length(EV_Distn_lognorm)>0) { Line_EV_Distn<-EV_Distn_lognorm } if (length(EV_Distn_norm)>0) { Line_EV_Distn<-EV_Distn_norm } EV_Distn<-scan(New_Input, what = "list", skip = (Line_EV_Distn[1]-1), nlines = 1, sep = ",") gc() # Varies K for each population # The K for each population were calculated in the RAMAS Spatial results (i.e. the replicate .ptc file). # The following code will sample new carrying capacity values # If K is not a function of habitat (or that relationship is unknown) and is therefore represented as a # constant (either a 0 or 1) in RAMAS Spatial, then the RAMAS Spatial results will indicate that the constant value is # that respective population's K. In this case, the K values will be sampled from the Metapop Link model. Under these # circumstances the user should have specified more than 1 population (with a separate K for each) in the Metapop Link file. # The new K values will be sampled using the mean and SD of K across all populations; a 10% CV is applied if there is SD = 0 # If K is a function of habitat, K values for each population were calculated by RAMAS Spatial based on this function and # new K values are sampled using the K as the mean and a 10% CV is applied # Reads in habitat specific carrying capacity function or constant value (line 16) k_link <- scan(ptcfile, what="numeric", skip=13,nlines=1) # Reads in the original number of populations in the Metapop Link model Original_Npops<-(Migration_linkmp-41)-4 # Reads in the population dataframe from the replicate .ptc file pop_link<-read.table(metapop_link,skip=44,sep=",",nrow=Original_Npops) # If K is not a function of habitat: calculates the mean and SD of K based on the link *.mp file; sets SD to 10% CV if SD = 0 # The objects Mean_K_link and SD_K_link are needed only if K is not a function of habitat, but objects are created regardless Mean_K_link<-mean(pop_link$V7, na.rm = TRUE) if (length(pop_link$V7)==1) { SD_K_link<-(Mean_K_link*0.1) } else { SD_K_link<-sd(pop_link$V7, na.rm = TRUE) if(SD_K_link==0) { SD_K_link<-(Mean_K_link*0.1) } } # NOTE: If no variability in carrying capacity, a 10% CV is assumed if (k_link==0|k_link==1) # if K is not a function of habitat { for (i in 1:Npops) { pop$V7[i]<-as.integer(rnorm(1, mean=Mean_K_link, sd=SD_K_link)) if(pop$V7[i]<0) { pop$V7[i]<-0 } # Sets constraints for carrying capacity } # Sets dispersal mortality new carrying capacity to 0 pop$V7[Npops]<-0 # Rounds initial abundances to 0 decimal places pop$V7<-round(pop$V7, digits=0) } else { K<-pop$V7[-Npops] #Creates a vector of K values stdK<-pop$V8[-Npops] #Creates a vector of SD of K K_new<-K #Creates a new empty vector with mode and length of vector 'K' #Samples a new K from a normal distribution with mean of K and CV =10% if no standard deviation was indicated for (i in 1:length(K_new)) { if (stdK[i]==0) #if SD of K is = 0 { K_new[i]<-rnorm(1,mean=K[i],sd=0.1*K[i]) if(K_new[i]<0) { K_new[i]<-0 } } else { K_new[i]<-rnorm(1,mean=K[i],sd=stdK[i]) #if SD of K is > 0 if(K_new[i]<0) { K_new[i]<-0 } } } # Adds new carrying capacities to population dataframe table pop$V7[-Npops]<-K_new # Ensures carrying capacities are rounded to 0 decimal places pop$V7<-round(pop$V7, digits=0) } # Calculates the number of populations with K = 0 N_sinks<-0 # Creates a new object for (i in 1:(Npops-1)) { if (pop$V7[i]==0) { N_sinks<-N_sinks+1 } } # Enforces constraint on the K of the Dispersal mortality population: no variability, and set to arbitrarily low number. pop$V7[Npops]<-2 pop$V8[Npops]<-0 #Some models may require a minimum K for all populations, including the sink population. #Randomizes Rmax if(SD_Rmax==0) { SD_Rmax<-Mean_Rmax*0.1 } Selected_Rmax<-rnorm(1,mean=Mean_Rmax,sd=SD_Rmax) if(Selected_Rmax<1) { Selected_Rmax<-1 } for(i in 1:Npops) { pop$V6[i]<-Selected_Rmax } # Varies dispersal rates and connection among populations # # Reads in the dispersal-distance parameters M_parameters<-scan(New_Input, skip = (Migration[1]+1), nlines = 1, sep=",") # Reads in the dispersal matrix Dispersal_matrix<-matrix(scan(New_Input, sep=",", nlines = (Correlation-Migration-3), skip = Migration+2), ncol = (Correlation-Migration-2), byrow = TRUE) Dispersal_matrix<-Dispersal_matrix[,-(length(Dispersal_matrix[1,]-1))] # Varies the number of connections among populations (i.e. population pairs connected through dispersal) in # the dispersal matrix (excluding possible connections with sink population), sampled from a uniform distribution N_elements<-runif(1,0,length(Dispersal_matrix[-length(Dispersal_matrix[1,]),-length(Dispersal_matrix[,1])])) N_elements<-round(N_elements) P_connection<-N_elements/length(Dispersal_matrix[-length(Dispersal_matrix[1,]),-length(Dispersal_matrix[,1])]) # Creates a matrix that randomly selects which elements will represent connections, and which will not. Logical_matrix<-matrix(, Npops[1], Npops[1]) for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { Logical_matrix[i,j]<-rbinom(1, 1, P_connection) } } # Converts the numeric matrix to a logical matrix where 'TRUE' signifies a connection Logical_matrix<-matrix(as.logical(Logical_matrix),Npops, Npops) # For population pairs not connected by dispersal in this replicate simulation, the corresponding element is set to 0 for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (any(Logical_matrix[i,j])) { Dispersal_matrix[i,j]<-Dispersal_matrix[i,j] } else { Dispersal_matrix[i,j]<-0 } } } # Applies a transformation (sampled from a normal distribution) to vary dispersal rates Transform_dispersal<-rnorm(1, 0, 0.1) for (i in 1:Npops) { for (j in 1:Npops) { Dispersal_matrix[i,j]<-(Dispersal_matrix[i,j]+(Dispersal_matrix[i,j]*Transform_dispersal)) } } # Code to sample dispersal rates from a random uniform distribution from 0 to 1, if no dispersal parameters were specified Unif_dispersal<-runif(1,0,1) if(sum(Dispersal_matrix)==0) { if(M_parameters[1]==0) { for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (any(Logical_matrix[i,j])) { Dispersal_matrix[i,j]<-Unif_dispersal } else { Dispersal_matrix[i,j]<-0 } } } } } # Ensures diagonal elements on the dispersal matrix are 0 diag(Dispersal_matrix)<-c(rep(0,length(Dispersal_matrix[1,]))) # Enforces a symmetrical matrix prior to adding the sink population, if distances among populations are calculated based # on center-center or edge-edge. If distance based on centre-edge, dispersal matrix remains assymetrical if ( ptcinput[27] == "Center to edge") { Dispersal_matrix[i,j]<-Dispersal_matrix[i,j] } else { Upper_logical<-upper.tri(Dispersal_matrix) for (i in 1:length(Dispersal_matrix[,1])) { for (j in 1:length(Dispersal_matrix[1,])) { if(any(Upper_logical[i,j])) { Dispersal_matrix[j,i]<-Dispersal_matrix[i,j] } } } #All dispersal rates are symmetric (same in both directions), except for sink population and if centre-edge is distance type } # Enforces constraints on individual elements (i.e. dispersal rate varies from 0 to 1) for (i in 1:length(Dispersal_matrix[,1])) { for (j in 1:length(Dispersal_matrix[1,])) { if (Dispersal_matrix[i,j]<0) { Dispersal_matrix[i,j]<-0 } if (Dispersal_matrix[i,j]>1) { Dispersal_matrix[i,j]<-1 } } } # Enforces a maximum total dispersal rate of 1 from each population Sum_Columns_Dispersal_matrix<-c(seq(1:Npops)) for (i in 1:Npops) { Sum_Columns_Dispersal_matrix[i]<-sum(Dispersal_matrix[,i]) } for(i in 1:Npops) { if (Sum_Columns_Dispersal_matrix[i]>1) { Dispersal_matrix[,i]<-Dispersal_matrix[,i]/Sum_Columns_Dispersal_matrix[i] } } # If total dispersal exceeds 100%, elements are reduced proportionally so that total dispersal is 1.00. # Varies dispersal mortality, sampled from a uniform distribution # Assumes same proportion of dispersers die regardless of source or target populations) Dispersal_surv<-runif(1, 0, 1) Total_mort<-matrix(0, Npops, Npops) for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { Dispersal_matrix[i,j]<-(Dispersal_matrix[i,j]*Dispersal_surv) Total_mort[i,j]<-(Dispersal_matrix[i,j]-(Dispersal_matrix[i,j]*Dispersal_surv)) } } # Adds all the dead dispersers from each population Disp_mort_sums<-apply(Total_mort, 2, sum, na.rm = TRUE) # Add dead dispersers to sink population Dispersal_matrix[length(Dispersal_matrix[,1]),]<-Disp_mort_sums #Keeps track of dispersal rates for calculating mean, standard deviation, minimum and maximum dispersal rates below Lower_logical<-lower.tri(Dispersal_matrix) Dispersal_rate<-Dispersal_matrix[-Npops,-Npops] for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if ( ptcinput[27] == "Center to edge") { Dispersal_rate[i,j]<-Dispersal_matrix[i,j] #if distance based on centre-edge, dispersal matrix remains assymetrical and a full matrix is kept (with no NAs) } else { if (any(Lower_logical[i,j])) { Dispersal_rate[i,j]<-Dispersal_matrix[i,j] } else { Dispersal_rate[i,j]<-NA # upper logical of matrix is coerced to NAs } } } } # Gives the line number in Spatial input file (*.ptc) that has "Landscape_indices" on it Landscape_indices<-grep("Landscape indices", ptcinput) # Creates an empty distance matrix Pairwise_Distance <- array(0, c(Npops-1,Npops-1)) # Reads in distance matrix; part of spatial .ptc results # These distances are based on the type of distance measure specified, ie centre-centre, centre-edge, edge-edge if ( ptcinput[27] == "Center to edge") { # To check Pairwise_Distance <- matrix(scan(ptcfile, nlines = (Npops-1), skip = (Landscape_indices-Npops)),ncol = (Npops-1), byrow = TRUE) } else { # Scans in pairwise distance matrix (only lower half exists) from varied Spatial Data file, which includes results #To delete Pairwise_Distance[upper.tri(Pairwise_Distance, diag = TRUE)] <- scan(ptcfile, nlines = Npops, skip = (Landscape_indices-Npops-1)) # modified (feb 2016 Pairwise_Distance[upper.tri(Pairwise_Distance, diag = TRUE)] <- scan(ptcfile, nlines = (Npops-1), skip = (Landscape_indices-Npops)) # Creates a full matrix of pairwise distances among populations, where the lower and upper halves of the matrix are symmetrical Pairwise_Distance <- Pairwise_Distance + t(Pairwise_Distance) } # Adds Dispersal mortality population rows and columns to distance matrix Pairwise_Distance<-rbind(Pairwise_Distance, matrix(0, 1, length(Pairwise_Distance[1,]))) Pairwise_Distance<-cbind(Pairwise_Distance, matrix(0, length(Pairwise_Distance[,1]),1)) # Keeps track of the distances among population pairs that are connected through dispersal Distances<-Dispersal_matrix[-Npops,-Npops] for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (ptcinput[27] == "Center to edge") { if(Dispersal_matrix[i,j]>0) { Distances[i,j]<-Pairwise_Distance[i,j] } else { Distances[i,j]<-0 } } else { if(any(Lower_logical[i,j])) { if(Dispersal_matrix[i,j]>0) { Distances[i,j]<-Pairwise_Distance[i,j] } else { Distances[i,j]<-0 } } else { Distances[i,j]<-NA #sets upper part of matrix to NAs } } } } # Ensures values in the dispersal matrix are rounded to 7 decimal places Dispersal_matrix<-round(Dispersal_matrix,7) # Add the extra column of "," at the end of each matrix row for writing to the replicate simulation input file Extra_NA<-c(rep(NA, length(Dispersal_matrix[1,]))) Dispersal_matrix<-cbind(Dispersal_matrix, Extra_NA, deparse.level=0) gc() # Varies coefficients of correlations in vital rates among population pairs # NOTE: Correlations among populations are assumed to range from 0 to 1, for simplicity. # Reads in the correlation-distance function parameters C_parameters<-scan(New_Input, skip = (Correlation[1]+1), nlines = 1, sep=",") # Reads in the correlation matrix Correlation_matrix<-matrix(scan(New_Input, sep = ",", nlines = Npops, skip =Correlation+2), nrow = Npops, byrow = TRUE) Correlation_matrix<-Correlation_matrix[,-(length(Correlation_matrix[1,]-1))] # Create a matrix of pairwise distances between geometric population centres Centre_Pairwise_Distance<-matrix(0,Npops,Npops) for (i in 1:Npops) { for (j in 1:Npops) { Centre_Pairwise_Distance[i,j]<-sqrt((pop[i,2]-pop[j,2])^2+(pop[i,3]-pop[j,3])^2) } } # Adjusts the correlations according to the distances among population pairs (centre-centre of populations) and for cases where autofill option was used for (i in 1:(Npops-1)) { for(j in 1:(Npops-1)) { Correlation_matrix[i,j]<-(C_parameters[1]*exp((-Centre_Pairwise_Distance[i,j]^C_parameters[3])/C_parameters[2])) } } # Applies a transformation (sampled from a normal distribution) to vary the magnitude of the correlations Transform_correlation<-rnorm(1,0,0.1) for (i in 1:Npops) { for (j in 1:Npops) { Correlation_matrix[i,j]<-(Correlation_matrix[i,j]+(Correlation_matrix[i,j]*Transform_correlation)) } } #Ensures diagonal elements on the correlation matrix are 1. diag(Correlation_matrix)<-c(rep(1,length(Correlation_matrix[1,]))) #Ensures a symmetrical matrix For_C_logical<-upper.tri(Correlation_matrix) for (i in 1:length(Correlation_matrix[,1])) { for (j in 1:length(Correlation_matrix[1,])) { if(any(For_C_logical[i,j])) { Correlation_matrix[i,j]<-Correlation_matrix[j,i] } } } #Enforces constraints on individual elements (elements range from 0 to 1) for (i in 1:length(Correlation_matrix[,1])) { for (j in 1:length(Correlation_matrix[1,])) { if (Correlation_matrix[i,j]<0) { Correlation_matrix[i,j]<-0 } if (Correlation_matrix[i,j]>1) { Correlation_matrix[i,j]<-1 } } } # Ensures correlation coefficients are rounded to 6 decimal places Correlation_matrix<-round(Correlation_matrix, digits = 6) # Models correlations as a constant value (no distance function) for input files where a constant or function was not specified. Cte<-runif(1,0,1) if(sum(Correlation_matrix)==length(Correlation_matrix[,1])) { if(sum(C_parameters)<=0.1) { for(i in 1: length(Correlation_matrix[,1])) { for(j in 1:length(Correlation_matrix[,1])) { if(any(Lower_logical[i,j])) { Correlation_matrix[i,j]<-Cte } } } } } #Code that will keep track of correlations for calculating mean, standard deviation, minimum and maximum correlations below Correlations<-Correlation_matrix[-Npops,-Npops] for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (any(Lower_logical[i,j])) { Correlations[i,j]<-Correlation_matrix[i,j] } else { Correlations[i,j]<-NA } } } #Code to add the extra column of "," at the end of each matrix row for writing input file Extra_NA_C<-c(rep(NA, length(Correlation_matrix[1,]))) Correlation_matrix<-cbind(Correlation_matrix, Extra_NA_C, deparse.level=0) #Varies stage and standard deviation matrices # Reads in the number of stage matrices Line_N_stage_matrices<-grep("stage matrix", InputFile) Number_stage_matrix<-scan(New_Input, what="list", skip=Line_N_stage_matrices-1, nlines=1) Number_matrices<-as(Number_stage_matrix[1], "numeric") # Reads in the stage matrix names and three associated lines of code for each of the stage matrices and saves them as a list Stage_matrix_name<-list(length=Number_matrices) for (i in 1:Number_matrices) { Stage_matrix_name[[i]]<-scan(New_Input, what="list", skip=(Line_N_stage_matrices+((i-1)*(4+N_stages))), nlines=4, sep=",") } # Reads in the number of standard deviation matrices Line_N_stdev_matrix<-grep("st.dev.", InputFile) Number_stdev_matrix<-scan(New_Input, what="list", skip=Line_N_stdev_matrix-1, nlines=1) Number_stdev_matrices<-as(Number_stdev_matrix[1], "numeric") # Reads in the names of each of the standard deviation matrix Stdev_matrix_name<-list(length=Number_stdev_matrices) for (i in 1:Number_stdev_matrices) { Stdev_matrix_name[[i]]<-scan(New_Input, what="list", skip=(Line_N_stdev_matrix+((i-1)*(1+N_stages))), nlines=1, sep=",") } # Reads in the stage matrices Stage_matrices<-list(length=Number_matrices) for (i in 1:Number_matrices) { Stage_matrices[[i]]<-matrix(scan(New_Input, skip =(Line_N_stage_matrices+4+(i-1)*(4+N_stages)), nlines=N_stages), byrow=TRUE, ncol=N_stages) } # Reads the corresponding standard deviation matrices Stdev_matrices<-list(length=Number_stdev_matrices) for (i in 1:Number_stdev_matrices) { Stdev_matrices[[i]]<-matrix(scan(New_Input, skip =(Line_N_stdev_matrix+1+(i-1)*(1+N_stages)), nlines=N_stages), byrow=TRUE, ncol=N_stages) } #Varies the stage matrix elements # Samples matrix element assuming lognormal distribution # As described in Burgman et al. (1993) Risk Assessment in Conservation Biology, Chapman and Hall, London, UK. if(EV_Distn[1]=="Lognormal") { lnorm_mean<-list() for(i in 1:Number_matrices) { lnorm_mean[[i]]<-matrix(0, N_stages, N_stages) } # Creates intermediate matrices of the value c = standard deviation/mean for(k in 1:Number_matrices) { for(i in 1:N_stages) { for(j in 1:N_stages) { lnorm_mean[[k]][i,j]<-Stdev_matrices[[k]][i,j]/Stage_matrices[[k]][i,j] } } } # Creates the new stage matrices for(k in 1:Number_matrices) { for(i in 1:N_stages) { for(j in 1:N_stages) { if(Stage_matrices[[k]][i,j]>0) { if(Stdev_matrices[[k]][i,j]>0) { Stage_matrices[[k]][i,j]<- exp(rnorm(1,mean=(log(Stage_matrices[[k]][i,j])-(0.5*log((lnorm_mean[[k]][i,j]^2)+1))), sd=sqrt(log((lnorm_mean[[k]][i,j]^2)+1)))) } else { Stage_matrices[[k]][i,j]<- exp(rnorm(1, mean=(log(Stage_matrices[[k]][i,j])-(0.5*log((0.1^2)+1))), sd= sqrt(log((0.1^2)+1)))) } } } } } } # Samples matrix elements assuming normal distribution if(EV_Distn[1]=="Normal") { for(k in 1:Number_matrices) { for(i in 1:N_stages) { for(j in 1:N_stages) { if(Stage_matrices[[k]][i,j]>0) { if(Stdev_matrices[[k]][i,j]>0) { Stage_matrices[[k]][i,j]<- rnorm(1,Stage_matrices[[k]][i,j], sd=Stdev_matrices[[k]][i,j]) } else { Stage_matrices[[k]][i,j]<- rnorm(1,Stage_matrices[[k]][i,j], sd=(Stage_matrices[[k]][i,j]*0.1)) } } } } } } # Enforces a constraint of 1 on survival rates # User-defined code will need to be added for PVAs that model both sexes separately. In these models, 2 lines of the matrix refer to fecundity values, which should not be constrained to <1.0. Please contact authors for sample code. if(N_stages[1]>1) { for(k in 1:Number_matrices) { for(i in 1:(N_stages-1)) { for(j in 1:N_stages) { if(Stage_matrices[[k]][i+1,j]>1) { (Stage_matrices[[k]][i+1,j]=0.99) } } } } } # Ensures matrix elements are rounded to 4 decimal places for(i in 1:Number_matrices) { Stage_matrices[[i]]<-round(Stage_matrices[[i]], digits = 4) } # Line_Constraints<-grep("Constraints", InputFile) #255 # Line_Constraints<-grep("Constraints Matrix",InputFile) # Reads in the constraints matrix Constraints <- matrix(scan(New_Input, skip=Line_Constraints, nlines=N_stages), N_stages, byrow=TRUE) Constraints<-round(Constraints, digits = 6) # Reads in the relative dispersal of stages Rel_Dispersal<-scan(New_Input, skip=Line_Constraints+N_stages, nlines=1, sep=" ") Rel_Dispersal<-round(Rel_Dispersal, digits = 6) # Reads in the initial abundances Initial_Abundances<-matrix(scan(New_Input, skip=(Line_Constraints+(N_stages[1]*3+3)), nlines=Npops), ncol=N_stages[1], byrow=TRUE) #Reads in the information from the "stages menu" Stages_data<-scan(New_Input, what="list", skip=145, nlines=35) #-------------------------------------# # Writes the new replicate simulation # # input file (rep_*.mp) # #-------------------------------------# #Writes model description in the first six lines of input file write(Description, file = filename, append = FALSE) #Writes the number of stochastic runs within each replicate simulation write(Nruns, file=filename, append=TRUE) #Writes number of time steps write("100",file=filename,append=TRUE) #Writes whether or not demographic stochasticity is assumed write(InputFile[9],file=filename,append=TRUE) #Writes number of stages and whether or not constraints are ignored write(InputFile[10],file=filename,append=TRUE,sep=",",ncolumns=2) #Writes the name of catastrophe 1 write(InputFile[11],file=filename,append=TRUE) #Writes the probability of catastrophe 1, if regional in scale write(line12, file=filename, append=TRUE) #Writes whether catastrophe 1 is local or regional write(line13_new,file=filename,append=TRUE) #Writes the parameters that catastrophe 1 affects write(InputFile[14], file=filename,append=TRUE) #Writes whether or not catastrophe 1 is spreading, and how write(InputFile[15], file=filename, append=TRUE) #Writes the probability of infection per disperser write(InputFile[16],file=filename,append=TRUE) #Writes the probability distance function if catastrophe 1 is 'spread by a vector' write(InputFile[17],file=filename,append=TRUE,ncolumns=4,sep=",") #Writes the name of catastrophe 2 write(InputFile[18],file=filename,append=TRUE) #Writes the probability of catastrophe 2, if regional in scale write(line19,file=filename,append=TRUE) #Writes whether catastrophe 2 is local or regional write(line20_new,file=filename,append=TRUE) #Writes the parameters that catastrophe 2 affects write(InputFile[21], file=filename,append=TRUE) #Writes whether or not catastrophe 2 is spreading, and how write(InputFile[22], file=filename, append=TRUE) #Writes the probability of infection per disperser write(InputFile[23], file=filename,append=TRUE) #Writes the probability distance function if catastrophe 2 is 'spread by a vector' write(InputFile[24], file=filename,append=TRUE) #Writes the temporal trend in vital rate after catastrophe, and the correlation between catastrophe 1 and 2 write(line25,file=filename,append=TRUE,sep=",",ncolumns=length(line25)) #Writes the parameters that density dependence affects write(InputFile[26],file=filename,append=TRUE) #Writes the assumed distribution of environment stochasticity in vital rates, and advanced settings write(InputFile[27], file=filename, append = TRUE) #Writes the CV in dispersal write(InputFile[28], file=filename, append = TRUE) #Writes what happens when population size falls below local threshold write("count in total",file=filename,append=TRUE) #Writes advanced settings for within-population correlation write(InputFile[30],file=filename,append=TRUE) #Writes whether or not dispersal depends on target population carrying capacity write(InputFile[31],file=filename,append=TRUE) #Writes advanced settings for density-dependence write(InputFile[32],file=filename,append=TRUE) #Writes whether or not density-dependence is population-specific or not write(InputFile[33], file=filename, append=TRUE) #Writes the type of density-dependence write(InputFile[34],file=filename,append=TRUE) #Writes the location of a file with a user-defined density-dependence function, if specified write(InputFile[35], file=filename,append=TRUE) #Writes the duration in units of a single time-step write(InputFile[36],file=filename,append=TRUE) #Writes the unit of the time step write(InputFile[37],file=filename,append=TRUE) #Writes the sex structure of the model write(InputFile[38],file=filename,append=TRUE) #Writes the number of stages for females write(InputFile[39],file=filename,append=TRUE) #Writes the mating system write(InputFile[40],file=filename,append=TRUE) #Writes the number of females a male can mate with write(InputFile[41],file=filename,append=TRUE) #Writes the number of males a female can mate with write(InputFile[42],file=filename,append=TRUE) #Writes the CV for sampling of population size write(InputFile[43],file=filename,append=TRUE) #Writes line 44 (related to dispersal?) write(InputFile[44],file=filename,append=TRUE) #Writes the population dataframe write.table(pop,file=filename,append=TRUE,sep=",",row.names=FALSE,col.names=FALSE,na="",quote=FALSE) #Writes "Migration" on a line write("Migration", file = filename, append = TRUE) #Writes FALSE for use of the dispersal-distance function write("FALSE", file = filename, append = TRUE) #Writes the dispersal-distance function parameters write(M_parameters, file = filename, sep = ",", append = TRUE) #Writes the dispersal rate matrix write.table(Dispersal_matrix, file=filename, append = TRUE, quote = FALSE, sep = ",", na="", row.names=FALSE, col.names=FALSE) #Writes "Correlation" on a line write("Correlation", file = filename, append = TRUE) #Writes FALSE for use of the correlation-distance function write("FALSE", file = filename, append = TRUE) #Writes the correlation-distance function parameters write(C_parameters, file = filename, sep = ",", append = TRUE) #Writes the correlation matrix write.table(Correlation_matrix, file=filename, append = TRUE, quote = FALSE, sep = ",", na="", row.names=FALSE, col.names=FALSE) #Writes the line that specifies number of stage matrices write(Number_stage_matrix, file = filename, ncolumns = 5, append = TRUE, sep = " ") #Writes the stage matrix names, the three lines of associated code, and the corresponding stage matrices for (i in 1:Number_matrices) { write(Stage_matrix_name[[i]], file = filename, append = TRUE) write.table(Stage_matrices[[i]], file = filename, row.names = FALSE, col.names = FALSE, append = TRUE, sep = " ") } #Writes the line that specifies the number of standard deviation matrices write(Number_stdev_matrix, file = filename, ncolumns = 5, append = TRUE, sep = " ") #Writes the standard deviation matrix names, the lines of associated code, and the corresponding matrices for (i in 1:Number_matrices) { write(Stdev_matrix_name[[i]], file = filename, append = TRUE) write.table(Stdev_matrices[[i]], file = filename, row.names = FALSE, col.names = FALSE, append = TRUE, sep = " ") } #Writes the line that says "Constraints Matrix write("Constraints Matrix", file = filename, append = TRUE) #Writes the constraints matrix write.table(Constraints, file = filename, sep = " ", append = TRUE, row.names=FALSE, col.names=FALSE) #Writes the line of stage-specific relative dispersals write.table(t(Rel_Dispersal), file=filename, append = TRUE, quote = FALSE, sep = " ", na="", row.names=FALSE, col.names=FALSE) #Writes the stage-specific multipliers for vital rates for catastrophe 1 write.table(Cat1mult_VR_new, file=filename, append = TRUE, quote = FALSE, sep = " ", na="", row.names=FALSE, col.names=FALSE) #Writes the stage-specific multipliers for abundances for catastrophe 1 write(Cat1mult_Abun_new,file=filename,append=TRUE,ncolumns=length(Cat1mult_Abun),sep=" ") #Writes the stage-specific multipliers for vital rates for catastrophe 2 write.table(Cat2mult_VR_new, file=filename, append = TRUE, quote = FALSE, sep = " ", na="", row.names=FALSE, col.names=FALSE) #Writes the stage-specific multipliers for abundances for catastrophe 2 write(Cat2mult_Abun_new,file=filename,append=TRUE,ncolumns=length(Cat1mult_Abun),sep=" ") #Write the stage-specific initial abundances matrix write.table(Initial_Abundances, file=filename, append = TRUE, col.names = FALSE, row.names=FALSE) #Reads and writes the information from the stages menu for (i in (Line_Constraints+1+(N_stages[1]*3+3)+Npops):(Line_Constraints+(N_stages[1]*3+3)+Npops+5*N_stages[1])) { write(InputFile[i], file =filename, append = TRUE) } #Writes that there is no population management write("0 (pop mgmnt)", file = filename, append = TRUE) #Writes the extinction and explosion thresholds and timestep datum #Writes a generic extinction threshold of zero write("0",file = filename, append = TRUE) #Writes a generic explosion threshold of zero write("0",file = filename, append = TRUE) #Writes a generic timestep datum equal to 3.000 (as required for 100 timesteps) write("3", file = filename, append = TRUE) #Writes end-of-file write ("-End of file-", file = filename, append = TRUE) #---END create rep_x.mp file-----------# #--------------------------------------# #--------------------------------------# # Writes values of input parameters # # for further analysis # #--------------------------------------# # Creates empty vectors to store # fecundities, survivals and # other input parameters variable<-as.data.frame(matrix(,1,42)) if(N_stages>1) { Vital_rates<-as.data.frame(matrix(,1,(N_stages*2))) } else { Vital_rates<-0 } # Calculates 42 variables of potential interest # Reads the number of populations included in replicate simulation variable[1]<-Npops-1 # Calculates the mean population-specific carrying capacity variable[2]<-mean(pop$V7[-Npops]) # Calculates the standard deviation of the population-specific carrying capacities variable[3]<-sd(pop$V7[-Npops]) # Reads the maximum population-specific carrying capacity variable[4]<-max(pop$V7[-Npops]) # Reads the minimum population-specific carrying capacity variable[5]<-min(pop$V7[-Npops]) # Calculates the mean geographic distances (centre-centre) among population pairs, excluding the sink population Centre_Pairwise_Distances<-Centre_Pairwise_Distance[-Npops,-Npops] variable[6]<-sum(Centre_Pairwise_Distances)/((Npops-1)*(Npops-1)-(Npops-1)) # Calculates the maximum distances among population pairs, excluding the sink population variable[7]<-max(Centre_Pairwise_Distances) # Calculates the minimum distances among population pairs, excluding the sink population for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (Centre_Pairwise_Distances[i,j]==0) { Centre_Pairwise_Distances[i,j]<-500000000 } } } variable[8]<-min(Centre_Pairwise_Distances) # Calculates the mean population-specific initial abundance variable[9]<-mean(pop$V4[-Npops]) # Calculates the standard deviation of the population-specific initial abundances variable[10]<-sd(pop$V4[-Npops]) # Reads the maximum population-specific initial abundance variable[11]<-max(pop$V4[-Npops]) # Reads the minimum population-specific initial abundance variable[12]<-min(pop$V4[-Npops]) # Reads the spatial extent of catastrophe 1 (writes NA if no catastrophe) if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"| line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { variable[13]<-ext_1 } else { variable[13]<-"NA" } # Reads the spatial extent of catastrophe 2 (writes NA if no catastrophe) if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { variable[14]<-ext_2 } else { variable[14]<-"NA" } # Reads the nature of the correlation between catastrophes 1 and 2 (writes NA if <2 catastrophes are modeled) if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"| line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { variable[15]<-corr_catas } } else { variable[15]<-"NA" } # Reads the probability of catastrophe 1 (writes NA if no catastrophe) if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"| line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { variable[16]<-prob1 } else { variable[16]<-"NA" } #Reads the population multiplier for catastrophe 1 (write NA if no catastrophe) if(Levels_12==0) { if (any(line14=="Dispersal Rates" | line14=="Carrying Capacities"| line14=="Vital Rates"| line14=="Abundances"|line14=="dispersal rates" | line14=="carrying capacities"| line14=="vital rates"|line14=="abundances")) { variable[17]<-mult1 } else { variable[17]<-"NA" } } #Reads the stage-specific multiplier if catastrophe 1 affects abundances (write NA if no catastrophe) if (any(line14=="Abundances"|line14=="abundances")) { variable[18]<-Cat1multA } else { variable[18]<-"NA" } #Reads the stage-specific multiplier if catastrophe 1 affects vital rates (write NA if no catastrophe) if (any(line14=="Vital Rates"|line14=="vital rates")) { variable[19]<-Cat1multVR } else { variable[19]<-"NA" } #Reads the probability of catastrophe 2 (writes NA if no catastrophe) if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { variable[20]<-prob2 } else { variable[20]<-"NA" } #Reads the population multiplier for catastrophe 2 (writes NA if no catastrophe) if (any(line21=="Dispersal Rates" | line21=="Carrying Capacities"| line21=="Vital Rates"| line21=="Abundances"| line21=="dispersal rates" | line21=="carrying capacities"| line21=="vital rates"|line21=="abundances")) { variable[21]<-mult2 } else { variable[21]<-"NA" } #Reads the stage-specific multiplier if catastrophe 2 affects abundances (write NA if no catastrophe) if (any( line21=="Abundances"|line21=="abundances")) { variable[22]<-Cat2multA } else { variable[22]<-"NA" } #Reads the stage-specific multiplier if catastrophe 2 affects vital rates (write NA if no catastrophe) if (any(line21=="Vital Rates"|line21=="vital rates" )) { variable[23]<-Cat2multVR } else { variable[23]<-"NA" } #Reads number by which dispersal rates were multiplied variable[24]<-Transform_dispersal #Calculates the number of populations connected by dispersal if distances were calculated using the centre-to-edge (2) or centre-to-centre distances (1) Connections<-0 for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (Dispersal_matrix[i,j]>0) { Connections<-Connections+1 } } } variable[25]<-Connections #Reads in the dispersal mortality variable[26]<-Dispersal_surv #Reads the number by which correlations were multiplied variable[27]<-Transform_correlation #Calculates the minimum convex polygon area occupied by the populations, excluding sink population if((Npops-1)>4) #MCP can only be calculated if 5 or more populations were modeled { variable[28]<-"NA" } else { variable[28]<-"NA" } #Calculates the number of sink populations (i.e. with carrying capacity or ceiling of zero) variable[29]<-N_sinks #coerces diagonal of matrix to NAs if (ptcinput[27] == "Center to edge") (2), if Center to center (1) or edge to edge (3), this is already NA diag(Dispersal_rate)<-c(rep(NA, length(Dispersal_rate[1,]))) #Calculates the mean rate of dispersal among all population pairs variable[30]<-mean(Dispersal_rate, na.rm=TRUE) #Calculates the minimum rate of dispersal among all population pairs variable[31]<-min(Dispersal_rate, na.rm=TRUE) #Calculates the maximum rate of dispersal among all population pairs variable[32]<-max(Dispersal_rate, na.rm=TRUE) #Calculates the mean correlation coefficients among all population pairs variable[33]<-mean(Correlations, na.rm=TRUE) #Calculates the minimum correlation coefficients among all population pairs variable[34]<-min(Correlations, na.rm=TRUE) # Calculates the maximum correlation coefficients among all population pairs variable[35]<-max(Correlations, na.rm=TRUE) # Coerces diagonal of matrix to NAs # object 'Distances': Keeps track of the distances among population pairs that are connected through dispersal diag(Distances)<-c(rep(NA, length(Distances[1,]))) Distances2<-Distances if (Connections>0) { Distances2[Distances2==0]<-NA # Calculates the mean geographic distances among population pairs connected by dispersal variable[36]<-mean(Distances2, na.rm=TRUE) # Calculates the minimum geographic distances among population pairs connected by dispersal variable[37]<-min(Distances2, na.rm=TRUE) # Calculates the maximum geographic distances among population pairs connected by dispersal variable[38]<-max(Distances2, na.rm=TRUE) } else { # If there are no connections then we do not want to coerce any 0s to NAs # These 3 variables will all be 0 in any case, as there will be no connected populations # Calculates the mean geographic distances among population pairs connected by dispersal variable[36]<-mean(Distances2, na.rm=TRUE) # Calculates the minimum geographic distances among population pairs connected by dispersal variable[37]<-min(Distances2, na.rm=TRUE) # Calculates the maximum geographic distances among population pairs connected by dispersal variable[38]<-max(Distances2, na.rm=TRUE) } # Records Rmax variable[39]<-Selected_Rmax # Calculates the mean geographic distances based on type of distance measure applied among population pairs, # excluding the sink population Pairwise_Distances<-Pairwise_Distance[-Npops,-Npops] variable[40]<-sum(Pairwise_Distances)/((Npops-1)*(Npops-1)-(Npops-1)) # Calculates the maximum distances among population pairs, excluding the sink population variable[41]<-max(Pairwise_Distances) # Calculates the minimum distances among population pairs, excluding the sink population for (i in 1:(Npops-1)) { for (j in 1:(Npops-1)) { if (Pairwise_Distances[i,j]==0) { Pairwise_Distances[i,j]<-500000000 } } } variable[42]<-min(Pairwise_Distances) gc() # Reads the fecundities and survival rates # This codes works for female-only models with one or more stage matrices, and where there are no usual transitions (e.g. multiple transitions from one stage to multiple other stages) Average<-matrix(0,N_stages,N_stages) if(N_stages>1) { if(Number_matrices==1) { if(any(InputFile[38]=="OnlyFemale" |InputFile[38]=="OnlyMale" | InputFile[38]=="AllIndividuals")) { for(i in 1:N_stages) { Vital_rates[i]<-Stage_matrices[[1]][1,i] } for (i in 1:(N_stages-1)) { Vital_rates[i+N_stages]<-Stage_matrices[[1]][i+1,i] } Vital_rates[N_stages*2]<-Stage_matrices[[1]][N_stages,N_stages] } else { for (i in 1:N_stages) { Vital_rates[i]<-NA } } } else { if(any(InputFile[38]=="OnlyFemale" |InputFile[38]=="OnlyMale" | InputFile[38]=="AllIndividuals")) { for(i in 1:N_stages) { for(j in 1:N_stages) { for(k in 1:Number_matrices) { Average[i,j]<-Average[i,j]+Stage_matrices[[k]][i,j] } } } Average<-Average/Number_matrices for(i in 1:N_stages) { Vital_rates[i]<-Average[1,i] } for (i in 1:(N_stages-1)) { Vital_rates[i+N_stages]<-Average[i+1,i] } Vital_rates[N_stages*2]<-Average[N_stages,N_stages] } else { for (i in 1:N_stages) { Vital_rates[i]<-NA } } } } # Calculates population growth rate for unstructured models (i.e. one stage) if(N_stages==1) { if(Number_matrices==1) { Vital_rates<-Stage_matrices[[1]] } else { for(k in 1:Number_matrices) { Average<-Average+Stage_matrices[[k]] } Average<-Average/Number_matrices Vital_rates<-Average } } # Ensures vital rates are rounded to 6 decimal places) if(mode(Vital_rates)=="numeric") { Vital_rates<-round(Vital_rates, digits = 6) } # Binds and writes variables to variables.txt and vital_rates.txt files data.frame(do.call("cbind",variable)) write.table(variable,file="variables.txt",append=TRUE,sep=",",col.names = FALSE, row.names= FALSE) write.table(Vital_rates, file="vital_rates.txt", append=TRUE, quote=FALSE, sep=",", row.names=FALSE,col.names=FALSE) } #---END loop-----------------------# # This brace encloses all the code within a "replicate simulation loop" #-----------------------------------------# #-----------------------------------------# #-----------------------------------------# # Initiates the batch file to run all # # replicate simulations in RAMAS Metapop # #-----------------------------------------# system("batch_file.bat", wait = TRUE) #-----------------------------------------# # Collates and writes predicted response # # variables of potential interest from # # the results files created by # # RAMAS Metapop # #-----------------------------------------# # Creates a table called predictions.txt Predictions <- list("Final_Meta_N_20_yrs","Final_Meta_N_100_yrs", "Final_Meta_Pext_20_yrs","Final_Meta_Pext_100_yrs") write.table(Predictions, file="predictions.txt", col.names = FALSE, row.names= FALSE, sep=",", quote=FALSE, append=FALSE) Prediction_values<-matrix(,1,4) #-----------------------------------------# # Extract results from all # # replicate simulations # #-----------------------------------------# for(y in 1:Nreps) { # Calculates the final number of individuals in the metapopulation after 20 years (excluding sink population) Abundance<-paste("Abund", "_",y,".txt", sep="") TempAbund_1<-scan(Abundance, skip=35,nlines=1) Prediction_values[1]<-TempAbund_1[4] # Calculates the final number of individuals in the metapopulation after 20 years (excluding sink population) TempAbund_2<-scan(Abundance, skip=115,nlines=1) Prediction_values[2]<-TempAbund_2[4] # Calculates the probability of metapopulation extinction after 20 years QuasiExtinction<-paste("QuasiExt", "_",y,".txt", sep="") TempQuasiExt_1<-scan(QuasiExtinction, skip=25,nlines=1) Prediction_values[3]<-TempQuasiExt_1[3] # Calculates the probability of metapopulation extinction after 100 years TempQuasiExt_2<-scan(QuasiExtinction, skip=65,nlines=1) Prediction_values[4]<-TempQuasiExt_2[3] #Writes the response variables for the replicate simulation to the file predictions.txt write.table(Prediction_values, file="predictions.txt", append=TRUE, quote=FALSE, sep=",", row.names=FALSE,col.names=FALSE) } #-----------------------------------------# # Collates the predictions, variables and # # vital rates into a single table # #-----------------------------------------# Table_a <- read.table("predictions.txt", header = TRUE, sep = ",", na.strings = "NA") Table_b <- read.table("variables.txt", header = TRUE, sep = ",", na.strings = "NA") Table_c <- read.table("vital_rates.txt", header = TRUE, sep = ",", na.strings = "NA") Table_d <- read.table("PTCvariables.txt", header = TRUE, sep = ",", na.strings = "NA") Table_e <- read.table("spatialvar.txt", header = TRUE, sep = ",", na.strings = "NA") # Bind all together Table_f <- cbind(Table_a, Table_b, Table_c, Table_d, Table_e) write.table(Table_f, file="data.csv", col.names=TRUE, row.names=FALSE, sep = ",") #---END GRIP 2.0 Code---------------------# #-----------------------------------------#