diff --git a/NAMESPACE b/NAMESPACE index 8f49d110b142b900f3a263957c2cb7f8984adc76..f5d9c3eea5d6cdbc33709ad4101b7832c18cadd5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,7 +17,7 @@ export(CreateCalibOptions) export(CreateInputsCrit) export(CreateInputsModel) export(CreateRunOptions) -export(Ginet) +export(Griwrm) export(RunModel) export(getNodeRanking) import(airGR) diff --git a/R/CreateInputsModel.Griwrm.R b/R/CreateInputsModel.Griwrm.R index 82c3f97ad2fc071fe4107cdb4f3218f87c4d753e..8395318780f5da1737de66467bda4bb8a9145d6b 100644 --- a/R/CreateInputsModel.Griwrm.R +++ b/R/CreateInputsModel.Griwrm.R @@ -1,6 +1,6 @@ #' Create InputsModel object for a GRIWRM network #' -#' @param x Ginet object describing the diagram of the semi-distributed model, see \code{[Ginet]}. +#' @param x Griwrm object describing the diagram of the semi-distributed model, see \code{[Griwrm]}. #' @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. @@ -38,19 +38,19 @@ CreateEmptyGriwrmInputsModel <- function() { #' Create one InputsModel for a GRIWRM node #' #' @param id string of the node identifier -#' @param ginet See \code{[Ginet]}. +#' @param griwrm See \code{[Griwrm]}. #' @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, DatesR, Precip, PotEvap, Qobs) { - node <- ginet[ginet$id == id,] - FUN_MOD <- ginet$model[ginet$id == id] +CreateOneGriwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs) { + node <- griwrm[griwrm$id == id,] + FUN_MOD <- griwrm$model[griwrm$id == id] # Set hydraulic parameters - UpstreamNodes <- ginet$id[ginet$down == id & !is.na(ginet$down)] + UpstreamNodes <- griwrm$id[griwrm$down == id & !is.na(griwrm$down)] Qupstream <- NULL LengthHydro <- NULL BasinAreas <- NULL @@ -65,10 +65,10 @@ CreateOneGriwrmInputsModel <- function(id, ginet, DatesR, Precip, PotEvap, Qobs) Qupstream <- cbind(Qupstream, Qupstream1) } } - LengthHydro <- ginet$length[ginet$id %in% UpstreamNodes] + LengthHydro <- griwrm$length[griwrm$id %in% UpstreamNodes] BasinAreas <- c( - ginet$area[ginet$id %in% UpstreamNodes], - node$area - sum(ginet$area[ginet$id %in% UpstreamNodes]) + griwrm$area[griwrm$id %in% UpstreamNodes], + node$area - sum(griwrm$area[griwrm$id %in% UpstreamNodes]) ) } diff --git a/R/Griwrm.R b/R/Griwrm.R index d3ba4bc35e089fd0f7a634ae27e498b7e5d28c04..28665823ee3adacefdc62f73a755d20a029eed49 100644 --- a/R/Griwrm.R +++ b/R/Griwrm.R @@ -5,35 +5,35 @@ #' @param cols named list or vector for matching columns of `db` parameter. By default, mandatory columns names are: `id`, `down`, `length`. 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 `Ginet` class object containing the description of diagram of the semi-distributed catchment model +#' @return `Griwrm` 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", model = "model"), keep_all = FALSE) { +Griwrm <- 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) { db <- dplyr::select(db, names(cols)) } - class(db) <- append(class(db), c("Griwrm", "Ginet")) + class(db) <- append(class(db), c("Griwrm", "Griwrm")) db } #' Sort the nodes from upstream to downstream. #' -#' @param ginet See \code{[Ginet]}. +#' @param griwrm See \code{[Griwrm]}. #' #' @return vector with the ordered node names. #' @export -getNodeRanking <- function(ginet) { - if(!is(ginet, "Ginet")) { - stop("getNodeRanking: ginet argument should be of class Ginet") +getNodeRanking <- function(griwrm) { + if(!is(griwrm, "Griwrm")) { + stop("getNodeRanking: griwrm argument should be of class Griwrm") } # Rank 1 - rank <- setdiff(ginet$id, ginet$down) + rank <- setdiff(griwrm$id, griwrm$down) ranking <- rank # Next ranks - while(any(ginet$id %in% rank)) { - rank <- ginet$down[ginet$id %in% rank] + while(any(griwrm$id %in% rank)) { + rank <- griwrm$down[griwrm$id %in% rank] ranking <- c(ranking, rank) } ranking <- unique(ranking, fromLast = TRUE) diff --git a/vignettes/V01_First_network.Rmd b/vignettes/V01_First_network.Rmd index a9833777ac43dd086c6971e43d4ef62a6fa802c0..a14ce2fb5512e0ddfb666aa0889dc17b929dbba1 100644 --- a/vignettes/V01_First_network.Rmd +++ b/vignettes/V01_First_network.Rmd @@ -22,13 +22,13 @@ List of nodes ```{r} seine_nodes <- readr::read_delim( - file = system.file("seine_data", "network_gauging_stations.txt", package = "griwrm"), + file = system.file("seine_data", "network_gauging_stations.txt", package = "griwrm"), delim = "\t" ) 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: +Create the griwrm object which lists the nodes and describes the network diagram. It's a dataframe of class `Griwrm` and `Griwrm` with specific column names: - `id`: the identifier of the node in the network. - `down`: the identifier of the next hydrological node downstream. @@ -36,25 +36,25 @@ Create the ginet object which lists the nodes and describes the network diagram. - `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. +`Griwrm` function helps to rename the columns of the dataframe and assign the variable classes. ```{r} seine_nodes$model <- "RunModel_GR4J" -# Generate the ginet object -ginet <- Ginet(seine_nodes, list(id = "id_sgl", down = "id_aval", length = "distance_aval")) -ginet +# Generate the griwrm object +griwrm <- Griwrm(seine_nodes, list(id = "id_sgl", down = "id_aval", length = "distance_aval")) +griwrm ``` -## Observation time series +## Observation time series -Loading hydrometeorological data on the Seine river basin from the ClimAware project: +Loading hydrometeorological data on the Seine river basin from the ClimAware project: ```{r, warning=FALSE, message=FALSE} Precip <- NULL PotEvap <- NULL -Qobs <- NULL +Qobs <- NULL MergeTS <- function(dfOld, id, dfNew) { names(dfNew) <- c("DatesR", id) # Renaming columns of the new input into date and sub-basin ID @@ -66,13 +66,13 @@ MergeTS <- function(dfOld, id, dfNew) { return(dfOut) } -for(id in ginet$id) { - url <- +for(id in griwrm$id) { + url <- file.path( - "https://stratus.irstea.fr/d/0b18e688851a45478f7a/files/?p=/climaware_hydro/Q_OBS_NAT", + "https://stratus.irstea.fr/d/0b18e688851a45478f7a/files/?p=/climaware_hydro/Q_OBS_NAT", paste0(id, "_NAT.txt&dl=1") ) - ts <- readr::read_delim(file = url, + ts <- readr::read_delim(file = url, delim = ";", skip = 16, trim_ws = TRUE) # Date conversion to POSIX ts$Date <- as.POSIXlt(lubridate::ymd(ts$Date)) @@ -81,7 +81,7 @@ for(id in ginet$id) { # ETP column is merged into PotEvap dataframe PotEvap <- MergeTS(PotEvap, id, ts[,c("Date", "ETP")]) # Convert Qobs from m3/s to mm/time step - ts$Qnat <- ts$Qnat * 86.4 / ginet$area[ginet$id == id] + ts$Qnat <- ts$Qnat * 86.4 / griwrm$area[griwrm$id == id] # Setting data gaps to NA ts$Qnat[ts$Qnat <= 0] <- NA # Qnat column is merged into Qobs dataframe @@ -97,11 +97,11 @@ Qobs$DateR <- NULL The GRIWRM InputsModel object is a list of airGR InputsModel. The identifier of the sub-basin is used as key in the list which is ordered from upstream to downstream. -The airGR CreateInputsModel function is extended in order to handle the ginet object which describe the basin diagram: +The airGR CreateInputsModel function is extended in order to handle the griwrm object which describe the basin diagram: ```{r} -InputsModel <- CreateInputsModel(ginet, DatesR, Precip, PotEvap, Qobs) +InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs) ``` @@ -111,6 +111,6 @@ InputsModel <- CreateInputsModel(ginet, DatesR, Precip, PotEvap, Qobs) ```{r} dir.create("_cache", showWarnings = FALSE) -save(ginet, Qobs, InputsModel, file = "_cache/V01.RData") +save(griwrm, Qobs, InputsModel, file = "_cache/V01.RData") ``` diff --git a/vignettes/V03_First_Calibration.Rmd b/vignettes/V03_First_Calibration.Rmd index fd3aba36c5fb22bed83e3c988d19ccbc06325f67..276f432e859c52327943f163bbbadf5e367c3093 100644 --- a/vignettes/V03_First_Calibration.Rmd +++ b/vignettes/V03_First_Calibration.Rmd @@ -25,7 +25,7 @@ load("_cache/V02.RData") ```{r} InputsCrit <- CreateInputsCrit( - InputsModel = InputsModel, + InputsModel = InputsModel, FUN_CRIT = airGR::ErrorCrit_KGE2, RunOptions = RunOptions, Qobs = Qobs ) @@ -53,7 +53,7 @@ save(OutputsCalib, file = "_cache/V03.RData") ## Run model with Michel calibration ```{r} -ParamMichel <- sapply(ginet$id, function(x) {OutputsCalib[[x]]$Param}) +ParamMichel <- sapply(griwrm$id, function(x) {OutputsCalib[[x]]$Param}) OutputsModels <- RunModel( InputsModel = InputsModel, @@ -73,7 +73,7 @@ save(ParamMichel, file = "_cache/V03.RData") ```{r, fig.height = 5, fig.width = 8} htmltools::tagList(lapply( - names(OutputsModels), + names(OutputsModels), function(x) { plot(OutputsModels[[x]], Qobs = Qobs[RunOptions[[x]]$IndPeriod_Run,x] , main = x) } diff --git a/vignettes/v02_First_run.Rmd b/vignettes/v02_First_run.Rmd index 0146449e318646b9b283c7d5e1f711f0226fb9b2..91712d0e7277e187215a32879f849446676c7a20 100644 --- a/vignettes/v02_First_run.Rmd +++ b/vignettes/v02_First_run.Rmd @@ -28,7 +28,7 @@ Run `vignette("01_First_network", package = "griwrm")` before this one in order load("_cache/V01.RData") ``` -### Loading +### Loading Data comes from calibration of ClimAware project with naturalised flows. @@ -40,25 +40,25 @@ ClimAwareParams <- readr::read_csv( ClimAwareParams ``` -The lag $\tau_0$ and route $K_0$ parameters of TGR are expressed as time delay in hours corresponding to the delay time between the farest upstream inlet and the outlet of the sub-basin. +The lag $\tau_0$ and route $K_0$ parameters of TGR are expressed as time delay in hours corresponding to the delay time between the farest upstream inlet and the outlet of the sub-basin. Almost all sub basin has only a lag parameter. The only exception is for La Marne à Noisiel (NOISI_17) that has a routing parameter which can be approximated to a single lag parameter equals to $\tau_0 + K_0$. This lag parameter has to be converted in a speed in m/s used in the airGR lag model: ```{r} # Convert TGR routing parameter into speed -params <- merge(ginet, ClimAwareParams, by.x = "id", by.y = "id_sgl") +params <- merge(griwrm, ClimAwareParams, by.x = "id", by.y = "id_sgl") -ParamClimAware <- sapply(ginet$id, function(id) { +ParamClimAware <- sapply(griwrm$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)) + UpstrNodes <- which(griwrm$down == id & !is.na(griwrm$down)) if(length(UpstrNodes) > 0) { - maxLength <- max(ginet$length[UpstrNodes]) + maxLength <- max(griwrm$length[UpstrNodes]) Param <- c( - Param, + Param, maxLength / ((nodeParam$Tau0 + nodeParam$K0) * 3600) ) } @@ -82,7 +82,7 @@ library(lubridate) IndPeriod_Run <- seq( which(InputsModel[[1]]$DatesR == (InputsModel[[1]]$DatesR[1] + months(12))), # Set aside warm-up period length(InputsModel[[1]]$DatesR) # Until the end of the time series -) +) ``` The warmup period could also be defined as is: @@ -95,8 +95,8 @@ IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1) ```{r} RunOptions <- CreateRunOptions( - InputsModel = InputsModel, - IndPeriod_WarmUp = IndPeriod_WarmUp, + InputsModel = InputsModel, + IndPeriod_WarmUp = IndPeriod_WarmUp, IndPeriod_Run = IndPeriod_Run ) ``` @@ -123,7 +123,7 @@ save(RunOptions, ParamClimAware, file = "_cache/V02.RData") ```{r, fig.height = 5, fig.width = 8} htmltools::tagList(lapply( - names(OutputsModelsClimAware), + names(OutputsModelsClimAware), function(x) { plot(OutputsModelsClimAware[[x]], Qobs = Qobs[IndPeriod_Run,x] , main = x) }