Simulation of Potential and Actual Yields

Abdou Diallo

2024-11-28

Load the necessary packages

Code

# Load the necessary packages
library(SticsRPacks)
library(tidyr)
library(dplyr)
library(ggplot2)

# Initialisation
javastics <- "/home/ardiallo/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0"
out_dir <- "/home/ardiallo/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0/STICS_homework"
workspace <- "/home/ardiallo/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0/STICS_homework"
Sys.setenv(javastics_path = javastics)

Introduction

Simulate the potential and actual yields of the three possible crops for the described area.

Crops

Spring barley
Spring durum wheat
Spring field pea

Variables to Check

masec_n:
mafruit: Biomass of harvested organs (t/ha)
totir: Cumulative amount of irrigation water (mm)
totapN: Cumulative amount of nitrogen (kg/ha)
turfac:
inns:

Data Preparation

Code

# get_usms_list(file = file.path(workspace, "/usms.xml"))
usm <- c("SpringBarley", "SpringFieldpea", "SpringDurumWheat")

Initial Parameters

Code

# Chargement des paramètres initiaux
file <- file.path(paste0(workspace, "/usms.xml"))
# get_param_xml(file, c("julres", "jultrav", "profres", "proftrav"))

Code

# D'abord modifier les option de simulation pour la question 1 puis remetre les options en 'yes' pour les autres questions
# "Model Inputs", "Global parameters", "General parameters", "Simulation options" and desactivate "Nitrogen stress activation" and "Water stress activation"


#A basic USM for each crop is available (spring barley, spring wheat and spring field pea).
file <- file.path(paste0(workspace, "/usms.xml"))
# get_param_xml(file)

The tec file

Use the standard cropping practices you are provided with (Excel file) to create the .tec file of each crop by modifying the base technical file (ORIGINAL_tec). Months are given for dates of practices, choose the days as you wish.

Code

### SpringBarley_tec.xml
param_df <- data.frame(tec_name="SpringBarley_tec.xml",
                       julres=214,#Organic inputs supply
                       jultrav_1 = 253, profres_1=10, proftrav_1=15, #Soil preparation
                       jultrav_2 = 375, profres_2=0, proftrav_2=10,# Second soil prep
                       iplt0 = 397, profsem = 2.5, densitesem = 350, variete = 1, #Sowing
                       julapN_or_sum_upvt = 434, 'absolute_value/%' = 60, engrais = 5) #fert

param_df <-rename(param_df, 'absolute_value/%' = "absolute_value..")

### Modify the tec file using the data.frame
gen_tec_xml(param_df = param_df, 
            file = file.path(workspace, "SpringBarley_tec.xml"), 
            out_dir = workspace)


### SpringFieldpea_tec.xml
param_df <- data.frame(tec_name="SpringFieldpea_tec.xml",
                       julres=214,#Organic inputs supply
                       jultrav_1 = 253, profres_1=10, proftrav_1=15, #Soil preparation
                       jultrav_2 = 375, profres_2=0, proftrav_2=10,# Second soil prep
                       iplt0 = 406, profsem = 3, densitesem = 100, variete = 1, #Sowing, 335?
                       julapN_or_sum_upvt = NA, 'absolute_value/%' = NA, engrais = NA) #fert


param_df <-rename(param_df, 'absolute_value/%' = "absolute_value..")

### Modify the tec file using the data.frame
gen_tec_xml(param_df = param_df, 
            file = file.path(workspace, "SpringFieldpea_tec.xml"), 
            out_dir = workspace)



### SpringDurumWheat
param_df <- data.frame(tec_name="SpringDurumWheat_tec.xml",
                       julres=214,#Organic inputs supply
                       jultrav_1 = 253, profres_1=10, proftrav_1=15, #Soil preparation
                       jultrav_2 = 0, profres_2=0, proftrav_2=0,# Second soil prep
                       iplt0 = 344, profsem = 2.5, densitesem = 350, variete = 1, #Sowing
                       julapN_or_sum_upvt = 434, 'absolute_value/%' = 90, engrais = 5) #fert

param_df <-rename(param_df, 'absolute_value/%' = "absolute_value..")

### Modify the tec file using the data.frame
gen_tec_xml(param_df = param_df, 
            file = file.path(workspace, "SpringDurumWheat_tec.xml"), 
            out_dir = workspace)

