R/get_input.R

In TTR.PGM: Thornley Transport Resistance Plant Growth Model

#' The default c4 photosynthesis parameters
#' @description The default parameters used by \code{get_input()} to calculate
#' C4 photosynthesis rates.
#' @export
p4 <- list(
  Vcmax25 = 33.73342,
  Vpmax25 = 60.96067,
  Jmax25 = 293.6919,
  ps2a25 = 0.03532988,
  ps2b25 = 0.9282173,
  ps2c25 = 448.2845,
  gbs25 = 0.003658694,
  Rd25 = 0.2419786,
  gstar25 = 0.0003816742,
  Sco25 = 1310.019,
  Kc25 = 145.8067,
  Ko25 = 29180.07,
  Kp25 = 9.05336,
  x = 0.4081943,
  alpha = 0.0001964351,
  a025 = 0.04700031,
  Vpr = 79.99999,
  Jscale = 1937.371,
  EVpmax = 42.55894,
  EJmax = 57.48121,
  ToptJmax = 35.77227,
  HJmax = 200,
  EVcmax = 41.86028,
  ToptVcmax = 34.49091,
  HVcmax = 200,
  Eps2a = 0,        # no temperature dependency
  Eps2b = 0,        # no temperature dependency
  Eps2c = 54.76275,
  Egbs = 116.0691,
  Toptgbs = 29.20352, 
  Hgbs = 200,         
  ERd = 67.38313,
  EKc = 64.73958,
  EKo = 10.30235,
  EKp = 32.80124,
  ESco = -30.871,
  Ea0 = 1.629807,
  z = 1,
  m = 4,             # Ball Berry parameters from aDGVM2 model description
  b = 0.04,          # Ball Berry parameters from aDGVM2 model description
  g0 = 0.01,         # Medlyn parameters from CABLE doi:10.5194/gmd-8-431-2015
  g1 = 1.62,         # Medlyn parameters from CABLE doi:10.5194/gmd-8-431-2015
  bb = 0             # If using Ball Berry, make this 1; else, use Medlyn
)


#' The default c3 photosynthesis parameters
#' @description The default parameters used by \code{get_input()} to calculate
#' C3 photosynthesis rates.
#' @export
p3 <- list(
  Vcmax25 = 48.32805,
  Jmax25 = 177.2075,
  ps2a25 = 0.0007222616,
  ps2b25 = 0.9991975,
  ps2c25 = 535.0539,
  Jscale = 502.591,
  Gstar25 = 3.734099,
  gm25 = 2.770173,
  Rd25 = 1.12025,
  Kc25 = 19.51134,
  Ko25 = 20833.4,
  Eps2a = 0,        # no temperature dependency
  Eps2b = 0,        # no temperature dependency
  Eps2c = 28.6007,
  EJmax = 15.40946,
  ToptJmax = 33.12833,
  HJmax = 200,
  EVcmax = 46.40474,
  ToptVcmax = 32.67514,
  HVcmax = 200,
  EGstar = 7.64469,
  Egm = 38.67419,
  Toptgm = 30.04408, 
  Hgm = 200,         
  ERd = 41.59605,
  EKc = 40.85566,
  EKo = 40.35529,
  m = 9,             # Ball Berry parameters from aDGVM2 model description
  b = 0.01,          # Ball Berry parameters from aDGVM2 model description
  g0 = 0.01,         # Medlyn parameters from CABLE doi:10.5194/gmd-8-431-2015
  g1 = 4,            # Medlyn parameters from CABLE doi:10.5194/gmd-8-431-2015
  bb = 0             # If using Ball Berry, make this 1; else, use Medlyn
)


