fea: calculate all Explore2070 indicators for influenced, observed and simulated flows

Refs #5, #12

NAMESPACE

@@ -23,11 +23,23 @@ export(addReservoirsQ)
 export(areaPoly)
 export(calcQA)
 export(calcQJXA)
+export(calcQJXA10)
+export(calcQJXA2)
+export(calcQJXA20)
 export(calcQJXAn)
 export(calcQLogNn)
 export(calcQMNA)
+export(calcQMNA10)
+export(calcQMNA2)
+export(calcQMNA5)
 export(calcQMNAn)
 export(calcVCN)
+export(calcVCN10_10)
+export(calcVCN10_2)
+export(calcVCN10_5)
+export(calcVCN30_10)
+export(calcVCN30_2)
+export(calcVCN30_5)
 export(calcVCNn)
 export(createBasinsObs)
 export(getAprioriIds)
@@ -44,6 +56,7 @@ export(max_na_rm)
 export(mean_na_rm)
 export(min_na_rm)
 export(not_na_count)
+export(plot_seine_map)
 export(pmedian)
 export(readDrias2020)
 export(readQsim)

R/calcQJXA.R

@@ -207,8 +207,14 @@ getGumbelParams <- function(x) {
 }
 
 
-#' Compute the flow for a given Gumbel law
+#' Yearly maximum daily flow for a return period
 #'
+#' @description
+#' The maximum daily flow is calculated from an adjustement on a Gumbel law.
+#'
+#' Functions named calcQJXA_\[return_period\] are also defined for commun indicators (2, 10, 20 years).
+#'
+#' @param x [data.frame] of flows (the first column should be a [POSIXct]), or a `QJXA` object (See [calcQJXA])
 #' @param gumbel A [list] with `a` the location parameter and `b` the scale parameter
 #' @param return_period [numeric] Return period. Its unit is given by the aggregation time step used for the Gumbel law.
 #'
