Getting started with apsimeval

The interface

Everything below can also be done from a browser interface:

library(apsimeval)
run_app()                                   # opens on the working directory
run_app(out_dir = "~/apsim/sundarbans/outputs")
run_app(out_dir = "outputs", project = "sundarbans.apsimproj")
run_app(out_dir = "outputs", poll_ms = 1000, launch.browser = FALSE)

run_app() needs the optional packages: install.packages(c("shiny", "DT", "bslib")). Leave it open while APSIM runs — with Live sync on it re-reads only the files that changed and redraws, keeping your current settings.

Reading output

sim <- read_out_dir("outputs")              # every .out in the folder
apsim_vars("outputs/N120.out")              # variables and units, header only
sim <- read_out_dir("outputs", vars = c("wagt", "wrr"))

Treatments from simulation names

detect_factors(sim$sim)                     # guesses the separator and fields
preview_factors(sim$sim, parts = c("site", "crop", "water", "N"))
sim <- set_factors(sim, parts = c("site", "crop", "water", "N"),
                   numeric_from = "N")

Irregular names take a regex, hopeless ones a lookup table:

sim <- set_factors(sim, regex = "(?<site>[A-Za-z]+)[-_]N(?<N>\\d+)")
map <- factor_template(sim$sim, factors = c("water", "N"))   # fill in and pass back
sim <- set_factors(sim, map = map)

Calibration and validation

obs <- map_obs(read_obs_raw("observed.xlsx"),
               date_col = "Sampling_Date", value_cols = c("Grain_yield", "Biomass"),
               treatment_col = "Treatment",
               # Spreadsheet headings are not APSIM variable names. Without
               # this, nothing pairs and nothing plots.
               var_map = c(Grain_yield = "wrr", Biomass = "wagt"))

check_variables(obs, sim)                   # catches variable name mismatches
check_levels(obs, sim)                      # catches treatment label mismatches

long <- as_long(sim, vars = "wrr")
pr   <- pair_obs(long, obs, tol_days = 3)

cut <- suggest_cutoff(pr)                   # first season calibrates
pr  <- split_phase(pr, "date", at = cut)

# Statistics are never pooled across variables: RMSE would mix units and R2
# would be inflated by the gap between the variables rather than model skill.
gof("obs", "pred", data = pr, by = "variable")
gof("obs", "pred", data = pr, by = c("variable", "phase"))
gof("obs", "pred", data = pr, by = c("variable", "phase", "N"))

The series stays continuous across the cutoff:

plot_ts(long, obs = obs, cutoff = cut)
# Colouring is a presentation choice and does not affect the statistics; the
# corner box follows `stats_by`, which defaults to the faceting column.
plot_one2one(pr, group = "N", stats_by = character(0))   # pooled box
plot_one2one(pr, group = "N", facet = "phase", sd_band = TRUE)

Scenarios and sensitivity

ag <- aggregate_sim(sim, by = c("sim", "year", "water", "N"),
                    vars = c(wrr = "max"), crop_only = TRUE)
plot_exceedance(ag, "wrr", group = "N", facet = "water", ref = 7000)
plot_bar(ag, "wrr", "N", letters = TRUE)

sens_anova(ag, "wrr", c("site", "water", "N"), order = 2)
sens_oat(ag, "wrr", c("site", "water", "N"))

Figures

set_labels(list(wrr = "Grain yield"))
set_palette(c(N0 = "#D55E00", N120 = "#0072B2"))
set_colour_mode("grey")                     # "colour", "grey", or "mono"

fig <- build_figure(list(
  panel_spec("one2one", group = "N", facet = "phase"),
  panel_spec("exceedance", var = "wrr", group = "N", ref = 7000),
  panel_spec("bar", var = "wrr", group = "N", letters = TRUE)
), sim = sim, paired = pr, ncol = 3)

save_pub(fig, "figure_2", width = "double", height = 80,
         dpi = 600, format = c("tiff", "pdf"))

Projects

