R/VerticalProfiles.R

In medfate: Mediterranean Forest Simulation

#' Vertical profiles
#' 
#' Functions to generate vertical profiles generated by an input \code{\link{forest}} object.
#' 
#' @param x An object of class \code{\link{forest}}
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#' @param z A numeric vector with height values.
#' @param d A numeric vector with soil layer widths.
#' @param gdd Growth degree days.
#' @param byCohorts Separate profiles for each cohort.
#' @param bySpecies Aggregate cohort profiles by species.
#' @param includeHerbs Include herbaceous layer in the profile.
#' @param u The value of measured wind speed (in m/s).
#' @param windMeasurementHeight Height corresponding to wind measurement (in cm over the canopy).
#' @param boundaryLayerSize Size of the boundary layer (in cm) over the canopy.
#' @param target Wind property to draw, either "windspeed", "kineticenergy" (turbulent kinetic energy) or "stress" (Reynold's stress).
#' @param draw Logical flag to indicate that a plot is desired.
#' @param xlim Limits of the x-axis.
#' 
#' @return If \code{draw = FALSE}, the functions return a numeric vector with values measured at each height. Units depend on the profile function:
#' \itemize{
#'   \item{\code{vprofile_leafAreaDensity}: Cumulative LAI (m2/m2) per height bin.}
#'   \item{\code{vprofile_fuelBulkDensity}: Fuel bulk density (kg/m3) per height bin.}
#'   \item{\code{vprofile_PARExtinction}: Percent of photosynthetically active radiation (\%) corresponding to each height.}
#'   \item{\code{vprofile_SWRExtinction}: Percent of shortwave radiation (\%) corresponding to each height.}
#'   \item{\code{vprofile_windExtinction}: Wind speed (m/s) corresponding to each height.}
#' }
#' 
#' If \code{draw = TRUE} the functions return a ggplot object, instead.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{forest}}, \code{\link{plot.forest}}, \code{\link{wind_canopyTurbulence}}
#' 
#' @examples 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Load example plot plant data
#' data(exampleforest)
#' 
#' vprofile_leafAreaDensity(exampleforest, SpParamsMED)
#' vprofile_fuelBulkDensity(exampleforest, SpParamsMED)
#' 
#' vprofile_PARExtinction(exampleforest, SpParamsMED)
#' vprofile_SWRExtinction(exampleforest, SpParamsMED)
#' 
#' vprofile_windExtinction(exampleforest, SpParamsMED)
#' 
#' @keywords internal
#' @name vprofile_leafAreaDensity
vprofile_leafAreaDensity<-function(x, SpParams = NULL, z = NULL, gdd = NA, 
                                   byCohorts = FALSE, bySpecies = FALSE, includeHerbs = FALSE,
                                   draw = TRUE, xlim = NULL) {
  if(!(inherits(x,"data.frame") || inherits(x, "forest"))) stop("'x' should be of class 'forest' or 'data.frame'")
  if(inherits(x, "forest")) {
    if(is.null(SpParams)) stop("Please, provide 'SpParams' to calculate leaf area.")
    spnames <- plant_speciesName(x, SpParams)
    herbHeight <- x$herbHeight
    herbLAI <- herb_LAI(x, SpParams)
    x <- forest2aboveground(x, SpParams, gdd)
  } else {
    if(any(!(c("LAI_expanded", "H", "CR", "SP") %in% names(x)))) {
      stop("Data frame should contain columns 'SP', 'LAI_expanded', 'H' and 'CR'")
    }
    spnames <- x$SP
  }
  cohortnames <- row.names(x)
  
  if(is.null(z)) z <- seq(0, ceiling(max(x$H)/100)*100 +10, by=1)
  w <- z[2:length(z)]- z[1:(length(z)-1)]
  
  lai_vect <- x$LAI_expanded
  h_vect <- x$H
  cr_vect <- x$CR
  if(includeHerbs) {
    lai_vect <- c(lai_vect, herbLAI)
    h_vect <- c(h_vect, herbHeight)
    cr_vect <- c(x$CR, 1.0)
    spnames <- c(spnames, "Herbaceous")
    cohortnames <- c(cohortnames, "Herbaceous")
  }
  if(!byCohorts) {
    lai <- .LAIprofileVectors(z, lai_vect, h_vect, cr_vect)
    lai <- 100*lai/w
    if(draw) {
      df <- data.frame("lai" = c(0,lai), "z" = z)
      g<-ggplot(df, aes(x=.data$lai, y=.data$z))+
        geom_path()+
        xlab("Leaf Area Density (m2/m3)")+
        ylab("Height (cm)")+
        theme_bw()
      if(!is.null(xlim)) g <- g + xlim(xlim)
    }
  } else {
    lai <- .LAIdistributionVectors(z, lai_vect, h_vect, cr_vect)
    lai <- 100*sweep(lai,1,w, "/")
    if(bySpecies) {
      lai <- t(apply(lai,1, tapply, spnames, sum, na.rm=T))
      cohortnames <- colnames(lai)
    } 
    if(draw) {
      lai <- rbind(rep(0, ncol(lai)),lai)
      df <- data.frame("lai" = as.vector(lai), "z" = z,
                      "Cohort" = gl(length(cohortnames), nrow(lai), labels=cohortnames))
      g<-ggplot(df, aes(x=.data$lai, y=.data$z))+
        geom_path(aes(col=.data$Cohort, linetype = .data$Cohort))+
        xlab("Leaf Area Density (m2/m3)")+
        ylab("Height (cm)")+
        scale_color_discrete(name="")+
        scale_linetype_discrete(name="")+
        theme_bw()
      if(!is.null(xlim)) g <- g + xlim(xlim)
    }
  }
  if(draw) return(g)
  else return(lai)
}

