In SticsRPacks/CroPloteR: A Package to Analyze Crop Model Simulations Outputs with Plots and Statistics

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

CroPlotR

CroPlotR aims at the standardization of the process of analyzing the outputs from crop models such as STICS, APSIM or really any model.

Its use does not need any particular adaptation if your model has been wrapped with the CroptimizR package.

If you want to be notified when a new release of this package is made, you can tick the Releases box in the "Watch / Unwatch => Custom" menu at the top right of this page.

Table of Contents

• 1. Installation
• 2. Examples
  • 2.1 Plotting
    • 2.1.1 Dynamic plots
    • 2.1.2 Scatter plots
    • 2.1.3 Group comparison
    • 2.1.4 Plot saving
    • 2.1.5 Plot extracting
  • 2.2 Statistics
    • 2.2.1 Dynamic plots
    • 2.2.2 Several groups
    • 2.2.3 Statistics plot
  • 2.3 Data manipulation
• 3. Tools
  • 3.1 ggplotly
  • 3.2 patchwork
• 4. Help
• 5. Citation

1. Installation

You can install the released version of CroPlotR from Github either using devtools or the lightweight remotes package:

• With devtools

devtools::install_github("SticsRPacks/CroPlotR@*release")

• With remotes

# install.packages("remotes")
remotes::install_github("SticsRPacks/CroPlotR@*release")

Normally, all the package dependencies will be installed for CRAN packages.

2. Examples

At the moment, only one function is exported for plots plot() (and its alias autoplot()), and one for the statistics summary(). These functions should be the only one you need for all your plots and summary statistics. Additional ones are provided to simplify the manipulation of simulated data (see 2.3 Data manipulation).

In the following, an example using the STICS crop model is presented. If you want to use another model for which a wrapper has been designed for the CroptimizR package, just consider defining the sim variable used in the examples below as sim <- result$sim_list, where result is the list returned by your model wrapper. Examples of use of CroPlotR with Stics and APSIM model wrappers can be found in CroptimizR's website (see Articles tab).

In the following example a simulation of three situations (called USM in STICS) with their observations is used:

• an intercrop of Wheat and pea
• a Pea in sole crop
• a Wheat in sole crop

Let's import the simulation and observation data:

library(CroPlotR)

# Importing an example with three situations with observation:
workspace <- system.file(
  file.path("extdata", "stics_example_1"),
  package = "CroPlotR"
)

situations <- SticsRFiles::get_usms_list(
  file = file.path(workspace, "usms.xml")
)

sim <- SticsRFiles::get_sim(
  workspace = workspace,
  usms_file = file.path(workspace, "usms.xml")
)

obs <- SticsRFiles::get_obs(
  workspace =  workspace,
  usm = situations,
  usms_file = file.path(workspace, "usms.xml")
)

2.1 Plotting

2.1.1 Dynamic plots

Here is an application of dynamic plots for the 3 situations:

p <- plot(sim, obs = obs)

