The SticsRpacks packages

SticsRPacks is a set of packages for managing the STICS soil-crop model from R.

It includes 4 packages:

• SticsRFiles: to manipulate STICS input and output files,
• SticsOnR: to run STICS simulations,
• CroPlotR: to plot models results / observations and compute stats,
• CroptimizR: to estimate model parameters.

SticsRFiles and SticsOnR are dedicated to the STICS model, while CroPlotR and CroptimizR are generic for crop models (i.e. they can be used with any crop model that has an adequate R wrapper).

The websites of the different packages can be found at:

• https://SticsRPacks.github.io/SticsRFiles/
• https://SticsRPacks.github.io/SticsOnR/
• https://SticsRPacks.github.io/CroPlotR/
• https://SticsRPacks.github.io/CroptimizR/

Finding documentation about the packages and their functions

There are different sources of documentation regarding the packages and their functions:

• the home pages of the packages websites (which links are given here-before), that introduce the packages main features,
• the vignettes (Articles tab at the top of the packages websites home pages), that describe in detail the main features of the packages and the use of the associated functions,
• the Reference tab at the top of the packages websites home pages, that lists and quickly describes all the exported functions of the packages, and gives links to their detailed documentation,
• the command ? functionName in the R console, which displays the detailed documentation of a given function.

Note also that completion can be used to list the different arguments of a function: if you type the name of a function in the R console followed by a bracket, e.g. read.csv( and then click on the Tab button of your keyboard, the list of arguments of the function will be displayed.

The SticsRpacks tutorial

This tutorial aims at introducing the use of these different packages with the STICS model. It was designed to be of increasing complexity, it is thus advised to follow it in order. The different sections are however independent and you can thus start it from any part using the menu on the left side (this may however take a little while since results from previous exercises must be automatically generated in this case).

The names and documentation of the functions to be used are given in the first exercises before letting you find them on your own as you go through the exercises.

Most of the exercises in the tutorial requires writing and running a few lines of R code through R console panels like this one:

# Insert your R code here!

# This is the first advise!

# Here's the SOLUTION:
get_var_info(keyword = "roots")

You just have to type the R code in the black part and click on the Run Code green button on the right side to run it. The results will be printed below the panel. You can use these panels to do whatever tests you want to do to explore the features of the packages.

The Start over button re-initialize the panel.

If you don’t find the solution, click on the button Hints, it will give you an advice. If you click again and again several advises will be given (dependgin on the exercise) and finally the expected solution of the exercise. To copy the solution in the panel, click on the Copy to Clipboard blue button on the right and then right-click in the black part of the panel and select “paste”. You can train on the panel above.

STICS training Practical session PW1

The practical sessions of the STICS training are traditionally performed using the JavaStics GUI and Excel.

Let’s see how the first sessions would look like in R using SticsRPacks.

The game is to run a single USM by activating and deactivating water and nitrogen stress in turns and to analyze the impact by looking at dynamic variables and stress indicators.

Let’s first initialize variables and define the setup:

# Let's define the USM and variables of interest to use in this example
usm_name <- "maize_ref_2013"
stress_ind_names <- c("turfac1moy", "turfac2moy", "swfac1moy", "swfac2moy", "inn1moy", "inn2moy")
dyn_var_names <- c("lai_n","zrac","resmes","azomes","QNplante","masec_n")
var_names <- c(stress_ind_names,dyn_var_names)

# Set-up the simulation environment:
wrapper_options = stics_wrapper_options(javastics = javastics_path, workspace = workspace_PW1)

# Generate the files for the simulation:
gen_usms_xml2txt(javastics_path, workspace_PW1, usm = usm_name)

Then, let’s run the simulations for the different configurations:

# Simulate a reference situation with water and nitrogen stress
sim_ws_ns = stics_wrapper(model_options = wrapper_options,
                          param_values = c(codeinnact=1, codeh2oact=1),
                          var = var_names)$sim_list

# Simulate with water stress only
sim_ws = stics_wrapper(model_options = wrapper_options,
                       param_values = c(codeinnact=2),
                       var = var_names)$sim_list

# Simulate with N stress only
sim_ns = stics_wrapper(model_options = wrapper_options,
                       param_values = c(codeh2oact=2),
                       var= var_names)$sim_list

# Simulate with water stress only
sim_without_stress = stics_wrapper(model_options = wrapper_options,
                                   param_values = c(codeinnact=2, codeh2oact=2),
                                   var = var_names)$sim_list

Now we can use the simulated results to plot the dynamic variables of interest for the different configurations:

plot(sim_ws_ns=sim_ws_ns, sim_ws=sim_ws, sim_ns=sim_ns, sim_without_stress=sim_without_stress, var = dyn_var_names )

## $maize_ref_2013

Finally, we can plot the stress indicators averaged over vegetative and reproductive phases:

# Get the final values
df <- CroPlotR::bind_rows(
  sim_ws_ns[[usm_name]] %>% summarise(across(.cols = stress_ind_names, last)) %>%
    bind_cols(Case="sim_ws_ns"),
  sim_ws[[usm_name]] %>% summarise(across(.cols = stress_ind_names, last)) %>%
    bind_cols(Case="sim_ws"),
  sim_ns[[usm_name]] %>% summarise(across(.cols = stress_ind_names, last)) %>%
    bind_cols(Case="sim_ns"),
  sim_without_stress[[usm_name]] %>% summarise(across(.cols = stress_ind_names, last)) %>%
    bind_cols(Case="sim_without_stress"))

# Plot bar graph with ggplot ... sorry, not yet available using CroPlotR ...
df_long <- tidyr::pivot_longer(df, stress_ind_names, names_to = "variables")
ggplot(df_long, aes(fill=Case, y=value, x=variables)) + geom_bar(position="dodge", stat="identity")

… and let the users interpret these results!

The last part may not be straightforward for R newcomers … however a function to filter values of simulated variables at given dates and/or stages will come soon :-).

Parameter optimization: A simple example

Here is an example use-case of an optimization of the stamflax parameter for a wheat crop using observations of Leaf Area Index for one situation (i.e. one USM):

