R/wasa_calibwrap.R
In tpilz/WasaEchseTools: Title: Handy tools for preprocessing and analysis of the hydrological models WASA-SED and ECHSE (limited to the WASA engine). Currently only meant for internal use! If you want to try it, handle with care!
#' Wrapper function for calibration of the WASA-SED model
#'
#' A wrapper function, executing a simulation of the WASA-SED model for a given
#' parameter realisation.
#'
#' @param pars A named vector of type numeric with the values of selected parameters
#' (directed to \code{\link[WasaEchseTools]{wasa_modify_pars}}).
#'
#' @param wasa_app Character string giving the system command of the WASA-SED application
#' (directed to \code{\link[WasaEchseTools]{wasa_run}}).
#'
#' @param sp_input_dir Character string of the directory containing the spatial WASA-SED
#' input. E.g. the output of \code{\link[lumpR]{db_wasa_input}} (argument 'dest_dir')
#' (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param meteo_dir Character string of the directory containing the time series input
#' files for the WASA-SED model. Requires the files temperature.dat, radiation.dat,
#' humidity.dat, and rain_daily.dat or rain_hourly.dat (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param dir_run Character specifying the directory for the model run with the current
#' parameter realisation (directed to \code{\link[WasaEchseTools]{wasa_run}}).
#' Default: A temporary directory created with \code{\link{tempfile}}.
#'
#' @param sim_start Object of class 'date' giving the start date of the simulation
#' (will be written into WASA-SED input file 'do.dat'; directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param sim_end Object of class 'date' giving the end date of the simulation
#' (will be written into WASA-SED input file 'do.dat' directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param resol Temporal resolution of the model simulations. Supported are: 'daily'
#' (default) and 'hourly'.
#'
#' @param warmup_start An object of class 'date' giving the start date of the warm-up period.
#' If \code{NULL} (default), the value in do.dat (lines 4 and 6) is used.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param radex_file Character string of the file (including path) of preprared
#' extraterrestrial radiation input (named 'extraterrestrial_radiation.dat')
#' (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}})).
#'
#' @param dat_streamflow OPTIONAL: Object of class 'xts' containing a time series of
#' streamflow at the catchment outlet in (m3/s) in the resolution of the model run.
#' NA values are allowed and will be discarded for goodness of fit calculations.
#' Only needed if \code{return_val = 'nse'}.
#'
#' @param dat_pr OPTIONAL: Object of class 'xts' containing a time series of catchment-wide
#' average precipitation in (m3) in the resolution of the model run. If not given (default),
#' it will be read (and calculated) from the WASA-SED input files if needed. However,
#' it is more efficient in terms of function execution time to specify it as input if needed.
#' Only needed if \code{return_val = 'hydInd'}.
#'
#' @param flood_thresh OPTIONAL: Numeric value giving the threshold in (m3/s) for the definition of
#' a flood event (directed to \code{\link[WasaEchseTools]{hydInd}}). Only needed if
#' \code{return_val = 'hydInd'}.
#'
#' @param thresh_zero OPTIONAL: Values of discharge in (m3/s) below this value will be treated as
#' zero flows (directed to \code{\link[WasaEchseTools]{hydInd}}). Only needed if
#' \code{return_val = 'hydInd'}.
#'
#' @param warmup_len Integer giving the length of the warm-up period in months. Default: 3.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param max_pre_runs Integer specifying the maximum number of warm-up iterations to be
#' applied. If the relative storage change is still larger than \code{storage_tolerance}
#' after \code{max_pre_runs} iterations, the warm-up will be aborted and the model be run
#' anyway. A warning will be issued. Default: 20. Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param storage_tolerance Numeric value giving the relative change of the model's water
#' storages between two connsecutive warm-up runs below which the warm-up will be
#' concluded and the actual model simulation be started. Default: 0.01.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param keep_warmup_states Logical value. Shall state storage files after finishing the warm-up
#' runs be stored for upcomming runs (i.e. *.stat files be copied into \code{sp_input_dir}/input)?
#' WARNING: Existing *.stat files will be overwritten!
#' Default: \code{FALSE}. This option might be useful for calibration runs as it
#' might reduce the number of necessary warm-up runs for each new parameter set.
#'
#' @param return_val Character vector specifying your choice of what this function
#' shall return. Default: 'river_flow'. See description of return value below.
#'
#' @param return_sp Logical flag specifying if certain output shall be given including
#' spatial variability at subbasin level instead of mere catchment specific values.
#' See description of return values below for spatial outputs. Default: \code{FALSE}.
#'
#' @param log_meta Character containg the name (and path) of a file into which information
#' about the function call will be written. See below for more information. Default:
#' \code{NULL}, i.e. no logile will be created. If the file already exists, the file
#' name to be created will be extended by a call to \code{\link{tempfile}}.
#'
#' @param keep_rundir Value of type \code{logical}. Shall directory \code{dir_run}
#' be retained (\code{TRUE}) or deleted (\code{FALSE}) after function execution?
#' Default: \code{FALSE}.
#'
#' @param keep_log Value of type \code{logical}. Shall the WASA-SED specific log
#' file (not to be confused with \code{log_meta}!) of the model run be written
#' (to \code{dir_run})? Default: \code{FALSE}. Will be ignored if \code{keep_rundir = FALSE}.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param error2warn Value of type \code{logical}. Shall runtime errors of the model be
#' reported as a warning instead of stopping this function with an error? If so, the
#' model run's log file will be saved and, in case of reasonable model output, this wrapper
#' function will proceed as usual. If no reasonable output could be found,
#' a value \code{NA} will be returned. Default: \code{FALSE}.
#' Also directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @details Function is a wrapper function, internally executing functions \code{\link[WasaEchseTools]{wasa_prep_runs}},
#' \code{\link[WasaEchseTools]{wasa_modify_pars}}, and \code{\link[WasaEchseTools]{wasa_run}}
#' and processing and returning the simulation output as specified. The function can
#' be employed by model calibration functions such as \code{\link[HydroBayes]{dream}}
#' or \code{\link[ppso]{optim_dds}}, or to execute single model runs within a single
#' function call.
#'
#' @return Function returns a named list or a single element (if \code{length(return_val) = 1})
#' of numeric value(s). Can be controlled by argument \code{return_val}. Currently
#' implemented are the options:
#'
#' river_flow_[mm|ts]: A single value of the catchment's simulated river outflow over
#' the specified simulation period in (mm) or a time series of class 'xts' in (m3/s).
#' Respects \code{return_sp} (only ts, river_flow_mm is only given for catchment outlet!).
#'
#' eta_mm[_ts]: A single value of the catchment's simulated amount of actual evapotranspiration
#' over the specified simulation period in (mm) or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' etp_mm[_ts]: A single value of the catchment's simulated amount of potential evapotranspiration
#' over the specified simulation period in (mm) or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' runoff_total_mm[_ts]: A single value of the catchment's simulated amount of total runoff,
#' i.e. runoff contribution into river(s), over the specified simulation period in (mm).
#' Respects \code{return_sp} or a time series of class 'xts' in (mm/timestep).
#'
#' runoff_gw_mm[_ts]: A single value of the catchment's simulated amount of groundwater runoff,
#' i.e. groundwater contribution into river(s), over the specified simulation period in (mm).
#' Respects \code{return_sp} or a time series of class 'xts' in (mm/timestep).
#'
#' runoff_sub_mm[_ts]: A single value of the catchment's simulated amount of subsurface runoff,
#' i.e. near-surface soil water (excl. groundwater!) contribution into river(s), over the specified
#' simulation period in (mm) or a time series of class 'xts' in (mm/timestep). Note: due to computational reasons, 'runoff_gw_mm' will
#' be given in addition if this value is set.
#' Respects \code{return_sp}.
#'
#' runoff_surf_mm[_ts]: A single value of the catchment's simulated amount of surface runoff,
#' i.e. surface water contribution into river(s), over the specified simulation period in (mm)
#' or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' hydInd: Named vector of hydrological indices calculated with function \code{\link[WasaEchseTools]{hydInd}}.
#' See function's doc for more information. This option requires the optional input arguments
#' \code{flood_thresh}, \code{thresh_zero}, and (optional) \code{dat_pr}.
#'
#' nse: A single numeric value giving the Nash-Sutcliffe efficiency of streamflow
#' simulation at the catchment outlet (a common goodness of fit measure of hydrological
#' model runs).
#'
#' kge: A single numeric value giving the Kling-Gupta efficiency of streamflow
#' simulation at the catchment outlet (see paper \url{https://doi.org/10.1016/j.jhydrol.2009.08.003}).
#'
#'
#' If argument \code{log_meta} was given, the logfile contains the following information:
#' parameter names and corresponding realisations, output element names and corresponding
#' values, \code{run_dir}: realisation of argument \code{dir_run}, \code{time_total}:
#' total time for function execution (i.e. until writing logfile, in secs.), \code{time_simrun}:
#' runtime of the simulation run (in secs.), \code{time_warmup}: total runtime for all
#' warm-up runs (in secs.), \code{warmup_iterations}: number of warm-up iterations,
#' \code{warmup_storchange}: relative storage change after the last warm-up iteration.
#'
#'
#' @note To avoid warm-up runs, set \code{max_pre_runs} or \code{warmup_len} to zero.
#'
#' \strong{Model's water balance}: runoff_total_mm = runoff_surf_mm + runoff_sub_mm + runoff_gw_mm.
#' Moreover, river_flow_mm should be slight less than (due to transmission losses and storage
#' changes; but all in all almost equal to) runoff_total_mm as long as there is no
#' artificial influence (e.g. from reservoirs).
#'
#' \strong{Initial states} can have a large influence, depending on setup and parameterisation.
#' Therefore, it is advisable to make some test runs with different warmup parameters (
#' \code{warmup_start}, \code{warmup_len}, \code{max_pre_runs}, \code{storage_tolerance},
#' \code{keep_warmup_states}) before starting a calibration run! Consecutive runs with
#' the same parameterisation should come up with the same results when \code{keep_warmup_states = TRUE}.
#'
#' If running this function in \strong{parallel} mode:
#' Setting argument \code{keep_warmup_states = TRUE} can have problematic side effects
#' due to parallel reading and (over-)writing of the same files. To accelerate
#' warmup anyway, you can generate default initial storage files by running a
#' default parametrisation over the warmup horizon until convergence is reached
#' and use the resulting storage files as initial storages for all calibration runs.
#' A similar problem arises When specifying argument \code{log_meta}. In that case,
#' create an empty initial logfile '\code{log_meta}' before starting the parallel
#' calibration routine to invoke the file name extension via \code{\link{tempfile}}
#' right away.
#'
#'
#' @author Tobias Pilz \email{tpilz@@uni-potsdam.de}
#'
#' @export
wasa_calibwrap <- function(
pars = NULL,
wasa_app = NULL,
sp_input_dir = NULL,
meteo_dir = NULL,
dir_run = paste0(tempfile(pattern = "wasa_calib_"), sep="/"),
sim_start = NULL,
sim_end = NULL,
resol = "daily",
warmup_start = NULL,
radex_file = NULL,
dat_streamflow = NULL,
dat_pr = NULL,
flood_thresh = NULL,
thresh_zero = NULL,
warmup_len = 3,
max_pre_runs = 20,
storage_tolerance = 0.01,
keep_warmup_states = FALSE,
return_val = "river_flow_ts",
return_sp = FALSE,
log_meta = NULL,
keep_rundir = FALSE,
keep_log = FALSE,
error2warn = FALSE
) {
time_start <- Sys.time()
return_na <- FALSE # set to true in tryCatch(...) if error2warn = T
# CHECKS #
if(resol == "daily") {
timestep <- 24
prec_file <- "rain_daily.dat"
} else if(resol == "hourly") {
timestep <- 1
prec_file <- "rain_hourly.dat"
}
if("hydInd" %in% return_val) {
if(!is.numeric(flood_thresh)) stop("Argument return_val = hydInd requires argument flood_thresh to be given!")
if(!is.numeric(thresh_zero)) stop("Argument return_val = hydInd requires argument thresh_zero to be given!")
}
if(any(c("nse", "kge") %in% return_val)) {
if(!is.xts(dat_streamflow)) stop("Argument return_val = {nse|kge} requires an object of class 'xts' for argument dat_streamflow!")
}
if(!is.null(log_meta)) {
f_log <- TRUE
logfile = log_meta
logdir <- sub(basename(logfile), "", logfile)
if(file.exists(logfile)) logfile <- tempfile(sub(".[a-z]+$", "", basename(logfile)), tmpdir = logdir, fileext = sub("^[a-zA-Z0-9_-]+", "", basename(logfile)))
} else {
f_log <- FALSE
logfile <- NULL
}
# determine output files to be produced by WASA
outfiles <- NULL; return_val_gr <- NULL
if(any(c("nse", "kge", "hydInd", "river_flow_ts", "river_flow_mm") %in% return_val))
outfiles <- c(outfiles, "river_flow")
if(any(str_detect(return_val, "runoff_surf"))) outfiles <- c(outfiles, "daily_total_overlandflow")
if(any(str_detect(return_val, "runoff_sub"))) {
outfiles <- c(outfiles, "daily_subsurface_runoff")
# needed here as in WASA subsurface runoff = lateral near-surface soil water runoff + groundwater runoff!
if(!any(str_detect(return_val, "runoff_gw"))) return_val <- c(return_val, "runoff_gw_mm")
}
if(any(str_detect(return_val, "runoff_gw"))) outfiles <- c(outfiles, "gw_discharge")
if(any(str_detect(return_val, "runoff_total|runoff_surf|runoff_sub|eta|etp")) && timestep < 24)
stop("Currently a 'return_val' of 'eta_mm[_ts]', 'etp_mm[_ts]', or 'runoff_total_mm[_ts]' can only be returned when running with daily resolution!")
if(any(str_detect(return_val, "eta"))) outfiles <- c(outfiles, "daily_actetranspiration")
return_val_gr <- c(return_val_gr, "daily_actetranspiration" = "eta")
if(any(str_detect(return_val, "etp"))) outfiles <- c(outfiles, "daily_potetranspiration")
return_val_gr <- c(return_val_gr, "daily_potetranspiration" = "etp")
if(any(str_detect(return_val, "runoff_total"))) outfiles <- c(outfiles, "daily_water_subbasin")
return_val_gr <- c(return_val_gr, "daily_water_subbasin" = "runoff_total")
# prepare run directory
wasa_prep_runs(sp_input_dir = sp_input_dir, meteo_dir = meteo_dir,
radex_file = radex_file, wasa_sim_dir = dir_run, proj_name = "",
sim_start = sim_start, sim_end = sim_end, timestep = timestep,
outfiles = outfiles)
# modify wasa input
if(!is.null(pars)) wasa_modify_pars(pars, paste(dir_run, "input", sep="/"))
# run wasa (including warmup)
wasa_run(dir_run, wasa_app, warmup_start, warmup_len, max_pre_runs, storage_tolerance,
keep_warmup_states = keep_warmup_states, log_meta = logfile, keep_log = keep_log,
error2warn = error2warn)
# check that everything went fine (in case of a WASA error there will be a file run_save* if error2warn = T)
if(length(dir(dir_run, pattern = "run_save_[a-z0-9]+.log", full.names = T))>0) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
unlink(logfile)
}
warning(paste0("There was a runtime error in model run in ", dir_run, ". NA will be returned!"))
return(NA)
}
# save warm-up state storage files
if(keep_warmup_states) file.copy(dir(paste(dir_run, "input", sep="/"), pattern = ".stat$|.stats", full.names = T), sp_input_dir, overwrite = T)
# get simulation data
dat_wasa <- NULL
if("river_flow" %in% outfiles) {
file_wasa <- paste(dir_run, "output/River_Flow.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
mutate(date = as.POSIXct(paste(year, day, sep="-"), "%Y-%j", tz ="UTC"), group = "river_flow") %>%
dplyr::select(-day, -year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("gw_discharge" %in% outfiles) {
file_wasa <- paste(dir_run, "output/gw_discharge.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "gw_discharge") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("daily_subsurface_runoff" %in% outfiles) {
file_wasa <- paste(dir_run, "output/daily_subsurface_runoff.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "runoff_sub") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("daily_total_overlandflow" %in% outfiles) {
file_wasa <- paste(dir_run, "output/daily_total_overlandflow.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "runoff_surf") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if(any(c("daily_actetranspiration", "daily_potetranspiration", "daily_water_subbasin") %in% outfiles)) {
outfiles_t <- outfiles[match(c("daily_actetranspiration", "daily_potetranspiration", "daily_water_subbasin"), outfiles, nomatch = 0)]
dat_wasa <- map_dfr(outfiles_t, function(f) {
file_wasa <- paste0(dir_run, "/output/", f, ".out")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_file %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = return_val_gr[grep(f, names(return_val_gr))]) %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas))
}) %>%
bind_rows(dat_wasa)
}
# ignore errors if desired
if(nrow(dat_wasa) == 0) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
}
if(!error2warn) {
stop(paste0("There could be no output extracted from a model run, dir_run = ", dir_run))
} else {
warning(paste0("No reasonable model output could be extracted from model run in ", dir_run, ". NA will be returned!"))
return(NA)
}
}
# prepare output according to specifications
out <- NULL
tryCatch({
if(any(str_detect(return_val, "hydInd"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_sim_xts <- xts(dat_tmp$value, dat_tmp$date)
# precipitation (model forcing); get catchment-wide value (area-weighted precipitation mean)
if(is.null(dat_pr)) {
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F)
dat_sub_area <- dat_sub_area[2, order(dat_sub_area[1,])]
sub_pars <- data.frame(object = paste0("sub_", 1:length(dat_sub_area)), area = dat_sub_area)
dat_prec <- left_join(read.table(paste(meteo_dir, prec_file, sep="/"), header=T, skip=2, check.names = F)[,-2] %>%
mutate(date = as.POSIXct(sprintf(.[[1]], fmt="%08d"), "%d%m%Y", tz="UTC")) %>%
dplyr::select(-1) %>%
melt(id.vars="date", variable.name = "object") %>%
mutate(object = paste0("sub_", object), variable = "precip"),
sub_pars %>%
mutate_if(is.factor, as.character),
by = "object") %>%
# in case of hourly resolution, calculate daily sums (WASA output River_Flow.out is only daily, even in case of hourly simulation resolution)
group_by(date, variable, object, area) %>%
summarise(value = sum(value)) %>%
group_by(date, variable) %>%
# sums over the catchment in m3
summarise(value = sum(value * area*1e3)) %>%
filter(date >= as.POSIXct(paste(sim_start, "00:00:00")) & date <= as.POSIXct(paste(sim_end, "23:59:59")))
dat_pr <- xts(dat_prec$value, dat_prec$date)
}
# calculate diagnostic values
out_vals <- suppressWarnings(hydInd(dat_sim_xts, dat_pr, na.rm = T, thresh.zero = thresh_zero, flood.thresh = flood_thresh))
out_vals[which(is.na(out_vals) | is.nan(out_vals))] <- 0
out[names(out_vals)] <- out_vals
out <- as.list(out)
}
if(any(str_detect(return_val, "river_flow"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow") %>%
spread(key = subbas, value = value) %>%
select(-group) %>%
rename_at(vars(-date), funs(paste0("river_flow_sub_", .)))
out_vals <- xts(select(dat_tmp, - date), dat_tmp$date)
if("river_flow_ts" %in% return_val) {
if(!return_sp) {
out_vals <- out_vals[,"river_flow_sub_1"]
colnames(out_vals) <- "river_flow"
}
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
}
if("river_flow_mm" %in% return_val) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
# get catchment area (m2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F)
dat_sub_area <- dat_sub_area[2, order(dat_sub_area[1,])]
dat_sub_area <- sum(dat_sub_area) * 1e6
# sum of catchment river outflow (mm)
outflow <- dat_tmp$value * timestep * 60*60 # m3/day
outflow <- sum(outflow) # m3
outflow <- outflow *1000 # L
outflow <- outflow / dat_sub_area # L/m2 = mm
out[["river_flow_mm"]] <- outflow
}
}
if(any(str_detect(return_val, "nse"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_nse <- left_join(select(dat_tmp, date, value),
data.frame(obs = dat_streamflow,
date = index(dat_streamflow)),
by = "date") %>%
drop_na(obs) %>%
mutate(diffsq = (value - obs)^2,
obsmean = mean(obs), diffsqobs = (obs - obsmean)^2) %>%
summarise(nse = 1 - sum(diffsq) / sum(diffsqobs))
out_vals <- as.numeric(dat_nse)
out[["nse"]] <- out_vals
}
if(any(str_detect(return_val, "kge"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_kge <- left_join(select(dat_tmp, date, value),
data.frame(obs = dat_streamflow,
date = index(dat_streamflow)),
by = "date") %>%
drop_na(obs) %>%
summarise(a = cor(obs, value) - 1, b = sd(value)/sd(obs) - 1, c = mean(value)/mean(obs) - 1,
kge = 1 - sqrt(a^2 + b^2 + c^2))
out_vals <- as.numeric(dat_kge$kge)
out[["kge"]] <- out_vals
}
if(any(str_detect(return_val, "runoff_gw"))) {
# raw: m3/timestep
dat_tmp <- filter(dat_wasa, group == "gw_discharge")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean for each timestep
dat_gw_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_gw_mm" %in% return_val) {
dat_gw_mm_sub <- dat_gw_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_gw_mm"]] <- unique(dat_gw_mm_sub$value_catch)
if(return_sp) out[paste("runoff_gw_mm_sub", dat_gw_mm$subbas, sep="_")] <- dat_gw_mm_sub$value
}
if("runoff_gw_mm_ts" %in% return_val) {
dat_gw_mm_ts_catch <- dat_gw_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_gw_mm_ts_catch$value, dat_gw_mm_ts_catch$date)
colnames(out_vals) <- "runoff_gw_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_gw_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_gw_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_sub"))) {
# raw: m3/day
dat_tmp <- filter(dat_wasa, group == "runoff_sub")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_sub_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
# NOTE: in WASA subsurface runoff = lateral near-surface soil water runoff + groundwater runoff!
dat_sub_mm_ts$value <- dat_sub_mm_ts$value - dat_gw_mm_ts$value
if("runoff_sub_mm" %in% return_val) {
dat_sub_mm_sub <- dat_sub_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_sub_mm"]] <- unique(dat_sub_mm_sub$value_catch)
if(return_sp) out[paste("runoff_sub_mm_sub", dat_sub_mm$subbas, sep="_")] <- dat_sub_mm_sub$value
}
if("runoff_sub_mm_ts" %in% return_val) {
dat_sub_mm_ts_catch <- dat_sub_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_sub_mm_ts_catch$value, dat_sub_mm_ts_catch$date)
colnames(out_vals) <- "runoff_sub_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_sub_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_sub_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_surf"))) {
# raw: m3/timestep
dat_tmp <- filter(dat_wasa, group == "runoff_surf")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_surf_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_surf_mm" %in% return_val) {
dat_surf_mm_sub <- dat_surf_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_surf_mm"]] <- unique(dat_surf_mm_sub$value_catch)
if(return_sp) out[paste("runoff_surf_mm_sub", dat_surf_mm_ts$subbas, sep="_")] <- dat_surf_mm_sub$value
}
if("runoff_surf_mm_ts" %in% return_val) {
dat_surf_mm_ts_catch <- dat_surf_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_surf_mm_ts_catch$value, dat_surf_mm_ts_catch$date)
colnames(out_vals) <- "runoff_surf_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_surf_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_surf_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_total"))) {
# raw: m3/s
dat_tmp <- filter(dat_wasa, group == "runoff_total")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# sum of catchment river outflow (mm)
outflow <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value*60*60*24 / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_total_mm" %in% return_val) {
outflow_sub <- outflow %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_total_mm"]] <- unique(outflow_sub$value_catch)
if(return_sp) out[paste("runoff_total_mm_sub", outflow$subbas, sep="_")] <- outflow_sub$value
}
if("runoff_total_mm_ts" %in% return_val) {
outflow_catch <- outflow %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(outflow_catch$value, outflow_catch$date)
colnames(out_vals) <- "runoff_total_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- outflow %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_total_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "eta|etp"))) {
# raw: mm/day
dat_tmp <- filter(dat_wasa, group %in% c("eta", "etp"))
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_sums <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area), area_weight = area / area_sum),
by = "subbas") %>%
arrange(group,subbas) # mm / day
if(any(c("eta_mm", "etp_mm") %in% return_val)) {
dat_sums_sub <- dat_sums %>%
group_by(group, subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
group_by(group) %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[unique(dat_sums_sub$group)] <- unique(dat_sums_sub$value_catch)
if(return_sp) out[paste0(dat_sums_sub$group, "_sub_", dat_sums_sub$subbas)] <- dat_sums_sub$value
}
if(any(str_detect(grep("eta|etp", return_val, value = T), "_ts"))) {
dat_sums_catch <- dat_sums %>%
group_by(group, date) %>%
summarise(value = sum(value*area / area_sum)) %>% # area-weighted catchment sum in (mm)
spread(key = group, value = value) %>%
rename_at(vars(-date), funs(paste0(., "_mm")))
out_vals <- xts(select(dat_sums_catch, -date), dat_sums_catch$date)
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_sums %>%
select(-area, -area_sum, -area_weight) %>%
unite(col = "var", group, subbas, sep="_mm_sub_") %>%
spread(key = var, value = value)
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(length(out) > 1 && class(out) != "list") out <- as.list(out)
}, error = function(e) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
}
if(!error2warn) {
stop(paste0("Error while compiling output for model run, dir_run = ", dir_run))
} else {
warning(paste0("Could not compile output for model run in ", dir_run, ". NA will be returned!"))
return_na <- TRUE
}
})
if(return_na) return(NA)
# clean up
if(!keep_rundir) unlink(dir_run, recursive = T)
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
if(any(names(out) == "xts")) {
logfile_xts <- str_replace(logfile, ".dat$", "_xts.dat")
write.zoo(out$xts, file=logfile_xts, index.name = "date", row.names = F, col.names = T, sep="\t", quote=F)
outdat_log <- out[-grep("xts", names(out))]
}
if(length(outdat_log) > 0) {
out_log_t <- data.frame(group = c(rep("pars", length(pars)), rep("output", length(outdat_log)), "meta"),
variable = c(names(pars), names(outdat_log), "time_total"),
value = c(round(pars, 4), round(unlist(outdat_log),4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F)
} else {
out_log_t <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F)
}
out_log <- out_log_t %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=logfile, sep="\t", quote=F, row.names=F, col.names=T)
}
# catch non-finite output
if(any(!is.finite(unlist(out)))) {
warning("Non-finite outut detected! Return NA. Consider argument 'log_meta'.")
return(NA)
}
# output
return(out)
} # EOF
tpilz/WasaEchseTools documentation built on May 5, 2019, 12:33 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Wrapper function for calibration of the WASA-SED model
#'
#' A wrapper function, executing a simulation of the WASA-SED model for a given
#' parameter realisation.
#'
#' @param pars A named vector of type numeric with the values of selected parameters
#' (directed to \code{\link[WasaEchseTools]{wasa_modify_pars}}).
#'
#' @param wasa_app Character string giving the system command of the WASA-SED application
#' (directed to \code{\link[WasaEchseTools]{wasa_run}}).
#'
#' @param sp_input_dir Character string of the directory containing the spatial WASA-SED
#' input. E.g. the output of \code{\link[lumpR]{db_wasa_input}} (argument 'dest_dir')
#' (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param meteo_dir Character string of the directory containing the time series input
#' files for the WASA-SED model. Requires the files temperature.dat, radiation.dat,
#' humidity.dat, and rain_daily.dat or rain_hourly.dat (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param dir_run Character specifying the directory for the model run with the current
#' parameter realisation (directed to \code{\link[WasaEchseTools]{wasa_run}}).
#' Default: A temporary directory created with \code{\link{tempfile}}.
#'
#' @param sim_start Object of class 'date' giving the start date of the simulation
#' (will be written into WASA-SED input file 'do.dat'; directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param sim_end Object of class 'date' giving the end date of the simulation
#' (will be written into WASA-SED input file 'do.dat' directed to \code{\link[WasaEchseTools]{wasa_prep_runs}}).
#'
#' @param resol Temporal resolution of the model simulations. Supported are: 'daily'
#' (default) and 'hourly'.
#'
#' @param warmup_start An object of class 'date' giving the start date of the warm-up period.
#' If \code{NULL} (default), the value in do.dat (lines 4 and 6) is used.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param radex_file Character string of the file (including path) of preprared
#' extraterrestrial radiation input (named 'extraterrestrial_radiation.dat')
#' (directed to \code{\link[WasaEchseTools]{wasa_prep_runs}})).
#'
#' @param dat_streamflow OPTIONAL: Object of class 'xts' containing a time series of
#' streamflow at the catchment outlet in (m3/s) in the resolution of the model run.
#' NA values are allowed and will be discarded for goodness of fit calculations.
#' Only needed if \code{return_val = 'nse'}.
#'
#' @param dat_pr OPTIONAL: Object of class 'xts' containing a time series of catchment-wide
#' average precipitation in (m3) in the resolution of the model run. If not given (default),
#' it will be read (and calculated) from the WASA-SED input files if needed. However,
#' it is more efficient in terms of function execution time to specify it as input if needed.
#' Only needed if \code{return_val = 'hydInd'}.
#'
#' @param flood_thresh OPTIONAL: Numeric value giving the threshold in (m3/s) for the definition of
#' a flood event (directed to \code{\link[WasaEchseTools]{hydInd}}). Only needed if
#' \code{return_val = 'hydInd'}.
#'
#' @param thresh_zero OPTIONAL: Values of discharge in (m3/s) below this value will be treated as
#' zero flows (directed to \code{\link[WasaEchseTools]{hydInd}}). Only needed if
#' \code{return_val = 'hydInd'}.
#'
#' @param warmup_len Integer giving the length of the warm-up period in months. Default: 3.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param max_pre_runs Integer specifying the maximum number of warm-up iterations to be
#' applied. If the relative storage change is still larger than \code{storage_tolerance}
#' after \code{max_pre_runs} iterations, the warm-up will be aborted and the model be run
#' anyway. A warning will be issued. Default: 20. Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param storage_tolerance Numeric value giving the relative change of the model's water
#' storages between two connsecutive warm-up runs below which the warm-up will be
#' concluded and the actual model simulation be started. Default: 0.01.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param keep_warmup_states Logical value. Shall state storage files after finishing the warm-up
#' runs be stored for upcomming runs (i.e. *.stat files be copied into \code{sp_input_dir}/input)?
#' WARNING: Existing *.stat files will be overwritten!
#' Default: \code{FALSE}. This option might be useful for calibration runs as it
#' might reduce the number of necessary warm-up runs for each new parameter set.
#'
#' @param return_val Character vector specifying your choice of what this function
#' shall return. Default: 'river_flow'. See description of return value below.
#'
#' @param return_sp Logical flag specifying if certain output shall be given including
#' spatial variability at subbasin level instead of mere catchment specific values.
#' See description of return values below for spatial outputs. Default: \code{FALSE}.
#'
#' @param log_meta Character containg the name (and path) of a file into which information
#' about the function call will be written. See below for more information. Default:
#' \code{NULL}, i.e. no logile will be created. If the file already exists, the file
#' name to be created will be extended by a call to \code{\link{tempfile}}.
#'
#' @param keep_rundir Value of type \code{logical}. Shall directory \code{dir_run}
#' be retained (\code{TRUE}) or deleted (\code{FALSE}) after function execution?
#' Default: \code{FALSE}.
#'
#' @param keep_log Value of type \code{logical}. Shall the WASA-SED specific log
#' file (not to be confused with \code{log_meta}!) of the model run be written
#' (to \code{dir_run})? Default: \code{FALSE}. Will be ignored if \code{keep_rundir = FALSE}.
#' Directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @param error2warn Value of type \code{logical}. Shall runtime errors of the model be
#' reported as a warning instead of stopping this function with an error? If so, the
#' model run's log file will be saved and, in case of reasonable model output, this wrapper
#' function will proceed as usual. If no reasonable output could be found,
#' a value \code{NA} will be returned. Default: \code{FALSE}.
#' Also directed to \code{\link[WasaEchseTools]{wasa_run}}.
#'
#' @details Function is a wrapper function, internally executing functions \code{\link[WasaEchseTools]{wasa_prep_runs}},
#' \code{\link[WasaEchseTools]{wasa_modify_pars}}, and \code{\link[WasaEchseTools]{wasa_run}}
#' and processing and returning the simulation output as specified. The function can
#' be employed by model calibration functions such as \code{\link[HydroBayes]{dream}}
#' or \code{\link[ppso]{optim_dds}}, or to execute single model runs within a single
#' function call.
#'
#' @return Function returns a named list or a single element (if \code{length(return_val) = 1})
#' of numeric value(s). Can be controlled by argument \code{return_val}. Currently
#' implemented are the options:
#'
#' river_flow_[mm|ts]: A single value of the catchment's simulated river outflow over
#' the specified simulation period in (mm) or a time series of class 'xts' in (m3/s).
#' Respects \code{return_sp} (only ts, river_flow_mm is only given for catchment outlet!).
#'
#' eta_mm[_ts]: A single value of the catchment's simulated amount of actual evapotranspiration
#' over the specified simulation period in (mm) or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' etp_mm[_ts]: A single value of the catchment's simulated amount of potential evapotranspiration
#' over the specified simulation period in (mm) or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' runoff_total_mm[_ts]: A single value of the catchment's simulated amount of total runoff,
#' i.e. runoff contribution into river(s), over the specified simulation period in (mm).
#' Respects \code{return_sp} or a time series of class 'xts' in (mm/timestep).
#'
#' runoff_gw_mm[_ts]: A single value of the catchment's simulated amount of groundwater runoff,
#' i.e. groundwater contribution into river(s), over the specified simulation period in (mm).
#' Respects \code{return_sp} or a time series of class 'xts' in (mm/timestep).
#'
#' runoff_sub_mm[_ts]: A single value of the catchment's simulated amount of subsurface runoff,
#' i.e. near-surface soil water (excl. groundwater!) contribution into river(s), over the specified
#' simulation period in (mm) or a time series of class 'xts' in (mm/timestep). Note: due to computational reasons, 'runoff_gw_mm' will
#' be given in addition if this value is set.
#' Respects \code{return_sp}.
#'
#' runoff_surf_mm[_ts]: A single value of the catchment's simulated amount of surface runoff,
#' i.e. surface water contribution into river(s), over the specified simulation period in (mm)
#' or a time series of class 'xts' in (mm/timestep).
#' Respects \code{return_sp}.
#'
#' hydInd: Named vector of hydrological indices calculated with function \code{\link[WasaEchseTools]{hydInd}}.
#' See function's doc for more information. This option requires the optional input arguments
#' \code{flood_thresh}, \code{thresh_zero}, and (optional) \code{dat_pr}.
#'
#' nse: A single numeric value giving the Nash-Sutcliffe efficiency of streamflow
#' simulation at the catchment outlet (a common goodness of fit measure of hydrological
#' model runs).
#'
#' kge: A single numeric value giving the Kling-Gupta efficiency of streamflow
#' simulation at the catchment outlet (see paper \url{https://doi.org/10.1016/j.jhydrol.2009.08.003}).
#'
#'
#' If argument \code{log_meta} was given, the logfile contains the following information:
#' parameter names and corresponding realisations, output element names and corresponding
#' values, \code{run_dir}: realisation of argument \code{dir_run}, \code{time_total}:
#' total time for function execution (i.e. until writing logfile, in secs.), \code{time_simrun}:
#' runtime of the simulation run (in secs.), \code{time_warmup}: total runtime for all
#' warm-up runs (in secs.), \code{warmup_iterations}: number of warm-up iterations,
#' \code{warmup_storchange}: relative storage change after the last warm-up iteration.
#'
#'
#' @note To avoid warm-up runs, set \code{max_pre_runs} or \code{warmup_len} to zero.
#'
#' \strong{Model's water balance}: runoff_total_mm = runoff_surf_mm + runoff_sub_mm + runoff_gw_mm.
#' Moreover, river_flow_mm should be slight less than (due to transmission losses and storage
#' changes; but all in all almost equal to) runoff_total_mm as long as there is no
#' artificial influence (e.g. from reservoirs).
#'
#' \strong{Initial states} can have a large influence, depending on setup and parameterisation.
#' Therefore, it is advisable to make some test runs with different warmup parameters (
#' \code{warmup_start}, \code{warmup_len}, \code{max_pre_runs}, \code{storage_tolerance},
#' \code{keep_warmup_states}) before starting a calibration run! Consecutive runs with
#' the same parameterisation should come up with the same results when \code{keep_warmup_states = TRUE}.
#'
#' If running this function in \strong{parallel} mode:
#' Setting argument \code{keep_warmup_states = TRUE} can have problematic side effects
#' due to parallel reading and (over-)writing of the same files. To accelerate
#' warmup anyway, you can generate default initial storage files by running a
#' default parametrisation over the warmup horizon until convergence is reached
#' and use the resulting storage files as initial storages for all calibration runs.
#' A similar problem arises When specifying argument \code{log_meta}. In that case,
#' create an empty initial logfile '\code{log_meta}' before starting the parallel
#' calibration routine to invoke the file name extension via \code{\link{tempfile}}
#' right away.
#'
#'
#' @author Tobias Pilz \email{tpilz@@uni-potsdam.de}
#'
#' @export
wasa_calibwrap <- function(
pars = NULL,
wasa_app = NULL,
sp_input_dir = NULL,
meteo_dir = NULL,
dir_run = paste0(tempfile(pattern = "wasa_calib_"), sep="/"),
sim_start = NULL,
sim_end = NULL,
resol = "daily",
warmup_start = NULL,
radex_file = NULL,
dat_streamflow = NULL,
dat_pr = NULL,
flood_thresh = NULL,
thresh_zero = NULL,
warmup_len = 3,
max_pre_runs = 20,
storage_tolerance = 0.01,
keep_warmup_states = FALSE,
return_val = "river_flow_ts",
return_sp = FALSE,
log_meta = NULL,
keep_rundir = FALSE,
keep_log = FALSE,
error2warn = FALSE
) {
time_start <- Sys.time()
return_na <- FALSE # set to true in tryCatch(...) if error2warn = T
# CHECKS #
if(resol == "daily") {
timestep <- 24
prec_file <- "rain_daily.dat"
} else if(resol == "hourly") {
timestep <- 1
prec_file <- "rain_hourly.dat"
}
if("hydInd" %in% return_val) {
if(!is.numeric(flood_thresh)) stop("Argument return_val = hydInd requires argument flood_thresh to be given!")
if(!is.numeric(thresh_zero)) stop("Argument return_val = hydInd requires argument thresh_zero to be given!")
}
if(any(c("nse", "kge") %in% return_val)) {
if(!is.xts(dat_streamflow)) stop("Argument return_val = {nse|kge} requires an object of class 'xts' for argument dat_streamflow!")
}
if(!is.null(log_meta)) {
f_log <- TRUE
logfile = log_meta
logdir <- sub(basename(logfile), "", logfile)
if(file.exists(logfile)) logfile <- tempfile(sub(".[a-z]+$", "", basename(logfile)), tmpdir = logdir, fileext = sub("^[a-zA-Z0-9_-]+", "", basename(logfile)))
} else {
f_log <- FALSE
logfile <- NULL
}
# determine output files to be produced by WASA
outfiles <- NULL; return_val_gr <- NULL
if(any(c("nse", "kge", "hydInd", "river_flow_ts", "river_flow_mm") %in% return_val))
outfiles <- c(outfiles, "river_flow")
if(any(str_detect(return_val, "runoff_surf"))) outfiles <- c(outfiles, "daily_total_overlandflow")
if(any(str_detect(return_val, "runoff_sub"))) {
outfiles <- c(outfiles, "daily_subsurface_runoff")
# needed here as in WASA subsurface runoff = lateral near-surface soil water runoff + groundwater runoff!
if(!any(str_detect(return_val, "runoff_gw"))) return_val <- c(return_val, "runoff_gw_mm")
}
if(any(str_detect(return_val, "runoff_gw"))) outfiles <- c(outfiles, "gw_discharge")
if(any(str_detect(return_val, "runoff_total|runoff_surf|runoff_sub|eta|etp")) && timestep < 24)
stop("Currently a 'return_val' of 'eta_mm[_ts]', 'etp_mm[_ts]', or 'runoff_total_mm[_ts]' can only be returned when running with daily resolution!")
if(any(str_detect(return_val, "eta"))) outfiles <- c(outfiles, "daily_actetranspiration")
return_val_gr <- c(return_val_gr, "daily_actetranspiration" = "eta")
if(any(str_detect(return_val, "etp"))) outfiles <- c(outfiles, "daily_potetranspiration")
return_val_gr <- c(return_val_gr, "daily_potetranspiration" = "etp")
if(any(str_detect(return_val, "runoff_total"))) outfiles <- c(outfiles, "daily_water_subbasin")
return_val_gr <- c(return_val_gr, "daily_water_subbasin" = "runoff_total")
# prepare run directory
wasa_prep_runs(sp_input_dir = sp_input_dir, meteo_dir = meteo_dir,
radex_file = radex_file, wasa_sim_dir = dir_run, proj_name = "",
sim_start = sim_start, sim_end = sim_end, timestep = timestep,
outfiles = outfiles)
# modify wasa input
if(!is.null(pars)) wasa_modify_pars(pars, paste(dir_run, "input", sep="/"))
# run wasa (including warmup)
wasa_run(dir_run, wasa_app, warmup_start, warmup_len, max_pre_runs, storage_tolerance,
keep_warmup_states = keep_warmup_states, log_meta = logfile, keep_log = keep_log,
error2warn = error2warn)
# check that everything went fine (in case of a WASA error there will be a file run_save* if error2warn = T)
if(length(dir(dir_run, pattern = "run_save_[a-z0-9]+.log", full.names = T))>0) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
unlink(logfile)
}
warning(paste0("There was a runtime error in model run in ", dir_run, ". NA will be returned!"))
return(NA)
}
# save warm-up state storage files
if(keep_warmup_states) file.copy(dir(paste(dir_run, "input", sep="/"), pattern = ".stat$|.stats", full.names = T), sp_input_dir, overwrite = T)
# get simulation data
dat_wasa <- NULL
if("river_flow" %in% outfiles) {
file_wasa <- paste(dir_run, "output/River_Flow.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
mutate(date = as.POSIXct(paste(year, day, sep="-"), "%Y-%j", tz ="UTC"), group = "river_flow") %>%
dplyr::select(-day, -year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("gw_discharge" %in% outfiles) {
file_wasa <- paste(dir_run, "output/gw_discharge.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "gw_discharge") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("daily_subsurface_runoff" %in% outfiles) {
file_wasa <- paste(dir_run, "output/daily_subsurface_runoff.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "runoff_sub") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if("daily_total_overlandflow" %in% outfiles) {
file_wasa <- paste(dir_run, "output/daily_total_overlandflow.out", sep="/")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_wasa <- dat_file %>%
select(Day, Year, num_range(prefix = "", range = 1:999)) %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = "runoff_surf") %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas)) %>%
bind_rows(dat_wasa)
}
if(any(c("daily_actetranspiration", "daily_potetranspiration", "daily_water_subbasin") %in% outfiles)) {
outfiles_t <- outfiles[match(c("daily_actetranspiration", "daily_potetranspiration", "daily_water_subbasin"), outfiles, nomatch = 0)]
dat_wasa <- map_dfr(outfiles_t, function(f) {
file_wasa <- paste0(dir_run, "/output/", f, ".out")
tryCatch(dat_file <- read.table(file_wasa, header=T, skip=1, check.names = F),
error = function(e) stop(paste("Error reading", file_wasa, ":", e)))
dat_file %>%
mutate(date = as.POSIXct(paste(Year, Day, sep="-"), "%Y-%j", tz ="UTC"), group = return_val_gr[grep(f, names(return_val_gr))]) %>%
select(-Day, -Year) %>%
gather(key = subbas, value = value, -group, -date) %>%
mutate(subbas = as.integer(subbas))
}) %>%
bind_rows(dat_wasa)
}
# ignore errors if desired
if(nrow(dat_wasa) == 0) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
}
if(!error2warn) {
stop(paste0("There could be no output extracted from a model run, dir_run = ", dir_run))
} else {
warning(paste0("No reasonable model output could be extracted from model run in ", dir_run, ". NA will be returned!"))
return(NA)
}
}
# prepare output according to specifications
out <- NULL
tryCatch({
if(any(str_detect(return_val, "hydInd"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_sim_xts <- xts(dat_tmp$value, dat_tmp$date)
# precipitation (model forcing); get catchment-wide value (area-weighted precipitation mean)
if(is.null(dat_pr)) {
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F)
dat_sub_area <- dat_sub_area[2, order(dat_sub_area[1,])]
sub_pars <- data.frame(object = paste0("sub_", 1:length(dat_sub_area)), area = dat_sub_area)
dat_prec <- left_join(read.table(paste(meteo_dir, prec_file, sep="/"), header=T, skip=2, check.names = F)[,-2] %>%
mutate(date = as.POSIXct(sprintf(.[[1]], fmt="%08d"), "%d%m%Y", tz="UTC")) %>%
dplyr::select(-1) %>%
melt(id.vars="date", variable.name = "object") %>%
mutate(object = paste0("sub_", object), variable = "precip"),
sub_pars %>%
mutate_if(is.factor, as.character),
by = "object") %>%
# in case of hourly resolution, calculate daily sums (WASA output River_Flow.out is only daily, even in case of hourly simulation resolution)
group_by(date, variable, object, area) %>%
summarise(value = sum(value)) %>%
group_by(date, variable) %>%
# sums over the catchment in m3
summarise(value = sum(value * area*1e3)) %>%
filter(date >= as.POSIXct(paste(sim_start, "00:00:00")) & date <= as.POSIXct(paste(sim_end, "23:59:59")))
dat_pr <- xts(dat_prec$value, dat_prec$date)
}
# calculate diagnostic values
out_vals <- suppressWarnings(hydInd(dat_sim_xts, dat_pr, na.rm = T, thresh.zero = thresh_zero, flood.thresh = flood_thresh))
out_vals[which(is.na(out_vals) | is.nan(out_vals))] <- 0
out[names(out_vals)] <- out_vals
out <- as.list(out)
}
if(any(str_detect(return_val, "river_flow"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow") %>%
spread(key = subbas, value = value) %>%
select(-group) %>%
rename_at(vars(-date), funs(paste0("river_flow_sub_", .)))
out_vals <- xts(select(dat_tmp, - date), dat_tmp$date)
if("river_flow_ts" %in% return_val) {
if(!return_sp) {
out_vals <- out_vals[,"river_flow_sub_1"]
colnames(out_vals) <- "river_flow"
}
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
}
if("river_flow_mm" %in% return_val) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
# get catchment area (m2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F)
dat_sub_area <- dat_sub_area[2, order(dat_sub_area[1,])]
dat_sub_area <- sum(dat_sub_area) * 1e6
# sum of catchment river outflow (mm)
outflow <- dat_tmp$value * timestep * 60*60 # m3/day
outflow <- sum(outflow) # m3
outflow <- outflow *1000 # L
outflow <- outflow / dat_sub_area # L/m2 = mm
out[["river_flow_mm"]] <- outflow
}
}
if(any(str_detect(return_val, "nse"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_nse <- left_join(select(dat_tmp, date, value),
data.frame(obs = dat_streamflow,
date = index(dat_streamflow)),
by = "date") %>%
drop_na(obs) %>%
mutate(diffsq = (value - obs)^2,
obsmean = mean(obs), diffsqobs = (obs - obsmean)^2) %>%
summarise(nse = 1 - sum(diffsq) / sum(diffsqobs))
out_vals <- as.numeric(dat_nse)
out[["nse"]] <- out_vals
}
if(any(str_detect(return_val, "kge"))) {
dat_tmp <- filter(dat_wasa, group == "river_flow" & subbas == 1)
dat_kge <- left_join(select(dat_tmp, date, value),
data.frame(obs = dat_streamflow,
date = index(dat_streamflow)),
by = "date") %>%
drop_na(obs) %>%
summarise(a = cor(obs, value) - 1, b = sd(value)/sd(obs) - 1, c = mean(value)/mean(obs) - 1,
kge = 1 - sqrt(a^2 + b^2 + c^2))
out_vals <- as.numeric(dat_kge$kge)
out[["kge"]] <- out_vals
}
if(any(str_detect(return_val, "runoff_gw"))) {
# raw: m3/timestep
dat_tmp <- filter(dat_wasa, group == "gw_discharge")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean for each timestep
dat_gw_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_gw_mm" %in% return_val) {
dat_gw_mm_sub <- dat_gw_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_gw_mm"]] <- unique(dat_gw_mm_sub$value_catch)
if(return_sp) out[paste("runoff_gw_mm_sub", dat_gw_mm$subbas, sep="_")] <- dat_gw_mm_sub$value
}
if("runoff_gw_mm_ts" %in% return_val) {
dat_gw_mm_ts_catch <- dat_gw_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_gw_mm_ts_catch$value, dat_gw_mm_ts_catch$date)
colnames(out_vals) <- "runoff_gw_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_gw_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_gw_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_sub"))) {
# raw: m3/day
dat_tmp <- filter(dat_wasa, group == "runoff_sub")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_sub_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
# NOTE: in WASA subsurface runoff = lateral near-surface soil water runoff + groundwater runoff!
dat_sub_mm_ts$value <- dat_sub_mm_ts$value - dat_gw_mm_ts$value
if("runoff_sub_mm" %in% return_val) {
dat_sub_mm_sub <- dat_sub_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_sub_mm"]] <- unique(dat_sub_mm_sub$value_catch)
if(return_sp) out[paste("runoff_sub_mm_sub", dat_sub_mm$subbas, sep="_")] <- dat_sub_mm_sub$value
}
if("runoff_sub_mm_ts" %in% return_val) {
dat_sub_mm_ts_catch <- dat_sub_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_sub_mm_ts_catch$value, dat_sub_mm_ts_catch$date)
colnames(out_vals) <- "runoff_sub_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_sub_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_sub_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_surf"))) {
# raw: m3/timestep
dat_tmp <- filter(dat_wasa, group == "runoff_surf")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_surf_mm_ts <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_surf_mm" %in% return_val) {
dat_surf_mm_sub <- dat_surf_mm_ts %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_surf_mm"]] <- unique(dat_surf_mm_sub$value_catch)
if(return_sp) out[paste("runoff_surf_mm_sub", dat_surf_mm_ts$subbas, sep="_")] <- dat_surf_mm_sub$value
}
if("runoff_surf_mm_ts" %in% return_val) {
dat_surf_mm_ts_catch <- dat_surf_mm_ts %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(dat_surf_mm_ts_catch$value, dat_surf_mm_ts_catch$date)
colnames(out_vals) <- "runoff_surf_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_surf_mm_ts %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_surf_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "runoff_total"))) {
# raw: m3/s
dat_tmp <- filter(dat_wasa, group == "runoff_total")
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# sum of catchment river outflow (mm)
outflow <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area)),
by = "subbas") %>%
mutate(value = value*60*60*24 / (area*1000)) %>% # mm / day
arrange(subbas)
if("runoff_total_mm" %in% return_val) {
outflow_sub <- outflow %>%
group_by(subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
ungroup() %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[["runoff_total_mm"]] <- unique(outflow_sub$value_catch)
if(return_sp) out[paste("runoff_total_mm_sub", outflow$subbas, sep="_")] <- outflow_sub$value
}
if("runoff_total_mm_ts" %in% return_val) {
outflow_catch <- outflow %>%
group_by(date) %>%
summarise(value = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out_vals <- xts(outflow_catch$value, outflow_catch$date)
colnames(out_vals) <- "runoff_total_mm"
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- outflow %>%
select(-area, -area_sum, -group) %>%
spread(key = subbas, value = value) %>%
rename_at(vars(-date), funs(paste0("runoff_total_mm_sub_", .)))
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(any(str_detect(return_val, "eta|etp"))) {
# raw: mm/day
dat_tmp <- filter(dat_wasa, group %in% c("eta", "etp"))
# get catchment area (km2)
dat_sub <- readLines(paste(dir_run, "input/Hillslope/hymo.dat", sep="/"))
dat_sub <- dat_sub[-c(1,2)]
dat_sub_area <- t(sapply(dat_sub, function(x) as.numeric(unlist(strsplit(x, "\t"))[c(1,2)]), USE.NAMES = F))
colnames(dat_sub_area) <- c("subbas", "area")
# merge data and calculate area-weighted mean and time period sum
dat_sums <- left_join(dat_tmp,
as.data.frame(dat_sub_area) %>%
mutate(subbas = as.integer(subbas), area_sum = sum(area), area_weight = area / area_sum),
by = "subbas") %>%
arrange(group,subbas) # mm / day
if(any(c("eta_mm", "etp_mm") %in% return_val)) {
dat_sums_sub <- dat_sums %>%
group_by(group, subbas, area, area_sum) %>%
summarise(value = sum(value)) %>% # sum over simulation period (mm)
group_by(group) %>%
mutate(value_catch = sum(value*area / area_sum)) # area-weighted catchment sum in (mm)
out[unique(dat_sums_sub$group)] <- unique(dat_sums_sub$value_catch)
if(return_sp) out[paste0(dat_sums_sub$group, "_sub_", dat_sums_sub$subbas)] <- dat_sums_sub$value
}
if(any(str_detect(grep("eta|etp", return_val, value = T), "_ts"))) {
dat_sums_catch <- dat_sums %>%
group_by(group, date) %>%
summarise(value = sum(value*area / area_sum)) %>% # area-weighted catchment sum in (mm)
spread(key = group, value = value) %>%
rename_at(vars(-date), funs(paste0(., "_mm")))
out_vals <- xts(select(dat_sums_catch, -date), dat_sums_catch$date)
if(any(names(out) == "xts")) {
out[["xts"]] <- cbind(out$xts, out_vals)
} else {
out[["xts"]] <- out_vals
}
if(return_sp) {
out_vals <- dat_sums %>%
select(-area, -area_sum, -area_weight) %>%
unite(col = "var", group, subbas, sep="_mm_sub_") %>%
spread(key = var, value = value)
out_vals <- xts(select(out_vals, -date), out_vals$date)
out[["xts"]] <- cbind(out$xts, out_vals)
}
}
}
if(length(out) > 1 && class(out) != "list") out <- as.list(out)
}, error = function(e) {
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
out_log <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F
) %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=paste0(logfile, ".err"), sep="\t", quote=F, row.names=F, col.names=T)
}
if(!error2warn) {
stop(paste0("Error while compiling output for model run, dir_run = ", dir_run))
} else {
warning(paste0("Could not compile output for model run in ", dir_run, ". NA will be returned!"))
return_na <- TRUE
}
})
if(return_na) return(NA)
# clean up
if(!keep_rundir) unlink(dir_run, recursive = T)
# write logfile
if(f_log) {
# read information from call to wasa_run()
dat_log <- read.table(logfile, header =T, stringsAsFactors = F)
time_end <- Sys.time()
if(any(names(out) == "xts")) {
logfile_xts <- str_replace(logfile, ".dat$", "_xts.dat")
write.zoo(out$xts, file=logfile_xts, index.name = "date", row.names = F, col.names = T, sep="\t", quote=F)
outdat_log <- out[-grep("xts", names(out))]
}
if(length(outdat_log) > 0) {
out_log_t <- data.frame(group = c(rep("pars", length(pars)), rep("output", length(outdat_log)), "meta"),
variable = c(names(pars), names(outdat_log), "time_total"),
value = c(round(pars, 4), round(unlist(outdat_log),4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F)
} else {
out_log_t <- data.frame(group = c(rep("pars", length(pars)), "meta"),
variable = c(names(pars), "time_total"),
value = c(round(pars, 4), round(difftime(time_end, time_start, units = "s"), 1)),
stringsAsFactors = F)
}
out_log <- out_log_t %>%
mutate(value = as.character(value)) %>%
bind_rows(.,dat_log)
write.table(out_log, file=logfile, sep="\t", quote=F, row.names=F, col.names=T)
}
# catch non-finite output
if(any(!is.finite(unlist(out)))) {
warning("Non-finite outut detected! Return NA. Consider argument 'log_meta'.")
return(NA)
}
# output
return(out)
} # EOF
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.