Note that the obs argument is explicitly named. This is because the first argument of the function is ... (we'll see why in a minute).

The plot function returns a named list of ggplot objects.

To plot all of them, just do

p

or simply

plot(sim, obs = obs)

In this case, the elements of the list take the name of the situations.

names(p)

To plot only one of the graph, access it using its name:

p$`IC_Wheat_Pea_2005-2006_N0`

or index:

p[[1]]

It is possible to aggregate plots of multiple situations on the same graph when situations follow one another over time. This can be done using the successive parameter.

workspace <- system.file(
  file.path("extdata", "stics_example_successive"),
  package = "CroPlotR"
)

situations <- SticsRFiles::get_usms_list(
  file = file.path(workspace, "usms.xml")
)

sim_rot <- SticsRFiles::get_sim(
  workspace = workspace,
  usm = situations,
  usms_file = file.path(workspace, "usms.xml")
)

plot(
  sim_rot,
  var = c("resmes", "masec_n"),
  successive = list(list("demo_Wheat1", "demo_BareSoil2", "demo_maize3"))
)

We can also overlay variables thanks to the "overlap" parameter with dynamic plots.

plot(sim, obs = obs, overlap = list(list("lai_n", "masec_n")))

Note that it is not possible to scale the variables right now from the plot function (see issue). If you want to do so, you are encouraged to scale before the plotting function, and to add a second axis using sec_axis on the resulting plot.

2.1.2 Scatter plots

Here are the same plots, but presented as scatter plots:

# Only plotting the first situation for this one:
plots <- plot(sim, obs = obs, type = "scatter", all_situations = FALSE)
plots$`IC_Wheat_Pea_2005-2006_N0`

Residues can also be represented against observations:

# Only plotting the first situation again:
plots <- plot(
  sim,
  obs = obs,
  type = "scatter",
  select_scat = "res",
  all_situations = FALSE
)

plots[[1]]

All these data can also be represented with a single graph for all situations:

plot(sim, obs = obs, type = "scatter", all_situations = TRUE)

When plotting residual scatter plots, reference_var allows to choose the reference variable on the x-axis. Thus, the observations or simulations of this reference variable (to be chosen by suffixing the variable name by "_obs" or "_sim") will be compared to the residuals of each of the variables.

plot(
  sim,
  obs = obs,
  type = "scatter",
  select_scat = "res",
  all_situations = TRUE,
  reference_var = "lai_n_sim"
)

The points on the graphs can be shown in different shapes to differentiate between situations when all_situations = TRUE. If desired, the names of the situations can be displayed.

plot(
  sim,
  obs = obs[c(2, 3)],
  type = "scatter",
  all_situations = TRUE,
  shape_sit = "txt"
)

As you can see, this can quickly become unreadable depending on the number of points and length of situation names; That is why you can simply assign a different symbol to each situation.

plot(
  sim,
  obs = obs,
  type = "scatter",
  all_situations = TRUE,
  shape_sit = "symbol"
)

It is also possible to represent a group of situations with the same symbol when, for example, clusters are identified.

plot(
  sim,
  obs = obs,
  type = "scatter",
  all_situations = TRUE,
  shape_sit = "group",
  situation_group = list(list("SC_Pea_2005-2006_N0", "SC_Wheat_2005-2006_N0"))
)

You can also name your situation_group list and thus customize (e.g shorten) the plot legend.

plot(
  sim,
  obs = obs,
  type = "scatter",
  all_situations = TRUE,
  shape_sit = "group",
  situation_group = list(
    "Two Single Crops" = list("SC_Pea_2005-2006_N0", "SC_Wheat_2005-2006_N0"))
)

By default, all variables are returned by plot(), but you can filter them using the var argument:

plot(sim, obs = obs, type = "scatter", all_situations = TRUE, var = c("lai_n"))

Error bars related to observations can also be added to the graph using the obs_sd parameter which must be of the same shape as obs. In our example, we will create a false data frame with the only purpose of having a preview of the result. To have 95% confidence, the error bar is equal to two standard deviations on each side of the point.

obs_sd <- obs
names_obs <- names(obs_sd$`SC_Pea_2005-2006_N0`)
obs_sd$`SC_Pea_2005-2006_N0`[, !(names_obs %in% c("Date", "Plant"))] <-
  0.05 * obs_sd$`SC_Pea_2005-2006_N0`[, !(names_obs %in% c("Date", "Plant"))]
obs_sd$`SC_Wheat_2005-2006_N0`[, !(names_obs %in% c("Date", "Plant"))] <-
  0.2 * obs_sd$`SC_Wheat_2005-2006_N0`[, !(names_obs %in% c("Date", "Plant"))]

plot(sim, obs = obs, obs_sd = obs_sd, type = "scatter", all_situations = TRUE)

2.1.3 Group comparison

We can compare groups of simulations alongside by simply adding the simulations objects one after the other (that is why the first argument of the function is ...). Group simulations can be the results of simulations from different model versions, or simulations with different parameter values.

workspace2 <- system.file(
  file.path("extdata", "stics_example_2"),
  package = "CroPlotR"
)

sim2 <- SticsRFiles::get_sim(
  workspace = workspace2,
  usms_file = file.path(workspace2, "usms.xml")
)

plot(sim, sim2, obs = obs, all_situations = FALSE)

Here only one plot is outputted because workspace2 only contains the intercrop situation.

We can also name the corresponding group in the plot by naming them while passing to the plot() function:

plot(
  "New version" = sim,
  original = sim2,
  obs = obs,
  type = "scatter",
  all_situations = FALSE
)

2.1.4 Plot saving

The plots can be saved to disk using the save_plot_png() function as follows:

plots <- plot("New version" = sim, original = sim2, obs = obs, type = "scatter")

save_plot_png(plot = plots, out_dir = "path/to/directory", suffix = "_scatter")

# or by piping:
plots <- plot(
  "New version" = sim,
  original = sim2,
  obs = obs,
  type = "scatter"
) %>%
save_plot_png(., out_dir = "path/to/directory", suffix = "_scatter")

They can also be saved using the save_plot_pdf() function that which, from a list of ggplots, generates a pdf file. If the file_per_var parameter is TRUE, in this case the function generates one pdf file per variable.

plots <- plot(sim, obs = obs)

save_plot_pdf(plot = plots, out_dir = "path/to/directory", file_per_var = FALSE)

2.1.5 Plot extracting

When we have plots with several variables and several situations, the extract_plot function allows to keep the situations and variables that we need.

In the following example, we want to extract the intercrop situation and the "masec_n" variable.

plots <- plot(sim, obs = obs, type = "scatter", all_situations = FALSE)

extract_plot(
  plots,
  situation = c("IC_Wheat_Pea_2005-2006_N0"), var = c("masec_n")
)

2.2 Statistics

2.2.1 Simple case

Here is an application of summary statistics for the 3 situations:

summary(sim, obs = obs, all_situations = FALSE)

s <- summary(sim, obs = obs, all_situations = FALSE)
knitr::kable(s)

Note that as for the plot() function the obs argument is explicitly named. This is because the first argument of the function is ... to be able to compare groups (i.e. model versions or simulation with different parameter values). In this example, a message warns the user because some observed values have a zero value which causes a division by zero in the calculation of certain statistical criteria, these values are therefore filtered for the calculation of these criteria.

And as for the plot() function again, it is possible to compute the statistical criteria for all situations at once.

summary(sim, obs = obs, all_situations = TRUE)

s <- summary(sim, obs = obs, all_situations = TRUE)
knitr::kable(s)

2.2.2 Several groups

We can get statistics for each group of simulations by simply adding the simulations objects one after the other (as for the plot() function).

summary(sim, sim2, obs = obs)

s <- summary(sim, sim2, obs = obs)
knitr::kable(s)

We can also name the corresponding group in the plot by naming them while passing to the summary() function:

summary("New version" = sim, original = sim2, obs = obs)

s <- summary("New version" = sim, original = sim2, obs = obs)
knitr::kable(s)

By default, all statistics are returned by summary, but you can filter them using the stat argument:

summary("New version" = sim, original = sim2, obs = obs,
        stats = c("R2", "nRMSE"))

s <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE")
)

