diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
index ef44dfd80d70d587df25f93b4aa58c10bcbcae24..8b39f0c3ca9085530b4daf04f83870ed976ed737 100644
--- a/R/Calibration.GRiwrmInputsModel.R
+++ b/R/Calibration.GRiwrmInputsModel.R
@@ -54,6 +54,11 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
IM <- l$InputsModel
IM$FUN_MOD <- "RunModel_Ungauged"
attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions
+ if(IM[[id]]$model$hasX4) {
+ subBasinAreas <- calcSubBasinAreas(IM)
+ donorArea <- subBasinAreas[id]
+ attr(RunOptions[[id]], "donorArea") <- donorArea
+ }
} else {
if (useUpstreamQsim && any(IM$UpstreamIsModeled)) {
# Update InputsModel$Qupstream with simulated upstream flows
@@ -79,10 +84,9 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
# Select nodes with model in the sub-network
g <- attr(IM, "GRiwrm")
Ids <- g$id[!is.na(g$donor) & g$donor == id]
- # Extract the X4 calibrated for the whole intermediate basin
if(IM[[id]]$model$hasX4) {
+ # Extract the X4 calibrated for the whole intermediate basin
X4 <- OutputsCalib[[id]]$ParamFinalR[IM[[id]]$model$iX4] # Global parameter
- subBasinAreas <- calcSubBasinAreas(IM)
}
for (uId in Ids) {
if(!IM[[uId]]$isReservoir) {
@@ -92,8 +96,10 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
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
+ OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$iX4] <- max(
+ X4 * (subBasinAreas[uId] / donorArea) ^ 0.3,
+ 0.5
+ )
}
} else {
OutputsCalib[[uId]] <- Calibration(
@@ -105,10 +111,15 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
)
}
}
+ if(useUpstreamQsim) {
+ OM_subnet <- RunModel_Ungauged(IM,
+ RunOptions[[id]],
+ OutputsCalib[[id]]$ParamFinalR,
+ output.all = TRUE)
+ OutputsModel <- c(OutputsModel, OM_subnet)
+ }
IM <- IM[[id]]
- }
-
- if(useUpstreamQsim) {
+ } else if(useUpstreamQsim) {
# Run the model for the sub-basin
OutputsModel[[id]] <- RunModel(
x = IM,
@@ -116,7 +127,6 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
Param = OutputsCalib[[id]]$ParamFinalR
)
}
-
}
return(OutputsCalib)
@@ -285,12 +295,14 @@ calcSubBasinAreas <- function(IM) {
#' <https://pastel.archives-ouvertes.fr/tel-01134990/document>
#'
#' @inheritParams airGR::RunModel
+#' @param ouput.all [logical] if `TRUE` returns the output of [RunModel.GRiwrm],
+#' returns the `OutputsModel` of the downstream node otherwise
#'
#' @inherit RunModel.GRiwrmInputsModel return return
#' @noRd
-RunModel_Ungauged <- function(InputsModel, RunOptions, Param) {
+RunModel_Ungauged <- function(InputsModel, RunOptions, Param, output.all = FALSE) {
InputsModel$FUN_MOD <- NULL
- SBVI <- sum(calcSubBasinAreas(InputsModel), na.rm = TRUE)
+ donorArea <- attr(RunOptions, "donorArea")
# Compute Param for each sub-basin
P <- lapply(InputsModel, function(IM) {
if (IM$isReservoir) {
@@ -298,12 +310,19 @@ RunModel_Ungauged <- function(InputsModel, RunOptions, Param) {
}
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
+ p[IM$model$iX4] <- max(
+ Param[IM$model$iX4] * (IM$BasinAreas[length(IM$BasinAreas)] / donorArea) ^ 0.3,
+ 0.5
+ )
}
return(p)
})
OM <- suppressMessages(
RunModel.GRiwrmInputsModel(InputsModel, attr(RunOptions, "GRiwrmRunOptions"), P)
)
- return(OM[[length(OM)]])
+ if (output.all) {
+ return(OM)
+ } else {
+ return(OM[[length(OM)]])
+ }
}
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 927862e874d466219859d60204c75b7c9e0bad2a..5a31d1c06113a3f60e9709ddfc0c3a04910ccdb4 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -422,8 +422,8 @@ getInputBV <- function(x, id, unset = NULL) {
#' @noRd
hasUngaugedNodes <- function(id, griwrm) {
nps <- getAllNodesProperties(griwrm)
- upIds <- griwrm$id[griwrm$down == id]
- upIds <- upIds[!is.na(upIds)]
+ upNodes <- griwrm[!is.na(griwrm$down) & griwrm$down == id, ]
+ upIds <- upNodes$id[upNodes$model != "Diversion"]
# No upstream nodes
if(length(upIds) == 0) return(FALSE)
# At least one upstream node is ungauged
diff --git a/tests/testthat/test-RunModel_Ungauged.R b/tests/testthat/test-RunModel_Ungauged.R
index 559053fc8351415096628f1c75c45295391d83b1..4d7a6da75f0a18da9cb160ee192a28d5563b0237 100644
--- a/tests/testthat/test-RunModel_Ungauged.R
+++ b/tests/testthat/test-RunModel_Ungauged.R
@@ -3,7 +3,7 @@ skip_on_cran()
# data set up
nodes <- loadSevernNodes()
-nodes <- nodes[!nodes$id %in% c("54002", "54057", "54095"), ]
+nodes <- nodes[nodes$id %in% c("54001", "54029", "54032"), ]
nodes[nodes$id == "54032", c("down", "length")] <- c(NA, NA)
nodes$model[nodes$id == "54029"] <- "Ungauged"
@@ -18,18 +18,30 @@ IC <- CreateInputsCrit(
FUN_CRIT = ErrorCrit_KGE2,
RunOptions = RunOptions,
Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"]],
- AprioriIds = c("54032" = "54001"),
- transfo = "sqrt",
- k = 0.15
)
CO <- CreateCalibOptions(InputsModel)
OC <- Calibration(InputsModel, RunOptions, IC, CO)
test_that("RunModel_Ungauged works for intermediate basin with ungauged station", {
- expect_true(all(sapply(OC, "[[", "CritFinal") > 0.96))
+ expect_true(all(sapply(OC, "[[", "CritFinal") > 0.95))
})
+Param <- sapply(OC, "[[", "ParamFinalR")
+OM <- RunModel(
+ InputsModel,
+ RunOptions = RunOptions,
+ Param = Param
+)
+CritValue <- ErrorCrit_KGE2(
+ InputsCrit = IC$`54032`,
+ OutputsModel = OM$`54032`
+)$CritValue
+
+# test_that("Ungauged node with gauged upstream node should works", {
+# expect_equal(OC$`54032`$CritFinal, CritValue)
+# })
+
test_that("RunModel_Ungauged works with a diversion as donor (#110)", {
nodes <- rbind(nodes,
data.frame(id = "54032", down = NA, length = NA, area = NA, model = "Diversion"))
@@ -45,9 +57,6 @@ test_that("RunModel_Ungauged works with a diversion as donor (#110)", {
FUN_CRIT = ErrorCrit_KGE2,
RunOptions = RunOptions,
Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
- AprioriIds = c("54032" = "54001"),
- transfo = "sqrt",
- k = 0.15
)
CO <- CreateCalibOptions(InputsModel)
@@ -70,8 +79,6 @@ IC <- CreateInputsCrit(
FUN_CRIT = ErrorCrit_KGE2,
RunOptions = RunOptions,
Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
- transfo = "sqrt",
- k = 0.15
)
CO <- CreateCalibOptions(InputsModel)
@@ -103,3 +110,85 @@ test_that("RunModel_Ungauged works with a diversion as upstream node (#113)", {
OCdiv <- Calibration(InputsModel, RunOptions, IC, CO)
expect_equal(OCdiv$`54032`$CritFinal, CritValue)
})
+
+test_that("RunModel_Ungauged works with a diversion as upstream node (#113)", {
+ nodes <- rbind(nodes,
+ data.frame(id = "54095", down = "54001", length = 100, area = NA, model = "Diversion"))
+ g <- CreateGRiwrm(nodes)
+ Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
+ colnames(Qobs2) <- "54095"
+ e <- setupRunModel(griwrm = g, runRunModel = FALSE, Qobs2 = Qobs2)
+ for(x in ls(e)) assign(x, get(x, e))
+ np <- getAllNodesProperties(griwrm)
+ OCdiv <- Calibration(InputsModel, RunOptions, IC, CO)
+ expect_equal(OCdiv$`54032`$CritFinal, CritValue)
+})
+
+test_that("Ungauged node with diversion outside the sub-network should work", {
+ nodes <- loadSevernNodes()
+ nodes <- nodes[!nodes$id %in% c("54002", "54057", "54029"), ]
+ nodes[nodes$id == "54032", c("down", "length")] <- c(NA, NA)
+ nodes$model[nodes$id == "54095"] <- "Ungauged"
+
+ # First without Diversion
+ g <- CreateGRiwrm(nodes)
+ e <- setupRunModel(griwrm = g, runRunModel = FALSE, Qobs2 = Qobs2)
+ for(x in ls(e)) assign(x, get(x, e))
+ np <- getAllNodesProperties(griwrm)
+
+ IC <- CreateInputsCrit(
+ InputsModel,
+ FUN_CRIT = ErrorCrit_KGE2,
+ RunOptions = RunOptions,
+ Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
+ )
+
+ CO <- CreateCalibOptions(InputsModel)
+ OC1 <- Calibration(InputsModel, RunOptions, IC, CO)
+ Param1 <- sapply(OC1, "[[", "ParamFinalR")
+ OM <- RunModel(
+ InputsModel,
+ RunOptions = RunOptions,
+ Param = Param1
+ )
+ CritValue <- ErrorCrit_KGE2(
+ InputsCrit = IC$`54032`,
+ OutputsModel = OM$`54032`
+ )$CritValue
+ # expect_equal(OC1$`54032`$CritFinal, CritValue)
+
+ # Second with Diversion with zero flow diverted for comparison
+ nodes <- rbind(nodes,
+ data.frame(id = "54095", down = "54032", length = 100,
+ area = NA, model = "Diversion"))
+ g <- CreateGRiwrm(nodes)
+ Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
+ colnames(Qobs2) <- "54095"
+ e <- setupRunModel(griwrm = g, runRunModel = FALSE, Qobs2 = Qobs2)
+ for(x in ls(e)) assign(x, get(x, e))
+ np <- getAllNodesProperties(griwrm)
+
+ IC <- CreateInputsCrit(
+ InputsModel,
+ FUN_CRIT = ErrorCrit_KGE2,
+ RunOptions = RunOptions,
+ Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
+ )
+
+ CO <- CreateCalibOptions(InputsModel)
+ OC2 <- Calibration(InputsModel, RunOptions, IC, CO)
+ expect_equal(OC2$`54001`$CritFinal, OC1$`54001`$CritFinal)
+ expect_equal(OC2$`54032`$CritFinal, OC1$`54032`$CritFinal)
+ Param2 <- sapply(OC2, "[[", "ParamFinalR")
+ expect_equal(Param2, Param1)
+ OM <- RunModel(
+ InputsModel,
+ RunOptions = RunOptions,
+ Param = Param2
+ )
+ CritValue <- ErrorCrit_KGE2(
+ InputsCrit = IC$`54032`,
+ OutputsModel = OM$`54032`
+ )$CritValue
+ # expect_equal(OC2$`54032`$CritFinal, CritValue)
+})