Commit 4d64bb1c authored by Dorchies David's avatar Dorchies David
Browse files

feat: Remove Gits class object

- remove Gits dependencies
- update vignettes
- correct management of data gaps in observed flow (previously interpolated and now removed from the calibration and graphics)

Refs #7
Showing with 98 additions and 96 deletions
+98 -96
Package: griwrm Package: griwrm
Title: GR based Integrated Water Resource Management Title: GR based Integrated Water Resource Management
Version: 0.2.0 Version: 0.2.1
Authors@R: Authors@R:
person(given = "David", person(given = "David",
family = "Dorchies", family = "Dorchies",
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
#' @param InputsModel object of class \emph{GriwrmInputsModel}, see \code{\link{CreateInputsModel.Griwrm}} for details. #' @param InputsModel object of class \emph{GriwrmInputsModel}, see \code{\link{CreateInputsModel.Griwrm}} for details.
#' @param FUN_CRIT \[function (atomic or list)\] error criterion function (e.g. \code{\link[airGR]{ErrorCrit_RMSE}}, \code{\link[airGR]{ErrorCrit_NSE}}) #' @param FUN_CRIT \[function (atomic or list)\] error criterion function (e.g. \code{\link[airGR]{ErrorCrit_RMSE}}, \code{\link[airGR]{ErrorCrit_NSE}})
#' @param RunOptions object of class \emph{GriwrmRunOptions}, see \code{[CreateRunOptions.Griwrm]} for details. #' @param RunOptions object of class \emph{GriwrmRunOptions}, see \code{[CreateRunOptions.Griwrm]} for details.
#' @param gits object of class \emph{Gits}, see [Gits]. #' @param Qobs matrix or data frame containing observed flows. Column names correspond to nodes ID
#' @param ... further arguments passed to \code{\link[airGR]{CreateInputsCrit}}. #' @param ... further arguments passed to \code{\link[airGR]{CreateInputsCrit}}.
#' #'
#' @return Object of class \emph{GriwrmInputsCrit} #' @return Object of class \emph{GriwrmInputsCrit}
...@@ -10,7 +10,7 @@ ...@@ -10,7 +10,7 @@
CreateInputsCrit.GriwrmInputsModel <- function(InputsModel, CreateInputsCrit.GriwrmInputsModel <- function(InputsModel,
FUN_CRIT = airGR::ErrorCrit_NSE, FUN_CRIT = airGR::ErrorCrit_NSE,
RunOptions, RunOptions,
gits, Qobs,
...) { ...) {
InputsCrit <- list() InputsCrit <- list()
class(InputsCrit) <- append(class(InputsCrit), "GriwrmInputsCrit") class(InputsCrit) <- append(class(InputsCrit), "GriwrmInputsCrit")
...@@ -20,7 +20,7 @@ CreateInputsCrit.GriwrmInputsModel <- function(InputsModel, ...@@ -20,7 +20,7 @@ CreateInputsCrit.GriwrmInputsModel <- function(InputsModel,
InputsModel = IM, InputsModel = IM,
FUN_CRIT = FUN_CRIT, FUN_CRIT = FUN_CRIT,
RunOptions = RunOptions[[IM$id]], RunOptions = RunOptions[[IM$id]],
Obs = gits[[IM$id]]$Qobs[RunOptions[[IM$id]]$IndPeriod_Run], Obs = Qobs[RunOptions[[IM$id]]$IndPeriod_Run, IM$id],
... ...
) )
} }
......
...@@ -2,19 +2,25 @@ ...@@ -2,19 +2,25 @@
#' #'
#' @param x Ginet object describing the diagram of the semi-distributed model, see \code{[Ginet]}. #' @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 girop Girop object giving the run-off model parameters, see \code{[Girop]}.
#' @param gits Gits object giving the observation time series, see \code{[Gits]}. #' @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.
#' @param Qobs Matrix or data frame of numeric containing potential observed flow in mm. Column names correspond to node IDs.
#' @param verbose (optional) boolean indicating if the function is run in verbose mode or not, default = \code{TRUE} #' @param verbose (optional) boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}
#' @param ... further arguments passed to \code{\link[airGR]{CreateInputsModel}}. #' @param ... further arguments passed to \code{\link[airGR]{CreateInputsModel}}.
#' #'
#' @return GriwrmInputsModel object equivalent to airGR InputsModel object for a semi-distributed model (See \code{\link[airGR]{CreateInputsModel}}) #' @return GriwrmInputsModel object equivalent to airGR InputsModel object for a semi-distributed model (See \code{\link[airGR]{CreateInputsModel}})
#' @export #' @export
CreateInputsModel.Griwrm <- function(x, girop, gits, verbose = TRUE, ...) { CreateInputsModel.Griwrm <- function(x, girop, DatesR, Precip, PotEvap, Qobs, verbose = TRUE, ...) {
InputsModel <- CreateEmptyGriwrmInputsModel() InputsModel <- CreateEmptyGriwrmInputsModel()
Qobs[is.na(Qobs)] <- -99 # airGRCreateInputsModel doesn't accept NA values
for(id in getNodeRanking(x)) { for(id in getNodeRanking(x)) {
if(verbose) cat("CreateInputsModel.griwrm: Treating sub-basin", id, "...\n") if(verbose) cat("CreateInputsModel.griwrm: Treating sub-basin", id, "...\n")
InputsModel[[id]] <- CreateOneGriwrmInputsModel(id, x, girop, gits, ...) InputsModel[[id]] <- CreateOneGriwrmInputsModel(
id, x, girop, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
)
} }
return(InputsModel) return(InputsModel)
} }
...@@ -35,10 +41,13 @@ CreateEmptyGriwrmInputsModel <- function() { ...@@ -35,10 +41,13 @@ CreateEmptyGriwrmInputsModel <- function() {
#' @param id string of the node identifier #' @param id string of the node identifier
#' @param ginet See \code{[Ginet]}. #' @param ginet See \code{[Ginet]}.
#' @param girop See \code{[Girop]}. #' @param girop See \code{[Girop]}.
#' @param gits See \code{[Gits]}. #' @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. #' @return \emph{InputsModel} object for one.
CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) { CreateOneGriwrmInputsModel <- function(id, ginet, girop, DatesR, Precip, PotEvap, Qobs) {
node <- ginet[ginet$id == id,] node <- ginet[ginet$id == id,]
FUN_MOD <- girop$model[girop$id == id] FUN_MOD <- girop$model[girop$id == id]
...@@ -51,7 +60,7 @@ CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) { ...@@ -51,7 +60,7 @@ CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) {
if(length(UpstreamNodes) > 0) { if(length(UpstreamNodes) > 0) {
# Sub-basin with hydraulic routing # Sub-basin with hydraulic routing
for(idUpstrNode in UpstreamNodes) { for(idUpstrNode in UpstreamNodes) {
Qupstream1 <- matrix(gits[[idUpstrNode]]$Qobs, ncol = 1) Qupstream1 <- matrix(Qobs[,idUpstrNode], ncol = 1)
if(is.null(Qupstream)) { if(is.null(Qupstream)) {
Qupstream <- Qupstream1 Qupstream <- Qupstream1
} else { } else {
...@@ -68,9 +77,9 @@ CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) { ...@@ -68,9 +77,9 @@ CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) {
# Set model inputs with the airGR function # Set model inputs with the airGR function
InputsModel <- CreateInputsModel( InputsModel <- CreateInputsModel(
FUN_MOD, FUN_MOD,
DatesR = gits$date, DatesR = DatesR,
Precip = gits[[id]]$Precip, Precip = Precip,
PotEvap = gits[[id]]$PotEvap, PotEvap = PotEvap,
Qupstream = Qupstream, Qupstream = Qupstream,
LengthHydro = LengthHydro, LengthHydro = LengthHydro,
BasinAreas = BasinAreas BasinAreas = BasinAreas
......
#' Title
#'
#' @param id string of the id of the node
#' @param ts numeric matrix or data frame containing 3 columns for precipitation, evaporation, and observed flow
#' @param cols named list or vector used for matching the columns of ts with the required columns names which are "Precip", "PotEvap", and "Qobs".
#'
#' @return \emph{Gits} class object which is a list containing a `date` element (Vector of PosiXlt timestamps) and an element named the id of the node containing a dataframe with observed data.
#' @export
Gits <- function(id, ts,
cols = list(date = "date", Precip = "Precip", PotEvap = "PotEvap", Qobs = "Qobs")) {
cols <- as.list(cols)
ts <- dplyr::rename(ts, unlist(cols))
if(any(is.na(ts$Qobs))) {
stop("Qobs should not contain any NA")
}
if(any(ts$Qobs < 0)) {
stop("Qobs should be strictly positive")
}
gitsOut <- list(date = ts$date)
cols$date <- NULL
gitsOut[[id]] <- dplyr::select(ts, names(cols))
class(gitsOut) <- append(class(gitsOut), "Gits")
gitsOut
}
#' Merge two gits objects with identical date time series.
#'
#' @param x Gits object to merge (See [Gits]).
#' @param y Gits object to merge (See [Gits]).
#' @param ... For merge generic function compatibility.
#'
#' @return Gits object merged with one item `Date` and Items corresponding to each node.
#' @export
merge.Gits <- function(x, y, ...) {
if(!is(y, "Gits")) {
stop("A Gits class object can only be merged with a Gits class object")
}
if(! identical(x$date, y$date)) {
stop("Time series dates are not identical")
}
y$date <- NULL
x <- utils::modifyList(x, y)
}
id_sgl id_hydro lambert2.x lambert2.y area description id_aval distance_aval
TRANN_01 766626.06 2369152.25 1557.06 l'Aube à Trannes ARCIS_24 68100
GURGY_02 H2221010 689713.00 2320549.00 3819.77 L'Yonne à Gurgy COURL_23 83612
BRIEN_03 H2482010 695313.00 2332549.00 2979.77 L'Armançon à Brienon-sur-Armançon COURL_23 84653
STDIZ_04 H5071010 791113.00 2407349.00 2347.53 La Marne à Saint-Dizier CHALO_21 85570
PARIS_05 H5920010 602213.00 2427449.00 43824.66 La Seine à Paris
BAR-S_06 H0400010 751913.00 2348349.00 2340.37 La Seine à Bar-sur-Seine MERY-_22 79766
CHAUM_07 716518.38 2241747.00 216.50 L'Yonne à Chaumard GURGY_02 153074
CUSSY_08 H2172310 726013.00 2275149.00 247.99 Le Cousin à Cussy-les-Forges GURGY_02 91378
STGER_09 718512.88 2266649.25 402.74 La Cure à St-Germain GURGY_02 94152
GUILL_10 H2322020 731013.00 2282349.00 488.71 Le Serein à Guillon CHABL_12 66026
AISY-_11 H2452020 742413.00 2298149.00 1349.51 L'Armançon à Aisy-sur-Armançon BRIEN_03 102428
CHABL_12 H2342010 709613.00 2314149.00 1116.27 Le Serein à Chablis COURL_23 111781
NOGEN_13 685912.88 2389349.25 9182.39 La Seine à Nogent-sur-Seine MONTE_15 63215
EPISY_14 H3621010 633413.00 2371049.00 3916.71 Le Loing à Épisy ALFOR_16 89196
MONTE_15 645112.88 2376849.25 21199.39 La Seine à Montereau ALFOR_16 94475
ALFOR_16 H4340020 606013.00 2420349.00 30784.71 La Seine à Alfortville PARIS_05 9263
NOISI_17 H5841010 620913.00 2428949.00 12547.72 La Marne à Noisiel PARIS_05 39384
MONTR_18 H5752020 638013.00 2431849.00 1184.81 Le Grand Morin à Montry NOISI_17 37915
LOUVE_19 H5083050 791613.00 2393949.00 461.74 La Blaise à Louvemont CHALO_21 86165
LASSI_20 H1362010 759513.00 2385549.00 876.53 La Voire à Lassicourt ARCIS_24 43618
CHALO_21 H5201010 747713.00 2441349.00 6291.55 La Marne à Châlons-sur-Marne NOISI_17 237937
MERY-_22 H0810010 714913.00 2390949.00 3899.62 La Seine à Méry-sur-Seine NOGEN_13 49933
COURL_23 H2721010 660813.00 2370449.00 10687.35 L'Yonne à Courlon-sur-Yonne MONTE_15 26159
ARCIS_24 H1501010 733313.00 2394749.00 3594.60 L'Aube à Arcis-sur-Aube NOGEN_13 70926
VITRY_25 H5172010 768513.00 2418849.00 2109.14 La Saulx à Vitry-en-Perthois CHALO_21 38047
...@@ -22,8 +22,8 @@ List of nodes ...@@ -22,8 +22,8 @@ List of nodes
```{r} ```{r}
seine_nodes <- readr::read_delim( seine_nodes <- readr::read_delim(
file = "https://stratus.irstea.fr/d/0b18e688851a45478f7a/files/?p=/climaware_hydro/griwrm_network.txt&dl=1", file = system.file("seine_data", "network_gauging_stations.txt", package = "griwrm"),
delim = ";" delim = "\t"
) )
seine_nodes seine_nodes
``` ```
...@@ -77,40 +77,45 @@ urls <- ...@@ -77,40 +77,45 @@ urls <-
) )
names(urls) <- ginet$id names(urls) <- ginet$id
load_ts <- function(x) { Precip <- NULL
ts <- readr::read_delim(file = x, PotEvap <- NULL
delim = ";", skip = 16, trim_ws = TRUE) Qobs <- NULL
ts$Date <- as.POSIXlt(lubridate::ymd(ts$Date))
# Interpolation of data gap in the discharge time serie MergeTS <- function(dfOld, id, dfNew) {
ts[ts$Qnat < 0, "Qnat"] <- NA names(dfNew) <- c("DatesR", id) # Renaming columns of the new input into date and sub-basin ID
if(is.na(ts$Qnat[nrow(ts)])) { if(is.null(dfOld)) {
ts$Qnat[nrow(ts)] <- 0 # No value at the end: time serie converge to zero dfOut <- dfNew # Creation of the first column
} else {
dfOut <- merge(dfOld, dfNew, by = "DatesR", all = TRUE) # Merge the new columns respecting to date column
} }
ts$Qnat <- zoo::na.approx(ts$Qnat) return(dfOut)
ts
} }
l <- lapply(urls, load_ts)
```
`Gits` object is a list containing a item named `date` with a timestamp vector of the time series and items named by the identifier of each node. These items contain a dataframe with the observations.
The Gits function creates a `Gits` object
```{r}
gits <- Gits(ginet$id[1], l[[ginet$id[1]]], cols = list(date = "Date", Precip = "Ptot", PotEvap = "ETP", Qobs = "Qnat"))
```
Copy the observations for each node the ginet network:
```{r}
for(id in ginet$id) { for(id in ginet$id) {
l[[id]]$Qnat url <-
l[[id]]$Qnat <- l[[id]]$Qnat * 86.4 / girop$area[girop$id == id] file.path(
l[[id]]$Qnat[l[[id]]$Qnat < 0] <- NA "https://stratus.irstea.fr/d/0b18e688851a45478f7a/files/?p=/climaware_hydro/Q_OBS_NAT",
gits <- merge(gits, Gits(id, l[[id]], cols = list(date = "Date", Precip = "Ptot", PotEvap = "ETP", Qobs = "Qnat"))) paste0(id, "_NAT.txt&dl=1")
)
ts <- readr::read_delim(file = url,
delim = ";", skip = 16, trim_ws = TRUE)
# Date conversion to POSIX
ts$Date <- as.POSIXlt(lubridate::ymd(ts$Date))
# Ptot column is merged into Precip dataframe
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]
# Setting data gaps to NA
ts$Qnat[ts$Qnat <= 0] <- NA
# Qnat column is merged into Qobs dataframe
Qobs <- MergeTS(Qobs, id, ts[,c("Date", "Qnat")])
} }
DatesR <- Precip$DatesR
Precip$DateR <- NULL
PotEvap$DateR <- NULL
Qobs$DateR <- NULL
``` ```
## Generate the GRIWRM InputsModel object ## Generate the GRIWRM InputsModel object
...@@ -121,7 +126,7 @@ The airGR CreateInputsModel function is extended in order to handle the ginet ob ...@@ -121,7 +126,7 @@ The airGR CreateInputsModel function is extended in order to handle the ginet ob
```{r} ```{r}
InputsModel <- CreateInputsModel(ginet, girop, gits) InputsModel <- CreateInputsModel(ginet, girop, DatesR, Precip, PotEvap, Qobs)
``` ```
...@@ -129,7 +134,8 @@ InputsModel <- CreateInputsModel(ginet, girop, gits) ...@@ -129,7 +134,8 @@ InputsModel <- CreateInputsModel(ginet, girop, gits)
## Save data for next vignettes ## Save data for next vignettes
```{r} ```{r}
dir.create("_cache", showWarnings = FALSE) dir.create("_cache", showWarnings = FALSE)
save(ginet, girop, gits, InputsModel, file = "_cache/V01.RData") save(ginet, girop, Qobs, InputsModel, file = "_cache/V01.RData")
``` ```
...@@ -27,7 +27,7 @@ load("_cache/V02.RData") ...@@ -27,7 +27,7 @@ load("_cache/V02.RData")
InputsCrit <- CreateInputsCrit( InputsCrit <- CreateInputsCrit(
InputsModel = InputsModel, InputsModel = InputsModel,
FUN_CRIT = airGR::ErrorCrit_KGE2, FUN_CRIT = airGR::ErrorCrit_KGE2,
RunOptions = RunOptions, gits = gits RunOptions = RunOptions, Qobs = Qobs
) )
str(InputsCrit) str(InputsCrit)
``` ```
...@@ -65,13 +65,20 @@ OutputsModels <- RunModel( ...@@ -65,13 +65,20 @@ OutputsModels <- RunModel(
) )
``` ```
## Save calibration data for next vignettes
```{r}
save(giropMichel, file = "_cache/V03.RData")
```
## Plot the result for each basin ## Plot the result for each basin
```{r, fig.height = 5, fig.width = 8} ```{r, fig.height = 5, fig.width = 8}
htmltools::tagList(lapply( htmltools::tagList(lapply(
names(OutputsModels), names(OutputsModels),
function(x) { function(x) {
plot(OutputsModels[[x]], Qobs = gits[[x]]$Qobs[RunOptions[[x]]$IndPeriod_Run] , main = x) plot(OutputsModels[[x]], Qobs = Qobs[RunOptions[[x]]$IndPeriod_Run,x] , main = x)
} }
)) ))
``` ```
......
...@@ -123,7 +123,7 @@ save(RunOptions, file = "_cache/V02.RData") ...@@ -123,7 +123,7 @@ save(RunOptions, file = "_cache/V02.RData")
htmltools::tagList(lapply( htmltools::tagList(lapply(
names(OutputsModelsClimAware), names(OutputsModelsClimAware),
function(x) { function(x) {
plot(OutputsModelsClimAware[[x]], Qobs = gits[[x]]$Qobs[IndPeriod_Run] , main = x) plot(OutputsModelsClimAware[[x]], Qobs = Qobs[IndPeriod_Run,x] , main = x)
} }
)) ))
``` ```
......
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