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

Merge branch '52-comptabilite-cemaneige' into 'dev'

Resolve "Handle CemaNeige compatibility"

Closes #52

See merge request !22
parents d1eb8a60 f795bb5b
Pipeline #24734 failed with stages
in 14 minutes and 6 seconds
......@@ -5,7 +5,18 @@
#' @param Precip Matrix or data frame of numeric containing precipitation in mm. Column names correspond to node IDs
#' @param PotEvap Matrix or data frame of numeric containing potential evaporation in mm. Column names correspond to node IDs
#' @param Qobs Matrix or data frame of numeric containing potential observed flow in mm. Column names correspond to node IDs
#' @param ... further arguments passed to [airGR::CreateInputsModel]
#' @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.
#'
#' See [airGR::CreateInputsModel] documentation for details concerning each input.
#'
#' @return GRiwrmInputsModel object equivalent to **airGR** InputsModel object for a semi-distributed model (See [airGR::CreateInputsModel])
#' @export
......@@ -53,16 +64,60 @@
#' Qobs = Qobs)
#' str(InputsModels)
#'
CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
CreateInputsModel.GRiwrm <- function(x, DatesR,
Precip,
PotEvap = NULL,
Qobs,
PrecipScale = TRUE,
TempMean = NULL, TempMin = NULL,
TempMax = NULL, ZInputs = NULL,
HypsoData = NULL, NLayers = 5, ...) {
# Check and format inputs
varNames <- c("Precip", "PotEvap", "TempMean",
"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))) {
stop(sprintf(
"'%s' must have column names",
varName
))
} else if(!all(colnames(v) %in% x$id)) {
stop(sprintf(
"'%s' column names must be included in 'id's of the GRiwrm object",
varName
))
}
} else if (!varName %in% c("ZInputs", "NLayers")) {
stop(sprintf("'%s' must be a matrix or a data.frame", varName))
}
}
})
InputsModel <- CreateEmptyGRiwrmInputsModel(x)
Qobs[is.na(Qobs)] <- -99 # airGR::CreateInputsModel doesn't accept NA values
for(id in getNodeRanking(x)) {
message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
InputsModel[[id]] <- CreateOneGRiwrmInputsModel(
id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
)
InputsModel[[id]] <-
CreateOneGRiwrmInputsModel(id = id,
griwrm = x,
DatesR = DatesR,
Precip = getInputBV(Precip, id),
PrecipScale,
PotEvap = getInputBV(PotEvap, id),
TempMean = getInputBV(TempMean, id),
TempMin = getInputBV(TempMin, id),
TempMax = getInputBV(TempMax, id),
ZInputs = getInputBV(ZInputs, id),
HypsoData = getInputBV(HypsoData, id),
NLayers = getInputBV(NLayers, id, 5),
Qobs = Qobs
)
}
attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
return(InputsModel)
......@@ -74,6 +129,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
#' @param griwrm a `GRiwrm` object (See [CreateGRiwrm])
#'
#' @return \emph{GRiwrmInputsModel} empty object
#' @noRd
CreateEmptyGRiwrmInputsModel <- function(griwrm) {
InputsModel <- list()
class(InputsModel) <- c("GRiwrmInputsModel", class(InputsModel))
......@@ -90,9 +146,10 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#' @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 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.
CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs) {
#' @noRd
CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
node <- griwrm[griwrm$id == id,]
FUN_MOD <- griwrm$model[griwrm$id == id]
......@@ -111,15 +168,19 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
griwrm$area[griwrm$id %in% UpstreamNodes],
node$area - sum(griwrm$area[griwrm$id %in% UpstreamNodes], na.rm = TRUE)
)
if (BasinAreas[length(BasinAreas)] < 0) {
stop(sprintf(
"Area of the catchment %s must be greater than the sum of the areas of its upstream catchments",
id
))
}
names(BasinAreas) <- c(UpstreamNodes, id)
}
# Set model inputs with the **airGR** function
InputsModel <- CreateInputsModel(
FUN_MOD,
DatesR = DatesR,
Precip = Precip,
PotEvap = PotEvap,
...,
Qupstream = Qupstream,
LengthHydro = LengthHydro,
BasinAreas = BasinAreas
......@@ -149,7 +210,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
#' @param InputsModel a `GRiwrmInputsModel`
#'
#' @return A [numeric] representing the time step in seconds
#'
#' @noRd
getModelTimeStep <- function(InputsModel) {
TS <- sapply(InputsModel, function(x) {
if (inherits(x, "hourly")) {
......@@ -165,3 +226,32 @@ getModelTimeStep <- function(InputsModel) {
}
return(unique(TS))
}
#' Select the node input for input arguments of [airGR::CreateInputsModel]
#'
#' @param x [matrix] [data.frame] or named [vector] the input argument
#' @param id [character] the id of the node
#' @param unset default value if the id is not found in `x`
#'
#' @return the selected column or value in respect to `id`
#' @noRd
getInputBV <- function(x, id, unset = NULL) {
if(is.null(x)) {
return(unset)
}
if (is.matrix(x) || is.data.frame(x)) {
if (!id %in% colnames(x)) {
return(unset)
}
} else {
# vector (for ZInputs and NLayers)
if (length(x) == 1 && is.null(names(x))) {
return(x)
} else if(!id %in% names(x)) {
return(unset)
} else {
return(x[id])
}
}
return(x[, id])
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/CreateInputsModel.GRiwrm.R
\name{CreateEmptyGRiwrmInputsModel}
\alias{CreateEmptyGRiwrmInputsModel}
\title{Create an empty InputsModel object for \strong{airGRiwrm} nodes}
\usage{
CreateEmptyGRiwrmInputsModel(griwrm)
}
\arguments{
\item{griwrm}{a \code{GRiwrm} object (See \link{CreateGRiwrm})}
}
\value{
\emph{GRiwrmInputsModel} empty object
}
\description{
Create an empty InputsModel object for \strong{airGRiwrm} nodes
}
......@@ -4,7 +4,21 @@
\alias{CreateInputsModel.GRiwrm}
\title{Create InputsModel object for a \strong{airGRiwrm} network}
\usage{
\method{CreateInputsModel}{GRiwrm}(x, DatesR, Precip, PotEvap, Qobs, ...)
\method{CreateInputsModel}{GRiwrm}(
x,
DatesR,
Precip,
PotEvap = NULL,
Qobs,
PrecipScale = TRUE,
TempMean = NULL,
TempMin = NULL,
TempMax = NULL,
ZInputs = NULL,
HypsoData = NULL,
NLayers = 5,
...
)
}
\arguments{
\item{x}{GRiwrm object describing the diagram of the semi-distributed model (See \link{CreateGRiwrm})}
......@@ -17,7 +31,21 @@
\item{Qobs}{Matrix or data frame of numeric containing potential observed flow in mm. Column names correspond to node IDs}
\item{...}{further arguments passed to \link[airGR:CreateInputsModel]{airGR::CreateInputsModel}}
\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{TempMean}{(optional) \link{matrix} or \link{data.frame} of time series of mean air temperature [°C], required 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{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{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{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{
GRiwrmInputsModel object equivalent to \strong{airGR} InputsModel object for a semi-distributed model (See \link[airGR:CreateInputsModel]{airGR::CreateInputsModel})
......@@ -25,6 +53,11 @@ GRiwrmInputsModel object equivalent to \strong{airGR} InputsModel object for a s
\description{
Create 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.
See \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} documentation for details concerning each input.
}
\examples{
#################################################################
# Run the `airGRRunModel_Lag` example in the GRiwrm fashion way #
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/CreateInputsModel.GRiwrm.R
\name{CreateOneGRiwrmInputsModel}
\alias{CreateOneGRiwrmInputsModel}
\title{Create one InputsModel for a \strong{airGRiwrm} node}
\usage{
CreateOneGRiwrmInputsModel(id, griwrm, DatesR, Precip, PotEvap, Qobs)
}
\arguments{
\item{id}{string of the node identifier}
\item{griwrm}{See \link{CreateGRiwrm})}
\item{DatesR}{vector of dates required to create the GR model and CemaNeige module inputs}
\item{Precip}{time series of potential evapotranspiration (catchment average) (mm/time step)}
\item{PotEvap}{time series of potential evapotranspiration (catchment average) (mm/time step)}
\item{Qobs}{Matrix or data frame of numeric containing observed flow (mm/time step). Column names correspond to node IDs}
}
\value{
\emph{InputsModel} object for one.
}
\description{
Create one InputsModel for a \strong{airGRiwrm} node
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/CreateInputsModel.GRiwrm.R
\name{getModelTimeStep}
\alias{getModelTimeStep}
\title{Check time steps of the model of all the nodes and return the time step in seconds}
\usage{
getModelTimeStep(InputsModel)
}
\arguments{
\item{InputsModel}{a \code{GRiwrmInputsModel}}
}
\value{
A \link{numeric} representing the time step in seconds
}
\description{
This function is called inside \link{CreateInputsModel.GRiwrm} for defining the time step of the big model.
}
#' Set up a data sample for Cemaneige tests
#'
#' A basin of 3 nodes with 2 identical upstream basins and a nul area downstream basin with a lag of 0 time step
#'
#' @return [list] with:
#' - `BasinInfo`
#' - `griwrm`
#' - `Precip`
#' - `TempMean`
#' - `PotEvap`
#' - `Qobs`
#' - `HypsoData`
#' @noRd
#'
#' @examples
setUpCemaNeigeData <- function() {
library(airGR)
data(L0123001)
detach("package:airGR")
# 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
ids <- c("Up1", "Up2", "Down")
l <- lapply(c("P", "T", "E", "Qmm"), function(x) {
m <- matrix(data = rep(BasinObs[[x]], 3), ncol = 3)
colnames(m) <- ids
return(m)
})
l$DatesR <- BasinObs$DatesR
names(l) <- c("Precip", "TempMean", "PotEvap", "Qobs")
l$HypsoData <- matrix(data = rep(BasinInfo$HypsoData, 3), ncol = 3)
colnames(l$HypsoData) <- ids
l$ZInputs <- rep(BasinInfo$HypsoData[51], 3)
names(l$ZInputs) <- ids
db <- data.frame(id = ids,
length = c(0, 0, NA),
down = c("Down", "Down", NA),
area = c(rep(BasinInfo$BasinArea, 2), BasinInfo$BasinArea * 2),
model = rep("RunModel_CemaNeigeGR4J", 3),
stringsAsFactors = FALSE)
# Create GRiwrm object from the data.frame
l$griwrm <- CreateGRiwrm(db)
l$DatesR <- BasinObs$DatesR
return(l)
}
context("CreateInputsModel")
l <- setUpCemaNeigeData()
test_that("CemaNeige data should be in InputsModel", {
InputsModels <- CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs)
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_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)))
})
})
test_that("downstream sub-catchment area should be positive", {
l$griwrm$area[3] <- 360
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
})
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$ZInputs <- l$ZInputs[1:2]
l$TempMean <- l$TempMean[,1:2]
l$HypsoData <- l$HypsoData[,1:2]
InputsModels <- CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs)
expect_false(inherits(InputsModels$Down, "CemaNeige"))
expect_null(InputsModels$Down$LayerPrecip)
})
test_that("throws error on wrong column name", {
colnames(l$Precip)[1] <- "Up0"
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
colnames(l$Precip) <- NULL
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs))
})
test_that("throws error when missing CemaNeige data", {
expect_error(CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
Qobs = l$Qobs))
})
......@@ -12,8 +12,9 @@ griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", leng
# InputsModel
DatesR <- Severn$BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
BasinsObs <- Severn$BasinsObs[griwrm$id]
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
......@@ -97,3 +98,41 @@ test_that("RunModelSupervisor with multi time steps controller, two regulations
)
expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
})
test_that("RunModel.GRiwrm handles CemaNeige", {
l <- setUpCemaNeigeData()
l$griwrm[l$griwrm$id == "Down", "model"] <- "RunModel_GR4J"
l$TempMean <- l$TempMean[,1:2]
l$ZInputs <- l$ZInputs[1:2]
l$TempMean <- l$TempMean[,1:2]
l$HypsoData <- l$HypsoData[,1:2]
InputsModels <- CreateInputsModel(l$griwrm,
DatesR = l$DatesR,
Precip = l$Precip,
PotEvap = l$PotEvap,
TempMean = l$TempMean,
ZInputs = l$ZInputs,
HypsoData = l$HypsoData,
Qobs = l$Qobs)
## 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"))
## preparation of the RunOptions object
RunOptions <- suppressWarnings(CreateRunOptions(InputsModels,
IndPeriod_Run = Ind_Run))
ids <- l$griwrm$id
names(ids) <- ids
Params <- lapply(ids, function(x) {
c(X1 = 408.774, X2 = 2.646, X3 = 131.264, X4 = 1.174,
CNX1 = 0.962, CNX2 = 2.249)
})
Params$Down <- c(1, Params$Down[1:4])
OutputsModel <- RunModel(
InputsModels,
RunOptions = RunOptions,
Param = Params
)
expect_named(OutputsModel, l$griwrm$id)
Qm3s <- attr(OutputsModel, "Qm3s")
expect_equal(Qm3s[,4], rowSums(Qm3s[,2:3]))
})
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