R/hydraulics_vulnerabilityCurvePlot.R

In medfate: Mediterranean Forest Simulation

#' @rdname hydraulics_conductancefunctions
#' 
#' @param x An object of class \code{\link{spwbInput}}.
#' @param soil A list containing the description of the soil (see \code{\link{soil}}).
#' @param type Plot type for \code{hydraulics_vulnerabilityCurvePlot}, either \code{"leaf"}, 
#'             \code{"stem"}, \code{"root"} or \code{"rhizosphere"}).
#' @param psiVec Vector of water potential values to evaluate for the vulnerability curve.
#' @param relative A flag to relativize vulnerability curves to the [0-1] interval.
#' @param speciesNames A flag to indicate the use of species names instead of cohort names in plots.
#' @param draw A flag to indicate whether the vulnerability curve should be drawn or just returned.
#' @param ylim,xlab,ylab Graphical parameters to override function defaults.
#' 
hydraulics_vulnerabilityCurvePlot<-function(x, soil = NULL, type="leaf",  psiVec =  seq(-0.1, -8.0, by =-0.01), 
                                            relative = FALSE, speciesNames = FALSE, 
                                            draw = TRUE,  ylim = NULL, xlab = NULL, ylab = NULL) {
  
  TYPES = c("leaf","stem","root","rhizosphere")
  type = match.arg(type,TYPES)  
  
  cohortnames = row.names(x$cohorts)
  spnames = x$cohorts$Name
  VCroot_kmax_layer = x$belowLayers$VCroot_kmax
  VGrhizo_kmax = x$belowLayers$VGrhizo_kmax
  nlayer = ncol(VCroot_kmax_layer)
  col = rainbow(nlayer, start = 0.8, end = 0.1)
  
  VCroot_kmax = x$paramsTransp$VCroot_kmax
  VCroot_c = x$paramsTransp$VCroot_c
  VCroot_d = x$paramsTransp$VCroot_d
  VCstem_kmax = x$paramsTransp$VCstem_kmax
  VCstem_c = x$paramsTransp$VCstem_c
  VCstem_d = x$paramsTransp$VCstem_d
  VCleaf_kmax = x$paramsTransp$VCleaf_kmax
  VCleaf_c = x$paramsTransp$VCleaf_c
  VCleaf_d = x$paramsTransp$VCleaf_d
  ncoh = nrow(x$above)

  if(relative) {
    VCstem_kmax = rep(1, ncoh)
    VCleaf_kmax = rep(1, ncoh)
    VCroot_kmax[] = 1
    VGrhizo_kmax[] = 1
  }
  res = NULL
  minPsi = min(psiVec)
  if(type=="leaf") {
    if(is.null(ylab)) ylab = expression(paste("Leaf conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
    if(is.null(xlab)) xlab = "Leaf pressure (-MPa)"
    res = matrix(0, ncol=ncoh, nrow = length(psiVec))
    rownames(res) = psiVec
    colnames(res) = cohortnames
    if(speciesNames) colnames(res) = spnames
    for(i in 1:ncoh) {
      res[,i] = unlist(lapply(psiVec, hydraulics_xylemConductance, VCleaf_kmax[i], VCleaf_c[i], VCleaf_d[i]))
    }
    if(draw) {
      return(.multiple_y(-psiVec, res, xlab=xlab, ylab=ylab, ylim=ylim))
    }
  } 
  else if(type=="stem") {
    res = matrix(0, ncol=ncoh, nrow = length(psiVec))
    rownames(res) = psiVec
    colnames(res) = cohortnames
    if(speciesNames) colnames(res) = spnames
    if(is.null(ylab)) ylab = expression(paste("Stem conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
    if(is.null(xlab)) xlab = "Stem pressure (-MPa)"

    for(i in 1:ncoh) {
      res[,i] = unlist(lapply(psiVec, hydraulics_xylemConductance, VCstem_kmax[i], VCstem_c[i], VCstem_d[i]))
    }
    if(draw) {
      return(.multiple_y(-psiVec, res, xlab=xlab, ylab=ylab, ylim=ylim))
    }
  } 
  else if(type=="root") {
    res = matrix(0, ncol=ncoh, nrow = length(psiVec))
    rownames(res) = psiVec
    colnames(res) = cohortnames
    if(speciesNames) colnames(res) = spnames
    if(is.null(ylab)) ylab = expression(paste("Root conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
    if(is.null(xlab)) xlab = "Root pressure (-MPa)"
    for(i in 1:ncoh) {
      res[,i] = unlist(lapply(psiVec, hydraulics_xylemConductance, VCroot_kmax[i], VCroot_c[i], VCstem_d[i]))
    }
    if(draw) {
      return(.multiple_y(-psiVec, res, xlab=xlab, ylab=ylab, ylim=ylim))
    }
  } 
  else if(type=="rootlayer") {
    if(is.null(ylab)) ylab = expression(paste("Root layer conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
    if(is.null(xlab)) xlab = "Root layer pressure (-MPa)"
    if(is.null(ylim)) ylim =c(0, max(VCroot_kmax_layer))
    
    for(i in 1:ncoh) {
      kroot = matrix(NA, ncol=length(psiVec), nrow = nlayer)
      for(j in 1:nlayer) {
        kroot[j,] = unlist(lapply(psiVec, hydraulics_xylemConductance, VCroot_kmax_layer[i,j], VCroot_c[i], VCroot_d[i]))
      }
      if(draw) {
        if(i==1) {
          matplot(-psiVec, t(kroot), type="l", xlim=c(0,-minPsi), ylim =ylim,
                  xlab = xlab, ylab = ylab, col=col)
        } else {
          matlines(-psiVec, t(kroot), lty=i, col=col)
        }
      }
    }
    if(draw) {
      if(!speciesNames) legend("topright", legend = cohortnames, lty=1:ncoh, col = 1:ncoh, bty="n")
      else legend("topright", legend = spnames, lty=1:ncoh, col = 1:ncoh, bty="n")
      legend("right", legend = paste("Layer", 1:nlayer), lty=1, col=col, bty="n")
    }
  } 
  else if(type=="rhizosphere") {
    if(is.null(soil)) stop("Please supply a non-null soil object!")
    if(is.null(ylab)) ylab = expression(paste("Rhizosphere conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
    if(is.null(xlab)) xlab = "Rhizosphere pressure (-MPa)"
    if(is.null(ylim)) ylim =c(0,100)
    VG_nc = soil$VG_n
    VG_alphac = soil$VG_alpha
    for(i in 1:ncoh) {
      krhizo = matrix(NA, ncol=length(psiVec), nrow = nlayer)
      for(j in 1:nlayer) {
        krhizo[j,] = unlist(lapply(psiVec, hydraulics_vanGenuchtenConductance, VGrhizo_kmax[i,j], VG_nc[j], VG_alphac[j]))
      }
      if(i==1) {
        matplot(-psiVec, t(krhizo), type="l", xlim=c(0,-minPsi), ylim =ylim,
                xlab = xlab, ylab = ylab, col=col)
      } else {
        matlines(-psiVec, t(krhizo), lty=i, col=col)
      }
    }
    if(!speciesNames) legend("topright", legend = cohortnames, lty=1:ncoh, col = 1:ncoh, bty="n")
    else legend("topright", legend = spnames, lty=1:ncoh, col = 1:ncoh, bty="n")
    legend("right", legend = paste("Layer", 1:nlayer), lty=1, col=col, bty="n")
  } 
  if(draw) invisible(res)
  else return(res)
}