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

refactor: Remove Girop object dependencies

Closes #9
Showing with 39 additions and 90 deletions
+39 -90
......@@ -18,7 +18,6 @@ export(CreateInputsCrit)
export(CreateInputsModel)
export(CreateRunOptions)
export(Ginet)
export(Girop)
export(RunModel)
export(getNodeRanking)
import(airGR)
#' Create InputsModel object for a GRIWRM network
#'
#' @param x Ginet object describing the diagram of the semi-distributed model, see \code{[Ginet]}.
#' @param girop Girop object giving the run-off model parameters, see \code{[Girop]}.
#' @param DateR Vector of POSIXlt observation time steps.
#' @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.
......@@ -11,7 +10,7 @@
#'
#' @return GriwrmInputsModel object equivalent to airGR InputsModel object for a semi-distributed model (See \code{\link[airGR]{CreateInputsModel}})
#' @export
CreateInputsModel.Griwrm <- function(x, girop, DatesR, Precip, PotEvap, Qobs, verbose = TRUE, ...) {
CreateInputsModel.Griwrm <- function(x, DatesR, Precip, PotEvap, Qobs, verbose = TRUE, ...) {
InputsModel <- CreateEmptyGriwrmInputsModel()
Qobs[is.na(Qobs)] <- -99 # airGRCreateInputsModel doesn't accept NA values
......@@ -19,7 +18,7 @@ CreateInputsModel.Griwrm <- function(x, girop, DatesR, Precip, PotEvap, Qobs, ve
for(id in getNodeRanking(x)) {
if(verbose) cat("CreateInputsModel.griwrm: Treating sub-basin", id, "...\n")
InputsModel[[id]] <- CreateOneGriwrmInputsModel(
id, x, girop, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
)
}
return(InputsModel)
......@@ -40,16 +39,15 @@ CreateEmptyGriwrmInputsModel <- function() {
#'
#' @param id string of the node identifier
#' @param ginet See \code{[Ginet]}.
#' @param girop See \code{[Girop]}.
#' @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 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, ginet, girop, DatesR, Precip, PotEvap, Qobs) {
CreateOneGriwrmInputsModel <- function(id, ginet, DatesR, Precip, PotEvap, Qobs) {
node <- ginet[ginet$id == id,]
FUN_MOD <- girop$model[girop$id == id]
FUN_MOD <- ginet$model[ginet$id == id]
# Set hydraulic parameters
UpstreamNodes <- ginet$id[ginet$down == id & !is.na(ginet$down)]
......@@ -67,10 +65,10 @@ CreateOneGriwrmInputsModel <- function(id, ginet, girop, DatesR, Precip, PotEvap
Qupstream <- cbind(Qupstream, Qupstream1)
}
}
LengthHydro <- ginet$length[girop$id %in% UpstreamNodes]
LengthHydro <- ginet$length[ginet$id %in% UpstreamNodes]
BasinAreas <- c(
girop$area[girop$id %in% UpstreamNodes],
girop$area[girop$id == id] - sum(girop$area[girop$id %in% UpstreamNodes])
ginet$area[ginet$id %in% UpstreamNodes],
node$area - sum(ginet$area[ginet$id %in% UpstreamNodes])
)
}
......
......@@ -2,13 +2,13 @@
#'
#' @param InputsModel object of class \emph{GriwrmInputsModel}, see \code{[CreateInputsModel.Griwrm]} for details.
#' @param RunOptions object of class \emph{GriwrmRunOptions}, see \code{[CreateRunOptions.Griwrm]} for details.
#' @param girop Girop object giving the run-off model parameters, see \code{[Girop]}.
#' @param Param list of parameter. The list item names are the IDs of the sub-basins. Each item is a vector of numerical parameters.
#' @param verbose (optional) boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}
#' @param ... Mandatory for S3 method signature function compatibility with generic.
#'
#' @return \emph{GriwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See \code{\link[airGR]{RunModel}}) for each node of the semi-distributed model.
#' @export
RunModel.GriwrmInputsModel <- function(InputsModel, RunOptions, girop, verbose = TRUE, ...) {
RunModel.GriwrmInputsModel <- function(InputsModel, RunOptions, Param, verbose = TRUE, ...) {
OutputsModel <- list()
class(OutputsModel) <- append(class(OutputsModel), "GriwrmOutputsModel")
......@@ -23,7 +23,7 @@ RunModel.GriwrmInputsModel <- function(InputsModel, RunOptions, girop, verbose =
OutputsModel[[IM$id]] <- RunModel(
InputsModel = IM,
RunOptions = RunOptions[[IM$id]],
Param = unlist(girop$params[girop$id == IM$id])
Param = Param[[IM$id]]
)
}
......
......@@ -7,8 +7,8 @@
#'
#' @return `Ginet` class object containing the description of diagram of the semi-distributed catchment model
#' @export
Ginet <- function(db, cols = list(id = "id", down = "down", length = "length", runoff = "runoff"), keep_all = FALSE) {
colsDefault <- list(id = "id", down = "down", length = "length", runoff = "runoff")
Ginet <- function(db, cols = list(id = "id", down = "down", length = "length", model = "model"), keep_all = FALSE) {
colsDefault <- list(id = "id", down = "down", length = "length", model = "model", area = "area")
cols <- utils::modifyList(colsDefault, as.list(cols))
db <- dplyr::rename(db, unlist(cols))
if(!keep_all) {
......
#' Generate the list of run-off models and their parameters
#'
#' @param db data frame containing at least the id the area and the model of the sub-basin.
#' @param cols named list or vector for matching columns of `db` parameter. By default, mandatory columns names are: `id`, `area`, `model`. But other names can be handled with a named list or vector containing items defined as `"required name" = "column name in db"`.
#' @param keep_all Keep all column of `db` or keep only columns defined in `cols`
#'
#' @return \emph{Girop} object.
#' @export
Girop <- function(db, cols = c(id = "id", area = "area", model = "model", params = "params"), keep_all = FALSE) {
colsDefault <- list(id = "id", area = "area", model = "model", params = "params")
cols <- utils::modifyList(colsDefault, as.list(cols))
if(!any(names(db) == cols$params)) {
# Add missing params column in the database
db[[cols$params]] = NA
}
db <- dplyr::rename(db, unlist(cols))
if(!keep_all) {
db <- dplyr::select(db, names(cols))
}
class(db) <- append(class(db), "Girop")
db
}
......@@ -30,52 +30,27 @@ seine_nodes
Create the ginet object which lists the nodes and describes the network diagram. It's a dataframe of class `Ginet` and `Griwrm` with specific column names:
- `id`: the identifier of the node in the network
- `down`: the identifier of the next node downstream
- `length`: hydraulic distance to the next downstream node
- `runoff`: does the node is a rainfall run-off model?
- `id`: the identifier of the node in the network.
- `down`: the identifier of the next hydrological node downstream.
- `length`: hydraulic distance to the next hydrological downstream node.
- `model`: Name of the hydrological model used (E.g. "RunModel_GR4J"). `NA` for other type of node.
- `area`: Area of the sub-catchment (km<sup>2</sup>). Used for hydrological model such as GR models. `NA` if not used.
`Ginet` function helps to rename the columns of the dataframe and assign the variable classes.
```{r}
# Specify that all nodes are of run-off type
seine_nodes$runoff <- TRUE
# Convert distance in km as it the unit used by airGR
seine_nodes$length <- seine_nodes$distance_aval / 1000
seine_nodes$model <- "RunModel_GR4J"
# Generate the ginet object
ginet <- Ginet(seine_nodes, list(id = "id_sgl", down = "id_aval"))
ginet <- Ginet(seine_nodes, list(id = "id_sgl", down = "id_aval", length = "distance_aval"))
ginet
```
Each line of the `Ginet` object having the `runoff` columns switched to `TRUE` should have a corresponding line in the `Girop` object which contains the parameters of the rainfall run-off models.
The `Girop` object is a dataframe of class `Girop` with specific columns:
- `id`: the identifier of the node in the network
- `area`: the total area of the basin (including upstream sub-basins) at the location of the node (km<sup>2</sup>)
- `model`: the name of the rainfall run-off model used (e.g. "RunModel_GR4J")
- `params`: a list containing the calibration parameters of the model
```{r}
# Specify which run-off model to use
seine_nodes$model = "RunModel_GR4J"
# Generate girop object
girop <- Girop(seine_nodes, list(id = "id_sgl", area = "area"))
girop
```
## Observation time series
Loading hydrometeorological data on the Seine river basin from the ClimAware project:
```{r, warning=FALSE, message=FALSE}
urls <-
file.path(
"https://stratus.irstea.fr/d/0b18e688851a45478f7a/files/?p=/climaware_hydro/Q_OBS_NAT",
paste0(ginet$id, "_NAT.txt&dl=1")
)
names(urls) <- ginet$id
Precip <- NULL
PotEvap <- NULL
......@@ -105,8 +80,8 @@ for(id in ginet$id) {
Precip <- MergeTS(Precip, id, ts[,c("Date", "Ptot")])
# ETP column is merged into PotEvap dataframe
PotEvap <- MergeTS(PotEvap, id, ts[,c("Date", "ETP")])
# Convert Qobs from m3/s to mm
ts$Qnat <- ts$Qnat * 86.4 / girop$area[girop$id == id]
# Convert Qobs from m3/s to mm/time step
ts$Qnat <- ts$Qnat * 86.4 / ginet$area[ginet$id == id]
# Setting data gaps to NA
ts$Qnat[ts$Qnat <= 0] <- NA
# Qnat column is merged into Qobs dataframe
......@@ -126,7 +101,7 @@ The airGR CreateInputsModel function is extended in order to handle the ginet ob
```{r}
InputsModel <- CreateInputsModel(ginet, girop, DatesR, Precip, PotEvap, Qobs)
InputsModel <- CreateInputsModel(ginet, DatesR, Precip, PotEvap, Qobs)
```
......@@ -136,6 +111,6 @@ InputsModel <- CreateInputsModel(ginet, girop, DatesR, Precip, PotEvap, Qobs)
```{r}
dir.create("_cache", showWarnings = FALSE)
save(ginet, girop, Qobs, InputsModel, file = "_cache/V01.RData")
save(ginet, Qobs, InputsModel, file = "_cache/V01.RData")
```
......@@ -53,22 +53,19 @@ save(OutputsCalib, file = "_cache/V03.RData")
## Run model with Michel calibration
```{r}
giropMichel <- girop
for(id in giropMichel$id) {
giropMichel$params[giropMichel$id == id] <- list(OutputsCalib[[id]]$Param)
}
ParamMichel <- sapply(ginet$id, function(x) {OutputsCalib[[x]]$Param})
OutputsModels <- RunModel(
InputsModel = InputsModel,
RunOptions = RunOptions,
girop = giropMichel
Param = ParamMichel
)
```
## Save calibration data for next vignettes
```{r}
save(giropMichel, file = "_cache/V03.RData")
save(ParamMichel, file = "_cache/V03.RData")
```
......
......@@ -49,20 +49,22 @@ This lag parameter has to be converted in a speed in m/s used in the airGR lag m
# Convert TGR routing parameter into speed
params <- merge(ginet, ClimAwareParams, by.x = "id", by.y = "id_sgl")
for(id in girop$id) {
# Maximum connection length with upstream nodes
ParamClimAware <- sapply(ginet$id, function(id) {
nodeParam <- ClimAwareParams[ClimAwareParams$id_sgl == id,]
# Record hydrological model parameters
Param <- unlist(nodeParam[c("S", "IGF", "KR", "T")])
# Add lag model parameter if upstream nodes exist
UpstrNodes <- which(ginet$down == id & !is.na(ginet$down))
if(length(UpstrNodes) > 0) {
maxLength <- max(ginet$length[UpstrNodes])
params$LAG <- maxLength * 1000 / ((params$Tau0 + params$K0) * 3600)
sLag <- "LAG"
} else {
sLag <- NULL
Param <- c(
Param,
maxLength / ((nodeParam$Tau0 + nodeParam$K0) * 3600)
)
}
girop$params[girop$id == id] <- list(params[params$id == id, c("S", "IGF", "KR", "T", sLag)])
}
return(Param)
})
dplyr::select(params, id, S, IGF, KR, T, LAG)
```
## GriwmRunOptions object
......@@ -107,14 +109,14 @@ RunOptions <- CreateRunOptions(
OutputsModelsClimAware <- RunModel(
InputsModel = InputsModel,
RunOptions = RunOptions,
girop
Param = ParamClimAware
)
```
## Save data for next vignettes
```{r}
save(RunOptions, file = "_cache/V02.RData")
save(RunOptions, ParamClimAware, file = "_cache/V02.RData")
```
## Plot the result for each basin
......
Supports Markdown
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