STICS-Pea

Appropriation of the STICS model

Abdou Diallo

2024-12-19

Setting

Packages

Code 
library(SticsRPacks)
library(readxl) # pour les fichier excel
library(DT) # pour datatable(), view df
library(hydroGOF) # pour calculer la RMSZ
library(htmltools) # pour les tag html
library(compareDF) # pour comparer graphiquement des df
library(dplyr) 
library(tidyr)
library(ggplot2)

Functions

Code 
#Getting information about parameters
get_param_info()
get_param_info(param = c("alb", "hum"))
get_param_info(keyword = "STATION")

Getting parameters values from XML files

Code 
#For one parameter and one occurrence
# Fixing files paths
sols <- file.path(workspace, "sols.xml")
par_gen <- file.path(workspace, "param_gen.xml")
station <- file.path(workspace, "climatex_sta.xml")

# A option parameter
get_param_xml(par_gen, param = "codeactimulch")

# A simple parameter
get_param_xml(station, param = "concrr")

# Using a conditional selection
get_param_xml(sols, param = "argi", select = "sol", select_value = "solcanne")

Getting variables informations

Code 
 # get_var_info(var = "msrac_n")
 # get_var_info(var = "masec_n")
 # get_var_info(var = "totir")
 # get_var_info(var = "totapN")
 # get_var_info(var = c("inn", "inn1moy","inn2moy"))
 # get_var_info(var = "totapN")
# get_var_info(var = "turfac")
# get_var_info(var = "inn")
# get_var_info(var = "leaching_from_lev")
# get_var_info(var = "N_volatilisation")

Introduction

Simulation of the actual yields of the pea crops:

  1. Spring Pea
  2. Winter Pea

Variables to check

  • masec_n: Biomass of aboveground plant (t.ha-1)
  • mafruit: Biomass of harvested organs (t.ha-1)
  • totir: Cumulative amount of irrigation water (mm)
  • totapN: Cumulative amount of nitrogen (kg.ha-1)
  • turfac: turgescence water stress index
  • inn: nitrogen nutrition index (NNI)

1 Data preparation

1.1 Data import

See the input Excel file : _inputs_stics_pea.xlsx.

Code 
# Getting the file in a specic directory
# usm_xl_file = download_usm_xl(out_dir = paste0(r_path,"/data"))

usm_xl_file = file.path(workspace, input_file)

1.2 USM file

Getting data

Code 
usms_param <- read_params_table(usm_xl_file, sheet_name = "USMs", 
                                num_na = "",
                                char_na = "")
tags$div(
  tags$h3(""),
  DT::datatable(usms_param, options = list(pageLength = 5))
)

Genarating the usms.xml file

Code 
# Generating a new usms.xml file, for all xl_param lines
out_file <- file.path(workspace, "usms.xml")
gen_usms_xml(file = out_file, param_df = usms_param)

The usms.xml file content

Focus of the first usm parameters:

## <?xml version="1.0" encoding="UTF-8" standalone="no"?>
## <usms version="10.1.0">
##   <usm nom="baccara_PP">
##     <datedebut>61</datedebut>
##     <datefin>201</datefin>
##     <finit>ACBB_Mons_T5_ini.xml</finit>
##     <nomsol>ACBB_Mons_T5</nomsol>
##     <fstation>ACBB_Mons_sta.xml</fstation>
##     <fclim1>ACBB_Mons.2010</fclim1>
##     <fclim2>ACBB_Mons.2010</fclim2>
##     <culturean>1</culturean>
##     <nbplantes>1</nbplantes>
##     <codesimul>0</codesimul>
##     <plante dominance="1">
##       <fplt>pea_v10_plt.xml</fplt>
##       <ftec>baccara_PP_tec.xml</ftec>
##       <flai>null</flai>
##       <fobs>03_T5_2010_b1.obs</fobs>
##     </plante>
##     <plante dominance="2">
##       <fplt>null</fplt>
##       <ftec>null</ftec>
##       <flai>null</flai>
##       <fobs>null</fobs>
## ...
## </usms>

1.3 Initialization

Getting data

Code 
ini_param <- read_params_table(usm_xl_file, sheet = "Ini")

