Azodyn-Pea 2019

Appropriation of the Azodyn-Pea 2019 model

Abdou Diallo

2024-12-19

Setting

Packages

Code

library(rvle)
library(hydroGOF) # pour calculer la RMSZ
library(htmltools) # pour les tag html
library(compareDF) # pour comparer graphiquement des df
library(DT) # pour l'affichage de df avec datable()

Functions

Code

# Fonction pour repérer les différences entre deux df
differences_summary <- function(df1, df2) {
  common_cols <- intersect(names(df1), names(df2)) # Colonnes communes
  differences <- lapply(common_cols, function(col) {
    rows_with_differences <- which(df1[[col]] != df2[[col]]) # Lignes avec différences
    if (length(rows_with_differences) > 0) {
      return(rows_with_differences) # Retourne les indices de lignes différentes
    } else {
      return(NULL) # Retourne NULL si pas de différences
    }
  })
  names(differences) <- common_cols # Associe les noms des colonnes
  return(differences)
}


#' transform a date in julian day
rvleExp.dateToNum_ad = function(dateStr) #' datStr, eg :  "2014-01-01"
{
  if (! is.character(dateStr)) {
    return(NA);
  } 
  return (as.numeric(as.Date(dateStr, format="%Y-%m-%d") + 2440588));
  #' the date in julian day: 2456659
}

# rvleExp.dateToNum_ad("2014-01-01")

#' transform a date of the form 2456659 into an object Date
rvleExp.dateFromNum_ad = function(dateNum, round=T) #' dateNum, eg :  2456659; #' round, are numerical rounded 
{
  if (round) {
    return(as.Date(as.Date(round(dateNum), origin="1970-01-01") - 2440588));
  } else {
    return(as.Date(as.Date(dateNum, origin="1970-01-01") - 2440588));
    #' return a date object, eg:  as.Date("2014-01-01")
  }
}
# rvleExp.dateFromNum_ad(2456659)

1 Case 1: Reference data data

Original data

1.1 Data import

Code

setwd("~/azodyn_check/AZODYN_optileg/R")
file_sim = "../data/pea_simul.csv";
file_obs = "../data/pea_observations.csv";
f = rvle.open(pkg="AZODYN", file="AzodynPea.vpz")
simData = rvleExp.parseSim(file_sim=file_sim, vleObj=f, id=NULL) 
# dim(simData) #  5 128; 5 137;  5 137; 5 136; 5 136
simData_cas1 = simData
obsData = rvleExp.parseObs(file_obs=file_obs, vleObj=f)

Code

# DT::datatable(simData, options = list(pageLength = 5))
tags$div(
  tags$h3("Overview of simulation data"),
  DT::datatable(simData, options = list(pageLength = 5))
)

Overview of simulation data

Code

tags$div(
  tags$h3("Overview of observation data"),
  DT::datatable(obsData, options = list(pageLength = 5))
)

Overview of observation data

1.2 Simulations

Code

res = rvle.plan_run(f)

[Manager] simulation mono nb simus: 5 
[Manager] end simulation and aggregation without parallelization

Code

res_cas1 = res
#rmse

Code

sim_vs_obs = rvleExp.compareSimObs(res=res, file_sim=simData, file_obs=obsData)

# dataframe vide pour stocker les résultats
rmse_cas1 <- data.frame(
  output = character(),
  RMSE = numeric(),
  stringsAsFactors = FALSE
)

# Boucle pour calculer la RMSE par groupe "output"
for (output in sort(unique(sim_vs_obs$output))) {
  sim_vs_obs_output <- sim_vs_obs[sim_vs_obs$output == output, ]
  
  # Calcul de la RMSE
  rmse_value <- rmse(sim_vs_obs_output$observed, sim_vs_obs_output$simulated)
  
  # Ajout des résultats au dataframe
  rmse_cas1 <- rbind(rmse_cas1, data.frame(output = output, RMSE = round (rmse_value,3)))
}

# rmse_cas1

1.3 Results case 1

ALL id

Code

rvleExp.plot(sim_vs_obs, file_sim=simData, file_obs=obsData)

[1] "rmse  LAI : 1.0569208998548"
[1] "rmse  MS : 122.076772209186"
[1] "rmse  NG : 1224.48033942446"
[1] "rmse  P1G : 0.0290073505444224"
[1] "rmse  Rdt : 13.2573989910306"

Focus on id 1

Code

rvleExp.plot(res, file_sim=simData, file_obs=obsData, 
             output_vars=NULL, id=1)

[1] "Rdt"
[1] "P1G"
[1] "NG"
[1] "MS"
[1] "LAI"

Focus on id 4

Code

rvleExp.plot(res, file_sim=simData, file_obs=obsData, 
             output_vars=NULL, id=4)

[1] "Rdt"
[1] "P1G"
[1] "NG"
[1] "MS"
[1] "LAI"

2 Case 2: Reference data data

Data modification

2.1 Data import

Code

