Merge branch '167-calibration-crash-with-sibling-receiver-node' into …

…'dev' Resolve "Calibration: crash with downstream receiver node with upstream reservoir" Closes #167 See merge request in-wop/airGRiwrm!102
inrae · Aug 22, 2024 · 9167c30 · 9167c30
2 parents f0df388 + 7b331fa
commit 9167c30
Show file tree

Hide file tree

Showing 12 changed files with 104 additions and 39 deletions.
diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
@@ -36,14 +36,13 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
 
   for (id in gaugedIds) {
     IM <- InputsModel[[id]]
-    message("Calibration.GRiwrmInputsModel: Processing sub-basin ", id, "...")
 
     if (inherits(InputsCrit[[id]], "InputsCritLavenneFunction")) {
       IC <- getInputsCrit_Lavenne(id, OutputsModel, InputsCrit)
     } else {
       IC <- InputsCrit[[id]]
     }
-    hasUngauged <- IM$hasUngauged
+    hasUngauged <- IM$hasUngaugedNodes
     if (hasUngauged) {
       l  <- updateParameters4Ungauged(id,
                                       InputsModel,
@@ -52,9 +51,12 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
                                       OutputsModel,
                                       useUpstreamQsim)
       IM <- l$InputsModel
+      message("Calibration.GRiwrmInputsModel: Processing sub-basins '",
+              paste(names(IM), collapse = "', '"), "' with '", id, "' as gauged donor...")
       IM$FUN_MOD <- "RunModel_Ungauged"
       attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions
     } else {
+      message("Calibration.GRiwrmInputsModel: Processing sub-basin '", id, "'...")
       if (useUpstreamQsim && any(IM$UpstreamIsModeled)) {
         # Update InputsModel$Qupstream with simulated upstream flows
         IM <- UpdateQsimUpstream(IM, RunOptions[[id]], OutputsModel)
@@ -74,7 +76,8 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
                                        OutputsCalib[[IM$gaugedId]]$ParamFinalR,
                                        IM$gaugedId,
                                        id,
-                                       CalibOptions[[id]]$FixedParam)
+                                       CalibOptions[[id]]$FixedParam,
+                                       verbose = TRUE)
       )
       class(OutputsCalib[[id]]) <- c("OutputsCalib", class(OutputsCalib[[id]]))
     } else {
@@ -99,7 +102,8 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
             ParamFinalR = transferGRparams(InputsModel,
                                            OutputsCalib[[id]]$ParamFinalR,
                                            id,
-                                           uId)
+                                           uId,
+                                           verbose = TRUE)
           )
           class(OutputsCalib[[uId]]) <- class(OutputsCalib[[id]])
         } else {

diff --git a/R/CreateCalibOptions.GRiwrmInputsModel.R b/R/CreateCalibOptions.GRiwrmInputsModel.R
@@ -8,7 +8,7 @@ CreateCalibOptions.GRiwrmInputsModel <- function(x, FixedParam = NULL, ...) {
   }
   np <- getAllNodesProperties(attr(x, "GRiwrm"))
   np <- np[!np$DirectInjection, ]
-  gaugedIds <- np$id[np$calibration == "Gauged"]
+  gaugedIds <- np$id[np$gauged]
 
   if (!is.null(FixedParam)) {
     if (!(is.list(FixedParam) || is.numeric(FixedParam))) {
@@ -50,7 +50,7 @@ CreateCalibOptions.GRiwrmInputsModel <- function(x, FixedParam = NULL, ...) {
     if (!is.null(FixedParam)) {
       FP <- FixedParam[[id]]
     }
-    if (np$calibration[np$id == id] == "Gauged") {
+    if (np$gauged[np$id == id]) {
       CalibOptions[[IM$id]] <- CreateCalibOptions(
         IM,
         FixedParam = FP,

diff --git a/R/CreateGRiwrm.R b/R/CreateGRiwrm.R
@@ -105,6 +105,7 @@ CreateGRiwrm <- function(db,
 
   # Set automatic downstream donors for ungauged nodes
   griwrm$donor <- setDonor(griwrm)
+  checkUngaugedCluster(griwrm)
 
   griwrm <- sort(griwrm)
 
@@ -291,7 +292,7 @@ setDonor <- function(griwrm) {
     }
     id <- griwrm$id[i]
     model <- griwrm$model[i]
-    if (is.na(model)) {
+    if (is.na(model) || model == "Diversion") {
       return(as.character(NA))
     }
     if (model == "RunModel_Reservoir" && is.na(griwrm$down[i])){
@@ -340,3 +341,22 @@ refineReservoirDonor <- function(i, griwrm) {
   # No upstream ungauged nodes: Reservoir is its own donor!
   return(griwrm$id[i])
 }
+
+checkUngaugedCluster <- function(griwrm) {
+  # Check presence of gauged nodes inside an ungauged cluster
+  clusters <- table(griwrm$donor)
+  clusters <- names(clusters[clusters > 1])
+  lapply(clusters, function(gaugedId) {
+    g <- getUngaugedCluster(griwrm, gaugedId)
+    p <- getAllNodesProperties(griwrm)
+    upstreamIdsInCluster <- unique(g$id[!g$id %in% g$down])
+    lapply(g$id, function(id) {
+      if (id != gaugedId) {
+        if (p$calibration[p$id == id] == "Gauged" && !id %in% upstreamIdsInCluster) {
+          stop("The gauged node '", id, "' is located in the cluster of the ungauged",
+               " nodes calibrated with the node '", gaugedId, "'")
+        }
+      }
+    })
+  })
+}
diff --git a/R/getNodeProperties.R b/R/getNodeProperties.R
@@ -16,7 +16,7 @@
 #' - "Upstream" ([logical]): is the node an upstream node?
 #' - "RunOff" ([logical]): is the node contains an hydrological model?
 #'
-#' `getAllNodeProperties` returns a [data.frame] constituted from the list returned
+#' `getAllNodesProperties` returns a [data.frame] constituted from the list returned
 #' by `getNodeProperties` for all nodes.
 #'
 #' @details
@@ -50,10 +50,23 @@ getNodeProperties <- function(id, griwrm) {
     Reservoir = !is.na(model) && model == "RunModel_Reservoir",
     airGR = grepl("RunModel_", donor_model)
   )
+  p$gauged <- isNodeGauged(id, griwrm)
   if (p$DirectInjection) {
     p$calibration <- "NA"
   } else {
-    p$calibration <- ifelse(isNodeGauged(id, griwrm), "Gauged", "Ungauged")
+    if (p$gauged) {
+      if (p$Reservoir) {
+        p$calibration <- "Reservoir"
+      } else {
+        p$calibration <- "Gauged"
+      }
+    } else {
+      if (is.na(griwrm$donor[id]) || isNodeDownstream(griwrm, id, griwrm$donor[id])) {
+        p$calibration <- "Ungauged"
+      } else {
+        p$calibration <- "Receiver"
+      }
+    }
   }
   p$Upstream <- p$position == "Upstream"
   p$RunOff <- !p$DirectInjection && !p$Reservoir && donor_model != "RunModel_Lag"

diff --git a/R/getNodeRanking.R b/R/getNodeRanking.R
@@ -36,7 +36,8 @@ getNodeRanking <- function(griwrm) {
       upId <- upIds[1]
       #Browse the ungauged sub-network until the donor
       upDonor <- unique(g$donor[g$id == upId])
-      g2 <- g[g$donor == upDonor, ]
+      cluster_nodes <- g$id[!is.na(g$donor) & g$donor == upDonor]
+      g2 <- g[g$id %in% cluster_nodes, ]
       # Check if upstream nodes have already been processed
       immediate_upstream_nodes <- g$id[!is.na(g$down) & g$down %in% g2$id]
       immediate_upstream_nodes <- immediate_upstream_nodes[!immediate_upstream_nodes %in% g2$id]

diff --git a/R/plot.GRiwrm.R b/R/plot.GRiwrm.R
@@ -144,7 +144,8 @@ getNodeClass <- function(id, griwrm) {
   }  else if (props$Reservoir) {
     nc <- "Reservoir"
   }  else {
-    nc <- paste0(props$position, props$calibration)
+    nc <- paste0(props$position,
+                 ifelse(props$gauged, "Gauged", "Ungauged"))
   }
   if (props$Diversion) nc <- paste0(nc, "Diversion")
   return(nc)

diff --git a/R/utils.Calibration.R b/R/utils.Calibration.R
@@ -56,7 +56,6 @@ reduceGRiwrmObj4Ungauged <- function(griwrm, obj) {
   return(obj)
 }
 
-
 #' Set a reduced GRiwrm network for calibration of a sub-network with ungauged
 #' hydrological nodes
 #'
@@ -78,20 +77,7 @@ updateParameters4Ungauged <- function(GaugedId,
                                       OutputsModel,
                                       useUpstreamQsim) {
 
-  ### Set the reduced network of the basin containing ungauged nodes ###
-  # Select nodes identified with the current node as donor gauged node
-  griwrm <- attr(InputsModel, "GRiwrm")
-  donorIds <- griwrm$id[!is.na(griwrm$donor) & griwrm$donor == GaugedId]
-  gDonor <- griwrm %>% dplyr::filter(.data$id %in% donorIds)
-  # Add upstream nodes for routing upstream flows
-  upNodes <- griwrm %>%
-    dplyr::filter(.data$down %in% gDonor$id,
-                  !.data$id %in% gDonor$id) %>%
-    dplyr::mutate(model = ifelse(!is.na(.data$model), NA, .data$model))
-  upIds <- upNodes$id
-  g <- rbind(upNodes, gDonor)
-  # Set downstream nodes
-  g$down[!g$down %in% g$id] <- NA
+  g <- getUngaugedCluster(attr(InputsModel, "GRiwrm"), GaugedId)
 
   ### Modify InputsModel for the reduced network ###
   # Remove nodes outside of reduced network
@@ -105,6 +91,7 @@ updateParameters4Ungauged <- function(GaugedId,
   # Update griwrm
   attr(InputsModel, "GRiwrm") <- g
   # Update Qupstream already modeled in the reduced network upstream nodes
+  upIds <- attr(g, "upIds")
   idIM <- unique(g$down[g$id %in% upIds])
   for (id in idIM) {
     if (useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsModeled)) {
@@ -211,13 +198,17 @@ RunModel_Ungauged <- function(InputsModel, RunOptions, Param, output.all = FALSE
 #' @param donor [character] id of the node which gives its parameters
 #' @param receiver [character] id of the node which receives the parameters from the donor
 #' @param default_param [numeric] vector of GR model parameters if parameters are missing from the donor
+#' @param verbose [logical] Add information message on donor and receiver
 #'
 #' @return A [numeric] [vector] with transferred parameters
 #' @export
 #'
-transferGRparams <- function(InputsModel, Param, donor, receiver, default_param = NULL) {
+transferGRparams <- function(InputsModel, Param, donor, receiver, default_param = NULL, verbose = FALSE) {
   missing_params <- setdiff(InputsModel[[receiver]]$model$indexParamUngauged,
                             InputsModel[[donor]]$model$indexParamUngauged)
+  if (verbose) {
+    message("Tranfering parameters from node '", donor, "' to node '", receiver, "'")
+  }
   if (length(missing_params) > 0) {
     if (is.null(default_param)) {
       stop("Missing parameters in transfer between nodes '",

diff --git a/R/utils.GRiwrm.R b/R/utils.GRiwrm.R
@@ -136,3 +136,37 @@ isNodeUpstream.GRiwrm <- function(x, current_node, candidate_node) {
 isNodeUpstream.GRiwrmInputsModel <- function(x, current_node, candidate_node) {
   isNodeUpstream(attr(x, "GRiwrm"), current_node, candidate_node)
 }
+
+#' Extract sub-network for calibration with ungauged nodes
+#'
+#' @inheritParams getNodeProperties
+#' @param GaugedId [character], the Id of the downstream gauged node in the
+#' ungauged cluster of sub-basins
+#'
+#' @return A [data.frame] of selected rows in `griwrm`.
+#' @noRd
+#'
+getUngaugedCluster <- function(griwrm, GaugedId) {
+  ### Set the reduced network of the basin containing ungauged nodes ###
+  # Select nodes identified with the current node as donor gauged node
+  g2 <- griwrm[getDiversionRows(griwrm, TRUE), ] # Remove duplicated by Diversions
+  donorIds <- g2$id[!is.na(g2$donor) & g2$donor == GaugedId]
+  # Remove receiver nodes that haven't GaugedId as downstream node
+  donorIds <- c(
+    GaugedId,
+    donorIds[sapply(donorIds, function(x) isNodeDownstream(griwrm, x, GaugedId))]
+  )
+  gDonor <- griwrm %>% dplyr::filter(.data$id %in% donorIds)
+  # Add upstream nodes for routing upstream flows
+  upNodes <- griwrm %>%
+    dplyr::filter(.data$down %in% gDonor$id,
+                  !.data$id %in% gDonor$id) %>%
+    dplyr::mutate(model = ifelse(!is.na(.data$model), NA, .data$model))
+  upIds <- upNodes$id
+  g <- rbind(upNodes, gDonor)
+  class(g) <- c("GRiwrm", class(g))
+  attr(g, "upIds") <- upIds
+  # Set downstream nodes
+  g$down[!g$down %in% g$id] <- NA
+  return(g)
+}
diff --git a/man/transferGRparams.Rd b/man/transferGRparams.Rd
diff --git a/tests/testthat/test-RunModel_Reservoir.R b/tests/testthat/test-RunModel_Reservoir.R
@@ -28,15 +28,6 @@ test_that("Calibration with Runmodel_Reservoir works!", {
                                  ErrorCrit_KGE2,
                                  RunOptions = RunOptions,
                                  Obs = Qobs[IndPeriod_Run,])
-  expect_message(
-    CreateInputsCrit(
-      InputsModel,
-      ErrorCrit_KGE2,
-      RunOptions = RunOptions,
-      Obs = Qobs[IndPeriod_Run,]
-    ),
-    regexp = "No observations"
-  )
 
   expect_warning(CreateCalibOptions(InputsModel), regexp = "FixedParam")
 

diff --git a/tests/testthat/test-createGRiwrm.R b/tests/testthat/test-createGRiwrm.R
@@ -68,7 +68,7 @@ test_that("Derivated ungauged node without downstream node should have derivated
   nodes <- rbind(nodes,
                  data.frame(id = "54001", down = "54032", length = 45, area = NA, model = "Diversion"))
   g <- CreateGRiwrm(nodes)
-  expect_equal(g$donor, rep("54032", 4))
+  expect_equal(g$donor, c("54032", "54032", NA, "54032"))
 })
 
 test_that("Reservoir between ungauged and gauged node should have the first downstream node as donor", {

diff --git a/tests/testthat/test-getNodeRanking.R b/tests/testthat/test-getNodeRanking.R
@@ -53,7 +53,8 @@ test_that("Impossible case detected: ungauged node with diversion to an upstream
                                            length = 20,
                                            model = "Diversion",
                                            area = NA))
-  expect_error(CreateGRiwrm(nodes_div))
+  expect_error(CreateGRiwrm(nodes_div),
+               regexp = "'54001' is located in the cluster")
 })
 
 test_that("donor of ungauged cluster is processed before sibling ungauged nodes (#155)", {