tags$div(
  tags$h3(""),
  DT::datatable(ini_param, options = list(pageLength = 5))
)

Generating multiple initialization files (*_ini.xml)

Code 
# *_ini.xml files, one for each xl_param line
gen_ini_xml(param_df = ini_param, out_dir = workspace)

Focus of a ini parameters file content sub-list:

## <?xml version="1.0" encoding="UTF-8"?>
## <initialisations version="10.1.0">
##   <nbplantes>1</nbplantes>
##   <plante dominance="1">
##     <stade0>snu</stade0>
##     <lai0>0</lai0>
##     <magrain0>0</magrain0>
##     <zrac0>0</zrac0>
##     <option choix="2" nom="Simulation of Nitrogen and Carbon reserves" nomParam="code_acti_reserve">
##       <choix code="1" nom="yes">
##         <maperenne0>0</maperenne0>
##         <QNperenne0>0</QNperenne0>
##         <masecnp0>0</masecnp0>
##         <QNplantenp0>0</QNplantenp0>
##       </choix>
##       <choix code="2" nom="no">
##         <masec0>0</masec0>
##         <QNplante0>0</QNplante0>
##         <restemp0>0</restemp0>
## ...
## </initialisations>

1.4 Soils

Getting data

Code 
soils_param <- read_params_table(usm_xl_file, sheet = "Soils")

tags$div(
  tags$h3(""),
  DT::datatable(soils_param, options = list(pageLength = 5))
)

Genarating the soils.xml file

Code 
out_file <- file.path(workspace, "sols.xml")
gen_sols_xml(file = out_file, param_df = soils_param)
Code 
# Verification
path_ref =  "/home/ardiallo/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0/stics_mod/example_pw"
df_ref = get_param_xml(file.path(path_ref, "sols.xml"))[[1]]
df = get_param_xml(file.path(workspace, "sols.xml"))[[1]]
df_ref = data.frame(df_ref)
df = data.frame(df)
all(names(df)==names(df_ref))

The sols.xml file content:

Focus of the first sol parameters sub-list:

## <?xml version="1.0" encoding="UTF-8"?>
## <sols version="10.1.0">
##   <sol nom="ACBB_Mons_T5">
##     <param format="real" max="60.0" min="0.0" nom="argi">17</param>
##     <param format="real" max="0.5" min="0.05" nom="norg">0.105</param>
##     <param format="real" max="60.0" min="10.0" nom="profhum">35</param>
##     <param format="real" max="100.0" min="0.0" nom="calc">0.53</param>
##     <param format="real" max="9.0" min="4.0" nom="pH">8.25</param>
##     <param format="real" max="0.5" min="0.0" nom="concseuil">0</param>
##     <param format="real" max="0.6" min="0.05" nom="albedo">0.22</param>
##     <param format="real" max="50.0" min="0.0" nom="q0">9.36</param>
##     <param format="real" max="1.0" min="0.0" nom="ruisolnu">0</param>
##     <param format="real" max="1000.0" min="10.0" nom="obstarac">210</param>
##     <param format="real" max="100.0" min="5.0" nom="pluiebat">50</param>
##     <param format="real" max="2.0" min="0.0" nom="mulchbat">1.5</param>
##     <param format="real" max="150.0" min="10.0" nom="zesx">60</param>
##     <param format="real" max="10.0" min="-10.0" nom="cfes">4</param>
##     <param format="real" max="0.2" min="0.01" nom="z0solnu">0.01</param>
##     <param format="real" max="20.0" min="8.0" nom="CsurNsol">9.46</param>
##     <param nom="finert">0.65</param>
##     <param format="real" max="5.0" min="0.0" nom="penterui">0.33</param>
##     <option choix="2" nom="pebbles" nomParam="codecailloux">
##       <choix code="1" nom="yes"/>
##       <choix code="2" nom="no"/>
##     </option>
##     <option choix="2" nom="macroporosity" nomParam="codemacropor">
##       <choix code="1" nom="yes"/>
##       <choix code="2" nom="no"/>
##     </option>
##     <option choix="2" nom="cracks (case of swelling clay soils)" nomParam="codefente">
## ...
##     </tableau>
## ...
## </sols>

1.5 Technical file

Getting data

