fix: calibration mixing ungauged and reservoir nodes

Refs #110

R/Calibration.GRiwrmInputsModel.R

@@ -48,6 +48,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       l  <- updateParameters4Ungauged(id,
                                       InputsModel,
                                       RunOptions,
+                                      CalibOptions,
                                       OutputsModel,
                                       useUpstreamQsim)
       IM <- l$InputsModel
@@ -84,14 +85,24 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
         subBasinAreas <- calcSubBasinAreas(IM)
       }
       for (uId in Ids) {
-        # Add OutputsCalib for ungauged nodes
-        OutputsCalib[[uId]] <- OutputsCalib[[id]]
-        # Copy parameters and transform X4 relatively to the sub-basin area
-        OutputsCalib[[uId]]$ParamFinalR <-
-          OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$indexParamUngauged]
-        if(IM[[id]]$model$hasX4) {
-          OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$iX4] <-
-            X4 * (subBasinAreas[uId] / sum(subBasinAreas)) ^ 0.3
+        if(!IM[[uId]]$isReservoir) {
+          # Add OutputsCalib for ungauged nodes
+          OutputsCalib[[uId]] <- OutputsCalib[[id]]
+          # Copy parameters and transform X4 relatively to the sub-basin area
+          OutputsCalib[[uId]]$ParamFinalR <-
+            OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$indexParamUngauged]
+          if(IM[[id]]$model$hasX4) {
+            OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$iX4] <-
+              X4 * (subBasinAreas[uId] / sum(subBasinAreas, na.rm = TRUE)) ^ 0.3
+          }
+        } else {
+          OutputsCalib[[uId]] <- Calibration(
+            InputsModel = IM[[uId]],
+            RunOptions = RunOptions[[uId]],
+            InputsCrit = IC,
+            CalibOptions = CalibOptions[[uId]],
+            ...
+          )
         }
       }
       IM <- IM[[id]]
@@ -171,6 +182,7 @@ reduceGRiwrmObj4Ungauged <- function(griwrm, obj) {
 updateParameters4Ungauged <- function(GaugedId,
                                       InputsModel,
                                       RunOptions,
+                                      CalibOptions,
                                       OutputsModel,
                                       useUpstreamQsim) {
 
@@ -188,9 +200,15 @@ updateParameters4Ungauged <- function(GaugedId,
   ### Modify InputsModel for the reduced network ###
   # Remove nodes outside of reduced network
   InputsModel <- reduceGRiwrmObj4Ungauged(g, InputsModel)
+  # Copy fixed parameters for Reservoirs
+  for (id in names(InputsModel)) {
+    if (InputsModel[[id]]$isReservoir) {
+      InputsModel[[id]]$FixedParam <- CalibOptions[[id]]$FixedParam
+    }
+  }
   # Update griwrm
   attr(InputsModel, "GRiwrm") <- g
-  # Update Qupstream already modelled in the reduced network upstream nodes
+  # Update Qupstream already modeled in the reduced network upstream nodes
   idIM <- unique(g$down[g$id %in% upIds])
   for (id in idIM) {
     if(useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsModeled)) {
@@ -258,9 +276,12 @@ calcSubBasinAreas <- function(IM) {
 #' @noRd
 RunModel_Ungauged <- function(InputsModel, RunOptions, Param) {
   InputsModel$FUN_MOD <- NULL
-  SBVI <- sum(calcSubBasinAreas(InputsModel))
+  SBVI <- sum(calcSubBasinAreas(InputsModel), na.rm = TRUE)
   # Compute Param for each sub-basin
   P <- lapply(InputsModel, function(IM) {
+    if (IM$isReservoir) {
+      return(IM$FixedParam)
+    }
     p <- Param[IM$model$indexParamUngauged]
     if(IM$model$hasX4) {
       p[IM$model$iX4] <- Param[IM$model$iX4] * (IM$BasinAreas[length(IM$BasinAreas)] / SBVI) ^ 0.3

tests/testthat/test-RunModel_Reservoir.R

@@ -79,4 +79,22 @@ test_that("Calibration with ungauged node and reservoir in the middle works", {
 
   expect_equal(g$donor[g$id == "54095"], "54001")
 
+  e <- setupRunModel(griwrm = g, runRunModel = FALSE, Qobs2 = Qobs2)
+  for(x in ls(e)) assign(x, get(x, e))
+
+  InputsCrit <- CreateInputsCrit(InputsModel,
+                                 ErrorCrit_KGE2,
+                                 RunOptions = RunOptions,
+                                 Obs = Qobs[IndPeriod_Run, ])
+  CalibOptions <- CreateCalibOptions(InputsModel)
+  CalibOptions[["Dam"]]$FixedParam <- c(650E6, 1)
+  OC <- Calibration(
+    InputsModel = InputsModel,
+    RunOptions = RunOptions,
+    InputsCrit = InputsCrit,
+    CalibOptions = CalibOptions
+  )
+  # X1, X2, X3 are identical
+  expect_equal(OC$`54001`$ParamFinalR[2:4], OC$`54095`$ParamFinalR[1:3])
+  expect_equal(OC$Dam$ParamFinalR, CalibOptions[["Dam"]]$FixedParam)
 })