Failed to fetch fork details. Try again later.
-
Delaigue Olivier authored2dec31be
Forked from
HYCAR-Hydro / airGR
Source project has a limited visibility.
context("Calibration")sModels <- c( "name IsHyst data aggreg ParamFinalR", "GR1A FALSE L0123001 %Y 0.91125", "GR2M FALSE L0123001 %Y%m 259.8228;0.9975", "GR4J FALSE L0123001 NA 223.6315877;0.5781516;97.5143942;2.2177177", "GR5J FALSE L0123001 NA 220.3863609;0.8944531;93.5640705;1.7628720;0.4846427", "GR6J FALSE L0123001 NA 192.8761657;0.6933087;49.1783293;2.2145422;0.5088240;6.8146261", "CemaNeigeGR4J FALSE L0123001 NA 2.043839e+02;5.781516e-01;1.025141e+02;2.217718e+00;1.501502e-03;1.432036e+01", "CemaNeigeGR5J FALSE L0123001 NA 1.983434e+02;8.747758e-01;9.849443e+01;1.768769e+00;4.829830e-01;1.501502e-03;1.432036e+01", "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", "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-04;4.869585e+01;1.111447e+01;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)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)) { # 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) if(model$data != "L0123003") TempMean <- BasinObs$T # 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 } # calibration period selection7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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)
)
)
# 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
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, ifelse(as.logical(model$IsHyst), "Hysteresis", ""), "works"), {
skip_on_cran()
ParamFinalR <- ModelCalibration(model)
expect_equal(ParamFinalR,
as.numeric(strsplit(model$ParamFinalR, ";")[[1]]),
tolerance = 1E-6)
})
}
#' 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
}
141142143
apply(dfModels, 1, TestModelCalibration)