Commit c506f0d7 authored by unknown's avatar unknown
Browse files

v0.0.1 first repository

parents
.Rproj.user
.Rhistory
.RData
.Ruserdata
Package: airGRscholar
Type: Package
Title: Tools to simplify the use of the airGR hydrological package by students
Version: 0.0.1
Date: 2016-12-05
Authors@R: c(person("Olivier", "Delaigue", role = c("aut", "cre"), email = "airGR@irstea.fr"), person("Laurent", "Coron", role = c("aut")))
Author: Olivier Delaigue and Laurent Coron
Maintainer: Olivier Delaigue <airGR@irstea.fr>
Depends: airGR, xts, dygraphs, shiny
Description: This package provides 3 functions which allows to run very simply the airGR hydrological modelling steps. It also provides plotting function to help student to explore observed data and to interprete the results of calibration and simulation of the GR hydrological models.
License: GPL-2
NeedsCompilation: no
URL: https://webgr.irstea.fr/airGR/?lang=en
Encoding: UTF-8
Modelling functions:
====================
ObsGR Creation of the inputs required to the CalGR
and SimGR functions
CalGR Calibration algorithm optimises the error
criterion selected as objective function
SimGR Run GR hydrological models
Plotting functions:
====================
dyplot Interactive time series plotting of ObsGR, CalGR
and SimGR objects
plot.ObsGR Time series plotting GR objects from inputs
plot.CalGR Plotting of the evolution of parameters and
objective function during the calibration step,
and time series plotting GR objects from calibration outputs
plot.SimGR Synthetic plotting of model outputs
Shiny interface:
====================
shiny.SimGR4J Interactive web applications to run manually
the GR4J hydrological model
#####################################
## S3 methods ##
#####################################
S3method("plot","ObsGR")
S3method("plot","CalGR")
S3method("plot","SimGR")
S3method("dyplot","default")
S3method("dyplot","ObsGR")
S3method("dyplot","CalGR")
S3method("dyplot","SimGR")
#####################################
## Export ##
#####################################
export(CalGR)
export(ObsGR)
export(SimGR)
export(plot.ObsGR)
export(plot.CalGR)
export(plot.SimGR)
export(dyplot)
export(dyplot.default)
export(dyplot.ObsGR)
export(dyplot.CalGR)
export(dyplot.SimGR)
export(shiny.SimGR4J)
#####################################
## Import ##
#####################################
import(stats)
import(graphics)
import(grDevices)
import(utils)
import(airGR)
import(xts)
import(dygraphs)
.onLoad <- function(libname, pkgname){
citation("airGR")
citation("airGRscholar")
}
CalGR <- function(ObsGR, CalCrit = c("NSE", "KGE", "KGE2", "RMSE"),
WupPer = NULL, CalPer, transfo = c("", "sqrt", "log", "inv", "sort"), verbose = TRUE) {
if (! any(class(ObsGR) %in% "ObsGR")) {
stop("Non convenient data for argument \"ObsGR\". Must be of class \"ObsGR\"")
}
WupInd <- NULL
if (!is.null(WupPer)) {
WupPer <- as.POSIXct(WupPer, tz = "UTC")
if (any(is.na(WupPer))) {
stop("Non convenient date format for the warm-up period \"WupPer\"")
} else {
if (! (any(ObsGR$InputsModel$DatesR == WupPer[1]) & any(ObsGR$InputsModel$DatesR == WupPer[2]))) {
stop("Non convenient date for the warm-up period \"WupPer\"")
} else {
WupInd <- which(ObsGR$InputsModel$DatesR == WupPer[1]):which(ObsGR$InputsModel$DatesR == WupPer[2])
}
}
}
CalPer <- as.POSIXct(CalPer, tz = "UTC")
if (any(is.na(CalPer))) {
stop("Non convenient date format for the calibration period \"CalPer\"")
} else {
if (! (any(ObsGR$InputsModel$DatesR == CalPer[1]) & any(ObsGR$InputsModel$DatesR == CalPer[2]))) {
stop("Non convenient date for the calibration period \"CalPer\"")
} else {
CalInd <- which(ObsGR$InputsModel$DatesR == CalPer[1]):which(ObsGR$InputsModel$DatesR == CalPer[2])
}
}
if (! any(CalCrit %in% c("NSE", "KGE", "KGE2", "RMSE"))) {
stop("Non convenient efficiency criteria \"EffCrit\"")
} else {
CalCrit <- CalCrit[1L]
CalCrit <- sprintf("ErrorCrit_%s", CalCrit)
FUN_CRIT <- get(CalCrit)
}
if (! any(transfo %in% c("", "sqrt", "log", "inv", "sort"))) {
stop("Non convenient transformation \"transfo\"")
} else {
transfo <- transfo[1L]
}
MOD_opt <- CreateRunOptions(FUN_MOD = get(ObsGR$TypeModel), InputsModel = ObsGR$InputsModel,
IndPeriod_WarmUp = WupInd, IndPeriod_Run = CalInd, verbose = FALSE)
MOD_crt <- CreateInputsCrit(FUN_CRIT = FUN_CRIT, InputsModel = ObsGR$InputsModel,
RunOptions = MOD_opt, Qobs = ObsGR$Qobs[CalInd], transfo = transfo)
CAL_opt <- CreateCalibOptions(FUN_MOD = get(ObsGR$TypeModel), FUN_CALIB = Calibration_Michel)
CAL <- Calibration(InputsModel = ObsGR$InputsModel, RunOptions = MOD_opt,
InputsCrit = MOD_crt, CalibOptions = CAL_opt,
FUN_MOD = get(ObsGR$TypeModel), FUN_CRIT = FUN_CRIT,
FUN_CALIB = Calibration_Michel)
SIM <- RunModel(InputsModel = ObsGR$InputsModel, RunOptions = MOD_opt,
Param = CAL$ParamFinalR, FUN_MOD = get(ObsGR$TypeModel))
CalGR <- list(OptionsCalib = MOD_opt, Qobs = ObsGR$Qobs[CalInd],
OutputsCalib = CAL, OutputsModel = SIM,
TypeModel = ObsGR$TypeModel, CalCrit = CalCrit,
PeriodModel = list(WarmUp = as.POSIXct(ObsGR$InputsModel$DatesR[range(MOD_opt$IndPeriod_WarmUp)]),
Run = CalPer))
class(CalGR) <- c("CalGR", "GR")
return(CalGR)
}
\ No newline at end of file
ObsGR <- function(ObsBV = NULL, DatesR = NULL, Precip = NULL, PotEvap = NULL, Qobs = NULL, TempMean = NULL,
ZInputs = NULL, HypsoData = NULL, NLayers = 5,
TypeModel, CemaNeige = FALSE) {
if (is.null(ObsBV) && (is.null(DatesR) | is.null(Precip) | is.null(PotEvap) | is.null(Qobs))) {
stop("Missing input data")
}
if (ncol(ObsBV) >= 5) {
TempMean <- ObsBV[, 5L]
}
if (!is.null(TempMean)) {
TempMean <- TempMean
} else {
TempMean <- NA
}
if (is.null(ObsBV)) {
ObsBV <- data.frame(DatesR = DatesR,
Precip = Precip,
PotEvap = PotEvap,
Qobs = Qobs,
TempMean = TempMean)
}
if (!is.null(ObsBV)) {
ObsBV <- data.frame(DatesR = ObsBV[, 1L],
Precip = ObsBV[, 2L],
PotEvap = ObsBV[, 3L],
Qobs = ObsBV[, 4L],
TempMean = TempMean)
}
SuiteGR <- paste("GR", c("1A", "2M", "4J", "5J", "6J", "4H"), sep = "")
if (! any(TypeModel %in% SuiteGR)) {
stop("Non convenient model")
} else {
if (! CemaNeige) {
TypeModel <- sprintf("RunModel_%s", TypeModel)
}
if (CemaNeige && grepl("J", TypeModel)) {
TypeModel <- sprintf("RunModel_CemaNeige%s", TypeModel)
}
if (CemaNeige && !grepl("J", TypeModel)) {
warning("CemaNeige can not be used with ", TypeModel)
TypeModel <- sprintf("RunModel_%s", TypeModel)
}
FUN_MOD <- get(TypeModel)
}
MOD_obs <- CreateInputsModel(FUN_MOD = FUN_MOD, DatesR = ObsBV$DatesR,
Precip = ObsBV$Precip, PotEvap = ObsBV$PotEvap, TempMean = ObsBV$TempMean,
ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers, verbose = FALSE)
ObsGR <- list(InputsModel = MOD_obs, Qobs = ObsBV$Qobs, TypeModel = TypeModel)
class(ObsGR) <- c("ObsGR", "GR")
return(ObsGR)
}
SimGR <- function(ObsGR, CalGR = NULL, Param = NULL, EffCrit = c("NSE", "KGE", "KGE2", "RMSE"),
WupPer = NULL, SimPer, transfo = c("", "sqrt", "log", "inv", "sort"), verbose = TRUE) {
if (! any(class(ObsGR) %in% "ObsGR")) {
stop("Non convenient data for argument \"ObsGR\". Must be of class \"ObsGR\"")
}
if (! any(class(CalGR) %in% "CalGR") & !is.null(CalGR)) {
stop("Non convenient data for argument \"CalGR\". Must be of class \"CalGR\"")
}
if (is.null(CalGR) & is.null(Param)) {
stop("Arguments \"CalGR\" and \"Param\" are missing, with no default. You must fill in one of these two arguments.")
}
if (is.null(Param)) {
Param <- CalGR$OutputsCalib$ParamFinalR
}
WupInd <- NULL
if (!is.null(WupPer)) {
WupPer <- as.POSIXct(WupPer, tz = "UTC")
if (any(is.na(WupPer))) {
stop("Non convenient date format for the warm-up period \"WupPer\"")
} else {
if (! (any(ObsGR$InputsModel$DatesR == WupPer[1]) & any(ObsGR$InputsModel$DatesR == WupPer[1]))) {
stop("Non convenient date for the warm-up period \"WupPer\"")
} else {
WupInd <- which(ObsGR$InputsModel$DatesR == WupPer[1]):which(ObsGR$InputsModel$DatesR == WupPer[2])
}
}
}
SimPer <- as.POSIXct(SimPer, tz = "UTC")
if (any(is.na(SimPer))) {
stop("Non convenient date format for the simulation period \"SimPer\"")
} else {
if (! (any(ObsGR$InputsModel$DatesR == SimPer[1]) & any(ObsGR$InputsModel$DatesR == SimPer[2]))) {
stop("Non convenient date for the simulation period \"SimPer\"")
} else {
SimInd <- which(ObsGR$InputsModel$DatesR == SimPer[1]):which(ObsGR$InputsModel$DatesR == SimPer[2])
}
}
if (! any(EffCrit %in% c("NSE", "KGE", "KGE2", "RMSE"))) {
stop("Non convenient efficiency criteria \"EffCrit\"")
} else {
EffCrit <- EffCrit[1L]
EffCrit <- sprintf("ErrorCrit_%s", EffCrit)
FUN_CRIT <- get(EffCrit)
}
if (! any(transfo %in% c("", "sqrt", "log", "inv", "sort"))) {
stop("Non convenient transformation \"transfo\"")
} else {
transfo <- transfo[1L]
}
MOD_opt <- CreateRunOptions(FUN_MOD = get(ObsGR$TypeModel), InputsModel = ObsGR$InputsModel,
IndPeriod_WarmUp = WupInd, IndPeriod_Run = SimInd, verbose = verbose)
MOD_crt <- CreateInputsCrit(FUN_CRIT = FUN_CRIT, InputsModel = ObsGR$InputsModel,
RunOptions = MOD_opt, Qobs = ObsGR$Qobs[SimInd], transfo = transfo)
SIM <- RunModel(InputsModel = ObsGR$InputsModel, RunOptions = MOD_opt,
Param = Param, FUN_MOD = get(ObsGR$TypeModel))
CRT <- ErrorCrit(InputsCrit = MOD_crt, OutputsModel = SIM, FUN_CRIT = FUN_CRIT, verbose = verbose)
SimGR <- list(OptionsSimul = MOD_opt, OutputsModel = SIM, Qobs = ObsGR$Qobs[SimInd],
TypeModel = ObsGR$TypeModel,
CalCrit = CalGR$CalCrit, EffCrit = CRT,
PeriodModel = list(WarmUp = as.POSIXct(ObsGR$InputsModel$DatesR[range(MOD_opt$IndPeriod_WarmUp)]),
Run = SimPer))
class(SimGR) <- c("SimGR", "GR")
return(SimGR)
}
\ No newline at end of file
.onLoad <- function(libname, pkgname){
# packageStartupMessage("\nTo cite package 'airGRscholar' in publications use :\n", appendLF = T)
# packageStartupMessage(print(citation("airGR"), style = "textVersion"), appendLF = T)
# packageStartupMessage("\nvv\n", appendLF = T)
# packageStartupMessage(print(citation("airGRscholar"), style = "textVersion"), appendLF = T)
print(citation("airGR"))
print(citation("airGRscholar"))
}
\ No newline at end of file
dyplot.CalGR <- function(x, ...) {
if (! any(class(x) %in% "CalGR")) {
stop("Non convenient data for x argument. Must be of class \"CalGR\"")
}
dyplot.default(x, ...)
}
dyplot.ObsGR <- function(x, ...) {
if (! any(class(x) %in% "ObsGR")) {
stop("Non convenient data for x argument. Must be of class \"ObsGR\"")
}
dyplot.default(x, ...)
}
dyplot <- function(x, ...) {
UseMethod("dyplot")
}
dyplot.SimGR <- function(x, ...) {
if (! any(class(x) %in% "SimGR")) {
stop("Non convenient data for x argument. Must be of class \"SimGR\"")
}
dyplot.default(x, ...)
}
dyplot.default <- function(x, col.Precip = "royalblue", col.Q = c("black", "red"),
xlab = "Time", ylab = c("Precipitation [mm]", "Flow [mm]"), main = NULL,
RangeSelector = TRUE, dyLegendShow = c("follow", "auto", "always", "onmouseover", "never"), ...) {
if (! any(class(x) %in% c("ObsGR", "CalGR", "SimGR"))) {
stop("Non convenient data for x argument. Must be of class \"ObsGR\" \"CalGR\", or \"SimGR\"")
}
if (any(class(x) %in% "ObsGR")) {
data <- data.frame(DatesR = x$InputsModel$DatesR,
Precip = x$InputsModel$Precip,
Qobs = x$Qobs)
} else {
data <- data.frame(DatesR = x$OutputsModel$DatesR,
Precip = x$OutputsModel$Precip,
Qobs = x$Qobs,
Qsim = x$OutputsModel$Qsim)
if (length(col.Q) < 2) {
col.Q <- rep(col.Q, 2)
}
}
data.xts <- xts(data[, -1L], order.by = data$DatesR)
graphOut <- dygraph(data.xts, main = main)
graphOut <- dyAxis(dygraph = graphOut, name = "y", label = "Flow",
valueRange = range(data.xts[, grep("^Q", colnames(data.xts))], na.rm = TRUE) * c(0.01, 1.59))
graphOut <- dyAxis(dygraph = graphOut, name = "y2", label = "Precip.", independentTicks = FALSE,
valueRange = rev(range(data.xts[, "Precip"], na.rm = TRUE) * c(0.01, 2.99)))
graphOut <- dySeries(dygraph = graphOut, "Qobs" , axis = 'y' , color = col.Q[seq_along(grep("^Q", colnames(data.xts)))])
graphOut <- dySeries(dygraph = graphOut, "Precip", axis = 'y2', stepPlot = TRUE, fillGraph = TRUE, color = col.Precip)
# if (ModelPeriod) {
# graphOut <- dyShading(dygraph = graphOut, from = x$PeriodModel$WarmUp[1L], to = x$PeriodModel$WarmUp[2L], color = "red")
# graphOut <- dyShading(dygraph = graphOut, from = x$PeriodModel$Run[1L] , to = x$PeriodModel$Run[2L] , color = "blue")
# }
if (RangeSelector) {
graphOut <- dyRangeSelector(dygraph = graphOut, height = 15)
}
if (any(dyLegendShow %in% c("follow", "auto", "always", "onmouseover", "never"))) {
graphOut <- dyLegend(dygraph = graphOut, show = dyLegendShow[1L])
}
graphOut <- dyOptions(dygraph = graphOut, useDataTimezone = TRUE)
return(graphOut)
}
plot.CalGR <- function(x, which = c("both", "hist", "time"), ...) {
if (! any(class(x) %in% "CalGR")) {
stop("Non convenient data for x argument. Must be of class \"CalGR\"")
}
nbParamX <- as.numeric(gsub("\\D", "", x$TypeModel))
nbParamC <- ifelse(grepl("CemaNeige", x$TypeModel), 2, 0)
nbParam <- nbParamX + nbParamC
nmParam <- c(sprintf("X%i", 1:nbParamX), sprintf("C%i", seq_len(nbParamC)))
if (any(which[1L] %in% c("both", "hist"))) {
layout.list <- list(matrix(c(1:2), ncol = 2),
matrix(c(1:3, 3), ncol = 2),
matrix(NA),
matrix(c(1:5, 5), ncol = 3),
matrix(c(1:5, 0, 6, 6), ncol = 4),
matrix(c(1:7, 7), ncol = 4),
matrix(c(1:5, 0, 6:8, 8), ncol = 5),
matrix(c(1:9, 9), ncol = 5))
layout.mat <- layout.list[[nbParam]]
layout.w <- c(rep(1, ncol(layout.mat)-1), 2)
layout.h <- rep(1, nrow(layout.mat))
layout(mat = layout.mat, widths = layout.w, heights = layout.h)
}
if (which[1L] %in% c("both")) {
layout.mat <- rbind(layout.mat, max(layout.mat)+1, max(layout.mat)+2)
layout(mat = layout.mat, widths = layout.w, heights = layout.h)
}
if (which[1L] %in% c("time")) {
layout(mat = matrix(1:2), widths = c(1, 2), heights = c(1, 2))
}
# if (nbParam == 1) {
# layout(mat = matrix(c(1:2), ncol = 2), widths = c(1, 1), heights = c(1))
# }
# if (nbParam == 2) {
# layout(mat = matrix(c(1:3, 3), ncol = 2), widths = c(1, 2), heights = c(1, 1))
# }
# if (nbParam == 4) {
# layout(mat = matrix(c(1:5, 5), ncol = 3), widths = c(1, 1, 2), heights = c(1, 1))
# }
# if (nbParam == 5) {
# layout(mat = matrix(c(1:5, 0, 6, 6), ncol = 4), widths = c(1, 1, 1, 2), heights = c(1, 1, 1))
# }
# if (nbParam == 6) {
# layout(mat = matrix(c(1:7, 7), ncol = 4), widths = c(1, 1, 1, 2), heights = c(1, 1, 1))
# }
# if (nbParam == 7) {
# layout(mat = matrix(c(1:5, 0, 6:8, 8), ncol = 5), widths = c(1, 1, 1, 1, 2), heights = c(1, 1, 1, 1))
# }
# if (nbParam == 8) {
# layout(mat = matrix(c(1:9, 9), ncol = 5), widths = c(1, 1, 1, 1, 2), heights = c(1, 1, 1, 1))
# }
if (any(which[1L] %in% c("both", "hist"))) {
par(mar = c(2.5, 2.5, 4.5, 1))
for (i in seq_len(ncol(x$OutputsCalib$HistParamR))) {
plot(x$OutputsCalib$HistParamR[, i], col = ifelse(grepl("X", nmParam[i]), "darkorchid3", "darkolivegreen3"),
type = "b", pch = 19, cex = 0.8,
xlab = "", ylab = "", main = nmParam[i])
}
plot(x$OutputsCalib$HistCrit, col = "tomato",
type = "b", pch = 19, cex = 1.0,
xlab = "", ylab = "", main = "Efficiency criteria\nduring the steepest-descent step")
}
if (any(which[1L] %in% c("both", "time"))) {
data <- data.frame(DatesR = x$OutputsModel$DatesR,
Precip = x$OutputsModel$Precip,
Qobs = x$Qobs,
Qsim = x$OutputsModel$Qsim)
par(mar = c(0.1, 4, 4, 1), xaxt = "n")
plot(Precip ~ DatesR, data = data, type = "h", col = "royalblue",
xlab = "", ylab = "Precip.", main = "", ylim = rev(range(data$Precip)))
par(mar = c(5, 4, 0.1, 1), xaxt = "s")
matplot(data$DatesR, data[, c("Qobs", "Qsim")], type = "l", col = c("black", "red"), lty = 1,
xlab = "Time", ylab = "Flow", main = "", xaxt = "n")
axis.POSIXct(side = 1, x = data$DatesR)
# if (plot.na) {
# axis(side = 1, at = as.POSIXct(data$DatesR[is.na(data$Qobs)]),
# labels = FALSE, lwd.ticks = 3, col.ticks = col.na, tck = 0.025, lend = "butt")
#
# legend("topright", legend = "NA", pch = 15, col = col.na, bty = "n", cex = 0.8)
# }
}
}
\ No newline at end of file
plot.ObsGR <- function(x, type = "l", col.Precip = "royalblue", col.Q = "black", col.na = "grey",
xlab = "Time", ylab = c("Precipitation [mm]", "Flow [mm]"), main = "", plot.na = TRUE, ...) {
if (! any(class(x) %in% "ObsGR")) {
stop("Non convenient data for x argument. Must be of class \"ObsGR\"")
}
data <- data.frame(DatesR = x$InputsModel$DatesR,
Precip = x$InputsModel$Precip,
Qobs = x$Qobs)
layout(mat = matrix(1:2), widths = c(1, 2), heights = c(1, 2))
par(mar = c(0.1, 4, 4, 2), xaxt = "n")
plot(Precip ~ DatesR, data = data, type = "h", col = col.Precip,
xlab = "", ylab = ylab[1L], main = "", ylim = rev(range(data$Precip)))
par(mar = c(5, 4, 0.1, 2), xaxt = "s")
plot(Qobs ~ DatesR, data = data, type = type, col = col.Q,
xlab = xlab, ylab = ylab[2L], main = main)
if (plot.na) {
axis(side = 1, at = as.POSIXct(data$DatesR[is.na(data$Qobs)]),
labels = FALSE, lwd.ticks = 3, col.ticks = col.na, tck = 0.025, lend = "butt")
legend("topright", legend = "NA", pch = 15, col = col.na, bty = "n", cex = 0.8)
}
}
\ No newline at end of file
plot.SimGR <- function(x, log_scale = FALSE) {
if (! any(class(x) %in% "SimGR")) {
stop("Non convenient data for x argument. Must be of class \"SimGR\"")
}
plot_OutputsModel(x$OutputsModel, Qobs = x$Qobs, log_scale = log_scale)
}
\ No newline at end of file
shiny.SimGR4J <- function(ObsGR, Param = c(200, 0, 100, 2), WupPer = NULL, SimPer) {
# assign("ObsGR" , ObsGR , envir = .GlobalEnv)
# assign("Param" , Param , envir = .GlobalEnv)
# assign("WupPer", WupPer, envir = .GlobalEnv)
# assign("SimPer", SimPer, envir = .GlobalEnv)
.GlobalEnv$.SimGR4J.args <- list(ObsGR = ObsGR, Param = Param, WupPer = WupPer, SimPer = SimPer)
on.exit(rm(.SimGR4J.args, envir = .GlobalEnv))
shiny::runApp(system.file("ManCalGR", package = "airGRscholar"), launch.browser = TRUE)
# rm("ObsGR" , envir = .GlobalEnv)
return(NULL)
}
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: knitr
LaTeX: pdfLaTeX
# server.R
shinyServer(function(input, output) {