The code to get this is the following:

usm = "wheat"

# First, import the observations and filter them:
obs = get_obs(workspace = workspace_path, usm = usm)
obs_filt = filter_obs(obs, "lai_n", include = TRUE)

# Set-up the simulation environment:
wrapper_options = stics_wrapper_options(javastics = javastics_path, workspace = workspace_path, parallel = TRUE)
optim_options = list(nb_rep = 3, out_dir = workspace_path)

# Define the parameter to optimize and the boundary values:
param_info = list(stamflax = list(lb = 200, ub = 400))

# Generate the files for the simulation:
gen_usms_xml2txt(javastics_path, workspace_path, usm = usm)

# Make the optimization:
optim_res = estim_param(obs_list = obs_filt,
                        model_function = stics_wrapper,
                        model_options = wrapper_options,
                        optim_options = optim_options,
                        param_info = param_info)

# Import the results before/after the simulation and plot them:
sim_after_optim = stics_wrapper(param_values = optim_res$final_values, model_options = wrapper_options)$sim_list
sim_before_optim = stics_wrapper(model_options = wrapper_options)$sim_list

plot(after = sim_after_optim, before = sim_before_optim, obs = obs_filt, var = "lai_n")

A lot is happening here!

But don’t worry if you don’t understand anything, you’ll learn everything you need to know in this tutorial. We just want you to see that you can make a parameter optimization with less than 30 lines of code, including the plots of the improvement of the resulting simulated variable of interest against observation. And the good news is when you understand these lines of code, you’ll be able to do almost anything with the packages, because it doesn’t get any longer or more complex than this.

OK now that we saw what we can do, let’s find out how step by step.

Get usm names

The function get_usms_list from the SticsRFiles package returns the list of USMs included in a usms.xml file.

Here’s its description as given in the R help panel (? get_usms_list) and on the SticsRFiles package website:

Description

Extracting a usm names list from an usms.xml file

Usage

get_usms_list(
  file,
  usm = NULL,
  usm_path = lifecycle::deprecated(),
  name = lifecycle::deprecated()
)

Arguments

file	Path (including name) of the USM xml file
usm	Vector of USM names (or partial names). Optional, if not provided, the function returns the names of all the USMs included in the given file.
usm_path	usm_path is no longer supported, use file instead.
name	name is no longer supported, use usm instead.

Value

A vector of usm names

Examples

path <- get_examples_path(file_type = "xml")

usms_list <- get_usms_list(file = file.path(path, "usms.xml"))

usms_list <- get_usms_list(file = file.path(path, "usms.xml"), usm = c("usm1", "usm2"))

• You will use in the following the get_usms_list function to get the list of USMs included in the “example” workspace of the JavaStics distribution installed by this tutorial. The path to this workspace is stored in the R object workspace_path (be careful, do not confuse the path of the JavaStics example workspace, stored in the R object workspace_path, and the path of the examples of the SticsRFiles package, retrieved using the get_examples_path function in the function documentation displayed above!).

First, display the path of the JavaStics “example” workspace, and list its content using the list.files function:

# Insert your R code here!

# print the object workspace_path to display the path

# apply the list.files function to workspace_path for listing its files

print(workspace_path)

list.files(workspace_path)

Now, display the list of USMs defined in the workspace workspace_path

# Insert your R code here!

# Use the file.path function to build full path to a file from directory names and filenames: file.path(path,filename)

# SOLUTION:
get_usms_list(file = file.path(workspace_path, "usms.xml"))

• Get the list of cover-crops USMs (their name starts with "cc_").

# Insert your R code here!

# The argument "usm" is particularly useful in this case ...

#SOLUTION:
get_usms_list(file = file.path(workspace_path, "usms.xml"), usm = "cc_")

Get parameters values from JavaStics input files

Parameters values included in JavaStics input files can be retrieved using the function get_param_xml from the SticsRFiles package:

Description

Extracting parameter values for a list of xml files and parameters

Usage

get_param_xml(
  file,
  param = NULL,
  select = NULL,
  select_value = NULL,
  xml_file = lifecycle::deprecated(),
  param_name = lifecycle::deprecated(),
  value = lifecycle::deprecated(),
  ...
)

Arguments

file	Vector of the xml file paths from which parameters values must be extracted
param	Vector of parameter names. Optional, if not provided, the function returns information for all parameters.
select	node name or attribute name to use for selection (optional, default to no selection)
select_value	Vector of values used for select (see examples). Optional, should be provided only if select is provided.
xml_file	xml_file is no longer supported, use file instead.
param_name	param_name is no longer supported, use param instead.
value	value is no longer supported, use select_value instead.
…	Pass further arguments to get_param_value()

Value

A list of parameter values for each xml_file (a list of list)

Examples

## Not run: 
# Soil file
file <- file.path(get_examples_path(file_type = "xml"), "sols.xml")

# For all soils
get_param_xml(file)
get_param_xml(file, "argi")
get_param_xml(file, c("argi", "norg"))

# For one soil selection
get_param_xml(file, "argi", select = "sol", select_value = "solcanne")

# For soils and parameters vectors
# scalar parameters per soil
get_param_xml(file, c("argi", "norg"),
  select = "sol", select_value = c("solcanne", "solbanane")
)

# vector parameters per soil (5 values, one per soil layer)
get_param_xml(file, c("epc", "HCCF"),
  select = "sol", select_value = c("solcanne", "solbanane")
)

# Crop management file
file <- file.path(get_examples_path(file_type = "xml"), "file_tec.xml")

get_param_xml(file)

# Getting parameters for irrigation (date and quantity)
get_param_xml(file, c("julapI_or_sum_upvt", "amount"))

# Getting all parameters for a given formalism: "irrigation"
get_param_xml(file, select = "formalisme", select_value = "irrigation")

## End(Not run)

• Get the values of clay content after decarbonation (argi) for all soils defined in the sols.xml file of the workspace_path workspace of the JavaStics distribution.

# Insert your R code here!

# The path of the "sols.xml" file of the JavaStics example folder can be defined using the function file.path(workspace_path, "sols.xml")

