R/IFNdynamics.R
In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis
Defines functions
IFNdynamics
Documented in
IFNdynamics
#' Forest dynamics
#'
#' Simulates one or several steps of forest dynamics, including growth, mortality and ingrowth processes.
#'
#' @param treeData A data frame with tree records in rows and columns 'ID', 'Species', 'DBH' (in cm) and 'N' (ha-1). Additionally, column 'OIFFIN' may be defined
#' to indicate trees that are signaled to be cut during the following step.
#' @param plotData A data frame with plot 'ID' as row names and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C), 'Prec' (mm), 'PET' (mm) and 'slope' (degrees).
#' An additional column 'Step' with integer values ('1', '2', etc.) may be used to indicate temporal variation in environmental conditions. In this case
#' 'ID' should also be provided as column, not row names, because plot IDs will not be unique. An additional column 'Province' may be used to specify Spanish province.
#' @param nb An object of class \code{listw} (see functions \code{\link{IFNknn}} and \code{\link{nb2listw}}).
#' @param NspNmatrix A numeric matrix with the density proportion of each species (columns) in each plot (rows).
#' @param numYears An integer or integer vector with the number of years corresponding to one time step (by default 10-yr time steps are assumed)
#' @param numSteps Number of time steps to simulate (normally, 1 step = 10 years). If temporal variation of environmental data is desired, values of column 'Step' of \code{plotData}
#' should be specified accordingly.
#' @param useProvince A flag to indicate that province information should be used in growth/height estimation.
#' @param useRP A flag to indicate that "region de procedencia" should be used in growth/height estimation
#' @param mergeTrees Flag to indicate that tree records of similar diameter, low density and the same species (not including ingrowth of the current step) should be merged to avoid accumulation of tree records.
#' @param killVeryLargeTrees Flag to indicate that trees with diameters larger than species-specific thresholds should be forced to die.
#' @param pastGrowthMortalityEffect Flag to indicate that the model should use previous diameter increment to determine mortality (whenever possible).
#' @param applyMortalityInManagedStands Flag to indicate that mortality should also be applied in stands with observed management (BAL.ext).
#' @param min.N Minimum density (individuals/ha) to keep a tree record
#' @param removeCutTrees Boolean flag to indicate that trees signaled as cut in 'OIFFIN' should be removed from the output
#' @param stochastic Boolean flag to indicate that simulation of growth, mortality and ingrowth should be stochastic (i.e. non-deterministic)
#' @param fire Boolean flag to indicate that simulations should include stochastic burning of plots (see \code{\link{IFNfire}}).
#' @param numRuns Number of stochastic runs (only applicable if \code{stochastic = TRUE})
#' @param ingrowth Boolean flag to simulate ingrowth during simulation
#' @param thresholdIngrowth Boolean flag to indicate that ingrowth should occur only when a given probability threshold is surpassed.
#' @param height Boolean flag to estimate height of trees after dynamic simulation
#' @param heightgrowth Boolean flag to estimate height increment whenever possible, instead of using static height models
#' @param sequence Boolean flag to indicate that the sequence of tree lists by step should be returned, instead of only the final tree list.
#' If \code{sequence = TRUE}, an additional column 'Step' is added to the output.
#' @param variable.radius A flag to indicate that sampling is with variable radius (this forces changes in density when
#' growth in tree diameters involves changes in sampling radius).
#' @param dead.table Boolean flag to indicate that a table of trees dead during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees died.
#' @param recr.labels Boolean flag to indicate whether ingrowth records should be labelled in the output.
#' \code{recr.labels = TRUE} results in a true/false column 'Ingrowth' being added to the output.
#' @param speciesNames Boolean flag to indicate that species names should be added to the output.
#' @param verbose Boolean flag to indicate console output of the process
#'
#' @return If \code{dead.table=FALSE} it returns a data frame with tree records in rows and columns 'ID', ('Step'), 'Species', ('Name'), 'DBH', 'N' , ('Ingrowth'), and ('H') depending on the selected output options
#' (i.e. \code{sequence}, \code{speciesNames}, \code{height} and \code{recr.labels}).
#' If \code{dead.table=FALSE}, it returns a list with two element, one named \code{out} with live trees (at the end of the simulation or, if \code{sequence=TRUE} the sequence of trees for all steps)
#' and another called \code{dead} with dead trees for each time step.
#'
#' @details The function checks whether plot IDs in \code{treeData} exist as row names of \code{plotData}
#' and whether species codes in \code{treeData} have models defined. If \code{plotData} contains field 'Province', this will
#' be used to overwrite any province especification in \code{treeData}. Option \code{stochastic = TRUE} causes the models
#' to be run in stochastic mode, which means that: (a) binomial draws are made for each tree to determine survival; (b)
#' predicted growth (e.g. diameter increment) is stochastic (Gamma distribution); (c) Ingrowth density is stochastic (Gamma distribution
#' coupled to a binomial draw); (d) Ingrowth diameter is stochastic (Gamma distribution).
#'
#' When running the function with spatial dependencies, local recruitment is still performed using submodel 'IFNingrowth' but
#' colonization of new species is determined using submodel 'IFNingrowthdisp' (because this latter submodel underestimates recruitment from species
#' already present in the plot).
#' @name IFNdynamics
#' @aliases IFNdynamics
#' @seealso \code{\link{IFNsubmodels}}
#'
#' @examples
#' # Load example tree data (two forest plots)
#' data(exampleTreeData)
#' head(exampleTreeData)
#'
#' # Load IFN environmental data
#' data(plotDataIFN)
#' head(plotDataIFN)
#'
#' # Load Regiones de Procedencia for IFN plots
#' data(plotRPIFN)
#' # Add Regiones de Procedencia to plot data
#' plotDataIFN = cbind(plotDataIFN, plotRPIFN)
#'
#' # Run
#' IFNdynamics(exampleTreeData, plotDataIFN)
#'
#' # Run with dispersal effects
#' data(examplePlotCoords)
#' nb = IFNknn(examplePlotCoords, k = 1)
#' Nmatrix = tapply(exampleTreeData$N, list(exampleTreeData$ID, exampleTreeData$Species),
#' sum, na.rm=T)
#' NspNmatrix = sweep(Nmatrix, 1, rowSums(Nmatrix, na.rm=T), "/")
#' NspNmatrix = rbind(NspNmatrix,c(0.5, NA, NA, NA))
#' rownames(NspNmatrix)[4] = "80003"
#' NspNmatrix = cbind(NspNmatrix, c(NA,NA,NA,0.5))
#' colnames(NspNmatrix)[5] = "25"
#' IFNdynamics(exampleTreeData, plotDataIFN, nb = nb, NspNmatrix = NspNmatrix)
IFNdynamics<-function(treeData, plotData,
nb = NULL, NspNmatrix = NULL,
numYears = 10, numSteps = 1, min.N = 1e-5,
ingrowth = TRUE, thresholdIngrowth = FALSE,
height = TRUE, heightgrowth = FALSE,
useProvince = TRUE, useRP = TRUE,
mergeTrees = FALSE,
killVeryLargeTrees = TRUE, removeCutTrees = TRUE, pastGrowthMortalityEffect = FALSE,
applyMortalityInManagedStands = TRUE,
stochastic = FALSE,
fire = FALSE,
numRuns = 1,
variable.radius = FALSE,
sequence = FALSE, dead.table = FALSE, recr.labels = TRUE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check species codes
.checkSpeciesCodes(treeData)
spatial = ((!is.null(nb)) && (!is.null(NspNmatrix)))
if(pastGrowthMortalityEffect) {
if("DBHincprev" %in% names(treeData)) DBHincprev = treeData$DBHincprev
else DBHincprev = rep(NA, nrow(treeData))
}
BALext = FALSE
if((numSteps==1) && ("OIFFIN" %in% names(treeData))) {
BALext = TRUE
if(verbose) {
cat("\n Calculating dynamics with BALext.\n")
}
}
if(!stochastic) numRuns = 1
treeDataOut = NULL
treeDataOutRun = NULL
treeDataDead = NULL
treeDataDeadRun = NULL
treeDataInitial = treeData
# Iterate over runs
for(r in 1:numRuns) {
# Recover initial data
treeData = treeDataInitial
# Store initial state as sequence
if(sequence) {
treeDataSeq = data.frame(Run = rep(r, nrow(treeData)),
Step = rep(0, nrow(treeData)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH,
N = treeData$N, stringsAsFactors = F)
if(height) { # Copy initial height or estimate it
if("H" %in% names(treeData)) treeDataSeq$H = treeData$H
else {
treeDataSeq$H = NA
}
}
if(pastGrowthMortalityEffect) treeDataSeq$DBHincprev = DBHincprev
if(recr.labels) treeDataSeq$Ingrowth = FALSE
}
# Iterate over steps
for(s in 1:numSteps) {
if(nrow(treeData)>0) {
if(verbose) {
cat(paste0("\n\n ---- Step ",s,"/",numSteps," ----\n"))
cat("\n Preparing tree & competition data...\n")
}
dataInput = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = BALext, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
#Subset environmental data corresponding to step c("ID","SWHC","Rad","Temp","Prec", "slope")
if("Step" %in% names(plotData)) {
plotDataStep = plotData[plotData$Step==s,, drop= FALSE]
if(nrow(plotDataStep)==0) stop(paste0("Environmental variables not defined for step ", s,"."))
# if("Province" %in% names(plotData)) plotDataStep$Province = plotData$Province[plotData$Step==s]
row.names(plotDataStep) = plotDataStep$ID
} else {
plotDataStep = plotData
# if("Province" %in% names(plotData)) plotDataStep$Province = plotData$Province
}
## Check IDs in environmental data
.checkIDs(treeData, plotDataStep)
if(useProvince) if("Province" %in% names(plotDataStep)) dataInput$Province = plotDataStep[as.character(dataInput$ID),"Province"]
dataInput$slope = plotDataStep[as.character(dataInput$ID),"slope"]
dataInput$SWHC = plotDataStep[as.character(dataInput$ID),"SWHC"]
dataInput$Rad = plotDataStep[as.character(dataInput$ID),"Rad"]
dataInput$Temp = plotDataStep[as.character(dataInput$ID),"Temp"]
dataInput$PET = plotDataStep[as.character(dataInput$ID),"PET"]
dataInput$Prec = plotDataStep[as.character(dataInput$ID),"Prec"]
dataInput$PPET = pmin(dataInput$Prec/dataInput$PET, 1.3)
dataInput$Temp15sq = (dataInput$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(dataInput$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
dataInput[[rp]] = plotDataStep[as.character(dataInput$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
dataInput[spns!=sp,rpmap] = NA
}
}
}
if(!spatial) { # Without dispersal
if(verbose) {
cat("\n Preparing species & competition data...\n")
}
spDataInput = prepareSpeciesCompetitionTable(treeData,
provinceFromID = TRUE, verbose = verbose)
} else { # With dispersal
if(verbose) {
cat("\n Preparing species & competition data (including dispersal)...\n")
}
spDataInput = prepareSpatialSpeciesCompetitionTable(treeData,
nb = nb, NspNmatrix = NspNmatrix,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotDataStep)) spDataInput$Province = plotDataStep[as.character(spDataInput$ID),"Province"]
spDataInput$slope = plotDataStep[as.character(spDataInput$ID),"slope"]
spDataInput$SWHC = plotDataStep[as.character(spDataInput$ID),"SWHC"]
spDataInput$Rad = plotDataStep[as.character(spDataInput$ID),"Rad"]
spDataInput$Temp = plotDataStep[as.character(spDataInput$ID),"Temp"]
spDataInput$PET = plotDataStep[as.character(spDataInput$ID),"PET"]
spDataInput$Prec = plotDataStep[as.character(spDataInput$ID),"Prec"]
spDataInput$PPET = pmin(spDataInput$Prec/spDataInput$PET, 1.3)
spDataInput$Temp15sq = (spDataInput$Temp-15)^2
# Substep 1: Growth
INC = .IFNgrowthInternal(dataInput, numYears = numYears, response = "DI",
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose=verbose)
if(heightgrowth) {
HINC = .IFNheightgrowthInternal(dataInput, numYears = numYears,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose=verbose)
newH = dataInput$h + HINC
}
# Substep 2: Survival
Psurv = .IFNsurvivalInternal(dataInput, numYears = numYears, killVeryLargeTrees = killVeryLargeTrees, verbose = verbose)
if(length(Psurv)!=nrow(dataInput)) stop("Unequal size for Psurv and dataInput")
if(pastGrowthMortalityEffect) {
if(nrow(dataInput)!=length(DBHincprev)) stop("Unequal size for DBHincprev and dataInput")
dataInput$ln.d.inc.prev = log(DBHincprev+1)
dataInput$ln.d.inc.prev.ln.d = dataInput$ln.d.inc.prev/dataInput$ln.d
selnoNA = !is.na(DBHincprev)
if(sum(selnoNA)>0) {
PsurvPG = .IFNsurvivalPGInternal(dataInput[selnoNA,,drop=FALSE], numYears = numYears, killVeryLargeTrees = killVeryLargeTrees, verbose = verbose)
if(length(PsurvPG)!=sum(selnoNA)) stop("Unequal size for PsurvPG and selnoNA")
if(sum(!is.na(PsurvPG))>0) Psurv[selnoNA][!is.na(PsurvPG)] = PsurvPG[!is.na(PsurvPG)]
}
}
if(!applyMortalityInManagedStands) {
if("BAL.ext" %in% names(dataInput)) {
ids_man = unique(dataInput$ID[dataInput$BAL.ext>0])
Psurv[dataInput$ID %in% ids_man] = 1
}
}
if(fire) { ## If stand burn then survival is not density-dependent but depends on fire
fireOutput = IFNfire(treeData, plotDataStep)
treePsurvfire = fireOutput$tree$Psurvfire
# if(sum(!is.na(treePsurvfire))>0) print(cbind(Psurv,treePsurvfire))
# Multiply survival probabilities for those stands that burned
Psurv[!is.na(treePsurvfire)] = Psurv[!is.na(treePsurvfire)]*treePsurvfire[!is.na(treePsurvfire)]
}
# Substep 3: Apply growth and survival
Nnew = treeData$N
if(stochastic) {
# Apply a binomial survival random variable
for(j in 1:length(Nnew)) Nnew[j] = min(Nnew[j], rbinom(1, prob = Psurv[j], size=round(Nnew[j])))
} else {
# Multiply survival probability by previous density estimate
Nnew = treeData$N*Psurv
}
treeDataOutStep = data.frame(Run = rep(r, length(treeData$ID)),
Step = rep(s, length(treeData$ID)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH + INC,
N = Nnew, stringsAsFactors = FALSE)
if(pastGrowthMortalityEffect) treeDataOutStep$DBHincprev = INC
if(BALext && removeCutTrees) { # If BALext was supplied remove cut trees
if(verbose) {
cat("\n Removing cut trees...\n")
}
if(heightgrowth) newH = newH[treeData$OIFFIN!="0",]
treeDataOutStep = treeDataOutStep[treeData$OIFFIN!="0",]
}
# Apply variable radius correction (if necessary)
if(variable.radius) {
dfr = .densityFactor(treeDataOutStep$DBH)/.densityFactor(treeData$DBH)
treeDataOutStep$N = treeDataOutStep$N*dfr
}
# Merge small trees (to avoid accumulation of records)
if(mergeTrees) {
treeDataOutStep = .mergeTrees(treeDataOutStep)
}
# Assumes dead trees did not grow in diameter
if(dead.table) {
treeDataDeadStep = data.frame(Run = rep(r, nrow(treeData)),
Step = rep(s, nrow(treeData)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH,
N = treeData$N - Nnew, stringsAsFactors = FALSE)
if(pastGrowthMortalityEffect) treeDataDeadStep$DBHincprev = INC
}
# Substep 4: Add ingrowth
if(ingrowth) {
#Local recruitment
col = (spDataInput$Nsp.N == 0) # Is new colonization?
recrData = .IFNingrowthInternal(spDataInput[!col,,drop=FALSE], numYears = numYears,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose=verbose)
if(sum(col)>0) { # If there is colonization apply models with dispersal for those records
recrData_disp = .IFNingrowthdispInternal(spDataInput[col,,drop=FALSE], numYears = numYears,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose=verbose)
ndisp = nrow(recrData_disp)
recrData = rbind(recrData, recrData_disp)
} else {
ndisp = 0
}
if(pastGrowthMortalityEffect) recrData$DBHincprev = NA
if(recr.labels) {
treeDataOutStep$Ingrowth = rep(FALSE,nrow(treeDataOutStep))
treeDataOutStep$Colonization = rep(FALSE,nrow(treeDataOutStep))
recrData$Run = rep(r, nrow(recrData))
recrData$Step = rep(s, nrow(recrData))
recrData$Ingrowth = rep(TRUE, nrow(recrData))
recrData$Colonization = c(rep(FALSE, nrow(recrData)-ndisp),rep(FALSE, ndisp))
}
if(fire) { ## Assume fire kills ingrowth
Psurv = rep(1, nrow(recrData))
Psurv[fireOutput$stand[recrData$ID,"fire"]==1] = 0
recrData$N = recrData$N*Psurv
recrData = recrData[recrData$N>0, , drop = FALSE]
# recrFire = .treeFire(recrData, fireOutput$stand)
# recrPsurvfire = recrFire$Psurvfire
# Psurv = rep(1, nrow(recrData))
# Psurv[!is.na(recrPsurvfire)] = recrPsurvfire[!is.na(recrPsurvfire)] # Replace survival probabilities for those stands that burned
# recrData$N = recrData$N*Psurv
}
treeDataOutStep = rbind(treeDataOutStep, recrData[,names(treeDataOutStep)])
}
# Substep 5: Height
if(height) {
treeDataOutStep$H = NA
#Replace height predictions derived from height increment model
if(heightgrowth) treeDataOutStep$H[1:length(newH)] = newH
hstatic = IFNheight(treeDataOutStep, plotDataStep,
useProvince = useProvince, useRP = useRP, verbose=verbose)
#Replace missing values with static predictions
naH = is.na(treeDataOutStep$H)
treeDataOutStep$H[naH] = hstatic[naH]
if(dead.table) {
treeDataDeadStep$H = IFNheight(treeDataDeadStep, plotDataStep,
useProvince = useProvince, useRP = useRP, verbose=verbose)
}
}
# Substep 6: Remove low density cohorts
treeDataOutStep = treeDataOutStep[treeDataOutStep$N>min.N, , drop = FALSE]
# For next step
if(s<numSteps) {
treeData = treeDataOutStep
if(pastGrowthMortalityEffect) DBHincprev = treeData$DBHincprev
}
if(sequence) {
treeDataSeq = rbind(treeDataSeq,treeDataOutStep[,names(treeDataSeq)])
}
if(dead.table) {
if(is.null(treeDataDeadRun)) {
treeDataDeadRun = treeDataDeadStep
} else {
treeDataDeadRun = rbind(treeDataDeadRun, treeDataDeadStep)
}
}
}
}
if(sequence) {
treeDataOutRun = treeDataSeq
} else {
treeDataOutRun = treeDataOutStep
}
if(!is.null(treeDataOut)) {
treeDataOut = rbind(treeDataOut, treeDataOutRun)
} else {
treeDataOut = treeDataOutRun
}
if(!is.null(treeDataDead)) {
treeDataDead = rbind(treeDataDead,treeDataDeadRun)
} else {
treeDataDead = treeDataDeadRun
}
}
# Species names
if(speciesNames) {
if(sequence) {
treeDataOut = data.frame(Run = treeDataOut$Run,
Step = treeDataOut$Step,
ID = treeDataOut$ID,
Species = treeDataOut$Species,
Name = speciesNames(treeDataOut$Species),
treeDataOut[,-c(1:4)], stringsAsFactors = FALSE)
} else {
treeDataOut = data.frame(Run = treeDataOut$Run,
ID = treeDataOut$ID,
Species = treeDataOut$Species,
Name = speciesNames(treeDataOut$Species),
treeDataOut[,-c(1:4)], stringsAsFactors = FALSE)
}
if(dead.table) {
treeDataDead = data.frame(Run = treeDataDead$Run,
Step = treeDataDead$Step,
ID = treeDataDead$ID,
Species = treeDataDead$Species,
Name = speciesNames(treeDataDead$Species),
treeDataDead[,-c(1:4)], stringsAsFactors = FALSE)
}
}
rownames(treeDataOut) = NULL
rownames(treeDataDead) = NULL
if(verbose) cat("\n")
if(!dead.table) {
return(treeDataOut)
}
return(list(out = treeDataOut, dead = treeDataDead))
}
miquelcaceres/IFNdyn documentation built on Feb. 1, 2021, 10:55 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions IFNdynamics
Documented in IFNdynamics
#' Forest dynamics
#'
#' Simulates one or several steps of forest dynamics, including growth, mortality and ingrowth processes.
#'
#' @param treeData A data frame with tree records in rows and columns 'ID', 'Species', 'DBH' (in cm) and 'N' (ha-1). Additionally, column 'OIFFIN' may be defined
#' to indicate trees that are signaled to be cut during the following step.
#' @param plotData A data frame with plot 'ID' as row names and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C), 'Prec' (mm), 'PET' (mm) and 'slope' (degrees).
#' An additional column 'Step' with integer values ('1', '2', etc.) may be used to indicate temporal variation in environmental conditions. In this case
#' 'ID' should also be provided as column, not row names, because plot IDs will not be unique. An additional column 'Province' may be used to specify Spanish province.
#' @param nb An object of class \code{listw} (see functions \code{\link{IFNknn}} and \code{\link{nb2listw}}).
#' @param NspNmatrix A numeric matrix with the density proportion of each species (columns) in each plot (rows).
#' @param numYears An integer or integer vector with the number of years corresponding to one time step (by default 10-yr time steps are assumed)
#' @param numSteps Number of time steps to simulate (normally, 1 step = 10 years). If temporal variation of environmental data is desired, values of column 'Step' of \code{plotData}
#' should be specified accordingly.
#' @param useProvince A flag to indicate that province information should be used in growth/height estimation.
#' @param useRP A flag to indicate that "region de procedencia" should be used in growth/height estimation
#' @param mergeTrees Flag to indicate that tree records of similar diameter, low density and the same species (not including ingrowth of the current step) should be merged to avoid accumulation of tree records.
#' @param killVeryLargeTrees Flag to indicate that trees with diameters larger than species-specific thresholds should be forced to die.
#' @param pastGrowthMortalityEffect Flag to indicate that the model should use previous diameter increment to determine mortality (whenever possible).
#' @param applyMortalityInManagedStands Flag to indicate that mortality should also be applied in stands with observed management (BAL.ext).
#' @param min.N Minimum density (individuals/ha) to keep a tree record
#' @param removeCutTrees Boolean flag to indicate that trees signaled as cut in 'OIFFIN' should be removed from the output
#' @param stochastic Boolean flag to indicate that simulation of growth, mortality and ingrowth should be stochastic (i.e. non-deterministic)
#' @param fire Boolean flag to indicate that simulations should include stochastic burning of plots (see \code{\link{IFNfire}}).
#' @param numRuns Number of stochastic runs (only applicable if \code{stochastic = TRUE})
#' @param ingrowth Boolean flag to simulate ingrowth during simulation
#' @param thresholdIngrowth Boolean flag to indicate that ingrowth should occur only when a given probability threshold is surpassed.
#' @param height Boolean flag to estimate height of trees after dynamic simulation
#' @param heightgrowth Boolean flag to estimate height increment whenever possible, instead of using static height models
#' @param sequence Boolean flag to indicate that the sequence of tree lists by step should be returned, instead of only the final tree list.
#' If \code{sequence = TRUE}, an additional column 'Step' is added to the output.
#' @param variable.radius A flag to indicate that sampling is with variable radius (this forces changes in density when
#' growth in tree diameters involves changes in sampling radius).
#' @param dead.table Boolean flag to indicate that a table of trees dead during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees died.
#' @param recr.labels Boolean flag to indicate whether ingrowth records should be labelled in the output.
#' \code{recr.labels = TRUE} results in a true/false column 'Ingrowth' being added to the output.
#' @param speciesNames Boolean flag to indicate that species names should be added to the output.
#' @param verbose Boolean flag to indicate console output of the process
#'
#' @return If \code{dead.table=FALSE} it returns a data frame with tree records in rows and columns 'ID', ('Step'), 'Species', ('Name'), 'DBH', 'N' , ('Ingrowth'), and ('H') depending on the selected output options
#' (i.e. \code{sequence}, \code{speciesNames}, \code{height} and \code{recr.labels}).
#' If \code{dead.table=FALSE}, it returns a list with two element, one named \code{out} with live trees (at the end of the simulation or, if \code{sequence=TRUE} the sequence of trees for all steps)
#' and another called \code{dead} with dead trees for each time step.
#'
#' @details The function checks whether plot IDs in \code{treeData} exist as row names of \code{plotData}
#' and whether species codes in \code{treeData} have models defined. If \code{plotData} contains field 'Province', this will
#' be used to overwrite any province especification in \code{treeData}. Option \code{stochastic = TRUE} causes the models
#' to be run in stochastic mode, which means that: (a) binomial draws are made for each tree to determine survival; (b)
#' predicted growth (e.g. diameter increment) is stochastic (Gamma distribution); (c) Ingrowth density is stochastic (Gamma distribution
#' coupled to a binomial draw); (d) Ingrowth diameter is stochastic (Gamma distribution).
#'
#' When running the function with spatial dependencies, local recruitment is still performed using submodel 'IFNingrowth' but
#' colonization of new species is determined using submodel 'IFNingrowthdisp' (because this latter submodel underestimates recruitment from species
#' already present in the plot).
#' @name IFNdynamics
#' @aliases IFNdynamics
#' @seealso \code{\link{IFNsubmodels}}
#'
#' @examples
#' # Load example tree data (two forest plots)
#' data(exampleTreeData)
#' head(exampleTreeData)
#'
#' # Load IFN environmental data
#' data(plotDataIFN)
#' head(plotDataIFN)
#'
#' # Load Regiones de Procedencia for IFN plots
#' data(plotRPIFN)
#' # Add Regiones de Procedencia to plot data
#' plotDataIFN = cbind(plotDataIFN, plotRPIFN)
#'
#' # Run
#' IFNdynamics(exampleTreeData, plotDataIFN)
#'
#' # Run with dispersal effects
#' data(examplePlotCoords)
#' nb = IFNknn(examplePlotCoords, k = 1)
#' Nmatrix = tapply(exampleTreeData$N, list(exampleTreeData$ID, exampleTreeData$Species),
#' sum, na.rm=T)
#' NspNmatrix = sweep(Nmatrix, 1, rowSums(Nmatrix, na.rm=T), "/")
#' NspNmatrix = rbind(NspNmatrix,c(0.5, NA, NA, NA))
#' rownames(NspNmatrix)[4] = "80003"
#' NspNmatrix = cbind(NspNmatrix, c(NA,NA,NA,0.5))
#' colnames(NspNmatrix)[5] = "25"
#' IFNdynamics(exampleTreeData, plotDataIFN, nb = nb, NspNmatrix = NspNmatrix)
IFNdynamics<-function(treeData, plotData,
nb = NULL, NspNmatrix = NULL,
numYears = 10, numSteps = 1, min.N = 1e-5,
ingrowth = TRUE, thresholdIngrowth = FALSE,
height = TRUE, heightgrowth = FALSE,
useProvince = TRUE, useRP = TRUE,
mergeTrees = FALSE,
killVeryLargeTrees = TRUE, removeCutTrees = TRUE, pastGrowthMortalityEffect = FALSE,
applyMortalityInManagedStands = TRUE,
stochastic = FALSE,
fire = FALSE,
numRuns = 1,
variable.radius = FALSE,
sequence = FALSE, dead.table = FALSE, recr.labels = TRUE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check species codes
.checkSpeciesCodes(treeData)
spatial = ((!is.null(nb)) && (!is.null(NspNmatrix)))
if(pastGrowthMortalityEffect) {
if("DBHincprev" %in% names(treeData)) DBHincprev = treeData$DBHincprev
else DBHincprev = rep(NA, nrow(treeData))
}
BALext = FALSE
if((numSteps==1) && ("OIFFIN" %in% names(treeData))) {
BALext = TRUE
if(verbose) {
cat("\n Calculating dynamics with BALext.\n")
}
}
if(!stochastic) numRuns = 1
treeDataOut = NULL
treeDataOutRun = NULL
treeDataDead = NULL
treeDataDeadRun = NULL
treeDataInitial = treeData
# Iterate over runs
for(r in 1:numRuns) {
# Recover initial data
treeData = treeDataInitial
# Store initial state as sequence
if(sequence) {
treeDataSeq = data.frame(Run = rep(r, nrow(treeData)),
Step = rep(0, nrow(treeData)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH,
N = treeData$N, stringsAsFactors = F)
if(height) { # Copy initial height or estimate it
if("H" %in% names(treeData)) treeDataSeq$H = treeData$H
else {
treeDataSeq$H = NA
}
}
if(pastGrowthMortalityEffect) treeDataSeq$DBHincprev = DBHincprev
if(recr.labels) treeDataSeq$Ingrowth = FALSE
}
# Iterate over steps
for(s in 1:numSteps) {
if(nrow(treeData)>0) {
if(verbose) {
cat(paste0("\n\n ---- Step ",s,"/",numSteps," ----\n"))
cat("\n Preparing tree & competition data...\n")
}
dataInput = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = BALext, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
#Subset environmental data corresponding to step c("ID","SWHC","Rad","Temp","Prec", "slope")
if("Step" %in% names(plotData)) {
plotDataStep = plotData[plotData$Step==s,, drop= FALSE]
if(nrow(plotDataStep)==0) stop(paste0("Environmental variables not defined for step ", s,"."))
# if("Province" %in% names(plotData)) plotDataStep$Province = plotData$Province[plotData$Step==s]
row.names(plotDataStep) = plotDataStep$ID
} else {
plotDataStep = plotData
# if("Province" %in% names(plotData)) plotDataStep$Province = plotData$Province
}
## Check IDs in environmental data
.checkIDs(treeData, plotDataStep)
if(useProvince) if("Province" %in% names(plotDataStep)) dataInput$Province = plotDataStep[as.character(dataInput$ID),"Province"]
dataInput$slope = plotDataStep[as.character(dataInput$ID),"slope"]
dataInput$SWHC = plotDataStep[as.character(dataInput$ID),"SWHC"]
dataInput$Rad = plotDataStep[as.character(dataInput$ID),"Rad"]
dataInput$Temp = plotDataStep[as.character(dataInput$ID),"Temp"]
dataInput$PET = plotDataStep[as.character(dataInput$ID),"PET"]
dataInput$Prec = plotDataStep[as.character(dataInput$ID),"Prec"]
dataInput$PPET = pmin(dataInput$Prec/dataInput$PET, 1.3)
dataInput$Temp15sq = (dataInput$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(dataInput$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
dataInput[[rp]] = plotDataStep[as.character(dataInput$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
dataInput[spns!=sp,rpmap] = NA
}
}
}
if(!spatial) { # Without dispersal
if(verbose) {
cat("\n Preparing species & competition data...\n")
}
spDataInput = prepareSpeciesCompetitionTable(treeData,
provinceFromID = TRUE, verbose = verbose)
} else { # With dispersal
if(verbose) {
cat("\n Preparing species & competition data (including dispersal)...\n")
}
spDataInput = prepareSpatialSpeciesCompetitionTable(treeData,
nb = nb, NspNmatrix = NspNmatrix,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotDataStep)) spDataInput$Province = plotDataStep[as.character(spDataInput$ID),"Province"]
spDataInput$slope = plotDataStep[as.character(spDataInput$ID),"slope"]
spDataInput$SWHC = plotDataStep[as.character(spDataInput$ID),"SWHC"]
spDataInput$Rad = plotDataStep[as.character(spDataInput$ID),"Rad"]
spDataInput$Temp = plotDataStep[as.character(spDataInput$ID),"Temp"]
spDataInput$PET = plotDataStep[as.character(spDataInput$ID),"PET"]
spDataInput$Prec = plotDataStep[as.character(spDataInput$ID),"Prec"]
spDataInput$PPET = pmin(spDataInput$Prec/spDataInput$PET, 1.3)
spDataInput$Temp15sq = (spDataInput$Temp-15)^2
# Substep 1: Growth
INC = .IFNgrowthInternal(dataInput, numYears = numYears, response = "DI",
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose=verbose)
if(heightgrowth) {
HINC = .IFNheightgrowthInternal(dataInput, numYears = numYears,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose=verbose)
newH = dataInput$h + HINC
}
# Substep 2: Survival
Psurv = .IFNsurvivalInternal(dataInput, numYears = numYears, killVeryLargeTrees = killVeryLargeTrees, verbose = verbose)
if(length(Psurv)!=nrow(dataInput)) stop("Unequal size for Psurv and dataInput")
if(pastGrowthMortalityEffect) {
if(nrow(dataInput)!=length(DBHincprev)) stop("Unequal size for DBHincprev and dataInput")
dataInput$ln.d.inc.prev = log(DBHincprev+1)
dataInput$ln.d.inc.prev.ln.d = dataInput$ln.d.inc.prev/dataInput$ln.d
selnoNA = !is.na(DBHincprev)
if(sum(selnoNA)>0) {
PsurvPG = .IFNsurvivalPGInternal(dataInput[selnoNA,,drop=FALSE], numYears = numYears, killVeryLargeTrees = killVeryLargeTrees, verbose = verbose)
if(length(PsurvPG)!=sum(selnoNA)) stop("Unequal size for PsurvPG and selnoNA")
if(sum(!is.na(PsurvPG))>0) Psurv[selnoNA][!is.na(PsurvPG)] = PsurvPG[!is.na(PsurvPG)]
}
}
if(!applyMortalityInManagedStands) {
if("BAL.ext" %in% names(dataInput)) {
ids_man = unique(dataInput$ID[dataInput$BAL.ext>0])
Psurv[dataInput$ID %in% ids_man] = 1
}
}
if(fire) { ## If stand burn then survival is not density-dependent but depends on fire
fireOutput = IFNfire(treeData, plotDataStep)
treePsurvfire = fireOutput$tree$Psurvfire
# if(sum(!is.na(treePsurvfire))>0) print(cbind(Psurv,treePsurvfire))
# Multiply survival probabilities for those stands that burned
Psurv[!is.na(treePsurvfire)] = Psurv[!is.na(treePsurvfire)]*treePsurvfire[!is.na(treePsurvfire)]
}
# Substep 3: Apply growth and survival
Nnew = treeData$N
if(stochastic) {
# Apply a binomial survival random variable
for(j in 1:length(Nnew)) Nnew[j] = min(Nnew[j], rbinom(1, prob = Psurv[j], size=round(Nnew[j])))
} else {
# Multiply survival probability by previous density estimate
Nnew = treeData$N*Psurv
}
treeDataOutStep = data.frame(Run = rep(r, length(treeData$ID)),
Step = rep(s, length(treeData$ID)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH + INC,
N = Nnew, stringsAsFactors = FALSE)
if(pastGrowthMortalityEffect) treeDataOutStep$DBHincprev = INC
if(BALext && removeCutTrees) { # If BALext was supplied remove cut trees
if(verbose) {
cat("\n Removing cut trees...\n")
}
if(heightgrowth) newH = newH[treeData$OIFFIN!="0",]
treeDataOutStep = treeDataOutStep[treeData$OIFFIN!="0",]
}
# Apply variable radius correction (if necessary)
if(variable.radius) {
dfr = .densityFactor(treeDataOutStep$DBH)/.densityFactor(treeData$DBH)
treeDataOutStep$N = treeDataOutStep$N*dfr
}
# Merge small trees (to avoid accumulation of records)
if(mergeTrees) {
treeDataOutStep = .mergeTrees(treeDataOutStep)
}
# Assumes dead trees did not grow in diameter
if(dead.table) {
treeDataDeadStep = data.frame(Run = rep(r, nrow(treeData)),
Step = rep(s, nrow(treeData)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH,
N = treeData$N - Nnew, stringsAsFactors = FALSE)
if(pastGrowthMortalityEffect) treeDataDeadStep$DBHincprev = INC
}
# Substep 4: Add ingrowth
if(ingrowth) {
#Local recruitment
col = (spDataInput$Nsp.N == 0) # Is new colonization?
recrData = .IFNingrowthInternal(spDataInput[!col,,drop=FALSE], numYears = numYears,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose=verbose)
if(sum(col)>0) { # If there is colonization apply models with dispersal for those records
recrData_disp = .IFNingrowthdispInternal(spDataInput[col,,drop=FALSE], numYears = numYears,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose=verbose)
ndisp = nrow(recrData_disp)
recrData = rbind(recrData, recrData_disp)
} else {
ndisp = 0
}
if(pastGrowthMortalityEffect) recrData$DBHincprev = NA
if(recr.labels) {
treeDataOutStep$Ingrowth = rep(FALSE,nrow(treeDataOutStep))
treeDataOutStep$Colonization = rep(FALSE,nrow(treeDataOutStep))
recrData$Run = rep(r, nrow(recrData))
recrData$Step = rep(s, nrow(recrData))
recrData$Ingrowth = rep(TRUE, nrow(recrData))
recrData$Colonization = c(rep(FALSE, nrow(recrData)-ndisp),rep(FALSE, ndisp))
}
if(fire) { ## Assume fire kills ingrowth
Psurv = rep(1, nrow(recrData))
Psurv[fireOutput$stand[recrData$ID,"fire"]==1] = 0
recrData$N = recrData$N*Psurv
recrData = recrData[recrData$N>0, , drop = FALSE]
# recrFire = .treeFire(recrData, fireOutput$stand)
# recrPsurvfire = recrFire$Psurvfire
# Psurv = rep(1, nrow(recrData))
# Psurv[!is.na(recrPsurvfire)] = recrPsurvfire[!is.na(recrPsurvfire)] # Replace survival probabilities for those stands that burned
# recrData$N = recrData$N*Psurv
}
treeDataOutStep = rbind(treeDataOutStep, recrData[,names(treeDataOutStep)])
}
# Substep 5: Height
if(height) {
treeDataOutStep$H = NA
#Replace height predictions derived from height increment model
if(heightgrowth) treeDataOutStep$H[1:length(newH)] = newH
hstatic = IFNheight(treeDataOutStep, plotDataStep,
useProvince = useProvince, useRP = useRP, verbose=verbose)
#Replace missing values with static predictions
naH = is.na(treeDataOutStep$H)
treeDataOutStep$H[naH] = hstatic[naH]
if(dead.table) {
treeDataDeadStep$H = IFNheight(treeDataDeadStep, plotDataStep,
useProvince = useProvince, useRP = useRP, verbose=verbose)
}
}
# Substep 6: Remove low density cohorts
treeDataOutStep = treeDataOutStep[treeDataOutStep$N>min.N, , drop = FALSE]
# For next step
if(s<numSteps) {
treeData = treeDataOutStep
if(pastGrowthMortalityEffect) DBHincprev = treeData$DBHincprev
}
if(sequence) {
treeDataSeq = rbind(treeDataSeq,treeDataOutStep[,names(treeDataSeq)])
}
if(dead.table) {
if(is.null(treeDataDeadRun)) {
treeDataDeadRun = treeDataDeadStep
} else {
treeDataDeadRun = rbind(treeDataDeadRun, treeDataDeadStep)
}
}
}
}
if(sequence) {
treeDataOutRun = treeDataSeq
} else {
treeDataOutRun = treeDataOutStep
}
if(!is.null(treeDataOut)) {
treeDataOut = rbind(treeDataOut, treeDataOutRun)
} else {
treeDataOut = treeDataOutRun
}
if(!is.null(treeDataDead)) {
treeDataDead = rbind(treeDataDead,treeDataDeadRun)
} else {
treeDataDead = treeDataDeadRun
}
}
# Species names
if(speciesNames) {
if(sequence) {
treeDataOut = data.frame(Run = treeDataOut$Run,
Step = treeDataOut$Step,
ID = treeDataOut$ID,
Species = treeDataOut$Species,
Name = speciesNames(treeDataOut$Species),
treeDataOut[,-c(1:4)], stringsAsFactors = FALSE)
} else {
treeDataOut = data.frame(Run = treeDataOut$Run,
ID = treeDataOut$ID,
Species = treeDataOut$Species,
Name = speciesNames(treeDataOut$Species),
treeDataOut[,-c(1:4)], stringsAsFactors = FALSE)
}
if(dead.table) {
treeDataDead = data.frame(Run = treeDataDead$Run,
Step = treeDataDead$Step,
ID = treeDataDead$ID,
Species = treeDataDead$Species,
Name = speciesNames(treeDataDead$Species),
treeDataDead[,-c(1:4)], stringsAsFactors = FALSE)
}
}
rownames(treeDataOut) = NULL
rownames(treeDataDead) = NULL
if(verbose) cat("\n")
if(!dead.table) {
return(treeDataOut)
}
return(list(out = treeDataOut, dead = treeDataDead))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.