Code 
tec_param <- read_params_table(usm_xl_file, sheet = "Tec")
tags$div(
  tags$h3(""),
  DT::datatable(tec_param, options = list(pageLength = 5))
)

Generating the crop management files (*_tec.xml)

Code 
# *_tec.xml files, one for each xl_param line
gen_tec_xml(param_df = tec_param, out_dir = workspace)

Focus of a tec parameters file content sub-list:

## <?xml version="1.0" encoding="UTF-8" standalone="no"?>
## <fichiertec version="10.1.0">
##   <formalisme nom="supply of organic residus">
##     <ta nb_interventions="0" nom="interventions">
##       <ta_entete nb_colonnes="7">
##         <colonne nom="julres"/>
##         <colonne nom="coderes"/>
##         <colonne nom="qres"/>
##         <colonne nom="Crespc"/>
##         <colonne nom="CsurNres"/>
##         <colonne nom="Nminres"/>
##         <colonne nom="eaures"/>
##       </ta_entete>
##     </ta>
##   </formalisme>
##   <formalisme nom="soil tillage">
##     <option choix="1" nom="Automatic calculation of the depth of residues incorporation in function of proftrav" nomParam="code_auto_profres">
##       <choix code="1" nom="yes">
##         <param format="real" max="0.0" min="1.0" nom="resk">0.14</param>
##         <param format="real" max="0.0" min="10.0" nom="resz">5</param>
##       </choix>
##       <choix code="2" nom="no"/>
##     </option>
##     <ta nb_interventions="1" nom="interventions">
##       <ta_entete nb_colonnes="3">
##         <colonne nom="jultrav"/>
##         <colonne nom="profres"/>
##         <colonne nom="proftrav"/>
##       </ta_entete>
##       <intervention nb_colonnes="3">
## ...
## </fichiertec>

** Cultivars**

Code 
# Read from a plant file (all cultivars available in a plant file)
get_cultivars_list(file = file.path(workspace,"plant", "pea_plt.xml"))
## [1] "baccara"

df = get_cultivars_param(file = file.path(workspace,"plant", "pea_plt.xml"))
tags$div(
  tags$h3("Cultivars list"),
  DT::datatable(df, options = list(pageLength = 5))
)

Cultivars list

1.6 Meteo station

Getting data

Code 
sta_param <- read_params_table(usm_xl_file, sheet = "Station")
tags$div(
  tags$h3(""),
  DT::datatable(sta_param, options = list(pageLength = 5))
)

Generating multiple meteorological station XML files (*_sta.xml)

Code 
# *_sta.xml files, one for each xl_param line
gen_sta_xml(param_df = sta_param, out_dir = workspace)
Code 
# Verification
path_ref =  "/home/ardiallo/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0/stics_mod/example_pw"
df_ref = get_param_xml(file.path(path_ref, "SUCS_sta.xml"))[[1]]
df = get_param_xml(file.path(workspace, "ACBB_Mons_sta.xml"))[[1]]
df_ref = data.frame(df_ref)
df = data.frame(df)
all(names(df)==names(df_ref))


get_param_xml(file.path(path_ref, "SUCS_sta.xml"), param = "concrr")
get_param_xml(file.path(workspace, "ACBB_Mons_sta.xml"), param = "concrr")

Focus of a sta parameters file content sub-list:

