R/IFNsubmodels.R
In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis
Defines functions
IFNfire
.treeFire
IFNheightgrowth
IFNheight
IFNingrowthdisp
IFNingrowth
IFNsurvivalPG
IFNsurvival
IFNgrowth
Documented in
IFNfire
IFNgrowth
IFNheight
IFNheightgrowth
IFNingrowth
IFNingrowthdisp
IFNsurvival
IFNsurvivalPG
#' Dynamic model components
#'
#' Set of basic model functions to simulate forest dynamics in IFN plots:
#' \itemize{
#' \item{\code{IFNgrowth} estimates growth as tree diameter increment or basal area increment.}
#' \item{\code{IFNsurvival} estimates survival probabilities.}
#' \item{\code{IFNsurvivalPG} estimates survival probabilities including previous diameter growth as predictor.}
#' \item{\code{IFNingrowth} estimates tree ingrowth densities and average diameters.}
#' \item{\code{IFNingrowthdisp} estimates tree ingrowth densities and average diameters taking into account spatial context.}
#' \item{\code{IFNheight} estimates height of trees (in m).}
#' \item{\code{IFNheightgrowth} estimates tree height increment (in m).}
#' \item{\code{IFNfire} estimates burning probability at the stand level and probability of death at the tree level.}
#' }
#'
#' @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' in rows and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C) and 'slope' (degrees).
#' @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 params Name of parameter data set or data frame with parameter values per species
#' @param response Type of output: "BAI" (basal area increment, cm2), "DI" (diameter increment, cm) or "DBH" (new diameter at breast height, cm)
#' @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 stochastic A flag to indicate that simulation of growth, mortality and ingrowth should be stochastic (i.e. non-deterministic)
#' @param prepared A flag to indicate that \code{treeData} is the result of calling function \code{\link{prepareTreeCompetitionTable}}
#' (for \code{IFNgrowth} and \code{IFNsurvival}), \code{\link{prepareSpeciesCompetitionTable}}
#' (for \code{IFNingrowth}) or \code{\link{prepareSpatialSpeciesCompetitionTable}}
#' (for \code{IFNingrowthdisp}). Otherwise these functions are called internally.
#' @param verbose A flag to indicate console output of the process
#'
#' @return Depends on the function:
#' \itemize{
#' \item{\code{IFNgrowth} returns a vector with basal area increment (in cm2), diameter increments (in cm) or new diameter values (in cm).}
#' \item{\code{IFNsurvival} and \code{IFNsurvivalPG} return a vector with survival probabilities (0-1).}
#' \item{\code{IFNingrowth} and \code{IFNingrowthdisp} return a data frame with new tree records in rows and columns 'ID', 'Species', ('Name'), 'N' and 'DBH'.}
#' \item{\code{IFNheight} return a vector with height values (in m).}
#' \item{\code{IFNheightgrowth} return a vector with height increment values (in m).}
#' \item{\code{IFNfire} returns a list with two elements: (1) 'stand' is a data frame with as many rows as plots and columns 'Pfire10' (probability of burning in 10 year interval), 'fire' stochastic trial of fire occurrence and 'Pdead' proportion of trees dead;
#' (2) 'tree' is a data frame with columns 'Psurv' (the probability of survival for each tree in each stand, following the order of 'treeData') and 'Psurv' (the probability of survival for each tree in each stand that burned in the stochastic trial, otherwise NA). }
#' }
#'
#' @details All functions check 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). Stochasticity is not available for \code{IFNsurvival} and
#' \code{IFNsurvivalPG} because these function return survival probabilities. Function \code{IFNheight} is only deterministic, because it is a static model.
#' Function \code{IFNfire} implements the fire probability model of González et al. (2006) and the fire damage model of González et al. (2007), both calibrated for
#' Catalonia. Height growht estimates of function \code{IFNheightgrowth} return \code{NA} values whenever the province is not included in the set of
#' provinces with enough data to define height growth models for the target species.
#'
#' @name IFNsubmodels
#' @aliases IFNgrowth IFNsurvival IFNingrowth IFNingrowthdisp IFNheight IFNheightgrowth IFNfire
#' @seealso \code{\link{IFNdynamics}}, \code{\link{prepareTreeCompetitionTable}}
#'
#' @references
#' González, J.R., Trasobares, A., Palahí, M., Pukkala, T., 2007. Predicting stand damage and tree survival in burned forests in Catalonia (North-East Spain). Ann. For. Sci. 64, 733–742. https://doi.org/10.1051/forest
#'
#' González, J.R., Palahí, M., Trasobares, A., Pukkala, T., 2006. A fire probability model for forest stands in Catalonia (north-east Spain). Ann. For. Sci. 63, 169–176. https://doi.org/10.1051/forest:2005109
#'
#' @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)
#'
#' # (1) Growth
#' IFNgrowth(exampleTreeData, plotDataIFN)
#'
#' # (2) Survival probabilities
#' IFNsurvival(exampleTreeData, plotDataIFN)
#'
#' # (3) Ingrowth model (local recruiment)
#' IFNingrowth(exampleTreeData, plotDataIFN)
#'
#' # (4) Ingrowth model (colonization)
#' 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"
#' IFNingrowthdisp(exampleTreeData, plotDataIFN, nb, NspNmatrix)
#'
#' # (5) Static height model
#' IFNheight(exampleTreeData, plotDataIFN)
#'
#'# (6) Height increment model
#' IFNheightgrowth(exampleTreeData, plotDataIFN)
#'
#' # (7) Stochastic fire model
#' IFNfire(exampleTreeData, plotDataIFN)
IFNgrowth<-function(treeData, plotData, numYears = 10,
params = "defaultGrowthParams", response = "DI", prepared = FALSE,
useProvince = TRUE, useRP = TRUE, stochastic = FALSE, verbose = FALSE) {
response = match.arg(response, c("BAI", "DI", "DBH"))
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs and species codes
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if(useProvince) if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(treeData$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
treeData[[rp]] = plotData[as.character(treeData$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
treeData[spns!=sp,rpmap] = NA
}
}
}
# print(head(treeData))
# print(apply(treeData[,RP_vars],2, function(x){sum(!is.na(x))}))
INC = .IFNgrowthInternal(treeData, numYears, params, response = response,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose = verbose)
if(verbose) cat("\n")
return(INC)
}
#' @rdname IFNsubmodels
#' @param killVeryLargeTrees Flag to indicate that trees with diameters larger than species-specific thresholds should be forced to die.
IFNsurvival<-function(treeData, plotData, numYears = 10,
params = "defaultSurvivalParams", prepared = FALSE,
killVeryLargeTrees = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
Psurv = .IFNsurvivalInternal(treeData, numYears, params, killVeryLargeTrees, verbose)
if(verbose) cat("\n")
return(Psurv)
}
#' @rdname IFNsubmodels
IFNsurvivalPG<-function(treeData, plotData, numYears = 10,
params = "defaultSurvivalPGParams", prepared = FALSE,
killVeryLargeTrees = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!("DBHincprev" %in% names(treeData))) stop("Column 'DBHincprev' needs to be defined in treeData.")
DBHincprev = treeData$DBHincprev
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
treeData$ln.d.inc.prev = log(DBHincprev+1)
treeData$ln.d.inc.prev.ln.d = treeData$ln.d.inc.prev/treeData$ln.d
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
Psurv = .IFNsurvivalPGInternal(treeData, numYears, params, killVeryLargeTrees, verbose)
if(verbose) cat("\n")
return(Psurv)
}
#' @rdname IFNsubmodels
#' @param paramsB Name of binomial model parameter data set or data frame with parameter values per species
#' @param paramsN Name of Gamma model parameter data set or data frame with parameter values per species
#' @param paramsDBH Name of Gamma model parameter data set or data frame with parameter values per species
#' @param speciesNames A flag to include species names in the output
#' @param thresholdIngrowth Boolean flag to indicate that ingrowth should occur only when a given probability threshold is surpassed.
IFNingrowth<-function(treeData, plotData, numYears = 10,
paramsB = "defaultIngrowthBParams",
paramsN = "defaultIngrowthNParams",
paramsDBH = "defaultIngrowthDBHParams",
prepared = FALSE, stochastic = FALSE,
thresholdIngrowth = FALSE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareSpeciesCompetitionTable(treeData,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
ing = .IFNingrowthInternal(treeData, numYears, paramsB, paramsN, paramsDBH,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose = verbose)
if(speciesNames) ing = data.frame(ing[,1:2], Name = speciesNames(ing$Species), ing[,3:4])
if(verbose) cat("\n")
return(ing)
}
#' @rdname IFNsubmodels
#' @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).
IFNingrowthdisp<-function(treeData, plotData,
nb, NspNmatrix, numYears = 10,
paramsB = "defaultIngrowthBParams_disp",
paramsN = "defaultIngrowthNParams_disp",
paramsDBH = "defaultIngrowthDBHParams_disp",
prepared = FALSE, stochastic = FALSE,
thresholdIngrowth = FALSE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareSpatialSpeciesCompetitionTable(treeData, nb = nb, NspNmatrix = NspNmatrix,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
ing = .IFNingrowthdispInternal(treeData,
numYears, paramsB, paramsN, paramsDBH,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose = verbose)
if(speciesNames) ing = data.frame(ing[,1:2], Name = speciesNames(ing$Species), ing[,3:4])
if(verbose) cat("\n")
return(ing[treeData$Nsp.N==0,, drop =FALSE])
}
#' @rdname IFNsubmodels
IFNheight<-function(treeData, plotData,
params = "defaultHeightParams",
useProvince = TRUE, useRP = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
H = numeric(0)
if(nrow(treeData)>0) {
dataInput = treeData
dataInput$Province = .getProvinceFromID(dataInput$ID)
if(useProvince) if("Province" %in% names(plotData)) dataInput$Province = plotData[as.character(dataInput$ID),"Province"]
dataInput$SWHC = plotData[as.character(dataInput$ID),"SWHC"]
dataInput$Rad = plotData[as.character(dataInput$ID),"Rad"]
dataInput$PET = plotData[as.character(dataInput$ID),"PET"]
dataInput$Prec = plotData[as.character(dataInput$ID),"Prec"]
dataInput$PPET = pmin(dataInput$Prec/dataInput$PET, 1.3)
dataInput$elevation = plotData[as.character(dataInput$ID),"elevation"]
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(dataInput$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
dataInput[[rp]] = plotData[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
}
}
}
H = .IFNheightInternal(dataInput, params, useProvince, useRP, verbose = verbose)
}
if(verbose) cat("\n")
return(H)
}
#' @rdname IFNsubmodels
IFNheightgrowth<-function(treeData, plotData, numYears = 10,
params = "defaultHeightGrowthParams", prepared = FALSE,
useProvince = TRUE, useRP = TRUE, stochastic = FALSE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs and species codes
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if(useProvince) if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(treeData$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
treeData[[rp]] = plotData[as.character(treeData$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
treeData[spns!=sp,rpmap] = NA
}
}
}
# print(head(treeData))
# print(apply(treeData[,RP_vars],2, function(x){sum(!is.na(x))}))
INC = .IFNheightgrowthInternal(treeData, numYears, params,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose = verbose)
if(verbose) cat("\n")
return(INC)
}
.treeFire<-function(treeData, standFire){
treeFire = data.frame(logitsurv = 2.224 + 0.110*treeData$DBH - 7.117*standFire[treeData$ID,"Pdead"])
treeFire$Psurv = 1/(1+exp(-treeFire$logitsurv))
treeFire$Psurvfire = treeFire$Psurv
treeFire$Psurvfire[standFire[treeData$ID,"fire"]==0] = NA
return(treeFire)
}
#' @rdname IFNsubmodels
IFNfire<-function(treeData, plotData) {
tBA = treeData$N*pi*(treeData$DBH/200)^2
fg = functionalGroups(treeData$Species)
tBAc = tBA*ifelse(fg=="conifer",1,0)
df = data.frame(BA = basalArea(treeData))
ids = row.names(df)
BAconifer = tapply(tBAc, treeData$ID, sum, na.rm=T)
df$BAconifer = BAconifer[ids]
df$BAhardwood = df$BA - df$BAconifer
N = tapply(treeData$N, treeData$ID, sum, na.rm=T)
df$N = N[ids]
Dg = tapply(treeData$DBH*tBA, treeData$ID, sum, na.rm=T)
df$Dg = Dg[ids]
df$Dg = df$Dg / df$BA
df$QMD = sqrt(df$BA/(0.0000785*df$N))
df$Sd = diameterSD(treeData)
df$elev = log(pmax(1,(plotData[ids, "elevation"]/100) - 6, na.rm=T))
df$slopeperc = 100*sin(plotData[ids, "slope"]*(pi/180))
df$Phard = df$BAhardwood/df$BA # Proportion in terms of BA
df$Pine = ifelse(df$BAconifer/df$BA > 0.5, 1, 0) # Is >50% of pines in terms of BA
df$logit12 = -1.925-2.256*df$elev-0.015*df$Dg+0.012*df$BA-1.763*df$Phard+2.081*(df$Sd/(df$Dg+0.01))
df$Pfire12 = 1/(1+exp(-df$logit12))
df$Pfire10 = 1-(1-df$Pfire12)^(10/12)
df$logitdead = -6.131 -0.329*df$BA + 0.060*df$slopeperc + 2.266*df$Pine + 4.319*(df$BA/(df$QMD+0.01)) + 6.718*(df$Sd/(df$QMD+0.01))
df$Pdead = 1/(1+exp(-df$logitdead))
# print(df$Pfire10)
df$fire = rbinom(nrow(df), 1, df$Pfire10)
# print(df)
treeFire = .treeFire(treeData, df)
return(list(stand = df[,c("Pfire10", "fire","Pdead")], tree = treeFire[,c("Psurv", "Psurvfire")]))
}
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 IFNfire .treeFire IFNheightgrowth IFNheight IFNingrowthdisp IFNingrowth IFNsurvivalPG IFNsurvival IFNgrowth
Documented in IFNfire IFNgrowth IFNheight IFNheightgrowth IFNingrowth IFNingrowthdisp IFNsurvival IFNsurvivalPG
#' Dynamic model components
#'
#' Set of basic model functions to simulate forest dynamics in IFN plots:
#' \itemize{
#' \item{\code{IFNgrowth} estimates growth as tree diameter increment or basal area increment.}
#' \item{\code{IFNsurvival} estimates survival probabilities.}
#' \item{\code{IFNsurvivalPG} estimates survival probabilities including previous diameter growth as predictor.}
#' \item{\code{IFNingrowth} estimates tree ingrowth densities and average diameters.}
#' \item{\code{IFNingrowthdisp} estimates tree ingrowth densities and average diameters taking into account spatial context.}
#' \item{\code{IFNheight} estimates height of trees (in m).}
#' \item{\code{IFNheightgrowth} estimates tree height increment (in m).}
#' \item{\code{IFNfire} estimates burning probability at the stand level and probability of death at the tree level.}
#' }
#'
#' @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' in rows and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C) and 'slope' (degrees).
#' @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 params Name of parameter data set or data frame with parameter values per species
#' @param response Type of output: "BAI" (basal area increment, cm2), "DI" (diameter increment, cm) or "DBH" (new diameter at breast height, cm)
#' @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 stochastic A flag to indicate that simulation of growth, mortality and ingrowth should be stochastic (i.e. non-deterministic)
#' @param prepared A flag to indicate that \code{treeData} is the result of calling function \code{\link{prepareTreeCompetitionTable}}
#' (for \code{IFNgrowth} and \code{IFNsurvival}), \code{\link{prepareSpeciesCompetitionTable}}
#' (for \code{IFNingrowth}) or \code{\link{prepareSpatialSpeciesCompetitionTable}}
#' (for \code{IFNingrowthdisp}). Otherwise these functions are called internally.
#' @param verbose A flag to indicate console output of the process
#'
#' @return Depends on the function:
#' \itemize{
#' \item{\code{IFNgrowth} returns a vector with basal area increment (in cm2), diameter increments (in cm) or new diameter values (in cm).}
#' \item{\code{IFNsurvival} and \code{IFNsurvivalPG} return a vector with survival probabilities (0-1).}
#' \item{\code{IFNingrowth} and \code{IFNingrowthdisp} return a data frame with new tree records in rows and columns 'ID', 'Species', ('Name'), 'N' and 'DBH'.}
#' \item{\code{IFNheight} return a vector with height values (in m).}
#' \item{\code{IFNheightgrowth} return a vector with height increment values (in m).}
#' \item{\code{IFNfire} returns a list with two elements: (1) 'stand' is a data frame with as many rows as plots and columns 'Pfire10' (probability of burning in 10 year interval), 'fire' stochastic trial of fire occurrence and 'Pdead' proportion of trees dead;
#' (2) 'tree' is a data frame with columns 'Psurv' (the probability of survival for each tree in each stand, following the order of 'treeData') and 'Psurv' (the probability of survival for each tree in each stand that burned in the stochastic trial, otherwise NA). }
#' }
#'
#' @details All functions check 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). Stochasticity is not available for \code{IFNsurvival} and
#' \code{IFNsurvivalPG} because these function return survival probabilities. Function \code{IFNheight} is only deterministic, because it is a static model.
#' Function \code{IFNfire} implements the fire probability model of González et al. (2006) and the fire damage model of González et al. (2007), both calibrated for
#' Catalonia. Height growht estimates of function \code{IFNheightgrowth} return \code{NA} values whenever the province is not included in the set of
#' provinces with enough data to define height growth models for the target species.
#'
#' @name IFNsubmodels
#' @aliases IFNgrowth IFNsurvival IFNingrowth IFNingrowthdisp IFNheight IFNheightgrowth IFNfire
#' @seealso \code{\link{IFNdynamics}}, \code{\link{prepareTreeCompetitionTable}}
#'
#' @references
#' González, J.R., Trasobares, A., Palahí, M., Pukkala, T., 2007. Predicting stand damage and tree survival in burned forests in Catalonia (North-East Spain). Ann. For. Sci. 64, 733–742. https://doi.org/10.1051/forest
#'
#' González, J.R., Palahí, M., Trasobares, A., Pukkala, T., 2006. A fire probability model for forest stands in Catalonia (north-east Spain). Ann. For. Sci. 63, 169–176. https://doi.org/10.1051/forest:2005109
#'
#' @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)
#'
#' # (1) Growth
#' IFNgrowth(exampleTreeData, plotDataIFN)
#'
#' # (2) Survival probabilities
#' IFNsurvival(exampleTreeData, plotDataIFN)
#'
#' # (3) Ingrowth model (local recruiment)
#' IFNingrowth(exampleTreeData, plotDataIFN)
#'
#' # (4) Ingrowth model (colonization)
#' 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"
#' IFNingrowthdisp(exampleTreeData, plotDataIFN, nb, NspNmatrix)
#'
#' # (5) Static height model
#' IFNheight(exampleTreeData, plotDataIFN)
#'
#'# (6) Height increment model
#' IFNheightgrowth(exampleTreeData, plotDataIFN)
#'
#' # (7) Stochastic fire model
#' IFNfire(exampleTreeData, plotDataIFN)
IFNgrowth<-function(treeData, plotData, numYears = 10,
params = "defaultGrowthParams", response = "DI", prepared = FALSE,
useProvince = TRUE, useRP = TRUE, stochastic = FALSE, verbose = FALSE) {
response = match.arg(response, c("BAI", "DI", "DBH"))
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs and species codes
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if(useProvince) if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(treeData$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
treeData[[rp]] = plotData[as.character(treeData$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
treeData[spns!=sp,rpmap] = NA
}
}
}
# print(head(treeData))
# print(apply(treeData[,RP_vars],2, function(x){sum(!is.na(x))}))
INC = .IFNgrowthInternal(treeData, numYears, params, response = response,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose = verbose)
if(verbose) cat("\n")
return(INC)
}
#' @rdname IFNsubmodels
#' @param killVeryLargeTrees Flag to indicate that trees with diameters larger than species-specific thresholds should be forced to die.
IFNsurvival<-function(treeData, plotData, numYears = 10,
params = "defaultSurvivalParams", prepared = FALSE,
killVeryLargeTrees = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
Psurv = .IFNsurvivalInternal(treeData, numYears, params, killVeryLargeTrees, verbose)
if(verbose) cat("\n")
return(Psurv)
}
#' @rdname IFNsubmodels
IFNsurvivalPG<-function(treeData, plotData, numYears = 10,
params = "defaultSurvivalPGParams", prepared = FALSE,
killVeryLargeTrees = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!("DBHincprev" %in% names(treeData))) stop("Column 'DBHincprev' needs to be defined in treeData.")
DBHincprev = treeData$DBHincprev
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
treeData$ln.d.inc.prev = log(DBHincprev+1)
treeData$ln.d.inc.prev.ln.d = treeData$ln.d.inc.prev/treeData$ln.d
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
Psurv = .IFNsurvivalPGInternal(treeData, numYears, params, killVeryLargeTrees, verbose)
if(verbose) cat("\n")
return(Psurv)
}
#' @rdname IFNsubmodels
#' @param paramsB Name of binomial model parameter data set or data frame with parameter values per species
#' @param paramsN Name of Gamma model parameter data set or data frame with parameter values per species
#' @param paramsDBH Name of Gamma model parameter data set or data frame with parameter values per species
#' @param speciesNames A flag to include species names in the output
#' @param thresholdIngrowth Boolean flag to indicate that ingrowth should occur only when a given probability threshold is surpassed.
IFNingrowth<-function(treeData, plotData, numYears = 10,
paramsB = "defaultIngrowthBParams",
paramsN = "defaultIngrowthNParams",
paramsDBH = "defaultIngrowthDBHParams",
prepared = FALSE, stochastic = FALSE,
thresholdIngrowth = FALSE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareSpeciesCompetitionTable(treeData,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
ing = .IFNingrowthInternal(treeData, numYears, paramsB, paramsN, paramsDBH,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose = verbose)
if(speciesNames) ing = data.frame(ing[,1:2], Name = speciesNames(ing$Species), ing[,3:4])
if(verbose) cat("\n")
return(ing)
}
#' @rdname IFNsubmodels
#' @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).
IFNingrowthdisp<-function(treeData, plotData,
nb, NspNmatrix, numYears = 10,
paramsB = "defaultIngrowthBParams_disp",
paramsN = "defaultIngrowthNParams_disp",
paramsDBH = "defaultIngrowthDBHParams_disp",
prepared = FALSE, stochastic = FALSE,
thresholdIngrowth = FALSE,
speciesNames = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareSpatialSpeciesCompetitionTable(treeData, nb = nb, NspNmatrix = NspNmatrix,
provinceFromID = TRUE, verbose = verbose)
}
if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
ing = .IFNingrowthdispInternal(treeData,
numYears, paramsB, paramsN, paramsDBH,
thresholdIngrowth = thresholdIngrowth,
stochastic = stochastic, verbose = verbose)
if(speciesNames) ing = data.frame(ing[,1:2], Name = speciesNames(ing$Species), ing[,3:4])
if(verbose) cat("\n")
return(ing[treeData$Nsp.N==0,, drop =FALSE])
}
#' @rdname IFNsubmodels
IFNheight<-function(treeData, plotData,
params = "defaultHeightParams",
useProvince = TRUE, useRP = TRUE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
H = numeric(0)
if(nrow(treeData)>0) {
dataInput = treeData
dataInput$Province = .getProvinceFromID(dataInput$ID)
if(useProvince) if("Province" %in% names(plotData)) dataInput$Province = plotData[as.character(dataInput$ID),"Province"]
dataInput$SWHC = plotData[as.character(dataInput$ID),"SWHC"]
dataInput$Rad = plotData[as.character(dataInput$ID),"Rad"]
dataInput$PET = plotData[as.character(dataInput$ID),"PET"]
dataInput$Prec = plotData[as.character(dataInput$ID),"Prec"]
dataInput$PPET = pmin(dataInput$Prec/dataInput$PET, 1.3)
dataInput$elevation = plotData[as.character(dataInput$ID),"elevation"]
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(dataInput$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
dataInput[[rp]] = plotData[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
}
}
}
H = .IFNheightInternal(dataInput, params, useProvince, useRP, verbose = verbose)
}
if(verbose) cat("\n")
return(H)
}
#' @rdname IFNsubmodels
IFNheightgrowth<-function(treeData, plotData, numYears = 10,
params = "defaultHeightGrowthParams", prepared = FALSE,
useProvince = TRUE, useRP = TRUE, stochastic = FALSE, verbose = FALSE) {
## Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
## Check IDs and species codes
.checkIDs(treeData, plotData)
.checkSpeciesCodes(treeData)
if(!prepared) {
if(verbose) {
cat("\n Preparing data...\n")
}
treeData = prepareTreeCompetitionTable(treeData, BAL = TRUE, BALext = FALSE, sortByDBH = FALSE,
provinceFromID = TRUE, verbose = verbose)
}
if(useProvince) if("Province" %in% names(plotData)) treeData$Province = plotData[as.character(treeData$ID),"Province"]
treeData$SWHC = plotData[as.character(treeData$ID),"SWHC"]
treeData$Rad = plotData[as.character(treeData$ID),"Rad"]
treeData$Temp = plotData[as.character(treeData$ID),"Temp"]
treeData$PET = plotData[as.character(treeData$ID),"PET"]
treeData$Prec = plotData[as.character(treeData$ID),"Prec"]
treeData$slope = plotData[as.character(treeData$ID),"slope"]
treeData$PPET = pmin(treeData$Prec/treeData$PET, 1.3)
treeData$Temp15sq = (treeData$Temp-15)^2
if(useRP) {
RP_vars = .getVarsRP()
spns = speciesNamesModels(treeData$Species)
RP_mapping = get("RP_mapping")
for(rp in RP_vars) {
treeData[[rp]] = plotData[as.character(treeData$ID),rp]
}
for(i in 1:length(RP_mapping)) {
sp = names(RP_mapping)[i]
rpmap = RP_mapping[[i]]
if(!is.na(rpmap[1])) {
treeData[spns!=sp,rpmap] = NA
}
}
}
# print(head(treeData))
# print(apply(treeData[,RP_vars],2, function(x){sum(!is.na(x))}))
INC = .IFNheightgrowthInternal(treeData, numYears, params,
useProvince = useProvince, useRP = useRP,
stochastic = stochastic, verbose = verbose)
if(verbose) cat("\n")
return(INC)
}
.treeFire<-function(treeData, standFire){
treeFire = data.frame(logitsurv = 2.224 + 0.110*treeData$DBH - 7.117*standFire[treeData$ID,"Pdead"])
treeFire$Psurv = 1/(1+exp(-treeFire$logitsurv))
treeFire$Psurvfire = treeFire$Psurv
treeFire$Psurvfire[standFire[treeData$ID,"fire"]==0] = NA
return(treeFire)
}
#' @rdname IFNsubmodels
IFNfire<-function(treeData, plotData) {
tBA = treeData$N*pi*(treeData$DBH/200)^2
fg = functionalGroups(treeData$Species)
tBAc = tBA*ifelse(fg=="conifer",1,0)
df = data.frame(BA = basalArea(treeData))
ids = row.names(df)
BAconifer = tapply(tBAc, treeData$ID, sum, na.rm=T)
df$BAconifer = BAconifer[ids]
df$BAhardwood = df$BA - df$BAconifer
N = tapply(treeData$N, treeData$ID, sum, na.rm=T)
df$N = N[ids]
Dg = tapply(treeData$DBH*tBA, treeData$ID, sum, na.rm=T)
df$Dg = Dg[ids]
df$Dg = df$Dg / df$BA
df$QMD = sqrt(df$BA/(0.0000785*df$N))
df$Sd = diameterSD(treeData)
df$elev = log(pmax(1,(plotData[ids, "elevation"]/100) - 6, na.rm=T))
df$slopeperc = 100*sin(plotData[ids, "slope"]*(pi/180))
df$Phard = df$BAhardwood/df$BA # Proportion in terms of BA
df$Pine = ifelse(df$BAconifer/df$BA > 0.5, 1, 0) # Is >50% of pines in terms of BA
df$logit12 = -1.925-2.256*df$elev-0.015*df$Dg+0.012*df$BA-1.763*df$Phard+2.081*(df$Sd/(df$Dg+0.01))
df$Pfire12 = 1/(1+exp(-df$logit12))
df$Pfire10 = 1-(1-df$Pfire12)^(10/12)
df$logitdead = -6.131 -0.329*df$BA + 0.060*df$slopeperc + 2.266*df$Pine + 4.319*(df$BA/(df$QMD+0.01)) + 6.718*(df$Sd/(df$QMD+0.01))
df$Pdead = 1/(1+exp(-df$logitdead))
# print(df$Pfire10)
df$fire = rbinom(nrow(df), 1, df$Pfire10)
# print(df)
treeFire = .treeFire(treeData, df)
return(list(stand = df[,c("Pfire10", "fire","Pdead")], tree = treeFire[,c("Psurv", "Psurvfire")]))
}
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.