#SOLUTION:
get_param_xml(file.path(workspace_path, "sols.xml"), c("argi"))

result <- get_param_xml(file.path(workspace_path, "sols.xml"), c("argi"))

• Knowing that the result of the previous exercise has been stored in an R object called result, plot an histogram of the values of argi.

# Insert your R code here!

# Use the hist function to plot an histogram

# Be careful, result is a list: argi is stored in result$sols.xml$argi

#SOLUTION:
hist(result$sols.xml$argi)

• Compute the Available Water Capacity of fine earth for soil "solmais" in mm. The equation is AWC=sum((HCCF-HMINF)*epc*DAF/10), where HCCF is the vector of STICS parameters for field capacity of the different soil layers, HMINF for wilting point, epc for thickness, DAF for soil density, and sum is over the different layers.

# Insert your R code here!

# Use select and value arguments of get_param_xml to extract information for a given soil

# SOLUTION:
result <- get_param_xml(file.path(workspace_path, "sols.xml"), c("HCCF", "HMINF", "epc", "DAF"), select = "sol", select_value = "solmais")
sum((result$sols.xml$HCCF - result$sols.xml$HMINF) * result$sols.xml$epc * result$sols.xml$DAF / 10)

Get weather variables

Weather variables defined in STICS input weather files can be retrieved using the function get_climate_txt from the SticsRFiles package:

Description

Read the meteorology input for STICS ("climat.txt")

Usage

get_climate_txt(
  workspace = getwd(),
  file_name = "climat.txt",
  preserve = TRUE,
  dirpath = lifecycle::deprecated(),
  filename = lifecycle::deprecated()
)

Arguments

workspace	Path of the workspace containing the Stics climate file to read
file_name	The meteorology file name (default to climat.txt).
preserve	Logical, TRUE`` for keeping the stics columns related to date (year, month, day, julian), or only keep the date as a POSIXctotherwise. Default toTRUE’.
dirpath	dirpath is no longer supported, use workspace instead.
filename	filename is no longer supported, use file_name instead.

Value

A data.frame of the input meteorological variables used as input for the STICS model.

Note

The time-related variables are summarised into one POSIXct column named Date.

Examples

library(SticsRFiles)
path <- get_examples_path(file_type = "txt")
Meteo <- get_climate_txt(path)

• Compute the total amount of precipitation that has fallen in Auzeville in 2013 from the file "Auzeville.2013" of the JavaStics Example workspace.

# Insert your R code here!

# Precipitation is stored in the column named ttrr

#SOLUTION:
result <- get_climate_txt(workspace = workspace_path, file_name = "Auzeville.2013")
sum(result$ttrr)

• Compute the amount of precipitation that has fallen between dates "2013-03-21" and "2013-09-20" using the result R object containing the data of the "Auzeville.2013" file.

file_path <- file.path(tmp_dir,"climate.RData")
if (!file.exists(file_path)) {
  result <- get_climate_txt(workspace = workspace_path, file_name = "Auzeville.2013")
  save("result", file=file_path)
} else {
  load(file = file_path)
}

# Insert your R code here!

# Use dplyr::filter(result,Date>"2013-03-21", Date<"2013-09-20") to filter the results wrt the Dates.

#SOLUTION:
filtered_result <- dplyr::filter(result, Date > "2013-03-21", Date < "2013-09-20")
sum(filtered_result$ttrr)

Getting simulations results

The get_sim() function can be used to get simulation results to R from a STICS workspace. Some options allows filtering data, e.g. for given USM names or output variables. Check the function help before trying to use it.

usms <- get_usms_list(file = file.path(workspace_path, "usms.xml"))
sim_done <- length(list.files(pattern = "^mod_s", workspace_path)) >=
  length(usms)

if (!sim_done) {
  print("Running simulations, please wait")
    run_javastics(javastics = javastics_path,
                  workspace = workspace_path, verbose = FALSE)
}

• Get the values of all the outputs simulated during the previous exercise in the example workspace

# Insert your R code here!

# The path of results is stored in the workspace_path object

get_sim(workspace = workspace_path)

• Get the values of all the outputs simulated for USMs banana and Turmeric

# Insert your R code here!

# The argument usm is particularly useful in this case ...

# SOLUTION
get_sim(workspace = workspace_path, usm = c("banana","Turmeric"))

• Get the values of the variable lai_n simulated for USMs wheat and maize

# Insert your R code here!

# Use the argument var for specifying which variable(s) to select ...

# SOLUTION
get_sim(workspace = workspace_path, usm = c("wheat", "maize"),
        var = "lai_n")

Data structures

The simulations (and observations) are stored using a list of data.frames. Each data.frame is used to store the data from a USM, with each row being a value for a given day and each column a variable (e.g. lai_n, the daily lai).

Each element (i.e. data.frame) of the list is named after the USM name. So we can get the simulation results (or the observations) either by indexing the list with an integer value, or by using the $ notation followed by the name of the USM, as you would do to get a variable with a data.frame.

usms <- c("wheat","DurumWheat", "maize", "sorghum")
sim <- get_sim(workspace = workspace_path, usm = usms)

usms <- get_usms_list(file = file.path(workspace_path, "usms.xml"))
sim_done <- length(list.files(pattern = "^mod_s", workspace_path)) >=
  length(usms)

if (!sim_done) {
  print("Running simulations, please wait")
    run_javastics(javastics = javastics_path,
                  workspace = workspace_path, verbose = FALSE)
}
usms <- c("wheat","DurumWheat", "maize", "sorghum")
sim <- get_sim(workspace = workspace_path, usm = usms, verbose = FALSE)
sim_all = CroPlotR::bind_rows(sim)

# Insert your R code here!

# use the names function on the `sim` object

# SOLUTION
names(sim)

# Insert your R code here!

# Index into the `sim` object using `[[`

# SOLUTION
sim[[2]]

# Insert your R code here!

# Index into the `sim` object using `$` and the name of the USM: maize

# SOLUTION
sim$maize

lapply(sim, function(x) max(x$lai_n))

# Insert your R code here!

# Use `bind_rows` on `sim`