setwd("~/azodyn_check/AZODYN_optileg/R")
file_sim = "../data/pea_simul_cas_4.csv"
file_obs = "../data/pea_observations.csv";
f = rvle.open(pkg="AZODYN", file="AzodynPea.vpz")
simData = rvleExp.parseSim(file_sim=file_sim, vleObj=f, id=NULL) 
# dim(simData) #  5 128; 5 137;  5 137; 5 136; 5 136
simData_cas2 = simData
obsData = rvleExp.parseObs(file_obs=file_obs, vleObj=f)

Code

# DT::datatable(simData, options = list(pageLength = 5))
tags$div(
  tags$h3("Overview of simulation data"),
  DT::datatable(simData, options = list(pageLength = 5))
)

Overview of simulation data

Code

tags$div(
  tags$h3("Overview of observation data"),
  DT::datatable(obsData, options = list(pageLength = 5))
)

Overview of observation data

2.2 Comparison of simulation data: Case1 vs Case2

1) Changing values

Code

df1 = simData_cas1
df2 = simData_cas2
# Comparaison

diffs <- differences_summary(df1, df2)

diffs <- diffs[!sapply(diffs, is.null)] # Retire les colonnes sans différences

# diffs


tags$div(
  tags$h3("Before modification"),
  DT::datatable(simData_cas1[, names(diffs)], options = list(pageLength = 5))
)

Before modification

Code

tags$div(
  tags$h3("After modification"),
  DT::datatable(simData_cas2[, names(diffs)], options = list(pageLength = 5))
)

After modification

2) Deleting conditions

Code

col_non_present2 = setdiff(names(df1),names(df2)) # colonne qui exite dans df1 mais qui n'existe pas dans df2
# df1[,col_non_present2]
# df2[,col_non_present2]
col_non_present2

[1] "cPlant.RDTmax_var"

3) Adding conditions

Code

col_non_present1 = setdiff(names(df2),names(df1)) # colonne qui exite dans df2 mais qui n'existe pas dans df1
# df1[,col_non_present1]
# df2[,col_non_present1]
# col_non_present1
tags$div(
  tags$h3(""),
  DT::datatable(df2[,col_non_present1], options = list(pageLength = 5))
)

2.3 Simulations

Code

res = rvle.plan_run(f)

[Manager] simulation mono nb simus: 5 
[Manager] end simulation and aggregation without parallelization

Code

res_cas2 = res
#rmse

Code

sim_vs_obs = rvleExp.compareSimObs(res=res, file_sim=simData, file_obs=obsData)

# dataframe vide pour stocker les résultats
rmse_cas2 <- data.frame(
  output = character(),
  RMSE = numeric(),
  stringsAsFactors = FALSE
)

# Boucle pour calculer la RMSE par groupe "output"
for (output in sort(unique(sim_vs_obs$output))) {
  sim_vs_obs_output <- sim_vs_obs[sim_vs_obs$output == output, ]
  
  # Calcul de la RMSE
  rmse_value <- rmse(sim_vs_obs_output$observed, sim_vs_obs_output$simulated)
  
  # Ajout des résultats au dataframe
  rmse_cas2 <- rbind(rmse_cas2, data.frame(output = output, RMSE = round (rmse_value,3)))
}

# rmse_cas2

2.4 Results case 2

ALL id

Code

rvleExp.plot(sim_vs_obs, file_sim=simData, file_obs=obsData)

[1] "rmse  LAI : 0.865361780796514"
[1] "rmse  MS : 62.6406879817228"
[1] "rmse  NG : 550.9981306973"
[1] "rmse  P1G : 0.0692214492904347"
[1] "rmse  Rdt : 21.4649640656617"

Focus on id 1

Code

rvleExp.plot(res, file_sim=simData, file_obs=obsData, 
             output_vars=NULL, id=1)

[1] "Rdt"
[1] "P1G"
[1] "NG"
[1] "MS"
[1] "LAI"

Focus on id 4

Code

rvleExp.plot(res, file_sim=simData, file_obs=obsData, 
             output_vars=NULL, id=4)

[1] "Rdt"
[1] "P1G"
[1] "NG"
[1] "MS"
[1] "LAI"

2.5 Comparison of RMSE: Case1 vs Case2

Code

# rmse_cas1
# rmse_cas2

rmse_cas_1_2 <- data.frame(
  Output = unique(rmse_cas1$output, rmse_cas2$output),
  RMSE_case_1 = rmse_cas1$RMSE,
  RMSE_case_2 = rmse_cas2$RMSE,
  stringsAsFactors = FALSE
)


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

RMSE

3 Case 3: OPTILEG data

Data preparation

3.1 Data import

Code

setwd("~/azodyn_check/AZODYN_optileg/R")
file_sim = "../data/pea_simul_test1.csv";
file_obs = "../data/pea_observations_test1.csv";
f = rvle.open(pkg="AZODYN", file="AzodynPea.vpz")
simData = rvleExp.parseSim(file_sim=file_sim, vleObj=f, id=NULL) 
# dim(simData) #  5 128; 5 137;  5 137; 5 136; 5 136
simData_cas3 = simData
obsData = rvleExp.parseObs(file_obs=file_obs, vleObj=f)

