################################################################################ # script to simulate stochastic spread of HWA under current climate # and associated loss of hemlock # # written by MC Fitzpatrick # at Harvard Forest, Petersham, MA and the Appalachian Lab, Frostburg, MD # 2009-2011 # # Model described in Fitzpatrick et al. (in press) Ecological Applications # # DESCRIPTION & NOTES # The code needs the following to run: # (1) map of hemlock abundance (available from Harvard Forest LTER archive) # (2) annual maps of winter temperature (also from Harvard Forest archive) # (3) a working directory (e.g., /SimOut) containing a folder named "modelOut" # # Code is provided as is, without support ################################################################################ ################################################################################ # CHUNK 1 - Helper functions for dispersal to speed up simulations ################################################################################ # progrediens generation long-distance dispersal dispersePR <- function(X){ props <- round(Npr[[X]] * dispEvents[[X]]) # number of dispersers direction <- runif(props, min=-pi, max=pi) # random direction for dispersal distance <- rlnorm(props, 8.46222, 1.183941) # Coords of dispersed propagules propCoords <- cbind(xFromCell(tsuga_abd, X) + distance * sin(direction), yFromCell(tsuga_abd, X) + distance * cos(direction)) # remove props that land outside of study area propCoords <- propCoords[propCoords[,1]>=xmin(tsuga_abd) & propCoords[,1]<=xmax(tsuga_abd) & propCoords[,2]>=ymin(tsuga_abd) & propCoords[,2]<=ymax(tsuga_abd) ,]} # sistens generation long-distance dispersal disperseSIS <- function(X){ props <- round(Nsis[[X]] * dispEvents[[X]]) # number of dispersers direction <- runif(props, min=-pi, max=pi) # random direction for dispersal distance <- rlnorm(props, 8.46222, 1.183941) # Coords of dispersed propagules propCoords <- cbind(xFromCell(tsuga_abd, X) + distance * sin(direction), yFromCell(tsuga_abd, X) + distance * cos(direction)) # remove props that land outside of study area propCoords <- propCoords[propCoords[,1]>=xmin(tsuga_abd) & propCoords[,1]<=xmax(tsuga_abd) & propCoords[,2]>=ymin(tsuga_abd) & propCoords[,2]<=ymax(tsuga_abd) ,]} ################################################################################ ################################################################################ # CHUNK 2 - Set working directory, load libraries & necessary data sets + code ################################################################################ source("....../spreadModel_auxiliary_functions.R") library(raster) library(maptools) #set working directory & load libraries root <- "/.../" setwd(root) #READ IN SPREAD SURFACES # Hemlock abundance in basal area: m2 ha-1, corrected for forest cover # Raster available from Harvard Forest LTER data archive tsuga_abd <- raster("tsuga_BA.img") tsuga_BA <- getValues(tsuga_abd) # List of temperature rasters for each year of simulation # Rasters available from Harvard Forest LTER data archive winterTempRasts <- list.files("/.../", pattern=".img$", full.names=T) nIter <- 1000 #number of stochastic simulations to perform iter=1 # This variable keeps track of the number of completed simulations # a csv file is written that records the number of cells # infested at each time steo, when nIter of these files have been written # the simulations stop simFiles <- length(list.files(".../SimOut/modelOut", pattern=".csv$")) ################################################################################ ################################################################################ # CHUNK 3 - Initialize simulation & convert basal area as needed ################################################################################ while(simFiles<=nIter){ setwd(paste(root, ".../SimOut", sep="")) # Directory where model output is written # length of simFiles is the number of successful simulations completed simFiles <- length(list.files("modelOut", pattern=".csv$")) ############################################################################ # Two lines of code allow simulations to run on multiple processors without # overwriting output from ongoing simulations on other processors; writes a # file to disk and assigns index to file names for each simlation iterFile <- paste(iter, ".aaa", sep="") if(!file.exists(iterFile)){write(iter, iterFile) ############################################################################ # Initial parameters startYear = 1951 endYear = 2009 nYears <- (endYear - startYear) + 1 # number of years # Matrices and vectors to hold output N <- matrix(0, nrow(tsuga_abd), ncol(tsuga_abd)) # populations n <- Npr <- Nsis <- N.out <- t(N) # matrices to hold populations n.cells <- NULL infestDate <- rep(0, length(tsuga_BA)) # vector to hold date cells # first become infested # HWA carrying capacity: assume K is related to the number of needles # amount of hemlock basal area (m2) in each 1 km2 cell (100 ha in 1 km2) tsuga_ba <- tsuga_BA*100 # leaf area (m2) from ba in cm2, as in Kenefic & Seymour 1999 LA <- ifelse(tsuga_ba > 0, (tsuga_ba*100*100)*0.1454+10.8264, 0) # number of needles, 25 mm2 average area per needle from # J. Turner REU study needlesMax <- needles <- round(LA / (25 / 1000 / 1000 )) # calculate area of cell that is covered by hemlock crown basalDiam <- 2*sqrt(tsuga_ba / pi)*100 # total basal diameter in cm # crown width (m), model from data in Santee & Monk 1981 CW <- ifelse(basalDiam == 0, 0, 1.9034 + 0.1962*basalDiam) # total crown area, crownArea <- (pi * (CW/2)^2) # infest cell near Richmond, VA with random initial population size N[rowFromCell(tsuga_abd, cellFromXY(tsuga_abd, c(1654065.944, 19831.545))), colFromCell(tsuga_abd, cellFromXY(tsuga_abd, c(1654065.944, 19831.545)))] <- rpois(1, lambda=1000) n.t <- t(N) # transpose N matrix to match raster indexing ################################################################################ ################################################################################ # CHUNK 4 - Main loop to perform population and dispersal dynamics ################################################################################ for(j in 2:nYears){ w=j-1 # index simYear <- j+1949 cells <- which(n.t > 0) # infested cells # coordinates of infested cells coords <- xyFromCell(tsuga_abd, cells) # if statement kills simulation if HWA fails to spread after 20 # generations if(length(cells) < 6 & j > 21){j=nYears} else{ # date location became infested infestDate[cells] <- ifelse(infestDate[cells]==0, simYear, infestDate[cells]) # track number of infested cells for writing to disk n.cells[j-1] <- length(cells) ncells <- as.matrix(n.cells) colnames(ncells) <- paste("sim", iter, sep="") # carrying capacity of progrediens Kpr <- needles ################################################################ # Tree death:calculate carrying capacity & reduce K as trees die ################################################################ if(length(cells) < 3){ newGrowth <- needles*runif(length(LA),0.01,0.05)} else{ # reduce trees as HWA increases # treeDeath = rate at which trees decline treeDeath <- ifelse(newGrowth == 0, 1, (1.5*N.out)/(Ksis + 1)) needles <- round(needles*(1-treeDeath)) newGrowth <- needles*runif(length(LA),0.01,0.05)} ################################################################ # carrying capacity of sistens # Assumes ~85% of sistens can settle on new growth # and the remaining ~15% on previous year's growth # A. Paradis diss. chapter II Ksis <- newGrowth*0.85 + needles*0.15 # percent of original basal area in stand that is alive perBAalive <- needles / (needlesMax+1) # total basal diameter in cm basalDiamAlive <- 2*sqrt((tsuga_ba*perBAalive) / pi)*100 CWalive <- ifelse(basalDiamAlive == 0, 0, 1.9034 + 0.1962*basalDiamAlive) crownArea <- (pi * (CWalive/2)^2) # crownArea that is alive ################################################################ # Progrediens generation ################################################################ # sistens fecundity, A. Paradis diss. eggsSis <- rpois(length(which(n.t>0)), 142.69) Npr <- n.t # population of progrediens generation # n.t is size of overwintering sistens population Npr[which(n.t>0)] <- n.t[which(n.t>0)]*eggsSis ################################## # Progrediens long-distance dispersal ################################# # proportion of individuals subject to LDD # when <3 cells are infested, typically need to send out more # dispersers or pop's rarely spread dispEvents <- if(length(cells) < 3){runif(length(Npr), 0, 0.0001)} else{runif(length(Npr), 0, 0.00000001)} # cells with enough individuals for dispersal to occur dispCells <- which(Npr*dispEvents>1) # write some text to screen to see how things are going cat("", fill=TRUE) cat("RUN", simFiles+1, fill=TRUE) cat("YEAR =", simYear, fill=TRUE) cat("NUMBER OF INFESTED STANDS (PR) =", length(which(Npr>0)), fill=TRUE) cat("NUMBER OF DISPERSING STANDS (PR) = ", length(dispCells), fill=TRUE) # disperse propagules from infested locations if(length(dispCells)>0){ # give list of propCoords for each dispersing cell propXY <- lapply(dispCells, dispersePR) # find cell numbers of where propagules landed infested <- as.vector(na.omit(unlist(lapply(propXY, cellInfest)))) # to speed things up, don't follow cells without hemlock infestedReduced <- infested[which(crownArea[infested]>0)] if(length(infestedReduced) != 0){ # number of times a propagule lands in each cell settleRate <- tabulate(infestedReduced) settleRate <- settleRate[which(settleRate>0)] # probability that HWA finds hemlock = fraction of # cell that is covered by hemlock crown probTsuga <- crownArea[unique(infestedReduced)]/1000000 estabCount <- rbinom(length(unique(infestedReduced)), settleRate, probTsuga) # assign number of dispersal events to cells in n n[unique(infestedReduced)] <- estabCount # new population size ='s Nsis + new arrivals # from dispersal (n) Npr <- n + Npr} } ################################# ################################## # Progrediens local diffusion ################################# diffEvents <- runif(length(Npr), 0, 0.000001) diffCells <- which(Npr*diffEvents>1) if(length(diffCells)>0){ infestableCells <- which(tsuga_BA>0) neighbors <- adjacency(tsuga_abd, diffCells, infestableCells, 8) Ndiff <- Npr[diffCells]*diffEvents[diffCells] estabCount <- round(Npr[neighbors[,1]]*diffEvents[neighbors[,1]]*crownArea[neighbors[,2]]/1000000) Npr[neighbors[,2]] <- Npr[neighbors[,2]] + estabCount} ################################# ################################## # Progrediens mortality ################################# #progrediens survival from A. Paradis diss. perKpr <- Npr[which(Npr>0)]/Kpr[which(Npr>0)] perKpr <- ifelse(perKpr>1,1,perKpr) Npr[which(Npr>0)] <- rbinom(length(Npr[which(Npr>0)]), Npr[which(Npr>0)], -0.9883333*perKpr+1) # Reduce populations to carrying capacity Npr[which(Npr>0)] <- round(ifelse(Npr[which(Npr>0)] > Kpr[which(Npr>0)], Kpr[which(Npr>0)], Npr[which(Npr>0)])) ################################# # make new, empty dispersal matrix n <- matrix(0, nrow=2571, ncol=2158) ################################################################ # Sistens generation ################################################################ # progrediens fecundity, A. Paradis diss. eggsPr <- rpois(length(which(Npr>0)), 22.21) Nsis <- Npr Nsis[which(Npr>0)] <- Npr[which(Npr>0)]*eggsPr ################################## # Sistens long-distance dispersal ################################# # proportion of individuals subject to LDD # when <3 cells are infested, typically need to send out more # dispersers or pop's rarely spread dispEvents <- if(length(cells) < 3){runif(length(Nsis), 0, 0.001)} else{runif(length(Nsis), 0, 0.00000001)} # cells with enough individuals for dispersal to occur dispCells <- which(Nsis*dispEvents>1) cat("NUMBER OF INFESTED STANDS (SIS) =", length(which(Nsis>0)), fill=TRUE) cat("NUMBER OF DISPERSING STANDS (SIS) = ", length(dispCells), fill=TRUE) # disperse propagules from infested locations if(length(dispCells)>0){ # give list of propCoords for each dispersing cell propXY <- lapply(dispCells, disperseSIS) # find cell numbers of where propagules landed infested <- as.vector(na.omit(unlist(lapply(propXY, cellInfest)))) # to speed things up, don't follow cells with no hemlock infestedReduced <- infested[which(crownArea[infested]>0)] if(length(infestedReduced) != 0){ # number of times a propagule lands in each cell settleRate <- tabulate(infestedReduced) settleRate <- settleRate[which(settleRate>0)] # probability that HWA finds hemlock = fraction of # cell that is covered by hemlock crown probTsuga <- crownArea[unique(infestedReduced)]/1000000 estabCount <- rbinom(length(unique(infestedReduced)), settleRate, probTsuga) # assign number of dispersal events to cells in n n[unique(infestedReduced)] <- estabCount # new population size ='s Nsis + new arrivals # from dispersal (n) Nsis <- n + Nsis} } ################################# ################################### ## Sistens local diffusion ################################## diffEvents <- runif(length(Nsis), 0, 0.000001) diffCells <- which(Nsis*diffEvents>1) if(length(diffCells)>0){ infestableCells <- which(tsuga_BA>0) neighbors <- adjacency(tsuga_abd, diffCells, infestableCells, 8) Ndiff <- Nsis[diffCells]*diffEvents[diffCells] estabCount <- round(Nsis[neighbors[,1]]*diffEvents[neighbors[,1]]*crownArea[neighbors[,2]]/1000000) Nsis[neighbors[,2]] <- Nsis[neighbors[,2]] + estabCount} ################################# ################################## # Sistens mortality ################################# # sistens survival from A. Paradis perKsis <- Nsis[which(Nsis>0)]/Ksis[which(Nsis>0)] perKsis <- ifelse(perKsis>1,1,perKsis) Nsis[which(Nsis>0)] <- rbinom(length(Nsis[which(Nsis>0)]), Nsis[which(Nsis>0)], -0.644*perKsis+1) # Reduce populations to carrying capacity Nsis[which(Nsis>0)] <- round(ifelse(Nsis[which(Nsis>0)] > Ksis[which(Nsis>0)], Ksis[which(Nsis>0)], Nsis[which(Nsis>0)])) # sistens aestivation survival, A. Paradis Nsis <- round(rbinom(length(Nsis),Nsis, 0.287666667)) ################################# # make new, empty dispersal matrix n <- matrix(0, nrow=2571, ncol=2158) N.out <- c(Nsis) ################################################################ # WINTER MORTALITY ################################################################ # winter mortality based on equation from Paradis et al 2008. # and winter temperature (DJFM), the model which had highest R2. # Winter (DJFM) temperature setwd(root) # read in temperature raster for year of sim # climate rasters available from Harvard Forest LTER data archive cold <- raster(winterTempRasts[w]) winterTemp <- getValues(cold) coeffs <- runif(length(which(N.out>0)), -0.104, -0.052) intc <- runif(length(which(N.out>0)), 0.402, 0.612) winterMort <- winterTemp[which(N.out>0)]*coeffs+intc N.out[which(N.out>0)] <- round(N.out[which(N.out>0)]*(1-winterMort)) ################################################################ setwd(paste(root, "/simOut", sep="")) n.t <- matrix(1,2571,2158)*N.out # write csv file recording infestation area (# of cells) if(j>30){write.csv(ncells, paste(getwd(), ".../modelOut/n.infested.cells.", iter, ".csv", sep=""), row.names=FALSE)} # WRITE OUTPUT TO FILE outYears <- c(2008, 2003, 1998, 1993, 1988, 1983, 1978, 1973) if(simYear %in% outYears){ #write population raster to file Nrast <- setValues(tsuga_abd, N.out) # assign population sizes to Nrast Nrast <- writeRaster(Nrast, filename=paste("Nrast", "_", simYear, "_", iter, ".tif", sep=""), format="GTiff", overwrite=TRUE) } if(j==nYears){ #write tree death to file newBA <- setValues(tsuga_abd, perBAalive*100) newBA <- writeRaster(newBA, filename=paste("perBAalive", "_", simYear, "_", iter, ".tif", sep=""), format="GTiff", overwrite=TRUE, datatype="INT1U") infestGen <- setValues(tsuga_abd, infestDate) infestGen <- writeRaster(infestGen, filename=paste("infestDate", "_", simYear, "_", iter, ".tif", sep=""), format="GTiff", overwrite=TRUE, datatype="INT2U") save.image(paste("hwa_spread_", iter, ".Rdata", sep=""))} } } } else{iter = iter + 1} # if another processors is on this interation, skip } # to next one and proceed... ################################################################################