# SOLUTION
CroPlotR::bind_rows(sim)

sim_all = CroPlotR::bind_rows(sim)

filter(sim_all, Date == as.POSIXct("2000-09-04", tz = "UTC"))

sim_all%>%
  group_by(situation)%>%
  summarise(max_lai = max(lai_n))

Plotting simulations results

The CroPlotR package can be used to plot simulated results as given by the get_sim function.

• Use the values of the variables lai_n and masec_n simulated for USM wheat already stored in the results R object and plot their dynamics

usms <- get_usms_list(file = file.path(workspace_path, "usms.xml"))
sim_done <- length(list.files(pattern = "^mod_s", workspace_path)) >=
  length(usms)

if (!sim_done) {
  print("Running simulations, please wait")
    run_javastics(javastics = javastics_path,
                  workspace = workspace_path, verbose = FALSE)
}
file_path <- file.path(tmp_dir,"sim-1.RData")
if (!file.exists(file_path)) {
  results <- get_sim(workspace = workspace_path, usm = "wheat",
                     var = c("lai_n", "masec_n"), verbose = FALSE)
  
  save("results", file=file_path)
} else {
  load(file=file_path)
}

# Insert your R code here!

# Use the results object containing get_sim output ...

# Then you can use the plot function as described in the CroPlotR documentation on the website <https://SticsRPacks.github.io/CroPlotR/reference/plot.cropr_simulation.html> or from R using ?plot.cropr_simulation

# SOLUTION
plot(results)

• Then, do the same but plotting both variables on the same graph

# Insert your R code here!

# Use the "overlap" argument ...

# try using : overlap = list(list("lai_n","masec_n"))

# SOLUTION
plot(results, overlap = list(list("lai_n","masec_n")))

Evaluating simulations results

Evaluation of simulation results is performed using observation data by plotting dynamic graphs, scatter plots and by calculating statistical criteria. Observed data can be loaded by using the SticsRFiles get_obs() function and provided to CroPlotR functions.

# Insert your R code here!

# Path of the observation files is the same as this of simulations outputs

# SOLUTION
get_obs(workspace = workspace_path)

# Insert your R code here!

# Use the same syntax as for the get_sim function

# SOLUTION
get_obs(workspace = workspace_path, usm = c("wheat", "maize"), var = c("lai_n","masec_n"))

usms <- get_usms_list(file = file.path(workspace_path, "usms.xml"))
sim_done <- length(list.files(pattern = "^mod_s", workspace_path)) >=
  length(usms)

if (!sim_done) {
  print("Running simulations, please wait")
    run_javastics(javastics = javastics_path,
                  workspace = workspace_path, verbose = FALSE)
}
file_path <- file.path(tmp_dir,"sim_obs1.RData")
if (!file.exists(file_path)) {
  results2 <- get_sim(workspace = workspace_path, usm = c("wheat", "maize"), verbose = FALSE)
  observations2 <- get_obs(workspace = workspace_path, usm = c("wheat", "maize"), verbose = FALSE)
  save(list=c("results2","observations2"), file=file_path)
} else {
  #print("loading")
  load(file=file_path)
}

# Insert your R code here!

# See previous exercises for dynamics plotting

# Use the "obs" argument of the CroPlotR plot function

# SOLUTION
plot(results2, var = c("lai_n","masec_n"), obs = observations2)

# Insert your R code here!

# Try type = "scatter" ;-)

# SOLUTION
plot(results2, var = c("lai_n","masec_n"), obs = observations2, type = "scatter")

# Insert your R code here!

# Look at the shape_sit argument ...

# SOLUTION
plot(results2, var = c("lai_n","masec_n"), obs = observations2, type = "scatter", shape_sit = "symbol")

usms <- get_usms_list(file = file.path(workspace_path, "usms.xml"))
sim_done <- length(list.files(pattern = "^mod_s", workspace_path)) >=
  length(usms)

if (!sim_done) {
  print("Running simulations, please wait")
    run_javastics(javastics = javastics_path,
                  workspace = workspace_path, verbose = FALSE)
}
usms <- c("wheat","DurumWheat", "maize", "sorghum")
file_path <- file.path(tmp_dir,"sim_obs2.RData")
if (!file.exists(file_path)) {
  results4 <- get_sim(workspace = workspace_path, usm = usms, verbose = FALSE)
  observations4 <- get_obs(workspace = workspace_path, usm = usms, verbose = FALSE)
  save(list=c("results4","observations4"), file=file_path)
} else {
  #print("loading")
  load(file=file_path)
}

# Insert your R code here!

# Set shape_sit = "group" and use the situation_group argument containing the groups definition (as a list of list)

# SOLUTION
groups <- list(list("wheat","DurumWheat"), list("maize", "sorghum"))
plot(results4, var = c("lai_n","masec_n"), obs = observations4, type = "scatter", shape_sit = "group", situation_group = groups)

# Insert your R code here!

# same use as for the plot function ...

# SOLUTION
summary(results4, obs = observations4)

# Insert your R code here!

# use the argument all_situations ...

# SOLUTION
summary(results4, obs = observations4, all_situations = FALSE)

# Insert your R code here!

# Just use the object returned by the summary function in the plot function ...

# SOLUTION
summ <- summary(results4, obs = observations4, stats = c("R2","nRMSE"))

plot(summ, title = "for all USMs")

# Insert your R code here!

# use the argument all_situations in summary ...

# SOLUTION
summ <- summary(results4, obs = observations4, all_situations = FALSE, stats = c("R2","nRMSE"))

plot(summ, title = "for each USM")

# Insert your R code here!

# Use the xvar argument of plot ...

# SOLUTION
summ <- summary(results4, obs = observations4, all_situations = FALSE, stats = c("R2","nRMSE"))

plot(summ, xvar = "situation", title = "USM names as X")

Generating JavaStics input files

The JavaStics Input files (XML format) can be generated from R using functions of the SticsRFiles package.

Generating observation files

Observation files can be generated from a data.frame using the function gen_obs of the SticsRFiles package.