Code

# DT::datatable(simData, options = list(pageLength = 5))
tags$div(
  tags$h3("Overview of simulation data"),
  DT::datatable(simData, options = list(pageLength = 5))
)

Overview of simulation data

Code

tags$div(
  tags$h3("Overview of observation data"),
  DT::datatable(obsData, options = list(pageLength = 5))
)

Overview of observation data

3.2 Simulations

Code

res = rvle.plan_run(f)

[Manager] simulation mono nb simus: 16 
[Manager] end simulation and aggregation without parallelization

Code

res_cas3 = res
#rmse

Code

sim_vs_obs = rvleExp.compareSimObs(res=res, file_sim=simData, file_obs=obsData)

  # conversion dates simul en dates normal
lignes_dates <- which(grepl("date_", sim_vs_obs$output))
df_res_date = data.frame ( droplevels(sim_vs_obs [lignes_dates,]));
rownames(df_res_date) = NULL

  # unique(sim_vs_obs$output)
  sim_vs_obs = droplevels(sim_vs_obs[!sim_vs_obs$output %in% c("date_MP","date_DRG"),])
  sim_vs_obs$observed = as.numeric(sim_vs_obs$observed)

Code

# dataframe vide pour stocker les résultats
rmse_cas3 <- data.frame(
  output = character(),
  RMSE = numeric(),
  stringsAsFactors = FALSE
)

# Boucle pour calculer la RMSE par groupe "output"
for (output in sort(unique(sim_vs_obs$output))) {
  sim_vs_obs_output <- sim_vs_obs[sim_vs_obs$output == output, ]
  
  # Calcul de la RMSE
  rmse_value <- rmse(sim_vs_obs_output$observed, sim_vs_obs_output$simulated)
  
  # Ajout des résultats au dataframe
  rmse_cas3 <- rbind(rmse_cas3, data.frame(output = output, RMSE = round (rmse_value,3)))
}

# rmse_cas3

3.3 Results case 3

ALL id

Code

rvleExp.plot(sim_vs_obs, file_sim=simData, file_obs=obsData)

[1] "rmse  LAI : 0.348544126896241"
[1] "rmse  MS : 268.121329318601"
[1] "rmse  TP : 13.4015025420531"
[1] "rmse  NG : 852.684094477781"
[1] "rmse  P1G14 : 0.0938011293135174"
[1] "rmse  Rdt14 : 32.4827506969135"
[1] "rmse  tNa : 1.0109520757031"

Focus on id 1

Code

  obsData_modif = obsData
  obsData_modif$date_MP = NULL;   obsData_modif$date_DRG = NULL
  rvleExp.plot(res, file_sim=simData, file_obs=obsData_modif, 
               output_vars=NULL, id=1)

[1] "tNveg"
[1] "tNa"
[1] "date_DRG"
[1] "TP"
[1] "date_MP"
[1] "QNG"
[1] "P1G14"
[1] "Rdt14"
[1] "NG"
[1] "MS"
[1] "QNveg"
[1] "LAI"

Focus on id 77

Code

  rvleExp.plot(res, file_sim=simData, file_obs=obsData_modif, 
               output_vars=NULL, id=77)

[1] "tNveg"
[1] "tNa"
[1] "date_DRG"
[1] "TP"
[1] "date_MP"
[1] "QNG"
[1] "P1G14"
[1] "Rdt14"
[1] "NG"
[1] "MS"
[1] "QNveg"
[1] "LAI"

Focus on the dates

Code

df_res_date = df_res_date[which(df_res_date$simulated > 0),] # a verifier est ce ok
rownames(df_res_date) = NULL
df_res_date$simulated_cor = rvleExp.dateFromNum_ad(df_res_date$simulated)
# df_res_date[!is.na(df_res_date$observed), ]
df_res_date = df_res_date[df_res_date$observed != "", ]


tags$div(
  tags$h3("Comparison of the dates"),
  DT::datatable(df_res_date, options = list(pageLength = 5)), options = list(pageLength = 5))

Comparison of the dates

3.4 Results all cases

Code

# rmse_cas1
# rmse_cas2
# rmse_cas3

rmse_cas_1_2_3 = rbind (rmse_cas_1_2,
             data.frame(Output = rep(NA, 2), RMSE_case_1 = rep(NA, 2), RMSE_case_2 = rep(NA, 2))) 

rmse_cas_1_2_3$Output_cas3 = rmse_cas3$output
rmse_cas_1_2_3$RMSE_case_3 = rmse_cas3$RMSE

tags$div(
  tags$h3("RMSE"),
  DT::datatable(rmse_cas_1_2_3, options = list(pageLength = 7))
)

RMSE

Questions

Additional information

Thank you for your attention