class: center, middle, inverse, title-slide # Découverte de
airGRiwrm
## Gestion intégrée de la ressource en eau avec
airGR
### David Dorchies ### UMR G-EAU ### Montpellier / Paris, le 30 septembre 2021 (MAJ 30/09/2021) ---
# Programme de la formation - 09h00: Introduction - 09h10: Installation du package **airGRiwrm** - 09h20: Le jeu de données de la formation - 09h30: Le modèle semi-distribué dans **airGR** - 10h00: **airGRiwrm**: une extension de **airGR** pour les modèles semi-distribués - 10h10: Créer un réseau semi-distribué de modèles avec **airGRiwrm** - 10h30: Calage d'un réseau semi-distribué de modèles hydrologiques - 11h00: Calage d'un réseau de modèles avec des débits influencés - 11h20: Régularisation des paramètres - 11h40: Régulation automatique centralisée - 11h50: Discussions autour de la famille des packages **airGR** --- # Introduction - Le package airGRiwrm est en développement depuis mai 2020 - Testé sur la Seine (stage 2021 Laura Nunez) avec 143 sous-bassins et l'influence des grands lacs - Forge : https://gitlab.irstea.fr/in-wop/airGRiwrm - Site internet : https://airgriwrm.g-eau.fr/dev/ - Lien de téléchargement de la formation : https://link.infini.fr/airgriwrm_2021-09 --- class: inverse, middle, center # Installation des packages *airGR* et *airGRiwrm* --- # Installation des packages *airGR* et *airGRiwrm* ### Sous Linux ou sous Windows avec les packages `devtools` et `Rtools` installés: ```r remotes::install_gitlab("in-wop/airgriwrm", host = "gitlab.irstea.fr", dependencies = TRUE) ``` ### Sinon, installation des packages locaux sous Windows Installation des dépendances : ```r install.packages(c("dplyr", "DiagrammeR")) ``` Installation de **airGR** et **airGRiwrm** version "dev" : ```r install.packages("../airGR_1.6.12.9000.zip", repo = NULL) install.packages("../airGRiwrm_0.5.0.9000.tar.gz", repo = NULL) ``` --- class: inverse, middle, center # Le jeu de données de la formation --- # Présentation du cas d'étude .pull-left[ ## 4 stations hydrométriques autour du lac "Seine" Données disponibles au pas de temps journalier: - Données météo des sous-bassins (P, ETP) (Source SAFRAN) - Débits observés des stations hydro (Source Banque hydro) - Débits observés en entrée (prise) et sortie (restitution) du lac "Seine" ] .pull-right[
] --- # Cas d'étude: chargement des données - Créer un R Notebook sous RStudio - Enregistrer le dossier `data` et le R Notebook dans le même dossier - Définir le working directory (Menu Session > Set working directory > To source file location) - Créer un premier chunk avec le code ci-dessous ```r load("data/Seine.RData") ls() ``` ``` ## [1] "BasinObs" "conn_seine" "DatesR" ## [4] "map" "Q_reservoir_seine" "reseau_reservoir_seine" ## [7] "sfSeine" ``` - `BasinObs` (*list*): Données hydroclimatiques de chaque bassin - `DatesR` (*POSIXct*): Dates des pas de temps d'observation - `Q_reservoir_seine` (*matrix*): Débits observés à la prise et à la restitution du lac Seine - `reseau_reservoir_seine` (*data.frame*): Description du réseau de stations --- ## Cas d'étude: description du réseau ```r reseau_reservoir_seine ``` |id | area|nom |id_aval | distance_aval| |:--------|-------:|:-----------------------------------------------------|:--------|-------------:| |H0800010 | 3499.50|La Seine [totale] à Troyes [après création grand lac] |NA | NA| |H0400010 | 2340.37|La Seine à Bar-sur-Seine |H0800010 | 45.81| |H0210010 | 1462.66|La Seine à Polisy |H0400010 | 7.86| |H0321030 | 548.93|L'Ource à Autricourt |H0400010 | 32.34| --- # Cas d'étude: description des obs ```r head(DatesR, 3) ``` ``` ## [1] "1958-08-01 UTC" "1958-08-02 UTC" "1958-08-03 UTC" ``` ```r tail(DatesR, 3) ``` ``` ## [1] "2019-07-29 UTC" "2019-07-30 UTC" "2019-07-31 UTC" ``` ```r str(BasinObs) ``` ``` ## List of 4 ## $ H0800010:'data.frame': 22280 obs. of 3 variables: ## ..$ P: num [1:22280] 1.481 0 0 0 0.367 ... ## ..$ E: num [1:22280] 4.24 3.34 3.07 3.29 3.92 ... ## ..$ Q: num [1:22280] 0.193 0.189 0.191 0.198 0.164 0.159 0.16 0.156 0.152 0.143 ... ## $ H0400010:'data.frame': 22280 obs. of 3 variables: ## ..$ P: num [1:22280] 5.92 0 0 0 1.23 ... ## ..$ E: num [1:22280] 4.27 3.34 2.96 3.24 3.83 ... ## ..$ Q: num [1:22280] 0.249 0.229 0.229 0.225 0.238 0.229 0.225 0.229 0.229 0.229 ... ## $ H0210010:'data.frame': 22280 obs. of 3 variables: ## ..$ P: num [1:22280] 6.16 0 0 0 1.24 ... ## ..$ E: num [1:22280] 4.3 3.33 3 3.23 3.88 ... ## ..$ Q: num [1:22280] 0.307 0.284 0.284 0.254 0.254 0.284 0.284 0.254 0.254 0.23 ... ## $ H0321030:'data.frame': 22280 obs. of 3 variables: ## ..$ P: num [1:22280] 0.877 0 0 0 1.061 ... ## ..$ E: num [1:22280] 4.43 3.46 3.09 3.32 3.8 ... ## ..$ Q: num [1:22280] NA NA NA NA NA NA NA NA NA NA ... ``` --- class: inverse, middle, center # Le modèle semi-distribué dans *airGR* --- # Modèle SD dans *airGR* ## Un modèle semi-distribué : c'est quoi ? .pull-left[ - Découpage du bassin versant en sous-bassins en fonction de la position des stations hydrométriques exutoires de chaque sous-bassin - Le ruissellement de chaque sous-bassin est modélisé par un modèle global (GR...) - Le modèle semi-distribué relie les stations hydrométrique grâce à un routage hydraulique de type "lag" ] .pull-right[ ![:scale 40%](img/Schema_BV_amont.png) ![:scale 40%](img/Schema_BVI_aval.png) ] --- # Modèle SD dans *airGR* ## La Seine à Bar-sur-Seine avec 2 apports amont .pull-left[
] .pull-right[ Modélisation du débit sur la station *H0400010* avec un modèle GR4J + Lag ayant comme apports amonts les débits observés aux stations *H0210010* et *H0321030.*
] --- # Modèle SD dans *airGR* ## La procédure de calage Pour effectuer un calage, il faut appeler les fonctions suivantes créant les objets ad-hoc:
--- # Modèle SD dans *airGR* ## `CreateInputsModel`: ```r library(airGR) ?airGR::CreateInputsModel ``` -- Il faut définir : -- - L'aire des bassins amont et du sous-bassin aval -- ```r BasinAreas <- reseau_reservoir_seine[c(3, 4, 2), "area"] BasinAreas[3] <- BasinAreas[3] - sum(BasinAreas[1:2]) ``` -- - La matrice des apports amont provenant des stations H0210010 et H0321030 -- ```r Qupstream <- cbind(BasinObs[["H0210010"]]$Q, BasinObs[["H0321030"]]$Q) ``` -- - La distance entre les apports amont et l'exutoire -- ```r # Distance entre les apports amont et l'exutoire LengthHydro <- reseau_reservoir_seine$distance_aval[3:4] ``` --- # Modèle SD dans *airGR* ## `CreateInputsModel`: .pull-left[ ```r IMairGR <- CreateInputsModel( FUN_MOD = RunModel_GR4J, DatesR = DatesR, Precip = BasinObs[["H0400010"]]$P, PotEvap = BasinObs[["H0400010"]]$E, * Qupstream = Qupstream, * LengthHydro = LengthHydro, * BasinAreas = BasinAreas ) ``` ``` ## Warning in CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = DatesR, Precip = ## BasinObs[["H0400010"]]$P, : 'Qupstream' contains NA values: model outputs will ## contain NAs ``` ] .pull-right[ ```r str(IMairGR) ``` ``` ## List of 6 ## $ DatesR : POSIXlt[1:22280], format: "1958-08-01" "1958-08-02" ... ## $ Precip : num [1:22280] 5.92 0 0 0 1.23 ... ## $ PotEvap : num [1:22280] 4.27 3.34 2.96 3.24 3.83 ... ## $ Qupstream : num [1:22280, 1:2] 449037 415395 415395 371516 371516 ... ## $ LengthHydro: num [1:2] 7.86 32.34 ## $ BasinAreas : num [1:3] 1463 549 329 ## - attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD" ``` ] --- # Modèle SD dans *airGR* ## `CreateRunOptions` ```r ?airGR::CreateRunOptions ``` -- ```r I_Run <- 366:length(DatesR) ROairGR <- CreateRunOptions( FUN_MOD = RunModel_GR4J, InputsModel = IMairGR, IndPeriod_WarmUp = 1:365, IndPeriod_Run = I_Run ) ``` --- # Modèle SD dans *airGR* ## `CreateInputsCrit`: ```r ?airGR::CreateInputsCrit ``` -- On va effectuer le calage sur la racine carrée des débits ```r ICairGR <- CreateInputsCrit( FUN_CRIT = ErrorCrit_KGE2, InputsModel = IMairGR, RunOptions = ROairGR, Obs = BasinObs[["H0400010"]]$Q[I_Run], * transfo = "sqrt" ) ``` --- # Modèle SD dans *airGR* ## `CreateCalibOptions`: ```r ?airGR::CreateCalibOptions ``` -- Il faut spécifier que le modéle est semi-distribué... -- ```r COairGR <- CreateCalibOptions( FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel, * IsSD = TRUE ) ``` --- # Modèle SD dans *airGR* ```r ?airGR::Calibration ``` La calibration a lieu comme pour n'importe quel modèle de *airGR*. -- On ajoute `suppressWarnings` pour éviter les avertissements liés au `NAs` dans les débits amont ```r OCairGR <- suppressWarnings( Calibration( InputsModel = IMairGR, RunOptions = ROairGR, InputsCrit = ICairGR, CalibOptions = COairGR, FUN_MOD = RunModel_GR4J )) ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 15.000, 432.681, -2.376, 20.697, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8303 ## Steepest-descent local search in progress ## Calibration completed (32 iterations, 532 runs) ## Param = 19.990, 371.074, -1002.215, 20.697, 4.493 ## Crit. KGE2[sqrt(Q)] = 0.8682 ``` --- # Modèle SD dans *airGR* .pull-left[ ## RunModel ```r OMairGR <- RunModel( InputsModel = IMairGR, RunOptions = ROairGR, Param = OCairGR$ParamFinalR, FUN_MOD = RunModel_GR4J ) ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 11553 ## time steps with NA values ``` ```r plot(OMairGR, BasinObs[["H0400010"]]$Q[I_Run]) ``` ] .pull-right[ ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-30-1.png)<!-- --> ] --- class: inverse, middle, center # *airGRiwrm*: extension de *airGR* pour un réseau de modèles semi-distribués --- class: hide_logo # *airGRiwrm*: extension de *airGR* pour un réseau de modèles semi-distribués **airGRiwrm** redéfinit les fonctions de **airGR** pour automatiser la gestion d'un réseau semi-distribué: ```r library(airGR) library(airGRiwrm) ``` ``` ## ## Attaching package: 'airGRiwrm' ``` ``` ## The following objects are masked from 'package:airGR': ## ## Calibration, CreateCalibOptions, CreateInputsCrit, ## CreateInputsModel, CreateRunOptions, RunModel ``` Les fonctions originales de **airGR** restent accessibles (Redirection avec les classes S3). --- class: inverse, middle, center # Créer un réseau de modèles semi-distribués avec *airGRiwrm* --- # Procédure de calage d'un réseau avec *airGRiwrm* Par rapport à *airGR*, on ajoute une étape préalable de description du réseau
--- # Créer un réseau de modèles semi-distribués avec *airGRiwrm* Le réseau est décrit par un objet `GRiwrm` qu'il faut créer : ```r ?CreateGRiwrm ``` -- ```r reseau_reservoir_seine ``` |id | area|nom |id_aval | distance_aval| |:--------|-------:|:-----------------------------------------------------|:--------|-------------:| |H0800010 | 3499.50|La Seine [totale] à Troyes [après création grand lac] |NA | NA| |H0400010 | 2340.37|La Seine à Bar-sur-Seine |H0800010 | 45.81| |H0210010 | 1462.66|La Seine à Polisy |H0400010 | 7.86| |H0321030 | 548.93|L'Ource à Autricourt |H0400010 | 32.34| -- Il faut définir quel modèle hydrologique on souhaite utilisé pour chaque noeud du réseau : -- ```r reseau_reservoir_seine$model <- "RunModel_GR4J" ``` --- # Création de l'objet `GRiwrm` On peut faire correspondre les noms de colonne avec le paramètre `cols` : -- ```r griwrmNat <- CreateGRiwrm( reseau_reservoir_seine, cols = list(down = "id_aval", length = "distance_aval")) str(griwrmNat) ``` ``` ## Classes 'GRiwrm' and 'data.frame': 4 obs. of 5 variables: ## $ id : chr "H0800010" "H0400010" "H0210010" "H0321030" ## $ down : chr NA "H0800010" "H0400010" "H0400010" ## $ length: num NA 45.81 7.86 32.34 ## $ model : chr "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" ## $ area : num 3500 2340 1463 549 ``` --- # Créer un réseau de modèles semi-distribués avec **airGRiwrm** La fonction `plot` permet de tracer une représentation schématique du réseau : ```r plot(griwrmNat) ``` .pull-left[
] .pull-right[ Code couleur: - rouge pour un noeud d'entrée directe (pas de simulation hydrologique) - bleu pour un noeud amont (modèle GR) - vert pour noeud intermédiaire/aval (modèle GR + LAG) ] --- class: inverse, middle, center # Calage d'un réseau de modèles hydrologiques semi-distribués avec *airGRiwrm* --- # `CreateInputsModel` dans *airGRiwrm* ```r ?CreateInputsModel.GRiwrm ``` -- Il faut rassembler les "inputs" dans des matrices avec le code de la station en nom de colonne : ```r P <- do.call(cbind, lapply(BasinObs, "[[", "P")) str(P) ``` ``` ## num [1:22280, 1:4] 1.481 0 0 0 0.367 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : NULL ## ..$ : chr [1:4] "H0800010" "H0400010" "H0210010" "H0321030" ``` -- ```r E <- do.call(cbind, lapply(BasinObs, "[[", "E")) Q <- do.call(cbind, lapply(BasinObs, "[[", "Q")) ``` --- # `CreateInputsModel` dans *airGRiwrm* ```r IMnoRes <- CreateInputsModel( griwrmNat, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = Q ) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0210010... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0321030... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0400010... ``` ``` ## Warning in airGR::CreateInputsModel(FUN_MOD = x, ...): 'Qupstream' contains NA ## values: model outputs will contain NAs ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0800010... ``` ``` ## Warning in airGR::CreateInputsModel(FUN_MOD = x, ...): 'Qupstream' contains NA ## values: model outputs will contain NAs ``` --- # L'objet `GRiwrmInputsModel` L'objet `GRiwrmInputsModel` est une liste d'`InputsModel` : 1 élément de liste = 1 modèle *GR ( + Lag )* ```r str(IMnoRes) ``` ``` ## List of 4 ## $ H0210010:List of 7 ## ..$ DatesR : POSIXlt[1:22280], format: "1958-08-01" "1958-08-02" ... ## ..$ Precip : num [1:22280] 6.16 0 0 0 1.24 ... ## ..$ PotEvap : num [1:22280] 4.3 3.33 3 3.23 3.88 ... ## ..$ id : chr "H0210010" ## ..$ down : chr "H0400010" ## ..$ BasinAreas: num 1463 ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR" ## $ H0321030:List of 7 ## ..$ DatesR : POSIXlt[1:22280], format: "1958-08-01" "1958-08-02" ... ## ..$ Precip : num [1:22280] 0.877 0 0 0 1.061 ... ## ..$ PotEvap : num [1:22280] 4.43 3.46 3.09 3.32 3.8 ... ## ..$ id : chr "H0321030" ## ..$ down : chr "H0400010" ## ..$ BasinAreas: num 549 ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR" ## $ H0400010:List of 11 ## ..$ DatesR : POSIXlt[1:22280], format: "1958-08-01" "1958-08-02" ... ## ..$ Precip : num [1:22280] 5.92 0 0 0 1.23 ... ## ..$ PotEvap : num [1:22280] 4.27 3.34 2.96 3.24 3.83 ... ## ..$ Qupstream : num [1:22280, 1:2] 449037 415395 415395 371516 371516 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "H0210010" "H0321030" ## ..$ LengthHydro : Named num [1:2] 7.86 32.34 ## .. ..- attr(*, "names")= chr [1:2] "H0210010" "H0321030" ## ..$ BasinAreas : Named num [1:3] 1463 549 329 ## .. ..- attr(*, "names")= chr [1:3] "H0210010" "H0321030" "H0400010" ## ..$ id : chr "H0400010" ## ..$ down : chr "H0800010" ## ..$ UpstreamNodes : chr [1:2] "H0210010" "H0321030" ## ..$ UpstreamIsRunoff: logi [1:2] TRUE TRUE ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD" ## $ H0800010:List of 11 ## ..$ DatesR : POSIXlt[1:22280], format: "1958-08-01" "1958-08-02" ... ## ..$ Precip : num [1:22280] 1.481 0 0 0 0.367 ... ## ..$ PotEvap : num [1:22280] 4.24 3.34 3.07 3.29 3.92 ... ## ..$ Qupstream : num [1:22280, 1] 582752 535945 535945 526583 557008 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr "H0400010" ## ..$ LengthHydro : Named num 45.8 ## .. ..- attr(*, "names")= chr "H0400010" ## ..$ BasinAreas : Named num [1:2] 2340 1159 ## .. ..- attr(*, "names")= chr [1:2] "H0400010" "H0800010" ## ..$ id : chr "H0800010" ## ..$ down : chr NA ## ..$ UpstreamNodes : chr "H0400010" ## ..$ UpstreamIsRunoff: logi TRUE ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD" ## - attr(*, "class")= chr [1:2] "GRiwrmInputsModel" "list" ## - attr(*, "GRiwrm")=Classes 'GRiwrm' and 'data.frame': 4 obs. of 5 variables: ## ..$ id : chr [1:4] "H0800010" "H0400010" "H0210010" "H0321030" ## ..$ down : chr [1:4] NA "H0800010" "H0400010" "H0400010" ## ..$ length: num [1:4] NA 45.81 7.86 32.34 ## ..$ model : chr [1:4] "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" ## ..$ area : num [1:4] 3500 2340 1463 549 ## - attr(*, "TimeStep")= num 86400 ``` --- # Calage avec *airGRiwrm* .pull-left[ ## CreateRunOptions ```r ?CreateRunOptions.GRiwrmInputsModel ``` ] -- .pull-right[ ```r ROnoRes <- CreateRunOptions( InputsModel = IMnoRes, IndPeriod_WarmUp = 1:365, IndPeriod_Run = I_Run ) ``` ] ```r str(ROnoRes) ``` ``` ## List of 4 ## $ H0210010:List of 8 ## ..$ IndPeriod_WarmUp: int [1:365] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ IndPeriod_Run : int [1:21915] 366 367 368 369 370 371 372 373 374 375 ... ## ..$ IniStates : num [1:67] 0 0 0 0 0 0 0 0 0 0 ... ## ..$ IniResLevels : num [1:4] 0.3 0.5 NA NA ## ..$ Outputs_Cal : chr [1:2] "Qsim" "Param" ## ..$ Outputs_Sim : Named chr [1:22] "DatesR" "PotEvap" "Precip" "Prod" ... ## .. ..- attr(*, "names")= chr [1:22] "" "GR1" "GR2" "GR3" ... ## ..$ FortranOutputs :List of 2 ## .. ..$ GR: chr [1:18] "PotEvap" "Precip" "Prod" "Pn" ... ## .. ..$ CN: NULL ## ..$ FeatFUN_MOD :List of 12 ## .. ..$ CodeMod : chr "GR4J" ## .. ..$ NameMod : chr "GR4J" ## .. ..$ NbParam : int 4 ## .. ..$ TimeUnit : chr "daily" ## .. ..$ Id : logi NA ## .. ..$ Class : chr [1:2] "daily" "GR" ## .. ..$ Pkg : chr "airGR" ## .. ..$ NameFunMod : chr "RunModel_GR4J" ## .. ..$ TimeStep : num 86400 ## .. ..$ TimeStepMean: int 86400 ## .. ..$ CodeModHydro: chr "GR4J" ## .. ..$ IsSD : logi FALSE ## ..- attr(*, "class")= chr [1:3] "RunOptions" "daily" "GR" ## $ H0321030:List of 8 ## ..$ IndPeriod_WarmUp: int [1:365] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ IndPeriod_Run : int [1:21915] 366 367 368 369 370 371 372 373 374 375 ... ## ..$ IniStates : num [1:67] 0 0 0 0 0 0 0 0 0 0 ... ## ..$ IniResLevels : num [1:4] 0.3 0.5 NA NA ## ..$ Outputs_Cal : chr [1:2] "Qsim" "Param" ## ..$ Outputs_Sim : Named chr [1:22] "DatesR" "PotEvap" "Precip" "Prod" ... ## .. ..- attr(*, "names")= chr [1:22] "" "GR1" "GR2" "GR3" ... ## ..$ FortranOutputs :List of 2 ## .. ..$ GR: chr [1:18] "PotEvap" "Precip" "Prod" "Pn" ... ## .. ..$ CN: NULL ## ..$ FeatFUN_MOD :List of 12 ## .. ..$ CodeMod : chr "GR4J" ## .. ..$ NameMod : chr "GR4J" ## .. ..$ NbParam : int 4 ## .. ..$ TimeUnit : chr "daily" ## .. ..$ Id : logi NA ## .. ..$ Class : chr [1:2] "daily" "GR" ## .. ..$ Pkg : chr "airGR" ## .. ..$ NameFunMod : chr "RunModel_GR4J" ## .. ..$ TimeStep : num 86400 ## .. ..$ TimeStepMean: int 86400 ## .. ..$ CodeModHydro: chr "GR4J" ## .. ..$ IsSD : logi FALSE ## ..- attr(*, "class")= chr [1:3] "RunOptions" "daily" "GR" ## $ H0400010:List of 8 ## ..$ IndPeriod_WarmUp: int [1:365] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ IndPeriod_Run : int [1:21915] 366 367 368 369 370 371 372 373 374 375 ... ## ..$ IniStates : num [1:67] 0 0 0 0 0 0 0 0 0 0 ... ## ..$ IniResLevels : num [1:4] 0.3 0.5 NA NA ## ..$ Outputs_Cal : chr [1:2] "Qsim" "Param" ## ..$ Outputs_Sim : Named chr [1:24] "DatesR" "PotEvap" "Precip" "Prod" ... ## .. ..- attr(*, "names")= chr [1:24] "" "GR1" "GR2" "GR3" ... ## ..$ FortranOutputs :List of 2 ## .. ..$ GR: chr [1:18] "PotEvap" "Precip" "Prod" "Pn" ... ## .. ..$ CN: NULL ## ..$ FeatFUN_MOD :List of 12 ## .. ..$ CodeMod : chr "GR4J" ## .. ..$ NameMod : chr "GR4J" ## .. ..$ NbParam : num 5 ## .. ..$ TimeUnit : chr "daily" ## .. ..$ Id : logi NA ## .. ..$ Class : chr [1:2] "daily" "GR" ## .. ..$ Pkg : chr "airGR" ## .. ..$ NameFunMod : chr "RunModel_GR4J" ## .. ..$ TimeStep : num 86400 ## .. ..$ TimeStepMean: int 86400 ## .. ..$ CodeModHydro: chr "GR4J" ## .. ..$ IsSD : logi TRUE ## ..- attr(*, "class")= chr [1:3] "RunOptions" "daily" "GR" ## $ H0800010:List of 8 ## ..$ IndPeriod_WarmUp: int [1:365] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ IndPeriod_Run : int [1:21915] 366 367 368 369 370 371 372 373 374 375 ... ## ..$ IniStates : num [1:67] 0 0 0 0 0 0 0 0 0 0 ... ## ..$ IniResLevels : num [1:4] 0.3 0.5 NA NA ## ..$ Outputs_Cal : chr [1:2] "Qsim" "Param" ## ..$ Outputs_Sim : Named chr [1:24] "DatesR" "PotEvap" "Precip" "Prod" ... ## .. ..- attr(*, "names")= chr [1:24] "" "GR1" "GR2" "GR3" ... ## ..$ FortranOutputs :List of 2 ## .. ..$ GR: chr [1:18] "PotEvap" "Precip" "Prod" "Pn" ... ## .. ..$ CN: NULL ## ..$ FeatFUN_MOD :List of 12 ## .. ..$ CodeMod : chr "GR4J" ## .. ..$ NameMod : chr "GR4J" ## .. ..$ NbParam : num 5 ## .. ..$ TimeUnit : chr "daily" ## .. ..$ Id : logi NA ## .. ..$ Class : chr [1:2] "daily" "GR" ## .. ..$ Pkg : chr "airGR" ## .. ..$ NameFunMod : chr "RunModel_GR4J" ## .. ..$ TimeStep : num 86400 ## .. ..$ TimeStepMean: int 86400 ## .. ..$ CodeModHydro: chr "GR4J" ## .. ..$ IsSD : logi TRUE ## ..- attr(*, "class")= chr [1:3] "RunOptions" "daily" "GR" ## - attr(*, "class")= chr [1:2] "list" "GRiwrmRunOptions" ``` --- # Calage avec *airGRiwrm* .pull-left[ ## CreateInputsCrit ```r ?CreateInputsCrit.GRiwrmInputsModel ``` ] -- .pull-right[ ```r ICnoRes <- CreateInputsCrit( InputsModel = IMnoRes, FUN_CRIT = airGR::ErrorCrit_KGE2, RunOptions = ROnoRes, Obs = Q[I_Run,], transfo = "sqrt" ) ``` ] ```r str(ICnoRes) ``` ``` ## List of 4 ## $ H0210010:List of 8 ## ..$ FUN_CRIT:function (InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 19 109 1 19 1 2628 2736 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003dd20f68> ## .. ..- attr(*, "class")= chr [1:2] "FUN_CRIT" "function" ## ..$ Obs : num [1:21915] 0.337 0.307 0.284 0.254 0.23 0.195 0.154 0.154 0.136 0.136 ... ## ..$ VarObs : chr "Q" ## ..$ BoolCrit: logi [1:21915] TRUE TRUE TRUE TRUE TRUE TRUE ... ## ..$ idLayer : logi NA ## ..$ transfo : chr "sqrt" ## ..$ epsilon : NULL ## ..$ Weights : NULL ## ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit" ## $ H0321030:List of 8 ## ..$ FUN_CRIT:function (InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 19 109 1 19 1 2628 2736 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003dd20f68> ## .. ..- attr(*, "class")= chr [1:2] "FUN_CRIT" "function" ## ..$ Obs : num [1:21915] NA NA NA NA NA NA NA NA NA NA ... ## ..$ VarObs : chr "Q" ## ..$ BoolCrit: logi [1:21915] TRUE TRUE TRUE TRUE TRUE TRUE ... ## ..$ idLayer : logi NA ## ..$ transfo : chr "sqrt" ## ..$ epsilon : NULL ## ..$ Weights : NULL ## ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit" ## $ H0400010:List of 8 ## ..$ FUN_CRIT:function (InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 19 109 1 19 1 2628 2736 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003dd20f68> ## .. ..- attr(*, "class")= chr [1:2] "FUN_CRIT" "function" ## ..$ Obs : num [1:21915] 0.118 0.111 0.111 0.118 0.111 0.111 0.107 0.107 0.107 0.103 ... ## ..$ VarObs : chr "Q" ## ..$ BoolCrit: logi [1:21915] TRUE TRUE TRUE TRUE TRUE TRUE ... ## ..$ idLayer : logi NA ## ..$ transfo : chr "sqrt" ## ..$ epsilon : NULL ## ..$ Weights : NULL ## ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit" ## $ H0800010:List of 8 ## ..$ FUN_CRIT:function (InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 19 109 1 19 1 2628 2736 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003dd20f68> ## .. ..- attr(*, "class")= chr [1:2] "FUN_CRIT" "function" ## ..$ Obs : num [1:21915] 0.081 0.077 0.074 0.079 0.077 0.074 0.072 0.067 0.064 0.067 ... ## ..$ VarObs : chr "Q" ## ..$ BoolCrit: logi [1:21915] TRUE TRUE TRUE TRUE TRUE TRUE ... ## ..$ idLayer : logi NA ## ..$ transfo : chr "sqrt" ## ..$ epsilon : NULL ## ..$ Weights : NULL ## ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit" ## - attr(*, "class")= chr [1:2] "GRiwrmInputsCrit" "list" ``` --- # Calage avec *airGRiwrm* .pull-left[ ## CreateCalibOptions ```r ?CreateInputsCrit.GRiwrmInputsModel ``` ] -- .pull-right[ ```r COnoRes <- CreateCalibOptions(IMnoRes) ``` ] ```r str(COnoRes) ``` ``` ## List of 4 ## $ H0210010:List of 4 ## ..$ FixedParam : logi [1:4] NA NA NA NA ## ..$ SearchRanges : num [1:2, 1:4] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ... ## ..$ FUN_TRANSFO :function (ParamIn, Direction) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 22 48 1 22 1 5580 5627 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003e7c1178> ## ..$ StartParamDistrib: num [1:3, 1:4] 169.017 247.151 432.681 -2.376 -0.649 ... ## ..- attr(*, "class")= chr [1:4] "CalibOptions" "daily" "GR" "HBAN" ## $ H0321030:List of 4 ## ..$ FixedParam : logi [1:4] NA NA NA NA ## ..$ SearchRanges : num [1:2, 1:4] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ... ## ..$ FUN_TRANSFO :function (ParamIn, Direction) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 1 22 48 1 22 1 5580 5627 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003e7c1178> ## ..$ StartParamDistrib: num [1:3, 1:4] 169.017 247.151 432.681 -2.376 -0.649 ... ## ..- attr(*, "class")= chr [1:4] "CalibOptions" "daily" "GR" "HBAN" ## $ H0400010:List of 4 ## ..$ FixedParam : logi [1:5] NA NA NA NA NA ## ..$ SearchRanges : num [1:2, 1:5] 1.00e-02 2.00e+01 4.59e-05 2.18e+04 -1.09e+04 ... ## ..$ FUN_TRANSFO :function (ParamIn, Direction) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 37 22 50 7 22 7 6123 6136 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003e6bfd58> ## ..$ StartParamDistrib: num [1:3, 1:5] 11.2 12.5 15 169 247.2 ... ## ..- attr(*, "class")= chr [1:5] "CalibOptions" "daily" "GR" "SD" ... ## $ H0800010:List of 4 ## ..$ FixedParam : logi [1:5] NA NA NA NA NA ## ..$ SearchRanges : num [1:2, 1:5] 1.00e-02 2.00e+01 4.59e-05 2.18e+04 -1.09e+04 ... ## ..$ FUN_TRANSFO :function (ParamIn, Direction) ## .. ..- attr(*, "srcref")= 'srcref' int [1:8] 37 22 50 7 22 7 6123 6136 ## .. .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x000000003e6bfd58> ## ..$ StartParamDistrib: num [1:3, 1:5] 11.2 12.5 15 169 247.2 ... ## ..- attr(*, "class")= chr [1:5] "CalibOptions" "daily" "GR" "SD" ... ## - attr(*, "class")= chr [1:2] "GRiwrmCalibOptions" "list" ``` --- # `Calibration` dans *airGRiwrm* ```r ?Calibration.GRiwrmInputsModel ``` -- ```r OCnoRes <- Calibration(InputsModel = IMnoRes, RunOptions = ROnoRes, InputsCrit = ICnoRes, CalibOptions = COnoRes) ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 432.681, -0.020, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8845 ## Steepest-descent local search in progress ## Calibration completed (44 iterations, 438 runs) ## Param = 275.925, 0.219, 99.248, 5.068 ## Crit. KGE2[sqrt(Q)] = 0.9339 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0321030... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 169.017, -0.649, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.9089 ## Steepest-descent local search in progress ## Calibration completed (25 iterations, 266 runs) ## Param = 177.977, -0.513, 59.607, 5.425 ## Crit. KGE2[sqrt(Q)] = 0.9380 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0400010... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 247.151, -2.376, 20.697, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8828 ## Steepest-descent local search in progress ## Calibration completed (57 iterations, 798 runs) ## Param = 2.839, 73.774, -29.099, 19.651, 3.003 ## Crit. KGE2[sqrt(Q)] = 0.8968 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0800010... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 432.681, -0.020, 83.096, 1.417 ## Crit. KGE2[sqrt(Q)] = 0.7232 ## Steepest-descent local search in progress ## Calibration completed (47 iterations, 695 runs) ## Param = 0.360, 2275.602, -0.120, 31.187, 2.462 ## Crit. KGE2[sqrt(Q)] = 0.7764 ``` --- # `Calibration` dans *airGRiwrm* ```r str(OCnoRes) ``` ``` ## List of 4 ## $ H0210010:List of 9 ## ..$ ParamFinalR : num [1:4] 275.925 0.219 99.248 5.068 ## ..$ CritFinal : num 0.934 ## ..$ NIter : num 44 ## ..$ NRuns : num 438 ## ..$ HistParamR : num [1:44, 1:4] 433 433 433 433 433 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:4] "Param1" "Param2" "Param3" "Param4" ## ..$ HistCrit : num [1:44, 1] 0.885 0.889 0.893 0.896 0.903 ... ## ..$ MatBoolCrit : logi [1:21915, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "KGE2[sqrt(Q)]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ H0321030:List of 9 ## ..$ ParamFinalR : num [1:4] 177.977 -0.513 59.607 5.425 ## ..$ CritFinal : num 0.938 ## ..$ NIter : num 25 ## ..$ NRuns : num 266 ## ..$ HistParamR : num [1:25, 1:4] 169 169 169 169 169 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:4] "Param1" "Param2" "Param3" "Param4" ## ..$ HistCrit : num [1:25, 1] 0.909 0.917 0.924 0.926 0.927 ... ## ..$ MatBoolCrit : logi [1:21915, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "KGE2[sqrt(Q)]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ H0400010:List of 9 ## ..$ ParamFinalR : num [1:5] 2.84 73.77 -29.1 19.65 3 ## ..$ CritFinal : num 0.897 ## ..$ NIter : num 57 ## ..$ NRuns : num 798 ## ..$ HistParamR : num [1:57, 1:5] 11.2 11.2 11.2 11.2 11.2 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:5] "Param1" "Param2" "Param3" "Param4" ... ## ..$ HistCrit : num [1:57, 1] 0.883 0.89 0.894 0.896 0.896 ... ## ..$ MatBoolCrit : logi [1:21915, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "KGE2[sqrt(Q)]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ H0800010:List of 9 ## ..$ ParamFinalR : num [1:5] 0.36 2275.6 -0.12 31.19 2.46 ## ..$ CritFinal : num 0.776 ## ..$ NIter : num 47 ## ..$ NRuns : num 695 ## ..$ HistParamR : num [1:47, 1:5] 11.2 11.2 11.2 11.2 10.6 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:5] "Param1" "Param2" "Param3" "Param4" ... ## ..$ HistCrit : num [1:47, 1] 0.723 0.756 0.772 0.772 0.772 ... ## ..$ MatBoolCrit : logi [1:21915, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "KGE2[sqrt(Q)]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## - attr(*, "class")= chr [1:2] "GRiwrmOutputsCalib" "list" ``` --- # `RunModel` dans *airGRiwrm* .pull-left[ Il faut une liste contenant les paramètres de chaque sous-bassin : ```r ParamNoRes <- lapply(OCnoRes, "[[", "ParamFinalR") str(ParamNoRes) ``` ``` ## List of 4 ## $ H0210010: num [1:4] 275.925 0.219 99.248 5.068 ## $ H0321030: num [1:4] 177.977 -0.513 59.607 5.425 ## $ H0400010: num [1:5] 2.84 73.77 -29.1 19.65 3 ## $ H0800010: num [1:5] 0.36 2275.6 -0.12 31.19 2.46 ``` ] .pull-right[ ```r ?RunModel.GRiwrmInputsModel ``` -- ```r OMnoRes <- RunModel(IMnoRes, RunOptions = ROnoRes, Param = ParamNoRes) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0321030... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0400010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0800010... ``` ] --- # Visualisation des simulations .pull-left[ ```r plot(OMnoRes, Q[I_Run,]) ``` ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-55-1.png)<!-- -->![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-55-2.png)<!-- -->![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-55-3.png)<!-- -->![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-55-4.png)<!-- --> ] -- .pull-right[ ```r Qm3sNoRes <- attr(OMnoRes, "Qm3s") plot(Qm3sNoRes[1:365,]) ``` ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-56-1.png)<!-- --> ] --- # Influence notable du réservoir sur le débit à l'aval .pull-left[ ```r plot(OMnoRes$H0800010, Q[I_Run, "H0800010"], which = "Regime") ``` ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-57-1.png)<!-- --> ] .pull-right[ Le régime est fortement influencé par le réservoir, le modèle ne peut le représenter correctement. Les six premiers mois de l'année correspondent à la période de remplissage du réservoir et les six derniers à la vidange pour le soutien d'étiage. Pour une modélisation correcte, il faut explicitement prendre en compte l'influence du réservoir. ] --- class: inverse, middle, center # Ajoutons l'influence des réservoirs... --- # Intégrer les influences dans le modèle Comment modéliser le lac réservoir Seine et ses connexions avec **airGRiwrm** ? .pull-left[
] -- .pull-right[
] --- # Intégrer les influences dans le modèle .pull-left[ Schéma du réseau...
] -- .pull-right[ ... et sa représentation dans le modèle
] --- # Intégrer les influences dans le modèle Ajout des noeuds de connexion au réservoir dans l'objet `GRiwrm` : -- ```r connections_reservoir <- data.frame( id = c("P_SEINE", "R_SEINE"), down = c("H0800010", "H0800010"), length = c(41, 11), model = NA, # Pas de modèle hydrologique = injection débit observé area = NA # Pas d'aire associée => Q fourni en m3/j ) griwrmInf <- rbind(griwrmNat, connections_reservoir) ``` |id |down | length|model | area| |:--------|:--------|------:|:-------------|-------:| |H0800010 |NA | NA|RunModel_GR4J | 3499.50| |H0400010 |H0800010 | 45.81|RunModel_GR4J | 2340.37| |H0210010 |H0400010 | 7.86|RunModel_GR4J | 1462.66| |H0321030 |H0400010 | 32.34|RunModel_GR4J | 548.93| |P_SEINE |H0800010 | 41.00|NA | NA| |R_SEINE |H0800010 | 11.00|NA | NA| --- # Intégrer les influences dans le modèle .pull-left[ Ajout des débits observés aux connexions : ```r summary(Q_reservoir_seine) ``` ``` ## P_SEINE R_SEINE ## Min. :-146.000 Min. : 0.000 ## 1st Qu.: -14.000 1st Qu.: 0.000 ## Median : -0.500 Median : 6.000 ## Mean : -8.533 Mean : 8.977 ## 3rd Qu.: 0.000 3rd Qu.:18.000 ## Max. : 0.010 Max. :45.000 ``` Les débits sont en m3/s... Les débits prélevés aux prises sont négatifs (débit retiré du modèle) et les débits restitués sont positifs (débit ajouté au modèle) ] .pull-right[ ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-65-1.png)<!-- --> ] --- # `CreateInputsModel` avec des débits influencés ## Conversion et fusion Par défaut, les débits des noeuds sans surface (`area=NA`) sont en m<sup>3</sup> par pas de temps: ```r Qinf <- cbind(Q, Q_reservoir_seine * 86400) IMinf <- CreateInputsModel(griwrmInf, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = Qinf) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0210010... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0321030... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0400010... ``` ``` ## Warning in airGR::CreateInputsModel(FUN_MOD = x, ...): 'Qupstream' contains NA ## values: model outputs will contain NAs ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin H0800010... ``` ``` ## Warning in airGR::CreateInputsModel(FUN_MOD = x, ...): 'Qupstream' contains NA ## values: model outputs will contain NAs ``` --- # Calage en débit influencé Appelons `CreateRunOptions`, `CreateInputsCrit` et `CreateCalibOptions`... -- ```r ROinf <- CreateRunOptions( InputsModel = IMinf, IndPeriod_WarmUp = 1:365, IndPeriod_Run = I_Run ) ICinf <- CreateInputsCrit( InputsModel = IMinf, FUN_CRIT = airGR::ErrorCrit_KGE2, RunOptions = ROinf, Obs = Qinf[I_Run,], transfo = "sqrt" ) COinf <- CreateCalibOptions(IMinf) ``` --- # Calage en débit influencé ## Calibration ```r OCinf <- suppressWarnings( Calibration(InputsModel = IMinf, RunOptions = ROinf, InputsCrit = ICinf, CalibOptions = COinf) ) ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 432.681, -0.020, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8845 ## Steepest-descent local search in progress ## Calibration completed (44 iterations, 438 runs) ## Param = 275.925, 0.219, 99.248, 5.068 ## Crit. KGE2[sqrt(Q)] = 0.9339 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0321030... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 169.017, -0.649, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.9089 ## Steepest-descent local search in progress ## Calibration completed (25 iterations, 266 runs) ## Param = 177.977, -0.513, 59.607, 5.425 ## Crit. KGE2[sqrt(Q)] = 0.9380 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0400010... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 247.151, -2.376, 20.697, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8828 ## Steepest-descent local search in progress ## Calibration completed (57 iterations, 798 runs) ## Param = 2.839, 73.774, -29.099, 19.651, 3.003 ## Crit. KGE2[sqrt(Q)] = 0.8968 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0800010... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 169.017, -0.649, 42.098, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.9261 ## Steepest-descent local search in progress ## Calibration completed (44 iterations, 667 runs) ## Param = 3.038, 176.367, -0.435, 24.229, 10.789 ## Crit. KGE2[sqrt(Q)] = 0.9356 ``` --- # Simulation en débit influencé Extraction des paramètres et lancement de `RunModel` : -- .pull-left[ ```r ParamInf <- lapply(OCinf, "[[", "ParamFinalR") OMinf <- RunModel(IMinf, RunOptions = ROinf, Param = ParamInf) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0321030... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0400010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0800010... ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 32 ## time steps with negative flow, set to zero. ``` Affichage du graphique *airGR* pour le bassin aval : -- ```r plot(OMinf$H0800010, Qinf[I_Run, "H0800010"]) ``` ] .pull-right[ ``` ## Warning in plot.OutputsModel(OMinf$H0800010, Qinf[I_Run, "H0800010"]): zeroes ## detected in 'Qsim': some plots in the log space will not be created using all ## time-steps ``` ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-71-1.png)<!-- --> ] --- class: inverse, middle, center # Régularisation des paramètres des modèles --- # Résultat du calage en débit influencé | | velocity| X1| X2| X3| X4| KGE2| |:--------|--------:|--------:|-----------:|--------:|---------:|---------:| |H0210010 | NA| 275.9254| 0.2192019| 99.24768| 5.068402| 0.9339461| |H0321030 | NA| 177.9773| -0.5131736| 59.60650| 5.424561| 0.9380038| |H0400010 | 2.839060| 73.7741| -29.0990363| 19.65126| 3.002584| 0.8967615| |H0800010 | 3.038152| 176.3665| -0.4351606| 24.22887| 10.788985| 0.9356111| Problème d'identification du paramètre `celerity` sur les bassins H0400010 et H0800010 : - des valeurs supérieures à 2 m/s sont peu réalistes - le pas de temps journalier ne permet pas d'identifier correctement ce paramètre sur des distances aussi courtes (paramètre peu sensible) --- # Régularisation avec la méthode de *Lavenne* Il faut définir des paramètres "a priori" à partir des paramètres d'un autre bassin ```r ?CreateInputsCrit.GRiwrmInputsModel ``` Dans *airGRiwrm*, on définit un bassin amont "donneur" pour chaque sous-bassin aval : ```r ICreg <- CreateInputsCrit( InputsModel = IMinf, FUN_CRIT = airGR::ErrorCrit_KGE2, RunOptions = ROinf, Obs = Qinf[I_Run,], transfo = "sqrt", * AprioriIds = c("H0400010" = "H0210010", "H0800010" = "H0400010"), * AprCelerity = 1 ) ``` Pour un sous-bassin aval "receveur" d'un bassin amont GR seul (pas de paramètre `celerity`), on définit une valeur par défaut `AprCelerity = 1` --- # Calage avec régularisation des paramètres ```r OCreg <- suppressWarnings( Calibration(InputsModel = IMinf, RunOptions = ROinf, InputsCrit = ICreg, CalibOptions = COinf) ) ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 432.681, -0.020, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.8845 ## Steepest-descent local search in progress ## Calibration completed (44 iterations, 438 runs) ## Param = 275.925, 0.219, 99.248, 5.068 ## Crit. KGE2[sqrt(Q)] = 0.9339 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0321030... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 169.017, -0.649, 83.096, 2.384 ## Crit. KGE2[sqrt(Q)] = 0.9089 ## Steepest-descent local search in progress ## Calibration completed (25 iterations, 266 runs) ## Param = 177.977, -0.513, 59.607, 5.425 ## Crit. KGE2[sqrt(Q)] = 0.9380 ## Calibration.GRiwrmInputsModel: Treating sub-basin H0400010... ## Crit. KGE2[sqrt(Q)] = 0.9339 ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9346 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 1.0073 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 1.0052 ## ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 247.151, -2.376, 20.697, 2.384 ## Crit. Composite = 0.8224 ## Steepest-descent local search in progress ## Calibration completed (38 iterations, 604 runs) ## Param = 0.990, 262.434, -5.810, 17.637, 4.765 ## Crit. Composite = 0.8845 ## Formula: sum(0.86 * KGE2[sqrt(Q)], 0.14 * GAPX[ParamT]) ## Calibration.GRiwrmInputsModel: Treating sub-basin H0800010... ## Crit. KGE2[sqrt(Q)] = 0.8919 ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9003 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 0.9989 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 1.0419 ## ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 169.017, -2.376, 83.096, 2.384 ## Crit. Composite = 0.8681 ## Steepest-descent local search in progress ## Calibration completed (39 iterations, 613 runs) ## Param = 0.830, 159.174, -1.992, 57.974, 4.375 ## Crit. Composite = 0.9307 ## Formula: sum(0.86 * KGE2[sqrt(Q)], 0.14 * GAPX[ParamT]) ``` --- # Résultat du calage avec régularisation Il faut calculer l'indicateur de performance KGE2 : ```r ParamReg <- lapply(OCreg, "[[", "ParamFinalR") OMreg <- RunModel(IMinf, ROinf, ParamReg) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0321030... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0400010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0800010... ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 34 ## time steps with negative flow, set to zero. ``` ```r KGEreg <- sapply(names(OMreg), function(id) { ErrorCrit_KGE2(ICinf[[id]], OMreg[[id]])$CritValue }) ``` ``` ## Crit. KGE2[sqrt(Q)] = 0.9339 ``` ``` ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9346 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 1.0073 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 1.0052 ``` ``` ## Crit. KGE2[sqrt(Q)] = 0.9380 ``` ``` ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9401 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 0.9854 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 0.9938 ``` ``` ## Crit. KGE2[sqrt(Q)] = 0.8919 ``` ``` ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9003 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 0.9989 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 1.0419 ``` ``` ## Crit. KGE2[sqrt(Q)] = 0.9325 ``` ``` ## SubCrit. KGE2[sqrt(Q)] cor(sim, obs, "pearson") = 0.9349 ## SubCrit. KGE2[sqrt(Q)] cv(sim)/cv(obs) = 0.9925 ## SubCrit. KGE2[sqrt(Q)] mean(sim)/mean(obs) = 0.9841 ``` --- # Résultat du calage avec régularisation | | velocity| X1| X2| X3| X4| KGE2 sans reg.| KGE2 avec reg.| |:--------|--------:|--------:|----------:|--------:|--------:|--------------:|--------------:| |H0210010 | NA| 275.9254| 0.2192019| 99.24768| 5.068402| 0.9339461| 0.9339461| |H0321030 | NA| 177.9773| -0.5131736| 59.60650| 5.424561| 0.9380038| 0.9380038| |H0400010 | 0.99| 262.4341| -5.8096883| 17.63702| 4.765015| 0.8967615| 0.8919017| |H0800010 | 0.83| 159.1743| -1.9918840| 57.97431| 4.374625| 0.9356111| 0.9325424| Le paramètre `velocity` a désormais une valeur "plus physique". La régularisation permet ainsi d'avoir des paramètres plus robustes pour une dégradation mineure de la performance. --- # Question bonus: comment simuler des débits naturalisés ? -- .pull-left[ ```r OMnat <- RunModel( IMnoRes, RunOptions = ROnoRes, Param = ParamReg ) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0210010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0321030... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0400010... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin H0800010... ``` ```r plot(OMnat$H0800010, Qinf[I_Run, "H0800010"], which = "Regime") ``` ] .pull-right[ ![](Formation_airGRiwrm_files/figure-html/unnamed-chunk-79-1.png)<!-- --> ] --- class: inverse, middle, center # La régulation automatique centralisée --- # La régulation automatique centralisée ## Exemples de régulations - Simuler le fonctionnement de réservoirs avec leurs règles de gestion (Suivi d'une courbe guide de remplissage, seuils de débits à respecter sur des points du réseau...) - Simuler des restrictions de prélèvements (AEP, irrigation...) pour la semaine en fonction du débit moyen simulé à l'aval la semaine précédente ## Principe de fonctionnement Un superviseur arrête la simulation à intervalles réguliers et effectue les opérations suivantes: - lit des "mesures" de débit simulés sur des noeud(s) au(x) pas de temps précédent(s) - calcule des commandes à partir d'un algorithmes fourni par l'utilisateur - applique les commandes sur des noeuds d'injection directe de débit jusqu'au prochain pas de supervision ## Exemple d'application Simulation d'une restriction de prélèvement de deux prises d'irrigation par "tour d'eau" : http://airgriwrm.g-eau.fr/dev/articles/V04_Closed-loop_regulated_withdrawal.html --- class: inverse, middle, center # Feuilles de route des modèles SD dans **airGR** et **airGRiwrm** --- class: hide_logo # Feuilles de route des modèles SD dans **airGR** et **airGRiwrm** ## airGR - Etendre le routage à d'autres fonctions que `Lag` - Simplifier le code en utilisant les classes S3 pour chaîner naturellement `CemaNeige + GR + SD` - Intégration facilitée des modèles **airGRplus** ## airGRiwrm - Capacité d'ajouter des stations non gaugées calées sur la station aval - Compatibilité de **airGRiwrm** avec **airGRPlus** (Après le passage des modèles **airGR** en classes S3) - Utilisation de modèles pour faire de la régulation locale (On remplace un GR par un modèle qui représente une influence anthropique): https://gitlab.irstea.fr/in-wop/airGRiwrm/-/issues/29 --- class: inverse, middle, center # Merci pour votre attention 😃