feat: Calibration with ungauged nodes

Refs #42

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)]])
+}

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)
 }

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