From e4d5e384300a3cbc0ad8475b5c737b50cad5744e Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@irstea.fr>
Date: Mon, 25 Jul 2022 09:28:21 +0200
Subject: [PATCH] feat: Calibration with ungauged nodes
Refs #42
---
R/Calibration.GRiwrmInputsModel.R | 191 +++++++++++++++++++--
R/CreateRunOptions.GRiwrmInputsModel.R | 1 +
vignettes/V05_Modelling_ungauged_nodes.Rmd | 93 ++++++++--
3 files changed, 257 insertions(+), 28 deletions(-)
diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
index e53eed8..0ea852e 100644
--- a/R/Calibration.GRiwrmInputsModel.R
+++ b/R/Calibration.GRiwrmInputsModel.R
@@ -31,34 +31,65 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
OutputsModel <- list()
class(OutputsModel) <- append("GRiwrmOutputsModel", class(OutputsModel))
- for(IM in InputsModel) {
- message("Calibration.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
+ b <- sapply(InputsModel, function(IM) !IM$isUngauged)
+ gaugedIds <- names(b[b])
- if(useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
- # Update InputsModel$Qupstream with simulated upstream flows
- IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]], OutputsModel)
- }
+ for(id in gaugedIds) {
+ IM <- InputsModel[[id]]
+ message("Calibration.GRiwrmInputsModel: Treating sub-basin ", id, "...")
- if (inherits(InputsCrit[[IM$id]], "InputsCritLavenneFunction")) {
- IC <- getInputsCrit_Lavenne(IM$id, OutputsModel, InputsCrit)
+ if (inherits(InputsCrit[[id]], "InputsCritLavenneFunction")) {
+ IC <- getInputsCrit_Lavenne(id, OutputsModel, InputsCrit)
+ } else {
+ IC <- InputsCrit[[id]]
+ }
+ hasUngauged <- IM$hasUngauged
+ if (hasUngauged) {
+ Sdown <- IM$BasinAreas[length(IM$BasinAreas)]
+ l <- updateParameters4Ungauged(id,
+ InputsModel,
+ RunOptions,
+ OutputsModel,
+ useUpstreamQsim)
+ IM <- l$InputsModel
+ IM$FUN_MOD <- "RunModel_Ungauged"
+ attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions
} else {
- IC <- InputsCrit[[IM$id]]
+ if(useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
+ # Update InputsModel$Qupstream with simulated upstream flows
+ IM <- UpdateQsimUpstream(IM, RunOptions[[id]], OutputsModel)
+ }
}
- OutputsCalib[[IM$id]] <- Calibration(
+ OutputsCalib[[id]] <- Calibration(
InputsModel = IM,
- RunOptions = RunOptions[[IM$id]],
+ RunOptions = RunOptions[[id]],
InputsCrit = IC,
- CalibOptions = CalibOptions[[IM$id]],
+ CalibOptions = CalibOptions[[id]],
...
)
+ if (hasUngauged) {
+ # Add OutputsCalib for ungauged nodes
+ g <- attr(InputsModel, "GRiwrm")
+ ungaugedIds <- g$id[g$gauged == id & g$id != id & !is.na(g$model)]
+ for (uId in ungaugedIds) {
+ OutputsCalib[[uId]] <- OutputsCalib[[id]]
+ PS <- attr(IM[[uId]], "ParamSettings")
+ OutputsCalib[[uId]]$ParamFinalR <- OutputsCalib[[uId]]$ParamFinalR[PS$Indexes]
+ if(PS$hasX4) {
+ OutputsCalib[[uId]]$ParamFinalR[PS$iX4] <-
+ OutputsCalib[[uId]]$ParamFinalR[PS$iX4] * (sum(IM[[uId]]$BasinAreas) / Sdown) ^ 0.3
+ }
+ }
+ }
+
if(useUpstreamQsim) {
# Run the model for the sub-basin
- OutputsModel[[IM$id]] <- RunModel(
+ OutputsModel[[id]] <- RunModel(
x = IM,
- RunOptions = RunOptions[[IM$id]],
- Param = OutputsCalib[[IM$id]]$ParamFinalR
+ RunOptions = RunOptions[[id]],
+ Param = OutputsCalib[[id]]$ParamFinalR
)
}
@@ -96,3 +127,133 @@ getInputsCrit_Lavenne <- function(id, OutputsModel, InputsCrit) {
AprCrit <- ErrorCrit(InputsCrit[[AprioriId]], OutputsModel[[AprioriId]])$CritValue
return(Lavenne_FUN(AprParamR, AprCrit))
}
+
+
+#' Reduce a GRiwrm list object (InputsModel, RunOptions...) for a reduced network
+#'
+#' @param griwrm See [CreateGRiwrm])
+#' @param obj Either a *GRiwrmInputsModel*, *GRiwrmOptions*... object
+#'
+#' @return The object containing only nodes of the reduced model
+#' @noRd
+reduceGRiwrmObj4Ungauged <- function(griwrm, obj) {
+ objAttributes <- attributes(obj)
+ obj <- lapply(obj, function(o) {
+ if(o$id %in% griwrm$id) {
+ o
+ } else {
+ NULL
+ }
+ })
+ obj[sapply(obj, is.null)] <- NULL
+ objAttributes$names <- names(obj)
+ attributes(obj) <- objAttributes
+ return(obj)
+}
+
+updateParameters4Ungauged <- function(GaugedId,
+ InputsModel,
+ RunOptions,
+ OutputsModel,
+ useUpstreamQsim) {
+
+ ### Set the reduced network of the basin containing ungauged nodes ###
+ # Select nodes identified with the current node as gauged node
+ griwrm <- attr(InputsModel, "GRiwrm")
+ g <- griwrm[griwrm$gauged == GaugedId, ]
+ # Add upstream nodes for routing upstream flows
+ upIds <- griwrm$id[griwrm$down %in% g$id & !griwrm$id %in% g$id]
+ g <- rbind(griwrm[griwrm$id %in% upIds, ], g)
+ g$model[g$id %in% upIds] <- NA
+ # Set downstream node
+ g$down[!g$down %in% g$id] <- NA
+
+ ### Modify InputsModel for the reduced network ###
+ # Remove nodes outside of reduced network
+ InputsModel <- reduceGRiwrmObj4Ungauged(g, InputsModel)
+ # Update griwrm
+ attr(InputsModel, "GRiwrm") <- g
+ # Update Qupstream of reduced network upstream nodes
+ g2 <- griwrm[griwrm$gauged == GaugedId,]
+ upIds2 <- g2$id[!g2$id %in% g2$down]
+ for (id in upIds2) {
+ if(useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsRunoff)) {
+ # Update InputsModel$Qupstream with simulated upstream flows
+ InputsModel[[id]] <- UpdateQsimUpstream(InputsModel[[id]],
+ RunOptions[[id]],
+ OutputsModel)
+ InputsModel[[id]]$UpstreamIsRunoff <-
+ rep(FALSE, length(InputsModel[[id]]$UpstreamIsRunoff))
+ }
+ }
+ # Add extra info for Param processing
+ nbParam <- RunOptions[[GaugedId]]$FeatFUN_MOD$NbParam
+ for (id in names(InputsModel)) {
+ attr(InputsModel[[id]], "ParamSettings") <-
+ list(Indexes = ifelse(inherits(InputsModel[[id]], "SD"), 1, 2):nbParam,
+ hasX4 = grepl("RunModel_GR[456][HJ]", InputsModel[[id]]$FUN_MOD),
+ iX4 = ifelse(inherits(InputsModel[[id]], "SD"), 5, 4))
+ }
+ # Add class InputsModel for airGR::Calibration checks
+ class(InputsModel) <- c("InputsModel", class(InputsModel))
+
+ ### Modify RunOptions for the reduced network ###
+ RunOptions <- reduceGRiwrmObj4Ungauged(g, RunOptions)
+ return(list(InputsModel = InputsModel, RunOptions = RunOptions))
+}
+
+updateRunOptions4Ungauged <- function(id, RunOptions) {
+ # Remove nodes outside of reduced network
+ RunOptions <- lapply(RunOptions, function(RO) {
+ if(RO$id %in% g$id) {
+ IM
+ } else {
+ NULL
+ }
+ })
+ return(RunOptions)
+}
+
+
+#' RunModel for a sub-network of ungauged nodes
+#'
+#' The function simulates a network with one set of parameters
+#' shared with ungauded nodes inside the basin.
+#'
+#' @details
+#' The network should contains only one gauged station at downstream and other
+#' nodes can be direct injection or ungauged nodes.
+#'
+#' This function works as functions similar to [airGR::RunModel_GR4J] except that
+#' `InputsModel` is a *GRiwrmInputsModel* containing the network of ungauged nodes
+#' and direct injection in the basin.
+#'
+#' `Param` is adjusted for each sub-basin using the method developped by
+#' Lobligeois (2014) for GR models.
+#'
+#' @references Lobligeois, Florent. Mieux connaître la distribution spatiale des
+#' pluies améliore-t-il la modélisation des crues ? Diagnostic sur 181 bassins
+#' versants français. Phdthesis, AgroParisTech, 2014.
+#' <https://pastel.archives-ouvertes.fr/tel-01134990/document>
+#'
+#' @inheritParams airGR::RunModel
+#'
+#' @inherit RunModel.GRiwrmInputsModel return return
+#' @noRd
+RunModel_Ungauged <- function(InputsModel, RunOptions, Param) {
+ InputsModel$FUN_MOD <- NULL
+ SBV <- sum(InputsModel[[length(InputsModel)]]$BasinAreas)
+ # Compute Param for each sub-basin
+ P <- lapply(InputsModel, function(IM) {
+ PS <- attr(IM, "ParamSettings")
+ p <- Param[PS$Indexes]
+ if(PS$hasX4) {
+ p[PS$iX4] <- Param[PS$iX4] * (IM$BasinAreas[length(IM$BasinAreas)] / SBV) ^ 0.3
+ }
+ return(p)
+ })
+ OM <- suppressMessages(
+ RunModel.GRiwrmInputsModel(InputsModel, attr(RunOptions, "GRiwrmRunOptions"), P)
+ )
+ return(OM[[length(OM)]])
+}
diff --git a/R/CreateRunOptions.GRiwrmInputsModel.R b/R/CreateRunOptions.GRiwrmInputsModel.R
index e7f503c..7b49fae 100644
--- a/R/CreateRunOptions.GRiwrmInputsModel.R
+++ b/R/CreateRunOptions.GRiwrmInputsModel.R
@@ -8,6 +8,7 @@ CreateRunOptions.GRiwrmInputsModel <- function(x, IniStates = NULL, ...) {
for(id in names(x)) {
RunOptions[[id]] <- CreateRunOptions(x[[id]], IniStates = IniStates[[id]], ...)
+ RunOptions[[id]]$id <- id
}
return(RunOptions)
}
diff --git a/vignettes/V05_Modelling_ungauged_nodes.Rmd b/vignettes/V05_Modelling_ungauged_nodes.Rmd
index 586b835..4aace25 100644
--- a/vignettes/V05_Modelling_ungauged_nodes.Rmd
+++ b/vignettes/V05_Modelling_ungauged_nodes.Rmd
@@ -25,20 +25,21 @@ library(airGRiwrm)
## Why modelling ungauged station in the semi-distributed model?
-This vignette introduces the implementation in airGRiwrm of the model developped by @lobligeoisMieuxConnaitreDistribution2014
-
-2 interests:
+Ungauged nodes in the semi-distributed model can be used to reach two different goals:
- increase spatial resolution of the rain fall to improve streamflow simulation [@lobligeoisWhenDoesHigher2014].
- simulate streamflows in location of interest for management purpose
+This vignette introduces the implementation in airGRiwrm of the method developped by @lobligeoisMieuxConnaitreDistribution2014 for calibrating ungauged nodes in a
+semi-distributed model.
+
## Presentation of the study case
-Using the study case of the vignette #1 and #2, we considere this time that nodes `54095`, `54001` and
+Using the study case of the vignette #1 and #2, we considere this time that nodes `54001` and
`54029` are ungauged. We simulate the streamflow at these locations by sharing
hydrological parameters of the gauged node `54032`.
-```{r, echo = FALSE}
+```{r network, echo = FALSE}
mmd <- function(x, ...) {
# For avoiding crash of R CMD build in console mode
if(Sys.getenv("RSTUDIO") == "1") {
@@ -52,9 +53,10 @@ id95[54095]
id01[54001]
id29[54029]
+id95 -->| 42 km| id01
+
subgraph Shared parameters from node 54032
id01 -->| 45 km| 54032
-id95 -->| 42 km| id01
id29 -->| 32 km| 54032
end
@@ -67,10 +69,10 @@ classDef IntUng fill:#efe
classDef IntGau fill:#afa
classDef DirInj fill:#faa
-class id95,id29 UpUng
+class id29 UpUng
class 54057,54032 IntGau
class id01 IntUng
-class 54002 UpGau
+class id95,54002 UpGau
")
```
@@ -85,11 +87,11 @@ With $X_4$ the unit hydrogram parameter for the entire basin at `54032` which as
Ungauged stations are specified by using the model "Ungauged" in the `model` column provided in the `CreateGRiwrm` function:
-```{r}
+```{r griwrm}
data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$model <- "RunModel_GR4J"
-nodes$model[nodes$gauge_id %in% c("54095", "54029", "54001")] <- "Ungauged"
+nodes$model[nodes$gauge_id %in% c("54029", "54001")] <- "Ungauged"
griwrmV05 <- CreateGRiwrm(
nodes,
list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
@@ -99,7 +101,7 @@ griwrmV05
On the following network scheme, the ungauged nodes are cleared than gauged ones with the same color (blue for upstream nodes and green for intermediate and downstream nodes)
-```{r}
+```{r plot_network}
plot(griwrmV05)
```
@@ -108,7 +110,7 @@ plot(griwrmV05)
The formatting of the input data is described in the vignette "V01_Structure_SD_model". The following code chunk resumes this formatting procedure:
-```{r}
+```{r obs}
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
@@ -120,8 +122,73 @@ Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
Then, the `GRiwrmInputsModel` object can be generated taking into account the new `GRiwrm` object:
-```{r}
+```{r InputsModel}
IM_U <- CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap, Qobs)
```
+## Calibration of the model integrating ungauged nodes
+
+Calibration options is detailed in vignette "V02_Calibration_SD_model".
+We also apply a parameter regularization here but only where an upstream simulated catchment is available.
+
+The following code chunk resumes this procedure:
+
+```{r RunOptions}
+IndPeriod_Run <- seq(
+ which(DatesR == (DatesR[1] + 365*24*60*60)), # Set aside warm-up period
+ length(DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)
+RunOptions <- CreateRunOptions(IM_U,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run)
+InputsCrit <- CreateInputsCrit(IM_U,
+ FUN_CRIT = ErrorCrit_KGE2,
+ RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,],
+ AprioriIds = c("54057" = "54032", "54032" = "54095"),
+ transfo = "sqrt", k = 0.15
+)
+CalibOptions <- CreateCalibOptions(IM_U)
+```
+
+The **airGR** calibration process is applied on each hydrological node of the `GRiwrm` network from upstream nodes to downstream nodes but this time the calibration of the sub-basin `54032` invokes a semi-distributed model composed of the nodes `54029`, `54001` and `54032`.
+
+```{r Calibration}
+OC_U <- suppressWarnings(
+ Calibration(IM_U, RunOptions, InputsCrit, CalibOptions))
+```
+
+Hydrological parameters for sub-basins
+
+## Run the model with the optimized model parameters
+
+The hydrologic model uses parameters herited from the calibration of the gauged sub-basin `54032` for the ungauged nodes `54001` and `54029`:
+
+```{r param}
+ParamV05 <- sapply(griwrmV05$id, function(x) {OC_U[[x]]$Param})
+dfParam <- do.call(
+ rbind,
+ lapply(ParamV05, function(x)
+ if(length(x)==4) {return(c(NA, x))} else return(x))
+)
+colnames(dfParam) <- c("velocity", paste0("X", 1:4))
+knitr::kable(round(dfParam, 3))
+```
+
+We can run the model with these calibrated parameters:
+
+```{r RunModel}
+OutputsModels <- RunModel(
+ IM_U,
+ RunOptions = RunOptions,
+ Param = ParamV05
+)
+```
+
+and plot the comparison of the modelled and the observed flows including on the so-called "ungauged" stations :
+
+```{r plot, fig.height = 5, fig.width = 8}
+plot(OutputsModels, Qobs = Qobs[IndPeriod_Run,], which = c("Regime", "CumFreq"))
+```
+
# References
--
GitLab