#' @rdname vprofile_leafAreaDensity
#' @keywords internal
vprofile_rootDistribution<-function(x, SpParams, d = NULL, bySpecies = FALSE, 
                                    draw = TRUE, xlim = NULL) {
  if(is.null(d)){
    zmax <- 0
    if(nrow(x$shrubData)>0) zmax <- max(zmax, ceiling(max(x$shrubData$Z95)/100)*100)
    if(nrow(x$treeData)>0) zmax <- max(zmax, ceiling(max(x$treeData$Z95)/100)*100)
    d <- rep(10,1+(zmax/10))
    z <- seq(0,zmax, by=10)
  } else {
    zmax <- sum(d)
    z <- numeric(length(d))
    z[1] <- 0
    for(i in 2:length(z)) z[i] <- z[i-1] + d[i-1]
  }
  cohortnames <- plant_ID(x, SpParams)
  rd <- .rootDistribution(d,x)
  rownames(rd) <- cohortnames
  if(bySpecies) {
    spnames <- plant_speciesName(x, SpParams)
    rd <- apply(rd,2, tapply, spnames, mean, na.rm=T)
  } 
  if(draw) {
    return(.multiple_x(x=t(rd)*100, y=z/10, xlab="% of fine roots/cm",
                       ylab="Depth (cm)", ylim=c(zmax/10,0), xlim = xlim))
  } else return(rd)
}

#' @rdname vprofile_leafAreaDensity
#' @keywords internal
vprofile_fuelBulkDensity<-function(x, SpParams, z = NULL, gdd = NA,
                                   draw = TRUE, xlim = NULL) {
  if(is.null(z)) z <- seq(0, ceiling(max(plant_height(x, SpParams))/100)*100 +10, by=1)
  wfp <- .woodyFuelProfile(z,x, SpParams, gdd)
  df <- data.frame("BD" = c(0,wfp), "Z" = z)
  if(draw) {
    g<-ggplot(df, aes(x=.data$BD, y=.data$Z))+
      geom_path()+
      xlab("Bulk density (kg/m3)")+
      ylab("Height (cm)")+
      theme_bw()
    if(!is.null(xlim)) g <- g + xlim(xlim)
  }
  if(draw) return(g)
  else return(wfp)
}