p <- new_project("sundarbans_rice", "outputs")
p <- set_project_factors(p, sim)
p$labels <- list(wrr = "Grain yield"); p$colour_mode <- "grey"
save_project(p, "sundarbans_rice.apsimproj")

p <- load_project("sundarbans_rice.apsimproj")
sim <- apply_factors(read_out_dir("outputs"), p)   # after a re-run

Replicated observed data

Field data usually arrives with one row per plot:

Sampling_Date	Site	Treatment	Replication	Biomass
2016-08-01	Gosaba	N120	1	1000
2016-08-01	Gosaba	N120	2	1100

Name the replication column and the replicates are averaged per date and treatment, keeping the spread for error bars:

raw <- read_obs_raw("observed.xlsx")
attr(raw, "guess")                # suggested date / treatment / rep / value columns

obs <- map_obs(raw, date_col = "Sampling_Date", value_cols = "Biomass",
               treatment_col = c("Site", "Treatment"),   # joined to match sim names
               rep_col = "Replication")

obs gains obs_sd and obs_n. A date with one replicate gets obs_sd = NA rather than zero, so no misleading zero-length bar is drawn.

plot_ts(long, obs = obs, cutoff = cut)          # error bars drawn automatically
plot_ts(long, obs = obs, cutoff = cut, obs_sd = NULL)   # suppress them
plot_one2one(pr, obs_sd = "obs_sd")             # horizontal bars on the scatter

If your sheet already holds means and a separate SD column (Biomass_SD), leave rep_col empty: the SD column is detected and paired automatically.

Timeline graphs

plot_series() handles several variables and several treatments at once.

Different units go on a secondary axis:

long <- as_long(sim, vars = c("wagt", "pond"), keep = c("water", "N"))
plot_series(long, vars = c("wagt", "pond"), secondary = "pond")

Or give each variable its own row, which avoids the dual-axis problem entirely:

plot_series(long, vars = c("wagt", "wrr", "pond"), layout = "stacked")

All treatments side by side — the usual way to judge a treatment effect, since the shapes are compared rather than just the endpoints:

plot_series(long, vars = c("wagt", "pond"), secondary = "pond",
            facet_by = "sim", ncol = 2)             # every simulation
plot_series(long, vars = "wagt", facet_by = "water") # or a factor column
plot_series(long, vars = "wagt", facet_by = "water",
            colour_by = "treatment")                # colour instead of panels

Combine layouts for a variable-by-treatment grid, and add the cutoff, observed points and replicate error bars as usual:

plot_series(long, vars = c("wagt", "wrr"), layout = "stacked", facet_by = "sim")
plot_series(long, vars = "wagt", facet_by = "sim", obs = obs, cutoff = cut)

Note keep = in as_long(): treatment columns have to survive the pivot or there is nothing left to facet by.

A caution on the secondary axis. Two series with different units can be made to look correlated by the scaling alone, so layout = "stacked" is the safer choice when the variables are not meant to be compared point by point. The right axis is always labelled in original units, and legend entries drawn against it are marked.

Observed data on a subset of the variables

Observations rarely cover every plotted variable — biomass and yield are measured, AWD ponding usually is not. Pass whatever you have; variables without observations are drawn as simulated-only lines and the omission is reported:

plot_series(long, vars = c("wagt", "wrr", "pond"), secondary = "pond", obs = obs)
#> No observed data for: pond (simulated only).

Each variable’s points carry their own marker shape and their own replicate error bars, computed exactly as elsewhere in the package (rep_col in map_obs(), or an existing obs_sd column). Points belonging to a secondary axis variable are rescaled with the line, and their error bars are rescaled as lengths rather than positions.

When panelling by treatment, observed records need to know which panel they belong to. If the observed data has no column matching facet_by, it is derived from the simulation label:

plot_series(long, vars = "wagt", facet_by = "water", obs = obs)  # obs has only `sim`

Records matching no simulation are dropped and counted in a message rather than being silently discarded.