knitr::kable(s)

Please read the help from summary.cropr_simulation() and predictor_assessment().

2.2.3 Statistics plot

It is also possible to plot the statistics:

In a rather obvious way, the resulting graph will take into account all the situations simultaneously or not according to the parameter given to summary. Here is an example with all_situations = FALSE.

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE"),
  all_situations = FALSE
)
plot(stats)

And here is an example with all_situations = TRUE.

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE"),
  all_situations = TRUE
)

plot(stats)

We can choose to plot either the group or the situation in x (and the other is used for grouping and colouring):

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE"),
  all_situations = FALSE
)

plot(stats, xvar = "situation", title = "Situation in X")

In the previous examples, each line corresponds to a statistical criterion. These can also be stacked.

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("pMSEs", "pMSEu"),
  all_situations = FALSE
)

plot(stats, xvar = "situation", title = "Stacked columns", group_bar = "stack")

Or put side by side.

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("pMSEs", "pMSEu"),
  all_situations = FALSE
)

plot(
  stats,
  xvar = "situation",
  title = "Side-by-side columns",
  group_bar = "dodge"
)

To compare different versions on a single criterion, the function produces a radar graph like the following one.

sim$`SC_Pea_2005-2006_N0`$mafruit <-
  (15 / 10) * sim$`SC_Pea_2005-2006_N0`$masec_n
