Commit 7a97389c authored by Dorchies David's avatar Dorchies David
Browse files

test(calibration): add tests for all hourly models

Refs #120
parent db303278
Pipeline #23644 failed with stages
in 31 minutes and 44 seconds
......@@ -12,91 +12,131 @@ sModels <- c(
"CemaNeigeGR6J FALSE L0123001 NA 184.9341841;0.5551637;59.7398917;2.2177177;0.4760000;6.0496475;0.0000000;14.4642868",
"CemaNeigeGR4J TRUE L0123001 NA 208.5127103;0.5781516;102.5140641;2.2274775;0.0000000;6.7644613;8.0000000;1.0000000",
"CemaNeigeGR5J TRUE L0123001 NA 202.350228;0.901525;98.494430;1.788288;0.483984;0.000000;7.401500;6.100000;1.000000",
"CemaNeigeGR6J TRUE L0123001 NA 188.67010241;0.56662930;60.34028760;2.22747748;0.47600000;5.98945247;0.03203203;7.93816892;10.80000000;1.00000000"
"CemaNeigeGR6J TRUE L0123001 NA 188.67010241;0.56662930;60.34028760;2.22747748;0.47600000;5.98945247;0.03203203;7.93816892;10.80000000;1.00000000",
"GR4H FALSE L0123003 NA 711.676649;-1.158469;150.556095;4.686093",
"GR5H FALSE L0123003 NA 804.0021672;-0.1898488;137.7524699;3.0436628;0.1951163",
"CemaNeigeGR4H FALSE L0123003 NA 1.595685e+03;-8.183484e-01;2.320697e+02;5.000000e-01;5.005005e-04;9.342369e+01",
"CemaNeigeGR5H FALSE L0123003 NA 33.34921883;-4.98925432;332.00673122;1.58534106;0.20792716;0.02214393;4.28498513",
"CemaNeigeGR4H TRUE L0123003 NA 1.766316e+03;-6.920667e-01;2.192034e+02;3.451688e+00;5.005005e-0;4.869585e+01;1.111447e+0;5.064090e-01",
"CemaNeigeGR5H TRUE L0123003 NA 66.6863310;-1.4558128;138.3795123;2.6499450;0.2325000;0.0000000;0.3017014;48.4000000;0.9914915"
)
dfModels <- read.table(text = paste(sModels, collapse = "\n"), header = TRUE)
dates <- c("1985-01-01", "1985-12-31", "1986-01-01", "2012-12-31")
TestModelCalibration <- function(model) {
model <- as.list(model)
test_that(paste(model$name, "works"), {
skip_on_cran()
sModel <- paste0("RunModel_", model$name)
ModelCalibration <- function(model) {
sModel <- paste0("RunModel_", model$name)
sIM_FUN_MOD <- sModel
if(model$data == "L0123003") {
# hourly time step database
dates <- c("2004-01-01 00:00", "2004-12-31 23:00", "2005-01-01 00:00", "2008-12-31 23:00")
date_format = "%Y-%m-%d %H:%M"
TempMean <- fakeHourlyTemp()
} else {
# yearly, monthly, daily time step databases
dates <- c("1985-01-01", "1985-12-31", "1986-01-01", "2012-12-31")
date_format <- "%Y-%m-%d"
if(!is.na(model$aggreg)) {
sIM_FUN_MOD <- "RunModel_GR4J"
} else {
sIM_FUN_MOD <- sModel
# Aggregation on monthly and yearly databases
sIM_FUN_MOD <- "RunModel_GR4J" # CreateInputsModel with daily data
date_format <- model$aggreg
}
}
## loading catchment data
data(list = model$data)
## loading catchment data
data(list = model$data)
if(model$data != "L0123003") TempMean <- BasinObs$T
# preparation of the InputsModel object with daily time step data
InputsModel <- CreateInputsModel(FUN_MOD = sIM_FUN_MOD,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
TempMean = BasinObs$T,
ZInputs = median(BasinInfo$HypsoData),
HypsoData = BasinInfo$HypsoData,
NLayers = 5)
# preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = sIM_FUN_MOD,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
TempMean = TempMean,
ZInputs = median(BasinInfo$HypsoData),
HypsoData = BasinInfo$HypsoData,
NLayers = 5)
if(!is.na(model$aggreg)) {
# conversion of InputsModel to target time step
InputsModel <- SeriesAggreg(InputsModel, Format = model$aggreg)
dfQobs <- SeriesAggreg(data.frame(DatesR = BasinObs$DatesR, Qmm = BasinObs$Qmm),
Format = model$aggreg, ConvertFun = "sum")
Obs <- dfQobs$Qmm
} else {
Obs <- BasinObs$Qmm
}
if(!is.na(model$aggreg)) {
# conversion of InputsModel to target time step
InputsModel <- SeriesAggreg(InputsModel, Format = model$aggreg)
# calibration period selection
if(is.na(model$aggreg)) {
date_format <- "%Y-%m-%d"
} else {
date_format <- model$aggreg
}
dates <- sapply(dates, function(x) format(as.Date(x), format = date_format))
Ind_WarmUp <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[1]),
which(format(InputsModel$DatesR, format = date_format)==dates[2])
)
Ind_Run <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[3]),
which(format(InputsModel$DatesR, format = date_format)==dates[4])
)
dfQobs <- SeriesAggreg(data.frame(DatesR = BasinObs$DatesR, Qmm = BasinObs$Qmm),
Format = model$aggreg, ConvertFun = "sum")
Obs <- dfQobs$Qmm
} else {
Obs <- BasinObs$Qmm
}
# calibration period selection
dates <- sapply(dates, function(x) format(as.Date(x), format = date_format))
Ind_WarmUp <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[1]),
which(format(InputsModel$DatesR, format = date_format)==dates[2])
)
Ind_Run <- seq(
which(format(InputsModel$DatesR, format = date_format)==dates[3]),
which(format(InputsModel$DatesR, format = date_format)==dates[4])
)
# preparation of the RunOptions object
suppressWarnings(
RunOptions <- CreateRunOptions(
FUN_MOD = sModel,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = Ind_WarmUp,
IsHyst = as.logical(model$IsHyst)
)
# preparation of the RunOptions object
suppressWarnings(
RunOptions <- CreateRunOptions(
FUN_MOD = sModel,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = Ind_WarmUp,
IsHyst = as.logical(model$IsHyst)
)
)
# calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = Obs[Ind_Run])
# preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(sModel, IsHyst = as.logical(model$IsHyst))
# calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = Obs[Ind_Run])
# preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(sModel, IsHyst = as.logical(model$IsHyst))
# calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = sModel)
expect_equal(OutputsCalib$ParamFinalR,
# calibration
suppressWarnings(OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = sModel))
OutputsCalib$ParamFinalR
}
TestModelCalibration <- function(model) {
model <- as.list(model)
test_that(paste(model$name, "works"), {
skip_on_cran()
ParamFinalR <- ModelCalibration(model)
expect_equal(ParamFinalR,
as.numeric(strsplit(model$ParamFinalR, ";")[[1]]),
tolerance = 1E-6)
})
# OutputsCalib$ParamFinalR
}
#' Create Fake hourly temperature from daily temperatures in L0123001
#'
#' @param start_date [character] start date in format "%Y-%m-%d"
#' @param end_date [character] end date in format "%Y-%m-%d"
#' @return [numeric] hourly temperature time series between `start_date` and `end_date`
fakeHourlyTemp <- function(start_date = "2004-01-01", end_date = "2008-12-31") {
dates <- as.POSIXct(c(start_date, end_date), tz = "UTC")
data(L0123002)
indJ <- seq.int(which(BasinObs$DatesR == as.POSIXct(dates[1])),
which(BasinObs$DatesR == as.POSIXct(dates[2])))
TJ <- BasinObs$T[indJ]
TH <- approx((seq.int(length(TJ)) - 1) * 24,TJ,
seq.int(length(TJ) * 24 ) - 1,
rule = 2)$y
varT_1J <- -sin(0:23/24 * 2 * pi) # Temp min at 6 and max at 18
varT <- rep(varT_1J, length(TJ))
TH <- TH + varT * 5 # For a mean daily amplitude of 10°
TH
}
apply(dfModels, 1, TestModelCalibration)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment