diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 04c4dec..3d68bbb 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -1,169 +1,167 @@ # Preface ---------------------------------------------------------------------- -remotes::install_github("kwb-r/kwb.hydrus1d@dev") -remotes::install_github("kwb-r/flextreat.hydrus1d@dev") +if (FALSE) { + remotes::install_github("kwb-r/kwb.hydrus1d@dev") + remotes::install_github("kwb-r/flextreat.hydrus1d@dev") +} + library(magrittr) library(flextreat.hydrus1d) - -create_temp_dir <- function(base_dir = NULL) { - # Erstellen des Basisverzeichnisses, falls nicht vorhanden - if (is.null(base_dir)) { - base_dir <- tempdir() # Erzeugt ein temporäres Verzeichnis - dir.create(base_dir, showWarnings = FALSE) +# get_atm ---------------------------------------------------------------------- +get_atm <- function(atm, extreme_rain = NULL) +{ + if (is.null(extreme_rain)) { + return(atm) } - # Generieren eines zufälligen Namens für das Unterverzeichnis - sub_dir <- file.path(base_dir, paste0("hydrus1d_", paste0(sample(letters, 10, replace = TRUE), collapse = ""))) - - # Erstellen des Unterverzeichnisses - dir.create(sub_dir) - - return(sub_dir) -} - -# Beispielaufrufe der Funktion -base_dir <- create_temp_dir() # Erstellt das Basisverzeichnis -cat("Basisverzeichnis:", base_dir, "\n") + if (!extreme_rain %in% c("dry", "wet")) { + stop("extreme_rain has to be either 'NULL', 'dry' or 'wet") + } -# Erstellen von zufälligen Unterverzeichnissen -sub_dir1 <- create_temp_dir(base_dir) -cat("Unterverzeichnis 1:", sub_dir1, "\n") + atm <- atm %>% + dplyr::mutate( + hydrologic_year = flextreat.hydrus1d::get_hydrologic_years(date), + year = as.integer(format(date, format = "%Y")), + day_of_year = lubridate::yday(date) + ) -sub_dir2 <- create_temp_dir(base_dir) -cat("Unterverzeichnis 2:", sub_dir2, "\n") + atm_stats <- atm %>% + dplyr::group_by(year) %>% + dplyr::summarise( + rain_mm = sum(rain_mm, na.rm = TRUE), + evapo_p_mean_mm = sum(evapo_p_mean_mm, na.rm = TRUE) + ) + if (extreme_rain == "dry") { + atm_dry <- atm_stats %>% + dplyr::filter(rain_mm == min(rain_mm)) -random_wait <- function() { - wait_time <- runif(1, min = 1, max = 10) # Generiert eine zufällige Zahl zwischen 1 und 10 - cat("Wartezeit:", wait_time, "Sekunden\n") # Ausgabe der Wartezeit zur Information - Sys.sleep(wait_time) # Warten für die zufällige Anzahl von Sekunden -} + atm_dry_sel <- atm %>% + dplyr::filter(year == atm_dry$year) %>% + dplyr::select(day_of_year, rain_mm) -get_atm <- function(atm, extreme_rain = NULL) { + atm_dry_input <- atm %>% + dplyr::select(- rain_mm, - hydrologic_year, - year) %>% + dplyr::left_join(atm_dry_sel) %>% + dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>% + dplyr::select(- day_of_year) %>% + dplyr::relocate(rain_mm, .after = clearwater.mmPerDay) - if(is.null(extreme_rain)) { - return(atm) + return(atm_dry_input) } - if (extreme_rain %in% c("dry", "wet")) { - atm <- atm %>% - dplyr::mutate(hydrologic_year = flextreat.hydrus1d::get_hydrologic_years(date), - year = as.integer(format(date, format = "%Y")), - day_of_year = lubridate::yday(date)) - - atm_stats <- atm %>% - dplyr::group_by(year) %>% - dplyr::summarise(rain_mm = sum(rain_mm, na.rm = TRUE), - evapo_p_mean_mm = sum(evapo_p_mean_mm, na.rm = TRUE)) - - if(extreme_rain == "dry") { - atm_dry <- atm_stats %>% - dplyr::filter(rain_mm == min(rain_mm)) - - atm_dry_sel <- atm %>% - dplyr::filter(year == atm_dry$year) %>% - dplyr::select(day_of_year, rain_mm) - - atm_dry_input <- atm %>% - dplyr::select(- rain_mm, - hydrologic_year, - year) %>% - dplyr::left_join(atm_dry_sel) %>% - dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>% - dplyr::select(- day_of_year) %>% - dplyr::relocate(rain_mm, .after = clearwater.mmPerDay) - - return(atm_dry_input) - } - - if (extreme_rain == "wet") { - atm_wet <- atm_stats %>% - dplyr::filter(rain_mm == max(rain_mm)) + if (extreme_rain == "wet") { + atm_wet <- atm_stats %>% + dplyr::filter(rain_mm == max(rain_mm)) + atm_wet_sel <- atm %>% + dplyr::filter(year == atm_wet$year) %>% + dplyr::select(day_of_year, rain_mm) - atm_wet_sel <- atm %>% - dplyr::filter(year == atm_wet$year) %>% - dplyr::select(day_of_year, rain_mm) + atm_wet_input <- atm %>% + dplyr::select(- rain_mm, - hydrologic_year, - year) %>% + dplyr::left_join(atm_wet_sel) %>% + dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>% + dplyr::select(- day_of_year) %>% + dplyr::relocate(rain_mm, .after = clearwater.mmPerDay) - atm_wet_input <- atm %>% - dplyr::select(- rain_mm, - hydrologic_year, - year) %>% - dplyr::left_join(atm_wet_sel) %>% - dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>% - dplyr::select(- day_of_year) %>% - dplyr::relocate(rain_mm, .after = clearwater.mmPerDay) - return(atm_wet_input) - } else { - stop("extreme_rain has to be either 'NULL', 'dry' or 'wet") - } - }} - - -prepare_solute_input <- function(dat, - selector, - Ks = NULL, - SnkL1 = NULL, - diff_w = 0, diff_g = 0, - kd = kd, - halftime_to_firstorderrate = halftime_to_firstorderrate) { + return(atm_wet_input) + } +} +# prepare_solute_input --------------------------------------------------------- +prepare_solute_input <- function( + dat, + selector, + Ks = NULL, + SnkL1 = NULL, + diff_w = 0, diff_g = 0, + kd = kd, + halftime_to_firstorderrate = halftime_to_firstorderrate +) +{ stopifnot(nrow(dat) <= 10) - if(selector$solute$No.Solutes != nrow(dat)) { + if (selector$solute$No.Solutes != nrow(dat)) { selector$solute$No.Solutes <- nrow(dat) } solute_names <- sprintf("solute_%d", seq_len(nrow(dat))) - solutes_new <- setNames(lapply(seq_len(nrow(dat)), function(i) { + solutes_new <- setNames(nm = solute_names, lapply(seq_len(nrow(dat)), function(i) { + dat_sel <- dat[i,] - ks <- kd(porosity = selector$waterflow$soil$ths - selector$waterflow$soil$thr, - retardation = dat_sel$retard, - bulk_density = selector$solute$transport$Bulk.d.) + ks <- kd( + porosity = selector$waterflow$soil$ths - selector$waterflow$soil$thr, + retardation = dat_sel$retard, + bulk_density = selector$solute$transport$Bulk.d. + ) - cols <- c("Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'", - "SnkS1'", "SnkG1'", "SnkL0", "SnkS0", "SnkG0", "Alfa") + cols <- c( + "Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'", + "SnkS1'", "SnkG1'", "SnkL0", "SnkS0", "SnkG0", "Alfa" + ) reaction <- matrix(data = 0, ncol = length(cols), nrow = length(ks)) %>% - as.data.frame() %>% tibble::as_tibble() + as.data.frame() %>% + tibble::as_tibble() + names(reaction) <- cols reaction$Beta <- 1 - reaction$Ks <- if(is.null(Ks)) {round(ks, 2)} else { Ks} - reaction$SnkL1 <- if(is.null(SnkL1)) { + + reaction$Ks <- if (is.null(Ks)) { + round(ks, 2) + } else { + Ks + } + + reaction$SnkL1 <- if (is.null(SnkL1)) { round(halftime_to_firstorderrate(dat_sel$half_life_days), 5) } else { SnkL1 } - list(diffusion = tibble::tibble(DifW = diff_w, DifG = diff_g), - reaction = reaction) - }), nm = solute_names) + list( + diffusion = tibble::tibble(DifW = diff_w, DifG = diff_g), + reaction = reaction + ) + + })) sel_tmp <- selector$solute[!names(selector$solute) %in% solute_names] - solutes_new_list <- c(sel_tmp[1:which(names(sel_tmp) == "transport")], - solutes_new, - sel_tmp[which(names(sel_tmp) == "kTopSolute"):length(sel_tmp)]) + solutes_new_list <- c( + sel_tmp[1:which(names(sel_tmp) == "transport")], + solutes_new, + sel_tmp[which(names(sel_tmp) == "kTopSolute"):length(sel_tmp)] + ) - c(selector[names(selector) != "solute"], - list(solute = solutes_new_list)) + c( + selector[names(selector) != "solute"], + list(solute = solutes_new_list) + ) } -get_mean <- function(col) { - x <- col %>% stringr::str_split_fixed("-", n = 2) +# get_mean --------------------------------------------------------------------- +get_mean <- function(col) +{ + x <- stringr::str_split_fixed(col, "-", n = 2) mode(x) <- "numeric" round(rowMeans(x), digits = 2) } - -kd <- function(porosity, retardation, bulk_density) { - +# kd --------------------------------------------------------------------------- +kd <- function(porosity, retardation, bulk_density) +{ #https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html (retardation - 1) * porosity / bulk_density - } -halftime_to_firstorderrate <- function(half_time) { - +# halftime_to_firstorderrate --------------------------------------------------- +halftime_to_firstorderrate <- function(half_time) +{ if(half_time != 0) { #https://chem.libretexts.org/Courses/Bellarmine_University/BU%3A_Chem_104_(Christianson)/Phase_2%3A_Understanding_Chemical_Reactions/4%3A_Kinetics%3A_How_Fast_Reactions_Go/4.5%3A_First_Order_Reaction_Half-Life#mjx-eqn-21.4.2 0.693 / half_time @@ -172,334 +170,382 @@ halftime_to_firstorderrate <- function(half_time) { } } +# provide_soil_columns --------------------------------------------------------- +provide_soil_columns <- function(path) +{ + kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE) %>% + janitor::clean_names() %>% + dplyr::mutate(half_life_days = dplyr::case_when( + grepl(">", hwz_tage) ~ hwz_tage %>% + stringr::str_remove(">") %>% + as.numeric() %>% + round(digits = 2), + grepl("<", hwz_tage) ~ hwz_tage %>% + stringr::str_remove("<") %>% + as.numeric() %>% + round(digits = 2), + grepl("-", hwz_tage) ~ get_mean(hwz_tage), + .default = as.numeric(hwz_tage) %>% + round(digits = 2)), + retard = dplyr::case_when( + grepl("-", retardation) ~ get_mean(retardation), + is.na(retardation) ~ 1L, + .default = as.numeric(retardation) %>% + round(digits = 2))) %>% + dplyr::filter(!is.na(half_life_days)) %>% + dplyr::mutate(id = 1:dplyr::n()) %>% + dplyr::relocate(id) # %>% dplyr::mutate(retard = 1, half_life_days = 0) +} + +# MAIN ------------------------------------------------------------------------- + #path <- "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_2_Boden-Grundwasser/daten_karten/Sickerwasserprognose/column-studies/Stoffeigenschaften_Säulen.xlsx" path <- "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/StofflicheModellrandbedingungen.xlsx" -soil_columns <- kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE) %>% - janitor::clean_names() %>% - dplyr::mutate(half_life_days = dplyr::case_when( - grepl(">", hwz_tage) ~ hwz_tage %>% stringr::str_remove(">") %>% as.numeric() %>% round(digits = 2), - grepl("<", hwz_tage) ~ hwz_tage %>% stringr::str_remove("<") %>% as.numeric() %>% round(digits = 2), - grepl("-", hwz_tage) ~ get_mean(hwz_tage), - .default = as.numeric(hwz_tage) %>% round(digits = 2)), - retard = dplyr::case_when( - grepl("-", retardation) ~ get_mean(retardation), - is.na(retardation) ~ 1L, - .default = as.numeric(retardation) %>% round(digits = 2))) %>% - dplyr::filter(!is.na(half_life_days)) %>% - dplyr::mutate(id = 1:dplyr::n()) %>% - dplyr::relocate(id) #%>% -#dplyr::mutate(retard = 1#, -#half_life_days = 0 -# ) - +soil_columns <- provide_soil_columns(path) ### Select 1 substance for 5 different half life classes defined in this table selected_substances <- readr::read_csv("inst/extdata/input-data/substance_classes.csv") -soil_columns_selected <- soil_columns %>% +soil_columns_selected <- soil_columns %>% dplyr::filter(substanz_nr %in% selected_substances$substance_id) %>% dplyr::mutate(id = 1:dplyr::n()) +arg_combis <- kwb.utils::expandGrid( + extreme_rain = c("", "wet", "dry"), # "" needs to be treated as NULL! + treatment = c("ka"), # "tracer" # c("ka", "o3") + scenario = unlist(lapply(c(1,10), function(x) { + paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x)) + })), + irrig_only_growing_season = FALSE, + duration_string = "long", #c("short", "long"), + retardation_scenario = c("retardation_yes", "retardation_no", "tracer") +) + +# generate_solute_ids ---------------------------------------------------------- +generate_solute_ids <- function(n) +{ + seq_start <- seq(1, n, 10) + + seq_end <- if (n > 10) { + seq(10, n, 10) + } else { + n + } -tracer <- FALSE -short <- FALSE -irrig_only_growing_season <- FALSE -retardation <- TRUE + if (length(seq_end) < length(seq_start)) { + seq_end <- c(seq_end, n) + } -soil_columns_selected <- if(retardation) { soil_columns_selected %>% - dplyr::mutate(retard = 1) -} else { - soil_columns_selected + tibble::tibble( + start = as.integer(seq_start), + end = as.integer(seq_end) + ) } -extreme_rains <- c("wet", "dry") # c(NULL, "wet", "dry") - -scenarios <- sapply(c(1,10), function(x) paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x))) %>% - as.vector() -#scenarios <- scenarios[2:7] +# generate_periods ------------------------------------------------------------- +generate_periods <- function(n) +{ + tibble::tibble( + start = seq(1, n, 10), + end = if (n %% 10 != 0) { + c(seq(10, n, 10), n) + } else { + seq(10, n, 10) + } + ) +} -#treatments <- "tracer" #c("ka", "o3") #"ka" #c("ka", "o3") #"ka" -treatments <- c("ka") +# prepare_files_for_irrig_int -------------------------------------------------- +prepare_files_for_irrig_int <- function(paths) +{ + p <- kwb.utils::createAccessor(paths) -str_final_dir <- if(treatments == "tracer") { "tracer"} else if (retardation) { - "retardation_yes"} else { - "retardation_no" + copy <- function(fun, from, to) { + kwb.utils::catAndRun( + sprintf("Copying from\n %s to\n %s", from, to), + fun(path = from, new_path = to) + ) } -# Main loop -------------------------------------------------------------------- -sapply(extreme_rains, function(extreme_rain) { + if (!fs::dir_exists(p("model_dir"))) { + copy(fs::dir_copy, from = p("model_dir_org"), to = p("model_dir")) + } - { + if (!fs::file_exists(p("model_gui_path"))) { + copy(fun = fs::file_copy, from = p("model_gui_path_org"), to = p("model_gui_path")) + } + model_out_files <- list.files( + p("model_dir"), + pattern = "\\.out$", + full.names = TRUE + ) + if (length(model_out_files) > 0L) { + fs::file_delete(model_out_files) } - sapply(treatments, function(treatment) { - - inner_function <- function(scenario, - treatment, - irrig_only_growing_season, - irrig_dir_string, - short, tracer, - duration_string, - extreme_rain, - extreme_rain_string, - str_final_dir, - soil_columns, - atm, - days_monthly, - get_atm, - kd, - halftime_to_firstorderrate, - periods, - prepare_solute_input, - check_hydrus_exe, - base_dir, - create_temp_dir, - `%>%`) { - - library(flextreat.hydrus1d) - - - atm <- get_atm(atm = atm, - extreme_rain = extreme_rain) - - if(irrig_only_growing_season) { - atm[which(!lubridate::month(atm$date) %in% 4:9), c("groundwater.mmPerDay", "clearwater.mmPerDay")] <- 0 - } - + soil_depth_cm <- 100 * + stringr::str_extract(p("model_name"), "soil-[0-9]+?m") %>% + stringr::str_extract("[0-9]") %>% + as.numeric() - atm <- if(short) { - atm %>% - dplyr::filter(date >= "2017-05-01" & date <= "2020-04-30") - } else { - atm %>% - dplyr::filter(date >= "2017-05-01" & date <= "2023-12-31") - } + if (soil_depth_cm != 200) { + soil_profile <- kwb.hydrus1d::read_profile(paths$profile) + profile_extended <- kwb.hydrus1d::extend_soil_profile( + soil_profile$profile, + x_end = -soil_depth_cm + ) + soil_profile_extended <- soil_profile + soil_profile_extended$profile <- profile_extended + kwb.hydrus1d::write_profile(soil_profile_extended, path = p("profile")) + } - days_monthly <- lubridate::days_in_month(seq.Date(from = min(atm$date), - to = max(atm$date), - by = "month")) + string_irrig_int <- stringr::str_extract(model_dir, "[0-9][0-9]?days") + # Return the string that is used as "irrig_interval" + paste( + as.integer(stringr::str_extract(string_irrig_int, "\\d+")), + stringr::str_extract(string_irrig_int, "[a-z]+") + ) +} - seq_start <- seq(1,nrow(soil_columns),10) +# sum_per_interval ------------------------------------------------------------- +sum_per_interval <- function(data, interval) +{ + data_org <- data + data <- dplyr::select(data, tidyselect::all_of( + c("date", "groundwater.mmPerDay", "clearwater.mmPerDay") + )) + cols_sum <- setdiff(names(data), "date") + data_summed <- data %>% + dplyr::mutate( + group = lubridate::floor_date(date, unit = interval) + ) %>% # Konvertiere date in datetime-Format + dplyr::group_by(group) %>% # Gruppiere nach Zeitintervallen + dplyr::summarise_at( + .vars = tidyselect::all_of(cols_sum), + .funs = sum + ) %>% # Berechne die Summe für jedes Intervall + dplyr::rename(date = group) + data_org[, cols_sum] <- 0 + data_org[data_org$date %in% data_summed$date, cols_sum] <- data_summed[, cols_sum] + data_org +} - seq_end <- if(nrow(soil_columns) > 10) { - seq(10,nrow(soil_columns),10) - } else { - nrow(soil_columns) - } +# get_valid_exe_path ----------------------------------------------------------- +get_valid_exe_path <- function(exe_dir) +{ + if (file.exists(file.path(exe_dir, "H1D_CALC.exe"))) { + file.path(exe_dir, "H1D_CALC.exe") + } else { + kwb.hydrus1d::check_hydrus_exe() + } +} +# provide_paths ---------------------------------------------------------------- +provide_paths <- function(config, start, end) +{ + #Y:\WWT_Department\Projects\FlexTreat\Work-packages\AP3\3_1_4_Prognosemodell\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no + + # Define a path grammar + PATH_GRAMMAR <- list( + #extdata = system.file("extdata", package = "flextreat.hydrus1d"), + + ###root_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no", + + ###root_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed////", + root_local = path.expand("~/Projekte/flextreat///"), + + #root_local = sprintf("D:/hydrus1d/irrig_fixed/%s/%s/retardation_yes", irrig_dir_string, sprintf("%s%s", duration_string, extreme_rain_string)), + #root_local = "C:/kwb/projects/flextreat/hydrus/Szenarien_10day", + #root_local = system.file("extdata/model", package = "flextreat.hydrus1d"), + exe_dir = "", + model_name_org = "model_to_copy", + model_name = "__soil-column_", + ###model_gui_path_org = "/.h1d", + model_gui_path_org = "/.h1d", + + model_gui_path = "/.h1d", + modelvs_gui_path = "/_vs.h1d", + + ###model_dir_org = "/", + model_dir_org = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/Hydrus1D/irrig_fixed///", + + model_dir = "/", + model_dir_vs = "/_vs", + atmosphere = "/ATMOSPH.IN", + atmosphere_vs = "/ATMOSPH.IN", + a_level = "/A_LEVEL.out", + hydrus1d = "/HYDRUS1D.DAT", + selector = "/SELECTOR.IN", + profile = "/PROFILE.DAT", + obs_node = "/Obs_Node.out", + balance = "/BALANCE.out", + t_level = "/T_LEVEL.out", + t_level_vs = "/T_LEVEL.out", + runinf = "/Run_Inf.out", + solute_id = "1", + solute = "/solute.out", + solute_vs = "/solute.out", + soil_data = "/input-data/soil/soil_geolog.csv", + solute_id_start = sprintf("%02d", start), + solute_id_end = sprintf("%02d", end), + location = if (tracer) { + "tracer" + } else { + sprintf("ablauf_%s_median", config$treatment) + } #"ablauf_ka_median" + ) + + tracer <- config$treatment == "tracer" + + # Resolve the path grammar by replacing the placeholders recursively + kwb.utils::resolve( + PATH_GRAMMAR, + irrig_dir_string = ifelse( + config$irrig_only_growing_season, + "irrig-period_growing-season", + "irrig-period_status-quo" + ), + duration_string = config$duration_string, + extreme_rain_string = if (config$extreme_rain != "") { + paste0("_", config$extreme_rain) + } else { + "" + }, + final_subdir = ifelse(tracer, "tracer", config$retardation_scenario), + scenario = config$scenario + ) +} - if(length(seq_end) < length(seq_start)) seq_end <- c(seq_end, nrow(soil_columns)) +# inner_function --------------------------------------------------------------- +inner_function <- function(atm_data, soil_columns, config, helper) +{ + { + # Define constants + IRRIGATION_COLUMNS <- c("groundwater.mmPerDay", "clearwater.mmPerDay") - solute_ids <- tibble::tibble(start = as.integer(seq_start), - end = as.integer(seq_end)) + # Provide variables from config + extreme_rain <- if (config$extreme_rain == "") { + NULL + } else { + config$extreme_rain + } + tracer <- config$treatment == "tracer" + retardation <- config$retardation == "retardation_yes" + if (retardation) { + soil_columns$retard <- 1 + } - periods <- tibble::tibble( - start = seq(1,length(days_monthly),10), - end = if(length(days_monthly) %% 10 != 0) { - c(seq(10,length(days_monthly),10), length(days_monthly)) - } else { - seq(10,length(days_monthly),10) - } - ) + atm <- get_atm(atm_data, extreme_rain) - loop_df <- if(tracer) { - periods - } else { - solute_ids - } + if (config$irrig_only_growing_season) { + atm[which(!lubridate::month(atm$date) %in% 4:9), IRRIGATION_COLUMNS] <- 0 + } + if (config$duration_string == "short") { + atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2020-04-30") + } else { + atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2023-12-31") + } + days_monthly <- lubridate::days_in_month( + seq.Date(from = min(atm$date), to = max(atm$date), by = "month") + ) - kwb.utils::catAndRun(messageText = sprintf("Running '%s' period for treatment '%s'", - extreme_rain, - treatment), - expr = { + loop_df <- if (tracer) { + helper("generate_periods")(n = length(days_monthly)) + } else { + helper("generate_solute_ids")(n = nrow(soil_columns)) + } + } + kwb.utils::catAndRun( + messageText = sprintf( + "Running '%s' period for treatment '%s'", + extreme_rain, + config$treatment + ), + expr = { sapply(seq_len(nrow(loop_df)), function(i) { - paths_list <- list( - #extdata = system.file("extdata", package = "flextreat.hydrus1d"), - #root_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no", - root_local = sprintf("C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/%s/%s/%s", - irrig_dir_string, - sprintf("%s%s", duration_string, extreme_rain_string), - str_final_dir), - #root_local = sprintf("D:/hydrus1d/irrig_fixed/%s/%s/retardation_yes", irrig_dir_string, sprintf("%s%s", duration_string, extreme_rain_string)), - #root_local = "C:/kwb/projects/flextreat/hydrus/Szenarien_10day", - #root_local = system.file("extdata/model", package = "flextreat.hydrus1d"), - exe_dir = "", - solute_id_start = sprintf("%02d", loop_df$start[i]), - solute_id_end = sprintf("%02d", loop_df$end[i]), - # months_start = periods$start[i], - # months_end = periods$end[i], - scenario = scenario, - location = if(tracer) {"tracer"} else {sprintf("ablauf_%s_median", treatment)}, #"ablauf_ka_median", - model_name_org = "model_to_copy", - model_name = "__soil-column_", - model_gui_path_org = "/.h1d", - model_gui_path = "/.h1d", - modelvs_gui_path = "/_vs.h1d", - model_dir_org = "/", - model_dir = "/", - model_dir_vs = "/_vs", - atmosphere = "/ATMOSPH.IN", - atmosphere_vs = "/ATMOSPH.IN", - a_level = "/A_LEVEL.out", - hydrus1d = "/HYDRUS1D.DAT", - selector = "/SELECTOR.IN", - profile = "/PROFILE.DAT", - obs_node = "/Obs_Node.out", - balance = "/BALANCE.out", - t_level = "/T_LEVEL.out", - t_level_vs = "/T_LEVEL.out", - runinf = "/Run_Inf.out", - solute_id = "1", - solute = "/solute.out", - solute_vs = "/solute.out", - soil_data = "/input-data/soil/soil_geolog.csv" - ) - - - paths <- kwb.utils::resolve(paths_list) - - solute_start_id <- as.numeric(paths$solute_id_start) - solute_end_id <- as.numeric(paths$solute_id_end) - n_solutes <- solute_end_id - (solute_start_id - 1) - - no_irrig <- stringr::str_detect(paths$model_dir, "no-irrig") - irrig_pattern <- "irrig-[0-9][0-9]?days" - irrig_int <- stringr::str_detect(paths$model_dir, irrig_pattern) - - if(irrig_int) { - - string_irrig <- stringr::str_extract(paths$model_dir, sprintf("_%s", irrig_pattern)) - - if(!fs::dir_exists(paths$model_dir)) { - - fs::dir_copy(paths$model_dir_org, paths$model_dir) - } - - if(!fs::file_exists(paths$model_gui_path)) { - fs::file_copy(paths$model_gui_path_org, paths$model_gui_path) - - } - - model_out_files <- list.files(paths$model_dir, pattern = "\\.out$", full.names = TRUE) - - if(length(model_out_files) > 0) { - fs::file_delete(model_out_files) - } - - - soil_depth_cm <- 100 *stringr::str_extract(paths$model_name, "soil-[0-9]+?m") %>% - stringr::str_extract("[0-9]") %>% as.numeric() - - if(soil_depth_cm != 200) { - soil_profile <- kwb.hydrus1d::read_profile(paths$profile) - profile_extended <- kwb.hydrus1d::extend_soil_profile(soil_profile$profile, - x_end = -soil_depth_cm) - - soil_profile_extended <- soil_profile - soil_profile_extended$profile <- profile_extended - - kwb.hydrus1d::write_profile(soil_profile_extended, - path = paths$profile) - } + #i <- 1L + { + paths <- helper("provide_paths")( + config, + start = loop_df$start[i], + end = loop_df$end[i] + ) - string_irrig_int <- stringr::str_extract(paths$model_dir, "[0-9][0-9]?days") + solute_start_id <- as.numeric(paths$solute_id_start) + solute_end_id <- as.numeric(paths$solute_id_end) + n_solutes <- solute_end_id - (solute_start_id - 1) - irrig_interval <- sprintf("%s %s", - string_irrig_int %>% - stringr::str_extract("\\d+") %>% - as.integer(), - string_irrig_int %>% - stringr::str_extract("[a-z]+")) + no_irrig <- stringr::str_detect(paths$model_dir, "no-irrig") + irrig_int <- stringr::str_detect(paths$model_dir, "irrig-[0-9][0-9]?days") } - - #no-irrigation - if(no_irrig) atm[,c("groundwater.mmPerDay", "clearwater.mmPerDay")] <- 0 - - - sum_per_interval <- function(data, interval) { - - data_org <- data - - data <- data %>% - dplyr::select(tidyselect::all_of(c("date", - "groundwater.mmPerDay", - "clearwater.mmPerDay"))) - - cols_sum <- names(data)[names(data) != "date"] - - data_summed <- data %>% - dplyr::mutate(group = lubridate::floor_date(date, unit = interval)) %>% # Konvertiere date in datetime-Format - dplyr::group_by(group) %>% # Gruppiere nach Zeitintervallen - dplyr::summarise_at(.vars = tidyselect::all_of(cols_sum), - .funs = sum) %>% # Berechne die Summe für jedes Intervall - dplyr::rename(date = group) - - data_org[, cols_sum] <- 0 - data_org[data_org$date %in% data_summed$date, cols_sum] <- data_summed[,cols_sum] - data_org + if (irrig_int) { + irrig_interval <- helper("prepare_files_for_irrig_int")(paths) } - if(irrig_int) { - - atm <- sum_per_interval(data = atm, interval = irrig_interval) + # no-irrigation + if (no_irrig) { + atm[, IRRIGATION_COLUMNS] <- 0 + } + if (irrig_int) { + atm <- helper("sum_per_interval")(data = atm, interval = irrig_interval) } days_total <- cumsum(days_monthly) - indeces <- as.integer(1:length(days_total)) + indices <- seq_along(days_total) - c_tops <- lapply(indeces, function(i) { + c_tops <- lapply(indices, function(i) { x <- rep(0, nrow(atm)) - if(i == 1) { - x_min = 1 - } else { - x_min = days_total[i - 1] + 1 - } + x_min <- ifelse(i == 1, 1, days_total[i - 1] + 1) x[x_min:days_total[i]] <- rep(100, days_monthly[i]) - tib <- data.frame(x) - colnames(tib) <- if(i == indeces[1]) { - "cTop"} else { - sprintf("cTop%d", which(indeces %in% i)) - } + + colnames(tib) <- if (i == indices[1]) { + "cTop" + } else { + sprintf("cTop%d", which(indices %in% i)) + } tib }) %>% dplyr::bind_cols() - c_tops_sel <- c_tops[,paths$solute_id_start:paths$solute_id_end] - + c_tops_sel <- c_tops[, paths$solute_id_start:paths$solute_id_end] - atm_prep <- if(tracer) { - flextreat.hydrus1d::prepare_atmosphere(atm = atm, - conc_irrig_clearwater = c_tops_sel, - conc_irrig_groundwater = c_tops_sel, - conc_rain = c_tops_sel) + atm_prep <- if (tracer) { + flextreat.hydrus1d::prepare_atmosphere( + atm = atm, + conc_irrig_clearwater = c_tops_sel, + conc_irrig_groundwater = c_tops_sel, + conc_rain = c_tops_sel + ) } else { - flextreat.hydrus1d::prepare_atmosphere(atm = atm, - conc_irrig_clearwater = soil_columns[solute_start_id:solute_end_id, paths$location], - conc_irrig_groundwater = 0, - conc_rain = 0 + flextreat.hydrus1d::prepare_atmosphere( + atm = atm, + conc_irrig_clearwater = soil_columns[solute_start_id:solute_end_id, paths$location], + conc_irrig_groundwater = 0, + conc_rain = 0 ) } n_tsteps <- nrow(atm_prep) - writeLines(kwb.hydrus1d::write_atmosphere(atm = atm_prep), - paths$atmosphere) - + writeLines( + kwb.hydrus1d::write_atmosphere(atm = atm_prep), + paths$atmosphere + ) selector <- kwb.hydrus1d::read_selector(path = paths$selector) @@ -512,39 +558,30 @@ sapply(extreme_rains, function(extreme_rain) { selector$time$TPrint[1] <- 1 } - solutes_parameters <- if(tracer) { - - } else { - soil_columns[solute_start_id:solute_end_id,] - } - - solutes_new <- prepare_solute_input(dat = soil_columns[solute_start_id:solute_end_id,], - Ks = if(tracer) {0} else {NULL}, - SnkL1 = if(tracer) {0} else {NULL}, - selector = selector, - diff_w = 0, - diff_g = 0, - kd = kd, - halftime_to_firstorderrate = halftime_to_firstorderrate) - + solutes_new <- helper("prepare_solute_input")( + dat = soil_columns[solute_start_id:solute_end_id,], + Ks = if (tracer) 0, # else NULL + SnkL1 = if (tracer) 0, # else NULL + selector = selector, + diff_w = 0, + diff_g = 0, + kd = kd, + halftime_to_firstorderrate = halftime_to_firstorderrate + ) kwb.hydrus1d::write_selector(solutes_new, paths$selector) - - hydrus1d <- kwb.hydrus1d::read_hydrus1d(paths$hydrus1d) hydrus1d$Main$NumberOfSolutes <- n_solutes kwb.hydrus1d::write_hydrus1d(hydrus1d, paths$hydrus1d) - exe_path <- if(file.exists(file.path(paths$exe_dir, "H1D_CALC.exe"))) { - file.path(paths$exe_dir, "H1D_CALC.exe") - } else { - kwb.hydrus1d::check_hydrus_exe() - } + exe_path <- helper("get_valid_exe_path")(paths$exe_dir) - kwb.hydrus1d::run_model(exe_path = exe_path, - model_path = paths$model_dir, - print_output = FALSE) + kwb.hydrus1d::run_model( + exe_path = exe_path, + model_path = paths$model_dir, + print_output = FALSE + ) atmos <- kwb.hydrus1d::read_atmosph(paths$atmosphere) @@ -567,499 +604,550 @@ sapply(extreme_rains, function(extreme_rain) { model_vs_out_files <- list.files(paths$model_dir_vs, pattern = "\\.out$", full.names = TRUE) - if(length(model_vs_out_files) > 0) { + if (length(model_vs_out_files) > 0) { fs::file_delete(model_vs_out_files) } - writeLines(kwb.hydrus1d::write_atmosphere(atm = atmos$data), - paths$atmosphere_vs) - - kwb.hydrus1d::run_model(exe_path = exe_path, - model_path = paths$model_dir_vs, - print_output = FALSE) - - })},dbg = TRUE)} - - - #sapply(scenarios, inner_function) - ncores <- parallel::detectCores() - 1L - cl <- parallel::makeCluster(ncores) - on.exit(parallel::stopCluster(cl)) - - parallel::parLapply( - cl = cl, - X = scenarios, - fun = inner_function, - treatment = treatment, - irrig_only_growing_season = irrig_only_growing_season, - irrig_dir_string = if(irrig_only_growing_season) { - "irrig-period_growing-season" - } else { - "irrig-period_status-quo" - }, - short = short, - tracer = treatment == "tracer", - duration_string = ifelse(short, "short", "long"), - extreme_rain = extreme_rain, - extreme_rain_string = if(any(c("dry", "wet") %in% extreme_rain)) { - sprintf("_%s", extreme_rain) - } else { - "" - }, - str_final_dir = str_final_dir, - soil_columns = soil_columns_selected, - atm = flextreat.hydrus1d::prepare_atmosphere_data(), - get_atm = get_atm, - kd = kd, - halftime_to_firstorderrate = halftime_to_firstorderrate, - prepare_solute_input = prepare_solute_input, - check_hydrus_exe = kwb.hydrus1d::check_hydrus_exe, - base_dir = create_temp_dir(), - create_temp_dir = create_temp_dir, - `%>%` = magrittr::`%>%` - ) - }) -}) - - -solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>% - dplyr::mutate(difftime = c(0,diff(time))) - -max(solute$sum_cv_top) + min(solute$sum_cv_bot) + min(solute$cv_ch1) - -plot(solute$time, solute$c_top) -points(solute$c_bot, col = "red") -(1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop)) * 100 - - -paths$solute_vs2 <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out" - -solute <- kwb.hydrus1d::read_solute(paths$solute) %>% - dplyr::mutate(difftime = c(0,diff(time))) - -(1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 - + writeLines( + kwb.hydrus1d::write_atmosphere(atm = atmos$data), + paths$atmosphere_vs + ) -sum(atmos$data$Prec) -max(tlevel$sum_infil) -max(tlevel$sum_evap) + kwb.hydrus1d::run_model( + exe_path = exe_path, + model_path = paths$model_dir_vs, + print_output = FALSE + ) -balance <- kwb.hydrus1d::read_balance(paths$balance) + }) + }, + dbg = TRUE + ) +} -balance %>% - dplyr::filter(subdomain_id == 0, - time < 400, - parameter == "CncBalT") %>% - ggplot2::ggplot(mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))) + ggplot2::geom_point() +# helper ----------------------------------------------------------------------- +helper <- kwb.utils::createAccessor(list( + generate_solute_ids = generate_solute_ids, + provide_paths = provide_paths, + prepare_files_for_irrig_int = prepare_files_for_irrig_int, + sum_per_interval = sum_per_interval, + prepare_solute_input = prepare_solute_input, + halftime_to_firstorderrate = halftime_to_firstorderrate, + get_valid_exe_path = get_valid_exe_path +)) +# Main loop -------------------------------------------------------------------- +if (FALSE) { -balance %>% - dplyr::filter(subdomain_id == 0, - time < 400, - parameter == "CncBalR") %>% - ggplot2::ggplot(mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))) + - ggplot2::geom_point() + atm_data <- flextreat.hydrus1d::prepare_atmosphere_data() -sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) + configs <- lapply(seq_len(nrow(arg_combis)), function(i) { + as.list(arg_combis[i, ]) + }) -solute_day <- flextreat.hydrus1d::aggregate_solute(solute) + # Sequential loop + for (config in configs) { + #config <- configs[[1L]] + inner_function( + atm_data = atm_data, + soil_columns = soil_columns_selected, + config = config, + helper = helper + ) + } + # Parallel loop + ncores <- parallel::detectCores() - 1L + cl <- parallel::makeCluster(ncores) + + parallel::parLapply( + cl = cl, + X = scenarios, + fun = inner_function, + atm_data = atm_data, + soil_columns = soil_columns_selected, + config = config, + helper = helper + ) + + on.exit(parallel::stopCluster(cl)) +} -plot(solute_day$date, solute_day$c_top) -abline(h = 7) +# After main loop -------------------------------------------------------------- +if (FALSE) +{ + solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>% + dplyr::mutate(difftime = c(0,diff(time))) -plot(atmos$data$tAtm, atmos$data$cTop) + max(solute$sum_cv_top) + min(solute$sum_cv_bot) + min(solute$cv_ch1) -sum(solute_day$cv_top[solute_day$cv_top > 0]) -sum(atmos$data$Prec*atmos$data$cTop) -sum(atmos$data$Prec*atmos$data$cTop)/sum(solute_day$cv_top[solute_day$cv_top > 0]) -sum(solute$cv_top) + plot(solute$time, solute$c_top) + points(solute$c_bot, col = "red") + (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop)) * 100 -sum((as.numeric(solute$cv_top) * solute$difftime)) -sum((solute$cv_ch1 * c(0,diff(solute$time)))) + paths$solute_vs2 <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out" -max(solute$sum_cv_top) + solute <- kwb.hydrus1d::read_solute(paths$solute) %>% + dplyr::mutate(difftime = c(0,diff(time))) -condition <- solute$cv_top > 0 -sum(solute$cv_top[condition]) -condition <- solute$cv_top < 0 -sum(solute$cv_top[condition]) + (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 -solute_aggr_date <- flextreat.hydrus1d::aggregate_solute(solute) + sum(atmos$data$Prec) + max(tlevel$sum_infil) + max(tlevel$sum_evap) -obsnode <- kwb.hydrus1d::read_obsnode(paths$obs_node) + balance <- kwb.hydrus1d::read_balance(paths$balance) -obsnode %>% - dplyr::filter(variable == "flux") %>% - dplyr::group_by(node_id) %>% - dplyr::summarise(sum = sum(value)) + balance %>% + dplyr::filter( + subdomain_id == 0, + time < 400, + parameter == "CncBalT" + ) %>% + ggplot2::ggplot( + mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id)) + ) + ggplot2::geom_point() -profile <- kwb.hydrus1d::read_profile(paths$profile) + balance %>% + dplyr::filter( + subdomain_id == 0, + time < 400, + parameter == "CncBalR" + ) %>% + ggplot2::ggplot(mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))) + + ggplot2::geom_point() -p <- obsnode %>% - dplyr::left_join(profile[,c("node_id", "x")]) %>% - dplyr::filter(variable == "conc1") %>% - ggplot2::ggplot(mapping = ggplot2::aes(x = time, - y = value, - col = as.factor(x))) + - ggplot2::geom_point() + - ggplot2::theme_bw() + sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) -plotly::ggplotly(p) + solute_day <- flextreat.hydrus1d::aggregate_solute(solute) -solute_aggr_date -View(tlevel_aggr_date) -View(solute_aggr_date) + plot(solute_day$date, solute_day$c_top) + abline(h = 7) + plot(atmos$data$tAtm, atmos$data$cTop) -t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) -t_level + sum(solute_day$cv_top[solute_day$cv_top > 0]) + sum(atmos$data$Prec*atmos$data$cTop) + sum(atmos$data$Prec*atmos$data$cTop)/sum(solute_day$cv_top[solute_day$cv_top > 0]) + sum(solute$cv_top) -## t_level aggregate -tlevel_aggr_date <- flextreat.hydrus1d::aggregate_tlevel(t_level) -tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel(t_level, - col_aggr = "yearmonth") -tlevel_aggr_year_hydrologic <- flextreat.hydrus1d::aggregate_tlevel(t_level, - col_aggr = "year_hydrologic") %>% - dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october + sum((as.numeric(solute$cv_top) * solute$difftime)) + sum((solute$cv_ch1 * c(0,diff(solute$time)))) + max(solute$sum_cv_top) -wb_date_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_date) -wb_yearmonth_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_yearmonth) -wb_yearhydrologic_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_year_hydrologic) + condition <- solute$cv_top > 0 + sum(solute$cv_top[condition]) + condition <- solute$cv_top < 0 + sum(solute$cv_top[condition]) -wb_date_plot -wb_yearmonth_plot -wb_yearhydrologic_plot + solute_aggr_date <- flextreat.hydrus1d::aggregate_solute(solute) + obsnode <- kwb.hydrus1d::read_obsnode(paths$obs_node) -solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))] + obsnode %>% + dplyr::filter(variable == "flux") %>% + dplyr::group_by(node_id) %>% + dplyr::summarise(sum = sum(value)) -paths$model_dir_vs + profile <- kwb.hydrus1d::read_profile(paths$profile) -#scenarios_solutes <- paste0(scenarios, "_soil-column") -scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median") + p <- obsnode %>% + dplyr::left_join(profile[,c("node_id", "x")]) %>% + dplyr::filter(variable == "conc1") %>% + ggplot2::ggplot(mapping = ggplot2::aes( + x = time, + y = value, + col = as.factor(x)) + ) + + ggplot2::geom_point() + + ggplot2::theme_bw() + plotly::ggplotly(p) -scenario_dirs <- fs::dir_ls(path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed", recurse = TRUE, regexp = "retardation.*vs$", type = "directory") + solute_aggr_date + View(tlevel_aggr_date) + View(solute_aggr_date) + t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) + t_level -sapply(scenario_dirs, function(scenario_dir) { + ## t_level aggregate + tlevel_aggr_date <- flextreat.hydrus1d::aggregate_tlevel(t_level) - solutes_list <- setNames(lapply(scenarios_solutes, function(scenario) { - solute_files <- fs::dir_ls(scenario_dir, - regexp = "solute\\d\\d?.out", - recurse = TRUE) + tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel( + t_level, + col_aggr = "yearmonth" + ) - stopifnot(length(solute_files) > 0) + tlevel_aggr_year_hydrologic <- flextreat.hydrus1d::aggregate_tlevel( + t_level, + col_aggr = "year_hydrologic" + ) %>% + dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october - solutes <- setNames(lapply(solute_files, function(file) { - kwb.hydrus1d::read_solute(file, dbg = TRUE) - }), nm = solute_files) %>% dplyr::bind_rows(.id = "path") + wb_date_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_date) + wb_yearmonth_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_yearmonth) + wb_yearhydrologic_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_year_hydrologic) + wb_date_plot + wb_yearmonth_plot + wb_yearhydrologic_plot - solute_files_df <- tibble::tibble(path = solute_files, - model_solute_id = path %>% basename() %>% stringr::str_extract(pattern = "[0-9][0-9]?") %>% as.integer(), - soilcolumn_id_start = path %>% dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(1,2) %>% as.integer(), - soilcolumn_id_end = path %>% dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(3,4) %>% as.integer(), - soil_column_id = soilcolumn_id_start + model_solute_id - 1) %>% - dplyr::left_join(soil_columns, by = c(soil_column_id = "id")) + solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))] + paths$model_dir_vs - dplyr::left_join(solutes, solute_files_df) - }), nm = scenarios_solutes) + #scenarios_solutes <- paste0(scenarios, "_soil-column") + scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median") - solutes_df <- solutes_list %>% - dplyr::bind_rows(.id = "scenario") + scenario_dirs <- fs::dir_ls( + path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed", + recurse = TRUE, + regexp = "retardation.*vs$", + type = "directory" + ) - solutes_df_stats <- solutes_df %>% - dplyr::bind_rows(.id = "scenario") %>% - dplyr::mutate(scen = stringr::str_remove(basename(dirname(path)), "_soil-column.*")) %>% - dplyr::group_by(path, scen,substanz_nr, substanz_name) %>% - dplyr::summarise(sum_cv_top = max(sum_cv_top), - sum_cv_bot = min(sum_cv_bot), - cv_ch1 = min(cv_ch1)) %>% - dplyr::mutate(mass_balance_error_percent = 100*(sum_cv_top + cv_ch1 + sum_cv_bot)/sum_cv_top) %>% - dplyr::arrange(mass_balance_error_percent) + sapply(scenario_dirs, function(scenario_dir) { + solutes_list <- lapply( + setNames(nm = scenarios_solutes), + function(scenario) { - solutes_df_stats$soil <- solutes_df_stats$sum_cv_top + solutes_df_stats$sum_cv_bot + solutes_df_stats$cv_ch1 + solute_files <- fs::dir_ls( + scenario_dir, + regexp = "solute\\d\\d?.out", + recurse = TRUE + ) - saveRDS(solutes_df, file = file.path(scenario_dir, "solutes.rds")) + stopifnot(length(solute_files) > 0) - openxlsx::write.xlsx(solutes_df_stats, file = file.path(scenario_dir, "hydrus_scenarios.xlsx")) -}) + solutes <- setNames(lapply(solute_files, function(file) { + kwb.hydrus1d::read_solute(file, dbg = TRUE) + }), nm = solute_files) %>% dplyr::bind_rows(.id = "path") + solute_files_df <- tibble::tibble( + path = solute_files, + model_solute_id = path %>% basename() %>% stringr::str_extract(pattern = "[0-9][0-9]?") %>% as.integer(), + soilcolumn_id_start = path %>% dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(1,2) %>% as.integer(), + soilcolumn_id_end = path %>% dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(3,4) %>% as.integer(), + soil_column_id = soilcolumn_id_start + model_solute_id - 1 + ) %>% + dplyr::left_join(soil_columns, by = c(soil_column_id = "id")) -scenario_dirs <- fs::dir_ls(path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", recurse = TRUE, regexp = "retardation_no/hydrus_scenarios.xlsx$", type = "file") + dplyr::left_join(solutes, solute_files_df) + } + ) -res_stats <- stats::setNames(lapply(scenario_dirs, function(scenario_dir) { - readxl::read_excel(scenario_dir) -}), nm = scenario_dirs) + solutes_df <- solutes_list %>% + dplyr::bind_rows(.id = "scenario") -View(solutes_list$`soil-1m_irrig-01days_soil-column`) + solutes_df_stats <- solutes_df %>% + dplyr::bind_rows(.id = "scenario") %>% + dplyr::mutate(scen = stringr::str_remove(basename(dirname(path)), "_soil-column.*")) %>% + dplyr::group_by(path, scen,substanz_nr, substanz_name) %>% + dplyr::summarise(sum_cv_top = max(sum_cv_top), + sum_cv_bot = min(sum_cv_bot), + cv_ch1 = min(cv_ch1)) %>% + dplyr::mutate(mass_balance_error_percent = 100*(sum_cv_top + cv_ch1 + sum_cv_bot)/sum_cv_top) %>% + dplyr::arrange(mass_balance_error_percent) -model_paths <- fs::dir_ls("C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", - recurse = TRUE, - regexp = "tracer$", - type = "directory") + solutes_df_stats$soil <- solutes_df_stats$sum_cv_top + solutes_df_stats$sum_cv_bot + solutes_df_stats$cv_ch1 -traveltimes_list <- lapply(model_paths, function(model_path) { - setNames(lapply(scenarios, function(scenario) { + saveRDS(solutes_df, file = file.path(scenario_dir, "solutes.rds")) - try({ + openxlsx::write.xlsx(solutes_df_stats, file = file.path(scenario_dir, "hydrus_scenarios.xlsx")) + }) - solute_files <- fs::dir_ls(path = model_path, - recurse = TRUE, - regexp = sprintf("tracer_%s_.*vs/solute\\d\\d?.out", scenario) - ) + scenario_dirs <- fs::dir_ls( + path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", + recurse = TRUE, + regexp = "retardation_no/hydrus_scenarios.xlsx$", + type = "file" + ) + + res_stats <- lapply( + stats::setNames(nm = scenario_dirs), + function(scenario_dir) { + readxl::read_excel(scenario_dir) + } + ) + + View(solutes_list$`soil-1m_irrig-01days_soil-column`) + + model_paths <- fs::dir_ls( + "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", + recurse = TRUE, + regexp = "tracer$", + type = "directory" + ) + + traveltimes_list <- lapply(model_paths, function(model_path) { + setNames(nm = (scenarios), lapply(scenarios, function(scenario) { + try({ + solute_files <- fs::dir_ls( + path = model_path, + recurse = TRUE, + regexp = sprintf("tracer_%s_.*vs/solute\\d\\d?.out", scenario) + ) + flextreat.hydrus1d::get_traveltimes(solute_files, dbg = TRUE) + }) + })) + }) - flextreat.hydrus1d::get_traveltimes(solute_files, dbg = TRUE) - })}), nm = (scenarios)) -}) + sapply(seq_along(traveltimes_list), function(i) { + label <- stringr::str_remove( + names(traveltimes_list)[i], + "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/" + ) -sapply(seq_along(traveltimes_list), function(i) { + htmlwidgets::saveWidget( + flextreat.hydrus1d::plot_traveltimes( + dplyr::bind_rows(traveltimes_list[[i]]), + title = label, + ylim = c(0, 650) + ), + file = sprintf("traveltimes_%s.html", label)) + }) - label <- stringr::str_remove(names(traveltimes_list)[i], "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/") - htmlwidgets::saveWidget(flextreat.hydrus1d::plot_traveltimes(dplyr::bind_rows(traveltimes_list[[i]]), - title = label, - ylim = c(0,650)), - file = sprintf("traveltimes_%s.html", label)) -}) + # extrahiere_letzte_drei_teile ------------------------------------------------- + extrahiere_letzte_drei_teile <- function(pfad) + { + sapply(pfad, function(pf) { + # Teile den Pfad anhand des Schrägstrichs auf + teile <- unlist(strsplit(pf, "/")) + # Waehle die letzten drei Teile aus + letzte_drei_teile <- tail(teile, 3) -extrahiere_letzte_drei_teile <- function(pfad) { + paste0(letzte_drei_teile, collapse = "_") + }) + } - sapply(pfad, function(pf) { - # Teile den Pfad anhand des Schrägstrichs auf - teile <- unlist(strsplit(pf, "/")) + # Plotting --------------------------------------------------------------------- + if (FALSE) + { + pdff <- "traveltimes_per-scenario.pdf" + kwb.utils::preparePdf(pdff) + sapply(seq_along(traveltimes_list), function(i) { + traveltimes_list_sel <- traveltimes_list[[i]] + + label <- extrahiere_letzte_drei_teile(names(traveltimes_list)[i]) + + traveltime_bp <- lapply(traveltimes_list_sel, function(x) { + x %>% + dplyr::filter(percentiles == 0.5) + }) %>% dplyr::bind_rows(.id = "scenario") %>% + dplyr::filter(!stringr::str_detect(scenario, "1.5")) %>% + dplyr::mutate(quarter = lubridate::quarter(date) %>% as.factor(), + soil_depth = stringr::str_extract(scenario, "soil-.*m") %>% + stringr::str_remove_all("soil-|m") %>% as.factor()) + + + scenario_by_median_traveltime <- traveltime_bp %>% + dplyr::group_by(scenario) %>% + dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>% + dplyr::arrange(median) + + traveltime_bp <- traveltime_bp %>% + dplyr::left_join(scenario_by_median_traveltime) + + + y_lim <- c(0,350) + + + tt_bp_total <- traveltime_bp %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff)) + + ggplot2::geom_boxplot(outliers = FALSE) + + ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1), + #col = "darkgrey", + alpha = 0.6) + + ggplot2::ylim(y_lim) + + ggplot2::labs(y = "Median Traveltime (days)", + x = sprintf("Scenario (%s)", label), + title = "Boxplot: median traveltime total") + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "top") + + print(tt_bp_total) + }) + + kwb.utils::finishAndShowPdf(pdff) + + pdff <- "traveltimes_all.pdf" + kwb.utils::preparePdf(pdff) + + traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% + dplyr::bind_rows(.id = "scenario_main_raw") %>% + dplyr::mutate( + scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"), + scenario_main = scenario_main_raw %>% extrahiere_letzte_drei_teile(), + quarter = lubridate::quarter(date) %>% as.factor(), + soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% + stringr::str_remove_all("soil-|m") %>% as.factor()) + + scenario_by_median_traveltime <- traveltimes_df %>% + dplyr::group_by(scenario_main, scenario_name) %>% + dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>% + dplyr::arrange(median) + + traveltimes_bp <- traveltimes_df %>% + dplyr::left_join(scenario_by_median_traveltime) + + y_lim <- c(0,350) + + tt_bp_total <- traveltimes_bp %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, median), + y = time_diff, + col = scenario_main)) + + ggplot2::geom_boxplot(outliers = FALSE) + + # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( + # jitter.width = 0.1, + # dodge.width = 0.75), + # alpha = 0.6) + + ggplot2::ylim(y_lim) + + ggplot2::labs(y = "Median Traveltime (days)", + x = "Scenario", + title = "Boxplot: median traveltime total") + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "top") + + print(tt_bp_total) + + kwb.utils::finishAndShowPdf(pdff) + + htmlwidgets::saveWidget( + widget = plotly::ggplotly(tt_bp_total), + file = "traveltimes_all.html" + ) - # Wähle die letzten drei Teile aus - letzte_drei_teile <- tail(teile, 3) + pdff <- "traveltimes_all_percent.pdf" + kwb.utils::preparePdf(pdff) + + traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% + dplyr::bind_rows(.id = "scenario_main_raw") %>% + dplyr::mutate( + scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"), + scenario_main = scenario_main_raw %>% extrahiere_letzte_drei_teile(), + quarter = lubridate::quarter(date) %>% as.factor(), + soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% + stringr::str_remove_all("soil-|m") %>% as.factor()) + + scenario_base_median <- traveltimes_df %>% + dplyr::filter( + scenario_name == "soil-2m_irrig-10days", + scenario_main == "irrig-period_status-quo_long_tracer", + percentiles == 0.5 + ) %>% + dplyr::select(- time_top, - time_bot) %>% + dplyr::rename(time_diff_base = time_diff) + + median_base_percent <- median(scenario_base_median$time_diff_base, na.rm = TRUE) + + scenario_by_median_traveltime <- traveltimes_df %>% + dplyr::group_by(scenario_main, scenario_name) %>% + dplyr::summarise(median = median(time_diff, na.rm = TRUE), + median_percent = round(100 * median / median_base_percent, 2)) %>% + dplyr::arrange(median_percent) + + traveltimes_bp <- traveltimes_df %>% + dplyr::left_join(scenario_base_median) %>% + dplyr::left_join(scenario_by_median_traveltime) - paste0(letzte_drei_teile, collapse = "_") - }) -} + } + # Plotting --------------------------------------------------------------------- -pdff <- "traveltimes_per-scenario.pdf" -kwb.utils::preparePdf(pdff) -sapply(seq_along(traveltimes_list), function(i) { - traveltimes_list_sel <- traveltimes_list[[i]] + y_lim <- c(0,350) - label <- extrahiere_letzte_drei_teile(names(traveltimes_list)[i]) + tt_bp_percent <- traveltimes_bp %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, median), + y = time_diff, + col = scenario_main)) + + ggplot2::geom_boxplot(outliers = FALSE) + + # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( + # jitter.width = 0.1, + # dodge.width = 0.75), + # alpha = 0.6) + + ggplot2::ylim(y_lim) + + ggplot2::labs(y = "Median Traveltime (%) compared to Status Quo", + x = "Scenario", + title = "Boxplot: median traveltime percent") + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "top") - traveltime_bp <- lapply(traveltimes_list_sel, function(x) { - x %>% - dplyr::filter(percentiles == 0.5) - }) %>% dplyr::bind_rows(.id = "scenario") %>% - dplyr::filter(!stringr::str_detect(scenario, "1.5")) %>% - dplyr::mutate(quarter = lubridate::quarter(date) %>% as.factor(), - soil_depth = stringr::str_extract(scenario, "soil-.*m") %>% - stringr::str_remove_all("soil-|m") %>% as.factor()) + print(tt_bp_percent) + kwb.utils::finishAndShowPdf(pdff) - scenario_by_median_traveltime <- traveltime_bp %>% - dplyr::group_by(scenario) %>% - dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>% - dplyr::arrange(median) + htmlwidgets::saveWidget( + widget = plotly::ggplotly(tt_bp_total), + file = "traveltimes_all.html" + ) - traveltime_bp <- traveltime_bp %>% - dplyr::left_join(scenario_by_median_traveltime) + tt_bp_total_soil <- traveltime_bp %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = soil_depth)) + + ggplot2::geom_boxplot(outliers = FALSE) + + ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1), + alpha = 0.6) + + ggplot2::ylim(y_lim) + + ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario", + col = "Soil Depth (m)", + title = "Boxplot: median traveltime total") + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "top") + tt_bp_total_soil - y_lim <- c(0,350) - tt_bp_total <- traveltime_bp %>% + tt_bp_total_quartal <- traveltime_bp %>% ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff)) + ggplot2::geom_boxplot(outliers = FALSE) + ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1), - #col = "darkgrey", + mapping = ggplot2::aes(col = quarter), alpha = 0.6) + ggplot2::ylim(y_lim) + - ggplot2::labs(y = "Median Traveltime (days)", - x = sprintf("Scenario (%s)", label), + ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario", + col = "Quartal", title = "Boxplot: median traveltime total") + ggplot2::theme_bw() + ggplot2::theme(legend.position = "top") - print(tt_bp_total) -}) -kwb.utils::finishAndShowPdfIf(pdff) - - -pdff <- "traveltimes_all.pdf" -kwb.utils::preparePdf(pdff) - -traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% - dplyr::bind_rows(.id = "scenario_main_raw") %>% - dplyr::mutate(scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"), - scenario_main = scenario_main_raw %>% extrahiere_letzte_drei_teile(), - quarter = lubridate::quarter(date) %>% as.factor(), - soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% - stringr::str_remove_all("soil-|m") %>% as.factor()) - - -scenario_by_median_traveltime <- traveltimes_df %>% - dplyr::group_by(scenario_main, scenario_name) %>% - dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>% - dplyr::arrange(median) - -traveltimes_bp <- traveltimes_df %>% - dplyr::left_join(scenario_by_median_traveltime) - - -y_lim <- c(0,350) - - -tt_bp_total <- traveltimes_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, median), - y = time_diff, - col = scenario_main)) + - ggplot2::geom_boxplot(outliers = FALSE) + - # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( - # jitter.width = 0.1, - # dodge.width = 0.75), - # alpha = 0.6) + - ggplot2::ylim(y_lim) + - ggplot2::labs(y = "Median Traveltime (days)", - x = "Scenario", - title = "Boxplot: median traveltime total") + - ggplot2::theme_bw() + - ggplot2::theme(legend.position = "top") - -print(tt_bp_total) - -kwb.utils::finishAndShowPdfIf(to.pdf = TRUE, pdff) - - - -htmlwidgets::saveWidget(widget = plotly::ggplotly(tt_bp_total), - file = "traveltimes_all.html") - - -pdff <- "traveltimes_all_percent.pdf" -kwb.utils::preparePdf(pdff) - -traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% - dplyr::bind_rows(.id = "scenario_main_raw") %>% - dplyr::mutate(scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"), - scenario_main = scenario_main_raw %>% extrahiere_letzte_drei_teile(), - quarter = lubridate::quarter(date) %>% as.factor(), - soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% - stringr::str_remove_all("soil-|m") %>% as.factor()) - -scenario_base_median <- traveltimes_df %>% - dplyr::filter(scenario_name == "soil-2m_irrig-10days", - scenario_main == "irrig-period_status-quo_long_tracer", - percentiles == 0.5 - ) %>% - dplyr::select(- time_top, - time_bot) %>% - dplyr::rename(time_diff_base = time_diff) - -median_base_percent <- median(scenario_base_median$time_diff_base, na.rm = TRUE) - -scenario_by_median_traveltime <- traveltimes_df %>% - dplyr::group_by(scenario_main, scenario_name) %>% - dplyr::summarise(median = median(time_diff, na.rm = TRUE), - median_percent = round(100 * median / median_base_percent, 2)) %>% - dplyr::arrange(median_percent) - -traveltimes_bp <- traveltimes_df %>% - dplyr::left_join(scenario_base_median) %>% - dplyr::left_join(scenario_by_median_traveltime) - - -y_lim <- c(0,350) - - -tt_bp_percent <- traveltimes_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, median), - y = time_diff, - col = scenario_main)) + - ggplot2::geom_boxplot(outliers = FALSE) + - # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( - # jitter.width = 0.1, - # dodge.width = 0.75), - # alpha = 0.6) + - ggplot2::ylim(y_lim) + - ggplot2::labs(y = "Median Traveltime (%) compared to Status Quo", - x = "Scenario", - title = "Boxplot: median traveltime percent") + - ggplot2::theme_bw() + - ggplot2::theme(legend.position = "top") - -print(tt_bp_percent) - -kwb.utils::finishAndShowPdfIf(to.pdf = TRUE, pdff) - - - -htmlwidgets::saveWidget(widget = plotly::ggplotly(tt_bp_total), - file = "traveltimes_all.html") - - + tt_bp_total_quartal -tt_bp_total_soil <- traveltime_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = soil_depth)) + - ggplot2::geom_boxplot(outliers = FALSE) + - ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1), - alpha = 0.6) + - ggplot2::ylim(y_lim) + - ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario", - col = "Soil Depth (m)", - title = "Boxplot: median traveltime total") + - ggplot2::theme_bw() + - ggplot2::theme(legend.position = "top") -tt_bp_total_soil - - - -tt_bp_total_quartal <- traveltime_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff)) + - ggplot2::geom_boxplot(outliers = FALSE) + - ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1), - mapping = ggplot2::aes(col = quarter), - alpha = 0.6) + - ggplot2::ylim(y_lim) + - ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario", - col = "Quartal", - title = "Boxplot: median traveltime total") + - ggplot2::theme_bw() + - ggplot2::theme(legend.position = "top") - -tt_bp_total_quartal - - - -tt_bp_quarter <- traveltime_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = quarter)) + - ggplot2::geom_boxplot(outliers = FALSE) + - ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( - jitter.width = 0.1, - dodge.width = 0.75), - alpha = 0.6) + - ggplot2::labs(y = "Median Traveltime (days)", - x = "Scenario", - col = "Quartal", - title = "Boxplot: median traveltime by quarter") + - ggplot2::ylim(y_lim) + - ggplot2::theme_bw() + - ggplot2::theme(legend.position = "top") - -tt_bp_quarter + tt_bp_quarter <- traveltime_bp %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = quarter)) + + ggplot2::geom_boxplot(outliers = FALSE) + + ggplot2::geom_jitter(position = ggplot2::position_jitterdodge( + jitter.width = 0.1, + dodge.width = 0.75), + alpha = 0.6) + + ggplot2::labs(y = "Median Traveltime (days)", + x = "Scenario", + col = "Quartal", + title = "Boxplot: median traveltime by quarter") + + ggplot2::ylim(y_lim) + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "top") -htmlwidgets::saveWidget(widget = plotly::ggplotly(tt_bp_total), - title = "Boxplot: median traveltime total", - file = "boxplot_traveltimes-median_total.html") + tt_bp_quarter + # Save HTML widgets ------------------------------------------------------------ + htmlwidgets::saveWidget( + widget = plotly::ggplotly(tt_bp_total), + title = "Boxplot: median traveltime total", + file = "boxplot_traveltimes-median_total.html" + ) -htmlwidgets::saveWidget(plotly::ggplotly(tt_bp_quarter), - title = "Boxplot: median traveltime by quarter", - file = "boxplot_traveltimes-median_quarter.html") + htmlwidgets::saveWidget( + plotly::ggplotly(tt_bp_quarter), + title = "Boxplot: median traveltime by quarter", + file = "boxplot_traveltimes-median_quarter.html" + ) +}