Skip to content

GitLab

Commit 54fb4300 authored by Dorchies David's avatar Dorchies David
Browse files

feat: migrate vgest related functions to rvgest package 

Refs rvgest#1
parent 24721744
Hide whitespace changes
Inline Side-by-side
Showing with 8 additions and 917 deletions
R/get_objectives.R deleted 100644 → 0
View file @ 24721744
#' Load objective data for \code{\link{vgest_run}}
#'
#' @param objective_thresholds File location of threshold table
#' @param objective_stations File location of objective station table
#' @param lakes File location of lake table
#'
#' @return A dataframe containing one line by objective and the following columns:
#'
#' - station (character): identifier of the objective station
#' - flood (boolean): TRUE for high flow mitigation objective, and FALSE for low flow support
#' - level (character): code composed with "l" for low-flow and "h" for high flow followed by the number of the level
#' - threshold (numeric): value of the threshold in m3/s
#' - lakes (dataframe): Dataframe containing lake details (name, min storage, max storage)
#'
#' @export
#'
#' @examples
#' # Get objectives stored in IRMaRA package
#' df <- get_objectives()
#' # Objectives at Paris
#' df[df[,"station"] == "PARIS_05",]
#' # Lake details concerning the 40th objective in the table
#' df[57, "lakes"]
get_objectives <- function(
objective_thresholds = app_sys("seine", "objective_thresholds.txt"),
objective_stations = app_sys("seine", "objective_stations.json"),
lakes = app_sys("seine", "lakes.txt")
) {
thresholds <- read.delim(objective_thresholds)
stations <- jsonlite::fromJSON(
readLines(objective_stations)
)
dfLakes <- read.delim(lakes)
row.names(dfLakes) <- dfLakes$name
bFirst = TRUE
for(iStation in 1:nrow(thresholds)) {
station <- thresholds[iStation, 1]
for(type in c("l", "h")) {
bFlood = (type == "h")
for(level in grep(paste0(type, ".*"), names(thresholds))) {
threshold <- thresholds[iStation, level]
df1 <- data.frame(
station = station,
flood = bFlood,
level = names(thresholds)[level],
threshold = threshold
)
df1$lakes <- list(dfLakes[stations[[station]]$lakes,])
if(bFirst) {
df <- df1
bFirst <- FALSE
} else {
df <- rbind(df, df1)
}
}
}
}
class(df) <- c("Objectives", class(df))
return(df)
}
\ No newline at end of file
R/get_vobj_ts.R deleted 100644 → 0
View file @ 24721744
#' Method to transform VOBJ data into time series
#'
#' Transform VOBJ data read by [vgest_read], [vgest_read_one] or [vgest_read_all].
#'
#' @param x
#'
#' @return
#' @export
#'
#' @examples
#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database"))
get_vobj_ts <- function(x) {
UseMethod("get_vobj_ts", x)
}
#' Transform a VOBJ data into time series
#'
#' @param dfVobj a dataframe provided by [vgest_read_one] or [vgest_read], or an item of the list provided by [vgest_read_all]
#'
#' @return A dataframe with one column `date` and one column `storage` containing the time series of the lake storage in m3.
#' @export
#'
#' @examples
#' # Get the first lake of the first objective
#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database")[[1]])
#' # Can also be done by
#' dfTs <- get_vobj_ts(vgest_read_one(1, get_objectives()[1,], "./database"))
get_vobj_ts.Vobj <- function(dfVobj) {
# Build the date vector
# Get year list
vYears <- as.numeric(
gsub("Y", "", names(dfVobj)[grep(pattern = "Y[0-9]{4}", names(dfVobj))])
)
# Combine dates dd/mm with years
vDates <- as.vector(sapply(vYears, function(x) {paste(dfVobj[,"DAY"], x, sep = "/")}))
# Convert dates to POSIX
vDates <- lubridate::dmy(vDates)
ts <- unlist(dfVobj[,grep(pattern = "Y[0-9]{4}", names(dfVobj))])
data.frame(date=vDates, storage=ts)
}
#' Transform list of VOBJ data into time series
#'
#' Use for binding results of several lakes for one objective at one station.
#'
#' @param lResultVobj list of VOBJ provided by [vgest_read_all]
#'
#' @return Dataframe with the time series with a column `date` and one column by lake with the storage of the lake in m3.
#' @export
#'
#' @examples
#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database"))
get_vobj_ts.ListVobj <-function(lResultVobj) {
lVobjTs <- lapply(lResultVobj, get_vobj_ts.Vobj)
df <- lVobjTs[[1]]
if(length(lVobjTs) > 1) {
for(i in 2:length(lVobjTs)) {
df <- cbind(df, lVobjTs[[i]]$storage)
}
}
names(df) <- c("date", names(lResultVobj))
return(df)
}
\ No newline at end of file
R/plot_isofrequency.R deleted 100644 → 0
View file @ 24721744
#' Plot one lake isofrequency annual curves for one objective at one station
#'
#' @param x [list] containing one [data.frame] per lake
#' @param freq [numeric] vector of frequencies to plot
#' @param result.dir [character] The path to the result database (default: "database")
#'
#' @return [NULL]
#' @export
#'
plot_isofrequency <- function(x, freq, result.dir = "database") {
lObj <- vgest_read_all(x, result.dir)
nLakes <- nrow(x$lakes[[1]])
if(nLakes > 1) {
par(mfrow=c(ceiling(nLakes/2),2))
}
lapply(1:nLakes, function(i) {plot_isofrequency_lake(lObj[[i]], freq, x$lakes[[1]][i,], ceiling(i/2))})
sLowHigh <- c("low", "high")
mtext(
paste("Objective", x$threshold, "m3/s", sLowHigh[x$flood + 1], "flows threshold at", x$station),
side = 3, col = "black", line = -2, cex = 1.2, font = 2, outer = TRUE
)
}
#' Plot one lake isofrequency annual curves for one objective at one station
#'
#' This function is mainly called by [plot_isofrequency()].
#'
#' @param vObj dataframe produced by [vgest_read_one()] and stored in a list by [vgest_read_all()]
#' @param frequencies vector of frequencies to plot
#' @param lake lake data extract from column `lakes` of objective data given by [get_objectives()]
#' @param top.margin top margin applied on the plot for the title
#'
#' @return [NULL]
#' @export
#'
plot_isofrequency_lake <- function(vObj, frequencies, lake, top.margin) {
frequencies <- paste0("F", format(frequencies, digits = 5, nsmall = 5))
vObj[,-1] <- (vObj[,-1]) / 1E6 + lake$min
# Contrasted color palette adapted from https://sashat.me/2017/01/11/list-of-20-simple-distinct-colors/
sColors <- c("#003A80", "#f58231", "#3cb44b", "#e6194B", "#911eb4", "#ffe119", "#f032e6", "#9A6324", "#000075", "#808000", "#42d4f4", "#a9a9a9", "#bfef45", "#469990")
par(mgp=c(2,0.5,0), mar = c(2.5,3.5,5.5 - 2.2*top.margin,0.5))
plot(
rep(lake$min, 365),
type = "l", lwd = 1.5, col = "grey", lty = 2,
xaxt="n",
main = NULL,
ylab = "Reservoir Storage (mcm)",
xlab = NULL,
ylim = c(0, lake$max),
cex.lab = 0.6, cex.axis = 0.6
)
lines(rep(lake$max, 365), lwd = 1.5, col = "grey", lty = 2)
i <- 0
for(frequency in frequencies) {
i <- i + 1
lines(vObj[,frequency], lwd = 2, col = sColors[i])
}
mtext(
paste(lake$name, "reservoir (Active cap.", lake$max - lake$min, "mcm)"),
side = 3, font = 2, cex = 0.75
)
xtick <- c(0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
xtick <- sapply(1:13, function(x) {sum(xtick[1:x])})
axis(side=1, at=xtick, labels = unlist(strsplit("JFMAMJJASONDJ", NULL)), cex.axis = 0.6)
}
R/save_isofrequency.R deleted 100644 → 0
View file @ 24721744
save_isofrequency <- function(x, freq, result.dir = "database", output.dir = "output") {
dir.create(output.dir, showWarnings = FALSE)
file <- file.path(
output.dir, paste0(
paste("isofreq", x$station, x$level, x$threshold, sep = "_"),
".png")
)
cat("Plotting", file, "...")
png(file = file,
width = 16, height = 10, units = "cm", res = 200
)
plot_isofrequency(x, freq, result.dir)
dev.off()
cat(" OK\n")
}
\ No newline at end of file
R/vgest_add_leap.R deleted 100644 → 0
View file @ 24721744
#' Add missing leap year data to naturalised flow file
#'
#' Update the file with missing leap year data
#'
#' @param Qfile Name of the file
#' @param vgest_location location of VGEST
#'
#' @return
#' @export
#'
vgest_add_leap <- function(Qfile, vgest_location = "../vgest") {
file <- file.path(vgest_location, "DONNEES", Qfile)
Qnat <- read.delim(file)
Qnat$Dates <- as.Date(as.character(Qnat$Dates), format = "%Y%m%d")
allDates <- seq(Qnat$Dates[1], tail(Qnat$Dates, 1), by="days")
missingDates <- allDates[!allDates %in% Qnat$Dates]
missingData <- cbind(data.frame(Dates = missingDates),
matrix(data = NA, nrow = length(missingDates), ncol = ncol(Qnat) - 1))
colnames(missingData) <- names(Qnat)
Qnat <- rbind(Qnat,
missingData)
Qnat <- dplyr::arrange(Qnat, Dates)
Qnat[,-1] <- zoo::na.approx(Qnat[,-1])
Qnat$Dates <- format(Qnat$Dates, "%Y%m%d")
columnNames <- gsub(pattern = "\\.", replacement = "-", names(Qnat))
write.table(x = Qnat, file = file, quote = FALSE,
sep = "\t", row.names = FALSE, col.names = columnNames)
}
\ No newline at end of file
R/vgest_cost.R deleted 100644 → 0
View file @ 24721744
#' Calculate the cost function for one simulation (One objective at one station)
#'
#' Uses ouput of PaChrono.txt or VObj\[1-4\].dat for the calculation.
#'
#' See [vgest_cost.ListVobj], [vgest_cost.Chrono] and [vgest_cost.ListChrono] for details.
#'
#' @param data Data source, see details
#' @param ... Other arguments passed to the methods
#'
#' @return
#' @export
#'
vgest_cost <- function(data, ...) {
UseMethod("vgest_cost", data)
}
#' Calculate the total cost for all the lakes
#'
#' For each lake, it's the mean daily storage for a low-flow support objective and the mean daily available storage capacity for a high flow mitigation objective.
#'
#' @param Vobj A [matrix] or a [data.frame] with one column by lake, in the same order as `objective$lakes`
#' @param objective one row of the dataframe given by [get_objectives()]
#'
#' @return total cost of all the lakes in m3/day
#'
vgest_cost_lakes <- function(Vobj, objective) {
lakes <- objective$lakes[[1]]
sum(
sapply(
1:nrow(lakes), function(i) {
V <- Vobj[,i]
if(objective$flood) {
V <- (lakes$max[i] - lakes$min[i]) * 1E6 - V
}
mean(V, na.rm = TRUE)
}
)
)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param data the list given by [vgest_read_all()]
#' @param objective one row of the dataframe given by [get_objectives()]
#'
#' @return the total cost for one objective at one station in m3/day
#' @export
#'
#' @examples
#' \donttest{
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lResultVobj <- vgest_read_all(objective)
#' vgest_cost(lResultVobj, objective)
#' }
vgest_cost.ListVobj <- function(data, objective) {
lVobj <- lapply(data, function(x) {unlist(x[,grep(pattern = "Y[0-9]{4}", names(x))])})
Vobj <- matrix(unlist(lVobj), ncol = length(lVobj))
vgest_cost_lakes(Vobj, objective)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param data the [data.frame] given by one item of the [list] provided by [vgest_read_chrono()]
#' @param objective one row of the [data.frame] given by [get_objectives()]
#'
#' @return the total cost for one objective at one station in m3/day
#' @export
#'
#' @examples
#' \donttest{
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lChrono <- vgest_read_chrono(objective)
#' vgest_cost(lChrono[[1]], objective)
#' }
vgest_cost.Chrono <- function(data, objective) {
vgest_cost_lakes(as.data.frame(data[,grep("Vobj\\.[1-9]\\.", names(data))]), objective)
}
#' Calculate the total cost for one objective at one station
#'
#'
#' @param data the [list] given by [vgest_read_chrono()]
#' @param objective [data.frame] given by [get_objectives()]
#'
#' @return A list with items named \[station\]_\[high/low\]_\[threshold\] containing the total cost for a list of objectives and stations in m3/day
#' @export
#'
#' @examples
#' \donttest{
#' # This should be done after the execution of vgest for the concerned objective
#' objective <- get_objectives()[1,]
#' lChronos <- vgest_read_chrono(objective, distributionType = 2)
#' vgest_cost(lChronos, objective)
#' }
vgest_cost.ListChrono <- function(data, objective) {
lCosts <- lapply(1:length(data), function(i) {vgest_cost.Chrono (data[[i]], objective[i,])})
names(lCosts) <- names(data)
return(lCosts)
}
\ No newline at end of file
R/vgest_read.R deleted 100644 → 0
View file @ 24721744
#' Read result file `VOBJi.DAT` for one lake of one objective at one station
#'
#' This function is called by [vgest_read_one()], please run this one instead.
#'
#' @param file complete path to the `VOBJi.DAT` file
#' @param bFlood boolean `TRUE` for high flow mitigation objective, `FALSE` for low flow support
#'
#' @return dataframe with the content of the `VOBJi.DAT` file
#' @export
#'
#' @examples
#' \donttest{
#' df <- vgest_read("database/PARIS_05/high_800/VOBJ1.dat")
#' }
vgest_read <- function(file, bFlood) {
# Managing headers
s <- readLines(file, n = 11)
headers <- unlist(strsplit(s[10], "\\s+"))
frequencies <- unlist(strsplit(s[11], "\\s+"))[-1]
if(bFlood) {
frequencies[-1] <- rev(frequencies[-1])
}
headers <- c(
"DAY",
paste0("Y", headers[2:(which(headers == "RANG-->")-1)]),
paste0("R", headers[which(headers == "RANG-->"):(which(headers == "RETOUR-->")-1)]),
paste0("F", frequencies)
)
# Managing column widths of fixed format
widths <- c(8, rep(16, length(headers) - 1))
# Read file in fixed format
df <- read.fwf(
file = file,
widths = widths,
header = FALSE,
skip = 11,
col.names = headers
)
# Remove redundancy columns
df$`FFREQUENCE...` <- NULL
df$`RRANG...` <- NULL
class(df) <- c("Vobj", class(df))
return(df)
}
#' Read result file `VOBJi.DAT` for one lake of one objective at one station
#'
#' This function is preferred to [vgest_read()] because it builds the path for the file to read.
#'
#' @param iLake Lake number for the current station
#' @param x one line of the dataframe produced by [get_objectives()]
#' @param result.dir path where results of VGEST runs are stored (See [vgest_run_store()])
#'
#' @return dataframe with the content of the `VOBJi.DAT` file
#' @export
#'
#' @examples
#' \donttest{
#' vgest_read_one(1, get_objectives()[1,], "./database")
#' }
vgest_read_one <- function(iLake, x, result.dir) {
sLowHigh <- c("low", "high")
file <- paste0(
file.path(
result.dir, x$station, paste0(sLowHigh[x$flood + 1], "_", x$threshold), "VOBJ"
),
iLake,
".dat"
)
vgest_read(file, x$flood)
}
#' Read all result files `VOBJi.DAT` for all the lakes of one objective at one station
#'
#' @param x one line of the dataframe produced by [get_objectives()]
#' @param result.dir path where results of VGEST runs are stored (See [vgest_run_store()])
#'
#' @return list with dataframes produced by
#' @export
#'
#' @examples
#' \donttest{
#' vgest_read_all(get_objectives()[1,], "./database")
#' }
vgest_read_all <- function(x, result.dir = "database") {
lObj <- lapply(1:nrow(x$lakes[[1]]), vgest_read_one, x, result.dir)
names(lObj) <- x$lakes[[1]]$name
class(lObj) <- c("ListVobj", class(lObj))
return(lObj)
}
\ No newline at end of file
R/vgest_read_chrono.R deleted 100644 → 0
View file @ 24721744
#' Read Chrono.txt or PaChrono.txt
#'
#' @description
#' Read the time series provided by VGEST for the forward calculation (Chrono.txt) and backward calculation (PaChrono.txt)
#'
#' @details
#' The format of the file is has follow. Headers are in line 53 followed by the complete time series which is backward in time for PaChrono.txt. The columns and respected widths are:
#' - Downstream date (10)
#' - QXsous (22)
#'
#' For each lake, we next have:
#' - Upstream date (14)
#' - 25 columns of flow and storage data (22)
#' - if `DistributionType %in% c(4,5)`, 1 extra column (22)
#' - 5 columns of codes (6)
#'
#' And then at the end:
#' - 7 columns of flow (22)
#' - 3 columns of codes (8)
#' - 1 extra column of code for `DistributionType==5` (8)
#'
#' The file is saved in RDS format for quicker reading the next time.
#'
#' @param x file to read with txt extension
#' @param ... further arguments passed to or from other methods
#'
#' @return
#' @export
#' @rdname vgest_read_chrono
#'
vgest_read_chrono <- function(x, ...) {
UseMethod("vgest_read_chrono", x)
}
#'
#' @param nLakes number of lakes in the file
#' @param distributionType Distribution type. See [vgest_write_batch] details
#' @return [data.frame] with the content of the file
#' @export
#' @rdname vgest_read_chrono
#'
vgest_read_chrono.default <- function(x, nLakes, distributionType) {
# Reading cached file if exists
rdsFile <- paste0(sub('\\..[^\\.]*$', '', x), ".rds")
if(file.exists(rdsFile)) {
df <- readRDS(rdsFile)
} else {
# Widths
width1Lake <- c(14, rep(22, 25))
if(distributionType %in% c(4,5)) { width1Lake <- c(width1Lake, 22)}
width1Lake <- c(width1Lake, rep(6,5))
widths <- c(10, 22, rep(width1Lake, nLakes), rep(22,7), rep(8,3))
if(distributionType == 5) { widths <- c(widths, 8) }
# Headers
if(distributionType %in% c(4,5)) {
headLine <- 54
} else {
headLine <- 53
}
s <- readLines(x, n = headLine)[headLine]
headers = trimws(unlist(read.fwf(textConnection(s), widths = widths, as.is = TRUE), use.names = FALSE))
# Read file in fixed format
df <- read.fwf(
file = x,
col.names = headers,
widths = widths,
skip = headLine,
strip.white = TRUE
)
df[,grep(pattern = "date", names(df))] <- lapply(df[,grep(pattern = "date", names(df))], as.Date, format = "%d/%m/%Y")
df <- dplyr::arrange(df, date.aval)
saveRDS(df, file = rdsFile)
}
class(df) <- c("Chrono", class(df))
return(df)
}
#'
#' @param result.dir path for storing the result of vgest run. The result is stored in a subfolder named high or low (depending on \code{bFlood}) followed by the threshold
#' @param distributionType Distribution type. See [vgest_write_batch] details
#' @param backward boolean `TRUE` for reading "PaChrono.txt", `FALSE` for reading "Chrono.txt"
#'
#' @return A list with items named \[station\]_\[high/low\]_\[threshold\] containing [data.frame] with the content of each file
#' @export
#' @rdname vgest_read_chrono
#'
#' @examples
#' \donttest{
#' objective <- get_objectives()[1,]
#' distributionType <- 2
#' vgest_run_store(objective,
#' 1, 1, "Q_NAT_1900-2009.txt",
#' "01/01/1900", "31/12/2009", distributionType)
#' lChrono <- vgest_read_chrono(objective, distributionType)
#' }
vgest_read_chrono.Objectives <- function(x, distributionType, result.dir = "database", backward = TRUE) {
if(backward) {
fileName <- "PaChrono.txt"
} else {
fileName <- "Chrono.txt"
}
sLowHigh <- c("low", "high")
objectivePaths <- file.path(result.dir, x$station, paste0(sLowHigh[x$flood + 1], "_", x$threshold))
l <- lapply(1:nrow(x), function(i) {
cat("Treating folder", objectivePaths[i], "...\n")
vgest_read_chrono.default(
file.path(objectivePaths[i], fileName),
nLakes = nrow(x[i, "lakes"][[1]]),
distributionType = distributionType
)
})
names(l) <- paste(x$station, sLowHigh[x$flood + 1], x$threshold, sep = "_")
class(l) <- c("ListChrono", class(l))
return(l)
}
R/vgest_read_qnat.R deleted 100644 → 0
View file @ 24721744
#' Read daily flow time series used by VGEST
#'
#' @details The file should be an ASCII file with tab separators, the first column contains the date in format `YYYYMMDD`.
#' The returned [data.frame] has the class attributes "Qnat" for further use with S3 methods.
#'
#' @param file Path to the file to read
#'
#' @return a [data.frame] with the content of the file
#' @export
#'
vgest_read_qnat <- function(file) {
Qnat <- read.delim(file)
Qnat$Dates <- as.POSIXct(as.character(Qnat$Dates), format = "%Y%m%d", tz = "UTC")
class(Qnat) <- c("Qnat", class(Qnat))
return(Qnat)
}
\ No newline at end of file
R/vgest_run.R deleted 100644 → 0
View file @ 24721744
#' Run VGEST and stop execution if an error is encountered during execution
#'
#' @details
#' Delete the content of the `RESULTAT` folder before running
#'
#' @param vgest_location Location of VGEST installation (Default "../vgest")
#'
#' @return Nil
#' @export
#'
#' @examples
#' \donttest{
#' # Run vgest with the current configuration
#' vgest_run()
#' }
vgest_run <- function(vgest_location = "../vgest") {
workingDir <- getwd()
setwd(vgest_location)