#' Calculate the composition (land cover) representing a
#' potential natural vegetation based on climate relationships
#'
#' The function returns relative abundance/land cover values that completely
#' cover the surface (i.e., they sum to 1) of a site specified by climate and/or
#' fixed input values.
#'
#' @section Details: Some of the land cover/vegetation types, i.e., trees,
#' annual grasses, and bare-ground are not estimated from climate
#' relationships; they are either set to 0, or alternatively fixed at the
#' value of the input argument(s).
#'
#' @section Details: The remaining vegetation types, i.e., shrubs, C3 grasses,
#' C4 grasses, forbs, and succulents, are estimated from climate relationships
#' using equations developed by Paruelo & Lauenroth 1996, or alternatively
#' fixed at the value of the input argument(s). If values for
#' \code{dailyC4vars} are provided, then equations developed by Teeri & Stowe
#' 1976 are used to limit the occurrence of C4 grasses.
#'
#' @section Details: The relative abundance values of the the vegetation types
#' that can be estimated and are not fixed by inputs, are estimated in two
#' steps: (i) as if they cover the entire surface; (ii) scaled to the
#' proportion of the surface that is not fixed by inputs.
#'
#' @section Notes: The equations developed Paruelo & Lauenroth 1996 are based
#' on sites with \var{MAT} from 2 C to 21.2 C and \var{MAP} from 117 to
#' 1011 mm. If \code{warn_extrapolation} is set to \code{TRUE}, then
#' inputs are checked against supported ranges, i.e., if \var{MAT} is below
#' 1 C, then it is reset to 1 C with a warning. If other inputs exceed their
#' ranges, then a warning is issued and the code proceeds.
#'
#' @param MAP_mm A numeric value. Mean annual precipitation in millimeter.
#' @param MAT_C A numeric value. Mean annual temperature in degree Celsius.
#' @param mean_monthly_ppt_mm A numeric vector of length 12. Mean monthly
#' precipitation in millimeter.
#' @param mean_monthly_Temp_C A numeric vector of length 12. Mean monthly
#' temperature in degree Celsius.
#' @param dailyC4vars A named list of length 3 or \code{NULL}.
#' If not \code{NULL}, then a correction for C4-grasses based on
#' Teeri & Stowe 1976 is applied based on the variables: \describe{
#' \item{\code{Month7th_NSadj_MinTemp_C}}{Mean minimum temperature of July on
#' the northern hemisphere and January on the southern hemisphere}
#' \item{\code{DegreeDaysAbove65F_NSadj_DaysC}}{Degree days above 65 F = 18.33
#' C in units of days x degree Celsius}
#' \item{\code{LengthFreezeFreeGrowingPeriod_NSadj_Days}}{Mean annual number
#' of days of the longest continuous period where minimum daily temperature
#' remain above freezing} }
#' @param isNorth A logical value. \code{TRUE} for locations on northern
#' hemisphere.
#' @param shrub_limit A numeric value. Default value is 0.2 based on page 1213
#' of Paruelo & Lauenroth 1996.
#' @param fix_succulents A logical value. If \code{TRUE}, then value for the
#' succulent component is fixed at \code{Succulents_Fraction} instead of
#' calculated from climatic relationships.
#' @param Succulents_Fraction A numeric value between 0 and 1. \code{NA} is
#' treated as if \code{fix_succulents} is \code{FALSE}.
#' @param fix_sumgrasses A logical value. If \code{TRUE}, then the summed
#' values of the three grass types, i.e., annuals, C3, and C4 grasses,
#' is fixed at \code{SumGrasses_Fraction}. If they are not fixed themselves,
#' then their values are estimated and scaled to sum to
#' \code{SumGrasses_Fraction}.
#' @param SumGrasses_Fraction A numeric value between 0 and 1.
#' \code{NA} is treated as if \code{fix_sumgrasses} is \code{FALSE}.
#' @param fix_annuals A logical value. If \code{TRUE}, then value for the annual
#' component is fixed at \code{Annuals_Fraction}.
#' @param Annuals_Fraction A numeric value between 0 and 1. Default value is 0.
#' @param fix_C4grasses A logical value. If \code{TRUE}, then value for the
#' C4-grass component is fixed at \code{C4_Fraction} instead of calculated
#' from climatic relationships.
#' @param C4_Fraction A numeric value between 0 and 1. \code{NA} is treated as
#' if \code{fix_C4grasses} is \code{FALSE}.
#' @param fix_C3grasses A logical value. If \code{TRUE}, then value for the
#' C3-grass component is fixed at \code{C3_Fraction} instead of calculated
#' from climatic relationships.
#' @param C3_Fraction A numeric value between 0 and 1. \code{NA} is treated as
#' if \code{fix_C3grasses} is \code{FALSE}.
#' @param fix_shrubs A logical value. If \code{TRUE}, then value for the shrub
#' component is fixed at \code{Shrubs_Fraction} instead of calculated from
#' climatic relationships.
#' @param Shrubs_Fraction A numeric value between 0 and 1. \code{NA} is treated
#' as if \code{fix_shrubs} is \code{FALSE}.
#' @param fix_forbs A logical value. If \code{TRUE}, then value for the forb
#' component is fixed at \code{Forbs_Fraction}.
#' @param Forbs_Fraction A numeric value between 0 and 1. Default value is 0.
#' @param fix_trees A logical value. If \code{TRUE}, then value for the tree
#' component is fixed at \code{Trees_Fraction}.
#' @param Trees_Fraction A numeric value between 0 and 1. Default value is 0.
#' @param fix_BareGround A logical value. If \code{TRUE}, then value for the
#' bare ground component is fixed at \code{BareGround_Fraction}.
#' @param BareGround_Fraction A numeric value between 0 and 1.
#' Default value is 0.
#' @param fill_empty_with_BareGround A logical value. If \code{TRUE}, then
#' incomplete land cover is considered (additional) bare-ground. If
#' \code{FALSE}, then some hacks are used to "fill in" incomplete land cover
#' with grasses and/or shrubs, and additionally,
#' if \code{fix_BareGround} is \code{FALSE} with bare-ground.
#' @param warn_extrapolation A logical value. If \code{TRUE}, then
#' warnings are issued if climate inputs \code{MAP_mm} and/or \code{MAT_C}
#' fall outside the range of supported values. See notes.
#'
#' @return A list with three named numeric vectors. \describe{
#' \item{Rel_Abundance_L0}{A numeric vector of length 8 with
#' relative abundance/cover [0-1] values of land cover types that sum to 1.
#' The names of the 8 types are: \var{Succulents}, \var{Forbs},
#' \var{Grasses_C3}, \var{Grasses_C4}, \var{Grasses_Annuals},
#' \var{Shrubs}, \var{Trees}, \var{BareGround}.
#' }
#' \item{Rel_Abundance_L1}{A numeric vector of length 5 with
#' relative abundance/cover [0-1] values of \pkg{rSOILWAT2} land cover
#' types that sum to 1.
#' The names of the 5 types are: \var{SW_TREES}, \var{SW_SHRUBS},
#' \var{SW_FORBS}, \var{SW_GRASS}, and \var{SW_BAREGROUND}.
#' }
#' \item{Grasses}{A numeric vector of length 3 with
#' relative abundance/cover [0-1] values of the grass types that sum to 1,
#' if there is any grass cover; otherwise, the values are 0.
#' The names of the 3 sub-types are: \var{Grasses_C3}, \var{Grasses_C4},
#' and \var{Grasses_Annuals}.
#' }
#' }
#'
#' @references Paruelo J.M., Lauenroth W.K. (1996) Relative abundance of plant
#' functional types in grasslands and shrublands of North America. Ecological
#' Applications, 6, 1212-1224.
#' @references Teeri J.A., Stowe L.G. (1976) Climatic patterns and the
#' distribution of C4 grasses in North America. Oecologia, 23, 1-12.
#'
#' @examples
#' ## Load weather dataset from rSOILWAT2
#' utils::data("weatherData", package = "rSOILWAT2")
#' clim1 <- calc_SiteClimate(weatherList = weatherData)
#' clim2 <- calc_SiteClimate(weatherList = weatherData, do_C4vars = TRUE)
#'
#' ## All estimable vegetation types are estimated:
#' estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim1[["MAP_cm"]],
#' MAT_C = clim1[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim1[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim1[["meanMonthlyTempC"]]
#' )
#'
#' ## Climate is outside supported range with MAT < 0 C:
#' estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim1[["MAP_cm"]],
#' MAT_C = clim1[["MAT_C"]] - clim1[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim1[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim1[["meanMonthlyTempC"]] - clim1[["MAT_C"]]
#' )
#'
#' ## Some land cover types are fixed and others are estimated, and
#' ## the C4-grass adjustment is used:
#' estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim2[["MAP_cm"]],
#' MAT_C = clim2[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim2[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim2[["meanMonthlyTempC"]],
#' dailyC4vars = clim2[["dailyC4vars"]],
#' fix_shrubs = TRUE, Shrubs_Fraction = 0.5,
#' fix_BareGround = TRUE, BareGround_Fraction = 0.25
#' )
#'
#' ## Fix total grass cover and annual grass cover,
#' ## but estimate relative proportions of C3 and C4 grasses:
#' estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim2[["MAP_cm"]],
#' MAT_C = clim2[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim2[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim2[["meanMonthlyTempC"]],
#' dailyC4vars = clim2[["dailyC4vars"]],
#' fix_sumgrasses = TRUE,
#' SumGrasses_Fraction = 0.8,
#' fix_annuals = TRUE,
#' Annuals_Fraction = 0.3
#' )
#'
#'
#' ## SOILWAT2 uses the same algorithm internally if requested to do so
#' # Obtain cover values from SOILWAT2 output
#' swin <- rSOILWAT2::sw_exampleData
#' swin@prod@veg_method <- 1L
#' swout <- sw_exec(swin)
#' tmp <- slot(slot(swout, "BIOMASS"), "Year")
#' pnvsim <- tmp[1, grep("fCover", colnames(tmp)), drop = TRUE]
#'
#' # Directly calculate cover values
#' climex <- calc_SiteClimate(weatherList = get_WeatherHistory(swin))
#' pnvex <- estimate_PotNatVeg_composition(
#' MAP_mm = 10 * climex[["MAP_cm"]],
#' MAT_C = climex[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * climex[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = climex[["meanMonthlyTempC"]]
#' )[["Rel_Abundance_L1"]]
#'
#' # They are identical
#' identical(pnvsim[["fCover_shrub"]], pnvex[["SW_SHRUB"]])
#' identical(pnvsim[["fCover_grass"]], pnvex[["SW_GRASS"]])
#' identical(pnvsim[["fCover_forbs"]], pnvex[["SW_FORBS"]])
#' identical(pnvsim[["fCover_tree"]], pnvex[["SW_TREES"]])
#' identical(pnvsim[["fCover_BareGround"]], pnvex[["SW_BAREGROUND"]])
#'
#'
#' @export
estimate_PotNatVeg_composition <- function(MAP_mm, MAT_C,
mean_monthly_ppt_mm, mean_monthly_Temp_C, dailyC4vars = NULL,
isNorth = TRUE, shrub_limit = 0.2,
fix_succulents = FALSE, Succulents_Fraction = NA,
fix_sumgrasses = FALSE, SumGrasses_Fraction = NA,
fix_annuals = TRUE, Annuals_Fraction = 0,
fix_C4grasses = FALSE, C4_Fraction = NA,
fix_C3grasses = FALSE, C3_Fraction = NA,
fix_shrubs = FALSE, Shrubs_Fraction = NA,
fix_forbs = FALSE, Forbs_Fraction = NA,
fix_trees = TRUE, Trees_Fraction = 0,
fix_BareGround = TRUE, BareGround_Fraction = 0,
fill_empty_with_BareGround = TRUE,
warn_extrapolation = TRUE) {
res <- .Call(C_rSW2_estimate_PotNatVeg_composition,
MAP_mm,
MAT_C,
mean_monthly_ppt_mm,
mean_monthly_Temp_C,
shrub_limit,
if (fix_sumgrasses) {
as.numeric(SumGrasses_Fraction)
} else {
NA_real_
},
as.logical(fill_empty_with_BareGround),
as.logical(warn_extrapolation),
if (is.null(dailyC4vars) || anyNA(dailyC4vars[1:3])) {
NULL
} else {
as.numeric(dailyC4vars[1:3])
},
as.logical(isNorth),
as.logical(fix_BareGround),
if (fix_succulents) {
as.numeric(Succulents_Fraction)
} else {
NA_real_
},
if (fix_annuals) {
as.numeric(Annuals_Fraction)
} else {
NA_real_
},
if (fix_C4grasses) {
as.numeric(C4_Fraction)
} else {
NA_real_
},
if (fix_C3grasses) {
as.numeric(C3_Fraction)
} else {
NA_real_
},
if (fix_shrubs) {
as.numeric(Shrubs_Fraction)
} else {
NA_real_
},
if (fix_forbs) {
as.numeric(Forbs_Fraction)
} else {
NA_real_
},
if (fix_trees) {
as.numeric(Trees_Fraction)
} else {
NA_real_
},
if (fix_BareGround) {
as.numeric(BareGround_Fraction)
} else {
NA_real_
}
)
res
}
#' Adjust a phenological pattern to a different temperature regime
#'
#' Extract the characteristics of seasonal (mean monthly) values of biomass,
#' activity, etc. (i.e., phenology) that were observed for a specific
#' (reference) climate (at a specific site), and project those characteristics
#' to the mean monthly temperature values of a different site and/or different
#' climate conditions.
#'
#' @param x A numeric vector of length 12. The values
#' reflect the phenological pattern that occur on average under
#' reference mean monthly temperature values \code{ref_temp}.
#' @param ref_temp A numeric vector of length 12. Reference mean monthly
#' temperature values in degree Celsius under which \code{x} was
#' determined / is valid.
#' @param target_temp A numeric vector of length 12. Mean monthly
#' temperature values in degree Celsius of a target site / condition
#' for which \code{x} is to be adjusted.
#' @param x_asif A numeric vector of length 12. If not \code{NULL}, then
#' the adjustments applied to \code{x} are based on values of \code{x_asif}.
#'
#' @return A copy of \code{x} with values that have been adjusted
#' to represent the phenology of a target climate described by
#' \code{target_temp}.
#'
#' @section Details: The implementation takes a geometrical interpretation to
#' the seasonal cycle.
#' The time axis is represented by a standardized twelve month circle;
#' temperature and vegetation values are standardized to lie within a band along
#' the time circle;
#' and warm- and cold-seasons are presumed to exist (i.e., temperate
#' climates will work best) and are represented by ellipses with centers at the
#' time of peak of temperature values (warm-season) below the standardized band
#' of values and half a year earlier (cold-season) above that band.
#' From these centers, radial expansion or contraction rates are calculated
#' as vectors from reference towards target temperatures and these vectors
#' are applied as adjustment rates to vegetation values.
#' Adjusted values are destandardized before returning.
#'
#' @section Details: The purpose of \code{x_asif} is to adjust \code{x} based
#' on the phenological pattern of \code{x_asif} instead of \code{x}. This
#' may be useful, for instance, if \code{x} represents total biomass and
#' \code{x_asif} fraction of live biomass and if we want the adjustments to
#' reflect that the live biomass fraction predominantly influences the
#' total biomass response to different temperatures.
#' @section Details: Adjustments based on \code{x_asif} reflect the
#' amount of variation in \code{x} relative to the amount of variation in
#' \code{x_asif}. Thus, the less variation in \code{x} relative to
#' \code{x_asif}, the less \code{x} is adjusted; and if \code{x} does not
#' vary across months, then \code{x} is not adjusted at all regardless of the
#' values of \code{x_asif}. This may not be desired in all cases.
#'
#' @section Notes: Phenological adjustments applied to \code{x} will not
#' necessarily maintain any specific characteristic of \code{x}, e.g.,
#' values in \code{[0, 1]}, \code{sum(x) == 1}, \code{max(x) == peak}.
#' Thus, post-adjustment scaling of the returned values may be required. See
#' examples.
#'
#'
#' @examples
#' sw_in <- rSOILWAT2::sw_exampleData
#' tmp <- swProd_MonProd_grass(sw_in)
#' phen_reference <- data.frame(
#' tmp,
#' Litter_pct = tmp[, "Litter"] / max(tmp[, "Litter"]),
#' Biomass_pct = tmp2 <- tmp[, "Biomass"] / max(tmp[, "Biomass"]),
#' Biomass2_pct = {tmp2[6:8] <- 1; tmp2}
#' )
#'
#' clim <- calc_SiteClimate(weatherList = rSOILWAT2::weatherData)
#' ref_temp <- clim[["meanMonthlyTempC"]]
#' temp_randomwarmer3C <- ref_temp + stats::rnorm(12, 3, 1)
#'
#' phen_adj <- sapply(
#' phen_reference,
#' adj_phenology_by_temp,
#' ref_temp = ref_temp,
#' target_temp = temp_randomwarmer3C,
#' x_asif = phen_reference[, "Live_pct"]
#' )
#'
#' ## Note: depending on what `x` represents, post-adjustment scaling may
#' ## be necessary
#'
#' # Maintain previous peak
#' rSW2utils::scale_to_reference_fun(
#' x = phen_adj[, "Litter_pct"],
#' x_ref = phen_reference[, "Litter_pct"],
#' fun = max
#' )
#'
#' # Maintain previous sum
#' rSW2utils::scale_to_reference_fun(
#' x = phen_adj[, "Biomass"],
#' x_ref = phen_reference[, "Biomass"],
#' fun = sum
#' )
#'
#' # Maintain previous frequency of peaks and cap the peaks at that value
#' rSW2utils::scale_to_reference_peak_frequency(
#' x = phen_adj[, "Biomass2_pct"],
#' x_ref = phen_reference[, "Biomass2_pct"],
#' cap_at_peak = TRUE
#' )
#'
#' ## Plot reference and adjusted monthly values
#' if (interactive()) {
#' par_prev <- par(mfrow = c(2, 1))
#'
#' plot(
#' 1:12,
#' phen_reference[, 1],
#' type = "l",
#' col = "red",
#' ylim = c(0, 300),
#' xlab = "Month"
#' )
#' lines(1:12, phen_adj[, 1])
#' legend(
#' "bottom",
#' ncol = 2,
#' legend = c("reference", "adjusted"),
#' lwd = 2, lty = 1, col = c("red", "black")
#' )
#'
#' plot(
#' 1:12,
#' phen_reference[, 2],
#' type = "l",
#' col = "red",
#' ylim = c(0, 1),
#' xlab = "Month"
#' )
#' lines(1:12, phen_adj[, 2])
#'
#' par(par_prev)
#' }
#'
#' @export
adj_phenology_by_temp <- function(x, ref_temp, target_temp, x_asif = NULL) {
mon <- rSW2_glovars[["st_mo"]]
nmon <- 12 # number of months `mon`
nwks <- 53
tol <- rSW2_glovars[["tol"]]
tol2 <- sqrt(tol)
if (is.null(x_asif)) {
x_asif <- x
}
stopifnot(
length(x) == nmon,
length(x_asif) == nmon,
length(ref_temp) == nmon,
length(target_temp) == nmon
)
#--- Free parameters
Nr <- 2 * nwks # number of steps for adjustments (1 full circle)
Ninfl <- 19 #17 # number of `nwks` units of seasonal influence
p_voffset <- 1 / 4 # vertical offset of season centers beyond 0-0.5 values
#------ Step 1) Define cold- and warm-season circles/ellipses
yr_std_d <- 1 # duration of 1 year = diameter
yr_std_r <- yr_std_d / 2 # radius of year
yr_ctrs <- c(0, yr_std_r) # time (x-axis) centers of cold-/warm-seasons
vals_std <- c(0, yr_std_r) # y-axis band of cold-/warm-season activity
vals_std_off <- vals_std[2] * p_voffset # y-axis offset of season centers
infl_warm <- Ninfl / nwks
infl_cold <- Ninfl / nwks
# Cold-season ellipse with center (0, 1), radii 0.5/1,
# anti-clockwise direction, and start at (-0.5, -0.5)
ctr_cold <- c(yr_ctrs[1], max(vals_std) + vals_std_off)
fun_cold <- rSW2utils::f_ellipse(
x0 = ctr_cold[1],
y0 = ctr_cold[2],
a = yr_std_r,
b = vals_std[2] + vals_std_off,
theta0 = pi
)
# Warm-season ellipse with center (0.5, -0.5), radii 0.5/1,
# anti-clockwise direction, and start at (0, 0)
ctr_warm <- c(yr_ctrs[2], min(vals_std) - vals_std_off)
fun_warm <- rSW2utils::f_ellipse(
x0 = ctr_warm[1],
y0 = ctr_warm[2],
a = yr_std_r,
b = vals_std[2] + vals_std_off,
theta0 = pi,
dir = -1
)
# Identify which points fall to the influence of cold- vs. warm-ellipses
# (with some overlap)
xrange_warm <- ctr_warm[1] + c(-1, 1) * infl_warm
xrange_warm2 <- ctr_warm[1] + c(-1, 1) * infl_warm / 2
xrange_cold <- ctr_cold[1] + c(-1, 1) * infl_cold
xrange_cold2 <- ctr_cold[1] + c(-1, 1) * infl_cold / 2
# Season sequences
# - must be continuous for "continuous adjustment sections"
# - make sure season ellipses' vertices are included in sequence
# - increase density of points near season ellipses' vertices
# (where uniroot tends to fail and vegetation curve often changes fast)
tmp <- pi / 2
seq_angle1 <- sort(unique(c(
c(-1, 1) * tmp,
seq_len(Nr) / Nr * 2 * pi
)))
seq_angle2 <- seq_len(2 * Nr) / Nr * pi
seq_cold1 <- fun_cold(seq_angle1)
ids_cold1 <- which(
seq_cold1[, 2] <= ctr_cold[2] &
seq_cold1[, 1] >= xrange_cold[1] &
seq_cold1[, 1] <= xrange_cold[2]
)
seq_cold2 <- fun_cold(seq_angle2)
ids_cold2 <- which(
seq_cold2[, 2] <= ctr_cold[2] &
seq_cold2[, 1] >= xrange_cold2[1] &
seq_cold2[, 1] <= xrange_cold2[2]
)
seq_cold <- unique(rbind(seq_cold1[ids_cold1, ], seq_cold2[ids_cold2, ]))
seq_cold <- seq_cold[order(seq_cold[, 1]), ]
seq_warm1 <- fun_warm(seq_angle1)
ids_warm1 <- which(
seq_warm1[, 2] >= ctr_warm[2] &
seq_warm1[, 1] >= xrange_warm[1] &
seq_warm1[, 1] <= xrange_warm[2]
)
seq_warm2 <- fun_warm(seq_angle2)
ids_warm2 <- which(
seq_warm2[, 2] >= ctr_warm[2] &
seq_warm2[, 1] >= xrange_warm2[1] &
seq_warm2[, 1] <= xrange_warm2[2]
)
seq_warm <- unique(rbind(seq_warm1[ids_warm1, ], seq_warm2[ids_warm2, ]))
seq_warm <- seq_warm[order(seq_warm[, 1]), ]
n_cold <- nrow(seq_cold)
n_warm <- nrow(seq_warm)
Nadj <- n_cold + n_warm
# Information about the points that define the lines which are used
# to intersect values for adjustments
# columns:
# 1:2 = x/y-coords of center of cold-/warm-season ellipses
# 3:4 = x/y-coords of points along cold-/warm-season ellipses
# 5:6 = min/max of x-coord under influence of cold-/warm-season ellipses
# 7 = x-extent of influence around center of cold-/warm-season ellipses
ptadj <- matrix(NA, nrow = Nadj, ncol = 7)
tmp <- seq_len(n_cold)
ptadj[tmp, 1:2] <- rep(ctr_cold, each = n_cold)
ptadj[tmp, 3:4] <- seq_cold
ptadj[tmp, 5:6] <- rep(xrange_cold, each = n_cold)
ptadj[tmp, 7] <- infl_cold
tmp <- n_cold + seq_len(n_warm)
ptadj[tmp, 1:2] <- rep(ctr_warm, each = n_warm)
ptadj[tmp, 3:4] <- seq_warm
ptadj[tmp, 5:6] <- rep(xrange_warm, each = n_warm)
ptadj[tmp, 7] <- infl_warm
#------ Step 2) Normalize (shift and scale) all values
# so that they fit into the season ellipses
# Normalize time of year values (on x-axis) into [0, 1]
mon_norm <- mon / nmon
# Normalize temperature values (on y-axis) into [0, 0.5]
temp_zero <- -min(ref_temp, target_temp)
tmp1 <- temp_zero + ref_temp
tmp2 <- temp_zero + target_temp
temp_2max <- max(tmp1, tmp2) / vals_std[2]
ref_temp_norm <- tmp1 / temp_2max
ref_temp_norm_mean <- mean(ref_temp_norm)
target_temp_norm <- tmp2 / temp_2max
# Normalize vegetation values (on y-axis) into [0, 0.5] so that mean = 0.25
veg_min <- c(0, 0)
veg_2max <- (c(max(x), max(x_asif)) - veg_min) / vals_std[2]
veg_4mean <- (c(mean(x), mean(x_asif)) - veg_min) / (vals_std[2] / 2)
veg_scale <- pmax(veg_2max, veg_4mean)
veg_norm <- (x - veg_min[1]) / veg_scale[1]
veg_asif_norm <- (x_asif - veg_min[2]) / veg_scale[2]
# Calculate factor to correct for different spreads between
# vegetation data and `x_asif`
tmp1 <- diff(range(x_asif))
if (abs(tmp1) < rSW2_glovars[["tol"]]) {
# if diff(range(x_asif)) -> 0,
# then max(x_asif) / diff(range(x_asif)) -> Inf
# ==> cap ratio at one to avoid overcorrection
cvtoasif <- 1
} else {
# if diff(range(x)) -> 0, then diff(range(x)) / max(x) -> 0
tmp2 <- max(x)
if (abs(tmp2) < rSW2_glovars[["tol"]]) {
cvtoasif <- 0
} else {
# cap at one to avoid overcorrection
cvtoasif <- rSW2utils::finite01(
diff(range(x)) / tmp2 * max(x_asif) / tmp1
)
}
}
# Shift temperature and vegetation values (on x-axis)
# so that their ellipses are centered on the warm- and cold-season centers
# i.e., shift time of year so that peaks lie at warm-season x-center of 0.5
# Original centers = time of peak season as mean month weighted by values
tmp <- order(ref_temp_norm, decreasing = TRUE)
tpeak_norm <- rSW2utils::circ_mean_weighted(
x = tmp,
w = ref_temp_norm[tmp],
int = nmon,
type = "ZeroPlus2Pi"
) / nmon
sx <- yr_std_r - tpeak_norm
mon_norm_center <- mon_norm + sx
#------ Step 2) Fit periodic cubic splines to normalized values
psps <- list(
t_ref = splines::periodicSpline(
ref_temp_norm ~ mon_norm_center,
period = yr_std_d
),
t_target = splines::periodicSpline(
target_temp_norm ~ mon_norm_center,
period = yr_std_d
),
veg = splines::periodicSpline(
veg_norm ~ mon_norm_center,
period = yr_std_d
),
veg_asif = splines::periodicSpline(
veg_asif_norm ~ mon_norm_center,
period = yr_std_d
)
)
itcurves <- 1:2
ivcurves <- 3:4
#------ Step 3) Correct vegetation
# based on radial temperature expansion/contraction
# Prepare container for adjusted (normalized) vegetation points (xy)
vadj <- vector("list", length = Nadj)
k_vadj <- 1
ksc <- 1 # counter of continous adjustment sections
kce <- c(1, 0, 0) # counter of continuous failure sections
# Loop through points along season-ellipses (`ptadj`)
for (k in seq_len(Nadj)) {
#--- Locate intersection between target temperature function and
# line connecting points (time of year, reference temperature function) and
# season center points
if (abs(ptadj[k, 3] - ptadj[k, 1]) > rSW2_glovars[["tol"]]) {
# Function describing line that is connecting the circle-center point
# with the k-th point on the season-circle
fline_radial <- rSW2utils::f_2pline(
x0 = ptadj[k, 1], y0 = ptadj[k, 2],
x1 = ptadj[k, 3], y1 = ptadj[k, 4]
)
# Search limits for locating intersections
if (ptadj[k, 3] > ptadj[k, 1]) {
# limits: ellipse center to max influence
olims <- c(ptadj[k, 1], ptadj[k, 6])
} else {
# limits: min influence to ellipse center
olims <- c(ptadj[k, 5], ptadj[k, 1])
}
olims <- rSW2utils::extend_range(olims, 1 + ptadj[k, 7] / 2)
# Locate the x-axis value where line crosses
# (i) normalized reference temperature curve
# (ii) normalized target temperature curve
xcrosses_tcurves <- lapply(
X = itcurves,
FUN = function(itp) {
rSW2utils::uniroots(
f = function(x) {
predict(psps[[itp]], x)[["y"]] - fline_radial(x)[, 2]
},
xlim = olims,
tol = tol,
expected_nroots = 1
)
}
)
# Success of detection of temperature curves intersections
good_tcurves <- sapply(
xcrosses_tcurves,
rSW2utils::has_uniroots,
tol = tol2
)
if (all(good_tcurves)) {
# Locate (all) x-axis values where line crosses
# (iii) normalized vegetation curves
xcrosses_vcurves <- lapply(
X = ivcurves,
FUN = function(ivp) {
rSW2utils::uniroots(
f = function(x) {
predict(psps[[ivp]], x)[["y"]] - fline_radial(x)[, 2]
},
xlim = olims,
tol = tol,
expected_nroots = 100
)
}
)
# Success of detection of vegetation curve intersections
good_vcurves <- lapply(
xcrosses_vcurves,
rSW2utils::has_uniroots,
tol = tol2
)
if (!identical(good_vcurves[[1]], good_vcurves[[2]])) {
good_vcurves <- FALSE
}
} else {
good_vcurves <- FALSE
}
} else {
# Vertical line (`f_2pline` fails because slope would be inf)
tmp <- list(list(root = ptadj[k, 3]))
xcrosses_tcurves <- lapply(seq_along(itcurves), function(...) tmp)
good_tcurves <- TRUE
xcrosses_vcurves <- lapply(seq_along(ivcurves), function(...) tmp)
good_vcurves <- lapply(seq_along(ivcurves), function(...) TRUE)
}
#--- Calculate and apply radial expansion/contraction
if (all(good_tcurves) && all(sapply(good_vcurves, any))) {
# Intersection points: reference temp, target temp
x_tvals <- sapply(xcrosses_tcurves, function(x) x[[1]][["root"]])
y_tvals <- sapply(
seq_along(x_tvals),
function(i) predict(psps[[i]], x_tvals[i])[["y"]]
)
# Calculate radial x- and y-expansion/contraction
# (from ref temps to target temps)
# as fraction of reference temps to season center
rx <- (x_tvals[2] - x_tvals[1]) / (x_tvals[1] - ptadj[k, 1])
rx <- if (is.finite(rx)) rx else 0
ry <- (y_tvals[2] - y_tvals[1]) / (y_tvals[1] - ptadj[k, 2])
ry <- if (is.finite(ry)) ry else 0
# Intersection points: vegetation(s)
x_vvals <- lapply(
seq_along(ivcurves),
function(kvc) {
sapply(
xcrosses_vcurves[[kvc]][good_vcurves[[kvc]]],
function(x) x[["root"]]
)
}
)
y_vvals <- lapply(
seq_along(ivcurves),
function(kvc) {
predict(psps[[ivcurves[kvc]]], x_vvals[[kvc]])[["y"]]
}
)
# Calculate adjusted vegetation point(s)
# based on `x_asif` (index 2) but applied to `x` (index 1)
for (kv in seq_along(x_vvals[[2]])) {
vadj[[k_vadj]] <- c(
# xy-coordinates
x_vvals[[1]][kv] + rx * cvtoasif * (x_vvals[[2]][kv] - ptadj[k, 1]),
y_vvals[[1]][kv] + ry * cvtoasif * (y_vvals[[2]][kv] - ptadj[k, 2]),
# season that adjusts
ptadj[k, 1],
# number of continous adjustment section
ksc
)
k_vadj <- k_vadj + 1
}
} else {
# discontinuation in adjustment sections
ksc <- ksc + 1
# count longest continuous failure section
if (kce[1] == k - 1) {
kce[1:2] <- c(k, kce[2] + 1)
} else {
if (kce[2] > kce[3]) {
kce[3] <- kce[2]
}
kce[1:2] <- c(k, 1)
}
}
}
if (ksc > 0.2 * Nadj) {
warning("More than 20% of adjustment points failed.")
}
if (max(kce[2:3]) > 0.1 * Nadj) {
warning(
"The longest contiguous section of failed adjustments ",
"is longer than 10% of all points.")
}
#--- Pre-process adjusted vegetation values separately for seasons
# (in case certain time points received both cold- and warm-season adj.)
# Clean, sort, and separate adjusted vegetation values
vadj <- vadj[seq_len(k_vadj - 1)]
vadj2 <- matrix(unlist(vadj), nrow = length(vadj), ncol = 4, byrow = TRUE)
vadj2 <- vadj2[complete.cases(vadj2), ]
# determine duplicated values: see `smooth.spline`
tmp <- vadj2[, 1]
isdup <- duplicated(round((tmp - mean(tmp)) / (1e-6 * stats::IQR(tmp))))
vadj2 <- vadj2[!isdup, ]
vadj2 <- vadj2[order(vadj2[, 1]), , drop = FALSE]
vadj2 <- lapply(
c(ctr_cold[1], ctr_warm[1]),
function(ctr) vadj2[abs(vadj2[, 3] - ctr) < tol, , drop = FALSE]
)
# Continous adjustment sections
scs <- lapply(vadj2, function(x) sort(unique(x[, 4])))
# Fit periodic linear spline to adjusted vegetation values per season
# and predict along regular sequence (for each continuous section)
fxpsp_adj_norm2 <- lapply(
vadj2,
function(x) {
if (nrow(x) >= 4) {
splines::interpSpline(x[, 2] ~ x[, 1])
}
}
)
itmp <- seq_len(Nadj) / Nadj - 1 / (2 * Nadj)
vadj3 <- lapply(
seq_along(vadj2),
function(k2) {
tmp <- matrix(NA, nrow = Nadj, ncol = 1 + length(scs[[k2]]))
tmp[, 1] <- itmp
tmp
}
)
# Columns: (1) x-axis location, (2) cold-season adj. veg,
# (3) warm-season adj. veg, (4) across season averaged adj. veg
vadj4 <- matrix(NA, nrow = Nadj, ncol = 4)
vadj4[, 1] <- itmp
# Predict seasonaly adjusted vegetation for each continous section
for (k2 in seq_along(vadj2)) {
if (!is.null(fxpsp_adj_norm2[[k2]])) {
for (k3 in seq_along(scs[[k2]])) {
ids <- vadj2[[k2]][, 4] %in% scs[[k2]][k3]
xl <- round(Nadj * range(vadj2[[k2]][ids, 1]))
xt2 <- seq(xl[1], xl[2])
tmp <- predict(
fxpsp_adj_norm2[[k2]],
xt2 / Nadj - 1 / (2 * Nadj)
)
tmp <- aggregate(
tmp[["y"]],
by = list(tmp[["x"]] %% yr_std_d),
FUN = mean
)
xt2p <- sapply(
tmp[, 1],
function(x) which.min(abs(x - vadj3[[k2]][, 1] %% yr_std_d))
)
vadj3[[k2]][xt2p, 1 + k3] <- tmp[, 2]
}
# Average across continous sections (within season)
vadj4[, 1 + k2] <- rowMeans(vadj3[[k2]][, -1, drop = FALSE], na.rm = TRUE)
}
}
# Average values across seasons
vadj4[, 4] <- rowMeans(vadj4[, -1], na.rm = TRUE)
isgood <- complete.cases(vadj4[, c(1, 4)])
vadj4 <- vadj4[isgood, , drop = FALSE]
#------ Step 4) Fit periodic linear spline to adjusted vegetation values
# De-normalize vegetation values back to original scale
veg_new <- veg_min[1] + veg_scale[1] * vadj4[, 4]
# Revert temporal shifts due to centering of refs, target temps and veg
toy_new <- (vadj4[, 1] - sx) %% yr_std_d
if (nrow(vadj4) >= 2) {
# Use linear (instead of cubic) spline to smooth across potential steps
# between cold- and warm-season adjustments
fxpsp_adj <- splines::periodicSpline(
veg_new ~ toy_new,
period = yr_std_d,
ord = 2 # linear spline
)
# Obtain monthly adjusted vegetation values
pres <- predict(fxpsp_adj, mon_norm)[["y"]]
} else {
warning("Less than 3 adjustment points succeeded: return mean value.")
pres <- rep(mean(veg_new), nmon)
}
pmax(0, pres)
}
#' Biomass equations
#'
#' @param MAP_mm A numeric vector. Mean annual precipitation in
#' millimeters (mm).
#' @references Milchunas & Lauenroth 1993 (Fig. 2):
#' Y [g/m2/yr] = c1 * MAP [mm/yr] + c2
#' @name biomass
NULL
#' Estimate shrub biomass density from mean annual precipitation
#' @export
#' @rdname biomass
Shrub_ANPP <- function(MAP_mm) 0.393 * MAP_mm - 10.2
#' Estimate grass biomass density from mean annual precipitation
#' @export
#' @rdname biomass
Grass_ANPP <- function(MAP_mm) 0.646 * MAP_mm - 102.5
#' Adjust mean monthly biomass values by precipitation
#' @section Details: Internally used by
#' \code{\link{estimate_PotNatVeg_biomass}}.
#' @noRd
adjBiom_by_ppt <- function(biom_shrubs, biom_C3, biom_C4, biom_annuals,
biom_maxs, map_mm_shrubs, map_mm_std_shrubs,
map_mm_grasses, map_mm_std_grasses,
vegcomp_std_shrubs, vegcomp_std_grass) {
# Intercepts to match outcomes of M & L 1993 equations under 'default'
# MAP with our previous default inputs for shrubs and sgs-grasslands
# Whereas these intercepts were introduced artificially, they could also be
# interpreted as perennial storage, e.g., Lauenroth & Whitman (1977) found
# "Accumulation in the standing dead was 63% of inputs, in the litter 8%,
# and belowground 37%.". Lauenroth, W.K. & Whitman, W.C. (1977)
# Dynamics of dry matter production in a mixed-grass prairie in western
# North Dakota. Oecologia, 27, 339-351.
Shrub_ANPPintercept <- (vegcomp_std_shrubs[1] * biom_maxs["Sh.Amount.Live"] +
vegcomp_std_shrubs[2] * biom_maxs["C3.Amount.Live"] +
vegcomp_std_shrubs[3] * biom_maxs["C4.Amount.Live"]) -
Shrub_ANPP(map_mm_std_shrubs)
Grasses_ANPPintercept <- (vegcomp_std_grass[1] * biom_maxs["Sh.Amount.Live"] +
vegcomp_std_grass[2] * biom_maxs["C3.Amount.Live"] +
vegcomp_std_grass[3] * biom_maxs["C4.Amount.Live"]) -
Grass_ANPP(map_mm_std_grasses)
# Get scaling values for scaled biomass; guarantee that > minimum.totalBiomass
minimum.totalBiomass <- 0 # This is a SOILWAT2 parameter
Shrub_BiomassScaler <- max(minimum.totalBiomass,
Shrub_ANPP(map_mm_shrubs) + Shrub_ANPPintercept)
Grass_BiomassScaler <- max(minimum.totalBiomass,
Grass_ANPP(map_mm_grasses) + Grasses_ANPPintercept)
# Scale live biomass amount by productivity; assumption:
# ANPP = peak standing live biomass
biom_shrubs$Sh.Amount.Live <- biom_shrubs$Sh.Amount.Live * Shrub_BiomassScaler
biom_C3$C3.Amount.Live <- biom_C3$C3.Amount.Live * Grass_BiomassScaler
biom_C4$C4.Amount.Live <- biom_C4$C4.Amount.Live * Grass_BiomassScaler
biom_annuals$Annual.Amount.Live <-
biom_annuals$Annual.Amount.Live * Grass_BiomassScaler
#Scale litter amount by productivity and adjust for ratio of litter/live
biom_shrubs$Sh.Litter <- biom_shrubs$Sh.Litter * Shrub_BiomassScaler *
biom_maxs["Sh.Litter"] / biom_maxs["Sh.Amount.Live"]
biom_C3$C3.Litter <- biom_C3$C3.Litter * Grass_BiomassScaler *
biom_maxs["C3.Litter"] / biom_maxs["C3.Amount.Live"]
biom_C4$C4.Litter <- biom_C4$C4.Litter * Grass_BiomassScaler *
biom_maxs["C4.Litter"] / biom_maxs["C4.Amount.Live"]
biom_annuals$Annual.Litter <- biom_annuals$Annual.Litter *
Grass_BiomassScaler *
biom_maxs["Annual.Litter"] / biom_maxs["Annual.Amount.Live"]
#Guarantee that live fraction = ]0, 1]
biom_shrubs$Sh.Perc.Live <- pmin(1, pmax(rSW2_glovars[["tol"]],
biom_shrubs$Sh.Perc.Live))
biom_C3$C3.Perc.Live <- pmin(1, pmax(rSW2_glovars[["tol"]],
biom_C3$C3.Perc.Live))
biom_C4$C4.Perc.Live <- pmin(1, pmax(rSW2_glovars[["tol"]],
biom_C4$C4.Perc.Live))
biom_annuals$Annual.Perc.Live <- pmin(1, pmax(rSW2_glovars[["tol"]],
biom_annuals$Annual.Perc.Live))
#Calculate total biomass based on scaled live biomass amount
biom_shrubs$Sh.Biomass <- biom_shrubs$Sh.Amount.Live /
biom_shrubs$Sh.Perc.Live
biom_C3$C3.Biomass <- biom_C3$C3.Amount.Live / biom_C3$C3.Perc.Live
biom_C4$C4.Biomass <- biom_C4$C4.Amount.Live / biom_C4$C4.Perc.Live
biom_annuals$Annual.Biomass <- biom_annuals$Annual.Amount.Live /
biom_annuals$Annual.Perc.Live
list(biom_shrubs = biom_shrubs,
biom_C3 = biom_C3,
biom_C4 = biom_C4,
biom_annuals = biom_annuals)
}
#' Adjust mean monthly biomass values of grass and shrub functional groups by
#' climate relationships
#'
#' @inheritParams adj_phenology_by_temp
#' @param target_MAP_mm A numeric value. Mean annual precipitation
#' in millimeter of the location.
#' @param tr_VegBiom A data.frame with 12 rows (one for each month) and columns
#' \code{X.Biomass}, \code{X.Amount.Live}, \code{X.Perc.Live}, and
#' \code{X.Litter} where \code{X} are for the functional groups shrubs,
#' \code{X = Sh}; C3-grasses, \code{X = C3}; C4-grasses, \code{X = C4}; and
#' annuals, \code{X = Annual} containing default input values.
#' Function default values are from Bradford et al. 2014, see
#' \code{\link{sw2_tr_VegBiom}}.
#' @param do_adjust_phenology A logical value. If \code{TRUE} then monthly
#' phenology is adjusted by temperature.
#' @param do_adjust_biomass A logical value. If \code{TRUE} then monthly biomass
#' is adjusted by precipitation.
#' @param fgrass_c3c4ann A numeric vector of length 3. Relative contribution
#' `[0-1]` of the C3-grasses, C4-grasses, and annuals functional groups.
#' The sum of `fgrass_c3c4ann` is 1.
#'
#' @section Default inputs: \itemize{
#' \item Shrubs are based on location \var{\sQuote{IM_USC00107648_Reynolds}}
#' which resulted in a vegetation composition of 70 % shrubs and 30 %
#' C3-grasses. Default monthly biomass values were estimated for
#' MAP = 450 mm yr-1.
#' \item Grasses are based on location \var{\sQuote{GP_SGSLTER}}
#' (shortgrass steppe) which resulted in 12 % shrubs, 22 % C3-grasses,
#' and 66 % C4-grasses. Default biomass values were estimated for
#' MAP = 340 mm yr-1.
#' \item Mean monthly reference temperature are the median values across
#' 898 big sagebrush sites
#' (see \url{https://github.com/DrylandEcology/rSFSTEP2/issues/195})
#' }
#'
#' @return A list with two elements \code{grass}, \code{shrub}. Each element is
#' a matrix with 12 rows (one for each month) and columns \code{Biomass},
#' \code{Amount.Live}, \code{Perc.Live}, and \code{Litter}.
#'
#' @section Details:
#' If `do_adjust_biomass`, then the internal function `adjBiom_by_ppt()` is
#' used to adjust biomass by annual precipitation amount.
#' If `do_adjust_phenology`, then the exported function
#' [adj_phenology_by_temp()] is used to adjust the seasonal pattern of biomass
#' (phenology) by monthly temperature.
#'
#' @references Bradford, J.B., Schlaepfer, D.R., Lauenroth, W.K. & Burke, I.C.
#' (2014). Shifts in plant functional types have time-dependent and regionally
#' variable impacts on dryland ecosystem water balance. J Ecol, 102,
#' 1408-1418.
#'
#' @examples
#' clim <- calc_SiteClimate(weatherList = rSOILWAT2::weatherData)
#'
#' veg_cover <- rSOILWAT2::estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim[["MAP_cm"]],
#' MAT_C = clim[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim[["meanMonthlyTempC"]],
#' dailyC4vars = clim[["dailyC4vars"]]
#' )
#'
#' rSOILWAT2::estimate_PotNatVeg_biomass(
#' target_temp = clim[["meanMonthlyTempC"]],
#' target_MAP_mm = 10 * clim[["MAP_cm"]],
#' do_adjust_phenology = TRUE,
#' do_adjust_biomass = TRUE,
#' fgrass_c3c4ann = veg_cover[["Grasses"]]
#' )
#'
#' @export
#' @md
estimate_PotNatVeg_biomass <- function(
target_temp,
target_MAP_mm,
ref_temp,
tr_VegBiom = rSOILWAT2::sw2_tr_VegBiom,
do_adjust_phenology = FALSE,
do_adjust_biomass = FALSE,
fgrass_c3c4ann = c(1, 0, 0)
) {
if (missing(ref_temp) || is.null(ref_temp)) {
# Mean monthly reference temperature
# corresponding to default phenology values:
# median values across 898 big sagebrush sites
# (see https://github.com/DrylandEcology/rSFSTEP2/issues/195)
ref_temp <- c(
-4.6768, -2.7282, 1.8257, 6.0538, 10.696, 15.3878,
19.7777, 18.8755, 13.7868, 7.2843, 0.4167, -4.6912
)
}
if (missing(target_MAP_mm) || is.null(target_MAP_mm)) {
target_MAP_mm <- 450
}
ns_VegBiom <- names(tr_VegBiom)
tmp <- strsplit(ns_VegBiom, split = ".", fixed = TRUE)
vtypes <- unique(sapply(tmp, function(x) x[[1]]))
stopifnot(identical(vtypes, c("Sh", "C3", "C4", "Annual")))
btypes <- unique(sapply(tmp, function(x) paste(x[-1], collapse = ".")))
stopifnot(identical(btypes, c("Litter", "Biomass", "Perc.Live")))
btypes2 <- c(btypes, "Amount.Live")
# Default shrub biomass input is at MAP = 450 mm/yr, and default grass
# biomass input is at MAP = 340 mm/yr
#Default site for the grass description is SGS LTER
StandardGrasses_MAP_mm <- 340
StandardGrasses_VegComposition <- c(0.12, 0.22, 0.66) # shrubs, C3, and C4
#Default site for the shrub description is Reynolds Creek, ID
StandardShrub_MAP_mm <- 250
StandardShrub_VegComposition <- c(0.7, 0.3, 0) # shrubs, C3, and C4
# Scale monthly values of litter and biomass amount by column-max
tmp <- apply(
tr_VegBiom[, grepl("Biomass", ns_VegBiom, fixed = TRUE)] *
tr_VegBiom[, grepl("Perc.Live", ns_VegBiom, fixed = TRUE)],
MARGIN = 2,
FUN = max
)
names(tmp) <- paste0(vtypes, ".Amount.Live")
colmax <- c(
apply(
tr_VegBiom,
MARGIN = 2,
FUN = max
),
tmp
)
# Pull vegetation types
x0 <- lapply(
vtypes,
function(vt) tr_VegBiom[, grepl(vt, ns_VegBiom)]
)
names(x0) <- paste0("biom_", vtypes)
# adjust phenology for mean monthly temperatures
if (do_adjust_phenology) {
x <- lapply(
X = x0,
FUN = function(vt_data) {
sapply(
X = vt_data,
adj_phenology_by_temp,
ref_temp = ref_temp,
target_temp = target_temp,
x_asif = vt_data[, grepl("Perc.Live", colnames(vt_data))]
)
}
)
# Scale to previous peak frequency (0-1)
for (ns in names(x0)) {
for (k in seq_len(ncol(x[[ns]]))) {
x[[ns]][, k] <- rSW2utils::scale_to_reference_peak_frequency(
x = x[[ns]][, k],
x_ref = x0[[ns]][, k],
cap_at_peak = FALSE
)
}
}
} else {
x <- x0
}
for (vt in vtypes) {
ns <- paste0("biom_", vt)
ntmp <- colnames(x[[ns]])
x[[ns]] <- as.data.frame(x[[ns]])
# Calculate 'live biomass amount'
x[[ns]][, paste0(vt, ".Amount.Live")] <-
x[[ns]][, grepl("Biomass", ntmp, fixed = TRUE)] *
x[[ns]][, grepl("Perc.Live", ntmp, fixed = TRUE)]
# Scale so that 1 = previous peak (as required by `adjBiom_by_ppt`)
for (bt in btypes2[-3]) {
ctmp <- paste0(vt, ".", bt)
x[[ns]][, ctmp] <- x[[ns]][, ctmp] / colmax[ctmp]
}
}
# if (do_adjBiom_by_ppt) then adjust biomass amounts by productivity
# relationship with MAP
# total biomass is back calculated
x <- adjBiom_by_ppt(
biom_shrubs = x[["biom_Sh"]],
biom_C3 = x[["biom_C3"]],
biom_C4 = x[["biom_C4"]],
biom_annuals = x[["biom_Annual"]],
biom_maxs = colmax,
map_mm_shrubs = if (do_adjust_biomass) {
target_MAP_mm
} else {
StandardShrub_MAP_mm
},
map_mm_std_shrubs = StandardShrub_MAP_mm,
map_mm_grasses = if (do_adjust_biomass) {
target_MAP_mm
} else {
StandardGrasses_MAP_mm
},
map_mm_std_grasses = StandardGrasses_MAP_mm,
vegcomp_std_shrubs = StandardShrub_VegComposition,
vegcomp_std_grass = StandardGrasses_VegComposition
)
biom_grasses <-
x[["biom_C3"]] * fgrass_c3c4ann[1] +
x[["biom_C4"]] * fgrass_c3c4ann[2] +
x[["biom_annuals"]] * fgrass_c3c4ann[3]
dimnames(biom_grasses)[[2]] <- btypes2
biom_shrubs <- x[["biom_shrubs"]]
dimnames(biom_shrubs)[[2]] <- btypes2
list(
grass = as.matrix(biom_grasses),
shrub = as.matrix(biom_shrubs)
)
}
#' Lookup transpiration coefficients
#'
#' Lookup transpiration coefficients for grasses, shrubs, and trees per soil
#' layer or per soil depth increment of 1 cm per distribution type for
#' each simulation run and copy values to \var{\sQuote{datafile.soils}}
#'
#' \itemize{
#' \item first row of datafile is label for per soil layer
#' \var{\dQuote{Layer}} or per soil depth increment of 1 cm
#' \var{\dQuote{DepthCM}}
#' \item second row of datafile is source of data
#' \item the other rows contain the data for each distribution type = columns
#' }
#'
#' @param tr_input_code
#' The `"desc"` component of [`rSOILWAT2::sw2_trco_table`].
#' @param tr_input_coeff
#' The `"data"` component of [`rSOILWAT2::sw2_trco_table`].
#' @param soillayer_no An integer value. The number of soil layers.
#' @param trco_type A character string. A column name of `tr_input_code`.
#' @param layers_depth An integer vector. The lower depths of soil layers [`cm`]
#' @param adjustType A character string. The method to adjust prescribed
#' coefficient profile onto provided depth profile `layers_depth`.
#'
#' @seealso [estimate_PotNatVeg_roots()] with example code
#'
#' @export
#' @md
TranspCoeffByVegType <- function(
tr_input_code,
tr_input_coeff,
soillayer_no,
trco_type,
layers_depth,
adjustType = c("positive", "inverse", "allToLast")
) {
#extract data from table by category
trco.code <- as.character(tr_input_code[,
which(colnames(tr_input_code) == trco_type)])
trco <- rep(0, times = soillayer_no)
trco.raw <- na.omit(tr_input_coeff[,
which(colnames(tr_input_coeff) == trco_type)])
if (trco.code == "DepthCM") {
temp <- sum(trco.raw, na.rm = TRUE)
trco_sum <- if (temp == 0 || is.na(temp)) 1L else temp
lup <- 1
for (l in 1:soillayer_no) {
llow <- as.numeric(layers_depth[l])
if (is.na(llow) || lup > length(trco.raw)) {
l <- l - 1
break
}
trco[l] <- sum(trco.raw[lup:llow], na.rm = TRUE) / trco_sum
lup <- llow + 1
}
usel <- l
} else if (trco.code == "Layer") {
usel <- if (length(trco.raw) < soillayer_no) {
length(trco.raw)
} else {
soillayer_no
}
ltemp <- seq_len(usel)
stemp <- sum(trco.raw[ltemp], na.rm = TRUE)
temp <- if (stemp == 0 && is.na(stemp)) 1 else stemp
trco[ltemp] <- trco.raw[ltemp] / temp
}
if (identical(adjustType, "positive")) {
# equivalent to: trco + (1 - sum(trco)) * trco / sum(trco)
trco <- trco / sum(trco)
} else if (identical(adjustType, "inverse")) {
irows <- 1:max(which(trco > 0))
# equivalent to: trco + (1 - sum(trco)) * rev(trco) / sum(trco)
trco[irows] <- trco[irows] + rev(trco[irows]) * (1 / sum(trco[irows]) - 1)
} else if (identical(adjustType, "allToLast")) {
irow <- max(which(trco > 0))
if (irow > 1) {
# adding all the missing roots because soil is too shallow to the
# deepest available layer
trco[irow] <- 1 - sum(trco[1:(irow - 1)])
} else {
trco[1] <- 1
}
}
trco
}
#' Calculate rooting profile for a soil profile from \var{lookup} tables
#'
#' @param layers_depth A numeric vector. Values describe
#' the lower soil layer depths [cm].
#' @param trco_type_by_veg A named list of character strings.
#' The rooting profiles, i.e., column names in the \var{lookup} table,
#' for each vegetation type. A \code{NA} indicates that a rooting profile
#' is not calculated for that vegetation type.
#' The values will be passed to argument \var{trco_type} of function
#' \code{\link{TranspCoeffByVegType}}.
#' @param trco_adj_by_veg A named list of character strings.
#' The type of adjustment from the full rooting profile
#' to the \code{layers_depth}.
#' The values will be passed to argument \var{adjustType} of function
#' \code{\link{TranspCoeffByVegType}}.
#' @param fgrass_c3c4ann A named, numeric vector of length 3.
#' Relative contribution [0-1] of the C3-grasses, C4-grasses, and
#' annuals functional groups. The sum of \code{fgrass_c3c4ann} is 1.
#' @param trco_table A named list with two elements. Default values
#' are taken from \code{\link{sw2_trco_table}}.
#'
#' @seealso \code{\link{TranspCoeffByVegType}}
#'
#' @examples
#' sw_in <- rSOILWAT2::sw_exampleData
#'
#' clim <- calc_SiteClimate(weatherList = rSOILWAT2::weatherData)
#'
#' veg_cover <- rSOILWAT2::estimate_PotNatVeg_composition(
#' MAP_mm = 10 * clim[["MAP_cm"]],
#' MAT_C = clim[["MAT_C"]],
#' mean_monthly_ppt_mm = 10 * clim[["meanMonthlyPPTcm"]],
#' mean_monthly_Temp_C = clim[["meanMonthlyTempC"]],
#' dailyC4vars = clim[["dailyC4vars"]]
#' )
#'
#' estimate_PotNatVeg_roots(
#' layers_depth = c(5, 10, 20, 30, 40, 50, 100, 200),
#' fgrass_c3c4ann = veg_cover[["Grasses"]]
#' )
#'
#' @export
estimate_PotNatVeg_roots <- function(
layers_depth,
trco_type_by_veg = list(
grass_C3 = "SchenkJackson2003_PCdry_grasses",
grass_C4 = "SchenkJackson2003_PCdry_grasses",
grass_annuals = "Jacksonetal1996_crops",
shrub = "SchenkJackson2003_PCdry_shrubs",
forb = "SchenkJackson2003_PCdry_forbs",
tree = "Bradfordetal2014_LodgepolePine"
),
trco_adj_by_veg = list(
grass_C3 = "positive",
grass_C4 = "positive",
grass_annuals = "positive",
shrub = "positive",
forb = "positive",
tree = "positive"
),
fgrass_c3c4ann = c(grass_C3 = NA, grass_C4 = NA, grass_annuals = NA),
trco_table = rSOILWAT2::sw2_trco_table
) {
n_slyrs <- length(layers_depth)
veg_types <- c("Grass", "Shrub", "Tree", "Forb")
res_trco <- array(
data = NA,
dim = c(n_slyrs, length(veg_types)),
dimnames = list(NULL, veg_types)
)
for (k1 in seq_along(veg_types)) {
tmp_type <- sort(grep(
veg_types[k1],
names(trco_type_by_veg),
ignore.case = TRUE,
value = TRUE
))
tmp_adj <- sort(grep(
veg_types[k1],
names(trco_adj_by_veg),
ignore.case = TRUE,
value = TRUE
))
if (!all(tmp_type == tmp_adj)) {
stop(
"Names of `trco_type_by_veg` do not match ",
"`trco_adj_by_veg` for ", shQuote(veg_types[k1])
)
}
if (
isTRUE(!anyNA(unlist(trco_type_by_veg[tmp_type]))) &&
isTRUE(!anyNA(unlist(trco_adj_by_veg[tmp_adj])))
) {
if (length(tmp_type) > 1 && veg_types[k1] == "Grass") {
tmp1 <- sapply(
strsplit(names(fgrass_c3c4ann), split = "_", fixed = TRUE),
FUN = function(x) x[[2]]
)
tmp2 <- sapply(
strsplit(tmp_type, split = "_", fixed = TRUE),
FUN = function(x) x[[2]]
)
tmp_igfcov <- match(tolower(tmp2), tolower(tmp1), nomatch = NA)
if (anyNA(tmp_igfcov)) {
stop(
"Names of `fgrass_c3c4ann` do not match ",
"`trco_adj_by_veg` for ", shQuote(veg_types[k1])
)
}
tmp_root <- rep(0, n_slyrs)
for (k2 in seq_along(tmp_type)) {
if (isTRUE(fgrass_c3c4ann[tmp_igfcov[k2]] > 0)) {
tmp <- TranspCoeffByVegType(
tr_input_code = trco_table[["desc"]],
tr_input_coeff = trco_table[["data"]],
soillayer_no = n_slyrs,
trco_type = trco_type_by_veg[[tmp_type[k2]]],
layers_depth = layers_depth,
adjustType = trco_adj_by_veg[[tmp_adj[k2]]]
)
tmp_root <- tmp_root + fgrass_c3c4ann[tmp_igfcov[k2]] * tmp
}
}
} else if (length(tmp_type) == 1) {
tmp_root <- TranspCoeffByVegType(
tr_input_code = trco_table[["desc"]],
tr_input_coeff = trco_table[["data"]],
soillayer_no = n_slyrs,
trco_type = trco_type_by_veg[[tmp_type]],
layers_depth = layers_depth,
adjustType = trco_adj_by_veg[[tmp_adj]]
)
} else {
stop(
"Root information for ",
shQuote(veg_types[k1]),
" incomplete."
)
}
# Check values
is_good <-
!anyNA(tmp_root) &&
all(tmp_root >= 0) &&
sum(tmp_root) - 1 <= sqrt(.Machine$double.eps)
if (!is_good) {
warning(
"Root information for ",
shQuote(veg_types[k1]),
" is problematic: ",
paste0(round(tmp_root, 4), collapse = " / "),
"; it was re-set to 0s."
)
tmp_root <- rep(0, n_slyrs)
}
res_trco[, veg_types[k1]] <- tmp_root
} else {
warning(
"No rooting profile selected for ",
shQuote(veg_types[k1]), "."
)
}
}
res_trco
}
#' Replace selected biomass values of a
#' \link[rSOILWAT2:swProd-class]{rSOILWAT2::swProd} object
#'
#' @param fg A character string. One of the functional groups represented by
#' \pkg{rSOILWAT2}
#' @param use A named logical vector. The names must represent the column
#' names of a `MonthlyVeg` element of an [`swProd-class`] object
#' @param prod_input A data frame. The values that replace the selected
#' biomass values.
#' @param prod_default A [`swProd-class`] object that contains
#' the `MonthlyVeg` element with biomass values to be updated.
#'
#' @return The requested `MonthlyVeg` element from `prod_default` with updated
#' values.
#'
#' @export
#' @md
update_biomass <- function(
fg = c("Grass", "Shrub", "Tree", "Forb"),
use, # nolint: function_argument_linter.
prod_input, # nolint: function_argument_linter.
prod_default # nolint: function_argument_linter.
) {
fg <- match.arg(fg)
comps <- c("_Litter", "_Biomass", "_FractionLive", "_LAIconv")
veg_ids <- lapply(
comps,
function(x) grep(paste0(fg, x), names(use))
)
veg_incl <- lapply(veg_ids, function(x) use[x])
temp <- swProd_MonProd_veg(prod_default, fg)
if (any(unlist(veg_incl))) {
for (k in seq_along(comps)) if (any(veg_incl[[k]]))
temp[veg_incl[[k]], k] <-
as.numeric(prod_input[, veg_ids[[k]][veg_incl[[k]]]])
}
temp
}
# Determine minimal number of rooted soil layers with veg > 0
get_min_rooted_soil_layers <- function(swInputData) {
veg_comp <- swProd_Composition(swInputData)
soils <- swSoils_Layers(swInputData)
var_veg1 <- c("Grass", "Shrub", "Tree", "Forb")
var_trco <- paste0("transp", var_veg1, "_frac")
var_comp <- sapply(var_veg1, grep, x = names(veg_comp), value = TRUE)
tmp <- apply(
soils[, var_trco, drop = FALSE],
MARGIN = 2,
FUN = function(x) sum(x > 0)
)
min(tmp[veg_comp[var_comp] > 0])
}
#' Check transpiration regions
#'
#' The transpiration regions are checked:
#' \enumerate{
#' \item There is least one transpiration region
#' \item All transpiration regions include at least one soil layer
#' \item Transpiration regions are strictly increasing
#' \item Transpiration regions go no deeper than the most shallow
#' rooting profile of any active vegetation type
#' }
#'
#' @param swInputData A \pkg{rSOILWAT2} input object of class
#' \code{\linkS4class{swInputData}}.
#'
#' @return A logical value.
#'
#' @examples
#' sw_in <- rSOILWAT2::sw_exampleData
#' check_TranspirationRegions(sw_in) ## Expected: TRUE
#'
#' # Make a mistake: set a transpiration region deeper than the rooting profile
#' swSite_TranspirationRegions(sw_in)[2, 2] <- 10
#' check_TranspirationRegions(sw_in) ## Expected: FALSE
#'
#' @export
check_TranspirationRegions <- function(swInputData) {
tr <- swSite_TranspirationRegions(swInputData)
# Checks
nrow(tr) > 0 &&
tr[1, 2] >= 1 &&
tr[nrow(tr), 2] <= get_min_rooted_soil_layers(swInputData) &&
!anyNA(
rSW2utils::check_monotonic_increase(
x = tr,
MARGIN = 2,
strictly = TRUE
)
)
}
#' Adjust transpiration regions for roots and soil layers
#'
#' @param swInputData A \pkg{rSOILWAT2} input object of class
#' \code{\linkS4class{swInputData}}.
#'
#' @return A transpiration region matrix.
#'
#' @examples
#' sw_in <- rSOILWAT2::sw_exampleData
#'
#' adjust_TranspirationRegions(sw_in)
#'
#' @export
adjust_TranspirationRegions <- function(swInputData) {
tr <- swSite_TranspirationRegions(swInputData)
n_tr <- min(4, nrow(tr))
tri_file <- cbind(Used_TF = 1, DeepestLayer = tr[, 2])
# adjust maximum transpiration region for minimum soil depth and rooting depth
n_max <- min(
nrow(swSoils_Layers(swInputData)),
get_min_rooted_soil_layers(swInputData)
)
if (max(tri_file[tri_file[, 1] > 0, 2], na.rm = TRUE) > n_max) {
for (k in rev(seq_len(n_tr))) {
if (tri_file[k, 1] > 0) {
if (tri_file[k, 2] > n_max) {
tri_file[k, 2] <- tr[k, 2] <- n_max
}
if (k > 1 && tri_file[k, 2] <= tri_file[k - 1, 2]) {
tr <- matrix(tr[-k, ], ncol = 2)
}
}
}
}
tr
}
#' Prepare transpiration regions
#'
#' Translate a soil vector of transpiration region values into a
#' transpiration region matrix.
#'
#' @param tr_lyrs An integer vector.
#' The transpiration region for each soil layer.
#'
#' @return A transpiration region matrix.
#'
#' @seealso
#' \code{\link{adjust_TranspirationRegions}} and
#' \code{\link{check_TranspirationRegions}}
#'
#' @examples
#' # Example values correspond to `CONUSSOIL_BSE_EVERY10cm` of \pkg{rSFSW2}:
#' prepare_TranspirationRegions(tr_lyrs = c(1, 1, 1, 2, 2, 2, 2, 3, 3, 3))
#'
#' @export
prepare_TranspirationRegions <- function(tr_lyrs) {
tmp <- cumsum(rle(tr_lyrs)[["lengths"]])
cbind(ndx = seq_along(tmp), layer = tmp)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.