## <?xml version="1.0" encoding="UTF-8" standalone="no"?>
## <fichiersta version="10.1.0">
##   <formalisme nom="Weather station">
##     <param format="real" max="10.0" min="2.0" nom="zr">2.5</param>
##     <param format="real" max="10.0" min="0.0" nom="NH3ref">0</param>
##     <param format="real" max="3.0" min="0.0" nom="concrr">0.02</param>
##     <param format="real" max="90.0" min="-90.0" nom="latitude">49</param>
##     <param format="real" max="1200.0" min="800.0" nom="patm">1000</param>
##     <param format="real" max="25.0" min="4.0" nom="aclim">12.6115643657415</param>
##   </formalisme>
##   <formalisme nom="climate">
##     <option choix="1" nom="reading OR calculation of PET" nomParam="codeetp">
##       <choix code="1" nom="PET-Penman_reading"/>
##       <choix code="2" nom="PET-Penman_calculation"/>
##       <choix code="3" nom="PET-Shuttleworth-Wallace_calculation"/>
##       <choix code="4" nom="PET-Priestley-Taylor_calculation">
##         <param format="real" max="2.0" min="1.0" nom="alphapt">1.26</param>
##       </choix>
##     </option>
##     <option choix="1" nom="climate change" nomParam="codeclichange">
##       <choix code="1" nom="no"/>
##       <choix code="2" nom="yes"/>
##     </option>
##     <option choix="1" nom="climate in altitude" nomParam="codaltitude">
##       <choix code="1" nom="no"/>
##       <choix code="2" nom="yes">
##         <param format="real" max="2000.0" min="0.0" nom="altistation">0</param>
##         <param format="real" max="2000.0" min="0.0" nom="altisimul">0</param>
##         <param format="real" max="3.0" min="0.1" nom="gradtn">-0.5</param>
##         <param format="real" max="3.0" min="0.1" nom="gradtx">-0.55</param>
##         <param format="real" max="2000.0" min="0.0" nom="altinversion">500</param>
##         <param format="real" max="3.0" min="0.1" nom="gradtninv">1.3</param>
##         <param format="real" max="1.0" min="0.0" nom="cielclair">0.8</param>
##         <option choix="1" nom="option.adret.or.ubac" nomParam="codadret">
##           <choix code="1" nom="adret(south)"/>
##           <choix code="2" nom="ubac(north)">
##             <param format="real" max="5.0" min="-5.0" nom="ombragetx">-1.4</param>
##           </choix>
##         </option>
##       </choix>
##     </option>
## ...
## </fichiersta>

1.7 Weather file

NOT OK : To be done manually on JavaSTICS

Getting data

Code 
xl_param <- read_params_table(usm_xl_file, sheet = "Weather")
tags$div(
  tags$h3(""),
  DT::datatable(xl_param, options = list(pageLength = 5))
)

1.8 Observations

Getting data

Code 
obs_param <- read_params_table(usm_xl_file, sheet = "Obs")
tags$div(
  tags$h3(""),
  DT::datatable(obs_param, options = list(pageLength = 5))
)

Generating the observations files (*_obs.xml)

Code 
if (nrow(obs_param)>0){ 
  gen_obs(df = obs_param, out_dir = workspace)
  }

1.9 Generating general parameters XML files

  • param_gen.xml;
  • param_newform.xml
Code 
# with the latest version
# gen_varmod(workspace,var = output)
# gen_general_param_xml(out_dir = workspace,overwrite = TRUE)

# or with a specific version
# gen_general_param_xml(out_dir = workspace_path, stics_version = "V10.2.0")


# get_usms_list(file = file.path(workspace, "/usms.xml"))

1.10 Generate txt files

Code 
# Generate the model input files (text files) from javaStics input files (XML files)
gen_usms_xml2txt(javastics = javastics,
                 out_dir = txt_dir,
                 workspace = workspace)

2 Simulations

2.1 Run simulations

Code 
## get the list of USMs
usm_list <- get_usms_list(file = file.path(workspace,"usms.xml"))

#### run the USMs
 ## define the model options
model_options <- stics_wrapper_options(javastics = javastics, 
                                        workspace = txt_dir,
                                        parallel = TRUE)


 ## run the USMs as they are defined
res <- stics_wrapper(model_options = model_options, var=output)

res_df = CroPlotR::bind_rows(res$sim_list)
 # names(res_df)
 # head(res_df)
 # tail(res_df)
 #

2.2 Results

Code 
plot(ori=res$sim_list,var=c("mafruit","masec_n","inn","turfac"))
## $baccara_PH

## 
## $baccara_PP

Code 

# plot(ori=res$sim_list,var=c("inn1moy","turfac1moy"))

# p <- plot(res$sim_list, overlap = list(list("mafruit","masec_n"),
#                                       list("turfac1moy","inn1moy")),
#           var = c("mafruit","masec_n", "inn1moy","turfac1moy"))
# p


# plot(res$sim_list,successive = list(usm_list))
# ggsave(file.path(eval_dir,"dyn.png"), width=30, units="cm")

2.3 Simulated vs Observed

A suivre !

Conclusions