diff --git a/DESCRIPTION b/DESCRIPTION index 0507aac19d3ec9c21b7c817a10147dce13b9a961..a3f44ba4304c2df39e5b26a1f950bb3aa509daae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,5 +37,5 @@ VignetteBuilder: Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.0 +RoxygenNote: 7.2.1 Language: en-US diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R index 6f16b6bf02eb0c03fc438e21af3681e96682402d..1924a775093a78c6d98b945ffc2c2af30f0ce171 100644 --- a/R/Calibration.GRiwrmInputsModel.R +++ b/R/Calibration.GRiwrmInputsModel.R @@ -54,7 +54,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel, IM$FUN_MOD <- "RunModel_Ungauged" attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions } else { - if(useUpstreamQsim && any(IM$UpstreamIsRunoff)) { + if (useUpstreamQsim && any(IM$UpstreamIsRunoff)) { # Update InputsModel$Qupstream with simulated upstream flows IM <- UpdateQsimUpstream(IM, RunOptions[[id]], OutputsModel) } @@ -73,20 +73,18 @@ Calibration.GRiwrmInputsModel <- function(InputsModel, g <- attr(IM, "GRiwrm") Ids <- g$id[g$donor == id & !is.na(g$model)] # Extract the X4 calibrated for the whole intermediate basin - PS <- attr(IM[[id]], "ParamSettings") - if(PS$hasX4) { - X4 <- OutputsCalib[[id]]$ParamFinalR[PS$iX4] # Global parameter + if(IM[[id]]$model$hasX4) { + X4 <- OutputsCalib[[id]]$ParamFinalR[IM[[id]]$model$iX4] # Global parameter 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 - PS <- attr(IM[[uId]], "ParamSettings") OutputsCalib[[uId]]$ParamFinalR <- - OutputsCalib[[uId]]$ParamFinalR[PS$Indexes] - if(PS$hasX4) { - OutputsCalib[[uId]]$ParamFinalR[PS$iX4] <- + 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 } } @@ -133,6 +131,10 @@ getInputsCrit_Lavenne <- function(id, OutputsModel, InputsCrit) { # Add default velocity parameter for a priori upstream catchment AprParamR <- c(AprCelerity, AprParamR) } + if (attr(InputsCrit[[id]], "model")$hasX4) { + featMod <- attr(InputsCrit[[id]], "model") + AprParamR[featMod$iX4] <- AprParamR[featMod$iX4] * featMod$X4Ratio + } AprCrit <- ErrorCrit(InputsCrit[[AprioriId]], OutputsModel[[AprioriId]])$CritValue return(Lavenne_FUN(AprParamR, AprCrit)) } @@ -195,14 +197,7 @@ updateParameters4Ungauged <- function(GaugedId, 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)) @@ -256,10 +251,9 @@ RunModel_Ungauged <- function(InputsModel, RunOptions, Param) { SBVI <- sum(calcSubBasinAreas(InputsModel)) # 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)] / SBVI) ^ 0.3 + 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 } return(p) }) diff --git a/R/CreateInputsCrit.GRiwrmInputsModel.R b/R/CreateInputsCrit.GRiwrmInputsModel.R index 1430274b2245739e360102c6eb7ba260c2fb56c8..55394a1e20449f443256372a1746cf1456918cd9 100644 --- a/R/CreateInputsCrit.GRiwrmInputsModel.R +++ b/R/CreateInputsCrit.GRiwrmInputsModel.R @@ -61,6 +61,10 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel, stop("'AprioriIds': the node \"", AprioriIds[id], "\" is an ungauged upstream node of the node \"", id,"\"") } + if (!identical(InputsModel[[id]]$FUN_MOD, InputsModel[[AprioriIds[id]]]$FUN_MOD)) { + stop("'AprioriIds': the node \"", AprioriIds[id], + "\" must use the same hydrological model as the node \"", id,"\"") + } }) } @@ -88,6 +92,11 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel, ) attr(InputsCrit[[IM$id]], "AprioriId") <- AprioriIds[IM$id] attr(InputsCrit[[IM$id]], "AprCelerity") <- AprCelerity + attr(InputsCrit[[IM$id]], "model") <- IM$model + if (IM$model$hasX4) { + attr(InputsCrit[[IM$id]], "model")$X4Ratio <- + (tail(IM$BasinAreas, 1) / tail(InputsModel[[AprioriIds[IM$id]]]$BasinAreas, 1))^0.3 + } class(InputsCrit[[IM$id]]) <- c("InputsCritLavenneFunction", class(InputsCrit[[IM$id]])) } } diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R index e316f7a2242547ea3702083ffefdf3e2a64e0966..db34acbbe7c3464f0eaec69fd20cddcf6f8ec031 100644 --- a/R/CreateInputsModel.GRiwrm.R +++ b/R/CreateInputsModel.GRiwrm.R @@ -194,10 +194,13 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) { # Add the model function InputsModel$FUN_MOD <- FUN_MOD + featModel <- .GetFeatModel(InputsModel) InputsModel$isUngauged <- griwrm$model[griwrm$id == id] == "Ungauged" InputsModel$gaugedId <- griwrm$donor[griwrm$id == id] InputsModel$hasUngaugedNodes <- hasUngaugedNodes(id, griwrm) - + InputsModel$model <- list(indexParamUngauged = ifelse(inherits(InputsModel, "SD"), 0, 1) + seq.int(featModel$NbParam), + hasX4 = grepl("RunModel_GR[456][HJ]", FUN_MOD), + iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4)) return(InputsModel) } @@ -282,3 +285,24 @@ hasUngaugedNodes <- function(id, griwrm) { } return(FALSE) } + + +#' function to extract model features partially copied from airGR:::.GetFeatModel +#' @noRd +.GetFeatModel <- function(InputsModel) { + path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR") + FeatMod <- read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE) + NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR", + yes = paste("RunModel", FeatMod$NameMod, sep = "_"), + no = FeatMod$NameMod) + IdMod <- which(sapply(NameFunMod, FUN = function(x) identical(InputsModel$FUN_MOD, x))) + if (length(IdMod) < 1) { + stop("'FUN_MOD' must be one of ", paste(NameFunMod, collapse = ", ")) + } + FeatMod <- as.list(FeatMod[IdMod, ]) + FeatMod$IsSD <- inherits(InputsModel, "SD") + if (FeatMod$IsSD) { + FeatMod$NbParam <- FeatMod$NbParam + 1 + } + return(FeatMod) +} diff --git a/tests/testthat/test-CreateInputsCrit.R b/tests/testthat/test-CreateInputsCrit.R index f5b55d689867bb9269ddbe3f14cd42d8379cc401..3d02d28056fedb515635e283026b6655c0f7b45b 100644 --- a/tests/testthat/test-CreateInputsCrit.R +++ b/tests/testthat/test-CreateInputsCrit.R @@ -109,6 +109,18 @@ test_that("Lavenne criterion: wrong sub-catchment order should throw error", { ) }) +test_that("Lavenne criterion: current node and a priori node must use the same model", { + InputsModel[["54032"]]$FUN_MOD <- RunModel_GR6J + expect_error( + CreateInputsCrit(InputsModel = InputsModel, + RunOptions = RunOptions, + Obs = Qobs[IndPeriod_Run,], + AprioriIds = AprioriIds, + transfo = "sqrt"), + regexp = "must use the same hydrological model" + ) +}) + test_that("Ungauged node as Apriori node should throw an error", { nodes$model[nodes$gauge_id == "54001"] <- "Ungauged" griwrm <- CreateGRiwrm(