They are not used by the STICS model but can be used to plot comparison between observed and simulated variables and to optimize STICS parameters both in JavaStics and in CroPlotR and CroptimizR packages as we have seen before.

• Generate observation files from the data stored in the df data.frame using the gen_obs function and check their content using the get_obs function.

df <- data.frame(usm_name = c("usm1","usm1","usm2"), ian = c(2020,2021,2021), mo = c(12,7,7), jo = c(25,27,20),jul = c(359,573,566),`lai_n` = c(0.5,4.5,4.2), mafruit = c(NA,6.5,6.1))

# Insert your R code here!

# Use the out_path argument of get_obs if you want to specify the path of the folder where the obs files will be created.
# By default they are created in the current working directory of the R process (given by getwd()).

# Use the workspace argument of get_obs to specify the path to the folder containing the observation files.
# if you did not use the out_path argument in gen_obs, then use workspace = getwd() in get_obs.

#SOLUTION:
gen_obs(df = df)
get_obs(workspace = getwd())

The need of a wrapper

We saw earlier how to run a simulation using run_javastics(). This function simply calls JavaStics from the command line to run simulations. This does the same as running an USM by opening the JavaStics user interface, but from R.

You are probably familiar with the xml files used as inputs for JavaStics (e.g. sols.xml or usms.xml). These files define all the parameters values for one or several USMs, with potentially several soils, management practices and climate. Then when we run a given USM, JavaStics reads the parameters values from the xml files, and writes new text files with only the values associated to this particular USM. Then it calls the STICS model that reads those text files, runs the simulation, and writes the outputs.

So when we run a USM from JavaStics we have four different steps:

1. JavaStics reads the parameters values for the USM from the xml files
2. JavaStics writes text files with the parameters values that only describe this USM
3. JavaStics calls the STICS model (i.e the STICS executable)
4. the STICS model reads the parameters values from the text files, runs a simulation, and writes the output files (i.e. mod_bXXX.sti and mod_sXXX.sti)

This method is simpler for the user because it hides the technical details, but it has some downsides: it is slow because we have to create the text files each time we make a simulation, and we cannot parallelize the simulations because we share the same text files.

But no worries, SticsRPacks provides a solution to this!

Instead of creating the text files one USM at a time, we can create them for all USMs once into one folder for each, and then we directly call the STICS model (i.e without using JavaStics) on these files. Using this method we are now able to run the simulations way faster, and in parallel!

So how do we do that? Well, we use three steps:

• create the text files (gen_usms_xml2txt())
• set-up our modeling environment (stics_wrapper_options())
• run the simulations and import the results (stics_wrapper())

Creating the text files

The first step is to create the text files for each USMs in a separate directory. To do so, we will use the function gen_usms_xml2txt().

Before generating the text files, we will first restrain the number of variables we want STICS to output using gen_varmod.

• Use the gen_varmod function to tell STICS we only want to simulate “lai_n”, “masec_n” and “mafruit”

# Insert your R code here!

# Use `workspace_path = workspace_path`

#SOLUTION:
gen_varmod(workspace_path, var = c("lai_n", "masec_n", "mafruit"))

It is generally useful (although not mandatory) to pre-filter the variables to avoid filtering the variables for the wrapper, the plots and the statistics.

Now we can generate our text files. Remember, the path of the JavaStics distribution downloaded by this tutorial is stored in the javastics_path R object, and the path to the example workspace in workspace_path.

• Use the function gen_usms_xml2txt() to generate the text files for the USMs “wheat”, “banana” and “soybean” in sub-folders of the workspace_path folder:

# Insert your R code here!

# Use `javastics_path = javastics_path` and `workspace_path = workspace_path`
# No need to use the `target_path` argument here, it will be set by default to `workspace_path`

#SOLUTION:
gen_usms_xml2txt(javastics_path, workspace_path, usm = c("wheat","banana","soybean"))

The text files have been generated in subfolders of workspace_path, i.e. /tmp/SticsRPacks/JavaSTICS-1.5.0-STICS-10.0.0/example, one folder per USM, each taking the name of the corresponding USM.

# generate missing dirs, if any
dirs <- c("wheat","banana","soybean")
dirs_ex <- dir.exists(file.path(workspace_path,dirs))
if (!all(dirs_ex)) gen_usms_xml2txt(javastics_path, workspace_path, usm = dirs[!dirs_ex], verbose = FALSE)

Now, print the files inside the created folders:

# Insert your R code here!

#SOLUTION:
dirs <- file.path(workspace_path, c("wheat","banana","soybean"))
list.files(dirs, full.names = TRUE)

Setting-up the modeling environment

Before running the simulations using STICS, we have to define some information about our modeling environment. This is done using stics_wrapper_options(). This function, defined in SticsOnR, helps us define where is our installation of JavaStics, which STICS executable to use (e.g. if you are using an executable from a research branch), which workspace to use, etc…

• Use stics_wrapper_options() to set-up our modeling environment using a parallel computation over the USMs for which we generated text files here-before:

# Insert your R code here!

# Use `javastics = javastics_path` and `workspace = workspace_path`
# (remember, the folders including the text files for each USMs have been created in the folder workspace_path)

#SOLUTION:
stics_wrapper_options(javastics = javastics_path, workspace = workspace_path, parallel = TRUE)

Running the simulations

Now everything is set-up to make our simulations! The model options we will use here are stored in the wrapper_options object with the options we just set up in the previous section.

# generate missing dirs, if any
dirs <- c("wheat","banana","soybean")
dirs_ex <- dir.exists(file.path(workspace_path,dirs))
if (!all(dirs_ex)) gen_usms_xml2txt(javastics_path, workspace_path, usm = dirs[!dirs_ex], verbose = FALSE)
wrapper_options_file <- file.path(tmp_dir,"wrapper_options.RData")
if (!file.exists(wrapper_options_file)) {
  wrapper_options = stics_wrapper_options(javastics = javastics_path, workspace = workspace_path, parallel = TRUE, verbose = FALSE)
  save(list = "wrapper_options", file = wrapper_options_file)
} else {
  #print("loading")
  load(file = wrapper_options_file)
}

