R/drawStandDynamics.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

#' 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)
  }
}