diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 599ec847b10e5b32ad7701752284926fa025ca7d..73d31f89ad7f752944b114e78ad61f76eab5b84a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -9,6 +9,7 @@ default:
- echo "setwd(\"$(pwd)\")" > .Rprofile
- PATH=~/R/sources/R-${R_VERSION}/bin:$PATH
- rename "s/${R_VERSION}.airGR/airGR/" *.tar.gz
+ - R -e 'chooseCRANmirror(graphics = FALSE, ind = 1); pkg <- "caRamel"; pkgInst <- installed.packages()[, "Package"]; pkgMiss <- setdiff(pkg, pkgInst); if (length(pkgMiss) > 0) install.packages(pkgMiss)'
.update_packages:
stage: update_packages
diff --git a/.regressionignore b/.regressionignore
index 1dd20c9de2d828cd149545d7fbc7da4326b40af2..401394a89412d17c6ec49dcd016599597a8dc38c 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -3,28 +3,6 @@
# The format of this file is: 5 lines of comments followed by one line by
# ignored variable : [Topic]<SPACE>[Variable].
# Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
-RunModel_GR1A BasinObs
-RunModel_GR1A ConvertFun
-RunModel_GR1A NewTabSeries
-RunModel_GR1A NewTimeFormat
-RunModel_GR1A OutputsModel
-RunModel_GR1A TabSeries
-RunModel_GR1A TimeFormat
-RunModel_GR1A YearFirstMonth
-RunModel_GR2M BasinObs
-RunModel_GR2M ConvertFun
-RunModel_GR2M NewTabSeries
-RunModel_GR2M NewTimeFormat
-RunModel_GR2M OutputsModel
-RunModel_GR2M RunOptions
-RunModel_GR1A OutputsModel
-Calibration_Michel CalibOptions
-Calibration CalibOptions
-CreateCalibOptions CalibOptions
-# New version of the SeriesAggreg function
-RunModel_GR2M TabSeries
-RunModel_GR2M TimeFormat
-SeriesAggreg BasinInfo
-SeriesAggreg BasinObs
-SeriesAggreg NewTabSeries
+
+
diff --git a/DESCRIPTION b/DESCRIPTION
index 33863672b2605bf5d9b0db4960404bc041052b1b..bc62189af5a743d8ff5a7fc95ef506146e5c48c7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.9.27
-Date: 2021-01-18
+Version: 1.6.10.4
+Date: 2021-01-29
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
@@ -28,7 +28,8 @@ Imports:
utils
Suggests:
knitr, rmarkdown,
- coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+ caRamel, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+ GGally, ggplot2,
testthat
Description: Hydrological modelling tools developed at INRAE-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR5H, GR4J, GR5J, GR6J, GR2M, GR1A), a snow accumulation and melt model (CemaNeige) and the associated functions for their calibration and evaluation. Use help(airGR) for package description and references.
License: GPL-2
diff --git a/NEWS.md b/NEWS.md
index bb266368a7502df97bc9525d3f4c593f56be260e..ba381fe7c89c0660fab470bc9c53885b50837a54 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,30 @@
+### 1.6.10.4 Release Notes (2021-01-29)
+
+#### New features
+
+- Added a section 'param_optim' vignette to explain how to manage with multiobjective optimization using the 'CaRamel' package. ([#61](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/61))
+
+
+#### Major user-visible changes
+
+- `Imax()` now returns an error message when `IndPeriod_Run` doesn't select 24 hours by day, instead of `numeric(0)`. ([#92](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/92))
+
+
+#### Minor user-visible changes
+
+- Fixed warning returned by GCC Fortran when compiling `frun_GR5H.f90`. ([#93](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/93))
+
+
+#### CRAN-compatibility updates
+
+- Coerce `POSIXlt` dates into character in `RunModel_GR1A()` example and in `SeriesAggreg()` tests in order to avoid bad subsetting on time series due to mixing UTC and local time on macOS flavors. ([#94](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/94))
+
+____________________________________________________________________________________
+
+
### 1.6.9.27 Release Notes (2021-01-18)
#### New features
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index cfa85ac63c9c52b86e691caa3f7dfa46b3e4fcba..69a3da22b2b312e01e80b37a0103d8fe1663dba8 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -219,10 +219,10 @@ Calibration_Michel <- function(InputsModel,
PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
if (PotentialCandidateT[1, I] < RangesT[1, I] ) {
- PotentialCandidateT[1,I] <- RangesT[1, I]
+ PotentialCandidateT[1, I] <- RangesT[1, I]
}
if (PotentialCandidateT[1, I] > RangesT[2, I]) {
- PotentialCandidateT[1,I] <- RangesT[2,I]
+ PotentialCandidateT[1, I] <- RangesT[2, I]
}
##We_check_the_set_is_not_outside_the_range_of_possible_values
if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 7b64ef2167fd56cf19479abe726c9e852218b183..78302745a0d96af3512c53068cacf531d964a566 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -21,7 +21,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (!is.logical(IsSD) | length(IsSD) != 1L) {
stop("'IsSD' must be a logical of length 1")
}
- ##check_FUN_MOD
+ ## check FUN_MOD
BOOL <- FALSE
if (identical(FUN_MOD, RunModel_GR4H)) {
@@ -87,7 +87,7 @@ CreateCalibOptions <- function(FUN_MOD,
return(NULL)
}
- ##check_FUN_CALIB
+ ## check FUN_CALIB
BOOL <- FALSE
if (identical(FUN_CALIB, Calibration_Michel)) {
@@ -100,9 +100,9 @@ CreateCalibOptions <- function(FUN_MOD,
}
- ##check_FUN_TRANSFO
+ ## check FUN_TRANSFO
if (is.null(FUN_TRANSFO)) {
- ##_set_FUN1
+ ## set FUN1
if (identical(FUN_MOD, RunModel_GR4H) |
identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
FUN_GR <- TransfoParam_GR4H
@@ -140,17 +140,17 @@ CreateCalibOptions <- function(FUN_MOD,
stop("'FUN_GR' was not found")
return(NULL)
}
- ##_set_FUN2
+ ## set FUN2
if (IsHyst) {
FUN_SNOW <- TransfoParam_CemaNeigeHyst
} else {
FUN_SNOW <- TransfoParam_CemaNeige
}
- ##_set_FUN_LAG
+ ## set FUN_LAG
if (IsSD) {
FUN_LAG <- TransfoParam_Lag
}
- ##_set_FUN_TRANSFO
+ ## set FUN_TRANSFO
if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
if (!IsSD) {
FUN_TRANSFO <- FUN_GR
@@ -179,7 +179,7 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4) ], Direction)
+ ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
@@ -198,7 +198,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (NParam <= 3) {
ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
} else {
- ParamOut[, 1:(NParam - 2)] <- FUN_GR( ParamIn[, 1:(NParam - 2)], Direction)
+ ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
if (!Bool) {
@@ -215,9 +215,9 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 2:(NParam - 4) ] <- FUN_GR( ParamIn[, 2:(NParam - 4) ], Direction)
- ParamOut[, (NParam - 3):NParam] <- FUN_SNOW( ParamIn[, (NParam - 3):NParam], Direction)
- ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
+ ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -235,9 +235,9 @@ CreateCalibOptions <- function(FUN_MOD,
if (NParam <= 3) {
ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
} else {
- ParamOut[, 2:(NParam - 2)] <- FUN_GR( ParamIn[, 2:(NParam - 2)], Direction)
+ ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
}
- ParamOut[, (NParam - 1):NParam] <- FUN_SNOW( ParamIn[, (NParam - 1):NParam], Direction)
+ ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
@@ -252,7 +252,7 @@ CreateCalibOptions <- function(FUN_MOD,
return(NULL)
}
- ##NParam
+ ## NParam
if ("GR4H" %in% ObjectClass) {
NParam <- 4
}
@@ -299,7 +299,7 @@ CreateCalibOptions <- function(FUN_MOD,
NParam <- NParam + 1
}
- ##check_FixedParam
+ ## check FixedParam
if (is.null(FixedParam)) {
FixedParam <- rep(NA, NParam)
} else {
@@ -317,10 +317,10 @@ CreateCalibOptions <- function(FUN_MOD,
}
}
- ##check_SearchRanges
+ ## check SearchRanges
if (is.null(SearchRanges)) {
- ParamT <- matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
- ncol = NParam, byrow = TRUE)
+ ParamT <- matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
+ ncol = NParam, byrow = TRUE)
SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
} else {
@@ -341,7 +341,7 @@ CreateCalibOptions <- function(FUN_MOD,
}
}
- ##check_StartParamList_and_StartParamDistrib__default_values
+ ## check StartParamList and StartParamDistrib default values
if (("HBAN" %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
if ("GR4H" %in% ObjectClass) {
ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
@@ -351,12 +351,12 @@ CreateCalibOptions <- function(FUN_MOD,
if (("GR5H" %in% ObjectClass) & ("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34,
+3.74, -0.41, +4.78, -8.94, -3.33,
- +4.29, +0.16, +5.39, -7.39, +3.33), ncol=5, byrow = TRUE);
+ +4.29, +0.16, +5.39, -7.39, +3.33), ncol = 5, byrow = TRUE)
}
if (("GR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
+3.62, -0.19, +4.80, -9.00, -6.31,
- +4.01, -0.04, +5.43, -7.53, -5.33), ncol=5, byrow = TRUE);
+ +4.01, -0.04, +5.43, -7.53, -5.33), ncol = 5, byrow = TRUE)
}
if ("GR4J" %in% ObjectClass) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
@@ -399,12 +399,12 @@ CreateCalibOptions <- function(FUN_MOD,
if (("CemaNeigeGR5H" %in% ObjectClass) & ("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
+3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
- +4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE);
+ +4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
}
if (("CemaNeigeGR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
+3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
- +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE);
+ +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
}
if ("CemaNeigeGR4J" %in% ObjectClass) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
@@ -440,7 +440,7 @@ CreateCalibOptions <- function(FUN_MOD,
}
- ##check_StartParamList_and_StartParamDistrib__format
+ ## check StartParamList and StartParamDistrib format
if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
if (!is.matrix(StartParamList)) {
stop("'StartParamList' must be a matrix")
diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 6fefda3fffb27b90dc7c23c947de8d7f26c85769..d09031ca95dc9833571fe46cb383673f10b816a8 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -8,333 +8,333 @@ CreateInputsModel <- function(FUN_MOD,
verbose = TRUE) {
- ObjectClass <- NULL
+ ObjectClass <- NULL
- FUN_MOD <- match.fun(FUN_MOD)
+ FUN_MOD <- match.fun(FUN_MOD)
- ##check_FUN_MOD
- BOOL <- FALSE
- if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H)) {
- ObjectClass <- c(ObjectClass, "hourly", "GR")
+ ##check_FUN_MOD
+ BOOL <- FALSE
+ if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H)) {
+ ObjectClass <- c(ObjectClass, "hourly", "GR")
- TimeStep <- as.integer(60 * 60)
+ TimeStep <- as.integer(60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR4J) |
- identical(FUN_MOD, RunModel_GR5J) |
- identical(FUN_MOD, RunModel_GR6J)) {
- ObjectClass <- c(ObjectClass, "daily", "GR")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_GR4J) |
+ identical(FUN_MOD, RunModel_GR5J) |
+ identical(FUN_MOD, RunModel_GR6J)) {
+ ObjectClass <- c(ObjectClass, "daily", "GR")
- TimeStep <- as.integer(24 * 60 * 60)
+ TimeStep <- as.integer(24 * 60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR2M)) {
- ObjectClass <- c(ObjectClass, "GR", "monthly")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_GR2M)) {
+ ObjectClass <- c(ObjectClass, "GR", "monthly")
- TimeStep <- as.integer(c(28, 29, 30, 31) * 24 * 60 * 60)
+ TimeStep <- as.integer(c(28, 29, 30, 31) * 24 * 60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_GR1A)) {
- ObjectClass <- c(ObjectClass, "GR", "yearly")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_GR1A)) {
+ ObjectClass <- c(ObjectClass, "GR", "yearly")
- TimeStep <- as.integer(c(365, 366) * 24 * 60 * 60)
+ TimeStep <- as.integer(c(365, 366) * 24 * 60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeige)) {
- ObjectClass <- c(ObjectClass, "daily", "CemaNeige")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_CemaNeige)) {
+ ObjectClass <- c(ObjectClass, "daily", "CemaNeige")
- TimeStep <- as.integer(24 * 60 * 60)
+ TimeStep <- as.integer(24 * 60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
- identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
- ObjectClass <- c(ObjectClass, "daily", "GR", "CemaNeige")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
+ identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+ identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+ ObjectClass <- c(ObjectClass, "daily", "GR", "CemaNeige")
- TimeStep <- as.integer(24 * 60 * 60)
+ TimeStep <- as.integer(24 * 60 * 60)
- BOOL <- TRUE
- }
- if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
- ObjectClass <- c(ObjectClass, "hourly", "GR", "CemaNeige")
+ BOOL <- TRUE
+ }
+ if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
+ ObjectClass <- c(ObjectClass, "hourly", "GR", "CemaNeige")
+
+ TimeStep <- as.integer(60 * 60)
- TimeStep <- as.integer(60 * 60)
+ BOOL <- TRUE
+ }
+ if (!BOOL) {
+ stop("incorrect 'FUN_MOD' for use in 'CreateInputsModel'")
+ }
- BOOL <- TRUE
+ ##check_arguments
+ if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
+ if (is.null(DatesR)) {
+ stop("'DatesR' is missing")
}
- if (!BOOL) {
- stop("incorrect 'FUN_MOD' for use in 'CreateInputsModel'")
+ if (!"POSIXlt" %in% class(DatesR) & !"POSIXct" %in% class(DatesR)) {
+ stop("'DatesR' must be defined as 'POSIXlt' or 'POSIXct'")
}
-
- ##check_arguments
- if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
- if (is.null(DatesR)) {
- stop("'DatesR' is missing")
- }
- if (!"POSIXlt" %in% class(DatesR) & !"POSIXct" %in% class(DatesR)) {
- stop("'DatesR' must be defined as 'POSIXlt' or 'POSIXct'")
- }
- if (!"POSIXlt" %in% class(DatesR)) {
- DatesR <- as.POSIXlt(DatesR)
- }
- if (!difftime(tail(DatesR, 1), tail(DatesR, 2), units = "secs")[[1]] %in% TimeStep) {
- TimeStepName <- grep("hourly|daily|monthly|yearly", ObjectClass, value = TRUE)
- stop(paste0("the time step of the model inputs must be ", TimeStepName, "\n"))
- }
- if (any(duplicated(DatesR))) {
- stop("'DatesR' must not include duplicated values")
- }
- LLL <- length(DatesR)
+ if (!"POSIXlt" %in% class(DatesR)) {
+ DatesR <- as.POSIXlt(DatesR)
}
- if ("GR" %in% ObjectClass) {
- if (is.null(Precip)) {
- stop("Precip is missing")
- }
- if (is.null(PotEvap)) {
- stop("'PotEvap' is missing")
- }
- if (!is.vector(Precip) | !is.vector(PotEvap)) {
- stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+ if (!difftime(tail(DatesR, 1), tail(DatesR, 2), units = "secs")[[1]] %in% TimeStep) {
+ TimeStepName <- grep("hourly|daily|monthly|yearly", ObjectClass, value = TRUE)
+ stop(paste0("the time step of the model inputs must be ", TimeStepName, "\n"))
+ }
+ if (any(duplicated(DatesR))) {
+ stop("'DatesR' must not include duplicated values")
+ }
+ LLL <- length(DatesR)
+ }
+ if ("GR" %in% ObjectClass) {
+ if (is.null(Precip)) {
+ stop("Precip is missing")
+ }
+ if (is.null(PotEvap)) {
+ stop("'PotEvap' is missing")
+ }
+ if (!is.vector(Precip) | !is.vector(PotEvap)) {
+ stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+ }
+ if (!is.numeric(Precip) | !is.numeric(PotEvap)) {
+ stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+ }
+ if (length(Precip) != LLL | length(PotEvap) != LLL) {
+ stop("'Precip', 'PotEvap' and 'DatesR' must have the same length")
+ }
+ }
+ if ("CemaNeige" %in% ObjectClass) {
+ if (is.null(Precip)) {
+ stop("'Precip' is missing")
+ }
+ if (is.null(TempMean)) {
+ stop("'TempMean' is missing")
+ }
+ if (!is.vector(Precip) | !is.vector(TempMean)) {
+ stop("'Precip' and 'TempMean' must be vectors of numeric values")
+ }
+ if (!is.numeric(Precip) | !is.numeric(TempMean)) {
+ stop("'Precip' and 'TempMean' must be vectors of numeric values")
+ }
+ if (length(Precip) != LLL | length(TempMean) != LLL) {
+ stop("'Precip', 'TempMean' and 'DatesR' must have the same length")
+ }
+ if (is.null(TempMin) != is.null(TempMax)) {
+ stop("'TempMin' and 'TempMax' must be both defined if not null")
+ }
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ if (!is.vector(TempMin) | !is.vector(TempMax)) {
+ stop("'TempMin' and 'TempMax' must be vectors of numeric values")
}
- if (!is.numeric(Precip) | !is.numeric(PotEvap)) {
- stop("'Precip' and 'PotEvap' must be vectors of numeric values")
+ if (!is.numeric(TempMin) | !is.numeric(TempMax)) {
+ stop("'TempMin' and 'TempMax' must be vectors of numeric values")
}
- if (length(Precip) != LLL | length(PotEvap) != LLL) {
- stop("'Precip', 'PotEvap' and 'DatesR' must have the same length")
+ if (length(TempMin) != LLL | length(TempMax) != LLL) {
+ stop("'TempMin', 'TempMax' and 'DatesR' must have the same length")
}
}
- if ("CemaNeige" %in% ObjectClass) {
- if (is.null(Precip)) {
- stop("'Precip' is missing")
- }
- if (is.null(TempMean)) {
- stop("'TempMean' is missing")
- }
- if (!is.vector(Precip) | !is.vector(TempMean)) {
- stop("'Precip' and 'TempMean' must be vectors of numeric values")
+ if (!is.null(HypsoData)) {
+ if (!is.vector(HypsoData)) {
+ stop("'HypsoData' must be a vector of numeric values if not null")
}
- if (!is.numeric(Precip) | !is.numeric(TempMean)) {
- stop("'Precip' and 'TempMean' must be vectors of numeric values")
+ if (!is.numeric(HypsoData)) {
+ stop("'HypsoData' must be a vector of numeric values if not null")
}
- if (length(Precip) != LLL | length(TempMean) != LLL) {
- stop("'Precip', 'TempMean' and 'DatesR' must have the same length")
+ if (length(HypsoData) != 101) {
+ stop("'HypsoData' must be of length 101 if not null")
}
- if (is.null(TempMin) != is.null(TempMax)) {
- stop("'TempMin' and 'TempMax' must be both defined if not null")
+ if (sum(is.na(HypsoData)) != 0 & sum(is.na(HypsoData)) != 101) {
+ stop("'HypsoData' must not contain any NA if not null")
}
- if (!is.null(TempMin) & !is.null(TempMax)) {
- if (!is.vector(TempMin) | !is.vector(TempMax)) {
- stop("'TempMin' and 'TempMax' must be vectors of numeric values")
- }
- if (!is.numeric(TempMin) | !is.numeric(TempMax)) {
- stop("'TempMin' and 'TempMax' must be vectors of numeric values")
- }
- if (length(TempMin) != LLL | length(TempMax) != LLL) {
- stop("'TempMin', 'TempMax' and 'DatesR' must have the same length")
- }
+ }
+ if (!is.null(ZInputs)) {
+ if (length(ZInputs) != 1) {
+ stop("'ZInputs' must be a single numeric value if not null")
}
- if (!is.null(HypsoData)) {
- if (!is.vector(HypsoData)) {
- stop("'HypsoData' must be a vector of numeric values if not null")
- }
- if (!is.numeric(HypsoData)) {
- stop("'HypsoData' must be a vector of numeric values if not null")
- }
- if (length(HypsoData) != 101) {
- stop("'HypsoData' must be of length 101 if not null")
- }
- if (sum(is.na(HypsoData)) != 0 & sum(is.na(HypsoData)) != 101) {
- stop("'HypsoData' must not contain any NA if not null")
- }
+ if (is.na(ZInputs) | !is.numeric(ZInputs)) {
+ stop("'ZInputs' must be a single numeric value if not null")
}
- if (!is.null(ZInputs)) {
- if (length(ZInputs) != 1) {
- stop("'ZInputs' must be a single numeric value if not null")
- }
- if (is.na(ZInputs) | !is.numeric(ZInputs)) {
- stop("'ZInputs' must be a single numeric value if not null")
- }
+ }
+ if (is.null(HypsoData)) {
+ if (verbose) {
+ warning("'HypsoData' is missing: a single layer is used and no extrapolation is made")
}
- if (is.null(HypsoData)) {
- if (verbose) {
- warning("'HypsoData' is missing: a single layer is used and no extrapolation is made")
- }
- HypsoData <- as.numeric(rep(NA, 101))
- ZInputs <- as.numeric(NA)
- NLayers <- as.integer(1)
+ HypsoData <- as.numeric(rep(NA, 101))
+ ZInputs <- as.numeric(NA)
+ NLayers <- as.integer(1)
+ }
+ if (is.null(ZInputs)) {
+ if (verbose & !identical(HypsoData, as.numeric(rep(NA, 101)))) {
+ warning("'ZInputs' is missing: HypsoData[51] is used")
}
- if (is.null(ZInputs)) {
- if (verbose & !identical(HypsoData, as.numeric(rep(NA, 101)))) {
- warning("'ZInputs' is missing: HypsoData[51] is used")
- }
- ZInputs <- HypsoData[51L]
- }
- if (NLayers <= 0) {
- stop("'NLayers' must be a positive integer value")
- }
- if (NLayers != as.integer(NLayers)) {
- warning("Coerce 'NLayers' to be of integer type (", NLayers, ": ", as.integer(NLayers), ")")
- NLayers <- as.integer(NLayers)
- }
+ ZInputs <- HypsoData[51L]
}
-
- ## check semi-distributed mode
- if (!is.null(Qupstream) & !is.null(LengthHydro) & !is.null(BasinAreas)) {
- ObjectClass <- c(ObjectClass, "SD")
- } else if (verbose & !all(c(is.null(Qupstream), is.null(LengthHydro), is.null(BasinAreas)))) {
- warning("Missing argument: 'Qupstream', 'LengthHydro' and 'BasinAreas' must all be set to run in a semi-distributed mode. The lumped mode will be used")
+ if (NLayers <= 0) {
+ stop("'NLayers' must be a positive integer value")
}
- if ("SD" %in% ObjectClass) {
- if (!("daily" %in% ObjectClass) & !("hourly" %in% ObjectClass)) {
- stop("Only daily and hourly time steps can be used in a semi-distributed mode")
- }
- if (!is.matrix(Qupstream) | !is.numeric(Qupstream)) {
- stop("'Qupstream' must be a matrice of numeric values")
- }
- if (!is.vector(LengthHydro) | !is.vector(BasinAreas) | !is.numeric(LengthHydro) | !is.numeric(BasinAreas)) {
- stop("'LengthHydro' and 'BasinAreas' must be vectors of numeric values")
- }
- if (ncol(Qupstream) != length(LengthHydro)) {
- stop("'Qupstream' number of columns and 'LengthHydro' length must be equal")
- }
- if (length(LengthHydro) + 1 != length(BasinAreas)) {
- stop("'BasinAreas' must have one more element than 'LengthHydro'")
- }
- if (nrow(Qupstream) != LLL) {
- stop("'Qupstream' must have same number of rows as 'DatesR' length")
- }
- if(any(is.na(Qupstream))) {
- stop("'Qupstream' cannot contain any NA value")
- }
+ if (NLayers != as.integer(NLayers)) {
+ warning("Coerce 'NLayers' to be of integer type (", NLayers, ": ", as.integer(NLayers), ")")
+ NLayers <- as.integer(NLayers)
+ }
+ }
+
+ ## check semi-distributed mode
+ if (!is.null(Qupstream) & !is.null(LengthHydro) & !is.null(BasinAreas)) {
+ ObjectClass <- c(ObjectClass, "SD")
+ } else if (verbose & !all(c(is.null(Qupstream), is.null(LengthHydro), is.null(BasinAreas)))) {
+ warning("Missing argument: 'Qupstream', 'LengthHydro' and 'BasinAreas' must all be set to run in a semi-distributed mode. The lumped mode will be used")
+ }
+ if ("SD" %in% ObjectClass) {
+ if (!("daily" %in% ObjectClass) & !("hourly" %in% ObjectClass)) {
+ stop("Only daily and hourly time steps can be used in a semi-distributed mode")
+ }
+ if (!is.matrix(Qupstream) | !is.numeric(Qupstream)) {
+ stop("'Qupstream' must be a matrice of numeric values")
+ }
+ if (!is.vector(LengthHydro) | !is.vector(BasinAreas) | !is.numeric(LengthHydro) | !is.numeric(BasinAreas)) {
+ stop("'LengthHydro' and 'BasinAreas' must be vectors of numeric values")
}
+ if (ncol(Qupstream) != length(LengthHydro)) {
+ stop("'Qupstream' number of columns and 'LengthHydro' length must be equal")
+ }
+ if (length(LengthHydro) + 1 != length(BasinAreas)) {
+ stop("'BasinAreas' must have one more element than 'LengthHydro'")
+ }
+ if (nrow(Qupstream) != LLL) {
+ stop("'Qupstream' must have same number of rows as 'DatesR' length")
+ }
+ if(any(is.na(Qupstream))) {
+ stop("'Qupstream' cannot contain any NA value")
+ }
+ }
- ##check_NA_values
- BOOL_NA <- rep(FALSE, length(DatesR))
+ ##check_NA_values
+ BOOL_NA <- rep(FALSE, length(DatesR))
- if ("GR" %in% ObjectClass) {
- BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
- if (sum(BOOL_NA_TMP) != 0) {
- BOOL_NA <- BOOL_NA | BOOL_NA_TMP
- if (verbose) {
- warning("Values < 0 or NA values detected in 'Precip' series")
- }
+ if ("GR" %in% ObjectClass) {
+ BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
+ if (sum(BOOL_NA_TMP) != 0) {
+ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+ if (verbose) {
+ warning("Values < 0 or NA values detected in 'Precip' series")
}
- BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap)
- if (sum(BOOL_NA_TMP) != 0) {
- BOOL_NA <- BOOL_NA | BOOL_NA_TMP
- if (verbose) {
- warning("Values < 0 or NA values detected in 'PotEvap' series")
- }
+ }
+ BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap)
+ if (sum(BOOL_NA_TMP) != 0) {
+ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+ if (verbose) {
+ warning("Values < 0 or NA values detected in 'PotEvap' series")
}
}
- if ("CemaNeige" %in% ObjectClass) {
- BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
+ }
+ if ("CemaNeige" %in% ObjectClass) {
+ BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
+ if (sum(BOOL_NA_TMP) != 0) {
+ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+ if (verbose) {
+ warning("Values < 0 or NA values detected in 'Precip' series")
+ }
+ }
+ BOOL_NA_TMP <- (TempMean < (-150)) | is.na(TempMean)
+ if (sum(BOOL_NA_TMP) != 0) {
+ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
+ if (verbose) {
+ warning("Values < -150 or NA values detected in 'TempMean' series")
+ }
+ }
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ BOOL_NA_TMP <- (TempMin < (-150)) | is.na(TempMin)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
- warning("Values < 0 or NA values detected in 'Precip' series")
+ warning("Values < -150 or NA values detected in 'TempMin' series")
}
}
- BOOL_NA_TMP <- (TempMean < (-150)) | is.na(TempMean)
+ BOOL_NA_TMP <- (TempMax < (-150)) | is.na(TempMax)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
- warning("Values < -150 or NA values detected in 'TempMean' series")
- }
- }
- if (!is.null(TempMin) & !is.null(TempMax)) {
- BOOL_NA_TMP <- (TempMin < (-150)) | is.na(TempMin)
- if (sum(BOOL_NA_TMP) != 0) {
- BOOL_NA <- BOOL_NA | BOOL_NA_TMP
- if (verbose) {
- warning("Values < -150 or NA values detected in 'TempMin' series")
- }
- }
- BOOL_NA_TMP <- (TempMax < (-150)) | is.na(TempMax)
- if (sum(BOOL_NA_TMP) != 0) {
- BOOL_NA <- BOOL_NA | BOOL_NA_TMP
- if (verbose) {
- warning("Values < -150 or NA values detected in 'TempMax' series")
- }
+ warning("Values < -150 or NA values detected in 'TempMax' series")
}
}
}
- if (sum(BOOL_NA) != 0) {
- WTxt <- NULL
- WTxt <- paste(WTxt, "\t Missing values are not allowed in 'InputsModel'", sep = "")
-
- Select <- (max(which(BOOL_NA)) + 1):length(BOOL_NA)
-
- if (Select[1L] > Select[2L]) {
- stop("time series could not be trunced since missing values were detected at the last time-step")
- }
- if ("GR" %in% ObjectClass) {
- Precip <- Precip[Select]
- PotEvap <- PotEvap[Select]
- }
- if ("CemaNeige" %in% ObjectClass) {
- Precip <- Precip[Select]
- TempMean <- TempMean[Select]
- if (!is.null(TempMin) & !is.null(TempMax)) {
- TempMin <- TempMin[Select]
- TempMax <- TempMax[Select]
- }
- }
+ }
+ if (sum(BOOL_NA) != 0) {
+ WTxt <- NULL
+ WTxt <- paste(WTxt, "\t Missing values are not allowed in 'InputsModel'", sep = "")
- DatesR <- DatesR[Select]
+ Select <- (max(which(BOOL_NA)) + 1):length(BOOL_NA)
- WTxt <- paste0(WTxt, "\t -> data were trunced to keep the most recent available time-steps")
- WTxt <- paste0(WTxt, "\t -> ", length(Select), " time-steps were kept")
-
- if (!is.null(WTxt) & verbose) {
- warning(WTxt)
- }
+ if (Select[1L] > Select[2L]) {
+ stop("time series could not be trunced since missing values were detected at the last time-step")
+ }
+ if ("GR" %in% ObjectClass) {
+ Precip <- Precip[Select]
+ PotEvap <- PotEvap[Select]
}
-
-
- ##DataAltiExtrapolation_Valery
if ("CemaNeige" %in% ObjectClass) {
- RESULT <- DataAltiExtrapolation_Valery(DatesR = DatesR,
- Precip = Precip, PrecipScale = PrecipScale,
- TempMean = TempMean, TempMin = TempMin, TempMax = TempMax,
- ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers,
- verbose = verbose)
- if (verbose) {
- if (NLayers == 1) {
- message("input series were successfully created on 1 elevation layer for use by CemaNeige")
- } else {
- message( "input series were successfully created on ", NLayers, " elevation layers for use by CemaNeige")
- }
+ Precip <- Precip[Select]
+ TempMean <- TempMean[Select]
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ TempMin <- TempMin[Select]
+ TempMax <- TempMax[Select]
}
}
+ DatesR <- DatesR[Select]
- ##Create_InputsModel
- InputsModel <- list(DatesR = DatesR)
- if ("GR" %in% ObjectClass) {
- InputsModel <- c(InputsModel, list(Precip = as.double(Precip), PotEvap = as.double(PotEvap)))
- }
- if ("CemaNeige" %in% ObjectClass) {
- InputsModel <- c(InputsModel, list(LayerPrecip = RESULT$LayerPrecip,
- LayerTempMean = RESULT$LayerTempMean,
- LayerFracSolidPrecip = RESULT$LayerFracSolidPrecip,
- ZLayers = RESULT$ZLayers))
+ WTxt <- paste0(WTxt, "\t -> data were trunced to keep the most recent available time-steps")
+ WTxt <- paste0(WTxt, "\t -> ", length(Select), " time-steps were kept")
+
+ if (!is.null(WTxt) & verbose) {
+ warning(WTxt)
}
- if ("SD" %in% ObjectClass) {
- InputsModel <- c(InputsModel, list(Qupstream = Qupstream,
- LengthHydro = LengthHydro,
- BasinAreas = BasinAreas))
+ }
+
+
+ ##DataAltiExtrapolation_Valery
+ if ("CemaNeige" %in% ObjectClass) {
+ RESULT <- DataAltiExtrapolation_Valery(DatesR = DatesR,
+ Precip = Precip, PrecipScale = PrecipScale,
+ TempMean = TempMean, TempMin = TempMin, TempMax = TempMax,
+ ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers,
+ verbose = verbose)
+ if (verbose) {
+ if (NLayers == 1) {
+ message("input series were successfully created on 1 elevation layer for use by CemaNeige")
+ } else {
+ message( "input series were successfully created on ", NLayers, " elevation layers for use by CemaNeige")
+ }
}
-
- class(InputsModel) <- c("InputsModel", ObjectClass)
-
- return(InputsModel)
+ }
+
+
+ ##Create_InputsModel
+ InputsModel <- list(DatesR = DatesR)
+ if ("GR" %in% ObjectClass) {
+ InputsModel <- c(InputsModel, list(Precip = as.double(Precip), PotEvap = as.double(PotEvap)))
+ }
+ if ("CemaNeige" %in% ObjectClass) {
+ InputsModel <- c(InputsModel, list(LayerPrecip = RESULT$LayerPrecip,
+ LayerTempMean = RESULT$LayerTempMean,
+ LayerFracSolidPrecip = RESULT$LayerFracSolidPrecip,
+ ZLayers = RESULT$ZLayers))
+ }
+ if ("SD" %in% ObjectClass) {
+ InputsModel <- c(InputsModel, list(Qupstream = Qupstream,
+ LengthHydro = LengthHydro,
+ BasinAreas = BasinAreas))
+ }
+
+ class(InputsModel) <- c("InputsModel", ObjectClass)
+
+ return(InputsModel)
diff --git a/R/DataAltiExtrapolation_Valery.R b/R/DataAltiExtrapolation_Valery.R
index d1e85cca51719f4e50dd12c3109e2e6780294637..24b1b48ed13e22130625c6b4b66d8136d1eccfec 100644
--- a/R/DataAltiExtrapolation_Valery.R
+++ b/R/DataAltiExtrapolation_Valery.R
@@ -3,549 +3,171 @@ DataAltiExtrapolation_Valery <- function(DatesR,
TempMean, TempMin = NULL, TempMax = NULL,
ZInputs, HypsoData, NLayers,
verbose = TRUE) {
-
- ##Altitudinal_gradient_functions_______________________________________________________________
- ##unique_gradient_for_precipitation
- GradP_Valery2010 <- function() {
- return(0.00041) ### value from Valery PhD thesis page 126
- }
- ##daily_gradients_for_mean_min_and_max_air_temperature
- GradT_Valery2010 <- function() {
- RESULT <- matrix(c(
- 01, 01, 0.434, 0.366, 0.498,
- 02, 01, 0.434, 0.366, 0.500,
- 03, 01, 0.435, 0.367, 0.501,
- 04, 01, 0.436, 0.367, 0.503,
- 05, 01, 0.437, 0.367, 0.504,
- 06, 01, 0.439, 0.367, 0.506,
- 07, 01, 0.440, 0.367, 0.508,
- 08, 01, 0.441, 0.368, 0.510,
- 09, 01, 0.442, 0.368, 0.512,
- 10, 01, 0.444, 0.368, 0.514,
- 11, 01, 0.445, 0.368, 0.517,
- 12, 01, 0.446, 0.368, 0.519,
- 13, 01, 0.448, 0.369, 0.522,
- 14, 01, 0.450, 0.369, 0.525,
- 15, 01, 0.451, 0.369, 0.527,
- 16, 01, 0.453, 0.370, 0.530,
- 17, 01, 0.455, 0.370, 0.533,
- 18, 01, 0.456, 0.370, 0.537,
- 19, 01, 0.458, 0.371, 0.540,
- 20, 01, 0.460, 0.371, 0.543,
- 21, 01, 0.462, 0.371, 0.547,
- 22, 01, 0.464, 0.372, 0.550,
- 23, 01, 0.466, 0.372, 0.554,
- 24, 01, 0.468, 0.373, 0.558,
- 25, 01, 0.470, 0.373, 0.561,
- 26, 01, 0.472, 0.374, 0.565,
- 27, 01, 0.474, 0.374, 0.569,
- 28, 01, 0.476, 0.375, 0.573,
- 29, 01, 0.478, 0.375, 0.577,
- 30, 01, 0.480, 0.376, 0.582,
- 31, 01, 0.483, 0.376, 0.586,
- 01, 02, 0.485, 0.377, 0.590,
- 02, 02, 0.487, 0.377, 0.594,
- 03, 02, 0.489, 0.378, 0.599,
- 04, 02, 0.492, 0.379, 0.603,
- 05, 02, 0.494, 0.379, 0.607,
- 06, 02, 0.496, 0.380, 0.612,
- 07, 02, 0.498, 0.381, 0.616,
- 08, 02, 0.501, 0.381, 0.621,
- 09, 02, 0.503, 0.382, 0.625,
- 10, 02, 0.505, 0.383, 0.630,
- 11, 02, 0.508, 0.384, 0.634,
- 12, 02, 0.510, 0.384, 0.639,
- 13, 02, 0.512, 0.385, 0.643,
- 14, 02, 0.515, 0.386, 0.648,
- 15, 02, 0.517, 0.387, 0.652,
- 16, 02, 0.519, 0.387, 0.657,
- 17, 02, 0.522, 0.388, 0.661,
- 18, 02, 0.524, 0.389, 0.666,
- 19, 02, 0.526, 0.390, 0.670,
- 20, 02, 0.528, 0.391, 0.674,
- 21, 02, 0.530, 0.392, 0.679,
- 22, 02, 0.533, 0.393, 0.683,
- 23, 02, 0.535, 0.393, 0.687,
- 24, 02, 0.537, 0.394, 0.691,
- 25, 02, 0.539, 0.395, 0.695,
- 26, 02, 0.541, 0.396, 0.699,
- 27, 02, 0.543, 0.397, 0.703,
- 28, 02, 0.545, 0.398, 0.707,
- 29, 02, 0.546, 0.399, 0.709,
- 01, 03, 0.547, 0.399, 0.711,
- 02, 03, 0.549, 0.400, 0.715,
- 03, 03, 0.551, 0.401, 0.718,
- 04, 03, 0.553, 0.402, 0.722,
- 05, 03, 0.555, 0.403, 0.726,
- 06, 03, 0.557, 0.404, 0.729,
- 07, 03, 0.559, 0.405, 0.732,
- 08, 03, 0.560, 0.406, 0.736,
- 09, 03, 0.562, 0.406, 0.739,
- 10, 03, 0.564, 0.407, 0.742,
- 11, 03, 0.566, 0.408, 0.745,
- 12, 03, 0.567, 0.409, 0.748,
- 13, 03, 0.569, 0.410, 0.750,
- 14, 03, 0.570, 0.411, 0.753,
- 15, 03, 0.572, 0.412, 0.756,
- 16, 03, 0.573, 0.413, 0.758,
- 17, 03, 0.575, 0.414, 0.761,
- 18, 03, 0.576, 0.415, 0.763,
- 19, 03, 0.577, 0.416, 0.765,
- 20, 03, 0.579, 0.417, 0.767,
- 21, 03, 0.580, 0.417, 0.769,
- 22, 03, 0.581, 0.418, 0.771,
- 23, 03, 0.582, 0.419, 0.773,
- 24, 03, 0.583, 0.420, 0.774,
- 25, 03, 0.584, 0.421, 0.776,
- 26, 03, 0.585, 0.422, 0.777,
- 27, 03, 0.586, 0.422, 0.779,
- 28, 03, 0.587, 0.423, 0.780,
- 29, 03, 0.588, 0.424, 0.781,
- 30, 03, 0.589, 0.425, 0.782,
- 31, 03, 0.590, 0.425, 0.783,
- 01, 04, 0.591, 0.426, 0.784,
- 02, 04, 0.591, 0.427, 0.785,
- 03, 04, 0.592, 0.427, 0.785,
- 04, 04, 0.593, 0.428, 0.786,
- 05, 04, 0.593, 0.429, 0.787,
- 06, 04, 0.594, 0.429, 0.787,
- 07, 04, 0.595, 0.430, 0.787,
- 08, 04, 0.595, 0.431, 0.788,
- 09, 04, 0.596, 0.431, 0.788,
- 10, 04, 0.596, 0.432, 0.788,
- 11, 04, 0.597, 0.432, 0.788,
- 12, 04, 0.597, 0.433, 0.788,
- 13, 04, 0.597, 0.433, 0.788,
- 14, 04, 0.598, 0.434, 0.788,
- 15, 04, 0.598, 0.434, 0.788,
- 16, 04, 0.598, 0.435, 0.787,
- 17, 04, 0.599, 0.435, 0.787,
- 18, 04, 0.599, 0.436, 0.787,
- 19, 04, 0.599, 0.436, 0.786,
- 20, 04, 0.599, 0.436, 0.786,
- 21, 04, 0.600, 0.437, 0.785,
- 22, 04, 0.600, 0.437, 0.785,
- 23, 04, 0.600, 0.437, 0.784,
- 24, 04, 0.600, 0.438, 0.784,
- 25, 04, 0.600, 0.438, 0.783,
- 26, 04, 0.601, 0.438, 0.783,
- 27, 04, 0.601, 0.438, 0.782,
- 28, 04, 0.601, 0.439, 0.781,
- 29, 04, 0.601, 0.439, 0.781,
- 30, 04, 0.601, 0.439, 0.780,
- 01, 05, 0.601, 0.439, 0.779,
- 02, 05, 0.601, 0.439, 0.778,
- 03, 05, 0.601, 0.439, 0.778,
- 04, 05, 0.601, 0.440, 0.777,
- 05, 05, 0.601, 0.440, 0.776,
- 06, 05, 0.601, 0.440, 0.775,
- 07, 05, 0.601, 0.440, 0.775,
- 08, 05, 0.601, 0.440, 0.774,
- 09, 05, 0.601, 0.440, 0.773,
- 10, 05, 0.602, 0.440, 0.772,
- 11, 05, 0.602, 0.440, 0.772,
- 12, 05, 0.602, 0.440, 0.771,
- 13, 05, 0.602, 0.440, 0.770,
- 14, 05, 0.602, 0.440, 0.770,
- 15, 05, 0.602, 0.440, 0.769,
- 16, 05, 0.602, 0.440, 0.768,
- 17, 05, 0.602, 0.440, 0.768,
- 18, 05, 0.602, 0.440, 0.767,
- 19, 05, 0.602, 0.440, 0.767,
- 20, 05, 0.602, 0.440, 0.766,
- 21, 05, 0.602, 0.440, 0.766,
- 22, 05, 0.602, 0.440, 0.765,
- 23, 05, 0.602, 0.440, 0.765,
- 24, 05, 0.602, 0.440, 0.764,
- 25, 05, 0.602, 0.440, 0.764,
- 26, 05, 0.602, 0.440, 0.764,
- 27, 05, 0.602, 0.439, 0.763,
- 28, 05, 0.602, 0.439, 0.763,
- 29, 05, 0.602, 0.439, 0.763,
- 30, 05, 0.602, 0.439, 0.762,
- 31, 05, 0.602, 0.439, 0.762,
- 01, 06, 0.602, 0.439, 0.762,
- 02, 06, 0.602, 0.439, 0.762,
- 03, 06, 0.602, 0.439, 0.762,
- 04, 06, 0.602, 0.439, 0.762,
- 05, 06, 0.602, 0.439, 0.762,
- 06, 06, 0.602, 0.438, 0.761,
- 07, 06, 0.602, 0.438, 0.761,
- 08, 06, 0.602, 0.438, 0.761,
- 09, 06, 0.602, 0.438, 0.761,
- 10, 06, 0.602, 0.438, 0.761,
- 11, 06, 0.602, 0.438, 0.762,
- 12, 06, 0.602, 0.438, 0.762,
- 13, 06, 0.602, 0.438, 0.762,
- 14, 06, 0.602, 0.438, 0.762,
- 15, 06, 0.602, 0.437, 0.762,
- 16, 06, 0.602, 0.437, 0.762,
- 17, 06, 0.602, 0.437, 0.762,
- 18, 06, 0.602, 0.437, 0.762,
- 19, 06, 0.602, 0.437, 0.763,
- 20, 06, 0.602, 0.437, 0.763,
- 21, 06, 0.602, 0.437, 0.763,
- 22, 06, 0.602, 0.436, 0.763,
- 23, 06, 0.602, 0.436, 0.763,
- 24, 06, 0.602, 0.436, 0.764,
- 25, 06, 0.602, 0.436, 0.764,
- 26, 06, 0.601, 0.436, 0.764,
- 27, 06, 0.601, 0.436, 0.764,
- 28, 06, 0.601, 0.436, 0.764,
- 29, 06, 0.601, 0.435, 0.765,
- 30, 06, 0.601, 0.435, 0.765,
- 01, 07, 0.601, 0.435, 0.765,
- 02, 07, 0.600, 0.435, 0.765,
- 03, 07, 0.600, 0.435, 0.765,
- 04, 07, 0.600, 0.434, 0.766,
- 05, 07, 0.600, 0.434, 0.766,
- 06, 07, 0.599, 0.434, 0.766,
- 07, 07, 0.599, 0.434, 0.766,
- 08, 07, 0.599, 0.434, 0.766,
- 09, 07, 0.598, 0.433, 0.766,
- 10, 07, 0.598, 0.433, 0.766,
- 11, 07, 0.598, 0.433, 0.766,
- 12, 07, 0.597, 0.433, 0.766,
- 13, 07, 0.597, 0.432, 0.767,
- 14, 07, 0.597, 0.432, 0.767,
- 15, 07, 0.596, 0.432, 0.767,
- 16, 07, 0.596, 0.432, 0.766,
- 17, 07, 0.595, 0.431, 0.766,
- 18, 07, 0.595, 0.431, 0.766,
- 19, 07, 0.594, 0.431, 0.766,
- 20, 07, 0.594, 0.430, 0.766,
- 21, 07, 0.593, 0.430, 0.766,
- 22, 07, 0.593, 0.430, 0.766,
- 23, 07, 0.592, 0.429, 0.765,
- 24, 07, 0.592, 0.429, 0.765,
- 25, 07, 0.591, 0.428, 0.765,
- 26, 07, 0.590, 0.428, 0.765,
- 27, 07, 0.590, 0.428, 0.764,
- 28, 07, 0.589, 0.427, 0.764,
- 29, 07, 0.588, 0.427, 0.764,
- 30, 07, 0.588, 0.426, 0.763,
- 31, 07, 0.587, 0.426, 0.763,
- 01, 08, 0.586, 0.425, 0.762,
- 02, 08, 0.586, 0.425, 0.762,
- 03, 08, 0.585, 0.424, 0.761,
- 04, 08, 0.584, 0.424, 0.761,
- 05, 08, 0.583, 0.423, 0.760,
- 06, 08, 0.583, 0.423, 0.760,
- 07, 08, 0.582, 0.422, 0.759,
- 08, 08, 0.581, 0.421, 0.758,
- 09, 08, 0.580, 0.421, 0.758,
- 10, 08, 0.579, 0.420, 0.757,
- 11, 08, 0.578, 0.420, 0.756,
- 12, 08, 0.578, 0.419, 0.755,
- 13, 08, 0.577, 0.418, 0.754,
- 14, 08, 0.576, 0.418, 0.754,
- 15, 08, 0.575, 0.417, 0.753,
- 16, 08, 0.574, 0.416, 0.752,
- 17, 08, 0.573, 0.415, 0.751,
- 18, 08, 0.572, 0.415, 0.750,
- 19, 08, 0.571, 0.414, 0.749,
- 20, 08, 0.570, 0.413, 0.748,
- 21, 08, 0.569, 0.413, 0.747,
- 22, 08, 0.569, 0.412, 0.746,
- 23, 08, 0.568, 0.411, 0.745,
- 24, 08, 0.567, 0.410, 0.744,
- 25, 08, 0.566, 0.409, 0.743,
- 26, 08, 0.565, 0.409, 0.742,
- 27, 08, 0.564, 0.408, 0.741,
- 28, 08, 0.563, 0.407, 0.740,
- 29, 08, 0.562, 0.406, 0.738,
- 30, 08, 0.561, 0.405, 0.737,
- 31, 08, 0.560, 0.405, 0.736,
- 01, 09, 0.558, 0.404, 0.735,
- 02, 09, 0.557, 0.403, 0.734,
- 03, 09, 0.556, 0.402, 0.732,
- 04, 09, 0.555, 0.401, 0.731,
- 05, 09, 0.554, 0.401, 0.730,
- 06, 09, 0.553, 0.400, 0.728,
- 07, 09, 0.552, 0.399, 0.727,
- 08, 09, 0.551, 0.398, 0.725,
- 09, 09, 0.550, 0.397, 0.724,
- 10, 09, 0.549, 0.396, 0.723,
- 11, 09, 0.548, 0.396, 0.721,
- 12, 09, 0.546, 0.395, 0.720,
- 13, 09, 0.545, 0.394, 0.718,
- 14, 09, 0.544, 0.393, 0.717,
- 15, 09, 0.543, 0.392, 0.715,
- 16, 09, 0.542, 0.391, 0.713,
- 17, 09, 0.541, 0.391, 0.712,
- 18, 09, 0.540, 0.390, 0.710,
- 19, 09, 0.538, 0.389, 0.709,
- 20, 09, 0.537, 0.388, 0.707,
- 21, 09, 0.536, 0.388, 0.705,
- 22, 09, 0.535, 0.387, 0.703,
- 23, 09, 0.533, 0.386, 0.702,
- 24, 09, 0.532, 0.385, 0.700,
- 25, 09, 0.531, 0.385, 0.698,
- 26, 09, 0.530, 0.384, 0.696,
- 27, 09, 0.528, 0.383, 0.694,
- 28, 09, 0.527, 0.383, 0.692,
- 29, 09, 0.526, 0.382, 0.690,
- 30, 09, 0.525, 0.381, 0.688,
- 01, 10, 0.523, 0.381, 0.686,
- 02, 10, 0.522, 0.380, 0.684,
- 03, 10, 0.521, 0.379, 0.682,
- 04, 10, 0.519, 0.379, 0.680,
- 05, 10, 0.518, 0.378, 0.678,
- 06, 10, 0.517, 0.377, 0.676,
- 07, 10, 0.515, 0.377, 0.674,
- 08, 10, 0.514, 0.376, 0.671,
- 09, 10, 0.512, 0.376, 0.669,
- 10, 10, 0.511, 0.375, 0.667,
- 11, 10, 0.510, 0.375, 0.664,
- 12, 10, 0.508, 0.374, 0.662,
- 13, 10, 0.507, 0.374, 0.659,
- 14, 10, 0.505, 0.373, 0.657,
- 15, 10, 0.504, 0.373, 0.654,
- 16, 10, 0.502, 0.372, 0.652,
- 17, 10, 0.501, 0.372, 0.649,
- 18, 10, 0.499, 0.372, 0.647,
- 19, 10, 0.498, 0.371, 0.644,
- 20, 10, 0.496, 0.371, 0.641,
- 21, 10, 0.495, 0.371, 0.639,
- 22, 10, 0.493, 0.370, 0.636,
- 23, 10, 0.492, 0.370, 0.633,
- 24, 10, 0.490, 0.370, 0.630,
- 25, 10, 0.489, 0.369, 0.628,
- 26, 10, 0.487, 0.369, 0.625,
- 27, 10, 0.485, 0.369, 0.622,
- 28, 10, 0.484, 0.368, 0.619,
- 29, 10, 0.482, 0.368, 0.616,
- 30, 10, 0.481, 0.368, 0.613,
- 31, 10, 0.479, 0.368, 0.610,
- 01, 11, 0.478, 0.368, 0.607,
- 02, 11, 0.476, 0.367, 0.604,
- 03, 11, 0.475, 0.367, 0.601,
- 04, 11, 0.473, 0.367, 0.598,
- 05, 11, 0.471, 0.367, 0.595,
- 06, 11, 0.470, 0.367, 0.592,
- 07, 11, 0.468, 0.367, 0.589,
- 08, 11, 0.467, 0.366, 0.586,
- 09, 11, 0.465, 0.366, 0.583,
- 10, 11, 0.464, 0.366, 0.580,
- 11, 11, 0.462, 0.366, 0.577,
- 12, 11, 0.461, 0.366, 0.574,
- 13, 11, 0.459, 0.366, 0.571,
- 14, 11, 0.458, 0.366, 0.568,
- 15, 11, 0.456, 0.366, 0.565,
- 16, 11, 0.455, 0.366, 0.562,
- 17, 11, 0.454, 0.366, 0.559,
- 18, 11, 0.452, 0.365, 0.556,
- 19, 11, 0.451, 0.365, 0.553,
- 20, 11, 0.450, 0.365, 0.550,
- 21, 11, 0.448, 0.365, 0.547,
- 22, 11, 0.447, 0.365, 0.544,
- 23, 11, 0.446, 0.365, 0.542,
- 24, 11, 0.445, 0.365, 0.539,
- 25, 11, 0.443, 0.365, 0.536,
- 26, 11, 0.442, 0.365, 0.533,
- 27, 11, 0.441, 0.365, 0.531,
- 28, 11, 0.440, 0.365, 0.528,
- 29, 11, 0.439, 0.365, 0.526,
- 30, 11, 0.438, 0.365, 0.523,
- 01, 12, 0.437, 0.365, 0.521,
- 02, 12, 0.436, 0.365, 0.519,
- 03, 12, 0.435, 0.365, 0.517,
- 04, 12, 0.434, 0.365, 0.515,
- 05, 12, 0.434, 0.365, 0.513,
- 06, 12, 0.433, 0.365, 0.511,
- 07, 12, 0.432, 0.365, 0.509,
- 08, 12, 0.431, 0.365, 0.507,
- 09, 12, 0.431, 0.365, 0.505,
- 10, 12, 0.430, 0.365, 0.504,
- 11, 12, 0.430, 0.365, 0.502,
- 12, 12, 0.429, 0.365, 0.501,
- 13, 12, 0.429, 0.365, 0.500,
- 14, 12, 0.429, 0.365, 0.498,
- 15, 12, 0.428, 0.365, 0.497,
- 16, 12, 0.428, 0.365, 0.496,
- 17, 12, 0.428, 0.365, 0.496,
- 18, 12, 0.428, 0.365, 0.495,
- 19, 12, 0.428, 0.365, 0.494,
- 20, 12, 0.428, 0.365, 0.494,
- 21, 12, 0.428, 0.365, 0.494,
- 22, 12, 0.428, 0.365, 0.493,
- 23, 12, 0.429, 0.365, 0.493,
- 24, 12, 0.429, 0.366, 0.493,
- 25, 12, 0.429, 0.366, 0.493,
- 26, 12, 0.430, 0.366, 0.494,
- 27, 12, 0.430, 0.366, 0.494,
- 28, 12, 0.431, 0.366, 0.495,
- 29, 12, 0.431, 0.366, 0.495,
- 30, 12, 0.432, 0.366, 0.496,
- 31, 12, 0.433, 0.366, 0.497), ncol = 5, byrow = TRUE)
-
- dimnames(RESULT) <- list(1:366, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))
-
- return(RESULT)
-
- }
-
-
-
+
+ ##Altitudinal_gradient_functions_______________________________________________________________
+ ##unique_gradient_for_precipitation
+ GradP_Valery2010 <- 0.00041 ### value from Valery PhD thesis page 126
+
+
+
##Format_______________________________________________________________________________________
- HypsoData <- as.double(HypsoData)
- ZInputs <- as.double(ZInputs)
-
-
-
+ HypsoData <- as.double(HypsoData)
+ ZInputs <- as.double(ZInputs)
+
+
+
##ElevationLayers_Creation_____________________________________________________________________
- ZLayers <- as.double(rep(NA, NLayers))
-
- if (!identical(HypsoData, as.double(rep(NA, 101)))) {
- nmoy <- 100 %/% NLayers
- nreste <- 100 %% NLayers
- ncont <- 0
-
- for (iLayer in 1:NLayers) {
- if (nreste > 0) {
- nn <- nmoy + 1
- nreste <- nreste - 1
- } else {
- nn <- nmoy
- }
- if (nn == 1) {
- ZLayers[iLayer] <- HypsoData[ncont + 1]
- }
- if (nn == 2) {
- ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
- }
- if (nn > 2) {
- ZLayers[iLayer] <- HypsoData[ncont + nn / 2]
- }
- ncont <- ncont + nn
+ ZLayers <- as.double(rep(NA, NLayers))
+
+ if (!identical(HypsoData, as.double(rep(NA, 101)))) {
+ nmoy <- 100 %/% NLayers
+ nreste <- 100 %% NLayers
+ ncont <- 0
+
+ for (iLayer in 1:NLayers) {
+ if (nreste > 0) {
+ nn <- nmoy + 1
+ nreste <- nreste - 1
+ } else {
+ nn <- nmoy
+ }
+ if (nn == 1) {
+ ZLayers[iLayer] <- HypsoData[ncont + 1]
+ }
+ if (nn == 2) {
+ ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
+ }
+ if (nn > 2) {
+ ZLayers[iLayer] <- HypsoData[ncont + nn / 2]
}
+ ncont <- ncont + nn
}
-
-
+ }
+
+
##Precipitation_extrapolation__________________________________________________________________
##Initialisation
- if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
- LayerPrecip <- list(as.double(Precip))
- } else {
- ##Elevation_gradients_for_daily_mean_precipitation
- GradP <- GradP_Valery2010() ### single value
- TabGradP <- rep(GradP, length(Precip))
- ##Extrapolation
- ##Thresold_of_inputs_median_elevation
- Zthreshold <- 4000
- LayerPrecip_mat <- sapply(1:NLayers, function(iLayer) {
- ##If_layer_elevation_smaller_than_Zthreshold
- if (ZLayers[iLayer] <= Zthreshold) {
- prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
- ##If_layer_elevation_greater_than_Zthreshold
+ if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
+ LayerPrecip <- list(as.double(Precip))
+ } else {
+ ##Elevation_gradients_for_daily_mean_precipitation
+ GradP <- GradP_Valery2010 ### single value
+ TabGradP <- rep(GradP, length(Precip))
+ ##Extrapolation
+ ##Thresold_of_inputs_median_elevation
+ Zthreshold <- 4000
+ LayerPrecip_mat <- sapply(1:NLayers, function(iLayer) {
+ ##If_layer_elevation_smaller_than_Zthreshold
+ if (ZLayers[iLayer] <= Zthreshold) {
+ prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
+ ##If_layer_elevation_greater_than_Zthreshold
+ } else {
+ ##If_inputs_median_elevation_smaller_than_Zthreshold
+ if (ZInputs <= Zthreshold) {
+ prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
+ ##If_inputs_median_elevation_greater_then_Zthreshold
} else {
- ##If_inputs_median_elevation_smaller_than_Zthreshold
- if (ZInputs <= Zthreshold) {
- prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
- ##If_inputs_median_elevation_greater_then_Zthreshold
- } else {
- prcp <- as.double(Precip)
- }
+ prcp <- as.double(Precip)
}
- return(prcp)
- })
- if (PrecipScale) {
- LayerPrecip_mat <- LayerPrecip_mat / rowMeans(LayerPrecip_mat) * Precip
- LayerPrecip_mat[is.nan(LayerPrecip_mat)] <- 0
}
- LayerPrecip <- as.list(as.data.frame(LayerPrecip_mat))
+ return(prcp)
+ })
+ if (PrecipScale) {
+ LayerPrecip_mat <- LayerPrecip_mat / rowMeans(LayerPrecip_mat) * Precip
+ LayerPrecip_mat[is.nan(LayerPrecip_mat)] <- 0
}
-
-
-
+ LayerPrecip <- as.list(as.data.frame(LayerPrecip_mat))
+ }
+
+
+
##Temperature_extrapolation____________________________________________________________________
##Initialisation
- LayerTempMean <- list()
- LayerTempMin <- list()
- LayerTempMax <- list()
-
- if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
- LayerTempMean[[1]] <- as.double(TempMean)
-
+ LayerTempMean <- list()
+ LayerTempMin <- list()
+ LayerTempMax <- list()
+
+ if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
+ LayerTempMean[[1]] <- as.double(TempMean)
+
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ LayerTempMin[[1]] <- as.double(TempMin)
+ LayerTempMax[[1]] <- as.double(TempMax)
+ }
+ } else {
+ ##Elevation_gradients_for_daily_mean_min_and_max_temperature
+ GradT <- .GradT_Valery2010
+ iday <- match(format(DatesR, format = "%d%m"),
+ sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
+ TabGradT <- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
+ ##Extrapolation
+ ##On_each_elevation_layer...
+ for (iLayer in 1:NLayers) {
+ LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
if (!is.null(TempMin) & !is.null(TempMax)) {
- LayerTempMin[[1]] <- as.double(TempMin)
- LayerTempMax[[1]] <- as.double(TempMax)
- }
- } else {
- ##Elevation_gradients_for_daily_mean_min_and_max_temperature
- GradT <- as.data.frame(GradT_Valery2010())
- iday <- match(format(DatesR, format = "%d%m"),
- sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
- TabGradT <-
- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
- ##Extrapolation
- ##On_each_elevation_layer...
- for (iLayer in 1:NLayers) {
- LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
- if (!is.null(TempMin) & !is.null(TempMax)) {
- LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
- LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
- }
+ LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
+ LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
}
}
-
-
-
+ }
+
+
+
##Solid_Fraction_for_each_elevation_layer______________________________________________________
- LayerFracSolidPrecip <- list()
-
- ##Thresold_of_inputs_median_elevation
- Zthreshold <- 1500
-
- ##Option
- Option <- "USACE"
- if (!is.na(ZInputs)) {
- if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) {
- Option <- "Hydrotel"
- }
+ LayerFracSolidPrecip <- list()
+
+ ##Thresold_of_inputs_median_elevation
+ Zthreshold <- 1500
+
+ ##Option
+ Option <- "USACE"
+ if (!is.na(ZInputs)) {
+ if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) {
+ Option <- "Hydrotel"
}
-
- ##On_each_elevation_layer...
- for (iLayer in 1:NLayers) {
+ }
- ##Turcotte_formula_from_Hydrotel
- if (Option == "Hydrotel") {
- TempMin <- LayerTempMin[[iLayer]]
- TempMax <- LayerTempMax[[iLayer]]
- SolidFraction <- 1 - TempMax / (TempMax - TempMin)
- SolidFraction[TempMin >= 0] <- 0
- SolidFraction[TempMax <= 0] <- 1
- }
- ##USACE_formula
- if (Option == "USACE") {
- USACE_Tmin <- -1.0
- USACE_Tmax <- 3.0
- TempMean <- LayerTempMean[[iLayer]]
- SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
- SolidFraction[TempMean > USACE_Tmax] <- 0
- SolidFraction[TempMean < USACE_Tmin] <- 1
- }
- LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
+ ##On_each_elevation_layer...
+ for (iLayer in 1:NLayers) {
+
+ ##Turcotte_formula_from_Hydrotel
+ if (Option == "Hydrotel") {
+ TempMin <- LayerTempMin[[iLayer]]
+ TempMax <- LayerTempMax[[iLayer]]
+ SolidFraction <- 1 - TempMax / (TempMax - TempMin)
+ SolidFraction[TempMin >= 0] <- 0
+ SolidFraction[TempMax <= 0] <- 1
}
- namesLayer <- sprintf("L%i", seq_along(LayerPrecip))
- names(LayerPrecip) <- namesLayer
- names(LayerTempMean) <- namesLayer
- if (!is.null(TempMin) & !is.null(TempMax)) {
- names(LayerTempMin) <- namesLayer
- names(LayerTempMax) <- namesLayer
- }
- names(LayerFracSolidPrecip) <- namesLayer
-
-
-
+ ##USACE_formula
+ if (Option == "USACE") {
+ USACE_Tmin <- -1.0
+ USACE_Tmax <- 3.0
+ TempMean <- LayerTempMean[[iLayer]]
+ SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
+ SolidFraction[TempMean > USACE_Tmax] <- 0
+ SolidFraction[TempMean < USACE_Tmin] <- 1
+ }
+ LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
+ }
+ namesLayer <- sprintf("L%i", seq_along(LayerPrecip))
+ names(LayerPrecip) <- namesLayer
+ names(LayerTempMean) <- namesLayer
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ names(LayerTempMin) <- namesLayer
+ names(LayerTempMax) <- namesLayer
+ }
+ names(LayerFracSolidPrecip) <- namesLayer
+
+
+
##END__________________________________________________________________________________________
- return(list(LayerPrecip = LayerPrecip,
- LayerTempMean = LayerTempMean,
- LayerTempMin = LayerTempMin,
- LayerTempMax = LayerTempMax,
- LayerFracSolidPrecip = LayerFracSolidPrecip,
- ZLayers = ZLayers))
-
-
+ return(list(LayerPrecip = LayerPrecip,
+ LayerTempMean = LayerTempMean,
+ LayerTempMin = LayerTempMin,
+ LayerTempMax = LayerTempMax,
+ LayerFracSolidPrecip = LayerFracSolidPrecip,
+ ZLayers = ZLayers))
+
+
}
diff --git a/R/Imax.R b/R/Imax.R
index 3c6058f1b1c4945ee734b2a50a185c5ae3d9f8c5..b783ce65939c9d93da948c9065676649c1467a66 100644
--- a/R/Imax.R
+++ b/R/Imax.R
@@ -1,17 +1,20 @@
-Imax <- function(InputsModel,
- IndPeriod_Run,
+Imax <- function(InputsModel,
+ IndPeriod_Run,
TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {
-
- ##_____Arguments_check_____________________________________________________________________
+
+
+ ## ---------- check arguments
+
+ ## InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
- }
+ }
if (!inherits(InputsModel, "hourly")) {
stop("'InputsModel' must be of class 'hourly'")
- }
-
- ##check_IndPeriod_Run
+ }
+
+ ## IndPeriod_Run
if (!is.vector(IndPeriod_Run)) {
stop("'IndPeriod_Run' must be a vector of numeric values")
}
@@ -21,24 +24,28 @@ Imax <- function(InputsModel,
if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
stop("'IndPeriod_Run' must be a continuous sequence of integers")
}
-
- ##TestedValues
+
+ ## TestedValues
if (!(is.numeric(TestedValues))) {
stop("'TestedValues' must be 'numeric'")
}
-
- ##aggregate data at the daily time step
- TabSeries <- data.frame(DatesR = InputsModel$DatesR[IndPeriod_Run],
- Precip = InputsModel$Precip[IndPeriod_Run],
- PotEvap = InputsModel$PotEvap[IndPeriod_Run])
- daily_data <- SeriesAggreg(TabSeries, Format = "%Y%m%d",
- ConvertFun = c("sum", "sum"))
-
- ##calculate total interception of daily GR models on the period
+
+ ## ---------- hourly inputs aggregation
+
+ ## aggregate data at the daily time step
+ daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")
+
+
+ ## ---------- calculate interception
+
+ ## calculate total interception of daily GR models on the period
cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
-
- ##calculate total interception of the GR5H interception store on the period
+ if (anyNA(cum_daily)) {
+ stop("'IndPeriod_Run' must be set to select 24 hours by day")
+ }
+
+ ## calculate total interception of the GR5H interception store on the period
## and compute difference with daily values
differences <- array(NA, c(length(TestedValues)))
for (Imax in TestedValues) {
@@ -52,8 +59,8 @@ Imax <- function(InputsModel,
}
differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
}
-
- ##return the Imax value that minimises the difference
+
+ ## return the Imax value that minimises the difference
return(TestedValues[which.min(differences)])
-
+
}
diff --git a/R/PE_Oudin.R b/R/PE_Oudin.R
index 2ce49b334736ff8ab69e7ac6823be62d7b4f3684..173c9e4fd812ba81d8fca76388b7eb959e1083fe 100644
--- a/R/PE_Oudin.R
+++ b/R/PE_Oudin.R
@@ -65,7 +65,7 @@ PE_Oudin <- function(JD, Temp,
LInputs = as.integer(length(Temp))
if (length(FI) == 1) {
- FI <- rep(FI, LInputs)
+ FI <- rep(FI, LInputs)
}
RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
@@ -96,7 +96,7 @@ PE_Oudin <- function(JD, Temp,
COSOM <- -1
}
if (COSOM > 1) {
- COSOM <- 1
+ COSOM <- 1
}
COSOM2 <- COSOM * COSOM
diff --git a/R/PEdaily_Oudin.R b/R/PEdaily_Oudin.R
index 3eee9283489c7d66f89f676f81f255c5b640cf4c..33472c886280417bb7ba0bf4619bfea29376f239 100644
--- a/R/PEdaily_Oudin.R
+++ b/R/PEdaily_Oudin.R
@@ -70,7 +70,7 @@ PEdaily_Oudin <- function(JD,
COSOM <- -1
}
if (COSOM > 1) {
- COSOM <- 1
+ COSOM <- 1
}
COSOM2 <- COSOM * COSOM
@@ -94,11 +94,11 @@ PEdaily_Oudin <- function(JD,
if (is.na(Temp[k])) {
PE_Oudin_D[k] <- NA
} else {
- if (Temp[k] >= -5.0) {
- PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
- } else {
- PE_Oudin_D[k] <- 0
- }
+ if (Temp[k] >= -5.0) {
+ PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
+ } else {
+ PE_Oudin_D[k] <- 0
+ }
}
}
diff --git a/R/RunModel.R b/R/RunModel.R
index a123e366705dd0de3b2cba80d5a3c70eef49e85e..1a6104d15fcceaa8195dc3944923d7110ace56dc 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -3,7 +3,7 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD) {
FUN_MOD <- match.fun(FUN_MOD)
if (inherits(InputsModel, "SD")) {
- # LAG Model take one parameter at the beginning of the vector
+ # Lag model take one parameter at the beginning of the vector
iFirstParamRunOffModel <- 2
} else {
# All parameters
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index 385c5f1a798de5ea4ac460fd3fb0d516c7c6430b..52e068c47e536527bbfbcbfef3e37251dc6befdf 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -59,8 +59,8 @@ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[2] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[2] ### routing store level (mm)
}
## Call GR model Fortan
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index 47c26ff968cce41eb982f897604c7d4b843d0162..749ae6a39ec4a34dcab587cbd552ff65b2ac481e 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -64,8 +64,8 @@ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
## Call GR model Fortan
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index fc647035051dc0f05fb741f0f43ae1f1945ded08..0be613be0d5fca604e72c9696e1d72bde93716ff 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -63,8 +63,8 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
## Call GR model Fortan
diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R
index 590a103291d3ea8e99f26681cd3eb1774ef05edb..3f985cb78b4956225029561409d4bff0ecdc80f0 100644
--- a/R/RunModel_GR5H.R
+++ b/R/RunModel_GR5H.R
@@ -1,6 +1,6 @@
RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
-
-
+
+
## Initialization of variables
NParam <- 5
FortranOutputs <- .FortranOutputs(GR = "GR5H")$GR
@@ -10,24 +10,24 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
} else {
Imax <- -99
}
-
-
+
+
## Arguments check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
- }
+ }
if (!inherits(InputsModel, "hourly" )) {
stop("'InputsModel' must be of class 'hourly' ")
- }
+ }
if (!inherits(InputsModel, "GR" )) {
stop("'InputsModel' must be of class 'GR' ")
- }
+ }
if (!inherits(RunOptions, "RunOptions" )) {
stop("'RunOptions' must be of class 'RunOptions' ")
- }
+ }
if (!inherits(RunOptions, "GR" )) {
stop("'RunOptions' must be of class 'GR' ")
- }
+ }
if (!is.vector(Param) | !is.numeric(Param)) {
stop("'Param' must be a numeric vector")
}
@@ -35,7 +35,7 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
stop(paste("'Param' must be a vector of length ", NParam, " and contain no NA", sep = ""))
}
Param <- as.double(Param)
-
+
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
@@ -49,35 +49,36 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
- }
-
+ }
+
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
- if ("all" %in% RunOptions$Outputs_Sim) { IndOutputs <- as.integer(1:length(FortranOutputs))
+ if ("all" %in% RunOptions$Outputs_Sim) {
+ IndOutputs <- as.integer(1:length(FortranOutputs))
} else {
IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
}
-
+
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-
+
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
if (IsIntStore) {
RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
}
}
-
+
## Call GR model Fortan
- RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
+ RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
@@ -97,14 +98,14 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
- RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
- ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
- IntStore = RESULTS$StateEnd[4L],
- UH1 = NULL, UH2 = RESULTS$StateEnd[(1:(40*24))+(7+20*24)],
- GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+ RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
+ ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+ IntStore = RESULTS$StateEnd[4L],
+ UH1 = NULL, UH2 = RESULTS$StateEnd[(1:(40*24))+(7+20*24)],
+ GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
-
+
## Output data preparation
## OutputsModel only
if (!ExportDatesR & !ExportStateEnd) {
@@ -113,30 +114,30 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
}
## DatesR and OutputsModel only
if (ExportDatesR & !ExportStateEnd) {
- OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
+ OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
}
## OutputsModel and StateEnd only
if (!ExportDatesR & ExportStateEnd) {
- OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
+ OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
list(RESULTS$StateEnd))
names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
}
## DatesR and OutputsModel and StateEnd
if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
- OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
+ OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
+ lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
list(RESULTS$StateEnd))
names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
}
-
+
## End
- rm(RESULTS)
+ rm(RESULTS)
class(OutputsModel) <- c("OutputsModel", "hourly", "GR")
if (IsIntStore) {
class(OutputsModel) <- c(class(OutputsModel), "interception")
}
return(OutputsModel)
-
+
}
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 33210729deffa802b0e51732e0659e27ee7802be..ee8f7756670dcafc0b3450f82cb29c5c536c07c0 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -64,8 +64,8 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
## Call GR model Fortan
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index c971056fbd21f4fd423a0b98abe0065d62c5090b..52df6553e18f143f379b4f29d764f4e929716abf 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -68,8 +68,8 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
- RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1] ### production store level (mm)
- RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3] ### routing store level (mm)
+ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
+ RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 8bfd1f568e73a5a259545047926aaa151a52cfc9..87f551f73f94e3ef91b93798b78ae35f2db5c0b6 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -50,8 +50,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
NbUpBasins <- length(InputsModel$LengthHydro)
LengthTs <- length(OutputsModel$QsimDown)
- OutputsModel$Qsim <-
- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+ OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
@@ -78,12 +77,10 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
}
for (upstream_basin in seq_len(NbUpBasins)) {
- Qupstream <-
- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+ Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
# Upstream flow with area needs to be converted to m3 by time step
- Qupstream <-
- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
+ Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
}
OutputsModel$Qsim <- OutputsModel$Qsim +
c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
diff --git a/R/TransfoParam_CemaNeige.R b/R/TransfoParam_CemaNeige.R
index add5b9192756d0b6e43c3fe45925a312d3f397a6..cd7e32af714493dfff8770a200b8c3c2297e7d39 100644
--- a/R/TransfoParam_CemaNeige.R
+++ b/R/TransfoParam_CemaNeige.R
@@ -28,7 +28,7 @@ TransfoParam_CemaNeige <- function(ParamIn, Direction) {
}
if (Direction == "RT") {
ParamOut <- ParamIn
- ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
+ ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
}
diff --git a/R/TransfoParam_CemaNeigeHyst.R b/R/TransfoParam_CemaNeigeHyst.R
index c4ee434b3f4bb13fd6b90ce389198f3f7348e365..6839edad549991d2a971c28d00dbbc62de37c9ea 100644
--- a/R/TransfoParam_CemaNeigeHyst.R
+++ b/R/TransfoParam_CemaNeigeHyst.R
@@ -22,15 +22,15 @@ TransfoParam_CemaNeigeHyst <- function(ParamIn, Direction) {
## transformation
if (Direction == "TR") {
- ParamOut <- ParamIn
+ ParamOut <- ParamIn
ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] * 5) + 50 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] / 19.98) + 0.5 ### Hyst CV
}
if (Direction == "RT") {
- ParamOut <- ParamIn
- ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
+ ParamOut <- ParamIn
+ ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] - 50) / 5 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] - 0.5) * 19.98 ### Hyst CV
diff --git a/R/TransfoParam_GR1A.R b/R/TransfoParam_GR1A.R
index 59b09d7d3ce8c9af29bac28f2e0115e1c7d761b4..b52e75a43bf798ccaa8bed3320b1e867476ed046 100644
--- a/R/TransfoParam_GR1A.R
+++ b/R/TransfoParam_GR1A.R
@@ -25,7 +25,7 @@ TransfoParam_GR1A <- function(ParamIn, Direction) {
ParamOut <- (ParamIn + 10.0) / 8
}
if (Direction == "RT") {
- ParamOut <- ParamIn * 8 - 10.0
+ ParamOut <- ParamIn * 8 - 10.0
}
diff --git a/R/TransfoParam_Lag.R b/R/TransfoParam_Lag.R
index 2294d0ab2442c860d8b6078580cc42583bc37130..0ee13b0462e535df8da47ca60452115c610bc1f7 100644
--- a/R/TransfoParam_Lag.R
+++ b/R/TransfoParam_Lag.R
@@ -16,7 +16,7 @@ TransfoParam_Lag <- function(ParamIn, Direction) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
- stop(sprintf("the LAG model requires %i parameters", NParam))
+ stop(sprintf("the Lag model requires %i parameters", NParam))
}
@@ -25,7 +25,7 @@ TransfoParam_Lag <- function(ParamIn, Direction) {
ParamOut <- 20 * (ParamIn + 10) / 20.0
}
if (Direction == "RT") {
- ParamOut <- ParamIn * 20.0 / 20 - 10
+ ParamOut <- ParamIn * 20.0 / 20 - 10
}
diff --git a/R/UtilsCemaNeige.R b/R/UtilsCemaNeige.R
new file mode 100644
index 0000000000000000000000000000000000000000..e492575c1e534939aa7de38711a09170573f05d7
--- /dev/null
+++ b/R/UtilsCemaNeige.R
@@ -0,0 +1,370 @@
+## daily gradients for mean, min and max air temperature
+.GradT_Valery2010 <- as.data.frame(matrix(c(
+ 01, 01, 0.434, 0.366, 0.498,
+ 02, 01, 0.434, 0.366, 0.500,
+ 03, 01, 0.435, 0.367, 0.501,
+ 04, 01, 0.436, 0.367, 0.503,
+ 05, 01, 0.437, 0.367, 0.504,
+ 06, 01, 0.439, 0.367, 0.506,
+ 07, 01, 0.440, 0.367, 0.508,
+ 08, 01, 0.441, 0.368, 0.510,
+ 09, 01, 0.442, 0.368, 0.512,
+ 10, 01, 0.444, 0.368, 0.514,
+ 11, 01, 0.445, 0.368, 0.517,
+ 12, 01, 0.446, 0.368, 0.519,
+ 13, 01, 0.448, 0.369, 0.522,
+ 14, 01, 0.450, 0.369, 0.525,
+ 15, 01, 0.451, 0.369, 0.527,
+ 16, 01, 0.453, 0.370, 0.530,
+ 17, 01, 0.455, 0.370, 0.533,
+ 18, 01, 0.456, 0.370, 0.537,
+ 19, 01, 0.458, 0.371, 0.540,
+ 20, 01, 0.460, 0.371, 0.543,
+ 21, 01, 0.462, 0.371, 0.547,
+ 22, 01, 0.464, 0.372, 0.550,
+ 23, 01, 0.466, 0.372, 0.554,
+ 24, 01, 0.468, 0.373, 0.558,
+ 25, 01, 0.470, 0.373, 0.561,
+ 26, 01, 0.472, 0.374, 0.565,
+ 27, 01, 0.474, 0.374, 0.569,
+ 28, 01, 0.476, 0.375, 0.573,
+ 29, 01, 0.478, 0.375, 0.577,
+ 30, 01, 0.480, 0.376, 0.582,
+ 31, 01, 0.483, 0.376, 0.586,
+ 01, 02, 0.485, 0.377, 0.590,
+ 02, 02, 0.487, 0.377, 0.594,
+ 03, 02, 0.489, 0.378, 0.599,
+ 04, 02, 0.492, 0.379, 0.603,
+ 05, 02, 0.494, 0.379, 0.607,
+ 06, 02, 0.496, 0.380, 0.612,
+ 07, 02, 0.498, 0.381, 0.616,
+ 08, 02, 0.501, 0.381, 0.621,
+ 09, 02, 0.503, 0.382, 0.625,
+ 10, 02, 0.505, 0.383, 0.630,
+ 11, 02, 0.508, 0.384, 0.634,
+ 12, 02, 0.510, 0.384, 0.639,
+ 13, 02, 0.512, 0.385, 0.643,
+ 14, 02, 0.515, 0.386, 0.648,
+ 15, 02, 0.517, 0.387, 0.652,
+ 16, 02, 0.519, 0.387, 0.657,
+ 17, 02, 0.522, 0.388, 0.661,
+ 18, 02, 0.524, 0.389, 0.666,
+ 19, 02, 0.526, 0.390, 0.670,
+ 20, 02, 0.528, 0.391, 0.674,
+ 21, 02, 0.530, 0.392, 0.679,
+ 22, 02, 0.533, 0.393, 0.683,
+ 23, 02, 0.535, 0.393, 0.687,
+ 24, 02, 0.537, 0.394, 0.691,
+ 25, 02, 0.539, 0.395, 0.695,
+ 26, 02, 0.541, 0.396, 0.699,
+ 27, 02, 0.543, 0.397, 0.703,
+ 28, 02, 0.545, 0.398, 0.707,
+ 29, 02, 0.546, 0.399, 0.709,
+ 01, 03, 0.547, 0.399, 0.711,
+ 02, 03, 0.549, 0.400, 0.715,
+ 03, 03, 0.551, 0.401, 0.718,
+ 04, 03, 0.553, 0.402, 0.722,
+ 05, 03, 0.555, 0.403, 0.726,
+ 06, 03, 0.557, 0.404, 0.729,
+ 07, 03, 0.559, 0.405, 0.732,
+ 08, 03, 0.560, 0.406, 0.736,
+ 09, 03, 0.562, 0.406, 0.739,
+ 10, 03, 0.564, 0.407, 0.742,
+ 11, 03, 0.566, 0.408, 0.745,
+ 12, 03, 0.567, 0.409, 0.748,
+ 13, 03, 0.569, 0.410, 0.750,
+ 14, 03, 0.570, 0.411, 0.753,
+ 15, 03, 0.572, 0.412, 0.756,
+ 16, 03, 0.573, 0.413, 0.758,
+ 17, 03, 0.575, 0.414, 0.761,
+ 18, 03, 0.576, 0.415, 0.763,
+ 19, 03, 0.577, 0.416, 0.765,
+ 20, 03, 0.579, 0.417, 0.767,
+ 21, 03, 0.580, 0.417, 0.769,
+ 22, 03, 0.581, 0.418, 0.771,
+ 23, 03, 0.582, 0.419, 0.773,
+ 24, 03, 0.583, 0.420, 0.774,
+ 25, 03, 0.584, 0.421, 0.776,
+ 26, 03, 0.585, 0.422, 0.777,
+ 27, 03, 0.586, 0.422, 0.779,
+ 28, 03, 0.587, 0.423, 0.780,
+ 29, 03, 0.588, 0.424, 0.781,
+ 30, 03, 0.589, 0.425, 0.782,
+ 31, 03, 0.590, 0.425, 0.783,
+ 01, 04, 0.591, 0.426, 0.784,
+ 02, 04, 0.591, 0.427, 0.785,
+ 03, 04, 0.592, 0.427, 0.785,
+ 04, 04, 0.593, 0.428, 0.786,
+ 05, 04, 0.593, 0.429, 0.787,
+ 06, 04, 0.594, 0.429, 0.787,
+ 07, 04, 0.595, 0.430, 0.787,
+ 08, 04, 0.595, 0.431, 0.788,
+ 09, 04, 0.596, 0.431, 0.788,
+ 10, 04, 0.596, 0.432, 0.788,
+ 11, 04, 0.597, 0.432, 0.788,
+ 12, 04, 0.597, 0.433, 0.788,
+ 13, 04, 0.597, 0.433, 0.788,
+ 14, 04, 0.598, 0.434, 0.788,
+ 15, 04, 0.598, 0.434, 0.788,
+ 16, 04, 0.598, 0.435, 0.787,
+ 17, 04, 0.599, 0.435, 0.787,
+ 18, 04, 0.599, 0.436, 0.787,
+ 19, 04, 0.599, 0.436, 0.786,
+ 20, 04, 0.599, 0.436, 0.786,
+ 21, 04, 0.600, 0.437, 0.785,
+ 22, 04, 0.600, 0.437, 0.785,
+ 23, 04, 0.600, 0.437, 0.784,
+ 24, 04, 0.600, 0.438, 0.784,
+ 25, 04, 0.600, 0.438, 0.783,
+ 26, 04, 0.601, 0.438, 0.783,
+ 27, 04, 0.601, 0.438, 0.782,
+ 28, 04, 0.601, 0.439, 0.781,
+ 29, 04, 0.601, 0.439, 0.781,
+ 30, 04, 0.601, 0.439, 0.780,
+ 01, 05, 0.601, 0.439, 0.779,
+ 02, 05, 0.601, 0.439, 0.778,
+ 03, 05, 0.601, 0.439, 0.778,
+ 04, 05, 0.601, 0.440, 0.777,
+ 05, 05, 0.601, 0.440, 0.776,
+ 06, 05, 0.601, 0.440, 0.775,
+ 07, 05, 0.601, 0.440, 0.775,
+ 08, 05, 0.601, 0.440, 0.774,
+ 09, 05, 0.601, 0.440, 0.773,
+ 10, 05, 0.602, 0.440, 0.772,
+ 11, 05, 0.602, 0.440, 0.772,
+ 12, 05, 0.602, 0.440, 0.771,
+ 13, 05, 0.602, 0.440, 0.770,
+ 14, 05, 0.602, 0.440, 0.770,
+ 15, 05, 0.602, 0.440, 0.769,
+ 16, 05, 0.602, 0.440, 0.768,
+ 17, 05, 0.602, 0.440, 0.768,
+ 18, 05, 0.602, 0.440, 0.767,
+ 19, 05, 0.602, 0.440, 0.767,
+ 20, 05, 0.602, 0.440, 0.766,
+ 21, 05, 0.602, 0.440, 0.766,
+ 22, 05, 0.602, 0.440, 0.765,
+ 23, 05, 0.602, 0.440, 0.765,
+ 24, 05, 0.602, 0.440, 0.764,
+ 25, 05, 0.602, 0.440, 0.764,
+ 26, 05, 0.602, 0.440, 0.764,
+ 27, 05, 0.602, 0.439, 0.763,
+ 28, 05, 0.602, 0.439, 0.763,
+ 29, 05, 0.602, 0.439, 0.763,
+ 30, 05, 0.602, 0.439, 0.762,
+ 31, 05, 0.602, 0.439, 0.762,
+ 01, 06, 0.602, 0.439, 0.762,
+ 02, 06, 0.602, 0.439, 0.762,
+ 03, 06, 0.602, 0.439, 0.762,
+ 04, 06, 0.602, 0.439, 0.762,
+ 05, 06, 0.602, 0.439, 0.762,
+ 06, 06, 0.602, 0.438, 0.761,
+ 07, 06, 0.602, 0.438, 0.761,
+ 08, 06, 0.602, 0.438, 0.761,
+ 09, 06, 0.602, 0.438, 0.761,
+ 10, 06, 0.602, 0.438, 0.761,
+ 11, 06, 0.602, 0.438, 0.762,
+ 12, 06, 0.602, 0.438, 0.762,
+ 13, 06, 0.602, 0.438, 0.762,
+ 14, 06, 0.602, 0.438, 0.762,
+ 15, 06, 0.602, 0.437, 0.762,
+ 16, 06, 0.602, 0.437, 0.762,
+ 17, 06, 0.602, 0.437, 0.762,
+ 18, 06, 0.602, 0.437, 0.762,
+ 19, 06, 0.602, 0.437, 0.763,
+ 20, 06, 0.602, 0.437, 0.763,
+ 21, 06, 0.602, 0.437, 0.763,
+ 22, 06, 0.602, 0.436, 0.763,
+ 23, 06, 0.602, 0.436, 0.763,
+ 24, 06, 0.602, 0.436, 0.764,
+ 25, 06, 0.602, 0.436, 0.764,
+ 26, 06, 0.601, 0.436, 0.764,
+ 27, 06, 0.601, 0.436, 0.764,
+ 28, 06, 0.601, 0.436, 0.764,
+ 29, 06, 0.601, 0.435, 0.765,
+ 30, 06, 0.601, 0.435, 0.765,
+ 01, 07, 0.601, 0.435, 0.765,
+ 02, 07, 0.600, 0.435, 0.765,
+ 03, 07, 0.600, 0.435, 0.765,
+ 04, 07, 0.600, 0.434, 0.766,
+ 05, 07, 0.600, 0.434, 0.766,
+ 06, 07, 0.599, 0.434, 0.766,
+ 07, 07, 0.599, 0.434, 0.766,
+ 08, 07, 0.599, 0.434, 0.766,
+ 09, 07, 0.598, 0.433, 0.766,
+ 10, 07, 0.598, 0.433, 0.766,
+ 11, 07, 0.598, 0.433, 0.766,
+ 12, 07, 0.597, 0.433, 0.766,
+ 13, 07, 0.597, 0.432, 0.767,
+ 14, 07, 0.597, 0.432, 0.767,
+ 15, 07, 0.596, 0.432, 0.767,
+ 16, 07, 0.596, 0.432, 0.766,
+ 17, 07, 0.595, 0.431, 0.766,
+ 18, 07, 0.595, 0.431, 0.766,
+ 19, 07, 0.594, 0.431, 0.766,
+ 20, 07, 0.594, 0.430, 0.766,
+ 21, 07, 0.593, 0.430, 0.766,
+ 22, 07, 0.593, 0.430, 0.766,
+ 23, 07, 0.592, 0.429, 0.765,
+ 24, 07, 0.592, 0.429, 0.765,
+ 25, 07, 0.591, 0.428, 0.765,
+ 26, 07, 0.590, 0.428, 0.765,
+ 27, 07, 0.590, 0.428, 0.764,
+ 28, 07, 0.589, 0.427, 0.764,
+ 29, 07, 0.588, 0.427, 0.764,
+ 30, 07, 0.588, 0.426, 0.763,
+ 31, 07, 0.587, 0.426, 0.763,
+ 01, 08, 0.586, 0.425, 0.762,
+ 02, 08, 0.586, 0.425, 0.762,
+ 03, 08, 0.585, 0.424, 0.761,
+ 04, 08, 0.584, 0.424, 0.761,
+ 05, 08, 0.583, 0.423, 0.760,
+ 06, 08, 0.583, 0.423, 0.760,
+ 07, 08, 0.582, 0.422, 0.759,
+ 08, 08, 0.581, 0.421, 0.758,
+ 09, 08, 0.580, 0.421, 0.758,
+ 10, 08, 0.579, 0.420, 0.757,
+ 11, 08, 0.578, 0.420, 0.756,
+ 12, 08, 0.578, 0.419, 0.755,
+ 13, 08, 0.577, 0.418, 0.754,
+ 14, 08, 0.576, 0.418, 0.754,
+ 15, 08, 0.575, 0.417, 0.753,
+ 16, 08, 0.574, 0.416, 0.752,
+ 17, 08, 0.573, 0.415, 0.751,
+ 18, 08, 0.572, 0.415, 0.750,
+ 19, 08, 0.571, 0.414, 0.749,
+ 20, 08, 0.570, 0.413, 0.748,
+ 21, 08, 0.569, 0.413, 0.747,
+ 22, 08, 0.569, 0.412, 0.746,
+ 23, 08, 0.568, 0.411, 0.745,
+ 24, 08, 0.567, 0.410, 0.744,
+ 25, 08, 0.566, 0.409, 0.743,
+ 26, 08, 0.565, 0.409, 0.742,
+ 27, 08, 0.564, 0.408, 0.741,
+ 28, 08, 0.563, 0.407, 0.740,
+ 29, 08, 0.562, 0.406, 0.738,
+ 30, 08, 0.561, 0.405, 0.737,
+ 31, 08, 0.560, 0.405, 0.736,
+ 01, 09, 0.558, 0.404, 0.735,
+ 02, 09, 0.557, 0.403, 0.734,
+ 03, 09, 0.556, 0.402, 0.732,
+ 04, 09, 0.555, 0.401, 0.731,
+ 05, 09, 0.554, 0.401, 0.730,
+ 06, 09, 0.553, 0.400, 0.728,
+ 07, 09, 0.552, 0.399, 0.727,
+ 08, 09, 0.551, 0.398, 0.725,
+ 09, 09, 0.550, 0.397, 0.724,
+ 10, 09, 0.549, 0.396, 0.723,
+ 11, 09, 0.548, 0.396, 0.721,
+ 12, 09, 0.546, 0.395, 0.720,
+ 13, 09, 0.545, 0.394, 0.718,
+ 14, 09, 0.544, 0.393, 0.717,
+ 15, 09, 0.543, 0.392, 0.715,
+ 16, 09, 0.542, 0.391, 0.713,
+ 17, 09, 0.541, 0.391, 0.712,
+ 18, 09, 0.540, 0.390, 0.710,
+ 19, 09, 0.538, 0.389, 0.709,
+ 20, 09, 0.537, 0.388, 0.707,
+ 21, 09, 0.536, 0.388, 0.705,
+ 22, 09, 0.535, 0.387, 0.703,
+ 23, 09, 0.533, 0.386, 0.702,
+ 24, 09, 0.532, 0.385, 0.700,
+ 25, 09, 0.531, 0.385, 0.698,
+ 26, 09, 0.530, 0.384, 0.696,
+ 27, 09, 0.528, 0.383, 0.694,
+ 28, 09, 0.527, 0.383, 0.692,
+ 29, 09, 0.526, 0.382, 0.690,
+ 30, 09, 0.525, 0.381, 0.688,
+ 01, 10, 0.523, 0.381, 0.686,
+ 02, 10, 0.522, 0.380, 0.684,
+ 03, 10, 0.521, 0.379, 0.682,
+ 04, 10, 0.519, 0.379, 0.680,
+ 05, 10, 0.518, 0.378, 0.678,
+ 06, 10, 0.517, 0.377, 0.676,
+ 07, 10, 0.515, 0.377, 0.674,
+ 08, 10, 0.514, 0.376, 0.671,
+ 09, 10, 0.512, 0.376, 0.669,
+ 10, 10, 0.511, 0.375, 0.667,
+ 11, 10, 0.510, 0.375, 0.664,
+ 12, 10, 0.508, 0.374, 0.662,
+ 13, 10, 0.507, 0.374, 0.659,
+ 14, 10, 0.505, 0.373, 0.657,
+ 15, 10, 0.504, 0.373, 0.654,
+ 16, 10, 0.502, 0.372, 0.652,
+ 17, 10, 0.501, 0.372, 0.649,
+ 18, 10, 0.499, 0.372, 0.647,
+ 19, 10, 0.498, 0.371, 0.644,
+ 20, 10, 0.496, 0.371, 0.641,
+ 21, 10, 0.495, 0.371, 0.639,
+ 22, 10, 0.493, 0.370, 0.636,
+ 23, 10, 0.492, 0.370, 0.633,
+ 24, 10, 0.490, 0.370, 0.630,
+ 25, 10, 0.489, 0.369, 0.628,
+ 26, 10, 0.487, 0.369, 0.625,
+ 27, 10, 0.485, 0.369, 0.622,
+ 28, 10, 0.484, 0.368, 0.619,
+ 29, 10, 0.482, 0.368, 0.616,
+ 30, 10, 0.481, 0.368, 0.613,
+ 31, 10, 0.479, 0.368, 0.610,
+ 01, 11, 0.478, 0.368, 0.607,
+ 02, 11, 0.476, 0.367, 0.604,
+ 03, 11, 0.475, 0.367, 0.601,
+ 04, 11, 0.473, 0.367, 0.598,
+ 05, 11, 0.471, 0.367, 0.595,
+ 06, 11, 0.470, 0.367, 0.592,
+ 07, 11, 0.468, 0.367, 0.589,
+ 08, 11, 0.467, 0.366, 0.586,
+ 09, 11, 0.465, 0.366, 0.583,
+ 10, 11, 0.464, 0.366, 0.580,
+ 11, 11, 0.462, 0.366, 0.577,
+ 12, 11, 0.461, 0.366, 0.574,
+ 13, 11, 0.459, 0.366, 0.571,
+ 14, 11, 0.458, 0.366, 0.568,
+ 15, 11, 0.456, 0.366, 0.565,
+ 16, 11, 0.455, 0.366, 0.562,
+ 17, 11, 0.454, 0.366, 0.559,
+ 18, 11, 0.452, 0.365, 0.556,
+ 19, 11, 0.451, 0.365, 0.553,
+ 20, 11, 0.450, 0.365, 0.550,
+ 21, 11, 0.448, 0.365, 0.547,
+ 22, 11, 0.447, 0.365, 0.544,
+ 23, 11, 0.446, 0.365, 0.542,
+ 24, 11, 0.445, 0.365, 0.539,
+ 25, 11, 0.443, 0.365, 0.536,
+ 26, 11, 0.442, 0.365, 0.533,
+ 27, 11, 0.441, 0.365, 0.531,
+ 28, 11, 0.440, 0.365, 0.528,
+ 29, 11, 0.439, 0.365, 0.526,
+ 30, 11, 0.438, 0.365, 0.523,
+ 01, 12, 0.437, 0.365, 0.521,
+ 02, 12, 0.436, 0.365, 0.519,
+ 03, 12, 0.435, 0.365, 0.517,
+ 04, 12, 0.434, 0.365, 0.515,
+ 05, 12, 0.434, 0.365, 0.513,
+ 06, 12, 0.433, 0.365, 0.511,
+ 07, 12, 0.432, 0.365, 0.509,
+ 08, 12, 0.431, 0.365, 0.507,
+ 09, 12, 0.431, 0.365, 0.505,
+ 10, 12, 0.430, 0.365, 0.504,
+ 11, 12, 0.430, 0.365, 0.502,
+ 12, 12, 0.429, 0.365, 0.501,
+ 13, 12, 0.429, 0.365, 0.500,
+ 14, 12, 0.429, 0.365, 0.498,
+ 15, 12, 0.428, 0.365, 0.497,
+ 16, 12, 0.428, 0.365, 0.496,
+ 17, 12, 0.428, 0.365, 0.496,
+ 18, 12, 0.428, 0.365, 0.495,
+ 19, 12, 0.428, 0.365, 0.494,
+ 20, 12, 0.428, 0.365, 0.494,
+ 21, 12, 0.428, 0.365, 0.494,
+ 22, 12, 0.428, 0.365, 0.493,
+ 23, 12, 0.429, 0.365, 0.493,
+ 24, 12, 0.429, 0.366, 0.493,
+ 25, 12, 0.429, 0.366, 0.493,
+ 26, 12, 0.430, 0.366, 0.494,
+ 27, 12, 0.430, 0.366, 0.494,
+ 28, 12, 0.431, 0.366, 0.495,
+ 29, 12, 0.431, 0.366, 0.495,
+ 30, 12, 0.432, 0.366, 0.496,
+ 31, 12, 0.433, 0.366, 0.497),
+ ncol = 5, byrow = TRUE,
+ dimnames = list(NULL, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))))
diff --git a/R/UtilsErrorCrit.R b/R/UtilsErrorCrit.R
index 95ca18fce0ec8bbd7bd9aa9ddf807a5250c4bb52..4e200bb28ec42c6378c4ab75ab7cfb59868a66ca 100644
--- a/R/UtilsErrorCrit.R
+++ b/R/UtilsErrorCrit.R
@@ -91,7 +91,7 @@
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
- InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
+ InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
}
if (InputsCrit$transfo == "boxcox") {
muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
diff --git a/inst/vignettesData/vignetteParamOptimCaramel.rda b/inst/vignettesData/vignetteParamOptimCaramel.rda
new file mode 100644
index 0000000000000000000000000000000000000000..7d2e9fcc41220cab22fa5b1e57e7c666ef74803e
Binary files /dev/null and b/inst/vignettesData/vignetteParamOptimCaramel.rda differ
diff --git a/man/CreateCalibOptions.Rd b/man/CreateCalibOptions.Rd
index b418995a04c37002f14ad8bc3bfb6332167f725f..c98ad748b70f67f09ada0082c21502b5358963ac 100644
--- a/man/CreateCalibOptions.Rd
+++ b/man/CreateCalibOptions.Rd
@@ -80,7 +80,7 @@ CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
Users wanting to use \code{FUN_MOD}, \code{FUN_CALIB} or \code{FUN_TRANSFO} functions that are not included in
the package must create their own \code{CalibOptions} object accordingly. \cr
-## ---- CemaNeige version
+## --- CemaNeige version
If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used. \cr
If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr
diff --git a/man/CreateInputsCrit.Rd b/man/CreateInputsCrit.Rd
index 848a35fd1d0a01e5eb3b3c4ed64008671e080d00..eed50f03672fc974f63264c022db962ef7a30433 100644
--- a/man/CreateInputsCrit.Rd
+++ b/man/CreateInputsCrit.Rd
@@ -73,11 +73,11 @@ To calculate composite or multiple criteria, it is necessary to use the \code{Er
\details{
Users wanting to use \code{FUN_CRIT} functions that are not included in the package must create their own InputsCrit object accordingly. \cr \cr
-## ---- Period of calculation
+## --- Period of calculation
Criteria can be calculated over discontinuous periods (i.e. only over winter periods, or when observed discharge is below a certain threshold). To do so, please indicate in \code{Bool_Crit} which indices must be used for calcullation. Discontinuous periods are allowed in the \code{Bool_Crit} argument.
-## ---- Transformations
+## --- Transformations
Transformations are simple functions applied to the observed and simulated variables used in order to change their distribution. Transformations are often used in hydrology for modifying the weight put on errors made for high flows or low flows. The following transformations are available: \cr \cr
\itemize{
@@ -93,11 +93,11 @@ We do not advise computing KGE or KGE' with log-transformation as it might be wr
In order to make sure that KGE and KGE2 remain dimensionless and are not impacted by zero values, the Box-Cox transformation (\code{transfo = "boxcox"}) uses the formulation given in Equation 10 of Santos et al. (2018). Lambda is set to 0.25 accordingly. \cr \cr
The syntax of the power transformation allows a numeric or a string of characters. For example for a squared transformation, the following can be used: \code{transfo = 2}, \code{transfo = "2"} or \code{transfo = "^2"}. Negative values are allowed. Fraction values are not allowed (e.g., \code{"-1/2"} must instead be written \code{"-0.5"}).\cr \cr
-## ---- The epsilon value
+## --- The epsilon value
The epsilon value is useful when \code{"log"} or \code{"inv"} transformations are used (to avoid calculation of the inverse or of the logarithm of zero). If an epsilon value is provided, then it is added to the observed and simulated variable time series at each time step and before the application of a transformation. The epsilon value has no effect when the \code{"boxcox"} transformation is used. The impact of this value and a recommendation about the epsilon value to use (usually one hundredth of average observation) are discussed in Pushpalatha et al. (2012) for NSE and in Santos et al. (2018) for KGE and KGE'. \cr \cr
-## ---- Single, multiple or composite criteria calculation
+## --- Single, multiple or composite criteria calculation
Users can set the following arguments as atomic or list: \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, \code{Weights}. If the list format is chosen, all the lists must have the same length. \cr
Calculation of a single criterion (e.g. NSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments atomics only. \cr
diff --git a/man/CreateRunOptions.Rd b/man/CreateRunOptions.Rd
index 6ed188c4670de2f1944cf2dc00ae4af3b979ce2b..20859d4c03a58c843d4c3b7c103d9d0dda56b0db 100644
--- a/man/CreateRunOptions.Rd
+++ b/man/CreateRunOptions.Rd
@@ -72,11 +72,11 @@ CreateRunOptions(FUN_MOD, InputsModel,
Users wanting to use \code{FUN_MOD} functions that are not included in
the package must create their own \code{RunOptions} object accordingly.
-## ---- IndPeriod_WarmUp and IndPeriod_Run
+## --- IndPeriod_WarmUp and IndPeriod_Run
Since the hydrological models included in airGR are continuous models, meaning that internal states of the models are propagated to the next time step, \code{IndPeriod_WarmUp} and \code{IndPeriod_Run} must be continuous periods, represented by continuous indices values; no gaps are allowed. To calculate criteria or to calibrate a model over discontinuous periods, please see the \code{Bool_Crit} argument of the \code{\link{CreateInputsCrit}} function.
-## ---- Initialisation options
+## --- Initialisation options
The model initialisation options can either be set to a default configuration or be defined by the user.
@@ -111,7 +111,7 @@ However, it is also possible to perform a long-term initialisation if other indi
}
}
-## ---- CemaNeige version
+## --- CemaNeige version
If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used. \cr
If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr
diff --git a/man/Param_Sets_GR4J.Rd b/man/Param_Sets_GR4J.Rd
index 7a3abb7cb111b98e0b79055443cf256aa96a8111..5798b2952e164a8e9ea759dca97d6c8bca928d24 100644
--- a/man/Param_Sets_GR4J.Rd
+++ b/man/Param_Sets_GR4J.Rd
@@ -56,7 +56,7 @@ Param_Sets_GR4J <- as.matrix(Param_Sets_GR4J)
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
-## ---- calibration step
+## --- calibration step
## short calibration period selection (< 6 months)
Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
@@ -82,7 +82,7 @@ OutputsCrit_Loop <- apply(Param_Sets_GR4J, 1, function(Param) {
Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ])
-## ---- validation step
+## --- validation step
## validation period selection
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-03-01"),
diff --git a/man/RunModel_CemaNeige.Rd b/man/RunModel_CemaNeige.Rd
index 2fc511d5641a82f178833e331779eb2b9cb2ff2a..685cb2a6d5eacb7961b52a0a4fdb3645e30cfb85 100644
--- a/man/RunModel_CemaNeige.Rd
+++ b/man/RunModel_CemaNeige.Rd
@@ -79,7 +79,7 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-0
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-## ---- original version of CemaNeige
+## --- original version of CemaNeige
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel,
@@ -94,7 +94,7 @@ OutputsModel <- RunModel_CemaNeige(InputsModel = InputsModel,
plot(OutputsModel)
-## ---- version of CemaNeige with the Linear Hysteresis
+## --- version of CemaNeige with the Linear Hysteresis
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel,
diff --git a/man/RunModel_CemaNeigeGR4H.Rd b/man/RunModel_CemaNeigeGR4H.Rd
index 7dbe52223adbac99b1419c6f246fbe48f1cded56..d2a499da73a6f7faf9e7b8fbff32f9ac4b3ec7df 100644
--- a/man/RunModel_CemaNeigeGR4H.Rd
+++ b/man/RunModel_CemaNeigeGR4H.Rd
@@ -108,7 +108,7 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H:\%M")=="2
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H:\%M")=="2008-12-31 23:00"))
-## ---- original version of CemaNeige
+## --- original version of CemaNeige
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4H, InputsModel = InputsModel,
diff --git a/man/RunModel_CemaNeigeGR4J.Rd b/man/RunModel_CemaNeigeGR4J.Rd
index c451c6a3f2be63b5ba5c29736b287358fe363b01..9ab766451012e0d66c42401884b4f54c795fd080 100644
--- a/man/RunModel_CemaNeigeGR4J.Rd
+++ b/man/RunModel_CemaNeigeGR4J.Rd
@@ -104,7 +104,7 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-0
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-## ---- original version of CemaNeige
+## --- original version of CemaNeige
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel,
@@ -125,7 +125,7 @@ InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsMo
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
-## ---- version of CemaNeige with the Linear Hysteresis
+## --- version of CemaNeige with the Linear Hysteresis
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel,
diff --git a/man/RunModel_CemaNeigeGR5H.Rd b/man/RunModel_CemaNeigeGR5H.Rd
index 5b9f07e5e0008631fe6a381e8f23bcf6faff6861..e7c3f8a9826af875b8064ba7884800926204410b 100644
--- a/man/RunModel_CemaNeigeGR5H.Rd
+++ b/man/RunModel_CemaNeigeGR5H.Rd
@@ -109,7 +109,7 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H:\%M")=="2
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H:\%M")=="2008-12-31 23:00"))
-## ---- original version of CemaNeige
+## --- original version of CemaNeige
## Imax computation
Imax <- Imax(InputsModel = InputsModel, IndPeriod_Run = Ind_Run,
diff --git a/man/RunModel_CemaNeigeGR6J.Rd b/man/RunModel_CemaNeigeGR6J.Rd
index 65e9635d50cf00bc89d0421c6830654f867ddddb..e9a37d4de87e101a6135302f1d0d97e33dfd808c 100644
--- a/man/RunModel_CemaNeigeGR6J.Rd
+++ b/man/RunModel_CemaNeigeGR6J.Rd
@@ -108,7 +108,7 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-0
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-## ---- original version of CemaNeige
+## --- original version of CemaNeige
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel,
@@ -129,7 +129,7 @@ InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsMo
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
-## ---- version of CemaNeige with the Linear Hysteresis
+## --- version of CemaNeige with the Linear Hysteresis
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel,
diff --git a/man/RunModel_GR1A.Rd b/man/RunModel_GR1A.Rd
index 56b2bce2c475ab9ab23219ba10580efe8a2ada9d..0ebba8426ced4aa870354ba128ca21c9d2839e11 100644
--- a/man/RunModel_GR1A.Rd
+++ b/man/RunModel_GR1A.Rd
@@ -60,7 +60,7 @@ TabSeries <- data.frame(DatesR = BasinObs$DatesR,
P = BasinObs$P,
E = BasinObs$E,
Qmm = BasinObs$Qmm)
-TabSeries <- TabSeries[TabSeries$DatesR < "2012-09-01", ]
+TabSeries <- TabSeries[TabSeries$DatesR < as.POSIXct("2012-09-01", tz = "UTC"), ]
BasinObs <- SeriesAggreg(TabSeries, Format = "\%Y",
YearFirstMonth = 09,
ConvertFun = c("sum", "sum", "sum"))
diff --git a/man/TransfoParam.Rd b/man/TransfoParam.Rd
index de58b74b8908342e621d973d77e5e7eaef589440..2029ed45c330e662c20d34a163f33e9a3febfe49 100644
--- a/man/TransfoParam.Rd
+++ b/man/TransfoParam.Rd
@@ -55,7 +55,7 @@ TransfoParam_CemaNeigeHyst(ParamIn, Direction)
\examples{
library(airGR)
-## ---- generic function
+## --- generic function
## transformation Raw -> Transformed for the GR4J model
Xraw <- matrix(c(+221.41, -3.63, +30.00, +1.37,
@@ -72,7 +72,7 @@ Xtran <- matrix(c(+3.60, -2.00, +3.40, -9.10,
Xraw <- TransfoParam(ParamIn = Xtran, Direction = "TR", FUN_TRANSFO = TransfoParam_GR4J)
-## ---- specific function
+## --- specific function
## transformation Raw -> Transformed for the GR4J model
Xraw <- matrix(c(+221.41, -3.63, +30.00, +1.37,
diff --git a/src/airGR.c b/src/airGR.c
index b4c256b242438470673b6a80c2f0b2ff26505b80..b1b549cc67dc039674ac8592d2662f1816111fef 100644
--- a/src/airGR.c
+++ b/src/airGR.c
@@ -2,7 +2,7 @@
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>
-/* FIXME:
+/* FIXME:
Check these declarations against the C/Fortran source code.
*/
@@ -26,7 +26,7 @@ static const R_FortranMethodDef FortranEntries[] = {
{"frun_gr4j", (DL_FUNC) &F77_NAME(frun_gr4j), 11},
{"frun_gr5j", (DL_FUNC) &F77_NAME(frun_gr5j), 11},
{"frun_gr6j", (DL_FUNC) &F77_NAME(frun_gr6j), 11},
- {"frun_pe_oudin", (DL_FUNC) &F77_NAME(frun_pe_oudin), 5},
+ {"frun_pe_oudin", (DL_FUNC) &F77_NAME(frun_pe_oudin), 5},
{NULL, NULL, 0}
};
diff --git a/src/frun_GR5H.f90 b/src/frun_GR5H.f90
index 2027f47cee451b0a239d868974bfffb6ff109357..e5ec1a1cbd991fe81b03c7cbb1fea59e60de3daf 100644
--- a/src/frun_GR5H.f90
+++ b/src/frun_GR5H.f90
@@ -76,8 +76,11 @@
doubleprecision, dimension(NMISC) :: MISC
doubleprecision :: D,P1,E,Q
- IF (Imax .LT. 0.) IsIntStore = .FALSE.
- IF (Imax .GE. 0.) IsIntStore = .TRUE.
+ IF (Imax .LT. 0.d0) THEN
+ IsIntStore = .FALSE.
+ ELSE
+ IsIntStore = .TRUE.
+ ENDIF
!--------------------------------------------------------------
! Initializations
diff --git a/tests/testthat/helper_regression.R b/tests/testthat/helper_regression.R
index cbcf9dbfee09f0309b429ad54636123ba988a34a..bc05c5521a19946aa2cef6170f7f28f830dfefe2 100644
--- a/tests/testthat/helper_regression.R
+++ b/tests/testthat/helper_regression.R
@@ -1,5 +1,5 @@
StoreStableExampleResults <- function(
- package = "airGR",
+ package = "airGR",
path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "stable"),
...) {
install.packages(package, repos = "http://cran.r-project.org")
@@ -7,8 +7,8 @@ StoreStableExampleResults <- function(
}
StoreDevExampleResults <- function(
- package = "airGR",
- path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "dev"),
+ package = "airGR",
+ path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "dev"),
...) {
StoreExampleResults(package = package, path = path, ...)
}
@@ -75,7 +75,9 @@ StoreTopicResults <- function(topic, package, path, run.dontrun = TRUE, run.dont
}
CompareStableDev <- function() {
- res = testthat::test_file("tests/testthat/regression.R")
- dRes = as.data.frame(res)
- if(any(dRes[,"failed"]>0) | any(dRes[,"error"])) quit(status = 1)
+ res <- testthat::test_file("tests/testthat/regression.R")
+ dRes <- as.data.frame(res)
+ if(any(dRes[, "failed"] > 0) | any(dRes[, "error"])) {
+ quit(status = 1)
+ }
}
diff --git a/tests/testthat/helper_vignettes.R b/tests/testthat/helper_vignettes.R
index 8dfe5e12d84bd9cf03acf5db2a37efbfbe638fc1..931f99a55ee90fa4417040a36b58f66220aaac71 100644
--- a/tests/testthat/helper_vignettes.R
+++ b/tests/testthat/helper_vignettes.R
@@ -3,12 +3,12 @@
#' @param fileRmd Rmd file to
#' @param tmpFolder Folder storing the script containing extracted chunks
#' @param force.eval Force execution of chunks with parameter eval=FALSE
-RunRmdChunks <- function(fileRmd,
- tmpFolder = "../tmp",
- force.eval = TRUE) {
+RunRmdChunks <- function(fileRmd,
+ tmpFolder = "../tmp",
+ force.eval = TRUE) {
dir.create(tmpFolder, showWarnings = FALSE)
output <- file.path(tmpFolder,
- gsub("\\.Rmd", "\\.R", basename(fileRmd), ignore.case = TRUE))
+ gsub("\\.Rmd", "\\.R", basename(fileRmd), ignore.case = TRUE))
knitr::purl(fileRmd, output = output, quiet = TRUE)
sTxt <- readLines(output)
if (force.eval) {
@@ -30,8 +30,8 @@ RunRmdChunks <- function(fileRmd,
for (i in 1:length(chunksEvalStart)) {
# Remove comments on eval=F chunk lines
sTxt[chunksEvalStart[i]:chunksEvalEnd[i]] <- gsub(pattern = "^## ",
- replace = "",
- x = sTxt[chunksEvalStart[i]:chunksEvalEnd[i]])
+ replace = "",
+ x = sTxt[chunksEvalStart[i]:chunksEvalEnd[i]])
}
}
@@ -70,12 +70,12 @@ RunRmdChunks <- function(fileRmd,
RunVignetteChunks <- function(vignette,
tmpFolder = "../tmp",
force.eval = TRUE) {
- if(file.exists(file.path("../../vignettes/", paste0(vignette, ".Rmd")))) {
+ if(file.exists(sprintf("../../vignettes/%s.Rmd", vignette))) {
# testthat context in development environnement
- RunRmdChunks(file.path("../../vignettes/", paste0(vignette, ".Rmd")), tmpFolder, force.eval)
+ RunRmdChunks(sprintf("../../vignettes/%s.Rmd", vignette), tmpFolder, force.eval)
} else {
# R CMD check context in package environnement
- RunRmdChunks(system.file(file.path("doc/", paste0(vignette, ".Rmd")), package = "airGR"), tmpFolder, force.eval)
+ RunRmdChunks(system.file(sprintf("doc/%s.Rmd", vignette), package = "airGR"), tmpFolder, force.eval)
}
return(TRUE)
}
@@ -92,4 +92,4 @@ TestQmmQlsConversion <- function(BasinObs, BasinArea, tolerance = 1E-7) {
Conversion <- Conversion / 86400 # Day -> seconds
notNA <- which(!is.na(BasinObs$Qmm))
expect_equal(BasinObs$Qmm[notNA] * Conversion, BasinObs$Qls[notNA], tolerance = tolerance)
-}
\ No newline at end of file
+}
diff --git a/tests/testthat/regression.R b/tests/testthat/regression.R
index 2ccadba442cf1ab03b3834b7e28a6533e60f21b7..1732b2c544430557ada828b159717e8e0d8e2179 100644
--- a/tests/testthat/regression.R
+++ b/tests/testthat/regression.R
@@ -35,7 +35,7 @@ if (dir.exists(file.path(tmp_path, "stable")) & dir.exists(file.path(tmp_path, "
message("File ", file.path(getwd(), regIgnoreFile), " not found")
regIgnore <- NULL
}
- lapply(X = refVarFiles, CompareWithStable, testDir = file.path(tmp_path, "dev"), regIgnore = regIgnore)
+ lapply(refVarFiles, FUN = CompareWithStable, testDir = file.path(tmp_path, "dev"), regIgnore = regIgnore)
} else {
stop("Regression tests compared to released version needs that you run the following instructions first:\n",
"Rscript tests/testthat/regression_tests.R stable\n",
diff --git a/tests/testthat/regression_tests.R b/tests/testthat/regression_tests.R
index 2b59829af88508a115abe5c20df9f0c8a7810107..4456e09edc85aaa6b06bd36f3b518d41c192e786 100644
--- a/tests/testthat/regression_tests.R
+++ b/tests/testthat/regression_tests.R
@@ -1,9 +1,9 @@
# Execute Regression test by comparing RD files stored in folders /tests/tmp/ref and /tests/tmp/test
-Args = commandArgs(trailingOnly=TRUE)
+Args <- commandArgs(trailingOnly = TRUE)
source("tests/testthat/helper_regression.R")
-lActions = list(
+lActions <- list(
stable = StoreStableExampleResults,
dev = StoreDevExampleResults,
compare = CompareStableDev
@@ -15,7 +15,7 @@ if(Args %in% names(lActions)) {
stop("This script should be run with one argument in the command line:\n",
"`Rscript tests/regression_tests.R [stable|dev|compare]`.\n",
"Available arguments are:\n",
- "- stable: install stable version from CRAN, run and store examples\n",
+ "- stable: install stable version from CRAN, run and store examples\n",
"- dev: install dev version from current directory, run and store examples\n",
"- compare: stored results of both versions")
}
diff --git a/tests/testthat/test-CreateRunOptions.R b/tests/testthat/test-CreateRunOptions.R
index 91723330598bf17fe8d5dca46a00074d62d92168..08c877dd62acedeceb0f30ad9a01ac5e4a582aae 100644
--- a/tests/testthat/test-CreateRunOptions.R
+++ b/tests/testthat/test-CreateRunOptions.R
@@ -2,26 +2,29 @@ context("CreateRunOptions")
test_that("Warm start of GR4J should give same result as warmed model", {
data(L0123001)
- InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+ InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
- Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
- which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
- Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"),
+ Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
+ which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
+ Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-12-31"))
# 1990-1991
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
+ InputsModel = InputsModel,
+ IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
+ InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run1))
OutputsModel1 <- RunModel_GR4J(InputsModel = InputsModel,
- RunOptions = RunOptions1, Param = Param)
+ RunOptions = RunOptions1, Param = Param)
# Warm start 1991
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel, IndPeriod_Run = Ind_Run2,
+ InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
OutputsModel2 <- RunModel_GR4J(InputsModel = InputsModel,
diff --git a/tests/testthat/test-CreateiniStates.R b/tests/testthat/test-CreateiniStates.R
index 58dbcf8a7ee9fd26967338907048774f82bd0ee5..17ac824be4cdb090bfdca6c7d2ddd3b8dc089a54 100644
--- a/tests/testthat/test-CreateiniStates.R
+++ b/tests/testthat/test-CreateiniStates.R
@@ -3,25 +3,23 @@ context("CreateIniStates on SD model")
data(L0123001)
test_that("Error: SD argument provided on non-SD 'InputsModel'", {
- InputsModel <-
- CreateInputsModel(
- FUN_MOD = RunModel_GR4J,
- DatesR = BasinObs$DatesR,
- Precip = BasinObs$P,
- PotEvap = BasinObs$E
- )
+ InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E
+ )
expect_error(
- IniStates <-
- CreateIniStates(
- FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel,
- ProdStore = 0,
- RoutStore = 0,
- ExpStore = NULL,
- UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
- UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
- SD = list(rep(0, 10))
- ),
+ IniStates <- CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = list(rep(0, 10))
+ ),
regexp = "'SD' argument provided and"
)
})
@@ -29,10 +27,9 @@ test_that("Error: SD argument provided on non-SD 'InputsModel'", {
BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
-Qupstream <-
- floor((sin((
- seq_along(BasinObs$Qmm) / 365 * 2 * 3.14
- )) + 1) * mean(BasinObs$Qmm, na.rm = TRUE))
+Qupstream <- floor((sin((
+ seq_along(BasinObs$Qmm) / 365 * 2 * 3.14
+)) + 1) * mean(BasinObs$Qmm, na.rm = TRUE))
InputsModel <- CreateInputsModel(
FUN_MOD = RunModel_GR4J,
@@ -46,8 +43,27 @@ InputsModel <- CreateInputsModel(
test_that("Error: Non-list 'SD' argument", {
expect_error(
- IniStates <-
- CreateIniStates(
+ IniStates <- CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = rep(0, 10)
+ ),
+ regexp = "'SD' argument must be a list"
+ )
+})
+
+test_that("Error: Non-numeric items in 'SD' list argument", {
+ lapply(list(list(list(rep(
+ 0, 10
+ ))), list(toto = NULL)),
+ function(x) {
+ expect_error(
+ IniStates <- CreateIniStates(
FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
ProdStore = 0,
@@ -55,47 +71,27 @@ test_that("Error: Non-list 'SD' argument", {
ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
- SD = rep(0, 10)
+ SD = x
),
- regexp = "'SD' argument must be a list"
- )
-})
-
-test_that("Error: Non-numeric items in 'SD' list argument", {
- lapply(list(list(list(rep(0, 10))), list(toto = NULL)),
- function(x) {
- expect_error(
- IniStates <-
- CreateIniStates(
- FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel,
- ProdStore = 0,
- RoutStore = 0,
- ExpStore = NULL,
- UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
- UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
- SD = x
- ),
- regexp = "Each item of 'SD' list argument must be numeric"
- )
- })
+ regexp = "Each item of 'SD' list argument must be numeric"
+ )
+ })
})
test_that("Error: Number of items not equal to number of upstream connections", {
lapply(list(list(), list(rep(0, 10), rep(0, 10))),
function(x) {
expect_error(
- IniStates <-
- CreateIniStates(
- FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel,
- ProdStore = 0,
- RoutStore = 0,
- ExpStore = NULL,
- UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
- UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
- SD = x
- ),
+ IniStates <- CreateIniStates(
+ FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ ProdStore = 0,
+ RoutStore = 0,
+ ExpStore = NULL,
+ UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+ UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+ SD = x
+ ),
regexp = "list argument must be the same as the number of upstream"
)
})
diff --git a/tests/testthat/test-RunModel_LAG.R b/tests/testthat/test-RunModel_Lag.R
similarity index 91%
rename from tests/testthat/test-RunModel_LAG.R
rename to tests/testthat/test-RunModel_Lag.R
index 00d0104afc576bf336f5730af9d11df00ceeff85..bfef15d2833076d4dc30005615c46d9ed00eb2fe 100644
--- a/tests/testthat/test-RunModel_LAG.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -51,8 +51,8 @@ Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModel,
- IndPeriod_Run = Ind_Run))
+ InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run))
test_that("InputsModel parameter should contain an OutputsModel key", {
expect_error(
@@ -140,10 +140,8 @@ test_that("Params from calibration with simulated data should be similar to init
InputsModel = InputsModel,
RunOptions = RunOptions,
VarObs = "Q",
- Obs = (
- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
- BasinObs$Qmm[Ind_Run] * BasinAreas[2L]
- ) / sum(BasinAreas)
+ Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
+ BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
)
CalibOptions <- CreateCalibOptions(
FUN_MOD = RunModel_GR4J,
@@ -193,14 +191,19 @@ Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"
# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = IM, IndPeriod_Run = Ind_Run1))
+ InputsModel = IM,
+ IndPeriod_Run = Ind_Run1))
OutputsModel1 <- RunModel(InputsModel = IM,
- RunOptions = RunOptions1, Param = PSDini, FUN_MOD = RunModel_GR4J)
+ RunOptions = RunOptions1, Param = PSDini,
+ FUN_MOD = RunModel_GR4J)
# 1990-1991
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = IM, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
+ InputsModel = IM,
+ IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
OutputsModel <- RunModel(InputsModel = IM,
- RunOptions = RunOptions, Param = PSDini, FUN_MOD = RunModel_GR4J)
+ RunOptions = RunOptions,
+ Param = PSDini,
+ FUN_MOD = RunModel_GR4J)
test_that("Warm start should give same result as warmed model", {
# Warm start 1991
@@ -209,7 +212,9 @@ test_that("Warm start should give same result as warmed model", {
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
OutputsModel2 <- RunModel(InputsModel = IM,
- RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
+ RunOptions = RunOptions2,
+ Param = PSDini,
+ FUN_MOD = RunModel_GR4J)
# Compare 1991 Qsim from warm started and from 1990-1991
names(OutputsModel2$Qsim) <- NULL
expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
diff --git a/tests/testthat/test-SeriesAggreg.R b/tests/testthat/test-SeriesAggreg.R
index c9edc42ee8d8ba01b8a4bb7f517f096dd0e5ade7..b65428c241067c26eca8557d376aa9025610945f 100644
--- a/tests/testthat/test-SeriesAggreg.R
+++ b/tests/testthat/test-SeriesAggreg.R
@@ -61,9 +61,10 @@ test_that("Check SeriesAggreg output values on yearly aggregation", {
E = BasinObs$E,
Qmm = BasinObs$Qmm
)
- GoodValues <- apply(BasinObs[BasinObs$DatesR >= "1984-09-01" &
- BasinObs$DatesR < "1985-09-01",
- c("P", "E", "Qmm")], 2, sum)
+ GoodValues <- apply(BasinObs[BasinObs$DatesR >= as.POSIXct("1984-09-01", tz = "UTC") &
+ BasinObs$DatesR < as.POSIXct("1985-09-01", tz = "UTC"),
+ c("P", "E", "Qmm")],
+ MARGIN = 2, FUN = sum)
TestedValues <- unlist(SeriesAggreg(TabSeries,
Format = "%Y",
YearFirstMonth = 9,
@@ -229,7 +230,7 @@ test_that("SeriesAggreg should work with ConvertFun 'min', 'max' and 'median'",
Qls <- BasinObs[, c("DatesR", "Qls")]
test_ConvertFunRegime <- function(x, ConvertFun, TimeFormat) {
expect_equal(nrow(SeriesAggreg(x, TimeFormat, ConvertFun = ConvertFun)),
- length(unique(format(BasinObs$DatesR, "%Y"))))
+ length(unique(format(BasinObs$DatesR, "%Y"))))
}
lapply(c("max", "min", "median"), function(x) {test_ConvertFunRegime(Qls, x, "%Y")})
})
diff --git a/tests/testthat/test-evap.R b/tests/testthat/test-evap.R
index aac8122d134e540315c3f4c986defad81de0d452..e68a47159c5e744cf5566b7b5e4241495ed99f61 100644
--- a/tests/testthat/test-evap.R
+++ b/tests/testthat/test-evap.R
@@ -6,7 +6,7 @@ comp_evap <- function(BasinObs,
TimeStepOut = "daily") {
PotEvap <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
Temp = BasinObs$T,
- Lat = Lat, LatUnit = LatUnit,
+ Lat = Lat, LatUnit = LatUnit,
TimeStepIn = TimeStepIn, TimeStepOut = TimeStepOut)
PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
Temp = BasinObs$T,
@@ -19,7 +19,7 @@ comp_evap <- function(BasinObs,
test_that("PE_Oudin works", {
skip_on_cran()
rm(list = ls())
-
+
data(L0123001); BasinObs_L0123001 <- BasinObs
data(L0123002); BasinObs_L0123002 <- BasinObs
@@ -30,14 +30,14 @@ test_that("PE_Oudin works", {
Lat = 0.8, LatUnit = "rad",
TimeStepIn = "daily", TimeStepOut = "hourly"))
expect_true(comp_evap(BasinObs = BasinObs_L0123002,
- Lat = 0.9, LatUnit = "rad",
+ Lat = 0.9, LatUnit = "rad",
TimeStepIn = "daily", TimeStepOut = "daily"))
expect_true(comp_evap(BasinObs = BasinObs_L0123002,
Lat = 0.9, LatUnit = "rad",
TimeStepIn = "daily", TimeStepOut = "hourly"))
-
+
## check with several catchments using different values for Lat
-
+
## one by one
PotEvapFor1 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123001$DatesR)$yday + 1,
Temp = BasinObs_L0123001$T,
@@ -47,7 +47,7 @@ test_that("PE_Oudin works", {
Temp = BasinObs_L0123002$T,
Lat = 0.9, LatUnit = "rad",
RunFortran = TRUE)
-
+
## all in one
BasinObs_L0123001$Lat <- 0.8
BasinObs_L0123002$Lat <- 0.9
@@ -56,7 +56,6 @@ test_that("PE_Oudin works", {
Temp = BasinObs$T,
Lat = BasinObs$Lat, LatUnit = "rad",
RunFortran = TRUE)
-
+
expect_equal(PotEvapFor, c(PotEvapFor1, PotEvapFor2))
-
})
diff --git a/tests/testthat/test-vignettes.R b/tests/testthat/test-vignettes.R
index 7d19c74f1ad2ecefb51218928d45128e19b733af..88de5d118cbdd54a94c6bd298065df54db0a664a 100644
--- a/tests/testthat/test-vignettes.R
+++ b/tests/testthat/test-vignettes.R
@@ -11,11 +11,14 @@ test_that("V02.1_param_optim works", {
skip_on_cran()
rm(list = ls())
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+ load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
rda_resGLOB <- resGLOB
rda_resPORT <- resPORT
+ rda_optMO <- optMO
expect_true(RunVignetteChunks("V02.1_param_optim"))
- expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1E-7)
- expect_equal(resGLOB[,-1], rda_resGLOB[,-1], tolerance = 1E-2) # High tolerance due to randomisation in optimisations
+ expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1e-7)
+ expect_equal(resGLOB[, -1], rda_resGLOB[, -1], tolerance = 1e-2) # High tolerance due to randomisation in optimisations
+ expect_equal(summary(optMO$parameters), summary(rda_optMO$parameters), tolerance = 1e-7)
})
test_that("V02.2_param_mcmc works", {
@@ -25,15 +28,14 @@ test_that("V02.2_param_mcmc works", {
rda_gelRub <- gelRub
rda_multDRAM <- multDRAM
expect_true(RunVignetteChunks("V02.2_param_mcmc"))
- expect_equal(gelRub, rda_gelRub, tolerance = 1E-7)
- expect_equal(multDRAM, rda_multDRAM, tolerance = 1E-7)
+ expect_equal(gelRub, rda_gelRub, tolerance = 1e-7)
+ expect_equal(multDRAM, rda_multDRAM, tolerance = 1e-7)
})
test_that("V03_param_sets_GR4J works", {
skip_on_cran()
rm(list = ls())
expect_true(RunVignetteChunks("V03_param_sets_GR4J"))
-
})
test_that("V04_cemaneige_hysteresis works", {
@@ -46,8 +48,8 @@ test_that("V04_cemaneige_hysteresis works", {
rda_OutputsCrit_Val_NoHyst <- OutputsCrit_Val_NoHyst
expect_true(RunVignetteChunks("V04_cemaneige_hysteresis"))
TestQmmQlsConversion(BasinObs, BasinInfo$BasinArea)
- expect_equal(OutputsCrit_Cal, rda_OutputsCrit_Cal, tolerance = 1E-7)
- expect_equal(OutputsCrit_Cal_NoHyst, rda_OutputsCrit_Cal_NoHyst, tolerance = 1E-7)
- expect_equal(OutputsCrit_Val, rda_OutputsCrit_Val, tolerance = 1E-7)
- expect_equal(OutputsCrit_Val_NoHyst, rda_OutputsCrit_Val_NoHyst, tolerance = 1E-7)
+ expect_equal(OutputsCrit_Cal, rda_OutputsCrit_Cal, tolerance = 1e-7)
+ expect_equal(OutputsCrit_Cal_NoHyst, rda_OutputsCrit_Cal_NoHyst, tolerance = 1e-7)
+ expect_equal(OutputsCrit_Val, rda_OutputsCrit_Val, tolerance = 1e-7)
+ expect_equal(OutputsCrit_Val_NoHyst, rda_OutputsCrit_Val_NoHyst, tolerance = 1e-7)
})
diff --git a/vignettes/V01_get_started.Rmd b/vignettes/V01_get_started.Rmd
index feecacbda9c5da715c3e7b324c33073550b2d095..3c5eb159f74fda55ecdf2de6231b9cef7297c3b8 100644
--- a/vignettes/V01_get_started.Rmd
+++ b/vignettes/V01_get_started.Rmd
@@ -1,5 +1,6 @@
---
title: "Get Started with airGR"
+author: "Guillaume Thirel, Olivier Delaigue, Laurent Coron"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index 0063d9151e5c4c90de00a1d61aaa0374ebac4e31..2dd121e92d89e6838e3c7db3af73322e74080400 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -1,6 +1,6 @@
---
title: "Plugging in new calibration algorithms in airGR"
-author: "François Bourgin"
+author: "François Bourgin, Guillaume Thirel"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
@@ -15,9 +15,13 @@ library(airGR)
library(DEoptim)
library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5
library(Rmalschains)
+library(caRamel)
+library(ggplot2)
+library(GGally)
# source("airGR.R")
set.seed(321)
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
```
@@ -158,7 +162,7 @@ As it can be seen in the table below, the four additional optimization strategie
```{r, warning=FALSE, echo=FALSE, eval=FALSE}
resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
round(rbind(
- OutputsCalib$ParamFinalR ,
+ OutputsCalib$ParamFinalR,
airGR::TransfoParam_GR4J(ParamIn = optPORT$par , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par) , Direction = "TR"),
@@ -172,4 +176,100 @@ resGLOB
<!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->
+# Multiobjective optimization
+
+Multiobjective optimization is used to explore possible trade-offs between different performances criteria.
+Here we use the following R implementation of an efficient strategy:
+* [caRamel: Automatic Calibration by Evolutionary Multi Objective Algorithm](https://cran.r-project.org/package=caRamel)
+
+Motivated by using the rainfall-runoff model for low flow simulation, we explore the trade-offs between the KGE values obtained without any data transformation and with the inverse transformation.
+
+First, the OptimGR4J function previously used is modified to return two values.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+InputsCrit_inv <- InputsCrit
+InputsCrit_inv$transfo <- "inv"
+
+MOptimGR4J <- function(i) {
+ if (algo == "caRamel") {
+ ParamOptim <- x[i, ]
+ }
+ ## Transformation of the parameter set to real space
+ RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
+ Direction = "TR")
+ ## Simulation given a parameter set
+ OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Param = RawParamOptim)
+ ## Computation of the value of the performance criteria
+ OutputsCrit1 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit,
+ OutputsModel = OutputsModel,
+ verbose = FALSE)
+ ## Computation of the value of the performance criteria
+ OutputsCrit2 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit_inv,
+ OutputsModel = OutputsModel,
+ verbose = FALSE)
+ return(c(OutputsCrit1$CritValue, OutputsCrit2$CritValue))
+}
+```
+
+## caRamel
+caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+algo <- "caRamel"
+optMO <- caRamel::caRamel(nobj = 2,
+ nvar = 4,
+ minmax = rep(TRUE, 2),
+ bounds = matrix(c(lowerGR4J, upperGR4J), ncol = 2),
+ func = MOptimGR4J,
+ popsize = 100,
+ archsize = 100,
+ maxrun = 15000,
+ prec = rep(1.e-3, 2),
+ carallel = FALSE,
+ graph = FALSE)
+```
+
+The algorithm returns parameter sets that describe the pareto front, illustrating the trade-off between overall good performance and good performance for low flow.
+
+```{r, fig.width=6, fig.height=6, warning=FALSE}
+ggplot() +
+ geom_point(aes(optMO$objectives[, 1], optMO$objectives[, 2])) +
+ coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) +
+ xlab("KGE") + ylab("KGE [1/Q]") +
+ theme_bw()
+```
+
+The pameter sets can be viewed in the parameter space, illustrating different populations.
+
+```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
+param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+ airGR::TransfoParam_GR4J(x, Direction = "TR")
+ })
+GGally::ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
+```
+
+```{r fig.height=6, fig.width=12, message=FALSE, warning=FALSE}
+RunOptions$Outputs_Sim <- "Qsim"
+run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+ airGR::RunModel_GR4J(InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ Param = x)
+ }$Qsim)
+run_optMO <- data.frame(run_optMO)
+
+ggplot() +
+ geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+ y = run_optMO$X1)) +
+ geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+ y = run_optMO$X54),
+ colour = "darkred") +
+ scale_x_datetime(limits = c(as.POSIXct("1998-01-01"), NA)) +
+ ylab("Discharge [mm/d]") + xlab("Date") +
+ scale_y_sqrt() +
+ theme_bw()
+```
+
+
diff --git a/vignettes/V03_param_sets_GR4J.Rmd b/vignettes/V03_param_sets_GR4J.Rmd
index 8f9897dd0da11359cedea62393753860458dc14b..0db718558de3dbe391bab579703ede5c6220700f 100644
--- a/vignettes/V03_param_sets_GR4J.Rmd
+++ b/vignettes/V03_param_sets_GR4J.Rmd
@@ -1,5 +1,6 @@
---
title: "Generalist parameter sets for the GR4J model"
+author: "Olivier Delaigue, Guillaume Thirel"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
diff --git a/vignettes/V04_cemaneige_hysteresis.Rmd b/vignettes/V04_cemaneige_hysteresis.Rmd
index 8c87387dd98abc9e93c6676c39ddfa2f082160ea..a100eb06b8acfb2c252f0fc042a980bc355ec8ed 100644
--- a/vignettes/V04_cemaneige_hysteresis.Rmd
+++ b/vignettes/V04_cemaneige_hysteresis.Rmd
@@ -1,5 +1,6 @@
---
title: "Using satellite snow cover area data for calibrating and improving CemaNeige"
+author: "Guillaume Thirel, Olivier Delaigue"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 5e814abc10a965416d4aef1a78e95bb6631893b1..f970f603124a8b78ba73fbe23f29d0d5db27fa16 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -69,10 +69,12 @@ InputsModelUp <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$Da
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModelUp, IndPeriod_Run = Ind_Run,
- IniStates = NULL, IniResLevels = NULL, IndPeriod_WarmUp = NULL)
+ InputsModel = InputsModelUp
+ , IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
+ IniStates = NULL, IniResLevels = NULL)
InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelUp,
- RunOptions = RunOptionsUp, VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
+ RunOptions = RunOptionsUp,
+ VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
@@ -101,9 +103,9 @@ we need to create the `InputsModel` object completed with upstream information:
InputsModelDown1 <- CreateInputsModel(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
- Qupstream = matrix(QObsUp, ncol = 1), # Upstream observed flow
- LengthHydro = 100 * 1000, # Distance between upstream catchment outlet and the downstream one in m
- BasinAreas = c(180, 180) # Upstream and downstream areas in km²
+ Qupstream = matrix(QObsUp, ncol = 1), # upstream observed flow
+ LengthHydro = 1e2 * 1e3, # distance between upstream catchment outlet & the downstream one [m]
+ BasinAreas = c(180, 180) # upstream and downstream areas [km²]
)
```
@@ -111,15 +113,19 @@ And then calibrate the combination of Lag model for upstream flow transfer and G
```{r}
RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
- InputsModel = InputsModelDown1, IndPeriod_Run = Ind_Run,
- IniStates = NULL, IniResLevels = NULL, IndPeriod_WarmUp = NULL)
+ InputsModel = InputsModelDown1,
+ IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
+ IniStates = NULL, IniResLevels = NULL)
InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelDown1,
- RunOptions = RunOptionsDown, VarObs = "Q", Obs = QObsDown[Ind_Run])
+ RunOptions = RunOptionsDown,
+ VarObs = "Q", Obs = QObsDown[Ind_Run])
CalibOptionsDown <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel,
- IsSD = TRUE) # Don't forget to specify that it's an SD model here
-OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1, RunOptions = RunOptionsDown,
- InputsCrit = InputsCritDown, CalibOptions = CalibOptionsDown,
+ IsSD = TRUE) # specify that it's a SD model
+OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
+ RunOptions = RunOptionsDown,
+ InputsCrit = InputsCritDown,
+ CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
```
@@ -151,8 +157,10 @@ CritDown1 <- ErrorCrit_NSE(InputsCritDown, OutputsModelDown1)
We calibrate the model with the `InputsModel` object previously created for substituting the observed upstream flow with the simulated one:
```{r}
-OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2, RunOptions = RunOptionsDown,
- InputsCrit = InputsCritDown, CalibOptions = CalibOptionsDown,
+OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
+ RunOptions = RunOptionsDown,
+ InputsCrit = InputsCritDown,
+ CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
ParamDown2 <- OutputsCalibDown2$ParamFinalR
```
@@ -162,7 +170,7 @@ ParamDown2 <- OutputsCalibDown2$ParamFinalR
## Identification of Lag parameter
-The theoretical LAG parameter should be equal to:
+The theoretical Lag parameter should be equal to:
```{r}
Lag <- InputsModelDown1$LengthHydro / (2 * 86400)
@@ -172,10 +180,14 @@ paste(format(Lag), "m/s")
Both calibrations overestimate this parameter:
```{r}
-mLag <- matrix(c(Lag, OutputsCalibDown1$ParamFinalR[1], OutputsCalibDown2$ParamFinalR[1]), ncol = 1)
-rownames(mLag) = c("theoretical", "calibrated with observed upstream flow",
- "calibrated with simulated upstream flow")
-colnames(mLag) = c("Lag parameter")
+mLag <- matrix(c(Lag,
+ OutputsCalibDown1$ParamFinalR[1],
+ OutputsCalibDown2$ParamFinalR[1]),
+ ncol = 1,
+ dimnames = list(c("theoretical",
+ "calibrated with observed upstream flow",
+ "calibrated with simulated upstream flow"),
+ c("Lag parameter")))
knitr::kable(mLag)
```
@@ -197,11 +209,17 @@ CritDownTheo <- ErrorCrit_NSE(InputsCritDown, OutputsModelDownTheo)
## Parameters and performance of each subcatchment for all calibrations
```{r}
-comp <- matrix(c(0, OutputsCalibUp$ParamFinalR, rep(OutputsCalibDown1$ParamFinalR, 2),
- OutputsCalibDown2$ParamFinalR, ParamDownTheo), ncol = 5, byrow = TRUE)
-comp <- cbind(comp, c(OutputsCalibUp$CritFinal, OutputsCalibDown1$CritFinal,
- CritDown1$CritValue, OutputsCalibDown2$CritFinal, CritDownTheo$CritValue))
-colnames(comp) <- c("Lag", paste0("x", 1:4), "NSE")
+comp <- matrix(c(0, OutputsCalibUp$ParamFinalR,
+ rep(OutputsCalibDown1$ParamFinalR, 2),
+ OutputsCalibDown2$ParamFinalR,
+ ParamDownTheo),
+ ncol = 5, byrow = TRUE)
+comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
+ OutputsCalibDown1$CritFinal,
+ CritDown1$CritValue,
+ OutputsCalibDown2$CritFinal,
+ CritDownTheo$CritValue))
+colnames(comp) <- c("Lag", paste0("X", 1:4), "NSE")
rownames(comp) <- c("Calibration of the upstream subcatchment",
"Calibration 1 with observed upstream flow",
"Validation 1 with simulated upstream flow",