• Now use stics_wrapper() to simulate leaf area index and above-ground biomass for the USM “wheat” only (not “banana” and “soybean”) using the model options stored in wrapper_options, and print the results:

# Insert your R code here!

# The situation argument is used to choose the USMs to simulate.

# The var_names argument is used to choose the variables to return.

# Use get_var_info() to search for a variable name.

#SOLUTION:
res1 <- stics_wrapper(model_options = wrapper_options, situation = "wheat", var = c("lai_n", "masec_n"))
print(res1)

Data structure

The data structure returned by stics_wrapper is a list of two:

• error: TRUE if there was an error during the simulation, FALSE otherwise
• sim_list: a list of the simulation results

The content of sim_list follows the same format as the outputs returned by get_obs and get_sim: it is a list of data.frames. Each data.frame corresponds to the simulation results of a USM, with each row being a value for a given day, and each column a variable. The data can be accessed using indexing or the $ syntax, i.e.: res$sim_list$wheat would return a data.frame with the simulation results for the USM called “wheat”.

Parameter forcing

Sometimes it is useful to compare the outputs of the model with a different value for one or several parameters to better understand its impact or check if it is used.

stics_wrapper() provides an easy way to do so with the param_values argument. This argument is given as a named vector.

• Make a simulation on the “wheat” USM and output just the lai, but this time with dlaimax = 0.005 and durvieF = 80:

# Insert your R code here!

# `param_values` is given as a named vector, e.g. c("bdens" = 5, "hautmax" = 1.1)

#SOLUTION:
res2 <- stics_wrapper(model_options = wrapper_options, situation = "wheat", var = "lai_n", param_values = c("dlaimax" = 0.005, "durvieF" = 80))
print(res2)

Plot the comparison

Now we can plot the outputs of the two simulations using the CroPlotR plot() function. The results of both simulations are accessible in the res1 and res2 objects. Remember, the simulated variables values are stored in the sim_list field of each object. For example, to get the values of the simulated variables for res1, we would do: res1$sim_list.

# generate missing dirs, if any
dirs <- c("wheat","banana","soybean")
dirs_ex <- dir.exists(file.path(workspace_path,dirs))
if (!all(dirs_ex)) gen_usms_xml2txt(javastics_path, workspace_path, usm = dirs[!dirs_ex], verbose = FALSE)
if (!file.exists(file.path(tmp_dir,"comp_res.RData"))) {
  wrapper_options = stics_wrapper_options(javastics = javastics_path,
                                          workspace = workspace_path, parallel = TRUE,
                                          verbose = FALSE)
  res1 = stics_wrapper(model_options = wrapper_options, situation = "wheat",
                       var = c("lai_n", "masec_n"))
  res2 = stics_wrapper(model_options = wrapper_options, situation = "wheat",
                       var = "lai_n",
                       param_values = c("dlaimax" = 0.005, "durvieF" = 80))
  
  save(list=c("res1", "res2"), file=file.path(tmp_dir,"comp_res.RData"))
} else {
  #print("loading")
  load(file=file.path(tmp_dir,"comp_res.RData"))
}

Then, to plot a simulation result, we would do:

# Execute the code below to create a plot of the results:
plot(res1$sim_list)

It is also possible to compare the results of several simulations by giving them as arguments to the function one after the other.

• Plot the comparison between the outputs of both simulations:

# Insert your R code here!

# The simulation results can be listed one after the other in `plot()` for comparison

#SOLUTION:
plot(res1$sim_list,res2$sim_list)

We can also name the simulations in our plot by naming them when passed as arguments to plot().

• Plot the comparison between the outputs of both simulations, and name the first simulation “original” and the second simulation “dlaimax=0.005, durvieF=80”:

# Insert your R code here!

# Name the arguments of the plot function, e.g. plot("sim1" = res1) would name the simulation "sim1" (but would not output it on the plot because there is only one)

#SOLUTION:
plot("original" = res1$sim_list,"dlaimax=0.005, durvieF=80" = res2$sim_list)

Workflow

Model calibration is an important part of most modeling work, but it can be complex and cumbersome. SticsRPacks provides a nice, standardized and easy way to estimate the values of parameters using observations through the CroptimizR estim_param() function.

The workflow for parameter estimation is made using several steps:

• Modeling environment set-up (stics_wrapper_options()). estim_param() uses stics_wrapper() under the hood to run each simulation, so we must set-up our modeling environment first as we saw before.
• Import the observations for each USM (get_obs() + filter_obs()).
• Choose the parameter(s) to optimize and their boundary values.
• Set-up options for the optimization method,
• run the parameter optimization (estim_param).

# generate missing dirs, if any
dirs <- c("wheat","banana","soybean")
dirs_ex <- dir.exists(file.path(workspace_path,dirs))
if (!all(dirs_ex)) gen_usms_xml2txt(javastics_path, workspace_path, usm = dirs[!dirs_ex], verbose = FALSE)
if (!file.exists(file.path(tmp_dir,"store-param-estim.RData"))) {
  wrapper_options = stics_wrapper_options(javastics = javastics_path, workspace = workspace_path, parallel = TRUE, verbose = FALSE)
  obs = get_obs(workspace = workspace_path, usm = "wheat", verbose = FALSE)
  obs_filt = filter_obs(obs, "lai_n", include = TRUE)
  param_info = list(stamflax = list(lb = 200, ub = 400))
  optim_options = list(nb_rep = 3, maxeval = 10, xtol_rel = 1e-03, out_dir = workspace_path, ranseed = 1234)
  
  save(list=c("wrapper_options", "obs", "obs_filt", "param_info", "optim_options"),
       file=file.path(tmp_dir,"store-param-estim.RData"))
} else {
  #print("loading")
  load(file=file.path(tmp_dir,"store-param-estim.RData"))
}

Modelling environment

We will simply re-use the modeling environment from the previous exercise: wrapper_options. This is how we did it:

wrapper_options = stics_wrapper_options(javastics = javastics_path, workspace = workspace_path, parallel = TRUE)

Observations

• Import the observations for the “wheat” USM, and put it in an object called obs:

# Insert your R code here!

