Nothing
In airGRteaching: Teaching Hydrological Modelling with the GR Rainfall-Runoff Models ('Shiny' Interface Included)
R_bib <- toBibtex(citation())
R_bib <- gsub("@Manual\\{", "@Manual{Rsoftware_man", R_bib)
airGR_bib <- toBibtex(citation("airGR"))
airGR_bib <- gsub("@Article\\{", "@Article{airGR_art", airGR_bib)
airGR_bib <- gsub("@Manual\\{", "@Manual{airGR_man", airGR_bib)
airGRteaching_bib <- toBibtex(citation("airGRteaching"))
airGRteaching_bib <- gsub("@Manual\\{", "@Manual{airGRteaching_man", airGRteaching_bib)
airGRteaching_bib <- gsub("@InProceedings\\{", "@InProceedings{airGRteaching_pcd", airGRteaching_bib)
airGRdatasets_bib <- toBibtex(citation("airGRdatasets"))
airGRdatasets_bib <- gsub("@Manual\\{", "@Manual{airGRdatasets_man", airGRdatasets_bib)
options(encoding = "UTF-8")
writeLines(text = c(R_bib, airGR_bib, airGRteaching_bib, airGRdatasets_bib), con = "airGR_galaxy.bib")
options(encoding = "native.enc")
formatGR <- '<strong><font color="#0BA6AA">%s</font></strong>'
GR <- sprintf(formatGR, "GR")
airGR <- sprintf(formatGR, "airGR")
airGRteaching <- sprintf(formatGR, "airGRteaching")
knitr::opts_chunk$set(echo = TRUE, fig.path = "figure/")
library(airGRteaching)
Sys.setlocale("LC_TIME", "English")
colorize <- function(x, color) {
if (knitr::is_latex_output()) {
sprintf("\\textcolor{%s}{%s}", color, x)
} else if (knitr::is_html_output()) {
sprintf("<span style='color: %s;'>%s</span>", color, x)
} else {
x
}
}
# Catchment data loading
library(airGRdatasets)
data("B222001001", package = "airGRdatasets")
# Ctachment area
area <- B222001001$Meta$Area
# Observed daily time series
ts_init <- B222001001$TS
# Add date information to the time series
ts_init$Year <- format(ts_init$Date, format = "%Y")
ts_init$MonthDay <- format(ts_init$Date, format = "%m-%d")
# Display of the 1st time steps of the time series
head(ts_init)
name_sta <- gsub("the ", "", B222001001$Meta$Name)
name_riv <- gsub("(the )(.*)( at.*)", "\\2", name_sta)
# Calibration period
per_cal_hist <- c("2000-09-01", "2018-08-31")
# Forecasting period
per_fcst <- c("2018-09-01", "2018-10-20")
# Warm-up period to simulate on the historical period
per_wup_hist <- c("1999-01-01", "2000-08-31")
# Warm-up period to simulate on the forcasting period
per_wup_fcst <- c(per_wup_hist[1], per_cal_hist[2])
# Forecasting dates
dates_fcst <- seq(from = as.POSIXct(per_fcst[1], tz = "UTC", format = "%Y-%m-%d"),
to = as.POSIXct(per_fcst[2], tz = "UTC", format = "%Y-%m-%d"),
by = "1 day")
head(dates_fcst)
# Formatting of forecast dates (Month-Day)
month_day_fcst <- format(dates_fcst, "%m-%d")
head(month_day_fcst)
# Set values of the last winter to missing data
ts_plot <- ts_init
isd_wint <- ts_plot$Date >= as.POSIXct(per_fcst[1], tz = "UTC", format = "%Y-%m-%d")
ts_plot[isd_wint, c("Ptot", "Temp", "Evap", "Qls", "Qmmd")] <- NA
# Display of the last time steps of the time series
tail(ts_plot)
per_wup_pct <- as.POSIXct(per_wup_hist, tz = "UTC", format = "%Y-%m-%d")
per_wup_for <- format(per_wup_pct, "%e %B %Y")
per_cal_pct <- as.POSIXct(per_cal_hist, tz = "UTC", format = "%Y-%m-%d")
per_cal_for <- format(per_cal_pct, "%e %B %Y")
per_fcst_pct <- as.POSIXct(per_fcst, tz = "UTC", format = "%Y-%m-%d")
per_fcst_for <- format(per_fcst_pct, "%e %B %Y")
per_fcst_for_no_y <- format(per_fcst_pct, "%e %B")
# Select the time series over the observed period (no "future" dates)
ts_hist <- ts_plot[ts_plot$Date < dates_fcst[1], ]
# Display of the last time steps of the time series
tail(ts_hist)
# Select the time series after the observed period (only "future" dates)
ts_fcst <- ts_plot[ts_plot$Date >= dates_fcst[1] & ts_plot$Date <= dates_fcst[length(dates_fcst)], ]
# Display of the 1st time steps of the time series
head(ts_fcst)
# Information for the vignette
year_per <- as.integer(format(range(ts_hist$Date), format = "%Y"))
year_min <- year_per[1L]
year_max <- year_per[2L]
n_year <- year_max - year_min + 1
Objective
Context
End of summer r year_max, the river r name_riv has reached its lowest level for several years and the rain is still awaited for on this catchment characterized by many water uses in summer. In charge of the monitoring of the streamflows in period of low-flow, you have been called by the local officials who are panicking about the drought situation and ask you to produce an estimate of the streamflows of the river for the next weeks.
Unfortunately, you have not received the weather forecast from the regional forecasting service for several weeks... You do not have any forecast of the coming precipitation, but you still have to quantify the possible streamflows on the catchment for the next weeks.
The exercise consists in using the available hydro-climatic data (at daily time step) on r name_sta (r round(area) km²), and a rainfall-runoff model to forecast the streamflows of the coming weeks by using (1) the last observed streamflow on the considered catchment, (2) the observed precipitation climatology and (3) the observed streamflow climatology. This information will provide ranges of potential streamflows for the coming weeks (see following figure).
This work will be completed in four steps:
- Statistical analysis of the streamflow climatology to obtain an order of magnitude of the streamflows historically observed during the study period, i.e. the daily streamflow regime.
- Calibration of the hydrological model on the historical period considering a criterion centered on low flows.
- Pessimistic scenario of zero precipitation: rainfall-runoff simulation over the summer studied from the last observed streamflow considering no future precipitation and a series of potential evapotranspiration (PE) constituted from the interannual average PE.
- Non-zero future precipitation scenarios: rainfall-runoff simulations for the studied summer based on the last observed discharge by hindcasting the precipitation climatology, i.e. the precipitation observed during the other summers. There will be as many simulations as there are historical years.
# Indices of the plotting period
ind_xlim <- c(nrow(ts_plot)-300, nrow(ts_plot))
# Plotting hyrdograph
plot(x = ts_plot$Date, y = ts_plot$Qmmd,
xlab = "time [days]", ylab = "flow [mm/day]",
main = paste("Year", year_max),
type = "l", lwd = 1,
log = "y",
xlim = ts_plot$Date[ind_xlim])
<<v02_fig_pattern>>
# Add a future shading period
rect(xleft = per_fcst_pct[1], ybottom = 1e-6,
xright = per_fcst_pct[2], ytop = 1e+6,
col = "lightgrey", border = NA)
box()
# Add quetion mark in the future period
text(x = mean(per_fcst_pct),
y = quantile(ts_plot$Qmmd, probs = 0.5, na.rm = TRUE),
labels = "?", cex = 8, font = 2, col = "orangered")
Instructions
This section aims to make explicit some of the expected tasks and to describe the hydrological model calibration and simulation conditions (parameter calibration period, reservoir initialization periods, calibration criteria, etc.).
Analysis of streamflow climatology
In a hydrological forecasting context, "flow climatology" refers to the streamflow regime, i.e. the series of interannual mean streamflows. If you want to forecast the streamflow of the day to come (for example on r per_fcst_for_no_y[1]), the simplest forecast from a technical point of view consists in calculating the average of the streamflows observed on this day during the past years. This average can also be framed thanks to the calculations of quantiles of these same historical streamflows (for example the quantiles 10 and 90 %), which makes it possible to have a first "probabilistic" estimate of the future streamflows, without carrying out rainfall-runoff modelling and without using a precipitation forecast.
In this exercise, you analyze all the streamflows observed during the forecast period (from r per_fcst_for[1] to r per_fcst_for[2]) during the previous years. r n_year years are available for this analysis, which allows to constitute r n_year series of historical daily streamflows for the forecast period. Once grouped in a single table, these streamflows can be summarized by calculating quantiles such as 10, 25, 50, 75 and 90 % quantiles for each day considered. These quantiles will then be used as a first estimate of future streamflows.
Rainfall-runoff model
You will use the GR6J model [@pushpalatha_downward_2011]. It is a lumped rainfall-runoff model, operating on a daily time step and having six parameters. It requires continuous time series of daily precipitation and potential evapotranspiration (PE) as input.
This model is easily usable thanks to the r airGRteaching package [@airGRteaching_man; @airGRteaching_pcd], developped for the R software by the Catchment Hydrology research group of the HYCAR reaserch unit (INRAE, France).
The time series of observed precipitation, PE and streamflow can be easily formatted using the PrepGR() function. A rainfall-runoff simulation can be performed with the SimGR() function, and a parameter calibration using the CalGR().
Calibration (and warm-up) period
Initialization is particularly important in a context of hydrological forecasting by modelling: this warm-up consists in estimating the saturation state of the basin (i.e. the internal states of the model) which will then be used as initial conditions for all the tested scenarios. These scenarios are all based on the same state at the date of launching the forecasts.
The period to be considered to calibrate the parameters of the model on r name_sta catchment starts on r per_cal_for[1] and ends on r per_cal_for[2]. A warm-up period of r as.numeric(round(diff(per_wup_pct)/365.2425*12)) months should also be considered from r per_wup_for[1] to r per_wup_for[2].
Calibration criterion
The calibration criterion to be considered in this exercise is the Nash and Sutcliffe criterion [@nash_river_1970] calculated on the natural logarithms of the streamflows, noted $NSE_{ln}$ hereafter (see following equation). This criterion is widely used in hydrological modelling.
The NSE criterion, bounded between $-\infty$ and $1$, allows to quantify the performance of a model in a relative way, by comparing a series of simulated streamflows with a so-called "naive" model, here the average of the observed streamflows (i.e. a series of streamflows constituted in each time step by the average of the observed streamflows). Thus, a NSE value equal to 1 indicates a perfect agreement between the series of observed and simulated streamflows (which is never the case), whereas a NSE value lower than 0 means that the simulation considered is less efficient than the simulation of the "naive" model. The calculation of $NSE_{ln}$ is detailed in the following equation, where $Q_{obs,t}$ is the observed streamflow at time step $t$, $Q_{sim,t}$ is the simulated streamflow at time step $t$, $\overline{Q_{obs}}$ is the average of the observed streamflows, $\epsilon$ is a constant defined below, and $n$ is the number of observations:
\begin{equation}
NSE_{ln} = 1-\frac{\sum_{t=1}^{n}(ln(Q_{obs,t}+\epsilon)-ln(Q_{sim,t}+\epsilon))^{2}}{\sum_{t=1}^{n}(ln(Q_{obs,t}+\epsilon)-\overline{ln(Q_{obs}+\epsilon)})^{2}}
\end{equation}
The calculation of the $\epsilon$ value is detailed in the following equation. The addition of this $\epsilon$ is necessary when zero streamflows are present in the observed streamflow series.
\begin{equation}
\epsilon=\frac{\overline{Q_{obs}}}{100}
\end{equation}
The logarithmic transformation allows to give more weight to the lowest streamflows, and thus to limit the errors made on the simulation of low flows [@oudin_dynamic_2006]. Considering a series constituted by the inverse of the streamflows during the calibration is also recommended to obtain a good performance during the simulation of very low flows [@pushpalatha_review_2012].
The elements necessary for the calculation of the calibration criterion have to be set as arguments of the CalGR() function, which allows to specify if the criterion must be calculated on "transformed" streamflows with the transfo argument (here equal to "log").
Automatic calibration of the model parameters
Automatic parameter estimation aims at using a search algorithm in the parameter space. This algorithm will automatically generate sets of parameters, test them, and generate others according to the performance of those already tested, until converging to an optimum. The algorithm developed by @michel_hydrologie_1991 will be used in this exercise.
Simulation period
The simulation period, which can also be called the forecast period in this exercise, runs from r per_fcst_for[1] to r per_fcst_for[2]. The set of time steps preceding this calibration period can be used as the warm-up period.
Pessimistic scenario of zero precipitation
Having a rainfall-runoff model calibrated on the studied catchment allows to use future precipitation scenarios and to transform them into streamflow scenarios. The simplest scenario to test (and the most pessimistic one) is the so-called "zero precipitation" scenario: no precipitation will be observed over the time steps of the forecast period. This will be the lower bound of the rainfall-runoff modelling forecast.
To implement this scenario in the r airGRteaching environment, it is necessary to create a "dummy" data table containing a time series of precipitation equal to 0 over the forecast period. For PE, a realistic assumption is to use the interannual regime of this variable. Thus, for each forecast day, the interannual average value of observed PE will be used: the PE of r per_fcst_for_no_y[1] for this scenario will be equal to the average of observed PE values for all past r per_fcst_for_no_y[1].
Non-zero future precipitationscenario
This last step aims at performing several rainfall-runoff simulations considering non-zero future precipitation scenarios. These precipitation scenarios will be constructed from historical precipitation data. r n_year years of past precipitation are available for this analysis, which allows to build r n_year scenarios of historical daily precipitation and PE for the forecast period. A rainfall-runoff simulation will be performed for each of these past years, starting the simulation with the same warm-up period. Once grouped in a single table, the simulated streamflows can be summarized by calculating quantiles such as 10, 25, 50, 75 and 90 % quantiles for each day considered. These quantiles will then be used as an indication of possible future streamflows.
Data available
The data set available to the rainfall-runoff modelling consists of:
- a time series of daily total precipitation (liquid + solid) [mm/day] (
Ptot);
- a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/day] (
Evap);
- a time series of daily streamflows expressed as a specific discharge [mm/day] (
Qmmd).
Command lines for the production of simulations
Loading and formatting of data
The following command lines are used to prepare the data required to calibrate the rainfall-runoff model, and to define the working periods (initialization period, calibration period and forecast period).
Initial time series
Ploted time series
Observed time series
Forecast time series
Analysis of streamflow climatology
The following command lines allow the analysis of the streamflow climatology by calculating the interannual streamflow quantiles over the forecasting period:
# Calculation of the historical streamflow quantiles
ts_qclim_quant <- aggregate(Qmmd ~ MonthDay,
data = ts_hist[ts_hist$MonthDay %in% month_day_fcst, ],
FUN = function(x) {
quantile(x, probs = c(0.10, 0.25, 0.50, 0.75, 0.90))
})
ts_qclim_quant <- as.data.frame(ts_qclim_quant$Qmmd)
colnames(ts_qclim_quant) <- paste0("Q", gsub("\\D", "", colnames(ts_qclim_quant)))
rownames(ts_qclim_quant) <- month_day_fcst
# Display of the 1st calculated quantiles
head(ts_qclim_quant)
The streamflow climatology allows to obtain a first order of magnitude of the possible streamflows of the next days (cf. following figure), but does not allow to account for the last observed streamflows.
<<v02_fig_pattern>>
polygon(x = c(dates_fcst, rev(dates_fcst)),
y = c(ts_qclim_quant$Q10, rev(ts_qclim_quant$Q90)),
col = "lightgrey", border = NA)
polygon(x = c(dates_fcst, rev(dates_fcst)),
y = c(ts_qclim_quant$Q25, rev(ts_qclim_quant$Q75)),
col = "darkgrey", border = NA)
lines(x = dates_fcst, y = ts_qclim_quant$Q50, lwd = 2, col = "grey50")
legend("topright",
legend = c("Qobs", "Qclim (10,25,50,75,90 %)"),
col = c("black", "grey50"),
pch = c(NA, 15), lwd = c(1, 2), pt.cex = 2,
lty = c(1, 1), pt.bg = c(NA, "lightgrey"),
bg = "white")
Data processing for GR6J
The following command lines aim to prepare the available data for use by GR6J, using the PrepGR() function of the r airGRteaching package.
# Adding an epsilon to observed streamflows for criterion calculation
epsilon_qobs <- mean(ts_hist$Qmmd, na.rm = TRUE) / 100
# Data processing for GR6J
prep_hist <- PrepGR(DatesR = ts_hist$Date,
Precip = ts_hist$Ptot,
PotEvap = ts_hist$Evap,
Qobs = ts_hist$Qmmd + epsilon_qobs,
HydroModel = "GR6J",
CemaNeige = FALSE)
GR6J calibration on the historical period
The following command lines are used to calibrate the GR6J model on the historical period.
# Calibration step
cal_hist <- CalGR(PrepGR = prep_hist,
CalCrit = "NSE",
WupPer = per_wup_hist,
CalPer = per_cal_hist,
transfo = "log",
verbose = TRUE)
# Get parameter and criterion values
param_cal_hist <- GetParam(cal_hist)
crit_cal_hist <- GetCrit(cal_hist)
# Graphical assessment of the calibration performance
plot(cal_hist)
The six parameters and the value of the calibration criterion ($NSE_{ln}$) obtained after the automatic calibration procedure are:
- X1 =
r round(param_cal_hist[1], digits = 0) [mm]
- X2 =
r round(param_cal_hist[2], digits = 3) [mm/day]
- X3 =
r round(param_cal_hist[3], digits = 0) [mm]
- X4 =
r round(param_cal_hist[4], digits = 3) [day]
- X5 =
r round(param_cal_hist[5], digits = 3) [-]
- X6 =
r round(param_cal_hist[6], digits = 2) [mm]
- $NSE_{ln}$ =
r round(crit_cal_hist, digits = 3) [-]
The performance obtained in calibration is very good, since the criterion $NSE_{ln}$ is equal to r round(crit_cal_hist, digits = 3).
The following command lines allow to store in the same table the observed streamflows and the simulated streamflows with the set of parameters obtained by automatic calibration, in order to compare them.
# Combination of observed and simulated streamflow time series on the calibration period
ts_cal_hist <- as.data.frame(cal_hist)
head(ts_cal_hist)
The following figure represents the observed and simulated streamflow series over the end of the calibration period.
# Plot framework
<<v02_fig_pattern>>
# sim
lines(x = ts_cal_hist$Date, y = ts_cal_hist$Qsim, col = "orangered")
legend("topright",
legend = c("Qobs", "Qsim"),
lty = 1, col = c("black", "orangered"),
bg = "white")
Pessimistic zero precipitation scenario
The following command lines are used to create a "dummy" time series that combines observed data (for model warm-up) and data generated for the "zero precipitation" scenario. This time series will then be used as input to the GR6J model to "transform" this meteorological scenario into a hydrological scenario.
# Duplicate the time series with future dates to fill the forecast period
ts_fcst_p0 <- ts_fcst
# Setting zero precipitation for the forecast period
ts_fcst_p0$Ptot <- 0
# Addition of interannual average PE
ts_eclim <- aggregate(Evap ~ MonthDay,
data = ts_hist[ts_hist$MonthDay %in% month_day_fcst, ],
FUN = mean)
ts_fcst_p0$Evap <- ts_eclim$Evap
# Combine historical and forecasting time series
ts_scen_p0 <- rbind(ts_hist, ts_fcst_p0)
# Display of the last lines
tail(ts_scen_p0)
The following command lines allow us to format the previously created time series as inputs to the GR6J model and to perform a rainfall-runoff simulation.
# Data processing for GR6J
prep_scen_p0 <- PrepGR(DatesR = ts_scen_p0$Date,
Precip = ts_scen_p0$Ptot,
PotEvap = ts_scen_p0$Evap,
Qobs = NULL,
HydroModel = "GR6J",
CemaNeige = FALSE)
# Hydrological forecast
sim_scen_p0 <- SimGR(PrepGR = prep_scen_p0,
WupPer = per_wup_fcst,
SimPer = per_fcst,
Param = param_cal_hist,
verbose = FALSE)
# Simulated streamflow time series
ts_sim_scen_p0 <- as.data.frame(sim_scen_p0)
The following figure represents the observed and simulated streamflow series over the end of the calibration period, and the result of the simulation of the "zero precipitation" scenario. The slow decay of the simulated streamflows, generated by the slow emptying of the GR6J reservoirs, is observed.
# Plot framework
<<v02_fig_step_cal_hist>>
# Zero precipitation
lines(x = ts_sim_scen_p0$Dates, y = ts_sim_scen_p0$Qsim, lwd = 2, col = "blue")
# Legend
legend("topright",
legend = c("Qobs", "Qsim", "QscenP0"),
lwd = 2, lty = 1, col = c("black", "orangered", "blue"),
bg = "white")
It is important to note that the simulations performed so far do not take into account the last observed streamflow. This can generate a discontinuity in the sequence of streamflows, between the observations available until the time of the forecast, and the simulations produced from this time. This "error" can be corrected in several ways, by "assimilating" the last observed streamflow. The simplest correction consists in calculating a simple ratio between the last observed streamflow and the corresponding simulated streamflow, and to use this ratio to correct all the following simulated streamflows. The following command lines allow to calculate such a ratio, which will be used thereafter to correct the streamflows simulated by GR6J, by dividing them by this ratio (e.g. following figure).
# Correction (~ assimilation)
corr_qsim <- ts_sim_scen_p0$Qsim[1] / ts_hist$Qmmd[nrow(ts_hist)]
# Ratio display
corr_qsim
# Plot framework
<<v02_fig_pattern>>
# Zero precipitation
lines(x = ts_sim_scen_p0$Date, y = ts_sim_scen_p0$Qsim / corr_qsim, lwd = 2, col = "cornflowerblue")
# Legend
legend("topright",
legend = c("Qobs", "QscenP0corr"),
lwd = 2, lty = 1, col = c("black", "cornflowerblue"),
bg = "white")
Non-zero future precipitation scenarios
The following command lines allow us to perform rainfall-runoff simulations for the studied summer using the last observed streamflow by "replaying" the precipitation climatology, i.e. the precipitation observed during other summers.
# Historical years
year_hist <- unique(ts_hist$Year)
year_hist <- setdiff(year_hist, unique(ts_fcst$Year))
year_hist
# Loop on the historical years
ts_qscen_year <- sapply(year_hist, FUN = function(i_year) {
i_ts_hist <- ts_hist[ts_hist$Year == i_year & ts_hist$MonthDay %in% month_day_fcst, ]
i_ts_hist$Date <- dates_fcst
# Combine historical and forecasting (based on precipitation climatology) time series
i_ts_scen_py <- rbind(ts_hist, i_ts_hist)
# Data processing for GR6J
prep_scen_py <- PrepGR(DatesR = i_ts_scen_py$Date,
Precip = i_ts_scen_py$Ptot,
PotEvap = i_ts_scen_py$Evap,
HydroModel = "GR6J",
CemaNeige = FALSE)
# Hydrological forecast
sim_scen_py <- SimGR(PrepGR = prep_scen_py,
WupPer = per_wup_fcst,
SimPer = per_fcst,
Param = param_cal_hist,
verbose = FALSE)
# Correction of the simulated streamflow
sim_scen_py$OutputsModel$Qsim / corr_qsim
})
# Compute future streamflow quantiles (based on precipitation climatology)
ts_qscen_quant <- t(apply(ts_qscen_year, MARGIN = 1, FUN = function(x) {
quantile(x, probs = c(0.10, 0.25, 0.50, 0.75, 0.90))
}))
ts_qscen_quant <- as.data.frame(ts_qscen_quant)
colnames(ts_qscen_quant) <- paste0("Q", gsub("\\D", "", colnames(ts_qscen_quant)))
rownames(ts_qscen_quant) <- month_day_fcst
# Display of the first calculated quantiles
head(ts_qscen_quant)
Streamflow forecasts based on the precipitation climatology suggest that the next few days will be dry, and that an increase in the streamflow of the river r name_riv is not expected until a few days.
# Plot framework
<<v02_fig_pattern>>
polygon(x = c(dates_fcst, rev(dates_fcst)),
y = c(ts_qscen_quant$Q10, rev(ts_qscen_quant$Q90)),
col = adjustcolor("green2", alpha.f = 0.3), border = NA)
polygon(x = c(dates_fcst, rev(dates_fcst)),
y = c(ts_qscen_quant$Q25, rev(ts_qscen_quant$Q75)),
col = adjustcolor("green2", alpha.f = 0.5), border = NA)
lines(x = dates_fcst, y = ts_qscen_quant$Q50, lwd = 2, col = "green4")
legend("topright",
legend = c("Qobs", "QscenPy (10,25,50,75,90 %)"),
col = c("black", "green2"),
pch = c(NA, 15), lwd = c(1, 2), pt.cex = 2,
lty = c(1, 1), pt.bg = c(NA, "green2"),
bg = "white")
The following figure allows to compare all the results on the same graph. It displays observed and simulated streamflows over the end of the calibration period, the streamflow climatology over the forecast period, the result of the simulation of the "zero precipitation" scenario, and the streamflows based on the precipitation climatology.
# Plot framework
<<v02_fig_pattern>>
par(new = TRUE)
<<v02_fig_qclim>>
par(new = TRUE)
<<v02_fig_cal_hist>>
par(new = TRUE)
<<v02_fig_qscen_py>>
par(new = TRUE)
<<v02_fig_qscen_p0corr>>
legend("topright",
legend = c("Qobs",
"Qsim",
"QscenP0corr",
"Qclim (10,25,50,75,90 %)",
"QscenPy (10,25,50,75,90 %)"),
col = c("black", "orangered", "cornflowerblue", "grey50", "green2"),
pch = c(NA, NA, NA, 15, 15),
lwd = c(1, 1, 2, 2, 2),
pt.cex = 2,
lty = c(1, 1, 1, 1, 1),
pt.bg = c(NA, NA, NA, "grey50", "green2"),
bg = "white")
References
Try the airGRteaching package in your browser
Any scripts or data that you put into this service are public.
airGRteaching documentation built on July 26, 2023, 5:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
R_bib <- toBibtex(citation()) R_bib <- gsub("@Manual\\{", "@Manual{Rsoftware_man", R_bib) airGR_bib <- toBibtex(citation("airGR")) airGR_bib <- gsub("@Article\\{", "@Article{airGR_art", airGR_bib) airGR_bib <- gsub("@Manual\\{", "@Manual{airGR_man", airGR_bib) airGRteaching_bib <- toBibtex(citation("airGRteaching")) airGRteaching_bib <- gsub("@Manual\\{", "@Manual{airGRteaching_man", airGRteaching_bib) airGRteaching_bib <- gsub("@InProceedings\\{", "@InProceedings{airGRteaching_pcd", airGRteaching_bib) airGRdatasets_bib <- toBibtex(citation("airGRdatasets")) airGRdatasets_bib <- gsub("@Manual\\{", "@Manual{airGRdatasets_man", airGRdatasets_bib) options(encoding = "UTF-8") writeLines(text = c(R_bib, airGR_bib, airGRteaching_bib, airGRdatasets_bib), con = "airGR_galaxy.bib") options(encoding = "native.enc")
formatGR <- '<strong><font color="#0BA6AA">%s</font></strong>' GR <- sprintf(formatGR, "GR") airGR <- sprintf(formatGR, "airGR") airGRteaching <- sprintf(formatGR, "airGRteaching")
knitr::opts_chunk$set(echo = TRUE, fig.path = "figure/") library(airGRteaching) Sys.setlocale("LC_TIME", "English") colorize <- function(x, color) { if (knitr::is_latex_output()) { sprintf("\\textcolor{%s}{%s}", color, x) } else if (knitr::is_html_output()) { sprintf("<span style='color: %s;'>%s</span>", color, x) } else { x } }
# Catchment data loading library(airGRdatasets) data("B222001001", package = "airGRdatasets") # Ctachment area area <- B222001001$Meta$Area # Observed daily time series ts_init <- B222001001$TS # Add date information to the time series ts_init$Year <- format(ts_init$Date, format = "%Y") ts_init$MonthDay <- format(ts_init$Date, format = "%m-%d") # Display of the 1st time steps of the time series head(ts_init)
name_sta <- gsub("the ", "", B222001001$Meta$Name) name_riv <- gsub("(the )(.*)( at.*)", "\\2", name_sta)
# Calibration period per_cal_hist <- c("2000-09-01", "2018-08-31") # Forecasting period per_fcst <- c("2018-09-01", "2018-10-20") # Warm-up period to simulate on the historical period per_wup_hist <- c("1999-01-01", "2000-08-31") # Warm-up period to simulate on the forcasting period per_wup_fcst <- c(per_wup_hist[1], per_cal_hist[2]) # Forecasting dates dates_fcst <- seq(from = as.POSIXct(per_fcst[1], tz = "UTC", format = "%Y-%m-%d"), to = as.POSIXct(per_fcst[2], tz = "UTC", format = "%Y-%m-%d"), by = "1 day") head(dates_fcst) # Formatting of forecast dates (Month-Day) month_day_fcst <- format(dates_fcst, "%m-%d") head(month_day_fcst)
# Set values of the last winter to missing data ts_plot <- ts_init isd_wint <- ts_plot$Date >= as.POSIXct(per_fcst[1], tz = "UTC", format = "%Y-%m-%d") ts_plot[isd_wint, c("Ptot", "Temp", "Evap", "Qls", "Qmmd")] <- NA # Display of the last time steps of the time series tail(ts_plot)
per_wup_pct <- as.POSIXct(per_wup_hist, tz = "UTC", format = "%Y-%m-%d") per_wup_for <- format(per_wup_pct, "%e %B %Y") per_cal_pct <- as.POSIXct(per_cal_hist, tz = "UTC", format = "%Y-%m-%d") per_cal_for <- format(per_cal_pct, "%e %B %Y") per_fcst_pct <- as.POSIXct(per_fcst, tz = "UTC", format = "%Y-%m-%d") per_fcst_for <- format(per_fcst_pct, "%e %B %Y") per_fcst_for_no_y <- format(per_fcst_pct, "%e %B")
# Select the time series over the observed period (no "future" dates) ts_hist <- ts_plot[ts_plot$Date < dates_fcst[1], ] # Display of the last time steps of the time series tail(ts_hist)
# Select the time series after the observed period (only "future" dates) ts_fcst <- ts_plot[ts_plot$Date >= dates_fcst[1] & ts_plot$Date <= dates_fcst[length(dates_fcst)], ] # Display of the 1st time steps of the time series head(ts_fcst)
# Information for the vignette year_per <- as.integer(format(range(ts_hist$Date), format = "%Y")) year_min <- year_per[1L] year_max <- year_per[2L] n_year <- year_max - year_min + 1
Objective
Context
End of summer r year_max, the river r name_riv has reached its lowest level for several years and the rain is still awaited for on this catchment characterized by many water uses in summer. In charge of the monitoring of the streamflows in period of low-flow, you have been called by the local officials who are panicking about the drought situation and ask you to produce an estimate of the streamflows of the river for the next weeks.
Unfortunately, you have not received the weather forecast from the regional forecasting service for several weeks... You do not have any forecast of the coming precipitation, but you still have to quantify the possible streamflows on the catchment for the next weeks.
The exercise consists in using the available hydro-climatic data (at daily time step) on r name_sta (r round(area) km²), and a rainfall-runoff model to forecast the streamflows of the coming weeks by using (1) the last observed streamflow on the considered catchment, (2) the observed precipitation climatology and (3) the observed streamflow climatology. This information will provide ranges of potential streamflows for the coming weeks (see following figure).
This work will be completed in four steps:
- Statistical analysis of the streamflow climatology to obtain an order of magnitude of the streamflows historically observed during the study period, i.e. the daily streamflow regime.
- Calibration of the hydrological model on the historical period considering a criterion centered on low flows.
- Pessimistic scenario of zero precipitation: rainfall-runoff simulation over the summer studied from the last observed streamflow considering no future precipitation and a series of potential evapotranspiration (PE) constituted from the interannual average PE.
- Non-zero future precipitation scenarios: rainfall-runoff simulations for the studied summer based on the last observed discharge by hindcasting the precipitation climatology, i.e. the precipitation observed during the other summers. There will be as many simulations as there are historical years.
# Indices of the plotting period ind_xlim <- c(nrow(ts_plot)-300, nrow(ts_plot)) # Plotting hyrdograph plot(x = ts_plot$Date, y = ts_plot$Qmmd, xlab = "time [days]", ylab = "flow [mm/day]", main = paste("Year", year_max), type = "l", lwd = 1, log = "y", xlim = ts_plot$Date[ind_xlim])
<<v02_fig_pattern>> # Add a future shading period rect(xleft = per_fcst_pct[1], ybottom = 1e-6, xright = per_fcst_pct[2], ytop = 1e+6, col = "lightgrey", border = NA) box() # Add quetion mark in the future period text(x = mean(per_fcst_pct), y = quantile(ts_plot$Qmmd, probs = 0.5, na.rm = TRUE), labels = "?", cex = 8, font = 2, col = "orangered")
Instructions
This section aims to make explicit some of the expected tasks and to describe the hydrological model calibration and simulation conditions (parameter calibration period, reservoir initialization periods, calibration criteria, etc.).
Analysis of streamflow climatology
In a hydrological forecasting context, "flow climatology" refers to the streamflow regime, i.e. the series of interannual mean streamflows. If you want to forecast the streamflow of the day to come (for example on r per_fcst_for_no_y[1]), the simplest forecast from a technical point of view consists in calculating the average of the streamflows observed on this day during the past years. This average can also be framed thanks to the calculations of quantiles of these same historical streamflows (for example the quantiles 10 and 90 %), which makes it possible to have a first "probabilistic" estimate of the future streamflows, without carrying out rainfall-runoff modelling and without using a precipitation forecast.
In this exercise, you analyze all the streamflows observed during the forecast period (from r per_fcst_for[1] to r per_fcst_for[2]) during the previous years. r n_year years are available for this analysis, which allows to constitute r n_year series of historical daily streamflows for the forecast period. Once grouped in a single table, these streamflows can be summarized by calculating quantiles such as 10, 25, 50, 75 and 90 % quantiles for each day considered. These quantiles will then be used as a first estimate of future streamflows.
Rainfall-runoff model
You will use the GR6J model [@pushpalatha_downward_2011]. It is a lumped rainfall-runoff model, operating on a daily time step and having six parameters. It requires continuous time series of daily precipitation and potential evapotranspiration (PE) as input.
This model is easily usable thanks to the r airGRteaching package [@airGRteaching_man; @airGRteaching_pcd], developped for the R software by the Catchment Hydrology research group of the HYCAR reaserch unit (INRAE, France).
The time series of observed precipitation, PE and streamflow can be easily formatted using the PrepGR() function. A rainfall-runoff simulation can be performed with the SimGR() function, and a parameter calibration using the CalGR().
Calibration (and warm-up) period
Initialization is particularly important in a context of hydrological forecasting by modelling: this warm-up consists in estimating the saturation state of the basin (i.e. the internal states of the model) which will then be used as initial conditions for all the tested scenarios. These scenarios are all based on the same state at the date of launching the forecasts.
The period to be considered to calibrate the parameters of the model on r name_sta catchment starts on r per_cal_for[1] and ends on r per_cal_for[2]. A warm-up period of r as.numeric(round(diff(per_wup_pct)/365.2425*12)) months should also be considered from r per_wup_for[1] to r per_wup_for[2].
Calibration criterion
The calibration criterion to be considered in this exercise is the Nash and Sutcliffe criterion [@nash_river_1970] calculated on the natural logarithms of the streamflows, noted $NSE_{ln}$ hereafter (see following equation). This criterion is widely used in hydrological modelling.
The NSE criterion, bounded between $-\infty$ and $1$, allows to quantify the performance of a model in a relative way, by comparing a series of simulated streamflows with a so-called "naive" model, here the average of the observed streamflows (i.e. a series of streamflows constituted in each time step by the average of the observed streamflows). Thus, a NSE value equal to 1 indicates a perfect agreement between the series of observed and simulated streamflows (which is never the case), whereas a NSE value lower than 0 means that the simulation considered is less efficient than the simulation of the "naive" model. The calculation of $NSE_{ln}$ is detailed in the following equation, where $Q_{obs,t}$ is the observed streamflow at time step $t$, $Q_{sim,t}$ is the simulated streamflow at time step $t$, $\overline{Q_{obs}}$ is the average of the observed streamflows, $\epsilon$ is a constant defined below, and $n$ is the number of observations:
\begin{equation} NSE_{ln} = 1-\frac{\sum_{t=1}^{n}(ln(Q_{obs,t}+\epsilon)-ln(Q_{sim,t}+\epsilon))^{2}}{\sum_{t=1}^{n}(ln(Q_{obs,t}+\epsilon)-\overline{ln(Q_{obs}+\epsilon)})^{2}} \end{equation}
The calculation of the $\epsilon$ value is detailed in the following equation. The addition of this $\epsilon$ is necessary when zero streamflows are present in the observed streamflow series.
\begin{equation} \epsilon=\frac{\overline{Q_{obs}}}{100} \end{equation}
The logarithmic transformation allows to give more weight to the lowest streamflows, and thus to limit the errors made on the simulation of low flows [@oudin_dynamic_2006]. Considering a series constituted by the inverse of the streamflows during the calibration is also recommended to obtain a good performance during the simulation of very low flows [@pushpalatha_review_2012].
The elements necessary for the calculation of the calibration criterion have to be set as arguments of the CalGR() function, which allows to specify if the criterion must be calculated on "transformed" streamflows with the transfo argument (here equal to "log").
Automatic calibration of the model parameters
Automatic parameter estimation aims at using a search algorithm in the parameter space. This algorithm will automatically generate sets of parameters, test them, and generate others according to the performance of those already tested, until converging to an optimum. The algorithm developed by @michel_hydrologie_1991 will be used in this exercise.
Simulation period
The simulation period, which can also be called the forecast period in this exercise, runs from r per_fcst_for[1] to r per_fcst_for[2]. The set of time steps preceding this calibration period can be used as the warm-up period.
Pessimistic scenario of zero precipitation
Having a rainfall-runoff model calibrated on the studied catchment allows to use future precipitation scenarios and to transform them into streamflow scenarios. The simplest scenario to test (and the most pessimistic one) is the so-called "zero precipitation" scenario: no precipitation will be observed over the time steps of the forecast period. This will be the lower bound of the rainfall-runoff modelling forecast.
To implement this scenario in the r airGRteaching environment, it is necessary to create a "dummy" data table containing a time series of precipitation equal to 0 over the forecast period. For PE, a realistic assumption is to use the interannual regime of this variable. Thus, for each forecast day, the interannual average value of observed PE will be used: the PE of r per_fcst_for_no_y[1] for this scenario will be equal to the average of observed PE values for all past r per_fcst_for_no_y[1].
Non-zero future precipitationscenario
This last step aims at performing several rainfall-runoff simulations considering non-zero future precipitation scenarios. These precipitation scenarios will be constructed from historical precipitation data. r n_year years of past precipitation are available for this analysis, which allows to build r n_year scenarios of historical daily precipitation and PE for the forecast period. A rainfall-runoff simulation will be performed for each of these past years, starting the simulation with the same warm-up period. Once grouped in a single table, the simulated streamflows can be summarized by calculating quantiles such as 10, 25, 50, 75 and 90 % quantiles for each day considered. These quantiles will then be used as an indication of possible future streamflows.
Data available
The data set available to the rainfall-runoff modelling consists of:
- a time series of daily total precipitation (liquid + solid) [mm/day] (
Ptot); - a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/day] (
Evap); - a time series of daily streamflows expressed as a specific discharge [mm/day] (
Qmmd).
Command lines for the production of simulations
Loading and formatting of data
The following command lines are used to prepare the data required to calibrate the rainfall-runoff model, and to define the working periods (initialization period, calibration period and forecast period).
Initial time series
Ploted time series
Observed time series
Forecast time series
Analysis of streamflow climatology
The following command lines allow the analysis of the streamflow climatology by calculating the interannual streamflow quantiles over the forecasting period:
# Calculation of the historical streamflow quantiles ts_qclim_quant <- aggregate(Qmmd ~ MonthDay, data = ts_hist[ts_hist$MonthDay %in% month_day_fcst, ], FUN = function(x) { quantile(x, probs = c(0.10, 0.25, 0.50, 0.75, 0.90)) }) ts_qclim_quant <- as.data.frame(ts_qclim_quant$Qmmd) colnames(ts_qclim_quant) <- paste0("Q", gsub("\\D", "", colnames(ts_qclim_quant))) rownames(ts_qclim_quant) <- month_day_fcst # Display of the 1st calculated quantiles head(ts_qclim_quant)
The streamflow climatology allows to obtain a first order of magnitude of the possible streamflows of the next days (cf. following figure), but does not allow to account for the last observed streamflows.
<<v02_fig_pattern>> polygon(x = c(dates_fcst, rev(dates_fcst)), y = c(ts_qclim_quant$Q10, rev(ts_qclim_quant$Q90)), col = "lightgrey", border = NA) polygon(x = c(dates_fcst, rev(dates_fcst)), y = c(ts_qclim_quant$Q25, rev(ts_qclim_quant$Q75)), col = "darkgrey", border = NA) lines(x = dates_fcst, y = ts_qclim_quant$Q50, lwd = 2, col = "grey50") legend("topright", legend = c("Qobs", "Qclim (10,25,50,75,90 %)"), col = c("black", "grey50"), pch = c(NA, 15), lwd = c(1, 2), pt.cex = 2, lty = c(1, 1), pt.bg = c(NA, "lightgrey"), bg = "white")
Data processing for GR6J
The following command lines aim to prepare the available data for use by GR6J, using the PrepGR() function of the r airGRteaching package.
# Adding an epsilon to observed streamflows for criterion calculation epsilon_qobs <- mean(ts_hist$Qmmd, na.rm = TRUE) / 100 # Data processing for GR6J prep_hist <- PrepGR(DatesR = ts_hist$Date, Precip = ts_hist$Ptot, PotEvap = ts_hist$Evap, Qobs = ts_hist$Qmmd + epsilon_qobs, HydroModel = "GR6J", CemaNeige = FALSE)
GR6J calibration on the historical period
The following command lines are used to calibrate the GR6J model on the historical period.
# Calibration step cal_hist <- CalGR(PrepGR = prep_hist, CalCrit = "NSE", WupPer = per_wup_hist, CalPer = per_cal_hist, transfo = "log", verbose = TRUE) # Get parameter and criterion values param_cal_hist <- GetParam(cal_hist) crit_cal_hist <- GetCrit(cal_hist) # Graphical assessment of the calibration performance plot(cal_hist)
The six parameters and the value of the calibration criterion ($NSE_{ln}$) obtained after the automatic calibration procedure are:
- X1 =
r round(param_cal_hist[1], digits = 0)[mm] - X2 =
r round(param_cal_hist[2], digits = 3)[mm/day] - X3 =
r round(param_cal_hist[3], digits = 0)[mm] - X4 =
r round(param_cal_hist[4], digits = 3)[day] - X5 =
r round(param_cal_hist[5], digits = 3)[-] - X6 =
r round(param_cal_hist[6], digits = 2)[mm] - $NSE_{ln}$ =
r round(crit_cal_hist, digits = 3)[-]
The performance obtained in calibration is very good, since the criterion $NSE_{ln}$ is equal to r round(crit_cal_hist, digits = 3).
The following command lines allow to store in the same table the observed streamflows and the simulated streamflows with the set of parameters obtained by automatic calibration, in order to compare them.
# Combination of observed and simulated streamflow time series on the calibration period ts_cal_hist <- as.data.frame(cal_hist) head(ts_cal_hist)
The following figure represents the observed and simulated streamflow series over the end of the calibration period.
# Plot framework <<v02_fig_pattern>> # sim lines(x = ts_cal_hist$Date, y = ts_cal_hist$Qsim, col = "orangered") legend("topright", legend = c("Qobs", "Qsim"), lty = 1, col = c("black", "orangered"), bg = "white")
Pessimistic zero precipitation scenario
The following command lines are used to create a "dummy" time series that combines observed data (for model warm-up) and data generated for the "zero precipitation" scenario. This time series will then be used as input to the GR6J model to "transform" this meteorological scenario into a hydrological scenario.
# Duplicate the time series with future dates to fill the forecast period ts_fcst_p0 <- ts_fcst # Setting zero precipitation for the forecast period ts_fcst_p0$Ptot <- 0 # Addition of interannual average PE ts_eclim <- aggregate(Evap ~ MonthDay, data = ts_hist[ts_hist$MonthDay %in% month_day_fcst, ], FUN = mean) ts_fcst_p0$Evap <- ts_eclim$Evap # Combine historical and forecasting time series ts_scen_p0 <- rbind(ts_hist, ts_fcst_p0) # Display of the last lines tail(ts_scen_p0)
The following command lines allow us to format the previously created time series as inputs to the GR6J model and to perform a rainfall-runoff simulation.
# Data processing for GR6J prep_scen_p0 <- PrepGR(DatesR = ts_scen_p0$Date, Precip = ts_scen_p0$Ptot, PotEvap = ts_scen_p0$Evap, Qobs = NULL, HydroModel = "GR6J", CemaNeige = FALSE) # Hydrological forecast sim_scen_p0 <- SimGR(PrepGR = prep_scen_p0, WupPer = per_wup_fcst, SimPer = per_fcst, Param = param_cal_hist, verbose = FALSE) # Simulated streamflow time series ts_sim_scen_p0 <- as.data.frame(sim_scen_p0)
The following figure represents the observed and simulated streamflow series over the end of the calibration period, and the result of the simulation of the "zero precipitation" scenario. The slow decay of the simulated streamflows, generated by the slow emptying of the GR6J reservoirs, is observed.
# Plot framework <<v02_fig_step_cal_hist>> # Zero precipitation lines(x = ts_sim_scen_p0$Dates, y = ts_sim_scen_p0$Qsim, lwd = 2, col = "blue") # Legend legend("topright", legend = c("Qobs", "Qsim", "QscenP0"), lwd = 2, lty = 1, col = c("black", "orangered", "blue"), bg = "white")
It is important to note that the simulations performed so far do not take into account the last observed streamflow. This can generate a discontinuity in the sequence of streamflows, between the observations available until the time of the forecast, and the simulations produced from this time. This "error" can be corrected in several ways, by "assimilating" the last observed streamflow. The simplest correction consists in calculating a simple ratio between the last observed streamflow and the corresponding simulated streamflow, and to use this ratio to correct all the following simulated streamflows. The following command lines allow to calculate such a ratio, which will be used thereafter to correct the streamflows simulated by GR6J, by dividing them by this ratio (e.g. following figure).
# Correction (~ assimilation) corr_qsim <- ts_sim_scen_p0$Qsim[1] / ts_hist$Qmmd[nrow(ts_hist)] # Ratio display corr_qsim
# Plot framework <<v02_fig_pattern>> # Zero precipitation lines(x = ts_sim_scen_p0$Date, y = ts_sim_scen_p0$Qsim / corr_qsim, lwd = 2, col = "cornflowerblue") # Legend legend("topright", legend = c("Qobs", "QscenP0corr"), lwd = 2, lty = 1, col = c("black", "cornflowerblue"), bg = "white")
Non-zero future precipitation scenarios
The following command lines allow us to perform rainfall-runoff simulations for the studied summer using the last observed streamflow by "replaying" the precipitation climatology, i.e. the precipitation observed during other summers.
# Historical years year_hist <- unique(ts_hist$Year) year_hist <- setdiff(year_hist, unique(ts_fcst$Year)) year_hist # Loop on the historical years ts_qscen_year <- sapply(year_hist, FUN = function(i_year) { i_ts_hist <- ts_hist[ts_hist$Year == i_year & ts_hist$MonthDay %in% month_day_fcst, ] i_ts_hist$Date <- dates_fcst # Combine historical and forecasting (based on precipitation climatology) time series i_ts_scen_py <- rbind(ts_hist, i_ts_hist) # Data processing for GR6J prep_scen_py <- PrepGR(DatesR = i_ts_scen_py$Date, Precip = i_ts_scen_py$Ptot, PotEvap = i_ts_scen_py$Evap, HydroModel = "GR6J", CemaNeige = FALSE) # Hydrological forecast sim_scen_py <- SimGR(PrepGR = prep_scen_py, WupPer = per_wup_fcst, SimPer = per_fcst, Param = param_cal_hist, verbose = FALSE) # Correction of the simulated streamflow sim_scen_py$OutputsModel$Qsim / corr_qsim }) # Compute future streamflow quantiles (based on precipitation climatology) ts_qscen_quant <- t(apply(ts_qscen_year, MARGIN = 1, FUN = function(x) { quantile(x, probs = c(0.10, 0.25, 0.50, 0.75, 0.90)) })) ts_qscen_quant <- as.data.frame(ts_qscen_quant) colnames(ts_qscen_quant) <- paste0("Q", gsub("\\D", "", colnames(ts_qscen_quant))) rownames(ts_qscen_quant) <- month_day_fcst # Display of the first calculated quantiles head(ts_qscen_quant)
Streamflow forecasts based on the precipitation climatology suggest that the next few days will be dry, and that an increase in the streamflow of the river r name_riv is not expected until a few days.
# Plot framework <<v02_fig_pattern>> polygon(x = c(dates_fcst, rev(dates_fcst)), y = c(ts_qscen_quant$Q10, rev(ts_qscen_quant$Q90)), col = adjustcolor("green2", alpha.f = 0.3), border = NA) polygon(x = c(dates_fcst, rev(dates_fcst)), y = c(ts_qscen_quant$Q25, rev(ts_qscen_quant$Q75)), col = adjustcolor("green2", alpha.f = 0.5), border = NA) lines(x = dates_fcst, y = ts_qscen_quant$Q50, lwd = 2, col = "green4") legend("topright", legend = c("Qobs", "QscenPy (10,25,50,75,90 %)"), col = c("black", "green2"), pch = c(NA, 15), lwd = c(1, 2), pt.cex = 2, lty = c(1, 1), pt.bg = c(NA, "green2"), bg = "white")
The following figure allows to compare all the results on the same graph. It displays observed and simulated streamflows over the end of the calibration period, the streamflow climatology over the forecast period, the result of the simulation of the "zero precipitation" scenario, and the streamflows based on the precipitation climatology.
# Plot framework <<v02_fig_pattern>> par(new = TRUE) <<v02_fig_qclim>> par(new = TRUE) <<v02_fig_cal_hist>> par(new = TRUE) <<v02_fig_qscen_py>> par(new = TRUE) <<v02_fig_qscen_p0corr>> legend("topright", legend = c("Qobs", "Qsim", "QscenP0corr", "Qclim (10,25,50,75,90 %)", "QscenPy (10,25,50,75,90 %)"), col = c("black", "orangered", "cornflowerblue", "grey50", "green2"), pch = c(NA, NA, NA, 15, 15), lwd = c(1, 1, 2, 2, 2), pt.cex = 2, lty = c(1, 1, 1, 1, 1), pt.bg = c(NA, NA, NA, "grey50", "green2"), bg = "white")
References
Try the airGRteaching package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.