diff --git a/DESCRIPTION b/DESCRIPTION
index a3f44ba4304c2df39e5b26a1f950bb3aa509daae..dc992881e5161950faf7456c729df520a7a232de 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,7 +19,7 @@ URL: https://airgriwrm.g-eau.fr/,
BugReports: https://github.com/inrae/airGRiwrm/issues
Depends:
airGR (>= 1.7.0),
- R (>= 2.10)
+ R (>= 3.5)
Imports:
DiagrammeR,
dplyr,
diff --git a/NAMESPACE b/NAMESPACE
index 2117129665dfac85f6d3978cda702541e67292c5..a6e3b9618b3e40ac965d2462cf710b9d5b24c35a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -38,9 +38,12 @@ export(RunModel)
export(getNoSD_Ids)
export(getSD_Ids)
export(isNodeDownstream)
+export(plot.Qm3s)
import(airGR)
importFrom(grDevices,rainbow)
importFrom(graphics,matplot)
importFrom(graphics,par)
importFrom(graphics,plot)
importFrom(graphics,title)
+importFrom(utils,read.table)
+importFrom(utils,tail)
diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
index 1924a775093a78c6d98b945ffc2c2af30f0ce171..f4671e11d5bcc9e480eb6d6501a5334bed607d36 100644
--- a/R/Calibration.GRiwrmInputsModel.R
+++ b/R/Calibration.GRiwrmInputsModel.R
@@ -36,7 +36,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
for(id in gaugedIds) {
IM <- InputsModel[[id]]
- message("Calibration.GRiwrmInputsModel: Treating sub-basin ", id, "...")
+ message("Calibration.GRiwrmInputsModel: Processing sub-basin ", id, "...")
if (inherits(InputsCrit[[id]], "InputsCritLavenneFunction")) {
IC <- getInputsCrit_Lavenne(id, OutputsModel, InputsCrit)
@@ -54,7 +54,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
IM$FUN_MOD <- "RunModel_Ungauged"
attr(RunOptions[[id]], "GRiwrmRunOptions") <- l$RunOptions
} else {
- if (useUpstreamQsim && any(IM$UpstreamIsRunoff)) {
+ if (useUpstreamQsim && any(IM$UpstreamIsModeled)) {
# Update InputsModel$Qupstream with simulated upstream flows
IM <- UpdateQsimUpstream(IM, RunOptions[[id]], OutputsModel)
}
@@ -188,13 +188,13 @@ updateParameters4Ungauged <- function(GaugedId,
g2 <- griwrm[griwrm$donor == GaugedId,]
upIds2 <- g2$id[!g2$id %in% g2$down]
for (id in upIds2) {
- if(useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsRunoff)) {
+ if(useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsModeled)) {
# Update InputsModel$Qupstream with simulated upstream flows
InputsModel[[id]] <- UpdateQsimUpstream(InputsModel[[id]],
RunOptions[[id]],
OutputsModel)
- InputsModel[[id]]$UpstreamIsRunoff <-
- rep(FALSE, length(InputsModel[[id]]$UpstreamIsRunoff))
+ InputsModel[[id]]$UpstreamIsModeled <-
+ rep(FALSE, length(InputsModel[[id]]$UpstreamIsModeled))
}
}
diff --git a/R/CreateController.R b/R/CreateController.R
index c1fe654023c6eb57a2ca05916b9c1a0f5ee0a0ff..4079240e35f642e580f9cdc11e2d7f1a743b17d3 100644
--- a/R/CreateController.R
+++ b/R/CreateController.R
@@ -26,40 +26,20 @@
#' - `FUN` [function]: controller logic which calculates `U` from `Y`
#' @export
#'
-#' @examples
-#' # First create a Supervisor from a model
-#' data(Severn)
-#' nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
-#' nodes$model <- "RunModel_GR4J"
-#' griwrm <- CreateGRiwrm(nodes,
-#' list(id = "gauge_id",
-#' down = "downstream_id",
-#' length = "distance_downstream"))
-#' BasinsObs <- Severn$BasinsObs
-#' DatesR <- BasinsObs[[1]]$DatesR
-#' PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
-#' PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
-#' Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
-#' Precip <- ConvertMeteoSD(griwrm, PrecipTot)
-#' PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-#' InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
-#' sv <- CreateSupervisor(InputsModel)
-#'
-#' # A controller which usually releases 0.1 m3/s and provides
-#' # extra release if the downstream flow is below 0.5 m3/s
-#' logicDamRelease <- function(Y) max(0.5 - Y[1], 0.1)
-#' CreateController(sv, "DamRelease", Y = c("54001"), U = c("54095"), FUN = logicDamRelease)
+#' @example man-examples/RunModel.Supervisor.R
CreateController <- function(supervisor, ctrl.id, Y, U, FUN){
- if(!is.character(ctrl.id)) stop("Parameter `ctrl.id` should be character")
+ if (!is.character(ctrl.id)) stop("Parameter `ctrl.id` should be character")
+ if (length(ctrl.id) != 1) stop("Parameter `ctrl.id` should be of length 1")
+ stopifnot(is.Supervisor(supervisor))
FUN <- match.fun(FUN)
ctrlr <- list(
id = ctrl.id,
- U = CreateControl(U),
+ U = CreateControl(U, supervisor, TRUE),
Unames = U,
- Y = CreateControl(Y),
+ Y = CreateControl(Y, supervisor, FALSE),
Ynames = Y,
FUN = FUN
)
@@ -87,9 +67,19 @@ CreateController <- function(supervisor, ctrl.id, Y, U, FUN){
#' @examples
#' # For pointing the discharge at the oulet of basins "54095" and "54002"
#' CreateControl(c("54095", "54002"))
-CreateControl <- function(locations) {
- if(!is.character(locations)) {
- stop("Parameter `locations` should be character")
+CreateControl <- function(locations, sv, isU) {
+ if (!is.character(locations)) {
+ stop("Parameters `Y` and `U` should be character")
+ }
+ if (isU) {
+ if (!all(locations %in% sv$griwrm4U$id)) {
+ stop("Ids defined in `U` must be chosen from DirectInjection and Diversion nodes: ",
+ paste(sv$griwrm4U$id, collapse = ", "))
+ }
+ } else {
+ if (!all(locations %in% sv$griwrm$id)) {
+ stop("Ids defined in `Y` must be chosen from available Ids in the GRiwrm object")
+ }
}
m <- matrix(NA, ncol = length(locations), nrow = 0)
return(m)
diff --git a/R/CreateGRiwrm.R b/R/CreateGRiwrm.R
index 83cf2672d10e128e72558c0440d090f24ea0db4d..2c46b0ef1f5a829609038ca96cdb78198c106022 100644
--- a/R/CreateGRiwrm.R
+++ b/R/CreateGRiwrm.R
@@ -17,11 +17,16 @@
#'
#' * One of the hydrological models available in the *airGR* package defined by its
#' `RunModel` function (i.e.: `RunModel_GR4J`, `RunModel_GR5HCemaneige`...)
-#' * `NA` for injecting (or abstracting) a flow time series at the location of the node
-#' (direct flow injection)
#' * `Ungauged` for an ungauged node. The sub-basin inherits hydrological model and
#' parameters from a "donor" sub-basin. By default the donor is the first gauged
#' node at downstream
+#' * `NA` for injecting (or abstracting) a flow time series at the location of the node
+#' (direct flow injection)
+#' * `Diversion` for abstracting a flow time series from an existing node transfer it
+#' to another node. As a `Diversion` is attached to an existing node, this node is
+#' then described with 2 lines: one for the hydrological model and another one for the
+#' diversion
+#'
#'
#' @param db [data.frame] description of the network (See details)
#' @param cols [list] or [vector] columns of `db`. By default, mandatory column
@@ -131,7 +136,7 @@ getNodeRanking <- function(griwrm) {
stop("getNodeRanking: griwrm argument should be of class GRiwrm")
}
# Remove upstream nodes without model (direct flow connections)
- griwrm <- griwrm[!is.na(griwrm$model),]
+ griwrm <- griwrm[!is.na(griwrm$model), ]
# Rank 1
rank <- setdiff(griwrm$id, griwrm$down)
ranking <- rank
@@ -149,19 +154,69 @@ getNodeRanking <- function(griwrm) {
checkNetworkConsistency <- function(db) {
- if(sum(is.na(db$down)) != 1 | sum(is.na(db$length)) != 1) {
- stop("One and only one node must have 'NA' in columns 'down' and 'length")
+ db2 <- db[getDiversionRows(db, TRUE), ]
+ if (any(duplicated(db2$id))) {
+ stop("Duplicated nodes detected: ",
+ paste(db2$id[duplicated(db2$id)], collapse = "\n"),
+ "\nNodes `id` must be unique (except for `Diversion` nodes)")
}
- if(which(is.na(db$down)) != which(is.na(db$length))) {
- stop("The node with 'down = NA' must be the same as the one with 'length = NA'")
+ if (sum(is.na(db$down)) == 0) {
+ stop("At least one node must be a network downstream node",
+ " specified by 'down = NA'")
}
sapply(db$down[!is.na(db$down)], function(x) {
- if(!(x %in% db$id)) {
+ if (!(x %in% db$id)) {
stop("The 'down' id ", x, " is not found in the 'id' column")
}
})
+ db3 <- db2[!is.na(db2$model), ]
+ sapply(db$id[getDiversionRows(db)], function(x) {
+ i <- which(db$id == x & db$model == "Diversion")[1]
+ if (length(which(db3$id == x)) != 1) {
+ nodeError(db[i, ],
+ "A Diversion node must have the same `id` of one (and only one) node with a model")
+ }
+ if (length(unique(db$down[db$id == x])) != 2) {
+ nodeError(db[i, ], paste(
+ "The downstream node of a Diversion node must be different",
+ "than the downstream node of the node is attached to"))
+ }
+ })
+ apply(db, 1, checkNode, simplify = FALSE)
}
+checkNode <- function(node) {
+ node <- as.list(node)
+ if (!is.na(node$model)) {
+ if (node$model == "Diversion") {
+ if (!is.na(node$area)) {
+ nodeError(node, "A Diversion node must have its area equal to `NA`")
+ }
+ } else if (length(grep("RunModel_GR", node$model)) > 0 & is.na(node$area)) {
+ # TODO This test should be extended to airGRplus models
+ nodeError(node, "A node using an hydrological model must have a numeric area")
+ }
+ }
+ if (is.na(node$down) & !is.na(node$length)) {
+ nodeError(node, "A downstream end node defined by `down=NA` must have `length=NA`")
+ }
+ if (is.na(node$length) & !is.na(node$down)) {
+ nodeError(node, "A node with a defined downstream node must have a numeric `length`")
+ }
+}
+
+displayNodeDetails <- function(node) {
+ s <- sapply(names(node), function(x) {
+ sprintf("%s: %s", x, node[x])
+ })
+ paste("Error on the node:",
+ paste(s, collapse = "\n"),
+ sep = "\n")
+}
+
+nodeError <- function(node, s) {
+ stop(displayNodeDetails(node), "\n", s)
+}
#' Get the Id of the nearest gauged model at downstream
#'
@@ -172,12 +227,44 @@ checkNetworkConsistency <- function(db) {
#'
#' @noRd
getGaugedId <- function(id, griwrm) {
- if(!is.na(griwrm$model[griwrm$id == id]) & griwrm$model[griwrm$id == id] != "Ungauged") {
+ griwrm <- griwrm[getDiversionRows(griwrm, TRUE), ]
+ if (!is.na(griwrm$model[griwrm$id == id]) &
+ griwrm$model[griwrm$id == id] != "Ungauged") {
return(id)
- } else if(!is.na(griwrm$down[griwrm$id == id])){
+ } else if (!is.na(griwrm$down[griwrm$id == id])) {
return(getGaugedId(griwrm$down[griwrm$id == id], griwrm))
} else {
- stop("The model of the downstream node of a network cannot be `NA` or \"Ungauged\"")
+ stop("The model of the downstream node of a network",
+ " cannot be `NA` or \"Ungauged\"")
}
}
+getDiversionRows <- function(griwrm, inverse = FALSE) {
+
+ rows <- which(!is.na(griwrm$model) & griwrm$model == "Diversion")
+ if (inverse) {
+ if(length(rows) == 0) {
+ rows <- seq.int(nrow(griwrm))
+ } else {
+ rows <- setdiff(seq.int(nrow(griwrm)), rows)
+ }
+ }
+ return(rows)
+}
+
+getNodeProperties <- function(id, griwrm) {
+ upstreamIds <- griwrm$id[!griwrm$id %in% griwrm$down]
+ gaugedIds <- griwrm$id[!is.na(griwrm$model) & griwrm$model != "Ungauged"]
+ divertedIds <- griwrm$id[!is.na(griwrm$model) & griwrm$model == "Diversion"]
+ p <- list(
+ position = ifelse(id %in% upstreamIds, "Upstream", "Intermediate"),
+ hydrology = ifelse(id %in% gaugedIds, "Gauged",
+ ifelse(is.na(griwrm$model[griwrm$id == id]),
+ "DirectInjection",
+ "Ungauged"))
+ )
+ p$Upstream <- p$position == "Upstream"
+ p$DirectInjection = p$hydrology == "DirectInjection"
+ p$Diversion = id %in% divertedIds
+ return(p)
+}
diff --git a/R/CreateInputsCrit.GRiwrmInputsModel.R b/R/CreateInputsCrit.GRiwrmInputsModel.R
index 55394a1e20449f443256372a1746cf1456918cd9..b0eff5ba007792f48a33da26e5fd67e67cae16bd 100644
--- a/R/CreateInputsCrit.GRiwrmInputsModel.R
+++ b/R/CreateInputsCrit.GRiwrmInputsModel.R
@@ -1,5 +1,6 @@
#' @rdname CreateInputsCrit
#' @import airGR
+#' @importFrom utils tail read.table
#' @export
CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
FUN_CRIT = ErrorCrit_NSE,
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index db34acbbe7c3464f0eaec69fd20cddcf6f8ec031..557ae31009ebe931e61a1c94609bc9bb413625ec 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -2,37 +2,78 @@
#'
#' @param x \[GRiwrm object\] diagram of the semi-distributed model (See [CreateGRiwrm])
#' @param DatesR [POSIXt] vector of dates
-#' @param Precip (optional) [matrix] or [data.frame] frame of numeric containing precipitation in \[mm per time step\]. Column names correspond to node IDs
-#' @param PotEvap (optional) [matrix] or [data.frame] frame of numeric containing potential evaporation \[mm per time step\]. Column names correspond to node IDs
-#' @param Qobs (optional) [matrix] or [data.frame] frame of numeric containing observed flows in \[mm per time step\]. Column names correspond to node IDs
-#' @param PrecipScale (optional) named [vector] of [logical] indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default `TRUE` (the mean of the precipitation is kept to the original value)
-#' @param TempMean (optional) [matrix] or [data.frame] of time series of mean air temperature \[°C\], required to create the CemaNeige module inputs
-#' @param TempMin (optional) [matrix] or [data.frame] of time series of minimum air temperature \[°C\], possibly used to create the CemaNeige module inputs
-#' @param TempMax (optional) [matrix] or [data.frame] of time series of maximum air temperature \[°C\], possibly used to create the CemaNeige module inputs
-#' @param ZInputs (optional) named [vector] of [numeric] giving the mean elevation of the Precip and Temp series (before extrapolation) \[m\], possibly used to create the CemaNeige module input
-#' @param HypsoData (optional) [matrix] or [data.frame] containing 101 [numeric] rows: min, q01 to q99 and max of catchment elevation distribution \[m\], if not defined a single elevation is used for CemaNeige
-#' @param NLayers (optional) named [vector] of [numeric] integer giving the number of elevation layers requested [-], required to create CemaNeige module inputs, default=5
+#' @param Precip (optional) [matrix] or [data.frame] of [numeric] containing
+#' precipitation in \[mm per time step\]. Column names correspond to node IDs
+#' @param PotEvap (optional) [matrix] or [data.frame] of [numeric] containing
+#' potential evaporation \[mm per time step\]. Column names correspond to node IDs
+#' @param Qobs (optional) [matrix] or [data.frame] of [numeric] containing
+#' observed flows. It must be provided only for nodes of type "Direct
+#' injection" and "Diversion". See [CreateGRiwrm] for
+#' details about these node types. Unit is \[mm per time step\] for nodes
+#' with an area, and \[m3 per time step\] for nodes with `area=NA`.
+#' Column names correspond to node IDs. Negative flows are abstracted from
+#' the model and positive flows are injected to the model
+#' @param Qmin (optional) [matrix] or [data.frame] of [numeric] containing
+#' minimum flows to let downstream of a node with a Diversion \[m3 per
+#' time step\]. Default is zero. Column names correspond to node IDs
+#' @param PrecipScale (optional) named [vector] of [logical] indicating if the
+#' mean of the precipitation interpolated on the elevation layers must be
+#' kept or not, required to create CemaNeige module inputs, default `TRUE`
+#' (the mean of the precipitation is kept to the original value)
+#' @param TempMean (optional) [matrix] or [data.frame] of time series of mean
+#' air temperature \[°C\], required to create the CemaNeige module inputs
+#' @param TempMin (optional) [matrix] or [data.frame] of time series of minimum
+#' air temperature \[°C\], possibly used to create the CemaNeige module inputs
+#' @param TempMax (optional) [matrix] or [data.frame] of time series of maximum
+#' air temperature \[°C\], possibly used to create the CemaNeige module inputs
+#' @param ZInputs (optional) named [vector] of [numeric] giving the mean
+#' elevation of the Precip and Temp series (before extrapolation) \[m\],
+#' possibly used to create the CemaNeige module input
+#' @param HypsoData (optional) [matrix] or [data.frame] containing 101 [numeric]
+#' rows: min, q01 to q99 and max of catchment elevation distribution \[m\],
+#' if not defined a single elevation is used for CemaNeige
+#' @param NLayers (optional) named [vector] of [numeric] integer giving the number
+#' of elevation layers requested [-], required to create CemaNeige module
+#' inputs, default=5
#' @param ... used for compatibility with S3 methods
#'
-#' @details Meteorological data are needed for the nodes of the network that represent a catchment simulated by a rainfall-runoff model. Instead of [airGR::CreateInputsModel] that has [numeric] [vector] as time series inputs, this function uses [matrix] or [data.frame] with the id of the sub-catchment as column names. For single values (`ZInputs` or `NLayers`), the function requires named [vector] with the id of the sub-catchment as name item. If an argument is optional, only the column or the named item has to be provided.
+#' @details Meteorological data are needed for the nodes of the network that
+#' represent a catchment simulated by a rainfall-runoff model. Instead of
+#' [airGR::CreateInputsModel] that has [numeric] [vector] as time series inputs,
+#' this function uses [matrix] or [data.frame] with the id of the sub-catchment
+#' as column names. For single values (`ZInputs` or `NLayers`), the function
+#' requires named [vector] with the id of the sub-catchment as name item. If an
+#' argument is optional, only the column or the named item has to be provided.
#'
#' See [airGR::CreateInputsModel] documentation for details concerning each input.
#'
-#' @return A \emph{GRiwrmInputsModel} object which is a list of \emph{InputsModel} objects created by [airGR::CreateInputsModel] with one item per modeled sub-catchment.
+#' Number of rows of `Precip`, `PotEvap`, `Qobs`, `Qmin`, `TempMean`, `TempMin`,
+#' `TempMax` must be the same of the length of `DatesR` (each row corresponds to
+#' a time step defined in `DatesR`).
+#'
+#' For example of use of Direct Injection nodes, see vignettes
+#' "V03_Open-loop_influenced_flow" and "V04_Closed-loop_regulated_withdrawal".
+#'
+#' For example of use of Diversion nodes, see example below and vignette
+#' "V06_Modelling_regulated_diversion".
+#'
+#' @return A \emph{GRiwrmInputsModel} object which is a list of \emph{InputsModel}
+#' objects created by [airGR::CreateInputsModel] with one item per modeled sub-catchment.
#' @export
-#' @inherit RunModel.GRiwrmInputsModel return examples
+#' @example man-examples/RunModel.GRiwrmInputsModel.R
#'
CreateInputsModel.GRiwrm <- function(x, DatesR,
Precip = NULL,
PotEvap = NULL,
Qobs = NULL,
+ Qmin = NULL,
PrecipScale = TRUE,
TempMean = NULL, TempMin = NULL,
TempMax = NULL, ZInputs = NULL,
HypsoData = NULL, NLayers = 5, ...) {
# Check and format inputs
- varNames <- c("Precip", "PotEvap", "TempMean", "Qobs",
+ varNames <- c("Precip", "PotEvap", "TempMean", "Qobs", "Qmin",
"TempMin", "TempMax", "ZInputs", "HypsoData", "NLayers")
names(varNames) <- varNames
lapply(varNames, function(varName) {
@@ -51,17 +92,32 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
))
}
if (!varName %in% c("ZInputs", "NLayers", "HypsoData") && nrow(v) != length(DatesR)) {
- stop("'%s' number of rows and the length of 'DatesR' must be equal",
- varName)
+ stop(sprintf(
+ "'%s' number of rows and the length of 'DatesR' must be equal",
+ varName
+ ))
+ }
+ if (varName %in% c("Precip", "PotEvap", "Qmin")) {
+ if (any(is.na(v))) {
+ stop(sprintf(
+ "`NA` values detected in '%s'. Missing values are not allowed in InputsModel",
+ varName
+ ))
+ }
+ if (any(v < 0)) {
+ stop(sprintf(
+ "'%s' values must be positive or nul. Missing values are not allowed in InputsModel",
+ varName
+ ))
+ }
}
-
} else if (!varName %in% c("ZInputs", "NLayers")) {
stop(sprintf("'%s' must be a matrix or a data.frame", varName))
}
}
})
- directFlowIds <- x$id[is.na(x$model)]
+ directFlowIds <- x$id[is.na(x$model) | x$model == "Diversion"]
if (length(directFlowIds) > 0) {
err <- FALSE
if (is.null(Qobs)) {
@@ -69,17 +125,63 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
} else {
Qobs <- as.matrix(Qobs)
if (is.null(colnames(Qobs))) {
- err <- TRUE
+ err <- TRUE
} else if (!all(directFlowIds %in% colnames(Qobs))) {
- err <- TRUE
+ err <- TRUE
}
}
if (err) stop(sprintf("'Qobs' column names must at least contain %s", paste(directFlowIds, collapse = ", ")))
}
+ if (!all(colnames(Qobs) %in% directFlowIds)) {
+ stop("The following columns in 'Qobs' don't match with ",
+ "Direction Injection or Diversion nodes: ",
+ paste(setdiff(colnames(Qobs), directFlowIds), collapse = ", "))
+ Qobs <- Qobs[directFlowIds, ]
+ }
+ diversionRows <- getDiversionRows(x)
+ if (length(diversionRows) > 0) {
+ warn <- FALSE
+ if (is.null(Qmin)) {
+ warn <- TRUE
+ } else {
+ Qmin <- as.matrix(Qmin)
+ if (!all(colnames(Qmin) %in% x$id[diversionRows])) {
+ stop(paste(
+ "'Qmin' contains columns that does not match with IDs of Diversion nodes:\n",
+ setdiff(colnames(Qmin), x$id[diversionRows])
+ ))
+ }
+ if (is.null(colnames(Qmin))) {
+ warn <- TRUE
+ } else if (!all(x$id[diversionRows] %in% colnames(Qmin))) {
+ warn <- TRUE
+ }
+ if (any(is.na(Qmin))) {
+ stop("`NA` values are note allowed in 'Qmin'")
+ }
+ }
+ if (warn) {
+ warning(
+ sprintf(
+ "'Qmin' would include the following columns %s.\n Zero values are applied by default.",
+ paste(directFlowIds, collapse = ", ")
+ )
+ )
+ }
+ # Qmin completion for Diversion nodes with default zero values
+ Qmin0 <- matrix(0, nrow = length(DatesR), ncol = length(diversionRows))
+ colnames(Qmin0) <- x$id[diversionRows]
+ if (is.null(Qmin)) {
+ Qmin <- Qmin0
+ } else {
+ Qmin0[, colnames(Qmin)] <- Qmin
+ Qmin <- Qmin0
+ }
+ }
InputsModel <- CreateEmptyGRiwrmInputsModel(x)
- # Qobs completion
+ # Qobs completion for at least filling Qupstream of all nodes by zeros
Qobs0 <- matrix(0, nrow = length(DatesR), ncol = nrow(x))
colnames(Qobs0) <- x$id
if (is.null(Qobs)) {
@@ -90,11 +192,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
}
for(id in getNodeRanking(x)) {
- message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
- if (x$area[x$id == id] > 0 && any(Qobs[, id] < 0, na.rm = TRUE)) {
- stop(sprintf("Negative flow found in 'Qobs[, \"%s\"]'. ", id),
- "Catchment flow can't be negative, use `NA` for flow data gaps.")
- }
+ message("CreateInputsModel.GRiwrm: Processing sub-basin ", id, "...")
InputsModel[[id]] <-
CreateOneGRiwrmInputsModel(id = id,
@@ -109,7 +207,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
ZInputs = getInputBV(ZInputs, id),
HypsoData = getInputBV(HypsoData, id),
NLayers = getInputBV(NLayers, id, 5),
- Qobs = Qobs
+ Qobs = Qobs,
+ Qmin = getInputBV(Qmin, id)
)
}
attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
@@ -145,24 +244,39 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#'
#' @return \emph{InputsModel} object for one.
#' @noRd
-CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
+CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
+ hasDiversion <- getNodeProperties(id, griwrm)$Diversion
+ if (hasDiversion) {
+ rowDiv <- which(griwrm$id == id & griwrm$model == "Diversion")
+ diversionOutlet <- griwrm$down[rowDiv]
+ griwrm <- griwrm[-rowDiv, ]
+ }
node <- griwrm[griwrm$id == id,]
FUN_MOD <- griwrm$model[griwrm$id == griwrm$donor[griwrm$id == id]]
# Set hydraulic parameters
- UpstreamNodes <- griwrm$id[griwrm$down == id & !is.na(griwrm$down)]
+ UpstreamNodeRows <- which(griwrm$down == id & !is.na(griwrm$down))
Qupstream <- NULL
LengthHydro <- NULL
BasinAreas <- NULL
- if(length(UpstreamNodes) > 0) {
+ if(length(UpstreamNodeRows) > 0) {
# Sub-basin with hydraulic routing
- Qupstream <- as.matrix(Qobs[ , UpstreamNodes, drop=FALSE])
- LengthHydro <- griwrm$length[griwrm$id %in% UpstreamNodes]
- names(LengthHydro) <- UpstreamNodes
+ Qupstream <- as.matrix(Qobs[ , griwrm$id[UpstreamNodeRows], drop=FALSE])
+ upstreamDiversion <- which(
+ sapply(griwrm$id[UpstreamNodeRows],
+ function(id) {
+ getNodeProperties(id, griwrm)$Diversion
+ })
+ )
+ if (length(upstreamDiversion) > 0) {
+ Qupstream[, upstreamDiversion] <- - Qupstream[, upstreamDiversion]
+ }
+ LengthHydro <- griwrm$length[UpstreamNodeRows]
+ names(LengthHydro) <- griwrm$id[UpstreamNodeRows]
BasinAreas <- c(
- griwrm$area[griwrm$id %in% UpstreamNodes],
- node$area - sum(griwrm$area[griwrm$id %in% UpstreamNodes], na.rm = TRUE)
+ griwrm$area[UpstreamNodeRows],
+ node$area - sum(griwrm$area[UpstreamNodeRows], na.rm = TRUE)
)
if (BasinAreas[length(BasinAreas)] < 0) {
stop(sprintf(
@@ -170,7 +284,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
id
))
}
- names(BasinAreas) <- c(UpstreamNodes, id)
+ names(BasinAreas) <- c(griwrm$id[UpstreamNodeRows], id)
}
# Set model inputs with the **airGR** function
@@ -185,9 +299,13 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
# Add Identifiers of connected nodes in order to be able to update them with simulated flows
InputsModel$id <- id
InputsModel$down <- node$down
- if(length(UpstreamNodes) > 0) {
- InputsModel$UpstreamNodes <- UpstreamNodes
- InputsModel$UpstreamIsRunoff <- !is.na(griwrm$model[match(UpstreamNodes, griwrm$id)])
+ if(length(UpstreamNodeRows) > 0) {
+ InputsModel$UpstreamNodes <- griwrm$id[UpstreamNodeRows]
+ InputsModel$UpstreamIsModeled <- !is.na(griwrm$model[UpstreamNodeRows])
+ InputsModel$UpstreamVarQ <- ifelse(!is.na(griwrm$model[UpstreamNodeRows]) &
+ griwrm$model[UpstreamNodeRows] == "Diversion",
+ "Qdiv_m3",
+ "Qsim_m3")
} else {
InputsModel$BasinAreas <- node$area
}
@@ -201,6 +319,12 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
InputsModel$model <- list(indexParamUngauged = ifelse(inherits(InputsModel, "SD"), 0, 1) + seq.int(featModel$NbParam),
hasX4 = grepl("RunModel_GR[456][HJ]", FUN_MOD),
iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4))
+ InputsModel$hasDiversion <- hasDiversion
+ if (hasDiversion) {
+ InputsModel$diversionOutlet <- diversionOutlet
+ InputsModel$Qdiv <- -Qobs[, id]
+ InputsModel$Qmin <- Qmin
+ }
return(InputsModel)
}
@@ -288,6 +412,7 @@ hasUngaugedNodes <- function(id, griwrm) {
#' function to extract model features partially copied from airGR:::.GetFeatModel
+#' @importFrom utils tail
#' @noRd
.GetFeatModel <- function(InputsModel) {
path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index 390aab4b82119585f071f477404597b4949cd640..48dfeca159a60f37504a9142d6577f60f1a61add 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -17,7 +17,7 @@
#'
#' @rdname CreateRunOptions
#' @export
-#' @inherit RunModel.GRiwrmInputsModel return examples
+#' @example man-examples/RunModel.GRiwrmInputsModel.R
CreateRunOptions <- function(x, ...) {
UseMethod("CreateRunOptions", x)
}
diff --git a/R/CreateSupervisor.R b/R/CreateSupervisor.R
index 5aec3fbeff0a48e79782d0dc6a93dee82718ccab..93b7fe29fabdba949f71e1b6eea648f319669707 100644
--- a/R/CreateSupervisor.R
+++ b/R/CreateSupervisor.R
@@ -11,28 +11,13 @@
#' - some internal state variables updated during simulation (`ts.index`, `ts.previous`, `ts.date`, `ts.index0`, `controller.id`)
#' @export
#'
-#' @examples
-#' data(Severn)
-#' nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
-#' nodes$model <- "RunModel_GR4J"
-#' griwrm <- CreateGRiwrm(nodes,
-#' list(id = "gauge_id",
-#' down = "downstream_id",
-#' length = "distance_downstream"))
-#' BasinsObs <- Severn$BasinsObs
-#' DatesR <- BasinsObs[[1]]$DatesR
-#' PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
-#' PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
-#' Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
-#' Precip <- ConvertMeteoSD(griwrm, PrecipTot)
-#' PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-#' InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
-#' sv <- CreateSupervisor(InputsModel)
+#' @example man-examples/RunModel.Supervisor.R
CreateSupervisor <- function(InputsModel, TimeStep = 1L) {
if(!inherits(InputsModel, "GRiwrmInputsModel")) {
stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
}
- if(!is.integer(TimeStep)) stop("`TimeStep` parameter must be an integer")
+ if (!is.integer(TimeStep)) stop("`TimeStep` parameter must be an integer")
+ if (TimeStep < 1) stop("`TimeStep` parameter must be strictly positive")
# Create Supervisor environment from the parent of GlobalEnv
e <- new.env(parent = parent.env(globalenv()))
@@ -48,6 +33,12 @@ CreateSupervisor <- function(InputsModel, TimeStep = 1L) {
e$DatesR <- InputsModel[[1]]$DatesR
e$InputsModel <- InputsModel
e$griwrm <- attr(InputsModel, "GRiwrm")
+ # Commands U are only applied on DirectInjection and Diversion
+ e$griwrm4U <- e$griwrm[is.na(e$griwrm$model) | e$griwrm$model == "Diversion", ]
+ e$nodeProperties <- lapply(e$griwrm$id[getDiversionRows(e$griwrm, TRUE)],
+ getNodeProperties,
+ griwrm = e$griwrm)
+ names(e$nodeProperties) <- e$griwrm$id[getDiversionRows(e$griwrm, TRUE)]
e$OutputsModel <- list()
e$.TimeStep <- TimeStep
@@ -74,3 +65,8 @@ CreateSupervisor <- function(InputsModel, TimeStep = 1L) {
return(e)
}
+
+
+is.Supervisor <- function(x) {
+ return(inherits(x, "Supervisor") && x$.isSupervisor == "3FJKmDcJ4snDbVBg")
+}
diff --git a/R/RunModel.GRiwrmInputsModel.R b/R/RunModel.GRiwrmInputsModel.R
index 439d6f2d40e0d4c5d2b9d2e29687b78bce1ed7ca..bd0f9833c69ca0be68aa9ad51006aa498ea41487 100644
--- a/R/RunModel.GRiwrmInputsModel.R
+++ b/R/RunModel.GRiwrmInputsModel.R
@@ -7,90 +7,7 @@
#'
#' @return An object of class \emph{GRiwrmOutputsModel}. This object is a [list] of *OutputsModel* objects produced by [RunModel.InputsModel] for each node of the semi-distributed model.
#' @export
-#' @examples
-#' ###################################################################
-#' # Run the `airGR::RunModel_Lag` example in the GRiwrm fashion way #
-#' # Simulation of a reservoir with a purpose of low-flow mitigation #
-#' ###################################################################
-#'
-#' ## ---- preparation of the InputsModel object
-#'
-#' ## loading package and catchment data
-#' library(airGRiwrm)
-#' data(L0123001)
-#'
-#' ## ---- specifications of the reservoir
-#'
-#' ## the reservoir withdraws 1 m3/s when it's possible considering the flow observed in the basin
-#' Qupstream <- matrix(-sapply(BasinObs$Qls / 1000 - 1, function(x) {
-#' min(1, max(0, x, na.rm = TRUE))
-#' }), ncol = 1)
-#'
-#' ## except between July and September when the reservoir releases 3 m3/s for low-flow mitigation
-#' month <- as.numeric(format(BasinObs$DatesR, "%m"))
-#' Qupstream[month >= 7 & month <= 9] <- 3
-#' Qupstream <- Qupstream * 86400 ## Conversion in m3/day
-#'
-#' ## the reservoir is not an upstream subcachment: its areas is NA
-#' BasinAreas <- c(NA, BasinInfo$BasinArea)
-#'
-#' ## delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
-#' LengthHydro <- 150
-#' ## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
-#' Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
-#'
-#' # This example is a network of 2 nodes which can be describe like this:
-#' db <- data.frame(id = c("Reservoir", "GaugingDown"),
-#' length = c(LengthHydro, NA),
-#' down = c("GaugingDown", NA),
-#' area = c(NA, BasinInfo$BasinArea),
-#' model = c(NA, "RunModel_GR4J"),
-#' stringsAsFactors = FALSE)
-#'
-#' # Create GRiwrm object from the data.frame
-#' griwrm <- CreateGRiwrm(db)
-#' str(griwrm)
-#'
-#' # Formatting observations for the hydrological models
-#' # Each input data should be a matrix or a data.frame with the good id in the name of the column
-#' Precip <- matrix(BasinObs$P, ncol = 1)
-#' colnames(Precip) <- "GaugingDown"
-#' PotEvap <- matrix(BasinObs$E, ncol = 1)
-#' colnames(PotEvap) <- "GaugingDown"
-#'
-#' # Observed flows contain flows that are directly injected in the model
-#' Qobs = matrix(Qupstream, ncol = 1)
-#' colnames(Qobs) <- "Reservoir"
-#'
-#' # Creation of the GRiwrmInputsModel object (= a named list of InputsModel objects)
-#' InputsModels <- CreateInputsModel(griwrm,
-#' DatesR = BasinObs$DatesR,
-#' Precip = Precip,
-#' PotEvap = PotEvap,
-#' Qobs = Qobs)
-#' str(InputsModels)
-#'
-#' ## run period selection
-#' Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
-#' which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
-#'
-#' # Creation of the GriwmRunOptions object
-#' RunOptions <- CreateRunOptions(InputsModels,
-#' IndPeriod_Run = Ind_Run)
-#' str(RunOptions)
-#'
-#' # Parameters of the SD models should be encapsulated in a named list
-#' ParamGR4J <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
-#' Param <- list(`GaugingDown` = c(Velocity, ParamGR4J))
-#'
-#' # RunModel for the whole network
-#' OutputsModels <- RunModel(InputsModels,
-#' RunOptions = RunOptions,
-#' Param = Param)
-#' str(OutputsModels)
-#'
-#' # Compare Simulation with reservoir and observation of natural flow
-#' plot(OutputsModels, data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]))
+#' @example man-examples/RunModel.GRiwrmInputsModel.R
RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
checkRunModelParameters(x, RunOptions, Param)
@@ -99,10 +16,10 @@ RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
for(id in names(x)) {
- message("RunModel.GRiwrmInputsModel: Treating sub-basin ", x[[id]]$id, "...")
+ message("RunModel.GRiwrmInputsModel: Processing sub-basin ", x[[id]]$id, "...")
# Update x[[id]]$Qupstream with simulated upstream flows
- if(any(x[[id]]$UpstreamIsRunoff)) {
+ if(any(x[[id]]$UpstreamIsModeled)) {
x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]], OutputsModel)
}
# Run the model for the sub-basin
diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index d64876a9886cc1fbbf754ac2c87f9ea1fa89c508..12727eb4adc6c6aa5dd404e4b654ca97368c970c 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -18,10 +18,76 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) {
OutputsModel <- airGR::RunModel(x, RunOptions, Param, FUN_MOD)
if (is.null(OutputsModel$Qsim_m3)) {
# Add Qsim_m3 in m3/timestep
- OutputsModel$Qsim_m3 <- OutputsModel$Qsim * sum(x$BasinAreas) * 1e3
+ OutputsModel$Qsim_m3 <-
+ OutputsModel$Qsim * sum(x$BasinAreas, na.rm = TRUE) * 1e3
}
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
- OutputsModel$RunOptions$WarmUpQsim_m3 <- OutputsModel$RunOptions$WarmUpQsim * sum(x$BasinAreas) * 1e3
+ OutputsModel$RunOptions$WarmUpQsim_m3 <-
+ OutputsModel$RunOptions$WarmUpQsim * sum(x$BasinAreas, na.rm = TRUE) * 1e3
+ }
+ if (x$hasDiversion) {
+ OutputsModel <- RunModel_Diversion(x, RunOptions, OutputsModel)
+ }
+ return(OutputsModel)
+}
+
+
+#' Model the diversion of a flow from an existing modeled node
+#'
+#' On a Diversion node, this function is called after `airGR::RunModel` to
+#' divert a part of the flow to another node than the original downstream one.
+#'
+#' @param InputsModel \[object of class \emph{InputsModel}\] see
+#' [airGR::CreateInputsModel] for details
+#' @param RunOptions Same parameter as in [RunModel.GRiwrmInputsModel]
+#' @param OutputsModel Output of [airGR::RunModel]
+#' @param updateQsim [logical] for updating Qsim after diversion in the output
+#'
+#' @return Updated `OutputsModel` object after diversion
+#' @noRd
+#'
+RunModel_Diversion <- function(InputsModel,
+ RunOptions,
+ OutputsModel,
+ updateQsim = TRUE) {
+ OutputsModel$Qnat <- OutputsModel$Qsim
+ lQ <- calc_Qdiv(OutputsModel$Qsim_m3,
+ InputsModel$Qdiv[RunOptions$IndPeriod_Run],
+ InputsModel$Qmin[RunOptions$IndPeriod_Run])
+ #message(paste(InputsModel$Qdiv[RunOptions$IndPeriod_Run], lQ$Qdiv, lQ$Qsim, InputsModel$Qmin[RunOptions$IndPeriod_Run], sep = ", "))
+ OutputsModel$Qdiv_m3 <- lQ$Qdiv
+ OutputsModel$Qsim_m3 <- lQ$Qsim
+ if (updateQsim) {
+ OutputsModel$Qsim <-
+ OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ }
+ if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+ lQ <- calc_Qdiv(OutputsModel$RunOptions$WarmUpQsim_m3,
+ InputsModel$Qdiv[RunOptions$IndPeriod_WarmUp],
+ InputsModel$Qmin[RunOptions$IndPeriod_WarmUp])
+ OutputsModel$RunOptions$WarmUpQdiv_m3 <- lQ$Qdiv
+ OutputsModel$RunOptions$WarmUpQsim_m3 <- lQ$Qsim
}
return(OutputsModel)
}
+
+
+#' Compute diverted and simulated flow at a diversion
+#'
+#' @param Qnat [numeric] time series of flow before diversion (m3/time step)
+#' @param Qdiv [numeric] time series of planned diverted flow (m3/time step)
+#' @param Qmin [numeric] time series of minimum flow after diversion (m3/time step)
+#'
+#' @return A [list] with items:
+#' - Qdiv, the diverted flow after limitation of minimum flow
+#' - Qsim, the simulated flow after diversion and limitation
+#' @noRd
+calc_Qdiv<- function(Qnat, Qdiv, Qmin) {
+ Qsim <- Qnat - Qdiv
+ indexQmin <- which(Qsim < Qmin & Qdiv > 0)
+ if (any(indexQmin)) {
+ #Qsim[indexQmin] <- sapply(indexQmin, function(i) min(Qnat[i], Qmin[i]))
+ Qsim[indexQmin] <- pmin(Qnat[indexQmin], Qmin[indexQmin])
+ }
+ return(list(Qsim = Qsim, Qdiv = Qnat - Qsim))
+}
diff --git a/R/RunModel.SD.R b/R/RunModel.SD.R
index 6d141d421caee7ba321f18cae3d28fca2e9f8c41..5498bcfa30a510832ac21f4f79d7d56673931668 100644
--- a/R/RunModel.SD.R
+++ b/R/RunModel.SD.R
@@ -8,5 +8,9 @@
#' @export
#'
RunModel.SD <- function(x, RunOptions, Param, QcontribDown, ...) {
- airGR::RunModel_Lag(x, RunOptions = RunOptions, Param = Param[1], QcontribDown = QcontribDown)
+ OutputsModel <- airGR::RunModel_Lag(x,
+ RunOptions = RunOptions,
+ Param = Param[1],
+ QcontribDown = QcontribDown)
+ return(OutputsModel)
}
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index 937dc40ea1f76bb9173aeea48ac04d8f43ea5501..3c9518fd7b6ac8a661c2235a5cab01fe9ca5fa29 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -7,8 +7,13 @@
#'
#' @return \emph{GRiwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See [airGR::RunModel]) for each node of the semi-distributed model
#' @export
+#'
+#' @example man-examples/RunModel.Supervisor.R
RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
+ stopifnot(is.Supervisor(x),
+ inherits(RunOptions, "GRiwrmRunOptions"))
+
# Save InputsModel for restoration at the end (Supervisor is an environment...)
InputsModelBackup <- x$InputsModel
@@ -28,32 +33,35 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
# Run runoff model for each sub-basin
x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
- RunModel.GR(IM,
- RunOptions = RunOptions[[IM$id]],
- Param = Param[[IM$id]])
+ OM_GR <- RunModel.GR(IM,
+ RunOptions = RunOptions[[IM$id]],
+ Param = Param[[IM$id]])
+ if (IM$hasDiversion) OM_GR$Qnat <- OM_GR$Qsim
+ return(OM_GR)
})
- class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
- # Copy simulated pure runoff flows (no SD nodes) to Qupstream in downstream SD nodes
- for(id in getNoSD_Ids(x$InputsModel)) {
- downId <- x$InputsModel[[id]]$down
- if(!is.null(x$InputsModel[[downId]])) {
- x$InputsModel[[downId]]$Qupstream[RunOptions[[downId]]$IndPeriod_Run, id] <-
- x$OutputsModel[[id]]$Qsim_m3
- }
+ class(x$OutputsModel) <- c("GRiwrmOutputsModel", class(x$OutputsModel))
+
+ # Copy simulated pure runoff flows (no SD nor Diversion nodes) to Qupstream
+ for(id in getNoSD_Ids(x$InputsModel, include_diversion = FALSE)) {
+ updateQupstream.Supervisor(x, id, IndPeriod_Run)
}
- # Save Qsim for step by step simulation
+ # Save OutputsModel for step by step simulation
QcontribDown <- do.call(
cbind,
lapply(x$OutputsModel, "[[", "Qsim")
)
-
Qsim_m3 <- do.call(
cbind,
lapply(x$OutputsModel, "[[", "Qsim_m3")
)
+ if (length(getDiversionRows(x$griwrm)) > 0) {
+ # Outputs of Diversion nodes
+ Qdiv_m3 <- Qsim_m3[, x$griwrm$id[getDiversionRows(x$griwrm)], drop = FALSE] * NA
+ Qnat <- Qdiv_m3
+ }
- # Initialisation of model states by running the model with no supervision on warm-up period
+ # Initialization of model states by running the model with no supervision on warm-up period
RunOptionsWarmUp <- RunOptions
for(id in names(x$InputsModel)) {
RunOptionsWarmUp[[id]]$IndPeriod_Run <- RunOptionsWarmUp[[id]]$IndPeriod_WarmUp
@@ -73,8 +81,16 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
x$OutputsModel[[id]]$StateEnd <- serializeIniStates(OM_WarmUp[[id]]$StateEnd)
}
+ message("Processing: 0%", appendLF = FALSE)
+ iProgressSteps <- round(length(lSuperTS) * seq(0.1, 0.9, 0.1))
+
# Loop over time steps with a step equal to the supervision time step
- for(iTS in lSuperTS) {
+ for(i in seq_along(lSuperTS)) {
+ iProgressMessage <- which(i == iProgressSteps)
+ if (length(iProgressMessage) == 1) {
+ message(" ", 10 * iProgressMessage, "%", appendLF = FALSE)
+ }
+ iTS <- lSuperTS[[i]]
# Run regulation on the whole basin for the current time step
x$ts.index <- iTS - x$ts.index0
x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
@@ -82,32 +98,49 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
if(iTS[1] > ts.start) {
doSupervision(x)
}
-
# Loop over sub-basin using SD model
- for(id in getSD_Ids(x$InputsModel)) {
- # Run the SD model for the sub-basin and one time step
- RunOptions[[id]]$IndPeriod_Run <- iTS
+ for(id in getSD_Ids(x$InputsModel, add_diversion = TRUE)) {
+ # Run model for the sub-basin and one time step
RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd)
- x$OutputsModel[[id]] <- RunModel.SD(
- x$InputsModel[[id]],
- RunOptions = RunOptions[[id]],
- Param = Param[[id]],
- QcontribDown = QcontribDown[x$ts.index, id]
- )
- # Storing Qsim in the data.frame Qsim
+ RunOptions[[id]]$IndPeriod_Run <- iTS
+ if (RunOptions[[id]]$FeatFUN_MOD$IsSD) {
+ # Route upstream flows for SD nodes
+ x$OutputsModel[[id]] <- RunModel.SD(
+ x$InputsModel[[id]],
+ RunOptions = RunOptions[[id]],
+ Param = Param[[id]],
+ QcontribDown = QcontribDown[x$ts.index, id]
+ )
+ } else {
+ x$OutputsModel[[id]]$Qsim_m3 <- Qsim_m3[x$ts.index, id]
+ }
+ if (x$InputsModel[[id]]$hasDiversion) {
+ # Compute diverted and simulated flows on Diversion nodes
+ x$OutputsModel[[id]] <-
+ RunModel_Diversion(x$InputsModel[[id]],
+ RunOptions = RunOptions[[id]],
+ OutputsModel = x$OutputsModel[[id]])
+ }
+ # Storing Qsim_m3 and Qdiv_m3 data.frames
Qsim_m3[x$ts.index, id] <- x$OutputsModel[[id]]$Qsim_m3
- # Routing Qsim to Qupstream of downstream nodes
- if(!is.na(x$InputsModel[[id]]$down)) {
- x$InputsModel[[x$InputsModel[[id]]$down]]$Qupstream[iTS, id] <-
- x$OutputsModel[[id]]$Qsim_m3
+ if (x$InputsModel[[id]]$hasDiversion) {
+ Qdiv_m3[x$ts.index, id] <- x$OutputsModel[[id]]$Qdiv_m3
}
+ # Routing Qsim_m3 and Qdiv_m3 to Qupstream of downstream nodes
+ updateQupstream.Supervisor(x, id, iTS)
}
x$ts.previous <- x$ts.index
}
+
+ message(" 100%")
+
for(id in getSD_Ids(x$InputsModel)) {
x$OutputsModel[[id]]$Qsim_m3 <- Qsim_m3[, id]
x$OutputsModel[[id]]$Qsim <-
Qsim_m3[, id] / sum(x$InputsModel[[id]]$BasinAreas, na.rm = TRUE) / 1e3
+ if (x$InputsModel[[id]]$hasDiversion) {
+ x$OutputsModel[[id]]$Qdiv_m3 <- Qdiv_m3[, id]
+ }
}
attr(x$OutputsModel, "Qm3s") <- OutputsModelQsim(x$InputsModel, x$OutputsModel, IndPeriod_Run)
@@ -116,3 +149,18 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
return(x$OutputsModel)
}
+
+updateQupstream.Supervisor <- function(x, id, iTS) {
+ downId <- x$InputsModel[[id]]$down
+ if(!is.null(x$InputsModel[[downId]])) {
+ x$InputsModel[[downId]]$Qupstream[iTS, id] <-
+ x$OutputsModel[[id]]$Qsim_m3
+ }
+ if (x$InputsModel[[id]]$hasDiversion) {
+ divOutId <- x$InputsModel[[id]]$diversionOutlet
+ if (!is.na(divOutId)) {
+ x$InputsModel[[divOutId]]$Qupstream[iTS, id] <-
+ x$OutputsModel[[id]]$Qdiv_m3
+ }
+ }
+}
diff --git a/R/UpdateQsimUpstream.R b/R/UpdateQsimUpstream.R
index 3afa52d8c94c02ddb7b716077973866eb1d404a6..9d0a012a2b4af6697babb3ddb838b9b2eb3eabd9 100644
--- a/R/UpdateQsimUpstream.R
+++ b/R/UpdateQsimUpstream.R
@@ -11,12 +11,13 @@
#' @noRd
#'
UpdateQsimUpstream <- function(InputsModel, Runoptions, OutputsModel) {
- iQ <- which(InputsModel$UpstreamIsRunoff)
+ iQ <- which(InputsModel$UpstreamIsModeled)
for(i in iQ) {
- InputsModel$Qupstream[Runoptions$IndPeriod_Run, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$Qsim_m3
- if (!is.null(OutputsModel[[InputsModel$UpstreamNodes[i]]]$RunOptions$WarmUpQsim_m3)) {
- InputsModel$Qupstream[Runoptions$IndPeriod_WarmUp, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$RunOptions$WarmUpQsim_m3
- }
+ InputsModel$Qupstream[Runoptions$IndPeriod_Run, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]][[InputsModel$UpstreamVarQ[i]]]
+ varWarmupQ_m3 <- paste0("WarmUp", InputsModel$UpstreamVarQ[i])
+ if (!is.null(OutputsModel[[InputsModel$UpstreamNodes[i]]]$RunOptions[[varWarmupQ_m3]])) {
+ InputsModel$Qupstream[Runoptions$IndPeriod_WarmUp, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]]$RunOptions[[varWarmupQ_m3]]
+ }
}
return(InputsModel)
}
diff --git a/R/plot.GRiwrm.R b/R/plot.GRiwrm.R
index 1698a13cd3d740783e063e7c8ed1baa354c38109..aa2ee083e7ff226b9230768e7e37cd4172d7b554 100644
--- a/R/plot.GRiwrm.R
+++ b/R/plot.GRiwrm.R
@@ -24,8 +24,8 @@ plot.GRiwrm <- function(x,
height = "100%",
box_colors = c(UpstreamUngauged = "#eef",
UpstreamGauged = "#aaf",
- IntermUngauged = "#efe",
- IntermGauged = "#afa",
+ IntermediateUngauged = "#efe",
+ IntermediateGauged = "#afa",
DirectInjection = "#faa"),
...) {
@@ -40,31 +40,53 @@ plot.GRiwrm <- function(x,
length(height) == 1,
is.character(box_colors),
length(setdiff(names(box_colors), c("UpstreamUngauged", "UpstreamGauged",
- "IntermUngauged", "IntermGauged",
+ "IntermediateUngauged", "IntermediateGauged",
"DirectInjection"))) == 0)
+ nodes <- sprintf("id_%1$s[%1$s]", x$id)
g2 <- x[!is.na(x$down),]
- nodes <- paste(
- sprintf("id_%1$s[%1$s]", g2$id),
+ links <- paste(
+ sprintf("id_%1$s", g2$id),
"-->|",
round(g2$length, digits = 0),
"km|",
- sprintf("id_%1$s[%1$s]", g2$down)
+ sprintf("id_%1$s", g2$down)
)
- node_class <- list(
- UpstreamUngauged = x$id[!x$id %in% x$down & x$model == "Ungauged"],
- UpstreamGauged = x$id[!x$id %in% x$down & x$model != "Ungauged" & !is.na(x$model)],
- IntermUngauged = x$id[x$id %in% x$down & x$model == "Ungauged"],
- IntermGauged = x$id[x$id %in% x$down & x$model != "Ungauged" & !is.na(x$model)],
- DirectInjection = x$id[is.na(x$model)]
- )
- node_class <- lapply(node_class, function(x) if(length(x) > 0) paste0("id_", x))
- node_class[sapply(node_class, is.null)] <- NULL
+ x$nodeclass <- sapply(x$id, getNodeClass, griwrm = x)
+ node_class <- lapply(unique(x$nodeclass), function(nc) {
+ x$id[x$nodeclass == nc]
+ })
+ names(node_class) <- unique(x$nodeclass)
+ node_class <- lapply(node_class, function(id) if(length(id) > 0) paste0("id_", id))
node_class <- paste("class", sapply(node_class, paste, collapse = ","), names(node_class))
- css <- paste("classDef", names(box_colors), paste0("fill:", box_colors))
- diagram <- paste(c(paste("graph", orientation), nodes, node_class, css), collapse = "\n\n")
+ css <- c(
+ paste("classDef", names(box_colors), paste0("fill:", box_colors)),
+ paste("classDef",
+ paste0(names(box_colors), "Diversion"),
+ sprintf("fill:%s, stroke:%s, stroke-width:3px", box_colors, box_colors["DirectInjection"]))
+ )
+ if (length(getDiversionRows(g2)) > 0) {
+ css <- c(css,
+ paste("linkStyle",
+ getDiversionRows(g2) - 1,
+ sprintf("stroke:%s, stroke-width:2px,stroke-dasharray: 5 5;",
+ box_colors["DirectInjection"])))
+ }
+ diagram <- paste(c(paste("graph", orientation), nodes, links, node_class, css),
+ collapse = "\n\n")
if (display) {
DiagrammeR::mermaid(diagram = diagram, width, height, ...)
} else {
return(diagram)
}
}
+
+getNodeClass <- function(id, griwrm) {
+ props <- getNodeProperties(id, griwrm)
+ if (props$DirectInjection) {
+ nc <- "DirectInjection"
+ } else {
+ nc <- paste0(props["position"], props["hydrology"])
+ }
+ if (props$Diversion) nc <- paste0(nc, "Diversion")
+ return(nc)
+}
diff --git a/R/plot.GRiwrmOutputsModel.R b/R/plot.GRiwrmOutputsModel.R
index 01d12f402d12bd0e91210b0dddc7821aeb0792ce..05081e5215e0e40acfa08bd9b144ddf44cdc57d9 100644
--- a/R/plot.GRiwrmOutputsModel.R
+++ b/R/plot.GRiwrmOutputsModel.R
@@ -11,7 +11,7 @@
#' @importFrom graphics plot par title
#' @export
#'
-#' @inherit RunModel.GRiwrmInputsModel return examples
+#' @example man-examples/RunModel.GRiwrmInputsModel.R
#'
plot.GRiwrmOutputsModel <- function(x, Qobs = NULL, ...) {
diff --git a/R/plot.Qm3s.R b/R/plot.Qm3s.R
index 65202859f2aefb4b4eb1ff9b1ab2ea450b3d095c..0b378cacdbfde6546b236de4378ddff7c77cc1d1 100644
--- a/R/plot.Qm3s.R
+++ b/R/plot.Qm3s.R
@@ -1,12 +1,18 @@
#' Plot of a `Qm3s` object (time series of simulated flows)
#'
-#' @param x [data.frame] with a first column with [POSIXt] dates and followings columns with flows at each node of the network
+#' This function plot time series of flow rate in m3/s. It's a method for object
+#' of class "Qm3s" which can be directly called by `plot`. It can also be called
+#' as a function `plot.Qm3s` if the first parameter has the good format.
+#'
+#' @param x [data.frame] with a first column with [POSIXt] dates and followings
+#' columns with flows at each node of the network
#' @param type [character] plot type (See [plot.default]), default "l"
#' @param xlab [character] label for the x axis, default to "Date"
#' @param ylab [character] label for the y axis, default to "Flow (m3/s)"
#' @param main [character] main title for the plot, default to "Simulated flows"
#' @param col [character] plotting color (See [par]), default to rainbow colors
-#' @param legend [character] see parameter `legend` of [legend]. Set to [NULL] if display of the legend is not wanted
+#' @param legend [character] see parameter `legend` of [legend]. Set it to [NULL]
+#' to hide the legend
#' @param legend.cex [character] `cex` parameter for the text of the legend (See [par])
#' @param legend.x,legend.y Legend position, see `x` and `y` parameters in [graphics::legend]
#' @param lty [character] or [numeric] The line type (See [par])
@@ -14,14 +20,17 @@
#'
#' @return Screen plot window.
#'
+#' @example man-examples/RunModel.GRiwrmInputsModel.R
+#'
#' @importFrom grDevices rainbow
#' @importFrom graphics matplot
+#' @export plot.Qm3s
#' @export
#'
plot.Qm3s <- function(x,
type = "l",
xlab = "Date",
- ylab = expression("Flow (m"^"3"*"/s)"),
+ ylab = expression("Flow rate (m"^"3"*"/s)"),
main = "Simulated flows",
col = rainbow(ncol(x) - 1),
legend = colnames(x)[-1],
@@ -30,8 +39,12 @@ plot.Qm3s <- function(x,
legend.y = NULL,
lty = 1,
...) {
+
+ stopifnot(is.data.frame(x),
+ inherits(x[, 1], "POSIXct"))
+
matplot(
- x$DatesR,
+ x[, 1],
x[, -1],
type = type,
lty = lty,
diff --git a/R/utils.R b/R/utils.R
index 2d641f622c78ae5bbbce745b72533829f863647b..a36aa1b1e22737806998118047014c3e9bc6ab33 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,15 +1,16 @@
#' Function to obtain the ID of sub-basins using SD model
#'
#' @param InputsModel \[`GRiwrmInputsModel` object\]
+#' @param add_diversion [logical] for adding upstream nodes with diversion
#'
#' @return [character] IDs of the sub-basins using SD model
#' @export
-getSD_Ids <- function(InputsModel) {
+getSD_Ids <- function(InputsModel, add_diversion = FALSE) {
if (!inherits(InputsModel, "GRiwrmInputsModel")) {
stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
}
bSDs <- sapply(InputsModel, function (IM) {
- inherits(IM, "SD")
+ inherits(IM, "SD") | IM$hasDiversion
})
names(InputsModel)[bSDs]
}
@@ -17,15 +18,16 @@ getSD_Ids <- function(InputsModel) {
#' Function to obtain the ID of sub-basins not using SD model
#'
#' @param InputsModel \[`GRiwrmInputsModel` object\]
+#' @param include_diversion [logical] for including diversion nodes
#'
#' @return [character] IDs of the sub-basins not using the SD model
#' @export
-getNoSD_Ids <- function(InputsModel) {
+getNoSD_Ids <- function(InputsModel, include_diversion = TRUE) {
if (!inherits(InputsModel, "GRiwrmInputsModel")) {
stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
}
bSDs <- sapply(InputsModel, function (IM) {
- !inherits(IM, "SD")
+ !inherits(IM, "SD") & (include_diversion | !IM$hasDiversion)
})
names(InputsModel)[bSDs]
}
@@ -44,11 +46,14 @@ getDataFromLocation <- function(loc, sv) {
if (length(grep("\\[[0-9]+\\]$", loc)) > 0) {
stop("Reaching output of other controller is not implemented yet")
} else {
- node <- sv$griwrm$down[sv$griwrm$id == loc]
- if(is.na(node)) {
- # Downstream node: simulated flow at last supervision time step (bug #40)
- sv$OutputsModel[[loc]]$Qsim_m3
+ if(sv$nodeProperties[[loc]]["hydrology"] != "DirectInjection") {
+ if (sv$nodeProperties[[loc]]$Upstream) {
+ sv$OutputsModel[[loc]]$Qsim_m3[sv$ts.previous]
+ } else {
+ sv$OutputsModel[[loc]]$Qsim_m3
+ }
} else {
+ node <- sv$griwrm$down[sv$griwrm$id == loc]
sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.previous, loc]
}
}
@@ -64,13 +69,22 @@ getDataFromLocation <- function(loc, sv) {
#' @noRd
setDataToLocation <- function(ctrlr, sv) {
l <- lapply(seq(length(ctrlr$Unames)), function(i) {
- node <- sv$griwrm$down[sv$griwrm$id == ctrlr$Unames[i]]
# limit U size to the number of simulation time steps of the current supervision time step
U <- ctrlr$U[seq.int(length(sv$ts.index)),i]
- # ! Qupstream contains warm up period and run period => the index is shifted
- if(!is.null(sv$InputsModel[[node]])) {
- sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index,
- ctrlr$Unames[i]] <- U
+
+ locU <- ctrlr$Unames[i]
+ if (sv$nodeProperties[[locU]]$DirectInjection) {
+ # Direct injection node => update Qusptream of downstream node
+ node <- sv$griwrm4U$down[sv$griwrm4U$id == locU]
+ # ! Qupstream contains warm up period and run period => the index is shifted
+ if(!is.null(sv$InputsModel[[node]])) {
+ sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, locU] <- U
+ }
+ } else if (sv$nodeProperties[[locU]]$Diversion){
+ # Diversion node => update Qdiv with -U
+ sv$InputsModel[[locU]]$Qdiv[sv$ts.index0 + sv$ts.index] <- -U
+ } else {
+ stop("Node ", locU, " must be a Direct Injection or a Diversion node")
}
})
}
@@ -94,6 +108,21 @@ doSupervision <- function(supervisor) {
if(is.vector(supervisor$controllers[[id]]$U)) {
supervisor$controllers[[id]]$U <- matrix(supervisor$controllers[[id]]$U, nrow = 1)
}
+ # Check U output
+ if (
+ ncol(supervisor$controllers[[id]]$U) != length(supervisor$controllers[[id]]$Unames) |
+ (!nrow(supervisor$controllers[[id]]$U) %in% c(supervisor$.TimeStep, length(supervisor$ts.index)))
+ ) {
+ stop("The logic function of the controller ",
+ supervisor$controllers[[id]]$name,
+ " should return a matrix of dimension ",
+ supervisor$.TimeStep, ", ", length(supervisor$controllers[[id]]$Unames))
+ }
+ # For the last supervisor time step which can be truncated
+ if (length(supervisor$ts.index) < supervisor$.TimeStep) {
+ supervisor$controllers[[id]]$U <-
+ supervisor$controllers[[id]]$U[seq(length(supervisor$ts.index)), ]
+ }
# Write U to locations in the model
setDataToLocation(supervisor$controllers[[id]], sv = supervisor)
}
@@ -133,23 +162,22 @@ checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
griwrm <- attr(InputsModel, "GRiwrm")
# Get simulated flow for each node
- # Flow for each node is available in InputsModel$Qupstream except for the downstream node
- upperNodes <- griwrm$id[!is.na(griwrm$down)]
+ # Flow for each node is available in OutputsModel except for Direct Injection
+ # nodes where it is stored in InputsModel$Qupstream of the downstream node
+ QsimRows <- getDiversionRows(griwrm, TRUE)
lQsim <- lapply(
- upperNodes,
- function(x, griwrm, IndPeriod_Run) {
- node <- griwrm$down[griwrm$id == x]
- InputsModel[[node]]$Qupstream[IndPeriod_Run, x]
- },
- griwrm = griwrm, IndPeriod_Run = IndPeriod_Run
+ QsimRows,
+ function(i) {
+ x <- griwrm[i, ]
+ if (is.na(x$model)) {
+ InputsModel[[x$down]]$Qupstream[IndPeriod_Run, x$id]
+ } else {
+ OutputsModel[[x$id]]$Qsim_m3
+ }
+ }
)
- names(lQsim) <- upperNodes
- # Flow of the downstream node is only available in OutputsModel[[node]]$Qsim
- downNode <- names(InputsModel)[length(InputsModel)]
- lQsim[[downNode]] <- OutputsModel[[downNode]]$Qsim_m3
-
- names(lQsim) <- c(upperNodes, downNode)
- dfQsim <- cbind(data.frame(DatesR = as.POSIXct(InputsModel[[1]]$DatesR[IndPeriod_Run])),
+ names(lQsim) <- griwrm$id[QsimRows]
+ dfQsim <- cbind(data.frame(DatesR = InputsModel[[1]]$DatesR[IndPeriod_Run]),
do.call(cbind,lQsim) / attr(InputsModel, "TimeStep"))
class(dfQsim) <- c("Qm3s", class(dfQsim)) # For S3 methods
return(dfQsim)
diff --git a/data-raw/Severn.R b/data-raw/Severn.R
index 5535ff91fc05663950a9fc17249e574eb632bdab..db720be1031c9481c2b04cb91611cba091d76c42 100644
--- a/data-raw/Severn.R
+++ b/data-raw/Severn.R
@@ -21,7 +21,7 @@ read_CAMELS_ts <- function(file, cut) {
df <- read.csv(file)
iCut <- which(df$date == cut)
df <- df[iCut:nrow(df),]
- df$DatesR <- as.POSIXct(df$date)
+ df$DatesR <- as.POSIXct(df$date, tz = "UTC")
return(df[,c("DatesR", "precipitation", "peti", "discharge_spec")])
}
diff --git a/data/Severn.rda b/data/Severn.rda
index 1d26cfd47d14b9716071c14f9346a4c05d61fbf5..1a6827656f3f7eeb2c655f4e5cc3cd5fdcab0aeb 100644
Binary files a/data/Severn.rda and b/data/Severn.rda differ
diff --git a/dev/pkgdown_build_site.R b/dev/pkgdown_build_site.R
new file mode 100644
index 0000000000000000000000000000000000000000..fdb3c170ad3ad204c84ecb3d5b8560dcca1bf99f
--- /dev/null
+++ b/dev/pkgdown_build_site.R
@@ -0,0 +1,2 @@
+remotes::install_gitlab("in-wop/seinebasin", host = "gitlab.irstea.fr")
+pkgdown::build_site("..")
diff --git a/man-examples/CreateGRiwrm.R b/man-examples/CreateGRiwrm.R
index 9ef97c074a3acb9b22aba2e8b03e361548e7b3a4..718deb6b1366893691883b1aa3635b43e3ccdc97 100644
--- a/man-examples/CreateGRiwrm.R
+++ b/man-examples/CreateGRiwrm.R
@@ -1,4 +1,6 @@
-# Network of 2 nodes distant of 150 km:
+#########################################
+# Network of 2 nodes distant of 150 km: #
+#########################################
# - an upstream reservoir modelled as a direct flow injection (no model)
# - a gauging station downstream a catchment of 360 km² modelled with GR4J
db <- data.frame(id = c("Reservoir", "GaugingDown"),
@@ -12,7 +14,9 @@ griwrm_basic
# Network diagram with direct flow node in red, intermediate sub-basin in green
plot(griwrm_basic)
-# GR4J semi-distributed model of the Severn River
+###################################################
+# GR4J semi-distributed model of the Severn River #
+###################################################
data(Severn)
nodes <- Severn$BasinsInfo
nodes$model <- "RunModel_GR4J"
@@ -26,7 +30,9 @@ griwrm_severn
# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
plot(griwrm_severn)
-# Same model with an ungauged station at nodes 54029 and 54001
+####################################################################
+# Severn network with an ungauged station at nodes 54029 and 54001 #
+####################################################################
# By default the first gauged node at downstream is used for parameter calibration (54032)
nodes_ungauged <- nodes
nodes_ungauged$model[nodes_ungauged$gauge_id %in% c("54029", "54001")] <- "Ungauged"
@@ -35,3 +41,17 @@ griwrm_ungauged <- CreateGRiwrm(nodes_ungauged, rename_columns)
griwrm_ungauged
# Network diagram with gauged nodes of vivid color, and ungauged nodes of dull color
plot(griwrm_ungauged)
+
+#######################################################
+# Severn network with a Diversion on the node "54029" #
+# which transfer flows to the node "54001" #
+#######################################################
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "model", "area")]
+nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54029",
+ downstream_id = "54001",
+ distance_downstream = 20,
+ model = "Diversion",
+ area = NA))
+griwrm_div <- CreateGRiwrm(nodes_div, rename_columns)
+# Network diagram figures Diversion node by a red frame and a red arrow
+plot(griwrm_div)
diff --git a/man-examples/RunModel.GRiwrmInputsModel.R b/man-examples/RunModel.GRiwrmInputsModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..c84a96a8048755b34dc6a9e936fef0184114e70b
--- /dev/null
+++ b/man-examples/RunModel.GRiwrmInputsModel.R
@@ -0,0 +1,210 @@
+###################################################################
+# Run the `airGR::RunModel_Lag` example in the GRiwrm fashion way #
+# Simulation of a reservoir with a purpose of low-flow mitigation #
+###################################################################
+
+## ---- preparation of the InputsModel object
+
+## loading package and catchment data
+library(airGRiwrm)
+data(L0123001)
+
+## ---- specifications of the reservoir
+
+## the reservoir withdraws 1 m3/s when it's possible considering the flow observed in the basin
+Qupstream <- matrix(-sapply(BasinObs$Qls / 1000 - 1, function(x) {
+ min(1, max(0, x, na.rm = TRUE))
+}), ncol = 1)
+
+## except between July and September when the reservoir releases 3 m3/s for low-flow mitigation
+month <- as.numeric(format(BasinObs$DatesR, "%m"))
+Qupstream[month >= 7 & month <= 9] <- 3
+Qupstream <- Qupstream * 86400 ## Conversion in m3/day
+
+## the reservoir is not an upstream subcachment: its areas is NA
+BasinAreas <- c(NA, BasinInfo$BasinArea)
+
+## delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
+LengthHydro <- 150
+## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
+Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
+
+# This example is a network of 2 nodes which can be describe like this:
+db <- data.frame(id = c("Reservoir", "GaugingDown"),
+ length = c(LengthHydro, NA),
+ down = c("GaugingDown", NA),
+ area = c(NA, BasinInfo$BasinArea),
+ model = c(NA, "RunModel_GR4J"),
+ stringsAsFactors = FALSE)
+
+# Create GRiwrm object from the data.frame
+griwrm <- CreateGRiwrm(db)
+str(griwrm)
+
+# Formatting observations for the hydrological models
+# Each input data should be a matrix or a data.frame with the good id in the name of the column
+Precip <- matrix(BasinObs$P, ncol = 1)
+colnames(Precip) <- "GaugingDown"
+PotEvap <- matrix(BasinObs$E, ncol = 1)
+colnames(PotEvap) <- "GaugingDown"
+
+# Observed flows contain flows that are directly injected in the model
+Qobs = matrix(Qupstream, ncol = 1)
+colnames(Qobs) <- "Reservoir"
+
+# Creation of the GRiwrmInputsModel object (= a named list of InputsModel objects)
+InputsModels <- CreateInputsModel(griwrm,
+ DatesR = BasinObs$DatesR,
+ Precip = Precip,
+ PotEvap = PotEvap,
+ Qobs = Qobs)
+str(InputsModels)
+
+## run period selection
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
+ which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
+
+# Creation of the GRiwmRunOptions object
+RunOptions <- CreateRunOptions(InputsModels,
+ IndPeriod_Run = Ind_Run)
+str(RunOptions)
+
+# Parameters of the SD models should be encapsulated in a named list
+ParamGR4J <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+Param <- list(`GaugingDown` = c(Velocity, ParamGR4J))
+
+# RunModel for the whole network
+OutputsModels <- RunModel(InputsModels,
+ RunOptions = RunOptions,
+ Param = Param)
+str(OutputsModels)
+
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
diff --git a/man-examples/RunModel.Supervisor.R b/man-examples/RunModel.Supervisor.R
new file mode 100644
index 0000000000000000000000000000000000000000..6eacf4f8e91e432492ae70a762af9265f93908ea
--- /dev/null
+++ b/man-examples/RunModel.Supervisor.R
@@ -0,0 +1,124 @@
+###############################################################################
+# An example of reservoir management on an hypothetical dam at station "54095"
+# on the Severn river build to support low-flows at "54057"
+###############################################################################
+# A minimum flow of 20 m3/s is maintained at the dam location and an extra-release
+# is provided when the flow at the downstream station "54057" cross a minimum
+# threshold of 45 m3/s. The dam has a storage capacity of 60 millions m3
+###############################################################################
+library(airGRiwrm)
+
+# Load Severn network information
+data(Severn)
+nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$model <- "RunModel_GR4J"
+
+# Insert a dam downstream the location the gauging station 54095
+# The dam is a direct injection node
+nodes$downstream_id[nodes$gauge_id == "54095"] <- "Dam"
+nodes$distance_downstream[nodes$gauge_id == "54095"] <- 0
+nodes <- rbind(nodes,
+ data.frame(gauge_id = "Dam",
+ downstream_id = "54001",
+ distance_downstream = 42,
+ area = NA,
+ model = NA))
+
+griwrm <- CreateGRiwrm(nodes,
+ list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream"))
+plot(griwrm)
+
+# Format meteorological inputs for CreateInputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+
+# Create a release flow time series for the dam
+# This release will be modified by the Supervisor
+# We initiate it with the natural flow for having a good initialisation of the
+# model at the first time step of the running period
+Qobs <- data.frame(
+ Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
+)
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
+
+# InputsModel object
+IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# Initialisation of the Supervisor
+sv <- CreateSupervisor(IM_severn)
+
+# States of the reservoir are stored in the Supervisor variable
+# The Supervisor variable is an environment which can be available in
+# the controller function for storing and exchange data during the simulation
+sv$Vres <- 0 # Reservoir storage time series
+
+# Dam management is modeled by a controller
+# This controller usually releases Qmin and provides
+# extra release if flow mesured somewhere is below Qthreshold
+# Flow is expressed in m3 / time step
+# Y[1] = runoff flow at gauging station 54095 filling the reservoir
+# Y[2] = flow at gauging station 54057, location of the low-flow objective
+# The returned value is the release calculated at the reservoir
+# We need to enclose the Supervisor variable and other parameters in
+# the environment of the function with a function returning the logic function
+factoryDamLogic <- function(sv, Vmin, Vmax, Qmin, Qthreshold) {
+ function(Y) {
+ # Filling of the reservoir
+ sv$Vres <- c(sv$Vres, tail(sv$Vres, 1) + Y[1])
+ # The release is the max between: overflow, low-flow support and minimum flow
+ U <- U <- max(tail(sv$Vres, 1) - Vmax, Qthreshold - Y[2], Qmin)
+ sv$Vres[length(sv$Vres)] <- tail(sv$Vres, 1) - U
+ if (tail(sv$Vres, 1) < Vmin) {
+ # Reservoir is empty
+ U <- U - (Vmin - tail(sv$Vres, 1))
+ sv$Vres[length(sv$Vres)] <- Vmin
+ }
+ return(U)
+ }
+}
+
+# And define a final function enclosing logic and parameters together
+funDamLogic <- factoryDamLogic(
+ sv = sv, # The Supervisor which store the states of reservoir storage
+ Vmin = 0, # Minimum volume in the reservoir (m3)
+ Vmax = 60 * 1E6, # Maximum volume in the reservoir (m3)
+ Qmin = 20 * 86400, # Min flow to maintain downstream the reservoir (m3/day)
+ Qthreshold = 45 * 86400 # Min flow threshold to support at station 54057 (m3/day)
+)
+
+CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN = funDamLogic)
+
+# GRiwrmRunOptions object simulation of the hydrological year 2002-2003
+IndPeriod_Run <- which(
+ DatesR >= as.POSIXct("2002-10-15", tz = "UTC") &
+ DatesR <= as.POSIXct("2003-10-15", tz = "UTC")
+)
+IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# The Supervisor is used instead of InputsModel for running the model
+OM_dam <- RunModel(sv,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plotting the time series of flows and reservoir storage
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
+ ylim = c(0, 200))
+Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
+ V_reservoir = sv$Vres)
+plot.Qm3s(Vres)
+par(oldpar)
diff --git a/man/CreateController.Rd b/man/CreateController.Rd
index 75eb2d4ed9bea474d55eee25a575e7aec291e558..93c1021881deccf5b70c0c29573c50d042492f12 100644
--- a/man/CreateController.Rd
+++ b/man/CreateController.Rd
@@ -43,26 +43,128 @@ for the previous time step and returns calculated \code{U}. These \code{U} will
at their location for the current time step of calculation of the model.
}
\examples{
-# First create a Supervisor from a model
+###############################################################################
+# An example of reservoir management on an hypothetical dam at station "54095"
+# on the Severn river build to support low-flows at "54057"
+###############################################################################
+# A minimum flow of 20 m3/s is maintained at the dam location and an extra-release
+# is provided when the flow at the downstream station "54057" cross a minimum
+# threshold of 45 m3/s. The dam has a storage capacity of 60 millions m3
+###############################################################################
+library(airGRiwrm)
+
+# Load Severn network information
data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$model <- "RunModel_GR4J"
+
+# Insert a dam downstream the location the gauging station 54095
+# The dam is a direct injection node
+nodes$downstream_id[nodes$gauge_id == "54095"] <- "Dam"
+nodes$distance_downstream[nodes$gauge_id == "54095"] <- 0
+nodes <- rbind(nodes,
+ data.frame(gauge_id = "Dam",
+ downstream_id = "54001",
+ distance_downstream = 42,
+ area = NA,
+ model = NA))
+
griwrm <- CreateGRiwrm(nodes,
list(id = "gauge_id",
down = "downstream_id",
length = "distance_downstream"))
+plot(griwrm)
+
+# Format meteorological inputs for CreateInputs
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
-Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
-sv <- CreateSupervisor(InputsModel)
-# A controller which usually releases 0.1 m3/s and provides
-# extra release if the downstream flow is below 0.5 m3/s
-logicDamRelease <- function(Y) max(0.5 - Y[1], 0.1)
-CreateController(sv, "DamRelease", Y = c("54001"), U = c("54095"), FUN = logicDamRelease)
+# Create a release flow time series for the dam
+# This release will be modified by the Supervisor
+# We initiate it with the natural flow for having a good initialisation of the
+# model at the first time step of the running period
+Qobs <- data.frame(
+ Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
+)
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
+
+# InputsModel object
+IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# Initialisation of the Supervisor
+sv <- CreateSupervisor(IM_severn)
+
+# States of the reservoir are stored in the Supervisor variable
+# The Supervisor variable is an environment which can be available in
+# the controller function for storing and exchange data during the simulation
+sv$Vres <- 0 # Reservoir storage time series
+
+# Dam management is modeled by a controller
+# This controller usually releases Qmin and provides
+# extra release if flow mesured somewhere is below Qthreshold
+# Flow is expressed in m3 / time step
+# Y[1] = runoff flow at gauging station 54095 filling the reservoir
+# Y[2] = flow at gauging station 54057, location of the low-flow objective
+# The returned value is the release calculated at the reservoir
+# We need to enclose the Supervisor variable and other parameters in
+# the environment of the function with a function returning the logic function
+factoryDamLogic <- function(sv, Vmin, Vmax, Qmin, Qthreshold) {
+ function(Y) {
+ # Filling of the reservoir
+ sv$Vres <- c(sv$Vres, tail(sv$Vres, 1) + Y[1])
+ # The release is the max between: overflow, low-flow support and minimum flow
+ U <- U <- max(tail(sv$Vres, 1) - Vmax, Qthreshold - Y[2], Qmin)
+ sv$Vres[length(sv$Vres)] <- tail(sv$Vres, 1) - U
+ if (tail(sv$Vres, 1) < Vmin) {
+ # Reservoir is empty
+ U <- U - (Vmin - tail(sv$Vres, 1))
+ sv$Vres[length(sv$Vres)] <- Vmin
+ }
+ return(U)
+ }
+}
+
+# And define a final function enclosing logic and parameters together
+funDamLogic <- factoryDamLogic(
+ sv = sv, # The Supervisor which store the states of reservoir storage
+ Vmin = 0, # Minimum volume in the reservoir (m3)
+ Vmax = 60 * 1E6, # Maximum volume in the reservoir (m3)
+ Qmin = 20 * 86400, # Min flow to maintain downstream the reservoir (m3/day)
+ Qthreshold = 45 * 86400 # Min flow threshold to support at station 54057 (m3/day)
+)
+
+CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN = funDamLogic)
+
+# GRiwrmRunOptions object simulation of the hydrological year 2002-2003
+IndPeriod_Run <- which(
+ DatesR >= as.POSIXct("2002-10-15", tz = "UTC") &
+ DatesR <= as.POSIXct("2003-10-15", tz = "UTC")
+)
+IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# The Supervisor is used instead of InputsModel for running the model
+OM_dam <- RunModel(sv,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plotting the time series of flows and reservoir storage
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
+ ylim = c(0, 200))
+Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
+ V_reservoir = sv$Vres)
+plot.Qm3s(Vres)
+par(oldpar)
}
diff --git a/man/CreateGRiwrm.Rd b/man/CreateGRiwrm.Rd
index 24d40194789be032f734d5425a6b9c16cc052407..f2ee4db078b466de81be4b225cb94036b87da3b3 100644
--- a/man/CreateGRiwrm.Rd
+++ b/man/CreateGRiwrm.Rd
@@ -62,16 +62,22 @@ The "model" column should be filled by one of the following:
\itemize{
\item One of the hydrological models available in the \emph{airGR} package defined by its
\code{RunModel} function (i.e.: \code{RunModel_GR4J}, \code{RunModel_GR5HCemaneige}...)
-\item \code{NA} for injecting (or abstracting) a flow time series at the location of the node
-(direct flow injection)
\item \code{Ungauged} for an ungauged node. The sub-basin inherits hydrological model and
parameters from a "donor" sub-basin. By default the donor is the first gauged
node at downstream
+\item \code{NA} for injecting (or abstracting) a flow time series at the location of the node
+(direct flow injection)
+\item \code{Diversion} for abstracting a flow time series from an existing node transfer it
+to another node. As a \code{Diversion} is attached to an existing node, this node is
+then described with 2 lines: one for the hydrological model and another one for the
+diversion
}
}
}
\examples{
-# Network of 2 nodes distant of 150 km:
+#########################################
+# Network of 2 nodes distant of 150 km: #
+#########################################
# - an upstream reservoir modelled as a direct flow injection (no model)
# - a gauging station downstream a catchment of 360 km² modelled with GR4J
db <- data.frame(id = c("Reservoir", "GaugingDown"),
@@ -85,7 +91,9 @@ griwrm_basic
# Network diagram with direct flow node in red, intermediate sub-basin in green
plot(griwrm_basic)
-# GR4J semi-distributed model of the Severn River
+###################################################
+# GR4J semi-distributed model of the Severn River #
+###################################################
data(Severn)
nodes <- Severn$BasinsInfo
nodes$model <- "RunModel_GR4J"
@@ -99,7 +107,9 @@ griwrm_severn
# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
plot(griwrm_severn)
-# Same model with an ungauged station at nodes 54029 and 54001
+####################################################################
+# Severn network with an ungauged station at nodes 54029 and 54001 #
+####################################################################
# By default the first gauged node at downstream is used for parameter calibration (54032)
nodes_ungauged <- nodes
nodes_ungauged$model[nodes_ungauged$gauge_id \%in\% c("54029", "54001")] <- "Ungauged"
@@ -108,4 +118,18 @@ griwrm_ungauged <- CreateGRiwrm(nodes_ungauged, rename_columns)
griwrm_ungauged
# Network diagram with gauged nodes of vivid color, and ungauged nodes of dull color
plot(griwrm_ungauged)
+
+#######################################################
+# Severn network with a Diversion on the node "54029" #
+# which transfer flows to the node "54001" #
+#######################################################
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "model", "area")]
+nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54029",
+ downstream_id = "54001",
+ distance_downstream = 20,
+ model = "Diversion",
+ area = NA))
+griwrm_div <- CreateGRiwrm(nodes_div, rename_columns)
+# Network diagram figures Diversion node by a red frame and a red arrow
+plot(griwrm_div)
}
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index 29aaf5501b2ac4fc76b713fae301ab854f271857..9a438b0061505b4c1de25efb53f0d3191ffb97bc 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -10,6 +10,7 @@
Precip = NULL,
PotEvap = NULL,
Qobs = NULL,
+ Qmin = NULL,
PrecipScale = TRUE,
TempMean = NULL,
TempMin = NULL,
@@ -25,38 +26,79 @@
\item{DatesR}{\link{POSIXt} vector of dates}
-\item{Precip}{(optional) \link{matrix} or \link{data.frame} frame of numeric containing precipitation in [mm per time step]. Column names correspond to node IDs}
+\item{Precip}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+precipitation in [mm per time step]. Column names correspond to node IDs}
-\item{PotEvap}{(optional) \link{matrix} or \link{data.frame} frame of numeric containing potential evaporation [mm per time step]. Column names correspond to node IDs}
+\item{PotEvap}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+potential evaporation [mm per time step]. Column names correspond to node IDs}
-\item{Qobs}{(optional) \link{matrix} or \link{data.frame} frame of numeric containing observed flows in [mm per time step]. Column names correspond to node IDs}
+\item{Qobs}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+observed flows. It must be provided only for nodes of type "Direct
+injection" and "Diversion". See \link{CreateGRiwrm} for
+details about these node types. Unit is [mm per time step] for nodes
+with an area, and [m3 per time step] for nodes with \code{area=NA}.
+Column names correspond to node IDs. Negative flows are abstracted from
+the model and positive flows are injected to the model}
-\item{PrecipScale}{(optional) named \link{vector} of \link{logical} indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default \code{TRUE} (the mean of the precipitation is kept to the original value)}
+\item{Qmin}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+minimum flows to let downstream of a node with a Diversion [m3 per
+time step]. Default is zero. Column names correspond to node IDs}
-\item{TempMean}{(optional) \link{matrix} or \link{data.frame} of time series of mean air temperature [°C], required to create the CemaNeige module inputs}
+\item{PrecipScale}{(optional) named \link{vector} of \link{logical} indicating if the
+mean of the precipitation interpolated on the elevation layers must be
+kept or not, required to create CemaNeige module inputs, default \code{TRUE}
+(the mean of the precipitation is kept to the original value)}
-\item{TempMin}{(optional) \link{matrix} or \link{data.frame} of time series of minimum air temperature [°C], possibly used to create the CemaNeige module inputs}
+\item{TempMean}{(optional) \link{matrix} or \link{data.frame} of time series of mean
+air temperature [°C], required to create the CemaNeige module inputs}
-\item{TempMax}{(optional) \link{matrix} or \link{data.frame} of time series of maximum air temperature [°C], possibly used to create the CemaNeige module inputs}
+\item{TempMin}{(optional) \link{matrix} or \link{data.frame} of time series of minimum
+air temperature [°C], possibly used to create the CemaNeige module inputs}
-\item{ZInputs}{(optional) named \link{vector} of \link{numeric} giving the mean elevation of the Precip and Temp series (before extrapolation) [m], possibly used to create the CemaNeige module input}
+\item{TempMax}{(optional) \link{matrix} or \link{data.frame} of time series of maximum
+air temperature [°C], possibly used to create the CemaNeige module inputs}
-\item{HypsoData}{(optional) \link{matrix} or \link{data.frame} containing 101 \link{numeric} rows: min, q01 to q99 and max of catchment elevation distribution [m], if not defined a single elevation is used for CemaNeige}
+\item{ZInputs}{(optional) named \link{vector} of \link{numeric} giving the mean
+elevation of the Precip and Temp series (before extrapolation) [m],
+possibly used to create the CemaNeige module input}
-\item{NLayers}{(optional) named \link{vector} of \link{numeric} integer giving the number of elevation layers requested \link{-}, required to create CemaNeige module inputs, default=5}
+\item{HypsoData}{(optional) \link{matrix} or \link{data.frame} containing 101 \link{numeric}
+rows: min, q01 to q99 and max of catchment elevation distribution [m],
+if not defined a single elevation is used for CemaNeige}
+
+\item{NLayers}{(optional) named \link{vector} of \link{numeric} integer giving the number
+of elevation layers requested \link{-}, required to create CemaNeige module
+inputs, default=5}
\item{...}{used for compatibility with S3 methods}
}
\value{
-A \emph{GRiwrmInputsModel} object which is a list of \emph{InputsModel} objects created by \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} with one item per modeled sub-catchment.
+A \emph{GRiwrmInputsModel} object which is a list of \emph{InputsModel}
+objects created by \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} with one item per modeled sub-catchment.
}
\description{
Creation of an InputsModel object for a \strong{airGRiwrm} network
}
\details{
-Meteorological data are needed for the nodes of the network that represent a catchment simulated by a rainfall-runoff model. Instead of \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} that has \link{numeric} \link{vector} as time series inputs, this function uses \link{matrix} or \link{data.frame} with the id of the sub-catchment as column names. For single values (\code{ZInputs} or \code{NLayers}), the function requires named \link{vector} with the id of the sub-catchment as name item. If an argument is optional, only the column or the named item has to be provided.
+Meteorological data are needed for the nodes of the network that
+represent a catchment simulated by a rainfall-runoff model. Instead of
+\link[airGR:CreateInputsModel]{airGR::CreateInputsModel} that has \link{numeric} \link{vector} as time series inputs,
+this function uses \link{matrix} or \link{data.frame} with the id of the sub-catchment
+as column names. For single values (\code{ZInputs} or \code{NLayers}), the function
+requires named \link{vector} with the id of the sub-catchment as name item. If an
+argument is optional, only the column or the named item has to be provided.
See \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} documentation for details concerning each input.
+
+Number of rows of \code{Precip}, \code{PotEvap}, \code{Qobs}, \code{Qmin}, \code{TempMean}, \code{TempMin},
+\code{TempMax} must be the same of the length of \code{DatesR} (each row corresponds to
+a time step defined in \code{DatesR}).
+
+For example of use of Direct Injection nodes, see vignettes
+"V03_Open-loop_influenced_flow" and "V04_Closed-loop_regulated_withdrawal".
+
+For example of use of Diversion nodes, see example below and vignette
+"V06_Modelling_regulated_diversion".
}
\examples{
###################################################################
@@ -125,7 +167,7 @@ str(InputsModels)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-# Creation of the GriwmRunOptions object
+# Creation of the GRiwmRunOptions object
RunOptions <- CreateRunOptions(InputsModels,
IndPeriod_Run = Ind_Run)
str(RunOptions)
@@ -140,6 +182,133 @@ OutputsModels <- RunModel(InputsModels,
Param = Param)
str(OutputsModels)
-# Compare Simulation with reservoir and observation of natural flow
-plot(OutputsModels, data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]))
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "\%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
}
diff --git a/man/CreateRunOptions.Rd b/man/CreateRunOptions.Rd
index d244dea0cb88e899186054ea87b6fe41ce6b0aef..51c0478531a6d9864f02d6f654f2b9dbce81ecdd 100644
--- a/man/CreateRunOptions.Rd
+++ b/man/CreateRunOptions.Rd
@@ -110,7 +110,7 @@ str(InputsModels)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-# Creation of the GriwmRunOptions object
+# Creation of the GRiwmRunOptions object
RunOptions <- CreateRunOptions(InputsModels,
IndPeriod_Run = Ind_Run)
str(RunOptions)
@@ -125,6 +125,133 @@ OutputsModels <- RunModel(InputsModels,
Param = Param)
str(OutputsModels)
-# Compare Simulation with reservoir and observation of natural flow
-plot(OutputsModels, data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]))
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "\%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
}
diff --git a/man/CreateSupervisor.Rd b/man/CreateSupervisor.Rd
index e3c95b7cb3b51c78b61d2eadeba5c1f823ca0e3c..053aa41c25ff3b4f114c217b7ff429ce03b292d1 100644
--- a/man/CreateSupervisor.Rd
+++ b/man/CreateSupervisor.Rd
@@ -25,20 +25,128 @@ A \code{Supervisor} object which is an \link{environment} containing all the nec
Creation of a Supervisor for handling regulation in a model
}
\examples{
+###############################################################################
+# An example of reservoir management on an hypothetical dam at station "54095"
+# on the Severn river build to support low-flows at "54057"
+###############################################################################
+# A minimum flow of 20 m3/s is maintained at the dam location and an extra-release
+# is provided when the flow at the downstream station "54057" cross a minimum
+# threshold of 45 m3/s. The dam has a storage capacity of 60 millions m3
+###############################################################################
+library(airGRiwrm)
+
+# Load Severn network information
data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$model <- "RunModel_GR4J"
+
+# Insert a dam downstream the location the gauging station 54095
+# The dam is a direct injection node
+nodes$downstream_id[nodes$gauge_id == "54095"] <- "Dam"
+nodes$distance_downstream[nodes$gauge_id == "54095"] <- 0
+nodes <- rbind(nodes,
+ data.frame(gauge_id = "Dam",
+ downstream_id = "54001",
+ distance_downstream = 42,
+ area = NA,
+ model = NA))
+
griwrm <- CreateGRiwrm(nodes,
list(id = "gauge_id",
down = "downstream_id",
length = "distance_downstream"))
+plot(griwrm)
+
+# Format meteorological inputs for CreateInputs
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
-Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
-InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
-sv <- CreateSupervisor(InputsModel)
+
+# Create a release flow time series for the dam
+# This release will be modified by the Supervisor
+# We initiate it with the natural flow for having a good initialisation of the
+# model at the first time step of the running period
+Qobs <- data.frame(
+ Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
+)
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
+
+# InputsModel object
+IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# Initialisation of the Supervisor
+sv <- CreateSupervisor(IM_severn)
+
+# States of the reservoir are stored in the Supervisor variable
+# The Supervisor variable is an environment which can be available in
+# the controller function for storing and exchange data during the simulation
+sv$Vres <- 0 # Reservoir storage time series
+
+# Dam management is modeled by a controller
+# This controller usually releases Qmin and provides
+# extra release if flow mesured somewhere is below Qthreshold
+# Flow is expressed in m3 / time step
+# Y[1] = runoff flow at gauging station 54095 filling the reservoir
+# Y[2] = flow at gauging station 54057, location of the low-flow objective
+# The returned value is the release calculated at the reservoir
+# We need to enclose the Supervisor variable and other parameters in
+# the environment of the function with a function returning the logic function
+factoryDamLogic <- function(sv, Vmin, Vmax, Qmin, Qthreshold) {
+ function(Y) {
+ # Filling of the reservoir
+ sv$Vres <- c(sv$Vres, tail(sv$Vres, 1) + Y[1])
+ # The release is the max between: overflow, low-flow support and minimum flow
+ U <- U <- max(tail(sv$Vres, 1) - Vmax, Qthreshold - Y[2], Qmin)
+ sv$Vres[length(sv$Vres)] <- tail(sv$Vres, 1) - U
+ if (tail(sv$Vres, 1) < Vmin) {
+ # Reservoir is empty
+ U <- U - (Vmin - tail(sv$Vres, 1))
+ sv$Vres[length(sv$Vres)] <- Vmin
+ }
+ return(U)
+ }
+}
+
+# And define a final function enclosing logic and parameters together
+funDamLogic <- factoryDamLogic(
+ sv = sv, # The Supervisor which store the states of reservoir storage
+ Vmin = 0, # Minimum volume in the reservoir (m3)
+ Vmax = 60 * 1E6, # Maximum volume in the reservoir (m3)
+ Qmin = 20 * 86400, # Min flow to maintain downstream the reservoir (m3/day)
+ Qthreshold = 45 * 86400 # Min flow threshold to support at station 54057 (m3/day)
+)
+
+CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN = funDamLogic)
+
+# GRiwrmRunOptions object simulation of the hydrological year 2002-2003
+IndPeriod_Run <- which(
+ DatesR >= as.POSIXct("2002-10-15", tz = "UTC") &
+ DatesR <= as.POSIXct("2003-10-15", tz = "UTC")
+)
+IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# The Supervisor is used instead of InputsModel for running the model
+OM_dam <- RunModel(sv,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plotting the time series of flows and reservoir storage
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
+ ylim = c(0, 200))
+Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
+ V_reservoir = sv$Vres)
+plot.Qm3s(Vres)
+par(oldpar)
}
diff --git a/man/RunModel.GRiwrmInputsModel.Rd b/man/RunModel.GRiwrmInputsModel.Rd
index 80cc2f5b8cf9e035908fc4c8b483f2cdcca8ccf9..a87734567e37ae701cdfbea76d0696324bba1bf0 100644
--- a/man/RunModel.GRiwrmInputsModel.Rd
+++ b/man/RunModel.GRiwrmInputsModel.Rd
@@ -88,7 +88,7 @@ str(InputsModels)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-# Creation of the GriwmRunOptions object
+# Creation of the GRiwmRunOptions object
RunOptions <- CreateRunOptions(InputsModels,
IndPeriod_Run = Ind_Run)
str(RunOptions)
@@ -103,6 +103,133 @@ OutputsModels <- RunModel(InputsModels,
Param = Param)
str(OutputsModels)
-# Compare Simulation with reservoir and observation of natural flow
-plot(OutputsModels, data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]))
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "\%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
}
diff --git a/man/RunModel.Supervisor.Rd b/man/RunModel.Supervisor.Rd
index 385707baec6f5a7c62c116e80aaac71c9af0d221..540bc339d495580bd6b64944710a3c773954b0ac 100644
--- a/man/RunModel.Supervisor.Rd
+++ b/man/RunModel.Supervisor.Rd
@@ -21,3 +21,129 @@
\description{
RunModel function for a GRiwrmInputsModel object
}
+\examples{
+###############################################################################
+# An example of reservoir management on an hypothetical dam at station "54095"
+# on the Severn river build to support low-flows at "54057"
+###############################################################################
+# A minimum flow of 20 m3/s is maintained at the dam location and an extra-release
+# is provided when the flow at the downstream station "54057" cross a minimum
+# threshold of 45 m3/s. The dam has a storage capacity of 60 millions m3
+###############################################################################
+library(airGRiwrm)
+
+# Load Severn network information
+data(Severn)
+nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$model <- "RunModel_GR4J"
+
+# Insert a dam downstream the location the gauging station 54095
+# The dam is a direct injection node
+nodes$downstream_id[nodes$gauge_id == "54095"] <- "Dam"
+nodes$distance_downstream[nodes$gauge_id == "54095"] <- 0
+nodes <- rbind(nodes,
+ data.frame(gauge_id = "Dam",
+ downstream_id = "54001",
+ distance_downstream = 42,
+ area = NA,
+ model = NA))
+
+griwrm <- CreateGRiwrm(nodes,
+ list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream"))
+plot(griwrm)
+
+# Format meteorological inputs for CreateInputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+
+# Create a release flow time series for the dam
+# This release will be modified by the Supervisor
+# We initiate it with the natural flow for having a good initialisation of the
+# model at the first time step of the running period
+Qobs <- data.frame(
+ Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
+)
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
+
+# InputsModel object
+IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# Initialisation of the Supervisor
+sv <- CreateSupervisor(IM_severn)
+
+# States of the reservoir are stored in the Supervisor variable
+# The Supervisor variable is an environment which can be available in
+# the controller function for storing and exchange data during the simulation
+sv$Vres <- 0 # Reservoir storage time series
+
+# Dam management is modeled by a controller
+# This controller usually releases Qmin and provides
+# extra release if flow mesured somewhere is below Qthreshold
+# Flow is expressed in m3 / time step
+# Y[1] = runoff flow at gauging station 54095 filling the reservoir
+# Y[2] = flow at gauging station 54057, location of the low-flow objective
+# The returned value is the release calculated at the reservoir
+# We need to enclose the Supervisor variable and other parameters in
+# the environment of the function with a function returning the logic function
+factoryDamLogic <- function(sv, Vmin, Vmax, Qmin, Qthreshold) {
+ function(Y) {
+ # Filling of the reservoir
+ sv$Vres <- c(sv$Vres, tail(sv$Vres, 1) + Y[1])
+ # The release is the max between: overflow, low-flow support and minimum flow
+ U <- U <- max(tail(sv$Vres, 1) - Vmax, Qthreshold - Y[2], Qmin)
+ sv$Vres[length(sv$Vres)] <- tail(sv$Vres, 1) - U
+ if (tail(sv$Vres, 1) < Vmin) {
+ # Reservoir is empty
+ U <- U - (Vmin - tail(sv$Vres, 1))
+ sv$Vres[length(sv$Vres)] <- Vmin
+ }
+ return(U)
+ }
+}
+
+# And define a final function enclosing logic and parameters together
+funDamLogic <- factoryDamLogic(
+ sv = sv, # The Supervisor which store the states of reservoir storage
+ Vmin = 0, # Minimum volume in the reservoir (m3)
+ Vmax = 60 * 1E6, # Maximum volume in the reservoir (m3)
+ Qmin = 20 * 86400, # Min flow to maintain downstream the reservoir (m3/day)
+ Qthreshold = 45 * 86400 # Min flow threshold to support at station 54057 (m3/day)
+)
+
+CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN = funDamLogic)
+
+# GRiwrmRunOptions object simulation of the hydrological year 2002-2003
+IndPeriod_Run <- which(
+ DatesR >= as.POSIXct("2002-10-15", tz = "UTC") &
+ DatesR <= as.POSIXct("2003-10-15", tz = "UTC")
+)
+IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# The Supervisor is used instead of InputsModel for running the model
+OM_dam <- RunModel(sv,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plotting the time series of flows and reservoir storage
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
+ ylim = c(0, 200))
+Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
+ V_reservoir = sv$Vres)
+plot.Qm3s(Vres)
+par(oldpar)
+}
diff --git a/man/getNoSD_Ids.Rd b/man/getNoSD_Ids.Rd
index 50750697f100619ec9e4813a74cc898fa21400aa..644d49eaf4ad79b18c723423d51f58d8ecf10e4d 100644
--- a/man/getNoSD_Ids.Rd
+++ b/man/getNoSD_Ids.Rd
@@ -4,10 +4,12 @@
\alias{getNoSD_Ids}
\title{Function to obtain the ID of sub-basins not using SD model}
\usage{
-getNoSD_Ids(InputsModel)
+getNoSD_Ids(InputsModel, include_diversion = TRUE)
}
\arguments{
\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+
+\item{include_diversion}{\link{logical} for including diversion nodes}
}
\value{
\link{character} IDs of the sub-basins not using the SD model
diff --git a/man/getSD_Ids.Rd b/man/getSD_Ids.Rd
index 1a863101f37fddc766c1d33d37152fd54350cf5a..a1987eca097fe7e47c1da1a9408e39dadae1ede8 100644
--- a/man/getSD_Ids.Rd
+++ b/man/getSD_Ids.Rd
@@ -4,10 +4,12 @@
\alias{getSD_Ids}
\title{Function to obtain the ID of sub-basins using SD model}
\usage{
-getSD_Ids(InputsModel)
+getSD_Ids(InputsModel, add_diversion = FALSE)
}
\arguments{
\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+
+\item{add_diversion}{\link{logical} for adding upstream nodes with diversion}
}
\value{
\link{character} IDs of the sub-basins using SD model
diff --git a/man/plot.GRiwrm.Rd b/man/plot.GRiwrm.Rd
index 280c22d5708d8c29b6db9a4ee4244f40b9d8e16b..ff59579bf370a7530878c2b8120e07f360a2a7c6 100644
--- a/man/plot.GRiwrm.Rd
+++ b/man/plot.GRiwrm.Rd
@@ -10,8 +10,8 @@
orientation = "LR",
width = "100\%",
height = "100\%",
- box_colors = c(UpstreamUngauged = "#eef", UpstreamGauged = "#aaf", IntermUngauged =
- "#efe", IntermGauged = "#afa", DirectInjection = "#faa"),
+ box_colors = c(UpstreamUngauged = "#eef", UpstreamGauged = "#aaf", IntermediateUngauged
+ = "#efe", IntermediateGauged = "#afa", DirectInjection = "#faa"),
...
)
}
@@ -41,7 +41,9 @@ This function only works inside RStudio because the HTMLwidget produced by Diagr
is not handled on some platforms
}
\examples{
-# Network of 2 nodes distant of 150 km:
+#########################################
+# Network of 2 nodes distant of 150 km: #
+#########################################
# - an upstream reservoir modelled as a direct flow injection (no model)
# - a gauging station downstream a catchment of 360 km² modelled with GR4J
db <- data.frame(id = c("Reservoir", "GaugingDown"),
@@ -55,7 +57,9 @@ griwrm_basic
# Network diagram with direct flow node in red, intermediate sub-basin in green
plot(griwrm_basic)
-# GR4J semi-distributed model of the Severn River
+###################################################
+# GR4J semi-distributed model of the Severn River #
+###################################################
data(Severn)
nodes <- Severn$BasinsInfo
nodes$model <- "RunModel_GR4J"
@@ -69,7 +73,9 @@ griwrm_severn
# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
plot(griwrm_severn)
-# Same model with an ungauged station at nodes 54029 and 54001
+####################################################################
+# Severn network with an ungauged station at nodes 54029 and 54001 #
+####################################################################
# By default the first gauged node at downstream is used for parameter calibration (54032)
nodes_ungauged <- nodes
nodes_ungauged$model[nodes_ungauged$gauge_id \%in\% c("54029", "54001")] <- "Ungauged"
@@ -78,4 +84,18 @@ griwrm_ungauged <- CreateGRiwrm(nodes_ungauged, rename_columns)
griwrm_ungauged
# Network diagram with gauged nodes of vivid color, and ungauged nodes of dull color
plot(griwrm_ungauged)
+
+#######################################################
+# Severn network with a Diversion on the node "54029" #
+# which transfer flows to the node "54001" #
+#######################################################
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "model", "area")]
+nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54029",
+ downstream_id = "54001",
+ distance_downstream = 20,
+ model = "Diversion",
+ area = NA))
+griwrm_div <- CreateGRiwrm(nodes_div, rename_columns)
+# Network diagram figures Diversion node by a red frame and a red arrow
+plot(griwrm_div)
}
diff --git a/man/plot.GRiwrmOutputsModel.Rd b/man/plot.GRiwrmOutputsModel.Rd
index f25a416c29c271b075685bf13ed5fd8e7a6b024e..a2f12a09d49e28228f10befe657c404c61ca4fa0 100644
--- a/man/plot.GRiwrmOutputsModel.Rd
+++ b/man/plot.GRiwrmOutputsModel.Rd
@@ -88,7 +88,7 @@ str(InputsModels)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
-# Creation of the GriwmRunOptions object
+# Creation of the GRiwmRunOptions object
RunOptions <- CreateRunOptions(InputsModels,
IndPeriod_Run = Ind_Run)
str(RunOptions)
@@ -103,6 +103,133 @@ OutputsModels <- RunModel(InputsModels,
Param = Param)
str(OutputsModels)
-# Compare Simulation with reservoir and observation of natural flow
-plot(OutputsModels, data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]))
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "\%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
}
diff --git a/man/plot.Qm3s.Rd b/man/plot.Qm3s.Rd
index 7c2ddcf3df9604fcd5ff5ed9494ab104719ba9fc..2bb899b633fa59113db059f8936f716c3e485eb2 100644
--- a/man/plot.Qm3s.Rd
+++ b/man/plot.Qm3s.Rd
@@ -8,7 +8,7 @@
x,
type = "l",
xlab = "Date",
- ylab = expression("Flow (m"^"3" * "/s)"),
+ ylab = expression("Flow rate (m"^"3" * "/s)"),
main = "Simulated flows",
col = rainbow(ncol(x) - 1),
legend = colnames(x)[-1],
@@ -20,7 +20,8 @@
)
}
\arguments{
-\item{x}{\link{data.frame} with a first column with \link{POSIXt} dates and followings columns with flows at each node of the network}
+\item{x}{\link{data.frame} with a first column with \link{POSIXt} dates and followings
+columns with flows at each node of the network}
\item{type}{\link{character} plot type (See \link{plot.default}), default "l"}
@@ -32,7 +33,8 @@
\item{col}{\link{character} plotting color (See \link{par}), default to rainbow colors}
-\item{legend}{\link{character} see parameter \code{legend} of \link{legend}. Set to \link{NULL} if display of the legend is not wanted}
+\item{legend}{\link{character} see parameter \code{legend} of \link{legend}. Set it to \link{NULL}
+to hide the legend}
\item{legend.cex}{\link{character} \code{cex} parameter for the text of the legend (See \link{par})}
@@ -46,5 +48,219 @@
Screen plot window.
}
\description{
-Plot of a \code{Qm3s} object (time series of simulated flows)
+This function plot time series of flow rate in m3/s. It's a method for object
+of class "Qm3s" which can be directly called by \code{plot}. It can also be called
+as a function \code{plot.Qm3s} if the first parameter has the good format.
+}
+\examples{
+###################################################################
+# Run the `airGR::RunModel_Lag` example in the GRiwrm fashion way #
+# Simulation of a reservoir with a purpose of low-flow mitigation #
+###################################################################
+
+## ---- preparation of the InputsModel object
+
+## loading package and catchment data
+library(airGRiwrm)
+data(L0123001)
+
+## ---- specifications of the reservoir
+
+## the reservoir withdraws 1 m3/s when it's possible considering the flow observed in the basin
+Qupstream <- matrix(-sapply(BasinObs$Qls / 1000 - 1, function(x) {
+ min(1, max(0, x, na.rm = TRUE))
+}), ncol = 1)
+
+## except between July and September when the reservoir releases 3 m3/s for low-flow mitigation
+month <- as.numeric(format(BasinObs$DatesR, "\%m"))
+Qupstream[month >= 7 & month <= 9] <- 3
+Qupstream <- Qupstream * 86400 ## Conversion in m3/day
+
+## the reservoir is not an upstream subcachment: its areas is NA
+BasinAreas <- c(NA, BasinInfo$BasinArea)
+
+## delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
+LengthHydro <- 150
+## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
+Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
+
+# This example is a network of 2 nodes which can be describe like this:
+db <- data.frame(id = c("Reservoir", "GaugingDown"),
+ length = c(LengthHydro, NA),
+ down = c("GaugingDown", NA),
+ area = c(NA, BasinInfo$BasinArea),
+ model = c(NA, "RunModel_GR4J"),
+ stringsAsFactors = FALSE)
+
+# Create GRiwrm object from the data.frame
+griwrm <- CreateGRiwrm(db)
+str(griwrm)
+
+# Formatting observations for the hydrological models
+# Each input data should be a matrix or a data.frame with the good id in the name of the column
+Precip <- matrix(BasinObs$P, ncol = 1)
+colnames(Precip) <- "GaugingDown"
+PotEvap <- matrix(BasinObs$E, ncol = 1)
+colnames(PotEvap) <- "GaugingDown"
+
+# Observed flows contain flows that are directly injected in the model
+Qobs = matrix(Qupstream, ncol = 1)
+colnames(Qobs) <- "Reservoir"
+
+# Creation of the GRiwrmInputsModel object (= a named list of InputsModel objects)
+InputsModels <- CreateInputsModel(griwrm,
+ DatesR = BasinObs$DatesR,
+ Precip = Precip,
+ PotEvap = PotEvap,
+ Qobs = Qobs)
+str(InputsModels)
+
+## run period selection
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1990-01-01"),
+ which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d")=="1999-12-31"))
+
+# Creation of the GRiwmRunOptions object
+RunOptions <- CreateRunOptions(InputsModels,
+ IndPeriod_Run = Ind_Run)
+str(RunOptions)
+
+# Parameters of the SD models should be encapsulated in a named list
+ParamGR4J <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
+Param <- list(`GaugingDown` = c(Velocity, ParamGR4J))
+
+# RunModel for the whole network
+OutputsModels <- RunModel(InputsModels,
+ RunOptions = RunOptions,
+ Param = Param)
+str(OutputsModels)
+
+# Compare regimes of the simulation with reservoir and observation of natural flow
+plot(OutputsModels,
+ data.frame(GaugingDown = BasinObs$Qmm[Ind_Run]),
+ which = "Regime")
+
+# Plot together simulated flows (m3/s) of the reservoir and the gauging station
+plot(attr(OutputsModels, "Qm3s"))
+
+
+########################################################
+# Run the Severn example provided with this package #
+# A natural catchment composed with 6 gauging stations #
+########################################################
+
+data(Severn)
+nodes <- Severn$BasinsInfo
+nodes$model <- "RunModel_GR4J"
+# Mismatch column names are renamed to stick with GRiwrm requirements
+rename_columns <- list(id = "gauge_id",
+ down = "downstream_id",
+ length = "distance_downstream")
+g_severn <- CreateGRiwrm(nodes, rename_columns)
+
+# Network diagram with upstream basin nodes in blue, intermediate sub-basin in green
+plot(g_severn)
+
+# Format CAMEL-GB meteorological dataset for airGRiwrm inputs
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+
+# Precipitation and Potential Evaporation are related to the whole catchment
+# at each gauging station. We need to compute them for intermediate catchments
+# for use in a semi-distributed model
+Precip <- ConvertMeteoSD(g_severn, PrecipTot)
+PotEvap <- ConvertMeteoSD(g_severn, PotEvapTot)
+
+# CreateInputsModel object
+IM_severn <- CreateInputsModel(g_severn, DatesR, Precip, PotEvap)
+
+# GRiwrmRunOptions object
+# Run period is set aside the one-year warm-up period
+IndPeriod_Run <- seq(
+ which(IM_severn[[1]]$DatesR == (IM_severn[[1]]$DatesR[1] + 365*24*60*60)),
+ length(IM_severn[[1]]$DatesR) # Until the end of the time series
+)
+IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
+
+RO_severn <- CreateRunOptions(
+ IM_severn,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+)
+
+# Load parameters of the model from Calibration in vignette V02
+P_severn <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+
+# Run the simulation
+OM_severn <- RunModel(IM_severn,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Plot results of simulated flows in m3/s
+Qm3s <- attr(OM_severn, "Qm3s")
+plot(Qm3s[1:150, ])
+
+
+##################################################################
+# An example of water withdrawal for irrigation with restriction #
+# modeled with a Diversion node on the Severn river #
+##################################################################
+
+# A diversion is added at gauging station "54001"
+nodes_div <- nodes[, c("gauge_id", "downstream_id", "distance_downstream", "area", "model")]
+names(nodes_div) <- c("id", "down", "length", "area", "model")
+nodes_div <- rbind(nodes_div,
+ data.frame(id = "54001", # location of the diversion
+ down = NA, # the abstracted flow goes outside
+ length = NA, # down=NA, so length=NA
+ area = NA, # no area, diverted flow is in m3/day
+ model = "Diversion"))
+
+g_div <- CreateGRiwrm(nodes_div)
+# The node "54001" is surrounded in red to show the diverted node
+plot(g_div)
+
+# Computation of the irrigation withdraw objective
+irrigMonthlyPlanning <- c(0.0, 0.0, 1.2, 2.4, 3.2, 3.6, 3.6, 2.8, 1.8, 0.0, 0.0, 0.0)
+names(irrigMonthlyPlanning) <- month.abb
+irrigMonthlyPlanning
+DatesR_month <- as.numeric(format(DatesR, "\%m"))
+# Withdrawn flow calculated for each day is negative
+Qirrig <- matrix(-irrigMonthlyPlanning[DatesR_month] * 86400, ncol = 1)
+colnames(Qirrig) <- "54001"
+
+# Minimum flow to remain downstream the diversion is 12 m3/s
+Qmin <- matrix(12 * 86400, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54001"
+
+# Creation of GRimwrInputsModel object
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qirrig, Qmin = Qmin)
+
+# RunOptions and parameters are unchanged, we can directly run the simulation
+OM_div <- RunModel(IM_div,
+ RunOptions = RO_severn,
+ Param = P_severn)
+
+# Retrieve diverted flow at "54001" and convert it from m3/day to m3/s
+Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
+
+# Plot the diverted flow for the year 2003
+Ind_Plot <- which(
+ OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+ OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
+)
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
+ Diverted_flow = Qdiv_m3s[Ind_Plot])
+
+oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001"
+df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
+ attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
+names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
+plot.Qm3s(df54001, ylim = c(0,70))
+abline(h = 12, col = "green")
+par(oldpar)
}
diff --git a/tests/testthat/helper_RunModel.R b/tests/testthat/helper_RunModel.R
index 2b5b500e8743b67327376ce9b2db44c80d9e6679..20ce78b829f5b9711f046abde6dbbe96bf4d365a 100644
--- a/tests/testthat/helper_RunModel.R
+++ b/tests/testthat/helper_RunModel.R
@@ -34,16 +34,9 @@ setupRunModel <-
# Set network
if(is.null(griwrm)) {
- nodes <-
- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
- nodes$model <- "RunModel_GR4J"
+ nodes <- loadSevernNodes()
griwrm <-
- CreateGRiwrm(nodes,
- list(
- id = "gauge_id",
- down = "downstream_id",
- length = "distance_downstream"
- ))
+ CreateGRiwrm(nodes)
}
# Convert meteo data to SD (remove upstream areas)
@@ -54,7 +47,7 @@ setupRunModel <-
if (!runInputsModel)
return(environment())
InputsModel <-
- suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs))
+ suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap))
# RunOptions
if (!runRunOptions)
@@ -126,3 +119,19 @@ setupRunOptions <- function(InputsModel) {
return(environment())
}
+
+#' Load Severn example and format BasinsInfo for use in CreateGRiwrm
+#'
+#' @return [data.frame] with columns "id", "down", "length", "area"
+#' @noRd
+#'
+#' @examples
+#' nodes <- loadNodes
+loadSevernNodes <- function() {
+ data(Severn)
+ nodes <-
+ Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+ names(nodes) <- c("id", "down", "length", "area")
+ nodes$model <- "RunModel_GR4J"
+ return(nodes)
+}
diff --git a/tests/testthat/test-Calibration.R b/tests/testthat/test-Calibration.R
index 8a27cb827a4fe5c197b24f71add625ed0b08f1f0..8ff1ccc36d38d22aa65b353dcdb3a0da720e2948 100644
--- a/tests/testthat/test-Calibration.R
+++ b/tests/testthat/test-Calibration.R
@@ -98,3 +98,29 @@ test_that("Calibration with regularization is OK", {
)
})
})
+
+test_that("Calibration with Diversion works", {
+ n_div <- rbind(nodes,
+ data.frame(id = "54029", down = "54002", length = 50, area = NA, model = "Diversion"))
+ g_div <- CreateGRiwrm(n_div)
+ Qmin = matrix(1E5, nrow = length(DatesR), ncol = 1)
+ colnames(Qmin) = "54029"
+ Qdiv <- -Qmin
+ IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qdiv, Qmin = Qmin)
+ RO_div <- setupRunOptions(IM_div)$RunOptions
+ P_div <- ParamMichel
+ P_div$`54002` <- c(1, ParamMichel$`54002`)
+ IC_div <- CreateInputsCrit(
+ InputsModel = IM_div,
+ RunOptions = RO_div,
+ Obs = Qobs[IndPeriod_Run,],
+ )
+ CO_div <- CreateCalibOptions(IM_div)
+ OC <- Calibration(
+ InputsModel = IM_div,
+ RunOptions = RO_div,
+ InputsCrit = IC_div,
+ CalibOptions = CO_div
+ )
+ expect_length(OC$`54002`$ParamFinalR, 5)
+})
diff --git a/tests/testthat/test-CreateInputsCrit.R b/tests/testthat/test-CreateInputsCrit.R
index 3d02d28056fedb515635e283026b6655c0f7b45b..3b419ae5b8401affc4142881864ab97e9fee5c81 100644
--- a/tests/testthat/test-CreateInputsCrit.R
+++ b/tests/testthat/test-CreateInputsCrit.R
@@ -122,13 +122,9 @@ test_that("Lavenne criterion: current node and a priori node must use the same m
})
test_that("Ungauged node as Apriori node should throw an error", {
- nodes$model[nodes$gauge_id == "54001"] <- "Ungauged"
- griwrm <- CreateGRiwrm(
- nodes,
- list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
- )
- InputsModel <-
- suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs))
+ nodes$model[nodes$id == "54001"] <- "Ungauged"
+ griwrm <- CreateGRiwrm(nodes)
+ InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap)
expect_error(
CreateInputsCrit(InputsModel = InputsModel,
RunOptions = RunOptions,
diff --git a/tests/testthat/test-CreateInputsModel.R b/tests/testthat/test-CreateInputsModel.R
index 4622038f1b169fe6e7352984f92559ebc5f0af0b..d383a5c856d22697c8e70c332a2e6a8a6be078c7 100644
--- a/tests/testthat/test-CreateInputsModel.R
+++ b/tests/testthat/test-CreateInputsModel.R
@@ -127,7 +127,7 @@ test_that("throws error when missing CemaNeige data", {
regexp = "'TempMean' is missing")
})
-test_that("throws error when missing Qobs on node not related to an hydrological model", {
+test_that("throws error when missing Qobs on node Direct Injection node", {
l$griwrm$model[1] <- NA
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
@@ -159,14 +159,16 @@ test_that("must works with node not related to an hydrological model", {
expect_equal(colnames(IM[[2]]$Qupstream), c("Up1", "Up2"))
})
-test_that("negative observed flow on catchment should throw error", {
- l$Qobs[100, 1] <- -99
+test_that("Qobs on hydrological nodes should throw en error", {
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
- Qobs = l$Qobs),
- regexp = "Negative flow found")
+ Qobs = l$Qobs,
+ TempMean = l$TempMean,
+ ZInputs = l$ZInputs,
+ HypsoData = l$HypsoData),
+ regexp = "columns in 'Qobs' don't match with")
l$griwrm$model[1] <- NA
expect_s3_class(suppressWarnings(
CreateInputsModel(
@@ -174,7 +176,7 @@ test_that("negative observed flow on catchment should throw error", {
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
- Qobs = l$Qobs,
+ Qobs = l$Qobs[,1, drop = F],
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData
@@ -191,13 +193,51 @@ for(x in ls(e)) assign(x, get(x, e))
test_that("Ungauged node should inherits its FUN_MOD from the downstream gauged node", {
- nodes$model[nodes$gauge_id == "54032"] <- "Ungauged"
- griwrmV05 <- CreateGRiwrm(
- nodes,
- list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
- )
+ nodes$model[nodes$id == "54032"] <- "Ungauged"
+ griwrmV05 <- CreateGRiwrm(nodes)
IM <- suppressWarnings(
- CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap, Qobs)
+ CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap)
)
expect_equal(IM[["54032"]]$FUN_MOD, "RunModel_GR4J")
})
+
+test_that("Network with Diversion works", {
+ n_div <- rbind(nodes, data.frame(id = "54029",
+ down = "54002",
+ length = 20,
+ model = "Diversion",
+ area = NA))
+ g <- CreateGRiwrm(n_div)
+ Qobs = matrix(-1, nrow = length(DatesR), ncol = 1)
+ colnames(Qobs) = "54029"
+ IM <- suppressWarnings(
+ CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs)
+ )
+ expect_equal(IM[["54032"]]$UpstreamNodes, c("54001", "54029"))
+ expect_equal(IM[["54032"]]$UpstreamVarQ , c("Qsim_m3", "Qsim_m3"))
+ expect_equal(IM[["54002"]]$UpstreamNodes, "54029")
+ expect_equal(IM[["54002"]]$UpstreamIsModeled , TRUE)
+ expect_equal(IM[["54002"]]$UpstreamVarQ , "Qdiv_m3")
+ expect_equivalent(IM$`54029`$Qmin, matrix(0, nrow = length(DatesR), ncol = 1))
+})
+
+test_that("Diversion node: checks about 'Qmin'", {
+ n_div <- rbind(nodes,
+ data.frame(id = "54029", down = "54002", length = 50, area = NA, model = "Diversion"))
+ g <- CreateGRiwrm(n_div)
+ Qobs = matrix(-1, nrow = length(DatesR), ncol = 1)
+ colnames(Qobs) = "54029"
+ expect_warning(CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs = Qobs),
+ regexp = "Zero values")
+ Qmin <- -Qobs
+ IM <- CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = Qmin)
+ expect_equivalent(IM$`54029`$Qmin, Qmin)
+ QminNA <- Qmin
+ QminNA[1] <- NA
+ expect_error(CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = QminNA),
+ regexp = "NA")
+ QminBadCol <- Qmin
+ colnames(QminBadCol) = "54002"
+ expect_error(CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = QminBadCol),
+ regexp = "columns that does not match with IDs of Diversion nodes")
+})
diff --git a/tests/testthat/test-CreateSupervisor.R b/tests/testthat/test-CreateSupervisor.R
new file mode 100644
index 0000000000000000000000000000000000000000..28bd79d2c75452ca22b218e40761dbc74e06eeca
--- /dev/null
+++ b/tests/testthat/test-CreateSupervisor.R
@@ -0,0 +1,45 @@
+# data set up
+e <- setupRunModel(runInputsModel = FALSE)
+# variables are copied from environment 'e' to the current environment
+# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
+for(x in ls(e)) assign(x, get(x, e))
+
+# Add 2 nodes to the network
+nodes2 <- rbind(nodes,
+ data.frame(
+ id = c("R1", "R2"),
+ down = "54057",
+ length = 100,
+ area = NA,
+ model = NA
+ ))
+griwrm2 <- CreateGRiwrm(nodes2)
+
+# Add Qobs for the 2 new nodes and create InputsModel
+Qobs <- matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2)
+colnames(Qobs) <- c("R1", "R2")
+InputsModel <-
+ CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs)
+sv <- CreateSupervisor(InputsModel)
+
+test_that("Checks in CreateSupervisor", {
+ expect_error(CreateSupervisor(InputsModel, 3.14), regexp = "integer")
+ expect_error(CreateSupervisor(InputsModel, 0L), regexp = "positive")
+ expect_s3_class(sv, "Supervisor")
+})
+
+test_that("Checks in CreateController", {
+ FUN <- function(Y) return(0)
+ expect_error(CreateController(sv, 0, "54057", "R1", FUN),
+ regexp = "character")
+ expect_error(CreateController(sv, c("toto1", "toto2"), "54057", "R1", FUN),
+ regexp = "length")
+ expect_error(CreateController(sv, "toto", "fake", "R1", FUN),
+ regexp = "GRiwrm")
+ expect_error(CreateController(sv, "toto", "54057", "54057", FUN),
+ regexp = "DirectInjection")
+ CreateController(sv, "toto", "54057", "R1", FUN)
+ expect_s3_class(sv$controllers, "Controllers")
+ expect_s3_class(sv$controllers$toto, "Controller")
+})
+
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
index 4b48b0e0abed15d40a68d00e277ae66fba6a8380..27868ca43cd226af01dd04738ea3acb5b5afc8bb 100644
--- a/tests/testthat/test-RunModel.R
+++ b/tests/testthat/test-RunModel.R
@@ -4,8 +4,6 @@ e <- setupRunModel()
# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
for(x in ls(e)) assign(x, get(x, e))
-context("RunModel.GRiwrmInputsModel")
-
test_that("RunModel.GRiwrmInputsModel should return same result with separated warm-up", {
RO_WarmUp <- CreateRunOptions(
InputsModel,
@@ -34,8 +32,6 @@ test_that("RunModel.GRiwrmInputsModel should return same result with separated w
})
})
-context("RunModel.Supervisor")
-
test_that("RunModel.Supervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
sv <- CreateSupervisor(InputsModel)
OM_Supervisor <- RunModel(
@@ -48,61 +44,6 @@ test_that("RunModel.Supervisor with no regulation should returns same results as
})
})
-# Add 2 nodes to the network
-nodes2 <- rbind(nodes,
- data.frame(
- gauge_id = c("R1", "R2"),
- downstream_id = "54057",
- distance_downstream = 100,
- area = NA,
- model = NA
- ))
-griwrm2 <- CreateGRiwrm(nodes2,
- list(id = "gauge_id",
- down = "downstream_id",
- length = "distance_downstream"))
-
-# Add Qobs for the 2 new nodes and create InputsModel
-Qobs2 <- cbind(Qobs, matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2))
-colnames(Qobs2) <- c(colnames(Qobs2)[1:6], "R1", "R2")
-InputsModel <- suppressWarnings(
- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs2)
-)
-
-test_that("RunModel.Supervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
- # Create Supervisor
- sv <- CreateSupervisor(InputsModel)
- # Function to withdraw half of the measured flow
- fWithdrawal <- function(y) { -y/2 }
- # Function to release half of the the measured flow
- fRelease <- function(y) { y/2 }
- # Controller that withdraw half of the flow measured at node "54002" at location "R1"
- CreateController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
- # Controller that release half of the flow measured at node "54002" at location "R2"
- CreateController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
-
- OM_Supervisor <- RunModel(
- sv,
- RunOptions = RunOptions,
- Param = ParamMichel
- )
- expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
-})
-
-test_that("RunModel.Supervisor with multi time steps controller, two regulations in 1 centralised controller that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
- sv <- CreateSupervisor(InputsModel, TimeStep = 10L)
- fEverything <- function(y) {
- matrix(c(y[,1]/2, -y[,1]/2), ncol = 2)
- }
- CreateController(sv, "Everything", Y = c("54002", "54032"), U = c("R1", "R2"), FUN = fEverything)
- OM_Supervisor <- RunModel(
- sv,
- RunOptions = RunOptions,
- Param = ParamMichel
- )
- expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
-})
-
test_that("RunModel.GRiwrmInputsModel handles CemaNeige", {
l <- setUpCemaNeigeData()
l$griwrm[l$griwrm$id == "Down", "model"] <- "RunModel_GR4J"
@@ -118,8 +59,7 @@ test_that("RunModel.GRiwrmInputsModel handles CemaNeige", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
- HypsoData = l$HypsoData,
- Qobs = l$Qobs
+ HypsoData = l$HypsoData
)
)
## run period selection
@@ -145,25 +85,45 @@ test_that("RunModel.GRiwrmInputsModel handles CemaNeige", {
expect_equal(Qm3s[,4], rowSums(Qm3s[,2:3]))
})
-test_that("RunModel.Supervisor with NA values in Qupstream", {
- # Create Supervisor
- InputsModel$`54057`$Qupstream[, c("R1", "R2")] <- NA
- sv <- CreateSupervisor(InputsModel)
- # Function to withdraw half of the measured flow
- fWithdrawal <- function(y) { -y/2 }
- # Function to release half of the the measured flow
- fRelease <- function(y) { y/2 }
- # Controller that withdraw half of the flow measured at node "54002" at location "R1"
- CreateController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
- # Controller that release half of the flow measured at node "54002" at location "R2"
- CreateController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
+n_div <- rbind(nodes,
+ data.frame(id = "54029", down = "54002", length = 50, area = NA, model = "Diversion"))
+g_div <- CreateGRiwrm(n_div)
+Qmin = matrix(1E5, nrow = length(DatesR), ncol = 1)
+colnames(Qmin) = "54029"
+Qobs <- -Qmin
+IM_div <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = Qmin)
+RO_div <- setupRunOptions(IM_div)$RunOptions
+P_div <- ParamMichel
+P_div$`54002` <- c(1, ParamMichel$`54002`)
- OM_Supervisor <- RunModel(
- sv,
- RunOptions = RunOptions,
- Param = ParamMichel
- )
- expect_equal(OM_Supervisor[["54057"]]$Qsim[1:3], rep(as.double(NA),3))
- expect_equal(OM_Supervisor[["54057"]]$Qsim[4:length(IndPeriod_Run)],
- OM_GriwrmInputs[["54057"]]$Qsim[4:length(IndPeriod_Run)])
+test_that("RunModel_Diversion with zero diversion equals no diversion", {
+ Qobs[, ] <- 0
+ IM <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = Qmin)
+ OM <- RunModel(IM, RunOptions = RO_div, Param = P_div)
+ expect_s3_class(OM, "GRiwrmOutputsModel")
+ lapply(names(OM), function(id) {
+ expect_equal(OM[[!!id]]$Qsim, OM_GriwrmInputs[[!!id]]$Qsim)
+ expect_equal(OM[[!!id]]$Qsim_m3, OM_GriwrmInputs[[!!id]]$Qsim_m3)
+ expect_equal(OM[[!!id]]$RunOptions$WarmUpQsim, OM_GriwrmInputs[[!!id]]$RunOptions$WarmUpQsim)
+ expect_equal(OM[[!!id]]$RunOptions$WarmUpQsim_m3, OM_GriwrmInputs[[!!id]]$RunOptions$WarmUpQsim_m3)
+ })
+})
+
+test_that("Huge diversion would result in Qsim_m3 == Qmin", {
+ Qobs[, ] <- -1E12
+ IM <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = Qmin)
+ OM <- RunModel(IM, RunOptions = RO_div, Param = P_div)
+ expect_equal(OM[["54029"]]$Qsim_m3, Qmin[RunOptions[[1]]$IndPeriod_Run])
+ expect_equal(OM[["54029"]]$Qsim,
+ Qmin[RunOptions[[1]]$IndPeriod_Run] /
+ g_div$area[g_div$id == "54029" & g_div$model != "Diversion"] / 1E3)
+})
+
+test_that("Huge minimum remaining flow results in Qdiv = 0", {
+ Qobs[, ] <- -1000
+ Qmin[, ] <- 1E12
+ IM <- CreateInputsModel(g_div, DatesR, Precip, PotEvap, Qobs = Qobs, Qmin = Qmin)
+ OM <- RunModel(IM, RunOptions = RO_div, Param = P_div)
+ expect_equal(OM[["54029"]]$Qsim, OM[["54029"]]$Qnat)
+ expect_equal(OM[["54029"]]$Qdiv_m3, rep(0, length(IndPeriod_Run)))
})
diff --git a/tests/testthat/test-RunModel.Supervisor.R b/tests/testthat/test-RunModel.Supervisor.R
new file mode 100644
index 0000000000000000000000000000000000000000..d1f1636ab00cfa8f87b802ab843eeb9388e0b1d1
--- /dev/null
+++ b/tests/testthat/test-RunModel.Supervisor.R
@@ -0,0 +1,82 @@
+# data set up
+e <- setupRunModel()
+# variables are copied from environment 'e' to the current environment
+# https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
+for(x in ls(e)) assign(x, get(x, e))
+
+# Add 2 nodes to the network
+nodes2 <- rbind(nodes,
+ data.frame(
+ id = c("R1", "R2"),
+ down = "54057",
+ length = 100,
+ area = NA,
+ model = NA
+ ))
+griwrm2 <- CreateGRiwrm(nodes2)
+
+# Add Qobs for the 2 new nodes and create InputsModel
+Qobs <- matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2)
+colnames(Qobs) <- c("R1", "R2")
+InputsModel <-
+ CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs)
+
+test_that("RunModel.Supervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+ # Create Supervisor
+ sv <- CreateSupervisor(InputsModel)
+ # Function to withdraw half of the measured flow
+ fWithdrawal <- function(y) { -y/2 }
+ # Function to release half of the the measured flow
+ fRelease <- function(y) { y/2 }
+ # Controller that withdraw half of the flow measured at node "54002" at location "R1"
+ CreateController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
+ # Controller that release half of the flow measured at node "54002" at location "R2"
+ CreateController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
+
+ OM_Supervisor <- RunModel(
+ sv,
+ RunOptions = RunOptions,
+ Param = ParamMichel
+ )
+ expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+test_that("RunModel.Supervisor with multi time steps controller, two regulations
+ in 1 centralised controller that cancel each other out should returns
+ same results as RunModel.GRiwrmInputsModel", {
+ sv <- CreateSupervisor(InputsModel, TimeStep = 10L)
+ fEverything <- function(y) {
+ m <- matrix(c(y[,1]/2, -y[,1]/2), ncol = 2)
+ }
+ CreateController(sv, "Everything", Y = c("54002", "54032"), U = c("R1", "R2"), FUN = fEverything)
+ OM_Supervisor <- RunModel(
+ sv,
+ RunOptions = RunOptions,
+ Param = ParamMichel
+ )
+ expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+test_that("RunModel.Supervisor with NA values in Qupstream", {
+ # Create Supervisor
+ InputsModel$`54057`$Qupstream[, c("R1", "R2")] <- NA
+ sv <- CreateSupervisor(InputsModel)
+ # Function to withdraw half of the measured flow
+ fWithdrawal <- function(y) { -y/2 }
+ # Function to release half of the the measured flow
+ fRelease <- function(y) { y/2 }
+ # Controller that withdraw half of the flow measured at node "54002" at location "R1"
+ CreateController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
+ # Controller that release half of the flow measured at node "54002" at location "R2"
+ CreateController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
+
+ OM_Supervisor <- RunModel(
+ sv,
+ RunOptions = RunOptions,
+ Param = ParamMichel
+ )
+ expect_equal(OM_Supervisor[["54057"]]$Qsim[1:3], rep(as.double(NA),3))
+ expect_equal(OM_Supervisor[["54057"]]$Qsim[4:length(IndPeriod_Run)],
+ OM_GriwrmInputs[["54057"]]$Qsim[4:length(IndPeriod_Run)])
+})
+
diff --git a/tests/testthat/test-calc_Qdiv.R b/tests/testthat/test-calc_Qdiv.R
new file mode 100644
index 0000000000000000000000000000000000000000..f4499f34f923e2dd989dbf3c25548037bfca0349
--- /dev/null
+++ b/tests/testthat/test-calc_Qdiv.R
@@ -0,0 +1,34 @@
+test_that("calc_Qdiv works", {
+ expect_equal(
+ calc_Qdiv(Qnat = 1, Qdiv = 1, Qmin = 0),
+ list(Qsim = 0, Qdiv = 1)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 5, Qdiv = 4, Qmin = 3),
+ list(Qsim = 3, Qdiv = 2)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 5, Qdiv = 10, Qmin = 3),
+ list(Qsim = 3, Qdiv = 2)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 5, Qdiv = 1, Qmin = 3),
+ list(Qsim = 4, Qdiv = 1)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 3, Qdiv = 1, Qmin = 4),
+ list(Qsim = 3, Qdiv = 0)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 3, Qdiv = -1, Qmin = 5),
+ list(Qsim = 4, Qdiv = -1)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = 3, Qdiv = -1, Qmin = 2),
+ list(Qsim = 4, Qdiv = -1)
+ )
+ expect_equal(
+ calc_Qdiv(Qnat = rep(1E6, 2), Qdiv = rep(1000, 2), Qmin = rep(1E12, 2)),
+ list(Qsim = rep(1E6, 2), Qdiv = rep(0, 2))
+ )
+})
diff --git a/tests/testthat/test-createGRiwrm.R b/tests/testthat/test-createGRiwrm.R
index 5b0e896fc0cb48171ec88138cfbfb973802402a2..8d34821d91cc64fd04649dd0e4e48e64fdff8773 100644
--- a/tests/testthat/test-createGRiwrm.R
+++ b/tests/testthat/test-createGRiwrm.R
@@ -19,29 +19,50 @@ H5920010 602213 2427449 43824.66 La Seine à Paris [Austerlitz après création
})
+# Setup a simple data.frame for GRiwrm
+nodes <- loadSevernNodes()
+
+test_that("Hydrological model nodes must have numeric area", {
+ nodes$area[nodes$id == "54057"] <- NA
+ expect_error(CreateGRiwrm(nodes),
+ regexp = "hydrological")
+})
+
+test_that("Duplicated nodes", {
+ nodes <- rbind(nodes, nodes[4,])
+ expect_error(CreateGRiwrm(nodes),
+ regexp = "Duplicated nodes detected")
+})
+
test_that("NA or Ungauged nodes at downstream should throw an error", {
- data(Severn)
- nodes <-
- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
- nodes$model <- "RunModel_GR4J"
- nodes$model[nodes$gauge_id == "54057"] <- "Ungauged"
- expect_error(CreateGRiwrm(
- nodes,
- list(
- id = "gauge_id",
- down = "downstream_id",
- length = "distance_downstream"
- )
- ),
- regexp = "downstream node")
+ nodes$model[nodes$id == "54057"] <- "Ungauged"
+ expect_error(CreateGRiwrm(nodes),
+ regexp = "downstream node")
nodes$model[nodes$gauge_id == "54057"] <- NA
- expect_error(CreateGRiwrm(
- nodes,
- list(
- id = "gauge_id",
- down = "downstream_id",
- length = "distance_downstream"
- )
- ),
- regexp = "downstream node")
+ expect_error(CreateGRiwrm(nodes),
+ regexp = "downstream node")
+})
+
+test_that("Diversion node", {
+ nodes <- rbind(nodes,
+ data.frame(id = "54029", down = "54002", length = 50, area = NA, model = "Diversion"))
+ expect_s3_class(CreateGRiwrm(nodes), "GRiwrm")
+ n99 <- nodes
+ n99$area[n99$model == "Diversion"] <- 99
+ expect_error(CreateGRiwrm(n99),
+ regexp = "Diversion node must have its area")
+ n_orphan <- nodes
+ n_orphan$id[n_orphan$model == "Diversion"] <- "54999"
+ expect_error(CreateGRiwrm(n_orphan),
+ regexp = "Diversion node must have the same `id` of")
+ n_samedown <- nodes
+ n_samedown$down[n_samedown$model == "Diversion"] <- "54032"
+ expect_error(CreateGRiwrm(n_samedown),
+ regexp = "downstream node of a Diversion node must be different")
+})
+
+test_that("Allow several downstream ends", {
+ nodes <- rbind(nodes,
+ data.frame(id = "54029", down = NA, length = NA, area = NA, model = "Diversion"))
+ expect_s3_class(CreateGRiwrm(nodes), "GRiwrm")
})
diff --git a/tests/testthat/test-getNodeRanking.R b/tests/testthat/test-getNodeRanking.R
new file mode 100644
index 0000000000000000000000000000000000000000..2fcf4dc1d072061d23b58329a4693bda7335b146
--- /dev/null
+++ b/tests/testthat/test-getNodeRanking.R
@@ -0,0 +1,26 @@
+nodes <- loadSevernNodes()
+
+test_that("Check ranking on Severn example", {
+ g <- CreateGRiwrm(nodes)
+ expect_equal(getNodeRanking(g),
+ c("54095", "54002", "54029", "54001", "54032", "54057"))
+})
+
+test_that("Check ranking with direct injection node", {
+ nodes$model[nodes$id == "54029"] <- NA
+ g <- CreateGRiwrm(nodes)
+ expect_equal(getNodeRanking(g),
+ c("54095", "54002", "54001", "54032", "54057"))
+})
+
+test_that("Check ranking with Diversion", {
+ n_div <- rbind(nodes, data.frame(id = "54029",
+ down = "54002",
+ length = 20,
+ model = "Diversion",
+ area = NA))
+ g <- CreateGRiwrm(n_div)
+ r <- getNodeRanking(g)
+ expect_lt(which(r == "54029"), which(r == "54002"))
+})
+
diff --git a/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
index e2e46c119142e3a9758df053890563b0bfeaf4f6..8325e86da6ed16f781a1c098624cd7c00c1a88b3 100644
--- a/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
+++ b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
@@ -255,10 +255,10 @@ for (x in 1:ncol(m))
As for the previous model, we need to set up an `GRiwrmInputsModel` object containing all the model inputs:
```{r}
-# A flow is needed for all direct injection nodes in the network
+# Flow time series are needed for all direct injection nodes in the network
# even if they may be overwritten after by a controller
-# The Qobs matrix is completed with 2 new columns for the new nodes
-QobsIrrig <- cbind(Qobs, Irrigation1 = 0, Irrigation2 = 0)
+QobsIrrig <- data.frame(Irrigation1 = rep(0, length(DatesR)),
+ Irrigation2 = rep(0, length(DatesR)))
# Creation of the GRiwrmInputsModel object
IM_Irrig <- CreateInputsModel(griwrmV04, DatesR, Precip, PotEvap, QobsIrrig)
diff --git a/vignettes/V05_Modelling_ungauged_nodes.Rmd b/vignettes/V05_Modelling_ungauged_nodes.Rmd
index 4aace25a8d5a3d82db4a3c42a14f2593b6f7c2cd..190fedfefd03af1d93c2f9d90d88933009655c99 100644
--- a/vignettes/V05_Modelling_ungauged_nodes.Rmd
+++ b/vignettes/V05_Modelling_ungauged_nodes.Rmd
@@ -117,13 +117,12 @@ PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
Precip <- ConvertMeteoSD(griwrmV05, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrmV05, PotEvapTot)
-Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
```
Then, the `GRiwrmInputsModel` object can be generated taking into account the new `GRiwrm` object:
```{r InputsModel}
-IM_U <- CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap, Qobs)
+IM_U <- CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap)
```
## Calibration of the model integrating ungauged nodes
@@ -141,6 +140,7 @@ IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)
RunOptions <- CreateRunOptions(IM_U,
IndPeriod_WarmUp = IndPeriod_WarmUp,
IndPeriod_Run = IndPeriod_Run)
+Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
InputsCrit <- CreateInputsCrit(IM_U,
FUN_CRIT = ErrorCrit_KGE2,
RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,],
diff --git a/vignettes/V06_Modelling_regulated_diversion.Rmd b/vignettes/V06_Modelling_regulated_diversion.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..d8a64b4bcc4d62835b73ae29e0c5a50dfa015b04
--- /dev/null
+++ b/vignettes/V06_Modelling_regulated_diversion.Rmd
@@ -0,0 +1,215 @@
+---
+title: "Severn_06: Modelling of a regulated diversion"
+output: rmarkdown::html_vignette
+vignette: >
+ %\VignetteIndexEntry{Severn_06: Modelling of a regulated diversion}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ fig.width = 6,
+ fig.asp = 0.68,
+ out.width = "70%",
+ fig.align = "center"
+)
+```
+
+```{r setup}
+library(airGRiwrm)
+data(Severn)
+```
+
+## The study case
+
+In this vignette, we are still using the study case of the vignette #1 and #2, and
+we are going to simulate a diversion at the node "54001" to provide flow to the node
+"54029" in order to support low flows in the Severn river.
+
+This objective is to maintain a minimum flow equals to the 3th decile of the flow
+at station "54029" without remaining less than the first decile of the flow downstream
+the station "54001".
+
+Let's compute the flow quantiles at each station in m3/s:
+
+```{r}
+Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+Qobs <- Qobs[, Severn$BasinsInfo$gauge_id]
+Qobs_m3s <- t(apply(Qobs, 1, function(r) r * Severn$BasinsInfo$area * 1E3 / 86400))
+apply(Qobs_m3s[, c("54029", "54001")], 2, quantile, probs = seq(0,1,0.1), na.rm = TRUE)
+```
+
+The rule to apply is expressed as follow:
+
+> The diversion from "54001" is computed to maintain a minimum flow of 5.3 m3/s at
+the gauging station "54029". The diversion is allowed as far as the flow at the
+gauging station "54001" is above 11.5 m3/s.
+
+## Network
+
+```{r}
+
+nodes_div <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes_div$model <- "RunModel_GR4J"
+nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54001",
+ downstream_id = "54029",
+ distance_downstream = 25,
+ model = "Diversion",
+ area = NA))
+renameCols <- list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
+griwrmV06 <- CreateGRiwrm(nodes_div, renameCols)
+plot(griwrmV06)
+```
+
+## GRiwrmInputsModel object
+
+We produce below the same operations as in the vignette "V02_Calibration_SD_model" to prepare the input data:
+
+```{r}
+data(Severn)
+nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$model <- "RunModel_GR4J"
+griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+```
+
+The main difference here is that we need to provide a diverted flow and a minimum
+remaining flow at station "54029". As, the diverted flow will be calculated during
+simulation, we provide a initial diverted flow equals to zero for all the time steps.
+
+```{r}
+Qdiv <- matrix(rep(0, length(DatesR)), ncol = 1)
+colnames(Qdiv) <- "54001"
+Qmin <- matrix(rep(11.5 * 86400, length(DatesR)), ncol = 1)
+colnames(Qmin) <- "54001"
+IM_div <- CreateInputsModel(griwrmV06, DatesR, Precip, PotEvap, Qobs = Qdiv, Qmin = Qmin)
+```
+
+## Implementation of the regulation controller
+
+### The supervisor
+
+The simulation is piloted through a `Supervisor` that can contain one or more
+`Controller`. The supervision time step will be here the same as the simulation
+time step: 1 day. Each day, the decision is taken for the current day from the
+measurement simulated the previous time step.
+
+```{r}
+sv <- CreateSupervisor(IM_div, TimeStep = 1L)
+```
+
+### The control logic function
+
+On a Diversion node, the minimum remaining flow is provided with the inputs and
+is automatically taken into account by the model. So the control logic function
+has only to take care about how much water to divert to the gauging station
+"54002".
+
+We need to enclose the logic function in a function factory providing the
+Supervisor in the environment of the function:
+
+```{r}
+#' @param sv the Supervisor environment
+logicFunFactory <- function(sv) {
+ #' @param Y Flow measured at "54002" the previous time step
+ function(Y) {
+ Qnat <- Y
+ # We need to remove the diverted flow to compute the natural flow at "54002"
+ lastU <- sv$controllers[[sv$controller.id]]$U
+ if (length(lastU) > 0) {
+ Qnat <- max(0, Y + lastU)
+ }
+ message(paste(round(c(Y, lastU, Qnat, -max(5.3 * 86400 - Qnat, 0))), collapse = ", "))
+ return(-max(5.3 * 86400 - Qnat, 0))
+ }
+}
+```
+
+### The controller
+
+We declare the controller which defines where is the measurement `Y` , where to
+apply the decision `U` with which logic function:
+
+```{r}
+CreateController(sv,
+ ctrl.id = "Low flow support",
+ Y = "54029",
+ U = "54001",
+ FUN = logicFunFactory(sv))
+```
+
+
+## Running the simulation
+
+First we need to create a `GRiwrmRunOptions` object and load the parameters calibrated in the vignette "V02_Calibration_SD_model":
+
+```{r}
+# Running simulation between 2002 and 2005
+IndPeriod_Run <- which(
+ DatesR >= as.POSIXct("2003-03-01", tz = "UTC") &
+ DatesR <= as.POSIXct("2004-01-01", tz = "UTC")
+)
+IndPeriod_WarmUp = seq(IndPeriod_Run[1] - 366,IndPeriod_Run[1] - 1)
+RunOptions <- CreateRunOptions(IM_div,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run)
+ParamV02 <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+```
+
+The node "54029" was initially an upstream node. As it receive water from "54001"
+it is no longer an upstream node and need a parameter for routing its upstream
+flows. We arbitrary say that the velocity in the channel between "54001" and
+"54029" is 1 m/s.
+
+```{r}
+ParamV02$`54029` <- c(1, ParamV02$`54029`)
+```
+
+And we run the supervised model:
+
+```{r}
+OM_div <- RunModel(sv, RunOptions = RunOptions, Param = ParamV02)
+```
+
+To compare results with diversion and without diversion, we run the model without
+the supervision (remember that we have set the diverted flow at zero in the inputs):
+
+```{r}
+OM_nat <- RunModel(IM_div, RunOptions = RunOptions, Param = ParamV02)
+```
+
+
+## Exploring results
+
+Let's plot the diverted flow, and compare the flow at stations 54029 and 54001,
+with and without low-flow support at station 54001:
+
+```{r, fig.asp=1}
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR,
+ Diverted_flow = OM_div$`54001`$Qdiv_m3 / 86400)
+
+oldpar <- par(mfrow=c(3,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001" and "54029"
+thresholds <- c("54001" = 11.5, "54029" = 5.3)
+lapply(names(thresholds), function(id) {
+ df <- cbind(attr(OM_div, "Qm3s")[, c("DatesR", id)],
+ attr(OM_nat, "Qm3s")[, id])
+ names(df) <- c("DatesR",
+ paste(id, "with low-flow support"),
+ paste(id, "natural flow"))
+ plot.Qm3s(df, ylim = c(0,50), lty = c("solid", "dashed"))
+ abline(h = thresholds[id], col = "green", lty = "dotted")
+})
+par(oldpar)
+```
+