# Use `get_obs()`, workspace_path is the workspace.

#SOLUTION:
obs = get_obs(workspace = workspace_path, usm = "wheat")

• Filter the observations to get only the lai using filter_obs(), and put it in an object called obs_filt:

# Insert your R code here!

# Use the argument var_names in filter_obs ...

#SOLUTION:
obs_filt = filter_obs(obs, "lai_n", include = TRUE)

Parameters of the optimization method

The optimization methods available in CroptimizR have some parameters too. For example the simplex method has the following:

• nb_rep: Number of times the minimization process is done, each time starting with a different set of initial parameters values. This one is used to control whether our initial values impact our resulting optimized value or not. The more is better, but the whole optimization process is repeated nb_rep so it takes more and more time.
• maxeval: Maximum number of evaluations of the minimized criteria. That is to say the algorithm is stopped if it did not converge to a value after maxeval runs.
• xtol_rel: Tolerance criterion between two iterations, i.e. stop the algorithm when an optimization step (or an estimate of the optimum) changes every parameter by less than xtol_rel multiplied by the absolute value of the parameter.
• out_dir: path to the directory on which the results will be written (graph and .Rdata)
• ranseed: set the random seed so that each execution give the same results. If you want full randomization, don’t set it.

Only out_dir is mandatory, all the others have default values.

These parameters are passed to estim_param as a named list.

• Make a named list to set values for these parameters and put the results in an object named optim_options. We require 3 repetitions and maxeval=10 (to make it quick … in real life you should set maxeval to a large value to let the process converge). You can use any directory to store the results. If you don’t particularly care where to put them, you can put them in workspace_path.

# Insert your R code here!

# A name list is constructed like so: list(name1 = value1, name2 = value2)

#SOLUTION:
optim_options = list(nb_rep = 3, maxeval=10, out_dir = workspace_path)

Choose parameter

We have to choose which parameter will be optimized, and what are their respective boundary values. To do so, we define a list of list for each parameter like so:

param_info = list()
param_info$dlaimax = list(lb = 0.0005, ub = 0.0025)
param_info$durvieF = list(lb = 100, ub = 450)

Let’s try to optimize the value of the stamflax parameter.

• Create the parameter list for stamflax and name it param_info.

# Insert your R code here!

# You can use 200 and 400 for boundary values

#SOLUTION:
param_info = list(stamflax = list(lb = 200, ub = 400))

Optimization

Now that we defined all information we needed, we can use estim_param() to make the optimization.

As a reminder, here’s the objects we previously defined:

• obs_filt: the observations to use in the optimization

• stics_wrapper: this is the wrapper function for stics

• wrapper_options: the options passed to stics_wrapper

• optim_options: the parameters for the optimization method

• param_info: the parameter(s) names and boundary values

• Now estimate the value of stamflax using estim_param() and the previous objects we defined, and put the results in an object called optim_res.

# Insert your R code here!

# Read the documentation of estim_param() and associate the aforementioned objects to the corresponding arguments.

#SOLUTION:
optim_res = estim_param(obs_list = obs_filt,
                        model_function = stics_wrapper,
                        model_options = wrapper_options,
                        optim_options = optim_options,
                        param_info = param_info)

# generate missing dirs, if any
dirs <- c("wheat","banana","soybean")
dirs_ex <- dir.exists(file.path(workspace_path,dirs))
if (!all(dirs_ex)) gen_usms_xml2txt(javastics_path, workspace_path, usm = dirs[!dirs_ex], verbose = FALSE)
if (!file.exists(file.path(tmp_dir,"prepare-param-estim-out.RData"))) {
  load(file = system.file("extdata/optim_results.Rdata", package = "SticsRPacks"))
  optim_res = res
  rm(res)
  wrapper_options = stics_wrapper_options(javastics = javastics_path,
                                          workspace = workspace_path, parallel = TRUE,
                                          verbose = FALSE)
  obs_filt = filter_obs(get_obs(workspace = workspace_path, usm = "wheat",
                                verbose = FALSE), "lai_n", include = TRUE)
  sim_after_optim = stics_wrapper(param_values = optim_res$final_values,
                                  model_options = wrapper_options, situation = "wheat")$sim_list
  sim_before_optim = stics_wrapper(model_options = wrapper_options, situation = "wheat")$sim_list
  
  save(list=c("optim_res", "wrapper_options", "obs_filt", "sim_after_optim", "sim_before_optim"),
       file=file.path(tmp_dir,"prepare-param-estim-out.RData"))
} else {
  #print("loading")
  load(file=file.path(tmp_dir,"prepare-param-estim-out.RData"))
}

Plots

Now that our parameter value is optimized, we can re-run the model with this value to check the results. To do so, we call again stics_wrapper with the optimized value passed as argument.

• Use stics_wrapper to make a new simulation of the wheat USM with our optimized value optim_res$final_values, and store the simulated results (i.e. the content of sim_list) in an object called sim_after_optim.

# Insert your R code here!

# Read the documentation of stics_wrapper() and use the `param_values` argument along with the resulting optimized value stored in `optim_res$final_values`.

# Use the situation argument of stics_wrapper() to select the USM to simulate.

# stics_wrapper() returns a list of two:
# - sim_list for the results
# - error for indications on errors
# We need only the results in our object `sim_after_optim`.

#SOLUTION:
sim_after_optim = stics_wrapper(param_values = optim_res$final_values, model_options = wrapper_options, situation = "wheat")$sim_list

Now we can plot the output using plot().

• Use plot() to make a dynamic plot of the simulations along with the observations (obs_filt) for the lai_n variable.

# Insert your R code here!

# Read the documentation of plot(). This one is tricky, you can either see it on the package website <https://SticsRPacks.github.io/CroPlotR/reference/plot.cropr_simulation.html> or from R using ?plot.cropr_simulation

#SOLUTION:
plot(sim = sim_after_optim, obs = obs_filt, var = "lai_n")

We can also make a scatter plot for a different display of the fit between observations and simulations.

• Make the same plot but as a scatter plot instead of a dynamic plot

# Insert your R code here!

# Read the documentation of plot(). This one is tricky, you can either see it on the package website <https://SticsRPacks.github.io/CroPlotR/reference/plot.cropr_simulation.html> or from R using ?plot.cropr_simulation

#SOLUTION:
plot(sim = sim_after_optim, obs = obs_filt, type = "scatter", var = "lai_n")

Before/after comparison

It is often interesting to compare the simulations before and after the optimization to investigate if it worked properly. To do so, we just have to simulate both, and to plot the results using plot().

• Make a simulation using stics_wrapper with the original values of the parameters, and store the simulated results (i.e. the content of sim_list) in an object called sim_before_optim.

# Insert your R code here!

# estim_param() does not modify the input files, so you can just re-simulate the USM "wheat" with stics_wrapper

#SOLUTION:
sim_before_optim = stics_wrapper(model_options = wrapper_options, situation = "wheat")$sim_list

Now we can compare the simulations using plot() and summary().

• Make a plot to compare both simulations for lai_n

# Insert your R code here!

# Read the documentation of plot(). This one is trick, you can either see it on the package website <https://SticsRPacks.github.io/CroPlotR/reference/plot.cropr_simulation.html> or from R using ?plot.cropr_simulation

#SOLUTION:
plot(before_optim=sim_before_optim, after_optim=sim_after_optim, obs= obs_filt, var = "lai_n")

• Compute the nRMSE and Bias of both simulations for lai_n using summary()

# Insert your R code here!

# Read the documentation of summary(). This one is trick, you can either see it on the package website <https://SticsRPacks.github.io/CroPlotR/reference/summary.cropr_simulation.html> or from R using ?summary.cropr_simulation

# use the "stat" argument of summary to select the criteria to compute ...

#SOLUTION:
summary(before_optim=sim_before_optim, after_optim=sim_after_optim, obs = obs_filt, stats=c("nRMSE", "Bias"))

Simulating successive crops

The run_javastics function CAN NOT be used to run successive USMs. In this case, it is advised to use the stics_wrapper() function.

• Have a look at the argument successive_usms of the function stics_wrapper_options(). Then, simulate USMs demo_Wheat1, demo_BareSoil2and demo_maize3 (defined in the JavaStics distribution example workspace) both in independent and successive ways using stics_wrapper. Check that the last values of the simulated mafruit for demo_maize3 differ between the two different simulation modes (independent or successive USMs).

# Insert your R code here!

# First generate the txt files using gen_usms_xml2txt

# successive_usms argument must be a list of vectors of successive USMs,
# thus in this case list(c("demo_Wheat1","demo_BareSoil2","demo_maize3"))

# Use the tail function with argument n=1 to extract the last element of a vector.

#SOLUTION:

# Generate the txt files
gen_usms_xml2txt(javastics_path, workspace_path, usm = c("demo_Wheat1","demo_BareSoil2","demo_maize3"))

# Run without succession
wrapper_options <- stics_wrapper_options(javastics = javastics_path,
                                         workspace = workspace_path)
without_succ <- stics_wrapper(model_options = wrapper_options,
                              situation = c("demo_Wheat1","demo_BareSoil2","demo_maize3"))

# Run with succession
wrapper_options <- stics_wrapper_options(javastics = javastics_path,
                                         workspace = workspace_path,
                                         successive = list(c("demo_Wheat1","demo_BareSoil2","demo_maize3")))
with_succ <- stics_wrapper(model_options = wrapper_options,
                           situation = c("demo_Wheat1","demo_BareSoil2","demo_maize3"))

# Compare both
print(paste("mafruit demo_maize3 without succession:",
            tail(without_succ$sim_list$demo_maize3$mafruit, n=1)))
print(paste("mafruit demo_maize3 with succession:",
            tail(with_succ$sim_list$demo_maize3$mafruit, n=1)))

Plotting successive simulations

Dynamic plots of variables can be aggregated to generate continuous plots.

• Have a look at the argument successive of the plot function of the CroPlotR package. Then plot the azomes variable (amount of NO3-N in soil) simulated in the previous exercise for both modes (use the without_succ and with_succ R variables that contains the results returned by the wrapper, as presented in the SOLUTION of the previous exercise). Check that azomes is more continuous in the second case.

if (!file.exists(file.path(tmp_dir,"prepare-plot-successive.RData"))) {
  gen_usms_xml2txt(javastics_path, workspace_path,
                   usm = c("demo_Wheat1","demo_BareSoil2","demo_maize3"),
                   verbose = FALSE)
  
  # Run without succession
  wrapper_options <- stics_wrapper_options(javastics = javastics_path,
                                           workspace = workspace_path,
                                           verbose = FALSE)
  without_succ <- stics_wrapper(model_options = wrapper_options,
                                situation = c("demo_Wheat1","demo_BareSoil2","demo_maize3"))
  
  # Run with succession
  wrapper_options <- stics_wrapper_options(javastics = javastics_path,
                                           workspace = workspace_path,
                                           successive = list(c("demo_Wheat1","demo_BareSoil2","demo_maize3")),
                                           verbose = FALSE)
  with_succ <- stics_wrapper(model_options = wrapper_options,
                             situation = c("demo_Wheat1","demo_BareSoil2","demo_maize3"))
  
  save(list=c("without_succ", "with_succ"),
       file=file.path(tmp_dir,"prepare-plot-successive.RData"))
} else {
  #print("loading")
  load(file=file.path(tmp_dir,"prepare-plot-successive.RData"))
}

# Insert your R code here!

#SOLUTION:
plot("Without successive"=without_succ$sim_list,
     "With successive"=with_succ$sim_list,
     var = "azomes",
     successive = list(list("demo_Wheat1","demo_BareSoil2","demo_maize3")))

Estimating parameters on successive simulations

Estimating parameters on successive simulations is also possible. Actually it may not require specific options in the estim_param function: if the successive_usms argument is defined in the STICS wrapper options, then the parameters will be estimated on successive simulations.

Note however that by default the same values of the parameters will be applied on all USMs. If one want to estimate the parameters for specific USMs only, or with different values depending on the USMs, one must use the sit_list element of the param_info argument.