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 vulnerabilityFunction String indicating the function used to represent vulnerability in plant segments, either "Weibull" or "Sigmoid".
#' @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.
#' 
#' @keywords internal
hydraulics_vulnerabilityCurvePlot<-function(x, soil = NULL, type="leaf",  vulnerabilityFunction = "Weibull",
                                            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","rootlayer","rhizosphere")
  type <- match.arg(type,TYPES)  
  
  TYPES_FUN <- c("Weibull", "Sigmoid")
  vulnerabilityFunction <- match.arg(vulnerabilityFunction,TYPES_FUN)  
  
  if(!(x$control$transpirationMode %in% c("Sperry", "Sureau"))) stop("x has to be initialized using 'Sperry' or 'Sureau' transpiration modes")

  
  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)
  
  ncoh <- nrow(x$above)

  if(relative) {
    VCroot_kmax_layer[] <- 100
    VGrhizo_kmax[] <- 100
  }
  res <- NULL
  minPsi <- min(psiVec)
  if(type=="leaf") {
    if(is.null(ylab)) {
      if(!relative) ylab <- expression(paste("Leaf conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
      else ylab <- "Leaf conductance (% of maximum)"
    }
    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
    if(!relative) VCleaf_kmax <- x$paramsTransp$VCleaf_kmax
    else VCleaf_kmax <- rep(100, ncoh)
    if(vulnerabilityFunction=="Weibull"){
      VCleaf_c <- x$paramsTransp$VCleaf_c
      VCleaf_d <- x$paramsTransp$VCleaf_d
      for(i in 1:ncoh) {
        res[,i] <- unlist(lapply(psiVec, hydraulics_xylemConductance, VCleaf_kmax[i], VCleaf_c[i], VCleaf_d[i]))
      }
    } else {
      VCleaf_P50 <- x$paramsTransp$VCleaf_P50
      VCleaf_slope <- x$paramsTransp$VCleaf_slope
      for(i in 1:ncoh) res[,i] <-VCleaf_kmax[i]*(1.0 - (1.0/(1.0 + exp(VCleaf_slope[i] / 25.0 * (psiVec - VCleaf_P50[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)) {
      if(!relative) ylab <- expression(paste("Stem conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
      else  ylab <- "Stem conductance (% of maximum)"
    }
    if(is.null(xlab)) xlab <- "Stem pressure (-MPa)"
    if(!relative) VCstem_kmax <- x$paramsTransp$VCstem_kmax
    else VCstem_kmax <- rep(100, ncoh)
    if(vulnerabilityFunction=="Weibull") {
      VCstem_c <- x$paramsTransp$VCstem_c
      VCstem_d <- x$paramsTransp$VCstem_d
      for(i in 1:ncoh) {
        res[,i] <- unlist(lapply(psiVec, hydraulics_xylemConductance, VCstem_kmax[i], VCstem_c[i], VCstem_d[i]))
      }
    } else {
      VCstem_P50 <- x$paramsTransp$VCstem_P50
      VCstem_slope <- x$paramsTransp$VCstem_slope
      for(i in 1:ncoh) res[,i] <-VCstem_kmax[i]*(1.0 - (1.0/(1.0 + exp(VCstem_slope[i] / 25.0 * (psiVec - VCstem_P50[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)) {
      if(!relative) ylab <- expression(paste("Root conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
      else ylab <- "Root conductance (% of maximum)"
    }
    if(is.null(xlab)) xlab <- "Root pressure (-MPa)"
    if(!relative) VCroot_kmax <- x$paramsTransp$VCroot_kmax
    else VCroot_kmax <- rep(100, ncoh)
    if(vulnerabilityFunction == "Weibull"){
      VCroot_c <- x$paramsTransp$VCroot_c
      VCroot_d <- x$paramsTransp$VCroot_d
      for(i in 1:ncoh) {
        res[,i] <- unlist(lapply(psiVec, hydraulics_xylemConductance, VCroot_kmax[i], VCroot_c[i], VCroot_d[i]))
      }
    } else {
      VCroot_P50 <- x$paramsTransp$VCroot_P50
      VCroot_slope <- x$paramsTransp$VCroot_slope
      for(i in 1:ncoh) res[,i] <-VCroot_kmax[i]*(1.0 - (1.0/(1.0 + exp(VCroot_slope[i] / 25.0 * (psiVec - VCroot_P50[i])))))
    }
    if(draw) {
      return(.multiple_y(-psiVec, res, xlab=xlab, ylab=ylab, ylim=ylim))
    }
  } 
  else if(type=="rootlayer") {
    if(is.null(ylab)) {
      if(!relative) ylab <- expression(paste("Root layer conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
      else ylab <- "Root layer conductance (% of maximum)"
    }
    if(is.null(xlab)) xlab <- "Root layer pressure (-MPa)"
    if(is.null(ylim)) ylim <- c(0, max(VCroot_kmax_layer))
    
    if(vulnerabilityFunction=="Weibull") {
      VCroot_c <- x$paramsTransp$VCroot_c
      VCroot_d <- x$paramsTransp$VCroot_d
      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)
          }
        }
      }
    } else {
      VCroot_P50 <- x$paramsTransp$VCroot_P50
      VCroot_slope <- x$paramsTransp$VCroot_slope
      for(i in 1:ncoh) {
        kroot <- matrix(NA, ncol=length(psiVec), nrow = nlayer)
        for(j in 1:nlayer) {
          kroot[j,] <- VCroot_kmax_layer[i,j]*(1.0 - (1.0/(1.0 + exp(VCroot_slope[i] / 25.0 * (psiVec - VCroot_P50[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(!inherits(soil, "soil")) {
      if(inherits(soil, "data.frame")){
        soil <- soil(soil)
      } else {
        stop("soil has to be of class 'data.frame' or 'soil'")
      }
    } 
    if(is.null(ylab)) {
      if(!relative) ylab <- expression(paste("Rhizosphere conductance  ", (mmol%.%s^{-1}%.%m^{-2}%.%MPa^{-1})))
      else ylab <- "Rhizosphere conductance (% of maximum)"
    }
    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)
}