#--- Helper functions ----------------------------------------------------------
check_all_output_available_of_run <- (
path_to_run,
N_scen,
zipped_runs = ,
check_input =
) {
if (zipped_runs) {
# data frame with file names and file sizes
<- ::(path_to_run, = )
fnum <- ()
fsizes <- , "Length", = ]
} {
<- ((path_to_run), full.names = )
fnum <- ()
fsizes <- (
,
(f) (f),
FUN.VALUE = ,
USE.NAMES =
)
}
# Expect: One output file per scenario plus one file for inputs
((
fnum == N_scen + if (check_input) 1 0,
(fsizes > 0)
))
}
#' Shorten full names of \pkg{rSFSW2} simulation runs
#'
#' Full names of \pkg{rSFSW2} simulation runs are composed of three parts
#' \itemize{
#' \item a run number
#' \item experimental treatment
#' \item site label
#' }
#' and concatenated by \var{"_"} as separator, for instance,
#' \var{"001_Default_SiteA"}, \var{"002_Default_SiteB"}, etc.
#'
#' It may be desirable to shorten the full names in cases where only none/one
#' experimental treatment was used, for instance,
#' \var{SiteA"}, \var{"SiteB"}, etc.
#'
#' @param run_names A vector of character strings.
#' @param element_sep A character value. The separating character.
#' @param N_discard An integer value. The number of elements separated by
#' \code{element_sep} to discard.
#'
#' @examples
#' shorten_run_names(c("001_Default_SiteA", "002_Default_SiteB"))
#'
#' @export
<- (run_names, element_sep = "_", N_discard = 2) {
discard <- (N_discard)
(
(run_names, = element_sep, fixed = ),
(x) (x- discard], collapse = element_sep),
FUN.VALUE =
)
}
prepare_soils_for_site <- (
path,
name_sw2_run,
name_sw2_run_soils = ,
zipped_runs = ,
soils = ,
soil_variables = ,
var_soilsite = "site"
) {
used_soil <-
if (!(soils) && !(soils) && !(name_sw2_run_soils)) {
#--- Soils are pre-extracted: subset `soils`
# Check that we got good soil variable names
names_soil_variables <- (soil_variables)
(
names_soil_variables %in% (soils),
names_soil_variables %in% (())
)
# Locate site
idss <- (
names_soil_variables,
(k) (name_sw2_run_soils, soils[k]], var_soilsite])
)
(idss) <- names_soil_variables
(
(idss, , FUN.VALUE = ),
(idss, , FUN.VALUE = ) == 1L
)
# Prepare soil variable values
used_soil <- (
names_soil_variables,
(k) {
ics <- ("_L[[:digit:]]+$", (soils[k]]))
(soils[k]][idss[k]], ics])
}
)
(used_soil) <- names_soil_variables
} {
#--- Soils are not pre-extracted: read values from files
if ((soil_variables)) {
soil_variables <- ()
}
nsv <- (soil_variables)
tmp_soils <- get_soillayers_variable(
path = path,
name_sw2_run = name_sw2_run,
id_scen = 1,
zipped_runs = zipped_runs,
sw2_soil_var = nsv
)
used_soil <- (
depth_cm = if ("depth_cm" %in% nsv) tmp_soils"depth_cm", ],
sand_frac = if ("sand_frac" %in% nsv) tmp_soils"sand_frac", ],
clay_frac = if ("clay_frac" %in% nsv) tmp_soils"clay_frac", ],
gravel_content = if ("gravel_content" %in% nsv) {
tmp_soils"gravel_content", ]
}
)
}
used_soil
}
#--- Obtain input values -------------------------------------------------------
#' Determine widths (as weights) of soil layers within a zone
#'
#' @param soil_depths_cm A numeric vector.
#' The lower depth limits of soil layers in \var{[cm]}.
#' @param used_depth_range_cm A numeric vector of length two.
#' The upper and lower depth limit of the zone (depth range) to consider.
#' @param n_slyrs_has An integer value. The number of simulated soil layers
#' (optional). The code throws an error if there are fewer soil layers
#' than selected by \code{used_depth_range_cm} from \code{soil_depths_cm}.
#'
#' @return A numeric vector of the same length as \code{soil_depths_cm}.
#' The values correspond to soil layer widths that fall within the
#' selected zone. Other values, outside the zone, are set to \code{NA}.
#'
#' @examples
#' calc_soillayer_weights(c(5, 10, 30, 50), NULL)
#' calc_soillayer_weights(c(5, 10, 30, 50), c(0, 100))
#' calc_soillayer_weights(c(5, 10, 30, 50), c(0, 40))
#' calc_soillayer_weights(c(5, 10, 30, 50), c(20, 40))
#' calc_soillayer_weights(c(5, 10, 30, 50), c(20, 50))
#' calc_soillayer_weights(c(5, 10, 30, 50), c(20, 30), 3)
#' \dontrun{
#' calc_soillayer_weights(c(5, 10, 30, 50), c(20, 50), 3)
#' }
#'
#' @export
<- (
soil_depths_cm,
used_depth_range_cm = ,
n_slyrs_has =
) {
# if used_depth_range_cm is NULL, then ids is `logical(0)`
ids <-
soil_depths_cm <= used_depth_range_cm[1]] |
soil_depths_cm > used_depth_range_cm[2]]
x <- ((0, soil_depths_cm))
x[ids] <-
if (!(n_slyrs_has)) {
if (((tmp <- ((!(x)))) > n_slyrs_has)) {
(
"Deeper soil layers requested than actually simulated:",
"\n * position of deepest requested soil layers = ", tmp,
"\n * simulated number of soil layers = ", n_slyrs_has
)
}
# Make sure that returned object is not longer than `n_slyrs_has`
if (n_slyrs_has < (x)) {
x <- x[seq_len(n_slyrs_has)]
}
}
x
}
#' Check whether soil layers are available in a specified zone
#'
#' @inheritParams calc_soillayer_weights
#' @param strict A logical vector either of length one or equal to the
#' length of \code{used_depth_range_cm}. \code{TRUE} indicates that
#' the corresponding value (or all) of \code{used_depth_range_cm} is required
#' to be either zero or one of the values of \code{soil_depths_cm}.
#' @param type A character string. If a value is not present, then either
#' a warning is issued or an error.
#'
#' @return (Invisibly) a logical value indicating whether at least one
#' required value is missing (\code{FALSE}) or all are present (\code{TRUE}).
#'
#' @examples
#' check_soillayer_availability(
#' c(5, 10, 30, 50), c(0, 40),
#' strict = TRUE,
#' type = "warn"
#' )
#' check_soillayer_availability(
#' c(5, 10, 30, 50), c(0, 40),
#' strict = c(TRUE, FALSE),
#' type = "warn"
#' )
#' check_soillayer_availability(
#' c(5, 10, 30, 50), c(20, 30, 40),
#' strict = TRUE,
#' type = "warn"
#' )
#'
#' @export
<- (
soil_depths_cm,
used_depth_range_cm = ,
strict = ,
= ("warn", "error")
) {
sim_soil_depths_cm <- (((0, soil_depths_cm)))
msg <- (0)
if ((strict) && !(used_depth_range_cm)) {
used_depth_range_cm <- ((used_depth_range_cm))
strict <- (strict, (used_depth_range_cm))
ids <- used_depth_range_cm[strict] %in% sim_soil_depths_cm
if (!(ids)) {
msg <- (
"Requested soil depth(s) at",
(used_depth_range_cm[strict]!ids]),
"are not available in the simulated set of soil layers of",
(sim_soil_depths_cm)
)
if (() == "error") (msg) (msg)
}
}
((msg) == 0)
}
#' Determine positions of soil layers within a zone
#'
#' @inheritParams calc_soillayer_weights
#'
#' @return An integer vector.
#'
#' @examples
#' determine_used_soillayers(c(5, 10, 30, 50))
#' determine_used_soillayers(c(5, 10, 30, 50), c(0, 100))
#' determine_used_soillayers(c(5, 10, 30, 50), c(0, 40))
#' determine_used_soillayers(c(5, 10, 30, 50), c(20, 40))
#' determine_used_soillayers(c(5, 10, 30, 50), c(20, 50))
#' determine_used_soillayers(c(5, 10, 30, 50), c(20, 30), 3)
#' \dontrun{
#' determine_used_soillayers(c(5, 10, 30, 50), c(20, 50), 3)
#' }
#'
#' @export
<- (
soil_depths_cm,
used_depth_range_cm = ,
n_slyrs_has =
) {
soil_depths_cm <- soil_depths_cm!(soil_depths_cm)]
x <- (soil_depths_cm)
# if used_depth_range_cm is NULL, then ids is `logical(0)`
ids <-
soil_depths_cm <= used_depth_range_cm[1]] |
soil_depths_cm > used_depth_range_cm[2]]
if ((ids)) {
x <- x[which(!ids)]
}
if (
!(n_slyrs_has) &&
(x) > 0 &&
((tmp <- (x)) > n_slyrs_has)
) {
(
"Deeper soil layers requested than actually simulated:",
"\n * position of deepest requested soil layers = ", tmp,
"\n * simulated number of soil layers = ", n_slyrs_has
)
}
x
}
#------ SEUG DEVELOPMENT ERA ------
groupid_by_days <- (
start_year,
end_year,
group_by_month,
first_month_of_year = 1L
) {
days <- (
from = (start_year, 1, 1, tz = "UTC"),
to = (end_year, 12, 31, tz = "UTC"),
= "1 day"
)
<- (days)$mon + 1 # nolint: extraction_operator_linter.
rlels <- ()"lengths"]]
(
groupid = ((group_by_month, (rlels)), times = rlels),
ids_adj_yrs = if (first_month_of_year > 1) {
>= first_month_of_year
}
)
}
weighted_mean_across_soillayers <- (x, w) {
ids_lyrs <- ((x), (!(w)))
(x[ids_lyrs], w[ids_lyrs])
}
#' Extract values, determine time series, and assign time steps to groups
#' based on an \pkg{rSOILWAT2} output as generated by \pkg{rSFSW2}
#'
#' @param id_scen An integer value. The climate scenario ID as generated by
#' \pkg{rSFSW2}.
#' @param path A character string. The directory path where \pkg{rSFSW2}
#' stored the \pkg{rSOILWAT2} output, by default \dQuote{3_Runs}.
#' @param name_sw2_run A character string. The \pkg{rSFSW2} name identifying
#' a run, i.e., the full experimental label x site name. The folder name
#' within \code{path}.
#' @param zipped_runs A logical value. Describes whether \code{name_sw2_run}
#' is a zip archive or a regular folder.
#' @param group_by_month An integer vector of length 12. The group IDs of each
#' month. If no grouping is requested, then provide,
#' e.g., \code{rep(0, 12)}.
#' @param first_month_of_year An integer value. The number of the
#' month `[1-12]` when the "year" should start,
#' e.g., `10` for a "water-year" that starts in October.
#' Note: the first and last year are incomplete and will be discarded
#' if `first_month_of_year` is not `1`.
#' @param sw2_tp A character string. One of the \pkg{rSOILWAT2} output time
#' slots. Note: currently, only \dQuote{Month} and \dQuote{Day} are
#' implemented.
#' @param sw2_out A character vector. One (or multiple) of the \pkg{rSOILWAT2}
#' output groups. These will be recycled if \code{sw_var} is longer.
#' @param sw2_var A character vector. One (or multiple) of the column names
#' of the output group(s) corresponding to the \code{sw2_out}.
#' @param varnames_are_fixed A logical value. If \code{TRUE}, then
#' \code{sw2_var} is taken literally. Otherwise, \code{sw2_var}
#' are interpreted as regular expressions
#' (and may thus return multiple columns).
#'
#' @return A list with three elements: \itemize{
#' \item \var{vals} A list with one element for each \code{sw2_var}. Each
#' element is a vector of values for that \code{sw2_var},
#' if \code{varnames_are_fixed} is \code{TRUE}; otherwise, a vector
#' or matrix where columns represent all the \pkg{SOILWAT2} output
#' corresponding to that \code{sw2_var} interpreted as
#' regular expression.
#' The length of vectors (or the number of rows) corresponds to the
#' number time steps requested by \code{sw2_tp}.
#' \item \var{year} A vector of calendar years for each time step.
#' \item \var{groups_by_time} A vector of group ID, based on
#' \code{group_by_month}, for each time step.
#' }
#'
#' @md
<- (id_scen, path, name_sw2_run,
group_by_month, first_month_of_year,
sw2_tp, sw2_out, sw2_var, varnames_are_fixed = ,
zipped_runs =
) {
x <- (
path = path,
name_sw2_run = name_sw2_run,
zipped_runs = zipped_runs,
id_scen = id_scen,
sw2_tp = sw2_tp,
sw2_outs = sw2_out,
sw2_vars = sw2_var,
varnames_are_fixed = varnames_are_fixed
)
# Grouping by seasons and
# year-adjustment for incomplete first/last occurrence of first group
x_years <- x"time"]], "Year"]
if (sw2_tp == "Month") {
x_groups <- (group_by_month, (x_years))
if (first_month_of_year > 1) {
ids <- x"time"]], "Month"] >= first_month_of_year
x_years[ids] <- x_years[ids] + 1
}
} if (sw2_tp == "Day") {
ids <- groupid_by_days(
start_year = x_years[1]],
end_year = x_years[length(x_years)],
group_by_month,
first_month_of_year
)
x_groups <- ids"groupid"]]
if (first_month_of_year > 1) {
x_years[ids"ids_adj_yrs"]]] <- x_years[ids"ids_adj_yrs"]]] + 1
}
}
(
vals = x"values"]],
years = x_years,
groups_by_time = x_groups
)
}
#------ NEWRR DEVELOPMENT ERA ------
get_swp_weighted <- (
path, name_sw2_run, id_scen,
years,
soils,
zipped_runs = ,
used_depth_range_cm = ,
) {
("`get_swp_weighted()` uses matric-VWC!")
vwc <- (
path = path,
name_sw2_run = name_sw2_run,
zipped_runs = zipped_runs,
id_scen = id_scen,
years = years,
sw2_tp = "Day",
sw2_outs = "VWCMATRIC",
sw2_vars = "Lyr",
varnames_are_fixed =
)
N_days <- (vwc"values"]][1]])
T_by_lyr <- (
path = path,
name_sw2_run = name_sw2_run,
zipped_runs = zipped_runs,
id_scen = id_scen,
years = years,
sw2_tp = "Day",
sw2_outs = "TRANSP",
sw2_vars = "transp_total_Lyr",
varnames_are_fixed =
)
widths_cm <- (
soil_depths_cm = soils"depth_cm"]],
used_depth_range_cm = used_depth_range_cm
)
id_slyrs <- (!(widths_cm))
widths_cm <- widths_cm[id_slyrs]
tmp <- (, N_days)
for (k (N_days)) {
tmp[k] <- if ((T_by_lyr"values"]][1]][k, id_slyrs]) > 0) {
# values weighted by transpiration per layer
rSOILWAT2::VWCtoSWP(
vwc = (
x = vwc"values"]][1]][k, id_slyrs],
w = T_by_lyr"values"]][1]][k, id_slyrs]
),
sand = (
x = soils"sand_frac"]][id_slyrs],
w = T_by_lyr"values"]][1]][k, id_slyrs]
),
clay = (
x = soils"clay_frac"]][id_slyrs],
w = T_by_lyr"values"]][1]][k, id_slyrs]
)
)
} {
# values weighted by layer width
rSOILWAT2::VWCtoSWP(
vwc = (
x = vwc"values"]][1]][k, id_slyrs],
w = widths_cm
),
sand = (
x = soils"sand_frac"]][id_slyrs],
w = widths_cm
),
clay = (
x = soils"clay_frac"]][id_slyrs],
w = widths_cm
)
)
}
}
(
= vwc"time"]],
values = (swp_weighted = tmp)
)
}
calc_new_yearly_aggregations <- (
x_daily,
temp_monthly,
fun_time, # sum
fun_extreme, # max
periods, # for instance: list(op = `>`, limit = 0),
output = (
"values", "5roll", "seasonal_variability", "seasonality",
"mean_day", "extreme_mean005day", "extreme_mean010day",
"extreme_value_consecutive_periods",
"mean_duration_consecutive_periods_days",
"extreme_duration_consecutive_periods_days"
),
include_year =
) {
# one value per year
# - mean/sum per year
# - mean/sum per month: CV across months, cor(., temperature)
# - moving mean across 5- and 10-day windows: extreme across smoothed days
output <- (output, several.ok = )
((x_daily"values"]]) == 1)
years <- (x_daily"time"]], "Year"])
res <- (
= (years),
= (output) + (include_year),
= (, if (include_year) ("Year", output) output)
)
if (include_year) {
res, "Year"] <- (x_daily"time"]], "Year"])
}
if ((("values", "5roll") %in% output)) {
(!(fun_time), !(fun_time))
x_yearly <- (
X = x_daily"values"]][1]],
INDEX = x_daily"time"]], "Year"],
FUN = fun_time
)
if ("values" %in% output) {
res, "values"] <- x_yearly
}
if ("5roll" %in% output && ("zoo", quietly = )) {
res, "5roll"] <- zoo::rollapply(
= x_yearly,
= 5,
FUN = ,
fill = ,
partial = ,
align = "center"
)
}
}
if ((("seasonal_variability", "seasonality") %in% output)) {
(!(fun_time), !(fun_time))
x_monthly <- (
x = x_daily"values"]][1]],
= (
Month = x_daily"time"]], "Month"],
Year = x_daily"time"]], "Year"]
),
FUN = fun_time
)
if ("seasonal_variability" %in% output) {
res, "seasonal_variability"] <- (
X = x_monthly"x"]],
INDEX = x_monthly"Year"]],
FUN = (x) (x) / (x)
)
}
if ("seasonality" %in% output) {
(
!(temp_monthly),
!(temp_monthly),
x_monthly"Year"]] == temp_monthly"time"]], "Year"]
)
res, "seasonality"] <- ((
= (x_monthly"x"]], temp_monthly"values"]][1]]),
INDICES = x_monthly"Year"]],
FUN = (x) (x, 1], x, 2])
))
}
}
if ("mean_day" %in% output) {
res, "mean_day"] <- (
X = x_daily"values"]][1]],
INDEX = x_daily"time"]], "Year"],
FUN =
)
}
if (
"extreme_mean005day" %in% output &&
("zoo", quietly = )
) {
(!(fun_extreme), !(fun_extreme))
res, "extreme_mean005day"] <- (
X = zoo::rollapply(
= x_daily"values"]][1]],
= 5,
FUN = ,
fill = ,
partial = ,
align = "center"
),
INDEX = x_daily"time"]], "Year"],
FUN = fun_extreme
)
}
if (
"extreme_mean010day" %in% output &&
("zoo", quietly = )
) {
(!(fun_extreme), !(fun_extreme))
res, "extreme_mean010day"] <- (
X = zoo::rollapply(
= x_daily"values"]][1]],
= 10,
FUN = ,
fill = ,
partial = ,
align = "center"
),
INDEX = x_daily"time"]], "Year"],
FUN = fun_extreme
)
}
var_durations <- (
"mean_duration_consecutive_periods_days",
"extreme_duration_consecutive_periods_days"
)
var_periods <- (
"extreme_value_consecutive_periods",
var_durations
)
if ((var_periods %in% output)) {
(!(periods), !(periods))
x_periods <- (
what = periods"op"]],
= (
x_daily"values"]][1]],
periods"limit"]]
)
)
if ("extreme_value_consecutive_periods" %in% output) {
(
!(fun_time), !(fun_time),
!(fun_extreme), !(fun_extreme)
)
res, "extreme_value_consecutive_periods"] <- (
X = (
= (x_periods, x_daily"values"]][1]]),
INDICES = x_daily"time"]], "Year"],
FUN = (x) {
tmp <- (x, 1])
if ((tmp"values"]])) {
# Propagate NAs
} if ((tmp"values"]] == 1, na.rm = )) {
# Create index of days for each spell when
# `x_periods` is TRUE (i.e., `tmp[["values"]] == 1`)
ips <- (tmp"lengths"]])
ids <- (ips, tmp"lengths"]])
ids[ids %in% ips[tmp"values"]] != 1]] <-
(X = x, 2], INDEX = ids, FUN = fun_time)
} {
# No days when `x_periods` is TRUE
0
}
}
),
FUN = fun_extreme,
FUN.VALUE =
)
}
if ((var_periods %in% output)) {
x_consecutive_periods_days <- (
X = x_periods,
INDEX = x_daily"time"]], "Year"],
FUN = (x) {
tmp <- (x)
if ((tmp"values"]])) {
# Propagate NAs
} if ((tmp"values"]] == 1, na.rm = )) {
# Select duration of spells when
# `x_periods` is TRUE (i.e., `tmp[["values"]] == 1`)
tmp"lengths"]][tmp"values"]]]
} {
# No days when `x_periods` is TRUE
0
}
}
)
if ("mean_duration_consecutive_periods_days" %in% output) {
res, "mean_duration_consecutive_periods_days"] <- (
X = x_consecutive_periods_days,
FUN = ,
FUN.VALUE =
)
}
if ("extreme_duration_consecutive_periods_days" %in% output) {
(!(fun_extreme), !(fun_extreme))
res, "extreme_duration_consecutive_periods_days"] <- (
X = x_consecutive_periods_days,
FUN = fun_extreme,
FUN.VALUE =
)
}
}
}
res
}
#--- DOY of 10%, 50%, and 90% percentile of cumulative transpiration
#--- Annual total transpiration (mm)
calc_transp_seasonality <- (x, , probs) {
(
x = x,
= (Year = ),
FUN = (x) {
tdc <- (x)
ttot <- tdc[length(tdc)]
(
ttot,
if ((ttot)) {
(, (probs))
} {
(
ttot * probs,
(v) (v >= tdc),
FUN.VALUE =
)
}
)
}
)
}
#--- DOY of two maxima of a smoothed transpiration
#' @examples
#' ## No or less than two plateaus -> central candidate from among all values
#' central_candidate(1)
#' central_candidate(1:4)
#' central_candidate(c(1, 20:29))
#'
#' ## More than one plateau -> central candidate from longest plateau
#' central_candidate(c(1, 20:30, 40:60))
#' central_candidate(c(1, 20:30, 40:60, 70:80, 100:110))
#' central_candidate(c(1, 20:30, 40:60, 70:90))
#'
#' @noRd
central_candidate <- (ids) {
# TODO: consider re-writing in cpp11 or Rcpp
if ((ids) > 1) {
# plateaus or individual peaks?
tmp <- ids-1] - ids-(ids)] == 1 # faster than `diff(ids) == 1`
if ((tmp)) {
# We have plateaus: now, count number of plateaus
tmpp <- ((tmp, ))
tmpps <- (tmpp"values"]])
if ((tmpps) > 1) {
# We have multiple plateaus: now, identify longest plateau(s)
tmppl <- tmpp"lengths"]][tmpps]
tmpi <- (tmppl == (tmppl))
tmppsi <- if ((tmpi) > 1) {
# We have multiple longest plateaus:
# now, select nearest even middle one among longest plateaus
tmpps[tmpi[round(0.5 + (tmpi) / 2)]]
} {
tmpps[tmpi]
}
# Extract start:end of identified longest plateau
tmpcs <- ((1, tmpp"lengths"]]))
ids <- ids[tmpcs[tmppsi]:tmpcs[tmppsi + 1]]
}
}
# take central value (from among all peaks or from selected plateau)
# here, the 50%-quantile of the positions under type 1
ids((ids) - 1) %/% 2 + 1]
} {
ids
}
}
identify_simple_peaks <- (
x,
,
fun = (x) (x) == ( + 1) / 2
) {
(("zoo", quietly = ))
(zoo::rollapply(
x,
= ,
FUN = fun,
align = "center",
fill =
))
}
#' @section Notes: Peaks closer than a half-\code{window} width from either
#' end of \code{x} will not be identified. Add an appropriate
#' buffer around \code{x} if such peaks are important.
identify_peaks <- (x, , = ("maxima", "minima")) {
( %% 2L == 1L) # check that window is odd
<- ()
if ( == "maxima") {
fpeak <- (x) {
tmpx <- (x, na.rm = )
if ((tmpx) && (tmpx > (x, na.rm = ))) {
central_candidate((x >= tmpx)) == ( + 1) / 2
} {
}
}
} {
fpeak <- (x) {
tmpx <- (x, na.rm = )
if ((tmpx) && (tmpx < (x, na.rm = ))) {
central_candidate((x <= tmpx)) == ( + 1) / 2
} {
}
}
}
res <- identify_simple_peaks(x, , fun = fpeak)
# We may get runs of peaks if a peak is an almost window-wide plateau
tmp <- (res) == 1
if ((tmp)) {
# We have a run of peaks: select central one
tmpp <- ((tmp, ))
tmp2 <- ((0, tmpp"lengths"]])) + 1
tmps <- tmp2-(tmp2)][tmpp"values"]]]
tmp <- (
= tmps,
= tmps + tmpp"lengths"]][tmpp"values"]]]
)
for (k1 (tmps)) {
ids <- tmp"start"]][k1]:tmp"end"]][k1]
res[ids]-((ids) / 2)] <-
}
res <- ((res))
}
res
}
# Identify valley bottoms between peaks
identify_valleys_between_peaks <- (x, peaks) {
ids_edges <- ((!(x)))
# Peaks located at start/end of non-NA value period
ids_peak_on_limits <- (peaks %in% ids_edges)
# Number of valleys
# (one on either side of a peak unless peak is on start/end of data)
n_valleys <- (peaks) + 1 - (ids_peak_on_limits)
valleys <- (, n_valleys)
iv <- 1
for (k (peaks)) {
# Skip to next if current peak on an edge and more than one peak
if (k %in% ids_peak_on_limits && n_valleys > 1)
#--- Identify valley before peak
if (!(peaks[k] %in% ids_edges[1]]) && (valleys[max(1, iv - 1)])) {
ids_inbetween <- (
from = if (k > 1) peaks[k - 1] 1,
to = peaks[k]
)
tmpx <- (x[ids_inbetween], na.rm = )
if ((tmpx) && (tmpx < (x, na.rm = ))) {
valleys[iv] <-
ids_inbetween[1]] - 1 +
central_candidate((x[ids_inbetween] <= tmpx))
iv <- iv + 1
}
}
#--- Identify valley after peak
if (!(peaks[k] %in% ids_edges[2]])) {
ids_inbetween <- (
from = peaks[k],
to = if (k == (peaks)) (x) peaks[k + 1]
)
tmpx <- (x[ids_inbetween], na.rm = )
if ((tmpx) && (tmpx < (x, na.rm = ))) {
valleys[iv] <-
ids_inbetween[1]] - 1 +
central_candidate((x[ids_inbetween] <= tmpx))
iv <- iv + 1
}
}
}
valleys
}
peak_size_v1 <- (pids, peak_type = ("value", "volume"), values) {
peak_type <- (peak_type)
if (peak_type == "value") {
# smoothed transpiration at time of peak
pvals <- values[pids]
} {
# summed smoothed transpiration around peaks
# Identify valley bottoms between peaks
valleys <- identify_valleys_between_peaks(x = values, peaks = pids)
vids <- ((1, valleys, (values)))
# Sum smoothed transpiration from left to right valley bottoms per peak
pvals <- (
pids,
(p) {
tmp <- (p, vids)
(values[vids[tmp]:vids[tmp + 1]], na.rm = )
},
FUN.VALUE =
)
}
pvals
}
peak_size_v2 <- (
pids, peak_type = ("value", "volume"), , values
) {
peak_type <- (peak_type)
if (peak_type == "value") {
# smoothed transpiration at time of peak
pvals <- values[pids]
} {
# summed smoothed transpiration around peaks
# Identify valley bottoms between peaks
valleys1 <- identify_valleys_between_peaks(x = values, peaks = pids)
# Identify negated peaks
valleys2 <- identify_peaks(x = values, = , = "minima")
vids <- (((1, valleys1, valleys2, (values))))
# Sum smoothed transpiration from left to right valley bottoms per peak
pvals <- (
pids,
(p) {
tmp <- (p, vids)
(values[vids[tmp]:vids[tmp + 1]], na.rm = )
},
FUN.VALUE =
)
}
pvals
}
#------ ATLINKAGES DEVELOPMENT ERA ------
get_variable_in_months <- (
path, name_sw2_run, id_scen_used,
list_years_scen_used,
sw2_out, sw2_var, var_label,
= 1:12,
fun_time = ,
var_scaler = 1,
include_year = ,
zipped_runs = ,
) {
res <- ()
for (k1 (id_scen_used)) {
tmp_mon <- collect_sw2_sim_data(
path = path,
name_sw2_run = name_sw2_run,
id_scen = id_scen_used[k1],
years = list_years_scen_used[k1]],
output_sets = (
(
sw2_tp = "Month",
sw2_outs = sw2_out,
sw2_vars = sw2_var,
varnames_are_fixed =
)
),
zipped_runs = zipped_runs
)
# Calculate and format
ts_years <- (tmp_mon[1]]"time"]], "Year"])
res[k1]] <- format_values_to_matrix(
x = var_scaler * ((
X = tmp_mon[1]]"values"]][sw2_var]],
INDEX = tmp_mon[1]]"time"]], "Year"],
FUN = (x) fun_time(x[months])
)),
ts_years = ts_years,
timestep = "yearly",
out_label = var_label
)
if (include_year) {
res[k1]] <- (
Year = ts_years,
res[k1]]
)
}
}
res
}
#------ FORMAT VALUES ------
#' Convert daily, monthly, or yearly values to a matrix
#'
#' @param x A numeric vector (or a list with a vector), data.frame, or a matrix.
#' Rows represent time and columns separate variables.
#' @param ts_years A numeric vector.
#' The time series of years that the rows of `x` represent.
#' If a climatology, then set `ts_years` to `NA`.
#' @param out_label A character vector. The labels representing columns of `x`.
#' @param timestep A character string.
#' @param include_year A logical value. If `TRUE` then adds a first row of
#' years.
#'
#' @return A matrix where columns represent years and
#' rows represent within-year values
#' (repeated blocks if multiple elements/columns in `x`), i.e.,
#' none for yearly,
#' month of year 1 through 12 for monthly,
#' or day of year 1 through 366 for daily time step.
#' Row names are constructed from the combinations of time step,
#' the names of `x` (if any), and `out_label` (if not `NULL`).
#'
#' @examples
#' format_values_to_matrix(
#' data.frame(a = 1:10, b = 20:29),
#' ts_years = 1:10,
#' timestep = "yearly",
#' out_label = NULL
#' )
#' format_values_to_matrix(
#' data.frame(a = 1:10, b = 20:29),
#' ts_years = NA,
#' timestep = "yearly",
#' out_label = LETTERS[1:10]
#' )
#'
#' @noRd
format_values_to_matrix <- (
x,
ts_years,
timestep = ("daily", "monthly", "yearly"),
out_label = ,
include_year =
) {
timestep <- (timestep)
years <- (ts_years)
is_clim <- ((years))
if (!(x)) {
if ((x)) {
ns_x <- (x)
x <- (x)
(x) <- ns_x
} {
x <- (x)
}
}
if (timestep == "daily") {
doys366 <- (366)
doys365 <- (365)
x_template <- ( = (366, (years)))
}
#--- Determine row names
tmp <- (
(
EXPR = timestep,
yearly = "annual",
monthly = ("mon", ((12), = 2, flag = "0")),
daily = ("doy", (doys366, = 3, flag = "0"))
),
(x),
out_label
)
ns_row <- (
expand.grid(tmp[lengths(tmp) > 0]),
MARGIN = 1,
FUN = (x) ((x), collapse = "_")
)
#--- Format as matrix (columns represent years)
if (timestep == "yearly") {
tmp <- (x)
res <- (
= tmp,
= if (is_clim) {
(tmp)
} if (((x) == (ts_years))) {
(x)
} {
1
},
= (ts_years),
= (ns_row, ),
byrow =
)
} {
tmp <- if (timestep == "monthly") {
(
(x),
(k) {
(
= x[k]],
= (12, (years))
)
}
)
} if (timestep == "daily") {
(
(x),
(k1) {
xtmp <- x_template
for (k2 (years)) {
if (is_clim) {
xtmp[seq_along(x[k1]]), k2] <- x[k1]]
} {
doys <- if (rSW2utils::isLeapYear(years[k2])) doys366 doys365
xtmp[doys, k2] <- x[k1]][ts_years == years[k2]]
}
}
xtmp
}
)
}
res <- (, tmp)
(res) <- ns_row
}
if (include_year) {
res <- (Year = years, res)
}
res
}
#------ SEPARATE DATA FROM CALCULATIONS DEVELOPMENT ERA ------
create_sw2simtime <- (n) {
(
Year = (, n),
Month = ,
Day = ,
=
)
}
#' Prepare specified vs. requested time step `data.frame`
#'
#' @param xt A two-dimensional object with columns
#' `Year`,
#' `Month` (if `sw2_tp` equals "Month"), and
#' `Day` (if `sw2_tp` equals "Day"). Each row represents one time step
#' according to `sw2_tp`.
#' @param req_years An integer vector of Calendar years. If missing or `NULL`,
#' then time steps contained in `xt` are used.
#' If not missing and not `NULL` and different than years contained in `xt`,
#' then time steps combined from requested and simulated are used.
#' @param sw2_tp A character string. The daily, monthly, or yearly time step
#' describing content of `xt` and determining the same output time step.
#'
#' @return A `data.frame` with four columns `Year`, `Month`, `Day`, `mode`.
#' Rows represent daily, monthly, or yearly time steps that combine
#' requested (if `req_years` was provided) and simulated (`xt`) time steps.
#' `mode` is `"nosim"` if a time step was not simulated (but requested),
#' `"sim_keep"` if a time step was simulated
#' (and requested or `req_years` is missing), or
#' `"sim_discard"` if a time step was simulated (but not requested).
#'
#' @section Notes:
#' The column `Day` contains day of year (and not day of month)!
#'
#' @examples
#' xty <- data.frame(Year = 1991:2020)
#' determine_sw2_sim_time(xty, sw2_tp = "Year")
#' determine_sw2_sim_time(xty, req_years = 1981:2010, sw2_tp = "Year")
#' determine_sw2_sim_time(xty, req_years = 2011:2030, sw2_tp = "Year")
#'
#' @export
<- (
xt,
req_years = ,
sw2_tp = ("Day", "Month", "Year")
) {
# nolint start: commented_code_linter, line_length_linter.
# debugging un-memoized version:
# assignInNamespace("determine_sw2_sim_time", environment(rSW2metrics::determine_sw2_sim_time)$`_f`, "rSW2metrics")
# nolint end
sw2_tp <- (sw2_tp)
years_sim <- (xt, "Year"])
has_req_yrs <-
!(req_years) && (req_years) > 0 &&
!(req_years, years_sim)
years_used <- if (has_req_yrs) {
(((req_years, years_sim)))
} {
years_sim
}
if (sw2_tp == "Year") {
x_time <- create_sw2simtime(n = (years_used))
x_time, "Year"] <- years_used
if (has_req_yrs) {
ids_sim <- x_time, "Year"] %in% years_sim
}
} if (sw2_tp == "Month") {
x_time <- create_sw2simtime(n = 12 * (years_used))
x_time, "Year"] <- (years_used, each = 12)
x_time, "Month"] <- (12)
if (has_req_yrs) {
ids_sim <-
(x_time, "Year"], "-", x_time, "Month"]) %in%
(xt, "Year"], "-", xt, "Month"])
}
} if (sw2_tp == "Day") {
# Create `x_time` for all used years
req_ts_days <- ((
from = ((years_used), 1, 1, tz = "UTC"),
to = ((years_used), 12, 31, tz = "UTC"),
= "1 day"
))
x_time <- create_sw2simtime(n = (req_ts_days))
# nolint start: extraction_operator_linter.
x_time, "Year"] <- 1900 + req_ts_days$year
x_time, "Month"] <- 1 + req_ts_days$mon
x_time, "Day"] <- 1 + req_ts_days$yday
# nolint end
# Apparently, SW2 output can be generated with incorrect leap/nonleap-years;
# if so, then the last simulated day of a nonleap year would
# incorrectly be the 366-th day (which in reality doesn't exist) and
# for which `as.POSIXlt()` correctly produces an NA (with a warning);
# similarly, the last simulated day of a leap year would
# incorrectly be the 365-th day (and the 366th-day would be missing)
# --> add simulated non-existing leap-days &
# remove not simulated existing leap years
tmp_sim_seq <- (xt, "Year"], "-", xt, "Day"])
sim_ts_days <- (tmp_sim_seq, = "%Y-%j", tz = "UTC")
has_leap_issues <- (sim_ts_days)
if (has_leap_issues) {
# Create `x_time` for simulated time period
x_time2 <- create_sw2simtime(n = (xt))
x_time2, "Year"] <- xt, "Year"]
x_time2, "Day"] <- xt, "Day"]
# nolint start: extraction_operator_linter.
x_time2, "Month"] <- sim_ts_days$mon + 1
# nolint end
x_time2[is.na(sim_ts_days), "Month"] <- 12
# Add requested but not simulated years from `x_time`
if (has_req_yrs) {
ids_notsim <- x_time, "Year"] %in% (req_years, years_sim)
tmp <- (
x_time[ids_notsim, , = ],
x_time2,
= ,
=
)
x_time <- tmp[order(tmp"Year"]], tmp"Day"]]), , = ]
(x_time) <-
} {
# All used years are simulated years
x_time <- x_time2
}
}
if (has_req_yrs) {
ids_sim <-
(x_time, "Year"], "-", x_time, "Day"]) %in% tmp_sim_seq
}
}
if (has_req_yrs) {
ids_req <- x_time, "Year"] %in% req_years
x_time[ids_sim & ids_req, "mode"] <- "sim_keep"
x_time[ids_sim & !ids_req, "mode"] <- "sim_discard"
x_time!ids_sim & ids_req, "mode"] <- "nosim"
} {
x_time, "mode"] <- "sim_keep"
}
(
!(x_time, "mode"]),
(x_time) >= (xt)
)
x_time
}
# Load rSOILWAT2 output
# A small timing example suggests that loading simulation data from zipped
# archives is about 10% slower (6.8 milliseconds) for a
# simulation output rda-file of 9 MB from an archive with 21 scenarios
load_sw2_rda <- (path, fname, zipped_runs = ) {
sw2_data <- (parent = ())
con <- if (zipped_runs) {
(path, filename = fname)
} {
(path, fname)
}
( = con, envir = sw2_data)
if (zipped_runs) {
(con) # make sure connection is closed
}
sw2_data
}
collect_sw2_sim_data <- (
path,
name_sw2_run,
id_scen,
years,
output_sets = (
(
sw2_tp = ("Day", "Month", "Year"),
sw2_outs = ,
sw2_vars = ,
varnames_are_fixed =
)
),
fail = ,
zipped_runs =
) {
n_sets <- (output_sets)
for (k (n_sets)) {
out <- output_sets[k]]
(
(out"sw2_tp"]]) == 1,
out"sw2_tp"]] %in% ("Day", "Month", "Year"),
(out"sw2_outs"]]) %in% (1L, (out"sw2_vars"]])),
(out"varnames_are_fixed"]]) %in% (1L, (out"sw2_vars"]]))
)
}
#--- Load rSOILWAT2 output object: `runDataSC`
sim_data <- load_sw2_rda(
path = (path, name_sw2_run),
fname = ("sw_output_sc", id_scen, ".RData"),
zipped_runs = zipped_runs
)
#--- Extract variables
res <- ( = "list", = n_sets)
if (((output_sets)) > 0) {
(res) <- (output_sets)
}
for (k1 (n_sets)) {
out <- output_sets[k1]]
nvars <- (out"sw2_vars"]])
sw2_outs <- (out"sw2_outs"]], nvars)
fixed <- (out"varnames_are_fixed"]], nvars)
# Extract slots
x <- (
sw2_outs,
(ko) ((sim_data"runDataSC"]], ko), out"sw2_tp"]])
)
x_vals <- ()
# Subset columns
for (k2 (nvars)) {
if (fixed[k2]) {
x_vals[k2]] <- (
x[k2]], out"sw2_vars"]][k2], = fixed[k2]],
silent =
)
} {
tmp <- (out"sw2_vars"]][k2], (x[k2]]), value = )
if ((tmp) > 0) {
x_vals[k2]] <- x[k2]], tmp, = fixed[k2]]
(x_vals[k2]]) <- tmp
} {
x_vals[k2]] <- (("variables not found"), silent = )
}
}
if ((x_vals[k2]], "try-error")) {
if (fail) {
(x_vals[k2]])
} {
x_vals[k2]] <- (
= ((x[k2]]), 1),
= (, out"sw2_vars"]][k2])
)
}
}
}
nv <- (out"sw2_vars"]])
(x_vals) <- if ((nv)) {
out"sw2_vars"]]
} {
tmp <- nv
ids <- (nv) == 0
if ((ids)) tmp[ids] <- out"sw2_vars"]][ids]
tmp
}
#--- Extract time
# `determine_sw2_sim_time()` is memoized and a missing `years` passed to
# argument `req_years` doesn't work correctly
# (see https://github.com/r-lib/memoise/issues/19)
x_time <- (
xt = x[1]],
req_years = if ((years)) years,
sw2_tp = out"sw2_tp"]]
)
n_nosim <- (x_time, "mode"] == "nosim")
#--- Add entries for requested but not simulated time steps "nosim"
if (n_nosim > 0) {
n_sim <- (x_time) - n_nosim
ids <- if ((x_time[1, "mode"] == "nosim")) {
# "nosim" occurs before "sim_keep"
((, n_nosim), (n_sim))
} {
# "nosim" occurs after "sim_keep"
((n_sim), (, n_nosim))
}
x_vals <- (
x_vals,
(x) if (((x))) x[ids] x[ids, , = ]
)
}
#--- Removes entries for un-requested but simulated time steps "sim_discard"
ids <- (x_time, "mode"] == "sim_discard")
if ((ids) > 0) {
x_time <- x_time-ids, , = ]
x_vals <- (
x_vals,
(x) if (((x))) x-ids] x-ids, , = ]
)
}
res[k1]] <- (
= x_time,
values = x_vals
)
}
res
}
