#' Example forest dynamics by species
#'
#' Simulates stand dynamics and plots the temporal series of a stand-level variable.
#'
#' @param species Character vector of species codes to be studied.
#' @param prov A character or integer of Spanish province between 1 and 50. Rellevant for regionalized species models.
#' @param initialDBH Numeric vector with initial diameter values (in cm). Initial density associated to each tree record will depend on diameter.
#' @param numSteps Number of 10-yr steps to be simulated.
#' @param SWHC Soil water holding capacity (mm).
#' @param Rad Average daily radiation (MJ).
#' @param Temp Mean annual temperature (degrees C).
#' @param Prec Mean annual precipitation (mm)
#' @param PET Mean annual potential evapotranspiration (mm)
#' @param slope Slope in degrees.
#' @param elevation Elevation (m.a.s.l.).
#' @param testClimateChange Boolean flag to include simulations with altered climate
#'
#' @return A list with two elements:
#' \itemize{
#' \item{\code{standDynamics}: A list with one data frame per species with the complete stand dynamics (tree data per step).}
#' \item{\code{finalStates}: A list with one data frame per species with the final tree data at the end of the simulation.}
#' }
#'
#' @examples
#'
#' pines = c("21", "22", "23", "24", "25", "26", "27","28","31","17")
#' oaks = c("45", "44", "41", "43", "71","72")
#' other = c("51", "55","60", "73", "3","37", "83", "80")
#' sdyn = exampleStandDynamics(species=pines, numSteps = 40)
#'
#'
exampleStandDynamics<-function(species, prov = 8, initialDBH = NULL, numSteps = 20,
SWHC = 150, Rad = 25, Temp= 13, Prec = 900, PET = 1000,
slope = 0, elevation = 100,
testClimateChange = FALSE, ...) {
params = get("defaultGrowthParams")
sp_vec = as.character(params$species)
s = strsplit(sp_vec, split=",")
getRow<-function(sp) {
for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(i)
stop("Species not found")
}
examplePlots = get("df_list")
standDynamics = vector("list", length(species))
finalStates = vector("list", length(species))
res<-NULL
for(i in 1:length(species)) {
spi = species[i]
names(standDynamics)[i] = speciesNamesModels(spi)
names(finalStates)[i] = speciesNamesModels(spi)
if(!is.null(initialDBH)) {
initialDBHsp = initialDBH
treeData = data.frame(ID=as.character(as.numeric(prov)*10000), Species=spi,
DBH = initialDBHsp,
N=.densityFactor(initialDBHsp), stringsAsFactors = F)
treeData$H = IFNheight(treeData, climateData)
climateData = data.frame(SWHC = SWHC,
Rad = Rad,
Temp= Temp,
Prec = Prec,
PET = PET,
slope = slope,
elevation = elevation,
row.names=as.character(as.numeric(prov)*10000))
} else {
spn = as.character(params$names[getRow(spi)])
dfsp = examplePlots[[spn]]
treeData = dfsp
climateData = plotDataIFN[unique(treeData$ID), c("SWHC", "Rad", "Temp", "Prec", "PET", "slope", "elevation", "Province"), drop = FALSE]
}
if(testClimateChange) {
td2 = treeData
td2$ID = paste0(td2$ID, "_CC")
cd2 = climateData
row.names(cd2) = paste0(row.names(cd2), "_CC")
cd2$Prec = cd2$Prec*0.75
cd2$Temp = cd2$Temp+2.0
cd2$PET = cd2$PET*1.25
climateData = rbind(climateData, cd2)
treeData = rbind(treeData, td2)
}
res_sp = IFNdynamics(treeData, climateData, numSteps=numSteps, sequence = TRUE, ...)
finalStates[[i]] = res_sp[res_sp$Step==numSteps, ]
standDynamics[[i]] = res_sp
}
return(list(standDynamics = standDynamics, finalStates = finalStates))
}
#' Draws stand dynamics
#'
#' Draw temporal changes of stand-level variables.
#'
#' @param x the output of function \code{IFNdynamics}.
#' @param standVariable Variable to be plotted:
#' \itemize{
#' \item{'BA' - Basal area}
#' \item{'N' - Stem density}
#' \item{'QMD' - Quadratic mean diameter}
#' \item{'MAXD' - Maximum diameter}
#' \item{'HDOM' - Dominant height}
#' \item{'DDOM' - Dominant diameter}
#' \item{'DSD' - Standard deviation of diameter}
#' \item{'VCC' - Sum of stem volume including bark}
#' \item{'Biomass' - Sum of biomass of stem+branches+roots}
#' }
#' @param draw Logical flag to indicate that draw is produced
#'
#' @return If \code{draw = TRUE} returns a plot. A data frame with stand statistics per simulation step and run.
#'
#' @examples
#'
#' # Load example tree data (two forest plots)
#' data(exampleTreeData)
#'
#' # Load IFN environmental data
#' data(plotDataIFN)
#'
#' # Simulate stand dynamics and store transient states
#' x = IFNdynamics(exampleTreeData, plotDataIFN, numSteps = 10, sequence = TRUE)
#'
#' # Draw basal area dynamics
#' standDynamicsSummary(x)
standDynamicsSummary<-function(x, standVariable="BA", draw = TRUE) {
standVariable = match.arg(standVariable, c("BA","N","QMD", "DSD","MAXD","HDOM","DDOM", "VCC", "Biomass"))
if(inherits(x, "list")) {
treeData = x$out
} else {
treeData = x
}
if(!("Run" %in% names(treeData))) stop("Variable 'Run' must be defined for forest dynamics in 'x'")
if(!("Step" %in% names(treeData))) stop("Variable 'Step' must be defined for forest dynamics in 'x'")
if(!("ID" %in% names(treeData))) stop("Variable 'ID' must be defined for forest dynamics in 'x'")
steps = sort(as.numeric(unique(treeData$Step)))
runs = sort(as.numeric(unique(treeData$Run)))
IDs = unique(as.character(treeData$ID))
res = NULL
for(r in runs) {
treeDataRun = treeData[(treeData$Run==r), ]
for(s in steps) {
treeDataStep = treeDataRun[(treeDataRun$Step==s), ]
res_rs<-data.frame(Run = rep(r, length(IDs)),
Step = rep(s, length(IDs)),
Plot = IDs,
Year = rep(10*s, length(IDs)),
BA = rep(NA, length(IDs)),
N = rep(NA, length(IDs)),
DSD = rep(NA, length(IDs)),
QMD = rep(NA, length(IDs)),
MAXD = rep(NA, length(IDs)),
HDOM = rep(NA, length(IDs)),
DDOM = rep(NA, length(IDs)),
VCC = rep(NA, length(IDs)),
Biomass = rep(NA, length(IDs)))
for(i in 1:length(IDs)) {
id = IDs[i]
treeDataStepID = treeDataStep[(treeDataStep$ID==id), ]
if(nrow(treeDataStepID)>0) {
if((!draw) || (standVariable %in% c("BA", "QMD"))) res_rs$BA[i] = as.numeric(basalArea(treeDataStepID))
if((!draw) || (standVariable %in% c("N", "QMD"))) res_rs$N[i] = sum(treeDataStepID$N)
if((!draw) || (standVariable=="DSD")) res_rs$DSD[i] = as.numeric(diameterSD(treeDataStepID))
if((!draw) || (standVariable=="QMD")) res_rs$QMD[i] = 100*sqrt(res_rs$BA[i]/res_rs$N[i])
if((!draw) || (standVariable=="MAXD")) res_rs$MAXD[i] = max(treeDataStepID$DBH, na.rm=T)
if((!draw) || (standVariable=="HDOM")) res_rs$HDOM[i] = as.numeric(dominantHeight(treeDataStepID))
if((!draw) || (standVariable=="DDOM")) res_rs$DDOM[i] = as.numeric(dominantDiameter(treeDataStepID))
if((!draw) || (standVariable=="VCC")) {
vol = IFNvolume(treeDataStepID)
res_rs$VCC[i] = sum(vol$VCC, na.rm=T)
}
if((!draw) || (standVariable=="Biomass")) {
bio = IFNbiomass(treeDataStepID)
res_rs$Biomass[i] = (1e-3)*(sum(bio$Stem, na.rm=T) + sum(bio$Branches, na.rm=T) + sum(bio$Roots, na.rm=T))
}
}
}
if(!is.null(res)) {
res = rbind(res, res_rs)
} else {
res = res_rs
}
}
}
if(draw) {
if(standVariable=="BA") {
p<-ggplot(res)+
ylab("Basal area (m2/ha)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,BA, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="N") {
p<-ggplot(res)+
ylab("Density (tree/ha)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,N, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="DSD") {
p<-ggplot(res)+
ylab("Diameter SD (cm)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,DSD, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="QMD") {
p<-ggplot(res)+
ylab("Quad. mean diameter (cm)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,QMD, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="MAXD") {
p<-ggplot(res)+
ylab("Max. diameter (cm)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,MAXD, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="HDOM") {
p<-ggplot(res)+
ylab("Dom. height (m)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,HDOM, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="DDOM") {
p<-ggplot(res)+
ylab("Dom. diameter (cm)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,DDOM, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="VCC") {
p<-ggplot(res)+
ylab("Volume (m3/ha)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,VCC, colour=Plot), data = res[res$Run==r,], size=1)
}
} else if (standVariable=="Biomass") {
p<-ggplot(res)+
ylab("Wood biomass (Mg/ha)")
for(r in runs) {
p<-p+geom_line(mapping=aes(Year,Biomass, colour=Plot), data = res[res$Run==r,], size=1)
}
}
return(p)
} else {
return(res)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.