Calage avec différentes variantes du Critère KGE2

library(seinebasin2)
cfg <- loadConfig()
reservoirs <- "with"
Crit <- c("ErrorCrit_KGE2sqrt", "ErrorCrit_KGE2", "ErrorCrit_KGE3", "ErrorCrit_KGE4", "ErrorCrit_KGE5")
BasinsObs <- seinebasin2::loadBasinsObs("historical/BasinsObs_observations_day_1958-2019.RDS", cfg = cfg)

Calage des paramètres en fonction des variantes de la fonction ErrorCritKGE2

Calage des paramètres grâce à 4 fonctions coût, ErrorCrit_KGE2 voir [airGRiwrm], ErrorCrit_KGE3 prenant en compte 50% de KGE2[Q] + 50% de KGE2[1/Q], ErrorCrit_KGE4 combiant le 90% du KGE2[Q] et 10% du KGE calculé sur les débits maximum annuel moyens entre Qsim et Qobs. Voir [Error_CRIT.R] in seinebasin2/dev et enfin ErrorCrit_KGE5 définis comme 50% de KGE2[Q²] + 50% de KGE2[1/Q]. Afin de comparer les différentes variantes nous réutilisons la fonction coût utilisé à l’origine pour SeineBasin2 soit ErrorCrit_KGE2 en utilisant sqrt(Q).

Calibration_multipleErrorCrit(BasinObs, Crit, reservoirs, transfo = "")

Simulation avec les paramètres optimisés obtenus

param_list.path <- list.files(getDataPath(cfg$calibration$path, cfg$Version))
param_list.path <- param_list.path[grep(paste0(reservoirs,".csv"), param_list.path)]

OutputsModels.list <- lapply(param_list.path, function(x) {
    BasinsObsM <- BasinsObs
    if (reservoirs == "with") {
        BasinsObsM$griwrm <- addReservoirsGRiwrm(BasinsObsM$griwrm)
        Crit_FUN <- substr(gsub("Calib_|_with.csv", "", x), 1, 14)
        transfo <- gsub(paste0("Calib_|", Crit_FUN, "|_with.csv"), "", x)
        BasinsObsM$Q <-
            addReservoirsQ(BasinsObsM, FUN = Crit_FUN, transfo = transfo)
    }
    RunModel(BasinsObsM, x)
})
names(OutputsModels.list) <- gsub("Calib_|_with.csv", "", param_list.path)

Sauvegarde des plots issue de la simulation

IndPeriod_Run <- seq(366, length(BasinsObs$DatesR))
lapply(names(OutputsModels.list), function(x) {
    saveGRiwrmOMplot(
        OutputsModels.list[[x]],
        Qobs <- BasinsObs[["Q"]][IndPeriod_Run, ],
        savedir <- getDataPath(cfg$calibration$path, cfg$Version, 
                              paste0("plots/OM_", reservoirs),"ModelOutputs", x)
    )
})

Extraction et calcul des indicateurs hydrologiques

Nous nous intéresserons dans un premier temps aux indicateurs hydrologiques caractérisant les forts débits, QA, QJXA2, QJXA10, QJXA20 puis aux indicateurs caractérisant les étiages, QMNA2, QMNA5, QMNA10, VCN10_2, VCN10_5, VCN10_10, VCN30_2, VCN30_5 et VCN30_10.

indicators <- c("QA",
    "VCN10_2","VCN30_2","QMNA2",
    "VCN10_5","VCN30_5","QMNA5",
    "VCN10_10","VCN30_10","QMNA10",
    "QJXA2","QJXA10","QJXA20"
)

Préparation des débit simulés influencés

Calcule des indicateurs hydrologiques et sauvegarde des indicateurs calculés dans une dataframe comportant pour chaqu’un, l’indicateur observé, simulé et le rapport Sim/Obs.