sim$`SC_Wheat_2005-2006_N0`$mafruit <-
  (15 / 20) * sim$`SC_Wheat_2005-2006_N0`$masec_n
sim2$`IC_Wheat_Pea_2005-2006_N0`$mafruit <-
  sim2$`IC_Wheat_Pea_2005-2006_N0`$masec_n
obs$`IC_Wheat_Pea_2005-2006_N0`$mafruit <-
  (12 / 10) * obs$`IC_Wheat_Pea_2005-2006_N0`$masec_n
obs$`SC_Pea_2005-2006_N0`$mafruit <-
  (18 / 10) * obs$`SC_Pea_2005-2006_N0`$masec_n
obs$`SC_Wheat_2005-2006_N0`$mafruit <-
  (15 / 12) * obs$`SC_Wheat_2005-2006_N0`$masec_n

stats <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE"),
  all_situations = TRUE
)

plot(
  stats,
  type = "radar",
  crit_radar = "nRMSE",
  title = "Radar chart : nRMSE"
)

2.3 Data manipulation

Observation lists can easily be handled using e.g. dplyr, tidyr or tibble packages.

The use of these packages on simulated data as returned by CroptimizR model wrappers is sometimes prevented by their attribute cropr_simulation. To easily manipulate simulated data we thus provide two functions for (i) binding rows of data simulated on different situations in a single data.frame or tibble and (ii) go back to the original (cropr) format by splitting this single data.frame or tibble.

df <- bind_rows(sim)
head(df)

The resulting data.frame/tibble can then easily be manipulated using standard R packages. The column situation contains the name of the corresponding situation (as given in the named list sim).

To go back to the original format of simulated data handled by CroPlotR, use the split_df2sim function:

sim_new <- split_df2sim(df)
lapply(sim_new, head)

3. Tools

3.1 ggplotly

The ggplotly function in plotly library makes it very easy to create interactive graphics from a ggplot. Do not hesitate to call it with your plot and move your mouse over the graph to discover the features of this function.

library(plotly)

ggplotly(plot(sim, obs = obs, type = "dynamic")[[1]])

3.2 patchwork

There is also the patchwork library that allows you to easily combine several ggplot into one.

library(patchwork)

plot1 <- plot(sim, obs = obs, type = "scatter", var = "lai_n")[[1]]
plot2 <- plot(sim, obs = obs, var = "lai_n")[[1]]
plot3 <- plot(sim, obs = obs, type = "scatter", var = "masec_n")[[1]]
plot4 <- plot(sim, obs = obs, var = "masec_n")[[1]]

plot1 + plot2 + plot3 + plot4 + plot_layout(ncol = 2)

4. Help

You can find help for the functions directly using the name of the function followed by the class of the object you need the method for:

• plot:

?plot.cropr_simulation

?plot.statistics

• statistics:

?summary.cropr_simulation

If you have any problem, please fill an issue on Github.

5. Citation

If you have used this package for a study that led to a publication or report, please cite us. You can either use the citation tool from Github if you used the last version, or use citation("CroPlotR") from R otherwise.

SticsRPacks/CroPloteR documentation built on April 1, 2024, 9:25 a.m.