#' Create an input data object for \code{make_data()}
#' @description This function takes the input forcing data and formats it for use by the \code{make_data()} function.
#' The input forcing data are expected to have the units specified below.
#' All of the forcing data should be in dataframe or matrix form.
#' The following constant dimension sizes are expected:
#'
#' nspecies : The number of species
#'
#' nsites : The number of sites
#'
#' nsteps : The number of timesteps
#'
#' Currently only one species is supported.
#' Depending on the model configuration different parameters need to be specified. tcur, tnur, tgrowth and tloss always need to be supplied. For the std variant of the process model, rsds has to be specified as well. For other variants  photosynthesis rates can be supplied or calculated from tcur, catm and rsds. Soil water content needs to be supplied or it will be calculated from wp, fc, prec and pet. nsoil is optional. fire is currently not supported.
#' 
#' @param observations The observation data that will be used in the statistical model, there are no checks performed, so make sure to format this in the way the user defined statistical model expects
#' @param seconds The number of seconds in a time step
#' @param p3 The parameters to be used to calculate C3 photosynthesis rates (see \code{p3})
#' @param p4 The parameters to be used to calculate C4 photosynthesis rates (see \code{p4})
#' @param lon vector of length nsite containing the Longitudes of the sites
#' @param lat vector of length nsite containing the Latitudes of the sites
#' @param tcur matrix with ntimesteps rows, nsites columns, containing the temperature associated with photosynthesis [°C]
#' @param tnur matrix with ntimesteps rows, nsites columns, containing the temperature associated with nitrogen uptake [°C]
#' @param tgrowth matrix with ntimesteps rows, nsites columns, containing the temperature associated with growth [°C]
#' @param tloss matrix with ntimesteps rows, nsites columns, containing the temperature associated with biomass loss [°C]
#' @param photoc3 matrix with ntimesteps rows, nsites columns, containing c3 photosynthesis rates [mumol*m^(-2)*s^(-1)]
#' @param photoc4 matrix with ntimesteps rows, nsites columns, containing c4 photosynthesis rates [mumol*m^(-2)*s^(-1)]
#' @param rsds matrix with ntimesteps rows, nsites columns, containing short wave downward solar radiation at the surface [W*m^(-2)]
#' @param catm matrix with ntimesteps rows, nsites columns, containing Partial pressure of CO2 in the atmosphere [Pa]
#' @param swc matrix with ntimesteps rows, nsites columns, containing soil water content scaled so that wilting point = 0 and field capacity = 100 [\%]
#' @param pet matrix with ntimesteps rows, nsites columns, containing potential evapotranspiration [kg*m^(-2)*s^(-1)]
#' @param rain matrix with ntimesteps rows, nsites columns, containing precipitation [kg*m^(-2)*s^(-1)]
#' @param wp vector of length nsites containing wilting points [cm^3/cm^3]
#' @param fc vector of length nsites containing field capacities [cm^3/cm^3]
#' @param fire matrix with ntimesteps rows, nsites columns, containing binary information on fire occurrence [no units]; currently not supported
#' @param nsoil vector of length nsites containing the soil nitrogen content [kg/kg]
#' @return A list containing the forcing data for the process model:
#' \item{timeseries}{
#' All forcing variables which change with time
#' }
#' \item{timeinvariant}{
#' All forcing variables which do not change with time
#' }
#' \item{observations}{
#' Observations as specified above
#' }
#' @export
get_input <- function(observations,
                      tcur,
                      tnur,
                      tgrowth,
                      tloss,
                      seconds,
                      p3 = TTR.PGM::p3,
                      p4 = TTR.PGM::p4,
                      lon = NULL,
                      lat = NULL,
                      rsds = NULL,
                      catm = NULL,
                      photoc3 = NULL,
                      photoc4 = NULL,
                      swc = NULL,
                      pet = NULL,
                      rain = NULL,
                      wp = NULL,
                      fc = NULL,
                      fire = NULL,
                      nsoil = NULL){
  #check if input data requirements for swc are given
  if(is.null(swc) & (is.null(wp) | is.null(fc) | is.null(rain) | is.null(pet))){
    stop("If no soil water content matrix is given, please provide wp, fc, pet and rain")
  }
  if((is.null(photoc3) | is.null(photoc4)) & (is.null(rsds) | is.null(catm))){
    stop("If photosynthesis parameters are not given, please provide rsds and catm")
  }
  to_column<- function(vec){
    if(any("data.frame" == class(vec))){
      if(dim(vec)[[2]] != 1){
        stop("More than one column given for column variable. Please specify a single column matrix, dataframe or a numeric vector")
      }
      col <- vec[[1]]
      return(as.array(col, dim=c(length(col),1)))
    }
    if(is.null(dim(vec))){
      return(as.array(vec, dim=c(length(vec),1)))
    }
    return(vec)
  }
  matrix_cast <- function(df){
    unname(as.matrix(df))
  }

  #cast observations to matrix
  observations <- as.matrix(observations)

  #we don't yet know the dimensionality
  nsteps <- NULL
  nsites <- NULL
  nspecies <- NULL

  #this is the minimum data required for all calls to ttr
  req <- list(tcur = tcur,
             tnur = tnur,
             tgrowth = tgrowth,
             tloss = tloss
  )
  req <- mapply(
    function(d,n){assertMatrix(matrix_cast(d),any.missing=FALSE,.var.name=n)},
    d = req,
    n = names(req),
    SIMPLIFY=FALSE)

  req_dim_one <- unique(vapply(req, function(d){dim(d)[1]}, FUN.VALUE=1))
  req_dim_two <- unique(vapply(req, function(d){dim(d)[2]}, FUN.VALUE=1))
  if(!testNumeric(req_dim_one, min.len=1, max.len=1, any.missing=FALSE) |
     !testNumeric(req_dim_two, min.len=1, max.len=1, any.missing=FALSE)){
    print(req_dim_one)
    print(req_dim_two)
    stop("Dimensions differ among input data.")
  }

  if(!is.null(nsteps)){
    if(nsteps != req_dim_one){
      stop("Input data dimension 1 (nsteps) does not add up with observational data.")
    }
  } else{
    nsteps <- req_dim_one
  }
  if(!is.null(nsites)){
    if(nsites != req_dim_two){
      stop("Input data dimension 2 (nsites) does not add up with observational data.")
    }
  } else{
    nsites <- req_dim_two
  }
  lon <- as.vector(lon)
  lat <- as.vector(lat)
  if(!is.null(lon) & !is.null(lat)){
    if(length(lon) != length(lat) || length(lon) != nsites || length(lat) != nsites){
      stop("Longitude or latitude vector's length does not fit the data.")
    }
  }
  lonlat <- matrix(data=c(lon,lat),ncol=2)
  colnames(lonlat) <- c("lon", "lat")

  #from this point, nsites and nsteps is known.
  per_site_default <- list(
    wp = 0,
    fc = 1,
    nsoil = 1
  )
  assert_per_site <- function(v, std = 0, n = ""){
    if(is.null(v)) rep(std,nsites) else assertNumeric(assertArray(to_column(v),d=1, .var.name=n),len=nsites,any.missing=FALSE, .var.name=n)
  }
  #TODO: is SIMPLIFY reliable here?
  cols <- mapply(assert_per_site,
                 v = list(wp = wp, fc = fc, nsoil = nsoil),
                 std = per_site_default,
                 n = names(per_site_default),
                 SIMPLIFY=FALSE)
  time.invariant <- aperm(abind(cols, along=2), c(2,1))

  assert_timeseries <- function(v, std = 0, n = ""){
    if(is.null(v)) matrix(data = std, nrow = nsteps, ncol = nsites) else assertMatrix(matrix_cast(v), nrows=nsteps, ncols=nsites, any.missing=FALSE, .var.name=n)
  }
  timeseries_default <- list(
    rain = 0,
    pet = 0,
    fire = 0,
    catm = 33.8,
    rsds = 0
  )
  ts <- mapply(assert_timeseries,
               v = list(rain = rain,
                        pet = pet,
                        fire = fire,
                        catm = catm,
                        rsds = rsds),
               std=timeseries_default,
               n = names(timeseries_default),
               SIMPLIFY=FALSE)
  swc <- if(is.null(swc)) make_swc(ts$rain, ts$pet, cols$fc, cols$wp, seconds, 3) else assertMatrix(matrix_cast(swc), nrows=nsteps, ncols=nsites, any.missing=FALSE,.var.name="swc")

  J_to_mol=4.6;
  frac_PAR=0.5;
  ppfd <- ts$rsds * frac_PAR * J_to_mol #bigleaf conversion from rsds W/m2 to PAR umol/m2/s
  tleaf <- req$tcur + 273.15
  photoc3 <- if(is.null(photoc3)) t(calc_photo(t(tleaf), t(ppfd), t(ts$catm), p3, "C3")) else assertMatrix(matrix_cast(photoc3), nrows=nsteps, ncols=nsites, any.missing=FALSE,.var.name="photoc3")
  photoc4 <- if(is.null(photoc4)) t(calc_photo(t(tleaf), t(ppfd), t(ts$catm), p4, "C4")) else assertMatrix(matrix_cast(photoc4), nrows=nsteps, ncols=nsites, any.missing=F,.var.name="photoc4")

  timeseries_matrices <- list(tcur = req$tcur,
                              tnur = req$tnur,
                              tgrowth = req$tgrowth,
                              tloss = req$tloss,
                              ppfd = ppfd,
                              swc = swc,
                              prec = ts$rain,
                              pet = ts$pet,
                              C3p = photoc3,
                              C4p = photoc4,
                              fire = ts$fire)

  timeseries <- aperm(abind(timeseries_matrices, along=3),c(3,1,2))

  list(dimensions = c(nspecies = nspecies, nsteps = nsteps, nsites = nsites), time.invariant = time.invariant, timeseries = timeseries, observations=observations, lon = lon, lat = lat, lonlat = lonlat)
}