Qobs <- cbind(DatesR = BasinsObs$DatesR, BasinsObs$Q)
lapply(names(OutputsModels.list), function(x){
    Qsim <-
        cbind(DatesR = OutputsModels.list[[x]][[1]]$DatesR, as.data.frame(do.call(
            cbind, lapply(OutputsModels.list[[x]], "[[", "Qsim")
        )))
    Qobs <- Qobs[Qobs$DatesR %in% Qsim$DatesR,]
    Qsim <- Qsim[, names(Qobs)]
    Qsim[is.na(Qobs)] <- NA
    
    l <- lapply(indicators, calcQIndicators, Qobs = Qobs, Qsim = Qsim)
    mQc <- do.call(cbind, l)
    dfQc <- cbind(data.frame(Id = rownames(mQc)), mQc)

    #save indicators data.frame
    write.table(
    dfQc,
    file <- getDataPath(cfg$calibration$path, cfg$Version, paste0("Indic_", x, "_", reservoirs, ".tsv")),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE)
})

Sauvegarde des cartes pour les différents indicateurs hydrologiques

lapply(Crit, function(x){
    saveMAPindicators(indicators, reservoirs, Crit=x, cfg=cfg)
})

Création des radials plot pour les différentes stations et fonctions coût

Stations = griwrm$id

Q_indi_Calib.list <- list.files(getDataPath(cfg$calibration$path, cfg$Version))
Q_indi_Calib.list <- Q_indi_Calib.list[grep(paste0("_",reservoirs, ".tsv"), Q_indi_Calib.list)]
Q_indi_Calib.list <- lapply(Q_indi_Calib.list, function(x){
    readr::read_tsv(file = getDataPath(cfg$calibration$path, cfg$Version, x))})

radial_plot <- lapply(Stations, function(y) {
    indiCalib_spider <- do.call(data.frame, lapply(Q_indi_Calib.list, function(x) {
        df <- x[x$Id == y, ]
        df <- t(df[, which(colnames(df) %in% paste0("R-", indicators))])
    }))
    
    desc_KGE <-
        c("KGE2: KGE[Q]",
          "KGE2: KGE[sqrtQ]",
          "KGE3: 0.5KGE[Q]+0.5KGE[1/Q]",
          "KGE4: 0.9KGE[Q]+0.1KGE[QJXAm]",
          "KGE5: 0.5KGE[Q²]+0.5KGE[1/Q]"
        )
    colnames(indiCalib_spider) <- desc_KGE
    return(indiCalib_spider)
})
names(radial_plot) <- Stations

Exemlple du basin H5920010, la Seine à Paris Austerlitz

Rdplot <- IndicatorsRadialPlot(
    radial_plot[["H5920010"]],
    axisName = indicators,
    axisMin  = 0,
    axisMax  = 1.1,
    title = title(sub = paste0("[Q] ", "H5920010", "...")),
    MAR = c(5, 5, 5, 5),
    legend = TRUE,
    LegName = "ErrorCrit"
)

Sauvegardes des radials plot

lapply(names(radial_plot), function(x) {
    if (cfg$data$write_results == TRUE) {
        file.path = file.path(
            getDataPath(cfg$calibration$path),
            "v1.0",
            "plots",
            paste0("OM_", reservoirs),
            "Indicators_ResumeSpider",
            paste0("Compared_indicators_ErrorCrits_",
                   x, ".png")
        )
        
        if (cfg$data$write_results == TRUE) {
            png(
                file.path,
                height = 2000,
                width = 2500,
                units = "px",
                pointsize = 50
            )
        }
        
        Rdplot <- IndicatorsRadialPlot(
            radial_plot[[x]],
            axisName = indicators,
            axisMin  = 0,
            axisMax  = 1.1,
            title = title(sub = paste0("[Q] ", x, "...")),
            MAR = c(5, 5, 5, 5),
            legend = TRUE,
            LegName = "ErrorCrit"
        )
        
        if (cfg$data$write_results == TRUE) {
            dev.off()
        }
    }
})