Input data
Precipitation (10-minute resolution) and evapotranspiration (daily
ET0) for 2011-2025 are provided by GeoSphere Austria (Österreichischer
Wetterdienst, formerly ZAMG) and ship with the package under
inst/extdata/models/wien/. The HDF5 model template
(base.h5) is produced with the Tandler
“Regenwasserbewirtschaftung” calculation engine; the engine itself is
downloaded from the KWB-R/kwb.raindrop.binaries GitHub
Release on demand via kwb.raindrop::download_engine().
Define Paths and Scenarios
library(kwb.raindrop)
path_list <- list(
modelname = "Wien",
root_path = file.path(tempdir(), "raindrop_wien"),
dir_input = "<root_path>/models/<modelname>/input",
dir_output = "<root_path>/models/<modelname>/output",
dir_target_output = "<dir_output>/<dir_target>",
file_errors_hdf5 = "Fehlerprotokoll.h5",
file_results_hdf5_element = "Mulde_Rigole.h5",
file_results_hdf5_flaeche = "Dach.h5",
file_results_hdf5_verschaltungen = "<dir_target>_Verschaltungen.h5",
file_results_txt = "Mulde_Rigole_RAINDROP.txt",
file_results_txt_multilayer = "Mulde_Rigole_RAINDROP_multi_layer.txt",
file_target = "<dir_target>.h5",
# Inputs (data source: GeoSphere Austria, period 2011-2025)
path_base = system.file("extdata/models/wien/base.h5", package = "kwb.raindrop"),
path_exe = if (is_windows && !is_ghactions) kwb.raindrop::download_engine() else NA_character_,
path_et = system.file("extdata/models/wien/et.csv", package = "kwb.raindrop"),
path_rain = system.file("extdata/models/wien/rain.csv.gz", package = "kwb.raindrop"),
# Output paths (still templated against tempdir scratch root)
path_errors_hdf5 = "<dir_target_output>/<file_errors_hdf5>",
path_results_hdf5_element = "<dir_target_output>/<file_results_hdf5_element>",
path_results_hdf5_flaeche = "<dir_target_output>/<file_results_hdf5_flaeche>",
path_results_hdf5_verschaltungen = "<dir_target_output>/<file_results_hdf5_verschaltungen>",
path_results_txt = "<dir_target_output>/<file_results_txt>",
path_results_txt_multilayer = "<dir_target_output>/<file_results_txt_multilayer>",
path_target_input = "<dir_input>/<file_target>"
)
parameters <- tibble::tibble(
para_nama_short = c(
"connected_area",
"mulde_area",
"mulde_height",
"filter_hydraulicconductivity",
# "filter_height",
"storage_height",
# "bottom_hydraulicconductivity",
"lai"
),
para_name_long = c(
"/Massnahmenelemente/Dach/Allgemein/Flaeche",
"/Massnahmenelemente/Mulde_Rigole/Allgemein/Flaeche",
"/Massnahmenelemente/Mulde_Rigole/Eigenschaften_Oberflaeche/Ueberlaufhoehe",
"Bodenarten/Bodenfilter/Ks_HydraulicConductivity",
#"/Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken",
"/Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken",
# "/Massnahmenelemente/Mulde_Rigole/Allgemein/Endversickerungsrate",
"/Massnahmenelemente/Mulde_Rigole/Parameter_Evapotranspiration/LAI_LeafAreaIndex"
),
index = c(1L,
1L,
1L,
1L,
2L,
1L)
)
DT::datatable(parameters,
filter = "top",
options = list(pageLength = 25,
autoWidth = TRUE))
connected_area <- 1000
mulde_area <- c(25, 50, 75, 100, 125, 150, 175, 200)
mulde_height <- c(100, 200, 300)
filter_hydraulicconductivity <- c(36, 180, 360)
filter_height <- 300
# Storage-layer (Speicher = 2nd Bodenschichtung layer) presets per storage type.
# Each type brings its own storage_height levels plus the Speicher soil
# parameters (theta*). They are written to the model in the run loop below:
# storage_height -> Schichtdicken[2]
# Startwerte_theta_ActualSoilMoisture -> Startwerte_theta_ActualSoilMoisture[2]
# thetaWP/thetaFC/thetaS -> //Bodenarten/Speicher/theta*
storage_type <- list(
"infiltration_box" = list(
storage_height = c(300, 600, 900, 1200),
Startwerte_theta_ActualSoilMoisture = 0,
thetaWP_MoistureAtWiltingPoint = 0,
thetaFC_MoistureAtFieldCapacity = 0,
thetaS_MoistureAtSaturation = 0.95
),
"gravel_trench" = list(
storage_height = 3 * c(300, 600, 900, 1200),
Startwerte_theta_ActualSoilMoisture = 0,
thetaWP_MoistureAtWiltingPoint = 0,
thetaFC_MoistureAtFieldCapacity = 0,
thetaS_MoistureAtSaturation = 0.3
)
)
rain_factor <- 1
bottom_hydraulicconductivity <- 12 #c(1,5,10,20,45,90,180,270,360,1860,3600)
# LAI for Mulde_Rigole only (Dach kept at H5 default).
# 8.5 = status-quo H5 default; 3.9 = grass per Hoernschemeyer et al.,
# Water 2023, 15, 2840, Tab. 6, plant type 5 (grasses/herbs).
#lai <- c(3.9, 8.5)
lai <- 3.9
# storage_height is coupled to storage_type (each type has its own levels), so
# build one row per (storage_type, storage_height) carrying the matching
# Speicher soil parameters, then cross-join with all remaining combinations.
storage_grid <- do.call(rbind, lapply(names(storage_type), function(type_name) {
spec <- storage_type[[type_name]]
data.frame(
storage_type = type_name,
storage_height = spec$storage_height,
storage_theta_start = spec$Startwerte_theta_ActualSoilMoisture,
storage_thetaWP = spec$thetaWP_MoistureAtWiltingPoint,
storage_thetaFC = spec$thetaFC_MoistureAtFieldCapacity,
storage_thetaS = spec$thetaS_MoistureAtSaturation,
stringsAsFactors = FALSE
)
}))
# Alle Kombinationen erzeugen
param_grid_all_combinations <- expand.grid(
connected_area = connected_area,
mulde_area = mulde_area,
mulde_height = mulde_height,
filter_hydraulicconductivity = filter_hydraulicconductivity,
filter_height = filter_height,
bottom_hydraulicconductivity = bottom_hydraulicconductivity,
rain_factor = rain_factor,
lai = lai,
stringsAsFactors = FALSE
)
# Cross-join the free parameters with the coupled storage grid.
param_grid_all_combinations <- merge(param_grid_all_combinations, storage_grid,
by = NULL)
param_grid_all_combinations <- param_grid_all_combinations %>%
dplyr::bind_cols(tibble::tibble(scenario_name = sprintf("s%05d",
seq_len(nrow(param_grid_all_combinations)))))
# Reference = first storage type at its smallest storage_height (storage_height
# is coupled to storage_type, so both are fixed together).
ref_storage_type <- names(storage_type)[1]
ref_storage_height <- min(storage_type[[ref_storage_type]]$storage_height)
ref_scenario <- param_grid_all_combinations %>%
dplyr::filter(connected_area == min(unique(param_grid_all_combinations$connected_area)),
mulde_area == min(unique(param_grid_all_combinations$mulde_area)),
filter_height == min(filter_height),
filter_hydraulicconductivity == min(param_grid_all_combinations$filter_hydraulicconductivity),
bottom_hydraulicconductivity == min(unique(param_grid_all_combinations$bottom_hydraulicconductivity)),
mulde_height == min(param_grid_all_combinations$mulde_height),
storage_type == ref_storage_type,
storage_height == ref_storage_height,
lai == max(param_grid_all_combinations$lai)) %>%
dplyr::pull(scenario_name)
stopifnot(length(ref_scenario)==1)
scenarios_with_single_parameter_variation <- kwb.raindrop::find_single_param_variations(
data = param_grid_all_combinations,
ref_scenario = ref_scenario
) %>%
dplyr::pull(scenario_name) %>% unique()
#> Rows with exactly one differing parameter: 14 of 576
#> Single-parameter variations per parameter: connected_area=0, mulde_area=7, mulde_height=2, filter_hydraulicconductivity=2, filter_height=0, bottom_hydraulicconductivity=0, rain_factor=0, lai=0, storage_type=0, storage_height=3, storage_theta_start=0, storage_thetaWP=0, storage_thetaFC=0, storage_thetaS=0
param_grid <- param_grid_all_combinations %>%
dplyr::filter(scenario_name %in% scenarios_with_single_parameter_variation)
param_grid <- param_grid_all_combinations
# Nutzbares Speichervolumen der Speicherschicht [m3] = Muldenflaeche x
# Speicherhoehe x nutzbare Porositaet (thetaS - thetaFC) des Speichertyps
# (Sickerbox 0.95, Schotterrigole 0.3). thetaFC statt thetaWP: nur das
# oberhalb der Feldkapazitaet entwaesserbare Porenvolumen leert sich zwischen
# den Ereignissen und steht als Retentionsvolumen erneut zur Verfuegung;
# Wasser zwischen WP und FC haelt die Schicht gegen die Schwerkraft (in den
# Presets sind thetaFC = thetaWP = 0, beide Definitionen also identisch).
# Erscheint in der Grid-Tabelle, den Ergebnis-CSVs und im Plot-Tooltip.
param_grid <- param_grid %>%
dplyr::mutate(
storage_volume_m3 = mulde_area * storage_height / 1000 *
(storage_thetaS - storage_thetaFC)
)
DT::datatable(param_grid,
filter = "top",
options = list(pageLength = 25,
autoWidth = TRUE))
htmlwidgets::saveWidget(DT::datatable(parameters,
filter = "top",
options = list(pageLength = 25,
autoWidth = TRUE)), "parameters.html")
htmlwidgets::saveWidget(DT::datatable(param_grid,
filter = "top",
options = list(pageLength = 25,
autoWidth = TRUE)), "param_grid.html")
psi_s_mm <- function(kf_mmh) (3.237 * (kf_mmh/25.4)^(-0.328)) * 25.4
paths <- kwb.utils::resolve(path_list)
timeseries_et <- readr::read_delim(paths$path_et, delim = ";", col_types = "cd") %>%
dplyr::mutate(date = lubridate::dmy(date),
time = difftime(date, min(date), units = "hours") %>% as.integer()) %>%
dplyr::select(-date) %>%
dplyr::filter(!is.na(value)) %>%
dplyr::relocate(time,.before = value)
timeseries_et$time[nrow(timeseries_et)] <- ceiling(timeseries_et$time[nrow(timeseries_et)])
timeseries_rain <- readr::read_csv(paths$path_rain) %>%
dplyr::rename(datetime = time,
value = rr) %>%
dplyr::mutate(time = difftime(datetime, min(datetime), units = "secs") %>% as.double()/3600) %>%
dplyr::filter(!is.na(value)) %>%
dplyr::select(-datetime, -station) %>%
dplyr::relocate(time,.before = value)
#> Rows: 788833 Columns: 3
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> dbl (2): station, rr
#> dttm (1): time
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
timeseries_rain$time[nrow(timeseries_rain)] <- ceiling(timeseries_rain$time[nrow(timeseries_rain)])
timeseries_et <- if(max(timeseries_rain$time) > max(timeseries_et$time)) {
rain_last_time <- max(timeseries_rain$time)
et_last_value <- timeseries_et$value[nrow(timeseries_et)]
msg <- sprintf("Rain time series is %f hours longer than ET time series. Adding final time %f with last value ('%f') in ET time series",
rain_last_time - max(timeseries_et$time),
rain_last_time,
et_last_value
)
kwb.utils::catAndRun(messageText = msg,
expr = {
timeseries_et %>%
dplyr::bind_rows(tibble::tibble(time = rain_last_time,
value = et_last_value))
})
} else {
timeseries_et
}
timeseries_rain <- if(max(timeseries_et$time) > max(timeseries_rain$time)) {
et_last_time <- max(timeseries_rain$time)
rain_last_value <- timeseries_rain$value[nrow(timeseries_rain)]
msg <- sprintf("ET time series is %f hours longer than rain time series. Adding final time %f with last value ('%f') in rain time series",
et_last_time - max(timeseries_rain$time),
et_last_time,
rain_last_value
)
kwb.utils::catAndRun(messageText = msg,
expr = {
timeseries_rain %>%
dplyr::bind_rows(tibble::tibble(time = et_last_time,
value = rain_last_value))
})
} else {
timeseries_rain
}
txt <- sprintf("F\u00fcr den Datensatz '%s' gibt es %.2f %% NA Werte und %.2f %% Werte die gleich Null sind. (Regenmenge: %f mm/a)\n",
paths$path_rain,
100*sum(is.na(timeseries_rain$value))/nrow(timeseries_rain),
100*sum(timeseries_rain$value == 0, na.rm = TRUE)/nrow(timeseries_rain),
sum(timeseries_rain$value, na.rm = TRUE)/((range(timeseries_rain$time)[2] - range(timeseries_rain$time)[1])/24/365))
message(txt)
#> Für den Datensatz 'D:/a/_temp/Library/kwb.raindrop/extdata/models/wien/rain.csv.gz' gibt es 0.00 % NA Werte und 96.04 % Werte die gleich Null sind. (Regenmenge: 647.465225 mm/a)
txt <- sprintf("F\u00fcr den Datensatz '%s' gibt es %.2f %% NA Werte und %.2f %% Werte die gleich Null sind. (Evapotranspirations: %f mm/a)\n",
paths$path_et,
100*sum(is.na(timeseries_et$value))/nrow(timeseries_et),
100*sum(timeseries_et$value == 0, na.rm = TRUE)/nrow(timeseries_et),
sum(timeseries_et$value, na.rm = TRUE)/((range(timeseries_rain$time)[2] - range(timeseries_rain$time)[1])/24/365))
message(txt)
#> Für den Datensatz 'D:/a/_temp/Library/kwb.raindrop/extdata/models/wien/et.csv' gibt es 0.00 % NA Werte und 0.00 % Werte die gleich Null sind. (Evapotranspirations: 888.173330 mm/a)
### Convert rain from mm to mm/h
period <- c(diff(timeseries_rain$time), mean(diff(timeseries_rain$time)))
timeseries_rain$value <- timeseries_rain$value / period
### Convert ET0 from mm/day to mm/h (the engine reads //Kurven/ET0 as a mm/h
### rate, exactly like rain; daily values must be divided by their interval in
### hours = 24, otherwise ET0 is integrated 24x too high)
period_et <- c(diff(timeseries_et$time), mean(diff(timeseries_et$time)))
timeseries_et$value <- timeseries_et$value / period_et
#openxlsx::write.xlsx(list(regen = timeseries_rain, et = timeseries_et), "timeseries.xlsx")
run_one <- function(i,
timestep_hours,
debug = FALSE,
...) {
param_grid_tmp <- param_grid[i, ]
paths <- kwb.utils::resolve(path_list,
dir_target = param_grid_tmp$scenario_name)
fs::dir_create(paths$dir_input, recurse = TRUE)
fs::dir_create(paths$dir_output, recurse = TRUE)
fs::dir_create(paths$dir_target_output, recurse = TRUE)
fs::file_copy(path = paths$path_base,
new_path = paths$path_target_input,
overwrite = TRUE)
h5 <- hdf5r::H5File$new(paths$path_target_input, mode = "a")
new_path <- stringr::str_c(normalizePath(fs::path_abs(paths$dir_target_output)),
"\\")
vals <- kwb.raindrop::h5_read_values(h5)
vals$`//Berechnungsparameter/Ergebnispfad` <- new_path
vals$`//Berechnungsparameter/Zeitschritt_Infiltration` <- timestep_hours
vals$`//Berechnungsparameter/Zeitschritt_ET` <- timestep_hours
vals$`//Berechnungsparameter/Zeitschritt_Verschaltungen` <- timestep_hours
vals$`//Berechnungsparameter/R-Plots` <- 0
vals$`//Berechnungsparameter/Ausgabemodus` <- "Optimierung"
vals$`//Berechnungsparameter/Evapotranspiration_aktiv` <- 1
vals$`//Massnahmenelemente/Dach/Berechnungsparameter/Evapotranspiration_aktiv` <- 1
vals$`//Massnahmenelemente/Dach/Allgemein/Flaeche` <- param_grid_tmp$connected_area
vals$`//Massnahmenelemente/Mulde_Rigole/Berechnungsparameter/Evapotranspiration_aktiv` <- 1
vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Regen-Skalierungsfaktor` <- 1
vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Flaeche` <- param_grid_tmp$mulde_area
vals$`//Massnahmenelemente/Mulde_Rigole/Eigenschaften_Oberflaeche/Ueberlaufhoehe` <- param_grid_tmp$mulde_height
vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Startwerte_theta_ActualSoilMoisture` <- c(0.3, param_grid_tmp$storage_theta_start)
vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken` <- c(param_grid_tmp$filter_height,
param_grid_tmp$storage_height)
# Speicher (2nd Bodenschichtung layer) soil parameters depend on the
# storage type (infiltration_box vs. gravel_trench); see `storage_type`.
vals$`//Bodenarten/Speicher/thetaWP_MoistureAtWiltingPoint` <- param_grid_tmp$storage_thetaWP
vals$`//Bodenarten/Speicher/thetaFC_MoistureAtFieldCapacity` <- param_grid_tmp$storage_thetaFC
vals$`//Bodenarten/Speicher/thetaS_MoistureAtSaturation` <- param_grid_tmp$storage_thetaS
vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Endversickerungsrate` <- param_grid_tmp$bottom_hydraulicconductivity
vals$`//Bodenarten/Bodenfilter/Ks_HydraulicConductivity` <- param_grid_tmp$filter_hydraulicconductivity
vals$`//Bodenarten/Bodenfilter/Psi_Saugspannung_CapillarySuction` <- psi_s_mm(param_grid_tmp$filter_hydraulicconductivity)
vals$`//Massnahmenelemente/Mulde_Rigole/Parameter_Evapotranspiration/LAI_LeafAreaIndex` <- param_grid_tmp$lai
vals$`//Kurven/ET0` <- timeseries_et
vals$`//Kurven/Regen` <- timeseries_rain
vals$`//Kurven/Growth_1`$time[2] <- max(timeseries_rain$time)
vals$`//Kurven/Shading_1`$time[2] <- max(timeseries_rain$time)
kwb.raindrop::h5_write_values(h5, vals, resize = TRUE,
scalar_strategy = "error",
verbose = FALSE)
h5$close_all()
kwb.raindrop::run_model(path_exe = paths$path_exe,
path_input = paths$path_target_input,
debug = debug)
# Thin immediately: read only this run's results (lean = water balance +
# overflow rates, no states/meta/connected-area rates) and reduce to the
# single optimisation row. This way we never hold all scenarios' full
# time series in memory at once; the full result HDF5 stays on disk for
# ad-hoc inspection.
sim_one <- kwb.raindrop::get_simulation_results_optim(
paths = paths,
path_list = path_list,
simulation_names = param_grid_tmp$scenario_name,
debug = debug,
lean = TRUE
)
kwb.raindrop::add_overflow_events_and_waterbalance(
simulation_results = sim_one,
event_separation_hours = 4,
canonical_variables = kwb.raindrop::default_canonical_wb_variables()
)
}
n_cores <- parallel::detectCores()
# run_one() now returns the thinned per-run optimisation row, so run_scenarios()
# yields a list of one-row tibbles we simply bind below.
scenario_rows <- NULL
system.time(expr = {
scenario_rows <- kwb.raindrop::run_scenarios(indices = seq_len(nrow(param_grid)),
run_one_scenario = run_one,
timestep_hours = 0.1,
debug = FALSE,
parallel = TRUE,
workers = n_cores,
show_progress = TRUE
)
}
)
### Read results for first run
if(FALSE) {
paths <- kwb.utils::resolve(path_list,
dir_target = sprintf("s%05d", i = 1))
#simulation_names <- basename(fs::dir_ls(paths$dir_output))
simulation_names <- param_grid$scenario_name
debug <- TRUEn_cores
errors_df <- kwb.raindrop::read_raindrop_errors(simulation_names, path_list)
x <- tidyr::unnest(errors_df, errors)
x$Fehlerbeschreibung
}Analyse Results
# Each run was already thinned to its optimisation row inside run_one(), so we
# just bind the per-run rows here instead of re-reading every run's full
# results into memory. (The previous get_simulation_results_optim_parallel() +
# add_overflow_events_and_waterbalance() pass loaded all runs at once.)
simulation_results_optimisation <- dplyr::bind_rows(scenario_rows)
simulation_results_optimisation <- param_grid %>%
dplyr::left_join(simulation_results_optimisation,
by = c("scenario_name" = "s_name")) %>%
dplyr::relocate(scenario_name, .before = connected_area) %>%
kwb.raindrop::compute_costs()
readr::write_csv(simulation_results_optimisation,
file = sprintf("simulation_results_optimisation_%s.csv",
paths$modelname)
)
htmlwidgets::saveWidget(DT::datatable(simulation_results_optimisation,
filter = "top",
options = list(pageLength = 25,
autoWidth = TRUE)),
file = sprintf("simulation_results_optimisation_%s.html",
paths$modelname),
title = sprintf("RAINDROP - '%s': Solution Space",
paths$modelname)
)
### Plot results
# Fuer Plots/Tooltips: die an storage_type gekoppelten Speicher-Bodenparameter
# (storage_theta*) sind durch den Typ bestimmt, also redundant - sie wuerden
# nur jeden Tooltip aufblaehen. Fuer den Modelllauf oben werden sie gebraucht,
# ab hier nicht mehr. storage_volume_m3 bekommt im Tooltip eine eigene Zeile
# (aus den Ergebnisdaten) und fliegt hier ebenfalls raus, sonst stuende es
# doppelt unter "Variierende Parameter".
param_grid <- param_grid %>%
dplyr::select(-dplyr::starts_with("storage_theta"),
-dplyr::any_of("storage_volume_m3"))
params <- c(
#"connected_area",
"mulde_area",
"mulde_height",
"filter_hydraulicconductivity",
#"filter_height",
"storage_height",
#"bottom_hydraulicconductivity",
#"rain_factor",
"lai",
"storage_type"
)
lang <- "de"
max_n_overflows <- 5
# Kostensaetze-Caption fuer die interaktiven Kostenplots: ggplotly verwirft
# ggplot-Captions, daher wird sie dort per plotly_add_caption() nachgeruestet;
# die PDFs bekommen sie automatisch ueber den caption-Default der Funktionen.
cost_caption <- kwb.raindrop::cost_rates_caption(lang)
pdff <- sprintf("simulation_results_optimisation_%s_main-effects.pdf",
paths$modelname)
gg <- kwb.raindrop::plot_main_effects(
df = simulation_results_optimisation,
y = "n_overflows",
params = params,
lang = lang,
ylim_lower = 0
)
grDevices::pdf(pdff, width = 9, height = 3, onefile = TRUE)
print(gg)
dev.off()
#kwb.utils::finishAndShowPdf(pdff)
pdff <- sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_mulde-height_single.pdf",
paths$modelname)
# kwb.utils::preparePdf(pdfFile = pdff)
# on.exit(dev.off(), add = TRUE)
grDevices::pdf(pdff, width = 9, height = 4, onefile = TRUE)
p <- kwb.raindrop::plot_valid_design_space(
param_grid = param_grid,
sim_results = simulation_results_optimisation,
x = "mulde_area",
y = "mulde_height",
valid_max = max_n_overflows,
jitter = TRUE,
alpha_mode = "duplicates",
alpha_min = 0.25,
alpha_max = 1,
drop_overflow_gt_valid_max = FALSE,
keep_param_grid_limits = TRUE,
facet_storage_type = TRUE,
lang = lang,
subtitle = ""
)
suppressWarnings(print(p))
dev.off()
# --- 2) Interaktiv als HTML: nach dem PDF
plotly_gg <- plotly::ggplotly(gg)
htmlwidgets::saveWidget(
widget = plotly_gg,
file = sprintf("simulation_results_optimisation_%s_main-effects.html",
paths$modelname),
selfcontained = TRUE,
title = sprintf("RAINDROP - '%s': Main Effects",
paths$modelname)
)
pdff <- sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_parameters.pdf",
paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)
for (y in c("mulde_height", "filter_hydraulicconductivity", "storage_height")) {
p <- kwb.raindrop::plot_valid_design_space(
param_grid = param_grid,
sim_results = simulation_results_optimisation,
x = "mulde_area",
y = y,
valid_max = max_n_overflows,
jitter = TRUE,
alpha_mode = "duplicates",
alpha_min = 0.25,
alpha_max = 1,
drop_overflow_gt_valid_max = FALSE,
keep_param_grid_limits = TRUE,
facet_storage_type = TRUE
)
# interaktiv als HTML; Farb-Legende je Ueberlaufklasse statt
# (Farbe, Form)-Tupeln - die Formen erklaeren die Panel-Beschriftungen.
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_split_legend(plotly_p, lang = lang,
add_shape_legend = FALSE)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_%s.html",
paths$modelname,
y),
selfcontained = TRUE,
title = sprintf("'%s' - Design Space: mulde_area vs. %s",
paths$modelname,
y)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
}
dev.off()
#kwb.utils::finishAndShowPdf(pdff)
pdff <- sprintf("simulation_results_optimisation_%s_water-balance.pdf",
paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)
p <- kwb.raindrop::plot_wb_tradeoff_overflows(
simulation_results_optimisation = simulation_results_optimisation,
param_grid = param_grid,
x = max_n_overflows,
filter_n_gtx = FALSE,
use_jitter = TRUE
)
# interaktiv als HTML; getrennte Legenden (Ueberlaufklassen + Speichertyp)
# statt der (Farbe, Form)-Tupel von ggplotly
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_split_legend(plotly_p, lang = lang)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf("simulation_results_optimisation_%s_water-balance.html",
paths$modelname),
selfcontained = TRUE,
title = sprintf("'%s' - Water balance vs overflows",
paths$modelname)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
dev.off()
#kwb.utils::finishAndShowPdf(pdff)
pdff <- sprintf("simulation_results_optimisation_%s_cost-vs-overflow-volume.pdf",
paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)
p <- kwb.raindrop::plot_cost_vs_overflow_volume(
simulation_results_optimisation = simulation_results_optimisation,
param_grid = param_grid,
x = max_n_overflows,
filter_n_gtx = FALSE,
use_jitter = TRUE
)
# interaktiv als HTML; getrennte Legenden (Ueberlaufklassen + Speichertyp)
# statt der (Farbe, Form)-Tupel von ggplotly
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_split_legend(plotly_p, lang = lang)
plotly_p <- kwb.raindrop::plotly_add_caption(plotly_p, cost_caption)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf("simulation_results_optimisation_%s_cost-vs-overflow-volume.html",
paths$modelname),
selfcontained = TRUE,
title = sprintf("'%s' - Cost vs. overflow volume",
paths$modelname)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
dev.off()
#kwb.utils::finishAndShowPdf(pdff)
# Drei Kosten-Boxplots mit unterschiedlichem Optimierungsziel je Kategorie
# (Kosten als Tie-Break): (i) guenstigste; (ii) geringstes Ueberlaufvolumen
# (Label m3 + %); (iii) hoechste Evapotranspiration (Label %, Punktgroesse =
# Evapotranspiration). x = max_n_overflows, Linie ueber alle Klassen. Die beiden
# Speichertypen liegen als zwei Panels untereinander (Sickerbox oben,
# Schotterrigol unten); die Punkte bleiben Kreise, da die Panel-Streifen
# den Typ bereits benennen.
cost_boxplots <- list(
list(suffix = "cheapest", best_by = "min_cost",
size_by = "overflow_volume", label_best = FALSE),
list(suffix = "min-overflow", best_by = "min_overflow",
size_by = "overflow_volume", label_best = TRUE),
list(suffix = "max-evap", best_by = "max_evapotranspiration",
size_by = "evapotranspiration", label_best = TRUE)
)
for (cb in cost_boxplots) {
pdff <- sprintf(
"simulation_results_optimisation_%s_cost-by-overflows-boxplot-%s.pdf",
paths$modelname, cb$suffix)
kwb.utils::preparePdf(pdfFile = pdff)
p <- kwb.raindrop::plot_cost_overflow_boxplot(
simulation_results_optimisation = simulation_results_optimisation,
param_grid = param_grid,
x = max_n_overflows,
filter_n_gtx = FALSE,
use_jitter = TRUE,
lang = lang,
facet_storage_type = TRUE,
size_by = cb$size_by,
best_by = cb$best_by,
label_best = cb$label_best
)
# interaktiv als HTML
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_add_caption(plotly_p, cost_caption)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf(
"simulation_results_optimisation_%s_cost-by-overflows-boxplot-%s.html",
paths$modelname, cb$suffix),
selfcontained = TRUE,
title = sprintf("'%s' - Cost boxplot (%s)", paths$modelname, cb$suffix)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
dev.off()
}
# Kosten vs. Evapotranspiration: Streudiagramm ueber den Design-Raum, Punktform
# kodiert den Speichertyp (Viereck = Sickerbox, Dreieck = Schotterrigol),
# Farbe die Anzahl Ueberlaufereignisse.
pdff <- sprintf("simulation_results_optimisation_%s_cost-vs-evaporation.pdf",
paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)
p <- kwb.raindrop::plot_cost_vs_evaporation(
simulation_results_optimisation = simulation_results_optimisation,
param_grid = param_grid,
x = max_n_overflows,
filter_n_gtx = FALSE,
use_jitter = TRUE,
lang = lang
)
# interaktiv als HTML; getrennte Legenden (Ueberlaufklassen + Speichertyp)
# statt der (Farbe, Form)-Tupel von ggplotly
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_split_legend(plotly_p, lang = lang)
plotly_p <- kwb.raindrop::plotly_add_caption(plotly_p, cost_caption)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf("simulation_results_optimisation_%s_cost-vs-evaporation.html",
paths$modelname),
selfcontained = TRUE,
title = sprintf("'%s' - Cost vs. evapotranspiration",
paths$modelname)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
dev.off()
# Kosten je Prozent Evapotranspiration [EUR/%]: Boxplot je Ueberlaufklasse mit den
# beiden Speichertypen als zwei Panels untereinander (Sickerbox oben,
# Schotterrigol unten); guenstigstes Szenario je Box markiert (Label EUR/%),
# Punktgroesse = Evapotranspiration.
pdff <- sprintf(
"simulation_results_optimisation_%s_cost-per-evap-boxplot.pdf",
paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)
p <- kwb.raindrop::plot_cost_overflow_boxplot(
simulation_results_optimisation = simulation_results_optimisation,
param_grid = param_grid,
x = max_n_overflows,
filter_n_gtx = FALSE,
use_jitter = TRUE,
lang = lang,
y_var = "cost_per_evap_pct",
facet_storage_type = TRUE,
size_by = "evapotranspiration",
best_by = "min_cost",
label_best = TRUE
)
# interaktiv als HTML
plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
plotly_p <- kwb.raindrop::plotly_add_caption(plotly_p, cost_caption)
htmlwidgets::saveWidget(
widget = plotly_p,
file = sprintf(
"simulation_results_optimisation_%s_cost-per-evap-boxplot.html",
paths$modelname),
selfcontained = TRUE,
title = sprintf("'%s' - Cost per %% evapotranspiration", paths$modelname)
)
# statisch ins PDF (WICHTIG!)
suppressWarnings(print(p))
dev.off()