#' @rdname vprofile_leafAreaDensity
#' @keywords internal
vprofile_PARExtinction<-function(x, SpParams, z = NULL, gdd = NA, includeHerbs = FALSE, 
                                 draw = TRUE, xlim = c(0,100)) {
  if(is.null(z)) z <- seq(0, ceiling(max(plant_height(x, SpParams), na.rm = TRUE)/100)*100 +10, by=1)
  pep <- .parExtinctionProfile(z,x, SpParams, gdd, includeHerbs)
  df <- data.frame("PEP" = pep, "Z" = z)
  if(draw) {
    g<-ggplot(df, aes(x=.data$PEP, y=.data$Z))+
      geom_path()+
      xlab("Available PAR (%)")+
      ylab("Height (cm)")+
      theme_bw()
    if(!is.null(xlim)) g <- g + xlim(xlim)
  }
  if(draw) return(g)
  else return(pep)
}

#' @rdname vprofile_leafAreaDensity
#' @keywords internal
vprofile_SWRExtinction<-function(x, SpParams, z = NULL, gdd = NA, includeHerbs = FALSE, 
                                 draw = TRUE, xlim = c(0,100)) {
  if(is.null(z)) z <- seq(0, ceiling(max(plant_height(x, SpParams))/100)*100 +10, by=1)
  swr <- .swrExtinctionProfile(z,x, SpParams, gdd, includeHerbs)
  df <- data.frame("SWR" = swr, "Z" = z)
  if(draw) {
    g<-ggplot(df, aes(x=.data$SWR, y=.data$Z))+
      geom_path()+
      xlab("Available SWR (%)")+
      ylab("Height (cm)")+
      theme_bw()
    if(!is.null(xlim)) g <- g + xlim(xlim)
  }
  if(draw) return(g)
  else return(swr)
}

#' @rdname vprofile_leafAreaDensity
#' @keywords internal
vprofile_windExtinction<-function(x, SpParams, u = 1, windMeasurementHeight = 200,
                                  boundaryLayerSize = 2000, target = "windspeed",
                                  z = NULL, gdd = NA, includeHerbs = FALSE,
                                  draw = TRUE, xlim = NULL) {
  if(is.null(z)) z <- seq(0, ceiling(max(plant_height(x, SpParams))/100)*100 +boundaryLayerSize, by=1)
  lad <- vprofile_leafAreaDensity(x=x, SpParams = SpParams, z = z, gdd = gdd,
                                  byCohorts = FALSE, bySpecies = FALSE, includeHerbs = includeHerbs,
                                  draw = FALSE, xlim = xlim)
  canopyHeight <- max(plant_height(x, SpParams), na.rm=T)
  zmid <- 0.5*(z[1:(length(z)-1)] + z[2:(length(z))])
  df <- wind_canopyTurbulence(zmid, lad, canopyHeight, u, windMeasurementHeight)
  target <- match.arg(target, c("windspeed", "kineticenergy", "stress"))
  if(draw) {
    if(target=="windspeed") {
      g<-ggplot(df, aes(x=.data$u, y=.data$zmid))+
        geom_path()+
        xlab("Wind speed (m/s)")
    } else if(target == "kineticenergy") {
      g<-ggplot(df, aes(x=.data$k, y=.data$zmid))+
        geom_path()+
        xlab("Turbulent kinetic energy (m2/s2)")
    } else if(target == "stress") {
      g<-ggplot(df, aes(x=.data$uw, y=.data$zmid))+
        geom_path()+
        xlab("Reynolds stress (m2/s2)")
    }
    g <- g + geom_hline(yintercept=canopyHeight, col="gray", linetype="dashed")+
      ylab("Height (cm)")+
      theme_bw()
    if(!is.null(xlim)) g <- g + xlim(xlim)
  }
  if(draw) return(g)
  else return(df)
}