diff --git a/tests/testthat/test-Calibration.R b/tests/testthat/test-Calibration.R index b3b2821aa909fccdaa8c6bdd1ec44d8b74ca6962..594fc078f8e7a7822dea0bff88fe1e8f7d4322d0 100644 --- a/tests/testthat/test-Calibration.R +++ b/tests/testthat/test-Calibration.R @@ -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)