Commit b80336bd authored by Dorchies David's avatar Dorchies David
Browse files

Merge branch '60-createinputsmodel-make-qobs-parameter-optional' into 'dev'

Resolve "[CreateInputsModel] Make `Qobs` parameter optional"

Closes #60

See merge request !28
parents 9c199919 32aa9e51
Pipeline #29996 passed with stage
in 5 minutes and 29 seconds
......@@ -2,9 +2,9 @@
#'
#' @param x \[GRiwrm object\] diagram of the semi-distributed model (See [CreateGRiwrm])
#' @param DatesR [POSIXt] vector of dates
#' @param Precip [matrix] or [data.frame] frame of numeric containing precipitation in \[mm per time step\]. Column names correspond to node IDs
#' @param PotEvap [matrix] or [data.frame] frame of numeric containing potential evaporation \[mm per time step\]. Column names correspond to node IDs
#' @param Qobs [matrix] or [data.frame] frame of numeric containing observed flows in \[mm per time step\]. Column names correspond to node IDs
#' @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
......@@ -50,11 +50,9 @@
#' PotEvap <- matrix(BasinObs$E, ncol = 1)
#' colnames(PotEvap) <- "GaugingDown"
#'
#' # Observed flows are integrated now because we mix:
#' # - flows that are directly injected in the model
#' # - flows that could be used for the calibration of the hydrological models
#' Qobs = matrix(c(Qupstream, BasinObs$Qmm), ncol = 2)
#' colnames(Qobs) <- griwrm$id
#' # Observed flows should at least contains flows that are directly injected in the model
#' Qobs = matrix(Qupstream, ncol = 1)
#' colnames(Qobs) <- "Reservoir"
#' str(Qobs)
#'
#' InputsModels <- CreateInputsModel(griwrm,
......@@ -65,41 +63,73 @@
#' str(InputsModels)
#'
CreateInputsModel.GRiwrm <- function(x, DatesR,
Precip,
Precip = NULL,
PotEvap = NULL,
Qobs,
Qobs = NULL,
PrecipScale = TRUE,
TempMean = NULL, TempMin = NULL,
TempMax = NULL, ZInputs = NULL,
HypsoData = NULL, NLayers = 5, ...) {
# Check and format inputs
varNames <- c("Precip", "PotEvap", "TempMean",
varNames <- c("Precip", "PotEvap", "TempMean", "Qobs",
"TempMin", "TempMax", "ZInputs", "HypsoData", "NLayers")
names(varNames) <- varNames
lapply(varNames, function(varName) {
v <- get(varName)
if(!is.null(v)) {
if(is.matrix(v) || is.data.frame(v)) {
if(is.null(colnames(v))) {
if (!is.null(v)) {
if (is.matrix(v) || is.data.frame(v)) {
if (is.null(colnames(v))) {
stop(sprintf(
"'%s' must have column names",
varName
))
} else if(!all(colnames(v) %in% x$id)) {
} else if (!all(colnames(v) %in% x$id)) {
stop(sprintf(
"'%s' column names must be included in 'id's of the GRiwrm object",
varName
))
}
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)
}
} 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)]
if (length(directFlowIds) > 0) {
err <- FALSE
if (is.null(Qobs)) {
err <- TRUE
} else {
Qobs <- as.matrix(Qobs)
if (is.null(colnames(Qobs))) {
err <- TRUE
} else if (!all(directFlowIds %in% colnames(Qobs))) {
err <- TRUE
}
}
if (err) stop(sprintf("'Qobs' column names must at least contain %s", paste(directFlowIds, collapse = ", ")))
}
InputsModel <- CreateEmptyGRiwrmInputsModel(x)
# Qobs completion
Qobs0 <- matrix(0, nrow = length(DatesR), ncol = nrow(x))
colnames(Qobs0) <- x$id
if (is.null(Qobs)) {
Qobs <- Qobs0
} else {
Qobs0[, colnames(Qobs)] <- Qobs
Qobs <- Qobs0
}
for(id in getNodeRanking(x)) {
message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
InputsModel[[id]] <-
......@@ -141,9 +171,10 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#'
#' @param id string of the node identifier
#' @param griwrm See [CreateGRiwrm])
#' @param DatesR vector of dates required to create the GR model and CemaNeige module inputs
#' @param Precip time series of potential evapotranspiration (catchment average) (mm/time step)
#' @param PotEvap time series of potential evapotranspiration (catchment average) (mm/time step)
#' @param ... parameters sent to [airGR::CreateInputsModel]:
#' - `DatesR` [vector] of dates required to create the GR model and CemaNeige module inputs
#' - `Precip` [vector] time series of potential evapotranspiration (catchment average) (mm/time step)
#' - `PotEvap` [vector] time series of potential evapotranspiration (catchment average) (mm/time step)
#' @param Qobs Matrix or data frame of numeric containing observed flow (mm/time step). Column names correspond to node IDs
#'
#' @return \emph{InputsModel} object for one.
......
......@@ -7,9 +7,9 @@
\method{CreateInputsModel}{GRiwrm}(
x,
DatesR,
Precip,
Precip = NULL,
PotEvap = NULL,
Qobs,
Qobs = NULL,
PrecipScale = TRUE,
TempMean = NULL,
TempMin = NULL,
......@@ -25,11 +25,11 @@
\item{DatesR}{\link{POSIXt} vector of dates}
\item{Precip}{\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} frame of numeric containing precipitation in [mm per time step]. Column names correspond to node IDs}
\item{PotEvap}{\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} frame of numeric containing potential evaporation [mm per time step]. Column names correspond to node IDs}
\item{Qobs}{\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} frame of numeric containing observed flows in [mm per time step]. Column names correspond to node IDs}
\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)}
......@@ -88,11 +88,9 @@ colnames(Precip) <- "GaugingDown"
PotEvap <- matrix(BasinObs$E, ncol = 1)
colnames(PotEvap) <- "GaugingDown"
# Observed flows are integrated now because we mix:
# - flows that are directly injected in the model
# - flows that could be used for the calibration of the hydrological models
Qobs = matrix(c(Qupstream, BasinObs$Qmm), ncol = 2)
colnames(Qobs) <- griwrm$id
# Observed flows should at least contains flows that are directly injected in the model
Qobs = matrix(Qupstream, ncol = 1)
colnames(Qobs) <- "Reservoir"
str(Qobs)
InputsModels <- CreateInputsModel(griwrm,
......
......@@ -10,17 +10,16 @@ test_that("CemaNeige data should be in InputsModel", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs)
HypsoData = l$HypsoData)
)
l$DatesR <- as.data.frame(l$DatesR)
lapply(InputsModels, function(IM) {
lapply(c("DatesR", "Precip", "PotEvap"), function(varName) {
expect_equal(IM[[varName]], l[[varName]][,1])
expect_equal(IM[[varName]], l[[varName]][, 1])
})
expect_named(IM$LayerPrecip, paste0("L", seq(1,5)))
expect_named(IM$LayerTempMean, paste0("L", seq(1,5)))
expect_named(IM$LayerFracSolidPrecip, paste0("L", seq(1,5)))
expect_named(IM$LayerPrecip, paste0("L", seq(1, 5)))
expect_named(IM$LayerTempMean, paste0("L", seq(1, 5)))
expect_named(IM$LayerFracSolidPrecip, paste0("L", seq(1, 5)))
})
})
......@@ -32,16 +31,16 @@ test_that("downstream sub-catchment area should be positive", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
HypsoData = l$HypsoData),
regexp = "must be greater than the sum of the areas")
})
test_that("handles mix of with and without CemaNeige nodes", {
l$griwrm[l$griwrm$id == "Down", "model"] <- "RunModel_GR4J"
l$TempMean <- l$TempMean[,1:2]
l$TempMean <- l$TempMean[, 1:2]
l$ZInputs <- l$ZInputs[1:2]
l$TempMean <- l$TempMean[,1:2]
l$HypsoData <- l$HypsoData[,1:2]
l$TempMean <- l$TempMean[, 1:2]
l$HypsoData <- l$HypsoData[, 1:2]
InputsModels <- suppressWarnings(
CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
......@@ -49,8 +48,7 @@ test_that("handles mix of with and without CemaNeige nodes", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs)
HypsoData = l$HypsoData)
)
expect_false(inherits(InputsModels$Down, "CemaNeige"))
expect_null(InputsModels$Down$LayerPrecip)
......@@ -64,8 +62,8 @@ test_that("throws error on wrong column name", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
HypsoData = l$HypsoData),
regexp = "column names must be included in")
colnames(l$Precip) <- NULL
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
......@@ -73,14 +71,72 @@ test_that("throws error on wrong column name", {
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
HypsoData = l$HypsoData),
regexp = "must have column names")
})
test_that("throw error on missing column in inputs", {
l$Precip <- l$Precip[, -1]
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData),
regexp = "Precip is missing")
})
test_that("throw error on wrong number of rows in inputs", {
l$Precip <- l$Precip[-1, ]
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData),
regexp = "number of rows and the length of 'DatesR' must be equal")
})
test_that("throws error when missing CemaNeige data", {
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap),
regexp = "'TempMean' is missing")
})
test_that("throws error when missing Qobs on node not related to an hydrological model", {
l$griwrm$model[1] <- NA
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap),
regexp = "'Qobs' column names must at least contain")
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
Qobs = l$Qobs))
Qobs = l$Qobs[, -1]),
regexp = "'Qobs' column names must at least contain")
})
test_that("must works with node not related to an hydrological model", {
l$griwrm$model[1] <- NA
IM <- suppressWarnings(CreateInputsModel(
l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
Qobs = l$Qobs[, 1, drop = FALSE],
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData
))
expect_equal(IM[[2]]$Qupstream[, "Up1"], l$Qobs[, "Up1"] * l$griwrm[1, "area"] * 1E3)
expect_equal(colnames(IM[[2]]$Qupstream), c("Up1", "Up2"))
})
......@@ -109,7 +109,7 @@ The **airGR** `CreateInputsModel` function is extended in order to handle the `G
```{r}
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap)
```
# References
......@@ -47,7 +47,7 @@ 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)
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap)
str(InputsModel)
```
......
......@@ -69,15 +69,21 @@ 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)
Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
```
This time, we need to provide observed flows as inputs for the nodes '54002' and '54095':
```{r}
QobsInputs <- Qobs[, c("54002", "54095")]
```
The `GRiwrmInputsModel` object can be generated taking into account the new `GRiwrm` object:
Then, the `GRiwrmInputsModel` object can be generated taking into account the new `GRiwrm` object:
```{r}
IM_OL <- CreateInputsModel(griwrmV03, DatesR, Precip, PotEvap, Qobs)
IM_OL <- CreateInputsModel(griwrmV03, DatesR, Precip, PotEvap, QobsInputs)
```
## Calibration of the new model
......
......@@ -231,10 +231,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 the nodes in the network
# even if it is overwritten after by a controller
# A flow is 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 <- cbind(Qobs[, c("54002", "54095")], Irrigation1 = 0, Irrigation2 = 0)
# Creation of the GRiwrmInputsModel object
IM_Irrig <- CreateInputsModel(griwrmV04, DatesR, Precip, PotEvap, QobsIrrig)
......
......@@ -80,7 +80,7 @@ The **airGR** CreateInputsModel function is extended in order to handle the griw
```{r CreateInputsModel}
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qnat)
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap)
```
......
......@@ -195,7 +195,6 @@ Then, we load observation data for each reservoir and each connection:
library(seinebasin)
data(QNAT)
data(Qreservoirs)
Qnat <- cbind(Qnat, Qreservoirs)
```
## How to handle online reservoir? The Pannecière lake case.
......@@ -208,15 +207,15 @@ There are two possibilities:
If we know in advance the flow released by the reservoir, upstream flow informations is not usefull for the simulation. But reservoir management simulation would need upstream flow informations for simulating the reservoir state during simulation. The second alternative will be useful for the next phases.
```{r Qnat}
Qnat <- cbind(Qnat, -Qnat[,"CHAUM_07"] * griwrm2$area[griwrm2$id == "CHAUM_07"] * 1e3)
colnames(Qnat)[ncol(Qnat)] <- "PANNEC_P"
Qreservoirs <- cbind(Qreservoirs, -Qnat[,"CHAUM_07"] * griwrm2$area[griwrm2$id == "CHAUM_07"] * 1e3)
colnames(Qreservoirs)[ncol(Qreservoirs)] <- "PANNEC_P"
```
## Create the InputsModel object
```{r CreateInputsModel}
InputsModel2 <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qnat)
InputsModel2 <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qreservoirs)
```
## Run simulation with naturalised flow parameters
......
......@@ -48,13 +48,13 @@ data(QOBS)
iEnd <-which(DatesR == as.POSIXct("2008-07-31", tz = "UTC"))
data(Qreservoirs)
Qobs3 <- cbind(Qobs, Qreservoirs[1:iEnd,])
QresMarne <- Qreservoirs[1:iEnd, grep("MARNE", colnames(Qreservoirs))]
id_GR_nodes <- griwrm3$id[!is.na(griwrm3$model)]
InputsModel3 <- CreateInputsModel(griwrm3,
DatesR[1:iEnd],
Precip[1:iEnd, id_GR_nodes],
PotEvap[1:iEnd, id_GR_nodes],
Qobs3[, griwrm3$id])
QresMarne)
```
## GriwmRunOptions object
......@@ -89,7 +89,7 @@ RunOptions <- CreateRunOptions(
InputsCrit <- CreateInputsCrit(
InputsModel = InputsModel3,
FUN_CRIT = airGR::ErrorCrit_KGE2,
RunOptions = RunOptions, Obs = Qobs3[IndPeriod_Run,]
RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,]
)
```
......@@ -124,7 +124,7 @@ OutputsModels3 <- RunModel(
htmltools::tagList(lapply(
griwrm3$id[!is.na(griwrm3$model)],
function(x) {
Q3 <- Qobs3[RunOptions[[1]]$IndPeriod_Run,x]
Q3 <- Qobs[RunOptions[[1]]$IndPeriod_Run,x]
iQ3 <- which(!is.na(Q3))
IndPeriod_Obs <- iQ3[1]:tail(iQ3,1)
OutputsModels <- ReduceOutputsModel(OutputsModels3[[x]], IndPeriod_Obs)
......
......@@ -49,13 +49,13 @@ data(QOBS)
iEnd <-which(DatesR == as.POSIXct("2008-07-31", tz = "UTC"))
data(Qreservoirs)
Qobs3 <- cbind(Qobs, Qreservoirs[1:iEnd,])
QresMarne <- Qreservoirs[1:iEnd, grep("MARNE", colnames(Qreservoirs))]
id_GR_nodes <- griwrm3$id[!is.na(griwrm3$model)]
InputsModel3 <- CreateInputsModel(griwrm3,
DatesR[1:iEnd],
Precip[1:iEnd, id_GR_nodes],
PotEvap[1:iEnd, id_GR_nodes],
Qobs3[, griwrm3$id])
InputsModel3 <- CreateInputsModel(griwrm3,
DatesR[1:iEnd],
Precip[1:iEnd, id_GR_nodes],
PotEvap[1:iEnd, id_GR_nodes],
QresMarne)
```
## GriwmRunOptions object
......@@ -90,7 +90,7 @@ RunOptions <- CreateRunOptions(
InputsCrit <- CreateInputsCrit(
InputsModel = InputsModel3,
FUN_CRIT = airGR::ErrorCrit_KGE2,
RunOptions = RunOptions, Obs = Qobs3[IndPeriod_Run,]
RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,]
)
```
......@@ -125,7 +125,7 @@ OutputsModels3 <- RunModel(
htmltools::tagList(lapply(
griwrm3$id[!is.na(griwrm3$model)],
function(x) {
Q3 <- Qobs3[RunOptions[[1]]$IndPeriod_Run,x]
Q3 <- Qobs[RunOptions[[1]]$IndPeriod_Run,x]
iQ3 <- which(!is.na(Q3))
IndPeriod_Obs <- iQ3[1]:tail(iQ3,1)
OutputsModels <- ReduceOutputsModel(OutputsModels3[[x]], IndPeriod_Obs)
......
......@@ -59,8 +59,7 @@ data(QNAT)
InputsModel4 <- CreateInputsModel(griwrm4,
DatesR,
Precip[, selectedNodes],
PotEvap[, selectedNodes],
Qnat[, selectedNodes])
PotEvap[, selectedNodes])
```
## GriwmRunOptions object
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment