README.Rmd

---
output: github_document
editor_options: 
  chunk_output_type: console
  markdown: 
    wrap: 72
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

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

# CroPlotR

<!-- badges: start -->

[![Project Status: WIP – Initial development is in progress, but there
has not yet been a stable, usable release suitable for the
public.](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip)
[![Codecov test
coverage](https://codecov.io/gh/SticsRPacks/CroPlotR/branch/master/graph/badge.svg)](https://app.codecov.io/gh/SticsRPacks/CroPlotR?branch=master)
[![R build
status](https://github.com/SticsRPacks/CroPlotR/workflows/R-CMD-check/badge.svg)](https://github.com/SticsRPacks/CroPlotR/actions)
[![DOI](https://zenodo.org/badge/263962392.svg)](https://zenodo.org/badge/latestdoi/263962392)

<!-- badges: end -->

`CroPlotR` aims at the standardization of the process of analyzing the
outputs from crop models such as
[STICS](https://www6.paca.inrae.fr/stics_eng/),
[APSIM](https://www.apsim.info/) or really any model.

Its use does not need any particular adaptation if your model has been
wrapped with the [CroptimizR](https://github.com/SticsRPacks/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](https://github.com/SticsRPacks/CroPlotR).

## Table of Contents

-   [1. Installation](#1-installation)
-   [2. Examples](#2-examples)
    -   [2.1 Plotting](#21-plotting)
        -   [2.1.1 Dynamic plots](#211-dynamic-plots)
        -   [2.1.2 Scatter plots](#212-scatter-plots)
        -   [2.1.3 Group comparison](#213-group-comparison)
        -   [2.1.4 Plot saving](#214-plot-saving)
        -   [2.1.5 Plot extracting](#215-plot-extracting)
    -   [2.2 Statistics](#22-statistics)
        -   [2.2.1 Dynamic plots](#221-simple-case)
        -   [2.2.2 Several groups](#222-several-groups)
        -   [2.2.3 Statistics plot](#223-statistics-plot)
    -   [2.3 Data manipulation](#23-data-manipulation)
-   [3. Tools](#3-tools)
    -   [3.1 ggplotly](#31-ggplotly)
    -   [3.2 patchwork](#32-patchwork)
-   [4. Help](#4-help)
-   [5. Citation](#5-Citation)

## 1. Installation

You can install the released version of CroPlotR from
[Github](https://github.com/SticsRPacks/CroPlotR) either using
`devtools` or the lightweight `remotes` package:

-   With `devtools`

```{r eval=FALSE}
devtools::install_github("SticsRPacks/CroPlotR@*release")
```

-   With `remotes`

```{r eval=FALSE}
# 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()`](https://sticsrpacks.github.io/CroPlotR/reference/plot.cropr_simulation.html)
(and its alias `autoplot()`), and one for the statistics
[`summary()`](https://sticsrpacks.github.io/CroPlotR/reference/summary.cropr_simulation.html).
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](#23-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](https://github.com/SticsRPacks/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](https://sticsrpacks.github.io/CroptimizR/) (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:

```{r}
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:

```{r}
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

```{r}
p
```

or simply

```{r, eval=FALSE}
plot(sim, obs = obs)
```

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

```{r}
names(p)
```

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

```{r}
p$`IC_Wheat_Pea_2005-2006_N0`
```

or index:

```{r, eval=FALSE}
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.

```{r}
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.

```{r}
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](https://github.com/SticsRPacks/CroPlotR/issues/2)). If you
> want to do so, you are encouraged to scale before the plotting
> function, and to add a second axis using
> [sec_axis](https://ggplot2.tidyverse.org/reference/sec_axis.html) on
> the resulting plot.

#### 2.1.2 Scatter plots

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

```{r}
# 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:

```{r}
# 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:

```{r}
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.

```{r}
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.

```{r}
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.

```{r}
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.

```{r}
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.

```{r}
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:

```{r}
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.

```{r}
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.

```{r}
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:

```{r}
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:

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r}
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:

```{r eval=FALSE}
summary(sim, obs = obs, all_situations = FALSE)
```

```{r echo=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.

```{r eval=FALSE}
summary(sim, obs = obs, all_situations = TRUE)
```

```{r echo=FALSE}
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).

```{r eval=FALSE}
summary(sim, sim2, obs = obs)
```

```{r echo=FALSE}
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:

```{r eval=FALSE}
summary("New version" = sim, original = sim2, obs = obs)
```

```{r echo=FALSE}
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:

```{r eval=FALSE}
summary(
  "New version" = sim, original = sim2, obs = obs,
  stats = c("R2", "nRMSE")
)
```

```{r echo=FALSE}
s <- summary(
  "New version" = sim,
  original = sim2,
  obs = obs,
  stats = c("R2", "nRMSE")
)

knitr::kable(s)
```

Please read the help from
[`summary.cropr_simulation()`](https://sticsrpacks.github.io/CroPlotR/reference/summary.cropr_simulation.html)
and
[`predictor_assessment()`](https://sticsrpacks.github.io/CroPlotR/reference/predictor_assessment.html).

#### 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`.

```{r}
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`.

```{r}
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):

```{r}
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.

```{r}
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.

```{r}
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.

```{r}
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](https://CRAN.R-project.org/package=dplyr),
[tidyr](https://CRAN.R-project.org/package=tidyr) or
[tibble](https://CRAN.R-project.org/package=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.

```{r}
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:

```{r}
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.

```{r, eval = FALSE}
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.

```{r}
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:

```{r eval=FALSE}
?plot.cropr_simulation

?plot.statistics
```

-   statistics:

```{r eval=FALSE}
?summary.cropr_simulation
```

If you have any problem, please [fill an
issue](https://github.com/SticsRPacks/CroPlotR/issues) 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.