class: center, middle, inverse, title-slide # Formation à <strong>airGRiwrm</strong> ## Gestion intégrée de la ressource en eau avec <strong>airGR</strong> ### David Dorchies ### UMR G-EAU ### Montpellier / Paris, le 20 avril 2021 --- <div> <style type="text/css">.xaringan-extra-logo { width: 110px; height: 128px; z-index: 0; background-image: url(img/logo_airGRiwrm.png); background-size: contain; background-repeat: no-repeat; position: absolute; top:1em;right:1em; } </style> <script>(function () { let tries = 0 function addLogo () { if (typeof slideshow === 'undefined') { tries += 1 if (tries < 10) { setTimeout(addLogo, 100) } } else { document.querySelectorAll('.remark-slide-content:not(.title-slide):not(.inverse):not(.hide_logo)') .forEach(function (slide) { const logo = document.createElement('div') logo.classList = 'xaringan-extra-logo' logo.href = null slide.appendChild(logo) }) } } document.addEventListener('DOMContentLoaded', addLogo) })()</script> </div> # Programme de la formation - 14h05: Tour de table... - 14h10: Installation du package **airGRiwrm** - 14h20: Le jeu de données de la formation - 14h30: Le modèle semi-distribué dans **airGR** - 15h00: **airGRiwrm**: une extension de **airGR** pour un réseau de modèles semi-distribués - 15h10: Créer un réseau de modèles semi-distribués avec **airGRiwrm** - 15h30: Calage d'un réseau de modèles hydrologiques semi-distribués - 15h50: Ajoutons des perturbations humaines... - 16h10: Calage d'un réseau de modèles avec des débits influencés - 16h30: Simulation d'un réservoir régulé - 16h50: Feuille de route des modèles SD dans **airGR** et **airGRiwrm** --- class: inverse, middle, center # Tour de table... ### Votre expérience avec *airGR* ### Vos besoins en modèle hydrologique semi-distribué ? --- class: inverse, middle, center # Installation des packages *airGR* et *airGRiwrm* --- # Installation des packages *airGR* et *airGRiwrm* ### Sous Linux ou sous Windows avec le package `devtools` et `Rtools` installés: ```r remotes::install_gitlab("in-wop/airgriwrm", host = "gitlab.irstea.fr") ``` ### Sinon, installation des packages locaux sous Windows Installation de **airGR** version "dev": ```r install.packages("packages/airGR_1.6.10.9000.zip", repo = NULL) ``` et **airGRiwrm** ```r install.packages("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 "Marne" Données disponibles au pas de temps journalier: - Données météo des sous-bassins (P, ETP) - Débit naturalisés ] .pull-right[ <div id="htmlwidget-c744053093bb3c895dfb" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-c744053093bb3c895dfb">{"x":{"options":{"minZoom":1,"maxZoom":52,"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}},"preferCanvas":false,"bounceAtZoomLimits":false,"maxBounds":[[[-90,-370]],[[90,370]]]},"calls":[{"method":"addProviderTiles","args":["CartoDB.Positron","CartoDB.Positron","CartoDB.Positron",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["CartoDB.DarkMatter","CartoDB.DarkMatter","CartoDB.DarkMatter",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenStreetMap","OpenStreetMap","OpenStreetMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["Esri.WorldImagery","Esri.WorldImagery","Esri.WorldImagery",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenTopoMap","OpenTopoMap","OpenTopoMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"createMapPane","args":["point",440]},{"method":"addFlatGeoBuf","args":["locations_sf","locations_sf",null,true,"mvFeatureId",{"radius":6,"stroke":true,"color":"#333333","weight":1,"opacity":0.9,"fill":true,"fillColor":"#6666FF","fillOpacity":0.6},{"className":"","pane":"point"},"mapview-popup",{"radius":{"to":[3,15],"from":[3,15]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addScaleBar","args":[{"maxWidth":100,"metric":true,"imperial":true,"updateWhenIdle":true,"position":"bottomleft"}]},{"method":"addHomeButton","args":[4.35234940338672,48.5158749743783,4.92973814401406,48.9533034808677,"locations_sf","Zoom to locations_sf","<strong> locations_sf <\/strong>","bottomright"]},{"method":"addLayersControl","args":[["CartoDB.Positron","CartoDB.DarkMatter","OpenStreetMap","Esri.WorldImagery","OpenTopoMap"],"locations_sf",{"collapsed":true,"autoZIndex":true,"position":"topleft"}]},{"method":"addLegend","args":[{"colors":["#6666FF"],"labels":["locations_sf"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"","extra":null,"layerId":null,"className":"info legend","group":"locations_sf"}]},{"method":"addMarkers","args":[[48.6364152102786,48.5158749743783,48.7460418340255,48.9533034808677],[4.92895392640666,4.92973814401406,4.62711339136132,4.35234940338672],{"iconUrl":{"data":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAEAAAABCAQAAAC1HAwCAAAAC0lEQVR4nGP6zwAAAgcBApocMXEAAAAASUVORK5CYII=","index":0},"iconWidth":1,"iconHeight":1},null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},null,null,null,null,["STDIZ_04","LOUVE_19","VITRY_25","CHALO_21"],{"permanent":true,"direction":"top","textOnly":true,"offset":[0,20]},null]}],"fitBounds":[48.5158749743783,4.35234940338672,48.9533034808677,4.92973814401406,[]],"limits":{"lat":[48.5158749743783,48.9533034808677],"lng":[4.35234940338672,4.92973814401406]}},"evals":[],"jsHooks":{"render":[{"code":"function(el, x, data) {\n return (\n function(el, x, data) {\n // get the leaflet map\n var map = this; //HTMLWidgets.find('#' + el.id);\n // we need a new div element because we have to handle\n // the mouseover output separately\n // debugger;\n function addElement () {\n // generate new div Element\n var newDiv = $(document.createElement('div'));\n // append at end of leaflet htmlwidget container\n $(el).append(newDiv);\n //provide ID and style\n newDiv.addClass('lnlt');\n newDiv.css({\n 'position': 'relative',\n 'bottomleft': '0px',\n 'background-color': 'rgba(255, 255, 255, 0.7)',\n 'box-shadow': '0 0 2px #bbb',\n 'background-clip': 'padding-box',\n 'margin': '0',\n 'padding-left': '5px',\n 'color': '#333',\n 'font': '9px/1.5 \"Helvetica Neue\", Arial, Helvetica, sans-serif',\n 'z-index': '700',\n });\n return newDiv;\n }\n\n\n // check for already existing lnlt class to not duplicate\n var lnlt = $(el).find('.lnlt');\n\n if(!lnlt.length) {\n lnlt = addElement();\n\n // grab the special div we generated in the beginning\n // and put the mousmove output there\n\n map.on('mousemove', function (e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() +\n ' | x: ' + L.CRS.EPSG3857.project(e.latlng).x.toFixed(0) +\n ' | y: ' + L.CRS.EPSG3857.project(e.latlng).y.toFixed(0) +\n ' | epsg: 3857 ' +\n ' | proj4: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ');\n } else {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n }\n });\n\n // remove the lnlt div when mouse leaves map\n map.on('mouseout', function (e) {\n var strip = document.querySelector('.lnlt');\n if( strip !==null) strip.remove();\n });\n\n };\n\n //$(el).keypress(67, function(e) {\n map.on('preclick', function(e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n var txt = document.querySelector('.lnlt').textContent;\n console.log(txt);\n //txt.innerText.focus();\n //txt.select();\n setClipboardText('\"' + txt + '\"');\n }\n });\n\n }\n ).call(this.getMap(), el, x, data);\n}","data":null},{"code":"function(el, x, data) {\n return (function(el,x,data){\n var map = this;\n\n map.on('keypress', function(e) {\n console.log(e.originalEvent.code);\n var key = e.originalEvent.code;\n if (key === 'KeyE') {\n var bb = this.getBounds();\n var txt = JSON.stringify(bb);\n console.log(txt);\n\n setClipboardText('\\'' + txt + '\\'');\n }\n })\n }).call(this.getMap(), el, x, data);\n}","data":null}]}}</script> ] --- # 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 > rm(list = ls()) > load("data/marne_naturalised.RData") > ls() [1] "DatesR" "E" "marne_nodes" "P" "Q" ``` - `DatesR` (*POSIXct*): Dates des pas de temps d'observation - `E`, `P`, `Q` (*matrix*): Données hydroclimatiques - `marne_nodes` (*data.frame*): description du réseau --- ## Cas d'étude: description du réseau ```r marne_nodes ``` | |id_sgl | area|id_aval | distance_aval|description | |:--|:--------|-------:|:--------|-------------:|:----------------------------| |4 |STDIZ_04 | 2347.53|CHALO_21 | 85.570|La Marne à Saint-Dizier | |19 |LOUVE_19 | 461.74|CHALO_21 | 86.165|La Blaise à Louvemont | |25 |VITRY_25 | 2109.14|CHALO_21 | 38.047|La Saulx à Vitry-en-Perthois | |21 |CHALO_21 | 6291.55|NA | NA|La Marne à Châlons-sur-Marne | --- # Cas d'étude: description des obs ```r head(DatesR, 3) ``` ``` ## [1] "1994-08-01 CEST" "1994-08-02 CEST" "1994-08-03 CEST" ``` ```r tail(DatesR, 3) ``` ``` ## [1] "2008-07-29 CEST" "2008-07-30 CEST" "2008-07-31 CEST" ``` ```r head(Q) ``` ``` ## STDIZ_04 LOUVE_19 VITRY_25 CHALO_21 ## 13150 0.4195729 0.2058301 0.1966299 0.2032440 ## 13151 0.4011706 0.1871183 0.2253051 0.1991242 ## 13152 0.3165199 0.1871183 0.1884370 0.2183500 ## 13153 0.3054785 0.1684065 0.1761476 0.2526818 ## 13154 0.2723543 0.1871183 0.1597618 0.2650412 ## 13155 0.2502716 0.2058301 0.1474724 0.2499352 ``` --- 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[   ] --- # Modèle SD dans *airGR* ## La Marne à Châlons-sur-Marne et ses 3 débits amont .pull-left[ <div id="htmlwidget-00e8b26ce2c45f4dc5da" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-00e8b26ce2c45f4dc5da">{"x":{"options":{"minZoom":1,"maxZoom":52,"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}},"preferCanvas":false,"bounceAtZoomLimits":false,"maxBounds":[[[-90,-370]],[[90,370]]]},"calls":[{"method":"addProviderTiles","args":["CartoDB.Positron","CartoDB.Positron","CartoDB.Positron",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["CartoDB.DarkMatter","CartoDB.DarkMatter","CartoDB.DarkMatter",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenStreetMap","OpenStreetMap","OpenStreetMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["Esri.WorldImagery","Esri.WorldImagery","Esri.WorldImagery",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenTopoMap","OpenTopoMap","OpenTopoMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"createMapPane","args":["point",440]},{"method":"addFlatGeoBuf","args":["locations_sf","locations_sf",null,true,"mvFeatureId",{"radius":6,"stroke":true,"color":"#333333","weight":1,"opacity":0.9,"fill":true,"fillColor":"#6666FF","fillOpacity":0.6},{"className":"","pane":"point"},"mapview-popup",{"radius":{"to":[3,15],"from":[3,15]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addScaleBar","args":[{"maxWidth":100,"metric":true,"imperial":true,"updateWhenIdle":true,"position":"bottomleft"}]},{"method":"addHomeButton","args":[4.35234940338672,48.5158749743783,4.92973814401406,48.9533034808677,"locations_sf","Zoom to locations_sf","<strong> locations_sf <\/strong>","bottomright"]},{"method":"addLayersControl","args":[["CartoDB.Positron","CartoDB.DarkMatter","OpenStreetMap","Esri.WorldImagery","OpenTopoMap"],"locations_sf",{"collapsed":true,"autoZIndex":true,"position":"topleft"}]},{"method":"addLegend","args":[{"colors":["#6666FF"],"labels":["locations_sf"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"","extra":null,"layerId":null,"className":"info legend","group":"locations_sf"}]}],"fitBounds":[48.5158749743783,4.35234940338672,48.9533034808677,4.92973814401406,[]]},"evals":[],"jsHooks":{"render":[{"code":"function(el, x, data) {\n return (\n function(el, x, data) {\n // get the leaflet map\n var map = this; //HTMLWidgets.find('#' + el.id);\n // we need a new div element because we have to handle\n // the mouseover output separately\n // debugger;\n function addElement () {\n // generate new div Element\n var newDiv = $(document.createElement('div'));\n // append at end of leaflet htmlwidget container\n $(el).append(newDiv);\n //provide ID and style\n newDiv.addClass('lnlt');\n newDiv.css({\n 'position': 'relative',\n 'bottomleft': '0px',\n 'background-color': 'rgba(255, 255, 255, 0.7)',\n 'box-shadow': '0 0 2px #bbb',\n 'background-clip': 'padding-box',\n 'margin': '0',\n 'padding-left': '5px',\n 'color': '#333',\n 'font': '9px/1.5 \"Helvetica Neue\", Arial, Helvetica, sans-serif',\n 'z-index': '700',\n });\n return newDiv;\n }\n\n\n // check for already existing lnlt class to not duplicate\n var lnlt = $(el).find('.lnlt');\n\n if(!lnlt.length) {\n lnlt = addElement();\n\n // grab the special div we generated in the beginning\n // and put the mousmove output there\n\n map.on('mousemove', function (e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() +\n ' | x: ' + L.CRS.EPSG3857.project(e.latlng).x.toFixed(0) +\n ' | y: ' + L.CRS.EPSG3857.project(e.latlng).y.toFixed(0) +\n ' | epsg: 3857 ' +\n ' | proj4: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ');\n } else {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n }\n });\n\n // remove the lnlt div when mouse leaves map\n map.on('mouseout', function (e) {\n var strip = document.querySelector('.lnlt');\n if( strip !==null) strip.remove();\n });\n\n };\n\n //$(el).keypress(67, function(e) {\n map.on('preclick', function(e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n var txt = document.querySelector('.lnlt').textContent;\n console.log(txt);\n //txt.innerText.focus();\n //txt.select();\n setClipboardText('\"' + txt + '\"');\n }\n });\n\n }\n ).call(this.getMap(), el, x, data);\n}","data":null},{"code":"function(el, x, data) {\n return (function(el,x,data){\n var map = this;\n\n map.on('keypress', function(e) {\n console.log(e.originalEvent.code);\n var key = e.originalEvent.code;\n if (key === 'KeyE') {\n var bb = this.getBounds();\n var txt = JSON.stringify(bb);\n console.log(txt);\n\n setClipboardText('\\'' + txt + '\\'');\n }\n })\n }).call(this.getMap(), el, x, data);\n}","data":null}]}}</script> ] .pull-right[ <div id="htmlwidget-b5a2c097beeb73e2c36b" style="width:504px;height:504px;" class="grViz html-widget"></div> <script type="application/json" data-for="htmlwidget-b5a2c097beeb73e2c36b">{"x":{"diagram":"\ndigraph MarneNoConnect {\nrankdir=BT\nedge [color = Navy]\n\nnode [shape = egg,\n margin=0,\n fontname = Helvetica,\n color = DarkGreen,\n style = filled,\n fillcolor = lightpink]\nCHALO_21 [fillcolor = palegreen]\nVITRY_25 -> CHALO_21 [label=\"38km\"]\nLOUVE_19 -> CHALO_21 [label=\"86km\"]\nSTDIZ_04 -> CHALO_21 [label=\"86km\"]\n}\n","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script> ] --- # Modèle SD dans *airGR* ## `CreateInputsModel`: ```r library(airGR) ?airGR::CreateInputsModel ``` -- ```r BasinAreas <- marne_nodes$area BasinAreas[4] <- BasinAreas[4] - sum(BasinAreas[1:3]) IMairGR <- CreateInputsModel( FUN_MOD = RunModel_GR4J, DatesR = DatesR, Precip = P[, "CHALO_21"], PotEvap = E[, "CHALO_21"], Qupstream = Q[, c("STDIZ_04", "LOUVE_19", "VITRY_25")], LengthHydro = marne_nodes$distance_aval[1:3], BasinAreas = BasinAreas ) ``` --- # Modèle SD dans *airGR* ## `CreateInputsModel`: ```r str(IMairGR) ``` ``` ## List of 6 ## $ DatesR : POSIXlt[1:5114], format: "1994-08-01" "1994-08-02" ... ## $ Precip : num [1:5114] 0.6 0 0 0 0 3.8 16.4 3.2 7.8 5.5 ... ## $ PotEvap : num [1:5114] 3.8 4 4.4 4.7 4.6 4.4 4.3 3.8 4 3.6 ... ## $ Qupstream : num [1:5114, 1:3] 984960 941760 743040 717120 639360 ... ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : chr [1:5114] "13150" "13151" "13152" "13153" ... ## .. ..$ : chr [1:3] "STDIZ_04" "LOUVE_19" "VITRY_25" ## $ LengthHydro: num [1:3] 85.6 86.2 38 ## $ BasinAreas : num [1:4] 2348 462 2109 1373 ## - 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 ``` -- ```r ICairGR <- CreateInputsCrit( FUN_CRIT = ErrorCrit_NSE, InputsModel = IMairGR, RunOptions = ROairGR, Obs = Q[I_Run, "CHALO_21"] ) ``` --- # Modèle SD dans *airGR* ## `CreateCalibOptions`: ```r ?airGR::CreateCalibOptions ``` -- ```r COairGR <- CreateCalibOptions( FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel, IsSD = TRUE ) ``` --- # Modèle SD dans *airGR* ```r ?airGR::Calibration ``` -- ```r OCairGR <- 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 = 11.250, 432.681, -2.376, 83.096, 2.384 ## Crit. NSE[Q] = 0.8609 ## Steepest-descent local search in progress ## Calibration completed (110 iterations, 1337 runs) ## Param = 0.446, 80.963, -30.754, 1039.537, 10.678 ## Crit. NSE[Q] = 0.9860 ``` --- # Modèle SD dans *airGR* .pull-left[ ## RunModel ```r OMairGR <- RunModel( InputsModel = IMairGR, RunOptions = ROairGR, Param = OCairGR$ParamFinalR, FUN_MOD = RunModel_GR4J ) ``` ```r plot(OMairGR, Q[I_Run, "CHALO_21"]) ``` ] .pull-right[ <!-- --> ] --- 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 .pull-left[ ## .center[Disclaimer]  ] .pull-right[ 🚧 Le package est en développement 🚧 Cette formation est l'occasion: - débusquer des bugs - lever des ambiguïtés - remettre en cause et simplifier des procédures - définir une feuille de route pour la suite des développements ] --- class: hide_logo # *airGRiwrm*: extension de *airGR* pour un réseau de modèles semi-distribués **airGRiwrm** étend 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* --- # Créer un réseau de modèles semi-distribués avec *airGRiwrm* L'extension à tout le réseau est défini à partir d'un objet décrivant ce dernier: ```r ?GRiwrm ``` -- Définition des noeuds modélisés par un modèle hydrologique: ```r marne_nodes$model <- "RunModel_GR4J" ``` -- Création de l'objet `GRiwrm` (correspondance des noms de colonnes) ```r griwrmNat <- GRiwrm(marne_nodes, cols = list(id = "id_sgl", down = "id_aval", length = "distance_aval")) ``` --- # Créer un réseau de modèles semi-distribués avec **airGRiwrm** ```r DiagramGRiwrm(griwrmNat) ``` .pull-left[ <div id="htmlwidget-9eb4a5ce4abbc617b288" style="width:100%;height:100%;" class="DiagrammeR html-widget"></div> <script type="application/json" data-for="htmlwidget-9eb4a5ce4abbc617b288">{"x":{"diagram":"graph LR\nSTDIZ_04 -->| 86 km| CHALO_21\nLOUVE_19 -->| 86 km| CHALO_21\nVITRY_25 -->| 38 km| CHALO_21\nstyle CHALO_21 fill:#cfc\nstyle fill:#fcc"},"evals":[],"jsHooks":[]}</script> ] .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) ] --- # `CreateInputsModel` dans *airGRiwrm* ```r ?CreateInputsModel ``` -- ```r IMnat <- CreateInputsModel( griwrmNat, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = Q ) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin STDIZ_04... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin LOUVE_19... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin VITRY_25... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin CHALO_21... ``` Création d'une liste d'InputsModel : 1 élément = 1 modèle *GR ( + Lag )* --- # L'objet `GRiwrmInputsModel` ```r str(IMnat) ``` ``` ## List of 4 ## $ STDIZ_04:List of 7 ## ..$ DatesR : POSIXlt[1:5114], format: "1994-08-01" "1994-08-02" ... ## ..$ Precip : num [1:5114] 0.7 0 0 0 0 6.1 17.1 2.2 10.8 8.6 ... ## ..$ PotEvap : num [1:5114] 3.8 4 4.4 4.7 4.7 4.4 4.3 3.8 4.1 3.4 ... ## ..$ id : chr "STDIZ_04" ## ..$ down : chr "CHALO_21" ## ..$ BasinAreas: num 2348 ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR" ## $ LOUVE_19:List of 7 ## ..$ DatesR : POSIXlt[1:5114], format: "1994-08-01" "1994-08-02" ... ## ..$ Precip : num [1:5114] 0.4 0 0 0 0.1 4.7 23.6 3.1 10.8 9.7 ... ## ..$ PotEvap : num [1:5114] 3.8 4 4.4 4.7 4.7 4.5 4.3 3.8 4.1 3.5 ... ## ..$ id : chr "LOUVE_19" ## ..$ down : chr "CHALO_21" ## ..$ BasinAreas: num 462 ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR" ## $ VITRY_25:List of 7 ## ..$ DatesR : POSIXlt[1:5114], format: "1994-08-01" "1994-08-02" ... ## ..$ Precip : num [1:5114] 0.9 0 0 0 0 1.6 17 6.1 6.9 4.1 ... ## ..$ PotEvap : num [1:5114] 3.8 3.9 4.3 4.6 4.6 4.4 4.2 3.8 3.9 3.5 ... ## ..$ id : chr "VITRY_25" ## ..$ down : chr "CHALO_21" ## ..$ BasinAreas: num 2109 ## ..$ FUN_MOD : chr "RunModel_GR4J" ## ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR" ## $ CHALO_21:List of 11 ## ..$ DatesR : POSIXlt[1:5114], format: "1994-08-01" "1994-08-02" ... ## ..$ Precip : num [1:5114] 0.6 0 0 0 0 3.8 16.4 3.2 7.8 5.5 ... ## ..$ PotEvap : num [1:5114] 3.8 4 4.4 4.7 4.6 4.4 4.3 3.8 4 3.6 ... ## ..$ Qupstream : num [1:5114, 1:3] 984960 941760 743040 717120 639360 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:5114] "13150" "13151" "13152" "13153" ... ## .. .. ..$ : chr [1:3] "STDIZ_04" "LOUVE_19" "VITRY_25" ## ..$ LengthHydro : Named num [1:3] 85.6 86.2 38 ## .. ..- attr(*, "names")= chr [1:3] "STDIZ_04" "LOUVE_19" "VITRY_25" ## ..$ BasinAreas : Named num [1:4] 2348 462 2109 1373 ## .. ..- attr(*, "names")= chr [1:4] "STDIZ_04" "LOUVE_19" "VITRY_25" "CHALO_21" ## ..$ id : chr "CHALO_21" ## ..$ down : chr NA ## ..$ UpstreamNodes : chr [1:3] "STDIZ_04" "LOUVE_19" "VITRY_25" ## ..$ UpstreamIsRunoff: logi [1:3] TRUE TRUE 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] "STDIZ_04" "LOUVE_19" "VITRY_25" "CHALO_21" ## ..$ down : chr [1:4] "CHALO_21" "CHALO_21" "CHALO_21" NA ## ..$ length: num [1:4] 85.6 86.2 38 NA ## ..$ model : chr [1:4] "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" ## ..$ area : num [1:4] 2348 462 2109 6292 ## - attr(*, "TimeStep")= num 86400 ``` --- class: inverse, middle, center # Calage d'un réseau de modèles hydrologiques semi-distribués avec *airGRiwrm* --- # Calage avec *airGRiwrm* ```r ROnat <- CreateRunOptions( InputsModel = IMnat, IndPeriod_WarmUp = 1:365, IndPeriod_Run = I_Run ) ``` -- ```r ICnat <- CreateInputsCrit( InputsModel = IMnat, FUN_CRIT = airGR::ErrorCrit_NSE, RunOptions = ROnat, Obs = Q[I_Run,] ) ``` ```r COnat <- CreateCalibOptions(IMnat) ``` --- # `Calibration` dans *airGRiwrm* ```r OCnat <- Calibration(InputsModel = IMnat, RunOptions = ROnat, InputsCrit = ICnat, CalibOptions = COnat, useUpstreamQsim = TRUE) ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 247.151, -0.020, 83.096, 2.384 ## Crit. NSE[Q] = 0.8732 ## Steepest-descent local search in progress ## Calibration completed (26 iterations, 273 runs) ## Param = 208.513, -0.130, 74.440, 3.506 ## Crit. NSE[Q] = 0.9236 ## Calibration.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 247.151, -2.376, 83.096, 2.384 ## Crit. NSE[Q] = 0.8788 ## Steepest-descent local search in progress ## Calibration completed (17 iterations, 200 runs) ## Param = 179.469, -2.703, 122.732, 2.267 ## Crit. NSE[Q] = 0.8923 ## Calibration.GRiwrmInputsModel: Treating sub-basin VITRY_25... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 432.681, -0.649, 83.096, 2.384 ## Crit. NSE[Q] = 0.8258 ## Steepest-descent local search in progress ## Calibration completed (26 iterations, 274 runs) ## Param = 257.238, -1.698, 123.965, 4.238 ## Crit. NSE[Q] = 0.9312 ## Calibration.GRiwrmInputsModel: Treating sub-basin CHALO_21... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 432.681, -2.376, 83.096, 2.384 ## Crit. NSE[Q] = 0.8372 ## Steepest-descent local search in progress ## Calibration completed (56 iterations, 788 runs) ## Param = 0.432, 678.901, -3.290, 71.127, 10.000 ## Crit. NSE[Q] = 0.9524 ``` --- # `Calibration` dans *airGRiwrm* ```r str(OCnat) ``` ``` ## List of 4 ## $ STDIZ_04:List of 9 ## ..$ ParamFinalR : num [1:4] 208.51 -0.13 74.44 3.51 ## ..$ CritFinal : num 0.924 ## ..$ NIter : num 26 ## ..$ NRuns : num 273 ## ..$ HistParamR : num [1:26, 1:4] 247 247 247 247 247 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:4] "Param1" "Param2" "Param3" "Param4" ## ..$ HistCrit : num [1:26, 1] 0.873 0.902 0.914 0.914 0.92 ... ## ..$ MatBoolCrit : logi [1:4749, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "NSE[Q]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ LOUVE_19:List of 9 ## ..$ ParamFinalR : num [1:4] 179.47 -2.7 122.73 2.27 ## ..$ CritFinal : num 0.892 ## ..$ NIter : num 17 ## ..$ NRuns : num 200 ## ..$ HistParamR : num [1:17, 1:4] 247 247 247 179 179 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:4] "Param1" "Param2" "Param3" "Param4" ## ..$ HistCrit : num [1:17, 1] 0.879 0.879 0.881 0.89 0.89 ... ## ..$ MatBoolCrit : logi [1:4749, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "NSE[Q]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ VITRY_25:List of 9 ## ..$ ParamFinalR : num [1:4] 257.24 -1.7 123.97 4.24 ## ..$ CritFinal : num 0.931 ## ..$ NIter : num 26 ## ..$ NRuns : num 274 ## ..$ HistParamR : num [1:26, 1:4] 433 433 433 433 433 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:4] "Param1" "Param2" "Param3" "Param4" ## ..$ HistCrit : num [1:26, 1] 0.826 0.86 0.903 0.914 0.914 ... ## ..$ MatBoolCrit : logi [1:4749, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "NSE[Q]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## $ CHALO_21:List of 9 ## ..$ ParamFinalR : num [1:5] 0.432 678.901 -3.29 71.127 10 ## ..$ CritFinal : num 0.952 ## ..$ NIter : num 56 ## ..$ NRuns : num 788 ## ..$ HistParamR : num [1:56, 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:56, 1] 0.837 0.85 0.857 0.861 0.863 ... ## ..$ MatBoolCrit : logi [1:4749, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:2] "BoolCrit_Requested" "BoolCrit_Actual" ## ..$ CritName : chr "NSE[Q]" ## ..$ CritBestValue: num 1 ## ..- attr(*, "class")= chr [1:2] "OutputsCalib" "HBAN" ## - attr(*, "class")= chr [1:2] "list" "GRiwrmOutputsCalib" ``` --- # `RunModel` dans *airGRiwrm* ```r ParamNat <- sapply(names(OCnat), function(x) {OCnat[[x]]$Param}) str(ParamNat) ``` ``` ## List of 4 ## $ STDIZ_04: num [1:4] 208.51 -0.13 74.44 3.51 ## $ LOUVE_19: num [1:4] 179.47 -2.7 122.73 2.27 ## $ VITRY_25: num [1:4] 257.24 -1.7 123.97 4.24 ## $ CHALO_21: num [1:5] 0.432 678.901 -3.29 71.127 10 ``` ```r OMnat <- RunModel(IMnat, RunOptions = ROnat, Param = ParamNat) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin VITRY_25... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin CHALO_21... ``` --- # Visualisation des simulations .pull-left[ ```r plot(OMnat, Q[I_Run,]) ``` <!-- --><!-- --><!-- --><!-- --> ``` ## [[1]] ## NULL ## ## [[2]] ## NULL ## ## [[3]] ## NULL ## ## [[4]] ## NULL ``` ] .pull-right[ ```r Qm3sNat <- attr(OMnat, "Qm3s") plot(Qm3sNat[1:365,]) ``` <!-- --> ] --- class: inverse, middle, center # Ajoutons des perturbations humaines... --- # Intégrer les influences dans le modèle Comment modéliser le lac réservoir Marne et ses connexions avec **airGRiwrm** ? .pull-left[ <div id="htmlwidget-b6ef9361d583746319f6" style="width:504px;height:504px;" class="leaflet html-widget"></div> <script type="application/json" data-for="htmlwidget-b6ef9361d583746319f6">{"x":{"options":{"minZoom":1,"maxZoom":52,"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}},"preferCanvas":false,"bounceAtZoomLimits":false,"maxBounds":[[[-90,-370]],[[90,370]]]},"calls":[{"method":"addProviderTiles","args":["CartoDB.Positron","CartoDB.Positron","CartoDB.Positron",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["CartoDB.DarkMatter","CartoDB.DarkMatter","CartoDB.DarkMatter",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenStreetMap","OpenStreetMap","OpenStreetMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["Esri.WorldImagery","Esri.WorldImagery","Esri.WorldImagery",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenTopoMap","OpenTopoMap","OpenTopoMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"createMapPane","args":["point",440]},{"method":"addFlatGeoBuf","args":["locations_sf","locations_sf",null,true,"mvFeatureId",{"radius":6,"stroke":true,"color":"#333333","weight":1,"opacity":0.9,"fill":true,"fillColor":"#6666FF","fillOpacity":0.6},{"className":"","pane":"point"},"mapview-popup",{"radius":{"to":[3,15],"from":[3,15]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addScaleBar","args":[{"maxWidth":100,"metric":true,"imperial":true,"updateWhenIdle":true,"position":"bottomleft"}]},{"method":"addHomeButton","args":[4.35234940338672,48.5158749743783,4.92973814401406,48.9533034808677,"locations_sf","Zoom to locations_sf","<strong> locations_sf <\/strong>","bottomright"]},{"method":"addLayersControl","args":[["CartoDB.Positron","CartoDB.DarkMatter","OpenStreetMap","Esri.WorldImagery","OpenTopoMap"],"locations_sf",{"collapsed":true,"autoZIndex":true,"position":"topleft"}]},{"method":"addLegend","args":[{"colors":["#6666FF"],"labels":["locations_sf"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"","extra":null,"layerId":null,"className":"info legend","group":"locations_sf"}]},{"method":"addMarkers","args":[[48.6364152102786,48.5158749743783,48.7460418340255,48.9533034808677],[4.92895392640666,4.92973814401406,4.62711339136132,4.35234940338672],{"iconUrl":{"data":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAEAAAABCAQAAAC1HAwCAAAAC0lEQVR4nGP6zwAAAgcBApocMXEAAAAASUVORK5CYII=","index":0},"iconWidth":1,"iconHeight":1},null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},null,null,null,null,["STDIZ_04","LOUVE_19","VITRY_25","CHALO_21"],{"permanent":true,"direction":"top","textOnly":true,"offset":[0,20]},null]}],"fitBounds":[48.5158749743783,4.35234940338672,48.9533034808677,4.92973814401406,[]],"limits":{"lat":[48.5158749743783,48.9533034808677],"lng":[4.35234940338672,4.92973814401406]}},"evals":[],"jsHooks":{"render":[{"code":"function(el, x, data) {\n return (\n function(el, x, data) {\n // get the leaflet map\n var map = this; //HTMLWidgets.find('#' + el.id);\n // we need a new div element because we have to handle\n // the mouseover output separately\n // debugger;\n function addElement () {\n // generate new div Element\n var newDiv = $(document.createElement('div'));\n // append at end of leaflet htmlwidget container\n $(el).append(newDiv);\n //provide ID and style\n newDiv.addClass('lnlt');\n newDiv.css({\n 'position': 'relative',\n 'bottomleft': '0px',\n 'background-color': 'rgba(255, 255, 255, 0.7)',\n 'box-shadow': '0 0 2px #bbb',\n 'background-clip': 'padding-box',\n 'margin': '0',\n 'padding-left': '5px',\n 'color': '#333',\n 'font': '9px/1.5 \"Helvetica Neue\", Arial, Helvetica, sans-serif',\n 'z-index': '700',\n });\n return newDiv;\n }\n\n\n // check for already existing lnlt class to not duplicate\n var lnlt = $(el).find('.lnlt');\n\n if(!lnlt.length) {\n lnlt = addElement();\n\n // grab the special div we generated in the beginning\n // and put the mousmove output there\n\n map.on('mousemove', function (e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() +\n ' | x: ' + L.CRS.EPSG3857.project(e.latlng).x.toFixed(0) +\n ' | y: ' + L.CRS.EPSG3857.project(e.latlng).y.toFixed(0) +\n ' | epsg: 3857 ' +\n ' | proj4: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ');\n } else {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n }\n });\n\n // remove the lnlt div when mouse leaves map\n map.on('mouseout', function (e) {\n var strip = document.querySelector('.lnlt');\n if( strip !==null) strip.remove();\n });\n\n };\n\n //$(el).keypress(67, function(e) {\n map.on('preclick', function(e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n var txt = document.querySelector('.lnlt').textContent;\n console.log(txt);\n //txt.innerText.focus();\n //txt.select();\n setClipboardText('\"' + txt + '\"');\n }\n });\n\n }\n ).call(this.getMap(), el, x, data);\n}","data":null},{"code":"function(el, x, data) {\n return (function(el,x,data){\n var map = this;\n\n map.on('keypress', function(e) {\n console.log(e.originalEvent.code);\n var key = e.originalEvent.code;\n if (key === 'KeyE') {\n var bb = this.getBounds();\n var txt = JSON.stringify(bb);\n console.log(txt);\n\n setClipboardText('\\'' + txt + '\\'');\n }\n })\n }).call(this.getMap(), el, x, data);\n}","data":null}]}}</script> ] .pull-right[ <div id="htmlwidget-a43f4322b6231b129b45" style="width:100%;height:504px;" class="DiagrammeR html-widget"></div> <script type="application/json" data-for="htmlwidget-a43f4322b6231b129b45">{"x":{"diagram":"\ngraph BT\nLOUVE_19 -->|0.5km| MARNE_P28(Prise Blaise)\nMARNE_P23(Prise Marne) -->|3km| STDIZ_04\nMARNE_P23 --> M(Lac Marne)\nMARNE_P28 --> M\nMARNE_P28 -->|25.6km| MARNE_R25(Restitution)\nM --> MARNE_R25\nSTDIZ_04 -->|28.7km| MARNE_R25\nMARNE_R25 -->|56.9km| CHALO_21\nVITRY_25 -->|38km| CHALO_21\nstyle MARNE_P28 fill:#fcc\nstyle MARNE_P23 fill:#fcc\nstyle MARNE_R25 fill:#fcc\nstyle M fill:#aaf\n"},"evals":[],"jsHooks":[]}</script> ] --- # Intégrer les influences dans le modèle .pull-left[ <div id="htmlwidget-e8fb97cf5c7560cc2d9d" style="width:100%;height:504px;" class="DiagrammeR html-widget"></div> <script type="application/json" data-for="htmlwidget-e8fb97cf5c7560cc2d9d">{"x":{"diagram":"\ngraph BT\nLOUVE_19 -->|0.5km| MARNE_P28(P_BLAISE)\nMARNE_P23(P_MARNE) -->|3km| STDIZ_04\nMARNE_P23 --> M(Lac Marne)\nMARNE_P28 --> M\nMARNE_P28 -->|25.6km| MARNE_R25(R_MARNE)\nM --> MARNE_R25\nSTDIZ_04 -->|28.7km| MARNE_R25\nMARNE_R25 -->|56.9km| CHALO_21\nVITRY_25 -->|38km| CHALO_21\nstyle MARNE_P28 fill:#fcc\nstyle MARNE_P23 fill:#fcc\nstyle MARNE_R25 fill:#fcc\nstyle M fill:#aaf\n"},"evals":[],"jsHooks":[]}</script> ] -- .pull-right[ <div id="htmlwidget-861c90a3f4abf9260be0" style="width:100%;height:100%;" class="DiagrammeR html-widget"></div> <script type="application/json" data-for="htmlwidget-861c90a3f4abf9260be0">{"x":{"diagram":"graph LR\nSTDIZ_04 -->| 86 km| CHALO_21\nLOUVE_19 -->| 86 km| CHALO_21\nVITRY_25 -->| 38 km| CHALO_21\nP_BLAISE -->| 82 km| CHALO_21\nP_MARNE -->| 3 km| STDIZ_04\nR_MARNE -->| 57 km| CHALO_21\nstyle CHALO_21 fill:#cfc\nstyle STDIZ_04 fill:#cfc\nstyle P_BLAISE fill:#fcc\nstyle P_MARNE fill:#fcc\nstyle R_MARNE fill:#fcc"},"evals":[],"jsHooks":[]}</script> ] --- # Intégrer les influences dans le modèle Ajout des noeuds dans l'objet `GRiwrm` ```r griwrmInf <- rbind(griwrmNat, data.frame(id = c("P_BLAISE", "P_MARNE", "R_MARNE"), down = c("CHALO_21", "STDIZ_04", "CHALO_21"), length = c(82.5, 3, 56.9), model = NA, area = NA)) griwrmInf$model[!is.na(griwrmInf$area)] <- "RunModel_GR4J" ``` | |id |down | length|model | area| |:--|:--------|:--------|------:|:-------------|-------:| |4 |STDIZ_04 |CHALO_21 | 85.570|RunModel_GR4J | 2347.53| |19 |LOUVE_19 |CHALO_21 | 86.165|RunModel_GR4J | 461.74| |25 |VITRY_25 |CHALO_21 | 38.047|RunModel_GR4J | 2109.14| |21 |CHALO_21 |NA | NA|RunModel_GR4J | 6291.55| |1 |P_BLAISE |CHALO_21 | 82.500|NA | NA| |2 |P_MARNE |STDIZ_04 | 3.000|NA | NA| |3 |R_MARNE |CHALO_21 | 56.900|NA | NA| --- # Intégrer les influences dans le modèle .pull-left[ Ajout des débits observés aux prises et restitution ```r load("data/marne_reservoir.RData") str(Qreservoir) ``` ``` ## 'data.frame': 5114 obs. of 3 variables: ## $ P_MARNE : num 0 0 0 0 0 0 0 0 0 0 ... ## $ P_BLAISE: num 0 0 0 0 0 0 0 0 0 0 ... ## $ R_MARNE : num 16.9 16.6 16.6 15.6 15.6 15.9 15.9 16.8 21.4 18.6 ... ``` Les débits sont en m3/s... Les débits prélevés aux prises sont négatifs et les débits restitués sont positifs ] .pull-right[ <!-- --> ] --- # `CreateInputsModel` avec des débits influencés ## Conversion et fusion Par défaut, les débits des noeuds sans surface (`area=NA`) sont en m3 par pas de temps: ```r Qnat2 <- cbind(Q, Qreservoir * 86400) IMinf <- CreateInputsModel(griwrmInf, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = Qnat2) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin STDIZ_04... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin LOUVE_19... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin VITRY_25... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin CHALO_21... ``` --- # `RunModel` avec des débits influencés .pull-left[ L'ajout d'un noeud à l'amont de `STDIZ_04` le transforme en un Modèle SD (*GR+Lag*). Il faut ajouter un paramètre `Velocity` aux paramètres du modèle: ```r ParamNat2 <- ParamNat ParamNat2$STDIZ_04 <- c(1, ParamNat$STDIZ_04) ``` ] .pull-right[ La configuration des noeuds modélisés ne change pas: on peut garder le même objet `RunOptions` ```r OMinf <- RunModel(IMinf, RunOptions = ROnat, Param = ParamNat2) ``` ] ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 46 ## time steps with negative flow, set to zero. ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin VITRY_25... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin CHALO_21... ``` --- # `RunModel` avec des débits influencés .pull-left[ ```r plot(OMinf$CHALO_21, Qinf[I_Run, "CHALO_21"]) ``` <!-- --> ] .pull-right[ ```r plot(attr(OMinf, "Qm3s")) ``` <!-- --> ] --- class: inverse, middle, center # Calage d'un réseau de modèles avec des débits influencés --- # Calage en débit influencé Il suffit de refaire le calage en configuration influencée avec des débits observés influencés: ```r str(Qinf) ``` ``` ## 'data.frame': 5114 obs. of 4 variables: ## $ STDIZ_04: num 0.418 0.402 0.316 0.305 0.272 ... ## $ LOUVE_19: num 0.211 0.19 0.182 0.178 0.195 ... ## $ VITRY_25: num 0.197 0.225 0.188 0.176 0.158 ... ## $ CHALO_21: num 0.42 0.451 0.435 0.402 0.367 ... ``` ```r QinfRes <- cbind(Qinf, Qreservoir * 86400) IMinf <- CreateInputsModel(griwrmInf, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = QinfRes) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin STDIZ_04... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin LOUVE_19... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin VITRY_25... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin CHALO_21... ``` --- # Calage en débit influencé ```r ROinf <- CreateRunOptions( InputsModel = IMinf, IndPeriod_WarmUp = 1:365, IndPeriod_Run = I_Run ) ICinf <- CreateInputsCrit( InputsModel = IMinf, FUN_CRIT = airGR::ErrorCrit_NSE, RunOptions = ROinf, Obs = QinfRes[I_Run,] ) COinf <- CreateCalibOptions(IMinf) ``` --- # Calage en débit influencé ```r OCinf <- Calibration(InputsModel = IMinf, RunOptions = ROinf, InputsCrit = ICinf, CalibOptions = COinf, useUpstreamQsim = TRUE) ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 15.000, 247.151, -0.020, 83.096, 2.384 ## Crit. NSE[Q] = 0.7374 ## Steepest-descent local search in progress ## Calibration completed (47 iterations, 696 runs) ## Param = 0.480, 208.513, -0.201, 73.700, 3.535 ## Crit. NSE[Q] = 0.8424 ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 53 ## time steps with negative flow, set to zero. ``` ``` ## Calibration.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 247.151, -2.376, 83.096, 2.384 ## Crit. NSE[Q] = 0.8780 ## Steepest-descent local search in progress ## Calibration completed (23 iterations, 249 runs) ## Param = 180.270, -2.700, 122.413, 2.264 ## Crit. NSE[Q] = 0.8914 ## Calibration.GRiwrmInputsModel: Treating sub-basin VITRY_25... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (81 runs) ## Param = 432.681, -0.649, 83.096, 2.384 ## Crit. NSE[Q] = 0.8225 ## Steepest-descent local search in progress ## Calibration completed (29 iterations, 298 runs) ## Param = 254.678, -1.738, 125.211, 4.238 ## Crit. NSE[Q] = 0.9318 ## Calibration.GRiwrmInputsModel: Treating sub-basin CHALO_21... ## Grid-Screening in progress (0% 20% 40% 60% 80% 100%) ## Screening completed (243 runs) ## Param = 11.250, 432.681, -2.376, 83.096, 2.384 ## Crit. NSE[Q] = 0.8516 ## Steepest-descent local search in progress ## Calibration completed (63 iterations, 861 runs) ## Param = 0.621, 454.535, -4.306, 53.953, 10.371 ## Crit. NSE[Q] = 0.9357 ``` --- # Simulation en débit influencé .pull-left[ ```r ParamInf <- sapply(names(OCinf), function(x) {OCinf[[x]]$Param}) OMinf2 <- RunModel(IMinf, RunOptions = ROinf, Param = ParamInf) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 53 ## time steps with negative flow, set to zero. ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin VITRY_25... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin CHALO_21... ``` ```r plot(OMinf2$CHALO_21, Qinf[I_Run, "CHALO_21"]) ``` ] .pull-right[ <!-- --> ] --- # Question bonus: comment simuler des débits naturalisés ? -- .pull-left[ ```r ParamInf2Nat <- ParamInf # St-Dizier est un noeud amont en Q nat ParamInf2Nat$STDIZ_04 <- ParamInf2Nat$STDIZ_04[2:5] OMinf2nat <- RunModel( IMnat, RunOptions = ROnat, Param = ParamInf2Nat) ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin STDIZ_04... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin LOUVE_19... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin VITRY_25... ``` ``` ## RunModel.GRiwrmInputsModel: Treating sub-basin CHALO_21... ``` ```r plot(OMinf2nat$CHALO_21, Qinf[I_Run, "CHALO_21"]) ``` ] .pull-right[ <!-- --> ] --- class: inverse, middle, center # Simulation d'un réservoir régulé --- # Présentation des règles de gestion du réservoir Dans cet exemple, le réservoir ne sert qu'à protéger l'aval contre les crues en limitant le débit à `CHALO_21` à 330&nbsp;m³/s. .pull-left[ Caractéristiques physiques du réservoir: - Capacité utiles du réservoir 350 millions de m³ - Limite débit canal amenée de la Marne: 375&nbsp;m³/s - Limite débit canal de restitution Marne: 50&nbsp;m³/s - Pour simplifier le problème, on condamne le canal d'amenée de la Blaise. Les obligations réglementaires sont&nbsp;: - Le débit réservé à la prise d'eau de la Marne est de 8&nbsp;m³/s La règle de gestion consiste à garder le réservoir aussi vide que possible tout en respectant le débit maximal à l'aval. ] .pull-right[ Le débit de la "prise Blaise" est fixé à zéro Les débits de la prise Marne et de la restitution seront calculés pendant la simulation ```r Qreg <- QinfRes Qreg$P_BLAISE <- 0 IMreg <- CreateInputsModel( griwrmInf, DatesR = DatesR, Precip = P, PotEvap = E, Qobs = Qreg ) ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin STDIZ_04... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin LOUVE_19... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin VITRY_25... ``` ``` ## CreateInputsModel.GRiwrm: Treating sub-basin CHALO_21... ``` ] --- # Logique de contrôle à la prise .pull-left[ Elle s'écrit sous la forme d'une fonction qui prend les données mesurées `Y` sur le réseau au pas de temps précédent(s) et renvoie les commandes `U` qui écraseront la valeur du débit au noeud au pas de temps en courant (et suivants). Ici, les noeuds concernées sont: -- - `CHALO_21` (`Y[1]`) pour limiter le débit aval - `STDIZ_04` (`Y[2]`) pour contrôler le débit réservé à la prise - `P_MARNE` (`Y[3]`) pour renaturaliser le débit aval pour la commande -- - `P_MARNE` (`U[1]`) pour commander la prise d'eau ] .pull-right[ ```r FactoryLogicPrise <- function(sv, Qmax, Qlim, Qres, Vmax) { function(Y) { # Flood control if(Y[1] - Y[3] > Qmax) { U <- Y[1] - Y[3] - Qmax if(Y[2] - Y[3] - U < Qres) { U <- Y[2] - Y[3] - Qres } # Physical flow limits U <- max(0, min(Qlim, U)) # Reservoir storage sv$Vres <- sv$Vres + U if(sv$Vres > Vmax) { U <- U - (sv$Vres - Vmax) sv$Vres <- Vmax } } else { U <- 0 } return(-U) # Flow withdrown from the network } } ``` ] --- # Logique de contrôle à la restitution .pull-left[ A la restitution, l'objectif est de vider le réservoir dans le respect du débit de seuil de crue à `CHALO_21`, on a : - `Y[1]`: le débit à `CHALO_21` - `U[1]`: le débit du canal de restitution Ce contrôleur enregistre la chronique de volume du réservoir pour pouvoir le tracé après la simulation. ] .pull-right[ ```r FactoryLogicRestit <- function(sv, Qmax, Qlim) { function(Y) { # Flood control if(Y[1] < Qmax) { # Vidange dans limite du seuil de crue # et du débit limite du canal U <- min(Qlim, Qmax - Y[1]) U <- max(0, U) # Reservoir storage sv$Vres <- sv$Vres - U if(sv$Vres < 0) { U <- U + sv$Vres sv$Vres <- 0 } } else { U <- 0 } # Record Volume storage time series sv$VresChro <- c(sv$VresChro, sv$Vres) return(U) } } ``` ] --- # Le superviseur .pull-left[ Le supervisor coordonne les interactions entre la simulation et les contrôleurs C'est un environnement qui contient toutes les données utiles à la simulation / régulation On peut y ajouter ces propres variables qui seront utilisées pendant la régulation ] .pull-right[ ```r sv <- CreateSupervisor(IMreg, TimeStep = 1L) # Initial reservoir storage sv$Vres <- 0 # Reservoir storage time series sv$VresChro <- sv$Vres ``` ```r names(sv) ``` ] ``` ## [1] "ts.index0" "CreateController" "ts.date" "supervisor" ## [5] "ts.previous" "DatesR" "OutputsModel" "griwrm" ## [9] "InputsModel" ".TimeStep" "controllers" ".isSupervisor" ## [13] "VresChro" "Vres" "controller.id" "ts.index" ``` --- # Les contrôleurs Les contrôleurs embarquent une logique de contrôle et définissent les noeuds "mesurés" `Y` et "commandés" `U` Les fonctions de logique de contrôle encapsulent les données nécessaires à leur fonctionnement ```r fprise <- FactoryLogicPrise(sv = sv, Qmax = 150 * 86400, Qlim = 375 * 86400, Qres = 8 * 86400, Vmax = 350 * 1e6) CreateController(sv, ctrl.id = "Prise", Y = c("CHALO_21", "STDIZ_04", "P_MARNE"), U = "P_MARNE", FUN = fprise) ``` ``` ## The controller has been added to the supervisor ``` ```r frestit <- FactoryLogicRestit(sv = sv, Qmax = 200 * 86400, Qlim = 50 * 86400) CreateController(sv, ctrl.id = "Restitution", Y = "CHALO_21", U = "R_MARNE", FUN = frestit) ``` ``` ## The controller has been added to the supervisor ``` --- # `RunModel` avec un `Supervisor` ```r # Initial reservoir storage sv$Vres <- 0 # Reservoir storage time series sv$VresChro <- sv$Vres time.start <- Sys.time() OMreg <- RunModel(sv, RunOptions = ROinf, Param = ParamInf) Sys.time() - time.start ``` ``` ## Time difference of 6.398789 secs ``` --- # Resultat de la simulation avec régulation .pull-left[ ```r I_Zoom <- 1180:1380 Qm3sReg <- attr(OMreg, "Qm3s") plot(Qm3sReg[I_Zoom,]) ``` <!-- --> ] .pull-right[ ```r plot(Qm3sReg$DatesR[I_Zoom], sv$VresChro[I_Zoom], type = "l") ``` <!-- --> ] --- 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** - Ajout d'une méthode de régalurisation pour les modèles semi-distribués ## 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 😃