The climate file

Choose a year among 2011 to 2020 as you wish and define the climate files for the USM. Since all crops are spring crops requiring a winter soil preparation the preceding year, you must choose two successive years for your simulations (e.g. for a crop harvested in 2011, choose clima_RAIMAT2010 for the climate file of the first year, and clima_RAIMAT2011 for the climate file of the second year). Remember to report the years you have used when presenting your results.

Changing climate stations to 2010 and 2011

Code

file <- file.path(paste0(workspace, "/usms.xml"))

aa <- get_param_xml(file)
# aa$usms.xml$fclim1 
# aa$usms.xml$fclim2 


set_param_xml(file = file.path(workspace, "usms.xml"),
              param = c("fclim1"), values= "clima_RAIMAT.2010", overwrite = T)
set_param_xml(file = file.path(workspace, "usms.xml"),
              param = c("fclim2"), values= "clima_RAIMAT.2011", overwrite = T)

The txt file

Code

# Generate the model input files (text files) from javaStics input files (XML files)
gen_usms_xml2txt(javastics = javastics,
                 out_dir = out_dir,
                 workspace = workspace,
                 usm=c("SpringBarley","SpringFieldpea","SpringDurumWheat"))

Question 1: Simulate the potential yield of the three crops

Simulation results

Code

chemin = "~/Bureau/PhD_l/STICS/JavaSTICS-1.5.2-STICS-10.2.0/STICS_homework/script_et_info/"
res = readRDS(file = paste0(chemin,"res.rds"))

sim_all = CroPlotR::bind_rows(res$sim_list)
# names(sim_all)

# get_var_info(var = "masec_n")
# get_var_info(var = "mafruit")
# get_var_info(var = "totir")


df = get_var_info(var = names(sim_all) )
# df[names(sim_all) %in% df$name, c("name","definition")]

plot(ori=res$sim_list,
     var=c("mafruit","masec_n"))

$SpringBarley


$SpringDurumWheat


$SpringFieldpea

Code

# plot(ori=res$sim_list,
#      var=c("inns","turfac"))

# p <- plot(res$sim_list, overlap = list(list("mafruit","masec_n")))
# print(p)

Maximum potentiel yield 2010 - 2011

Code

#The maxi
df = sim_all %>%
  group_by(situation) %>%
  summarise(max_masec_n = round(max(masec_n),2),
            max_mafruit = round(max(mafruit),2)) %>%
  left_join(
    sim_all %>%
      group_by(situation) %>%
      filter(masec_n == max(masec_n) |  mafruit == max(mafruit)) %>%
      summarise(Date = first(Date)),
    by = "situation"
  )
df

# A tibble: 3 × 4
  situation        max_masec_n max_mafruit Date               
  <chr>                  <dbl>       <dbl> <dttm>             
1 SpringBarley            21.5        8.51 2011-06-15 00:00:00
2 SpringDurumWheat        21.1        6.07 2011-06-16 00:00:00
3 SpringFieldpea          11          5.36 2011-05-28 00:00:00

Code

# A tibble: 3 × 4

# situation
# <chr>
# max_masec_n
# <dbl>
# max_mafruit
# <dbl>
# Date
# <S3: POSIXct>
# SpringBarley  21.50   8.51    2011-06-15  
# SpringDurumWheat  21.06   6.07    2011-06-16  
# SpringFieldpea    11.00   5.36    2011-05-28  

# head(sim_all)
# tail(sim_all)

Code

# At 2011-06-30
# df = filter(sim_all, Date == as.POSIXct("2011-06-30", tz = "UTC"))
# df$Date = as.POSIXct("2011-06-30", tz = "UTC")
# df[, c("situation","masec_n","mafruit", "Date")]
# result = get_sim(workspace = workspace)
# p <- plot(result)
# print(p)

# MODIFICATION DES OPTION EN YES

Question 2: Actual Yields Over 10 Years

Years: 2010 to 2020

Code

output_2 =c("masec_n","mafruit","totir","totapN","turfac","inn", "turfac1moy","turfac2moy",
            "inn1moy","inn2moy","leaching_from_lev","leaching_from_plt","N_volatilisation")

# Simulate the actual yields obtained over 10 years, analyse the potential causes of the variability,
# in particular looking at yield components and environmental stresses; you can focus on one or two species out of the three.