@@ -286,3 +292,17 @@ calcQJXAn.GumbelParams <- function(x, return_period) {
 max_na_rm <- function(x) {
     suppressWarnings(max(x, na.rm = TRUE))
 }
+
+
+#' @export
+#' @rdname calcQJXAn
+calcQJXA2 <- function(x) calcQJXAn(x, 2)
+
+#' @export
+#' @rdname calcQJXAn
+calcQJXA10 <- function(x) calcQJXAn(x, 10)
+
+#' @export
+#' @rdname calcQJXAn
+calcQJXA20 <- function(x) calcQJXAn(x, 20)
+

R/calcQMNA.R

@@ -156,7 +156,9 @@ calcQMNA.POSIXt <- function(x, flows) {
 #'
 #' Calculation of QMNA according to a return period. An adjustment to a Galton distribution (log-normal distribution) is realised.
 #'
-#' @param x [numeric] vector of mean monthly flows
+#' Functions named calcQJMNA_\[return_period\] are also defined for commun indicators (2, 5, 10 years).
+#'
+#' @param x [data.frame] of flows (the first column should be a [POSIXct]), or a `QMNA` object (See [calcQMNA])
 #' @param return_period [numeric] number of years
 #' @return [numeric] QMNA for the given return period
 #' @rdname calcQMNAn
@@ -231,3 +233,15 @@ calcQLogNn <- function(x, return_period) {
 min_na_rm <- function(x) {
     suppressWarnings(min(x, na.rm = TRUE))
 }
+
+#' @rdname calcQMNAn
+#' @export
+calcQMNA2 <- function(x) calcQMNAn(x, 2)
+
+#' @rdname calcQMNAn
+#' @export
+calcQMNA5 <- function(x) calcQMNAn(x, 5)
+
+#' @rdname calcQMNAn
+#' @export
+calcQMNA10 <- function(x) calcQMNAn(x, 10)

R/calcVCN.R

@@ -32,10 +32,13 @@ calcVCN.numeric <- function(x, k, ...) {
 
 #' VCNn for a given return period
 #'
+#' @description
 #' Calculation of VCN according to a return period. An adjustment to a Galton distribution (log-normal distribution) is realised.
 #'
+#' Functions named calcVCN\[k\]_\[return_period\] are also defined for commun indicators.
+#'
 #' @param return_period [numeric] number of years
-#' @param x [data.frame] of flows. The first column should be a [POSIXct]
+#' @param x [data.frame] of flows (the first column should be a [POSIXct])
 #' @return [numeric] VCN for the given return period
 #' @rdname calcVCNn
 #' @export
@@ -65,3 +68,27 @@ calcVCNn.data.frame <- function(x, k, return_period, ...) {
     VCNn <- apply(VCN_years[, -1, drop = FALSE], 2, calcQLogNn, return_period = return_period)
     return(VCNn)
 }
+
+#' @rdname calcVCNn
+#' @export
+calcVCN10_2 <- function(x) calcVCNn(x, 10, 2)
+
+#' @rdname calcVCNn
+#' @export
+calcVCN10_5 <- function(x) calcVCNn(x, 10, 5)
+
+#' @rdname calcVCNn
+#' @export
+calcVCN10_10 <- function(x) calcVCNn(x, 10, 10)
+
+#' @rdname calcVCNn
+#' @export
+calcVCN30_2 <- function(x) calcVCNn(x, 30, 2)
+
+#' @rdname calcVCNn
+#' @export
+calcVCN30_5 <- function(x) calcVCNn(x, 30, 5)
+
+#' @rdname calcVCNn
+#' @export
+calcVCN30_10 <- function(x) calcVCNn(x, 30, 10)

R/plot_seine_map.R

+#' Plot choroplethe map of the Seine River Basin at BVI scale
+#'
+#' @param r named [numeric] [vector], the indicator to display with the names corresponding to gauging station IDs
+#' @param breaks see the parameter in [mapsf::mf_map]
+#' @param title see the parameter in [mapsf::mf_layout]
+#' @param credits see the parameter in [mapsf::mf_layout]
+#' @param ... parameters sent to [mapsf::mf_map]
+#'
+#' @return
+#' @export
+#'
+#' @examples
+plot_seine_map <- function(r, breaks, title, credits = "IGN, Inrae", ...) {
+
+    library(mapsf)
+    data("gis_bvi")
+    sf_bvi <- st_as_sf(gis_bvi)
+    sf_bvi$ratio <- r[sf_bvi$CODE]
+
+    #mf_init(x = st_union(sf_bvi), "brutal")
+    mf_map(x = sf_bvi,
+           var = "ratio",
+           type = "choro",
+           pal = "RdYlBu",
+           breaks = breaks,
+           ...)
+
+    # Add river network
+    mf_map(gis_rivers, type = "base", col = "royalblue", add = TRUE)
+
+    # Add main cities
+    data(gis_villes)
+    mf_label(x = gis_villes, var = "X1",
+             cex = 0.5, halo = TRUE, r = 0.10)
+
+    # Add label for lakes
+    lac <- data.frame(x = c(780355.2, 760259.1, 749809.1, 719263.0),
+                      y = c(2401371, 2376452, 2366806, 2238994),
+                      name = c("Lac Marne", "Lac Aube", "Lac Seine", "Lac Yonne"),
+                      zz = 1)
+    sf_lac <- st_as_sf(lac, coords = c("x", "y"), crs = 27572)
+    mf_map(sf_lac, var = "zz", type = "symb", pal = "red", leg_pos = NA)
+    mf_label(x = sf_lac, var = "name", col = "navyblue",
+             cex = 0.65, halo = TRUE, r = 0.10, pos = 4)
+
+    legend(
+        "bottomleft",
+        legend = "Lacs reservoirs",
+        pch = 22,
+        pt.lwd = 1,
+        pt.cex = 1.2,
+        pt.bg = "lightgrey",
+        col = "red",
+        cex = 0.8,
+        bty = "n",
+        inset = c(0.2, 0.06)
+    )
+    legend(
+        "bottomleft",
+        legend = "Hydrographie",
+        lty = 1,
+        seg.len = 1,
+        col = "royalblue",
+        cex = 0.8,
+        bty = "n",
+        inset = c(0.18, 0)
+    )
+
+    mf_layout(title = title,
+              credits = credits)
+}

bookdown/02-calage.Rmd

@@ -3,6 +3,9 @@
 ```{r setup, include=FALSE}
 library(seinebasin2)
 cfg <- loadConfig()
+knitr::opts_chunk$set(
+  fig.width = 8
+)
 ```
 
 Le calage du modèle hydrologique semi-distribué utilise en entrée les précipitations, températures et évapotranspiration extraites de la base SAFRAN [@vidal50yearHighresolutionAtmospheric2010] et aggrégée à l'échelle des sous-bassins versants intermédiaires définis par @nuneztorresSimulationBassinVersant2021.
@@ -34,7 +37,7 @@ OC <-
 
 ## Calage des stations du bassin avec prise en compte des réservoirs {#calage-avec-reservoirs}
 
-Pour caler le modèle avec l'influence des réservoirs, on ajoute les prises et les restitutions des réservoirs à la structure du modèle ainsi que les débits observés à ces prises et restitutions. 
+Pour caler le modèle avec l'influence des réservoirs, on ajoute les prises et les restitutions des réservoirs à la structure du modèle ainsi que les débits observés à ces prises et restitutions.
 
 La "prise d'eau" du lac de Pannecière est calculée en fonction de l'apport naturel arrivant dans le lac. Pour cela, on effectue un calage et un simulation des débits à l'amont du lac de Pannecière à l'aide de la fonction `addReservoirsQ`.
 
@@ -87,15 +90,26 @@ plot(OM_with[[station]], Qobs = BasinsObs$Q[IndPeriod_Run, station], which = "Re
 Les indicateurs utilisés pour l'évaluation sont ceux utilisés lors du projet Explore 2070 [@chazotExplore2070Lot2012].
 
 ```{r}
-getMapRatio <- function(Qsim, Qobs, FUN, varTitle) {
+# Prepare influenced simulated flows for indicator calculations
+Qobs <- cbind(DatesR = BasinsObs$DatesR, BasinsObs$Q)
+Qsim <- cbind(DatesR = OM_with[[1]]$DatesR, as.data.frame(do.call(cbind, lapply(OM_with, "[[", "Qsim"))))
+# Qobs should have same rows as Qsim
+Qobs <- Qobs[Qobs$DatesR %in% Qsim$DatesR, ]
+# Qsim should have same columns as Qobs
+Qsim <- Qsim[, names(Qobs)]
+# Qsim should have same gaps as Qobs
+Qsim[is.na(Qobs)] <- NA
+```
+
+```{r}
+calcRatio <- function(Qsim, Qobs, FUN) {
     pSim <- FUN(Qsim)
     pObs <- FUN(Qobs)
-    r <- pSim / pObs
-    
-    library(mapsf)
-    data("gis_bvi")
-    sf_bvi <- st_as_sf(gis_bvi)
-    sf_bvi$ratio <- r[sf_bvi$CODE]
+    pSim / pObs
+}
+
+plotRatio <- function(Qsim, Qobs, FUN, title) {
+    r <- calcRatio(Qsim, Qobs, FUN)
     
     # Define scale centered around 1
     decVal <- r[r < 1]
@@ -103,97 +117,67 @@ getMapRatio <- function(Qsim, Qobs, FUN, varTitle) {
     decQuant <- quantile(decVal, probs = c(0, 0.1, 0.4, 0.6, 0.9))
     incQuant <- quantile(incVal, probs = c(0.1, 0.4, 0.6, 0.9, 1))
     breaks <- c(decQuant, incQuant)
-    #mf_init(x = st_union(sf_bvi), "brutal")
-    mf_map(x = sf_bvi, 
-           var = "ratio", 
-           type = "choro", 
-           pal = "RdYlBu", 
-           breaks = breaks)
-    
-    # Add river network
-    mf_map(gis_rivers, type = "base", col = "royalblue", add = TRUE)
     
-    # Add main cities
-    data(gis_villes)
-    mf_label(x = gis_villes, var = "X1",
-         cex = 0.5, halo = TRUE, r = 0.10)
-    
-    # Add label for lakes
-    lac <- data.frame(x = c(780355.2, 760259.1, 749809.1, 719263.0),
-                      y = c(2401371, 2376452, 2366806, 2238994),
-                      name = c("Lac Marne", "Lac Aube", "Lac Seine", "Lac Yonne"),
-                      zz = 1)
-    sf_lac <- st_as_sf(lac, coords = c("x", "y"), crs = 27572)
-    mf_map(sf_lac, var = "zz", type = "symb", pal = "red", leg_pos = NA)
-    mf_label(x = sf_lac, var = "name", col = "navyblue",
-             cex = 0.65, halo = TRUE, r = 0.10, pos = 4)
-    
-    legend(
-        "bottomleft",
-        legend = "Lacs reservoirs",
-        pch = 22,
-        pt.lwd = 1,
-        pt.cex = 1.2,
-        pt.bg = "lightgrey",
-        col = "red",
-        cex = 0.8,
-        bty = "n",
-        inset = c(0.2, 0.06)
-    )
-    legend(
-        "bottomleft",
-        legend = "Hydrographie",
-        lty = 1,
-        seg.len = 1,
-        col = "royalblue",
-        cex = 0.8,
-        bty = "n",
-        inset = c(0.18, 0)
-    )
-    
-    mf_layout(title = varTitle,
-              credits = "Inrae")
+    plot_seine_map(r, breaks, title)
+
 }
-Qobs <- cbind(DatesR = BasinsObs$DatesR, BasinsObs$Q)
-Qsim <- cbind(DatesR = OM_with[[1]]$DatesR, as.data.frame(do.call(cbind, lapply(OM_with, "[[", "Qsim"))))
-# Qobs should have same rows as Qsim
-Qobs <- Qobs[Qobs$DatesR %in% Qsim$DatesR, ]
-# Qsim should have same columns as Qobs
-Qsim <- Qsim[, names(Qobs)]
-# Qsim should have same gaps as Qobs
-Qsim[is.na(Qobs)] <- NA
 ```
 
 
+
 ### Ratio débit moyen simulé et observé
 
 ```{r, fig.cap="Carte des R-QA des Bassins versants intermédiaires"}
-getMapRatio(Qsim, Qobs, calcQA, "R-QA")
+plotRatio(Qsim, Qobs, calcQA, "R-QA")
 ```
 
 ### VCN30 de période de retour 2 ans
 
 ```{r, fig.cap="Carte des R-VCN30-2 des Bassins versants intermédiaires"}
-calcVCN30Y2 <- function(x) {
-    calcVCNn(x, 30, 2)
-}
-getMapRatio(Qsim, Qobs, calcVCN30Y2, "R-VCN30-2")
+plotRatio(Qsim, Qobs, calcVCN30_2, "R-VCN30-2")
 ```
 
 ### QMNA5&nbsp;: Débit moyen mensuel minimum annuel de période de retour 5 ans
 
 ```{r, fig.cap="Carte des R-QMNA5 des Bassins versants intermédiaires"}
-calcQMNA5 <- function(x) {
-    calcQMNAn(x, 5)
-}
-getMapRatio(Qsim, Qobs, calcQMNA5, "R-QMNA5")
+plotRatio(Qsim, Qobs, calcQMNA5, "R-QMNA5")
 ```
 
-### QJXA10&nbsp;: Débit moyen journalier maximum annuel de période de retour 10 ans
+### QJXA10&nbsp;: Débit moyen journalier maximum annuel de période de retour 10 ans☺
 
-```{r, fig.cap="Carte des R-QMNA5 des Bassins versants intermédiaires"}
-calcQJXA10 <- function(x) {
-    calcQJXAn(x, 10)
+```{r, fig.cap="Carte des R-QJXA10 des Bassins versants intermédiaires"}
+plotRatio(Qsim, Qobs, calcQJXA10, "R-QJXA10")
+```
+
+## Débits caractéristiques
+
+Nous calculons ici l'ensemble des indicateurs hydrologiques utilisés dans les fiches Explore 2070:
+
+```{r}
+indicators <- c("QA", 
+                "VCN10_2", "VCN30_2", "QMNA2",
+                "VCN10_5", "VCN30_5", "QMNA5",
+                "VCN10_10", "VCN30_10", "QMNA10",
+                "QJXA2", "QJXA10", "QJXA20")
+
+calcQcaract <- function(indicator, Qobs, Qsim) {
+    FUN <- match.fun(paste0("calc", indicator))
+    QcObs <- FUN(Qobs)
+    QcSim <- FUN(Qsim)
+    m <- cbind(QcObs, QcSim, QcSim / QcObs)
+    colnames(m) <- c(paste(indicator, c("obs", "sim"), sep = "."), 
+                     paste0("R-", indicator))
+    return(m)
 }
-getMapRatio(Qsim, Qobs, calcQJXA10, "R-QJXA10")
+l <- lapply(indicators, calcQcaract, Qobs = Qobs, Qsim = Qsim)
+mQc <- do.call(cbind, l)
+dfQc <- cbind(data.frame(Id = rownames(mQc)), mQc)
+write.table(dfQc, 
+          file = getDataPath(cfg$calibration$path, "Q_indicators.tsv"), 
+          sep = "\t", 
+          quote = FALSE,
+          row.names = FALSE)
 ```
+
+Ces indicateurs sont téléchargeable à l'adresse&nbsp;:
+https://owncloud.dorch.fr/index.php/s/pCPZvY4lk6xGC8m/download?path=%2F02-calibration&files=Q_indicators.tsv

man/calcQJXAn.Rd

@@ -5,7 +5,10 @@
 \alias{calcQJXAn.QJXA}
 \alias{calcQJXAn.data.frame}
 \alias{calcQJXAn.GumbelParams}
-\title{Compute the flow for a given Gumbel law}
+\alias{calcQJXA2}
+\alias{calcQJXA10}
+\alias{calcQJXA20}
+\title{Yearly maximum daily flow for a return period}
 \usage{
 calcQJXAn(x, return_period)
 
@@ -14,8 +17,16 @@ calcQJXAn(x, return_period)
 \method{calcQJXAn}{data.frame}(x, return_period)
 
 \method{calcQJXAn}{GumbelParams}(x, return_period)
+
+calcQJXA2(x)
+
+calcQJXA10(x)
+
+calcQJXA20(x)
 }
 \arguments{
+\item{x}{\link{data.frame} of flows (the first column should be a \link{POSIXct}), or a \code{QJXA} object (See \link{calcQJXA})}
+
 \item{return_period}{\link{numeric} Return period. Its unit is given by the aggregation time step used for the Gumbel law.}
 
 \item{gumbel}{A \link{list} with \code{a} the location parameter and \code{b} the scale parameter}
@@ -24,7 +35,9 @@ calcQJXAn(x, return_period)
 A \link{numeric} with the calculated flow
 }
 \description{
-Compute the flow for a given Gumbel law
+The maximum daily flow is calculated from an adjustement on a Gumbel law.
+
+Functions named calcQJXA_[return_period] are also defined for commun indicators (2, 10, 20 years).
 }
 \examples{
 #! load data

man/calcQMNAn.Rd

@@ -4,6 +4,9 @@
 \alias{calcQMNAn}
 \alias{calcQMNAn.QMNA}
 \alias{calcQMNAn.data.frame}
+\alias{calcQMNA2}
+\alias{calcQMNA5}
+\alias{calcQMNA10}
 \title{QMNA for a given return period}
 \usage{
 calcQMNAn(x, return_period)
@@ -11,9 +14,15 @@ calcQMNAn(x, return_period)
 \method{calcQMNAn}{QMNA}(x, return_period)
 
 \method{calcQMNAn}{data.frame}(x, return_period)
+
+calcQMNA2(x)
+
+calcQMNA5(x)
+
+calcQMNA10(x)
 }
 \arguments{
-\item{x}{\link{numeric} vector of mean monthly flows}
+\item{x}{\link{data.frame} of flows (the first column should be a \link{POSIXct}), or a \code{QMNA} object (See \link{calcQMNA})}
 
 \item{return_period}{\link{numeric} number of years}
 }
@@ -23,6 +32,9 @@ calcQMNAn(x, return_period)
 \description{
 Calculation of QMNA according to a return period. An adjustment to a Galton distribution (log-normal distribution) is realised.
 }
+\details{
+Functions named calcQJMNA_[return_period] are also defined for commun indicators (2, 5, 10 years).
+}
 \examples{
 #! load data
 data(L0123001, package = "airGR")

man/calcVCNn.Rd

@@ -3,11 +3,29 @@
 \name{calcVCNn}
 \alias{calcVCNn}
 \alias{calcVCNn.data.frame}
+\alias{calcVCN10_2}
+\alias{calcVCN10_5}
+\alias{calcVCN10_10}
+\alias{calcVCN30_2}
+\alias{calcVCN30_5}
+\alias{calcVCN30_10}
 \title{VCNn for a given return period}
 \usage{
 calcVCNn(x, k, return_period, ...)
 
 \method{calcVCNn}{data.frame}(x, k, return_period, ...)
+
+calcVCN10_2(x)
+
+calcVCN10_5(x)
+
+calcVCN10_10(x)
+
+calcVCN30_2(x)
+
+calcVCN30_5(x)
+
+calcVCN30_10(x)
 }
 \arguments{
 \item{x}{\link{numeric} vector of VCN of one year}
@@ -19,6 +37,8 @@ calcVCNn(x, k, return_period, ...)
 }
 \description{
 Calculation of VCN according to a return period. An adjustment to a Galton distribution (log-normal distribution) is realised.
+
+Functions named calcVCN[k]_[return_period] are also defined for commun indicators.
 }
 \examples{
 #! load data

man/plot_seine_map.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_seine_map.R
+\name{plot_seine_map}
+\alias{plot_seine_map}
+\title{Plot choroplethe map of the Seine River Basin at BVI scale}
+\usage{
+plot_seine_map(r, breaks, title, credits = "IGN, Inrae", ...)
+}
+\arguments{
+\item{r}{named \link{numeric} \link{vector}, the indicator to display with the names corresponding to gauging station IDs}
+
+\item{breaks}{see the parameter in \link[mapsf:mf_map]{mapsf::mf_map}}
+
+\item{title}{see the parameter in \link[mapsf:mf_layout]{mapsf::mf_layout}}
+
+\item{credits}{see the parameter in \link[mapsf:mf_layout]{mapsf::mf_layout}}
+
+\item{...}{parameters sent to \link[mapsf:mf_map]{mapsf::mf_map}}
+}
+\value{
+
+}
+\description{
+Plot choroplethe map of the Seine River Basin at BVI scale
+}