file <- file.path(paste0(workspace, "/usms.xml"))

 # get_param_xml(file)#I nned to know the second day of simulation

 
 
 file <- file.path(paste0(workspace, "/usms.xml"))
 aa <- get_param_xml(file)
 # aa$usms.xml$fclim1 
 # aa$usms.xml$fclim2 
 
 set_param_xml(file = file.path(workspace, "usms.xml"),
               param = c("fclim1"), values= "clima_RAIMAT.2010", overwrite = T)
 set_param_xml(file = file.path(workspace, "usms.xml"),
               param = c("fclim2"), values= "clima_RAIMAT.2020", overwrite = T)
 
 
 ### SpringBarley_tec.xml
param_df <- data.frame(tec_name="SpringBarley_tec.xml",
                       julres=214,#Organic inputs supply
                       jultrav_1 = 253, profres_1=10, proftrav_1=15, #Soil preparation
                       jultrav_2 = 375, profres_2=0, proftrav_2=10,# Second soil prep
                       iplt0 = 397, profsem = 2.5, densitesem = 350, variete = 1, #Sowing
                       julapN_or_sum_upvt = 434, 'absolute_value/%' = 60, engrais = 5) #fert

param_df <-rename(param_df, 'absolute_value/%' = "absolute_value..")

### Modify the tec file using the data.frame
gen_tec_xml(param_df = param_df, 
            file = file.path(workspace, "SpringBarley_tec.xml"), 
            out_dir = workspace)
 
 # Generate the model input files (text files) from javaStics input files (XML files)
 gen_usms_xml2txt(javastics = javastics,
                  out_dir = out_dir,
                  workspace = workspace,
                  usm=usm)
 
 
 #### run the USMs
 ## define the model options
 model_options <- stics_wrapper_options(javastics = javastics, 
                                        workspace = workspace,
                                        parallel = TRUE)
 
 ## run the USMs as they are defined
res_2 <- stics_wrapper(model_options = model_options, var=output_2, 
                       situation = usm)
 
# plot(ori=res_2$sim_list,
#      var=c("mafruit","masec_n"))
# 
# plot(ori=res_2$sim_list,
#      var=c("inns","turfac"))

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

Code

sim_all_2 = CroPlotR::bind_rows(res_2$sim_list)
 # names(sim_all_2)
 # head(sim_all_2)
 # tail(sim_all_2)
 #

For SpringBarley

Code

### For SpringBarley
 SpringBarley_yield <- result_yield %>% filter(situation == "SpringBarley")
 SpringBarley_stress <- result_stress %>% filter(situation == "SpringBarley")



 # combinaison des donnée
 SpringBarley_combined <- bind_rows(
   SpringBarley_yield %>% select(years, masec_n, mafruit) %>% pivot_longer(cols = -years, names_to = "variable", values_to = "value"),
   SpringBarley_stress %>% select(years, turfac1moy, inn1moy) %>% pivot_longer(cols = -years, names_to = "variable", values_to = "value")
 )
 
 ordonnee =  c("masec_n", "mafruit", "inn", "inn1moy", "inn2moy","turfac","turfac1moy","turfac2moy")
 SpringBarley_combined$variable <- factor(SpringBarley_combined$variable, levels = ordonnee)
 ggplot(SpringBarley_combined, aes(x = years, y = value, color = variable, group = variable)) +
   geom_line(size = 1) +
   geom_point(size = 2) +
   facet_wrap(~ variable, scales = "free_y") + # Une facette par variable
   labs(
     title = "SpringBarley",
     x = "Years",
     y = "Value",
     color = "Variable"
   ) +
   theme_bw() +
   theme(
     plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
     legend.position = "bottom"
   )

Code

 #+scale_y_continuous(limits = c(0, NA))

The mean average yield between 2011 to 2020 = 2.41 ± 0.4

Question 3: Yield response to irrigation rate

SpringBarley: Years 2010 - 2020

The scenario 1 is better: doseI_1=60, julapI_1=510

Question 4: Yield response to irrigation and N fertiliser rate

SpringBarley: Years 2010 - 2020

Optimal irrigation and fertilisation rates increased the yield but also increased N leaching and volatilisation (see fig below)

Conclusions

Fertilisation = doseI_1=60, julapI_1=510

Irrigation = doseN_1=350, julapN_1=435