From 92a291254c652eaff4759f5409826a51919243e4 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 06:54:49 +0200 Subject: [PATCH 01/29] Add/improve section label --- R/.scenarios_parallel.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index f9c814b..cc5fa98 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -1,9 +1,10 @@ -# Preface ---------------------------------------------------------------------- +# Install packages ------------------------------------------------------------- if (FALSE) { remotes::install_github("kwb-r/kwb.hydrus1d@dev") remotes::install_github("kwb-r/flextreat.hydrus1d@dev") } +# Load packages ---------------------------------------------------------------- library(magrittr) library(flextreat.hydrus1d) From 0f8c88c3b776be612addd65445a76f738483f406 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:16:23 +0200 Subject: [PATCH 02/29] Reindent --- R/.scenarios_parallel.R | 586 +++++++++++++++++++++++----------------- 1 file changed, 340 insertions(+), 246 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index cc5fa98..ddf835c 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -831,49 +831,75 @@ if (FALSE) future.apply::future_sapply(scenario_dirs, function(scenario_dir) { - solutes_list <- setNames(lapply(scenarios_solutes, function(scenario) { - solute_files <- fs::dir_ls(scenario_dir, - regexp = "solute\\d\\d?.out", - recurse = TRUE) - - stopifnot(length(solute_files) > 0) - - 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")) + solutes_list <- setNames( + lapply(scenarios_solutes, function(scenario) { + solute_files <- fs::dir_ls( + scenario_dir, + regexp = "solute\\d\\d?.out", + recurse = TRUE + ) - dplyr::left_join(solutes, solute_files_df) - }), nm = scenarios_solutes) + stopifnot(length(solute_files) > 0) + + 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")) + + dplyr::left_join(solutes, solute_files_df) + }), + nm = scenarios_solutes + ) solutes_df <- solutes_list %>% dplyr::bind_rows(.id = "scenario") solutes_df_stats <- solutes_df %>% dplyr::bind_rows(.id = "scenario") %>% - dplyr::mutate(scen = stringr::str_remove(basename(dirname(path)), "_soil-column.*")) %>% + 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::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) - solutes_df_stats$soil <- solutes_df_stats$sum_cv_top + solutes_df_stats$sum_cv_bot + solutes_df_stats$cv_ch1 saveRDS(solutes_df, file = file.path(scenario_dir, "solutes.rds")) openxlsx::write.xlsx(solutes_df_stats, file = file.path(scenario_dir, "hydrus_scenarios.xlsx")) }) + # Close the parallel plan future::plan(future::sequential) }) @@ -892,38 +918,41 @@ if (FALSE) } ) - # load_default <- res_stats$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/retardation_no/ablauf_ka_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx` - load_default <- res_stats$`C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/irrig-period_status-quo/long/retardation_no/ablauf_o3_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx` %>% - # dplyr::mutate(retardation = basename(dirname(dirname(path))), - # duration = basename(dirname(dirname(dirname(path)))), - # irrigation_period = basename(dirname(dirname(dirname(dirname((path))))))) %>% + # load_default <- res_stats$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/retardation_no/ablauf_ka_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx` + name <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/irrig-period_status-quo/long/retardation_no/ablauf_o3_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx" + load_default <- res_stats[[name]] %>% + # dplyr::mutate( + # retardation = basename(dirname(dirname(path))), + # duration = basename(dirname(dirname(dirname(path)))), + # irrigation_period = basename(dirname(dirname(dirname(dirname((path)))))) + # ) %>% dplyr::select(- path, - scen, - mass_balance_error_percent, - soil) - res_stats_df <- dplyr::bind_rows(res_stats) %>% - dplyr::mutate(retardation = basename(dirname(dirname(path))), - duration = basename(dirname(dirname(dirname(path)))), - irrigation_period = basename(dirname(dirname(dirname(dirname((path))))))) %>% + dplyr::mutate( + retardation = basename(dirname(dirname(path))), + duration = basename(dirname(dirname(dirname(path)))), + irrigation_period = basename(dirname(dirname(dirname(dirname((path)))))) + ) %>% dplyr::select(- path, - mass_balance_error_percent, - soil) - - names(res_stats_df)[4:6] <- paste0("default_", names(res_stats_df)[4:6]) + names(res_stats_df)[4:6] <- paste0("default_", names(res_stats_df)[4:6]) res_stats_df <- res_stats_df %>% dplyr::left_join(load_default, by = c("substanz_nr", "substanz_name")) %>% - dplyr::mutate(percental_load_gw = dplyr::if_else(abs(default_sum_cv_bot) < 10000 | abs(sum_cv_bot) < 10000, - NA_real_, - 100 + 100 * (abs(sum_cv_bot) - abs(default_sum_cv_bot)) / abs(default_sum_cv_bot))) + dplyr::mutate(percental_load_gw = dplyr::if_else( + abs(default_sum_cv_bot) < 10000 | abs(sum_cv_bot) < 10000, + NA_real_, + 100 + 100 * (abs(sum_cv_bot) - abs(default_sum_cv_bot)) / abs(default_sum_cv_bot)) + ) View(res_stats_df) - - - #root_path <- "D:/hydrus1d/irrig_fixed_01" - root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" + #root_path <- "D:/hydrus1d/irrig_fixed_01" + root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" model_paths <- fs::dir_ls( - path = root_path, + path = root_path, recurse = TRUE, regexp = "tracer$", type = "directory" @@ -931,7 +960,10 @@ if (FALSE) #model_paths <- model_paths[stringr::str_detect(model_paths, "wet|dry", negate = TRUE)] - scenarios <- sapply(c(1,10), function(x) paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x))) %>% + scenarios <- sapply( + c(1,10), + function(x) paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x)) + ) %>% as.vector() #unique(configs$scenario) @@ -940,40 +972,40 @@ if (FALSE) ### read traveltimes sequential system.time( - traveltimes_list <- lapply(model_paths, function(model_path) { - setNames(nm = (scenarios), lapply(scenarios, function(scenario) { - try({ - message(sprintf("Scenario: %s", scenario)) - 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) - }) + traveltimes_list <- lapply(model_paths, function(model_path) { + setNames(nm = (scenarios), lapply(scenarios, function(scenario) { + try({ + message(sprintf("Scenario: %s", scenario)) + 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) + }) + })) })) - })) ### read traveltimes in parallel library(future.apply) # Set up parallel plan system.time(expr = { - future::plan(future::multisession) - - traveltimes_list <- future.apply::future_lapply(model_paths, function(model_path) { - setNames(nm = (scenarios), future.apply::future_lapply(scenarios, function(scenario) { - try({ - message(sprintf("Scenario: %s", scenario)) - 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) - }) - })) - }) + future::plan(future::multisession) + + traveltimes_list <- future.apply::future_lapply(model_paths, function(model_path) { + setNames(nm = (scenarios), future.apply::future_lapply(scenarios, function(scenario) { + try({ + message(sprintf("Scenario: %s", scenario)) + 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) + }) + })) + }) } ) @@ -983,13 +1015,16 @@ if (FALSE) sapply(seq_along(traveltimes_list), function(i) { - htmlwidgets::saveWidget(flextreat.hydrus1d::plot_traveltimes(traveltimes_list[[i]] %>% dplyr::bind_rows(), - title = sprintf("%s", extrahiere_letzte_drei_teile(names(traveltimes_list)[i])), - ylim = c(0,650)), - file = sprintf("traveltimes_%s.html", names(traveltimes_list)[i])) + htmlwidgets::saveWidget( + flextreat.hydrus1d::plot_traveltimes( + traveltimes_list[[i]] %>% + dplyr::bind_rows(), + title = sprintf("%s", extrahiere_letzte_drei_teile(names(traveltimes_list)[i])), + ylim = c(0, 650) + ), + file = sprintf("traveltimes_%s.html", names(traveltimes_list)[i])) }) - # extrahiere_letzte_drei_teile ------------------------------------------------- extrahiere_letzte_drei_teile <- function(pfad) { @@ -1004,62 +1039,70 @@ if (FALSE) }) } - # Plotting --------------------------------------------------------------------- + # Plotting if (FALSE) { pdff <- "traveltimes_per-scenario.pdf" kwb.utils::preparePdf(pdff) # sapply(names(traveltimes_list), function(path) { - # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows() - # - # label <- extrahiere_letzte_drei_teile(path) - - # traveltime_bp <- traveltimes_sel %>% - # dplyr::bind_rows() %>% - # 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()) - - traveltime_bp <- lapply(traveltimes_list, function(x) { - x %>% - dplyr::bind_rows() %>% - 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 %>% - dplyr::mutate(scenario_short = extrahiere_letzte_drei_teile(scenario)) %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, 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 = "Scenario", - title = "Boxplot: median traveltime total") + - ggplot2::theme_bw() + - ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust=1)) - - print(tt_bp_total) + # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows() + # + # label <- extrahiere_letzte_drei_teile(path) + + # traveltime_bp <- traveltimes_sel %>% + # dplyr::bind_rows() %>% + # 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()) + + traveltime_bp <- lapply(traveltimes_list, function(x) { + x %>% + dplyr::bind_rows() %>% + 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 %>% + dplyr::mutate(scenario_short = extrahiere_letzte_drei_teile(scenario)) %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, 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 = "Scenario", + title = "Boxplot: median traveltime total" + ) + + ggplot2::theme_bw() + + ggplot2::theme( + axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1) + ) + + print(tt_bp_total) kwb.utils::finishAndShowPdf(pdff) @@ -1067,12 +1110,15 @@ if (FALSE) pdff <- "traveltimes_all.pdf" kwb.utils::preparePdf(pdff) - traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% + traveltimes_df <- lapply(traveltimes_list, 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(), + 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()) @@ -1087,18 +1133,22 @@ if (FALSE) 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::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::labs( + y = "Median Traveltime (days)", + x = "Scenario", + title = "Boxplot: median traveltime total" + ) + ggplot2::theme_bw() + ggplot2::theme(legend.position = "top") @@ -1114,14 +1164,18 @@ if (FALSE) pdff <- "traveltimes_all_percent.pdf" kwb.utils::preparePdf(pdff) - traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>% + traveltimes_df <- lapply(traveltimes_list, 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_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( @@ -1134,123 +1188,163 @@ if (FALSE) traveltimes_bp <- traveltimes_df %>% dplyr::filter(percentiles == 0.5) %>% - dplyr::left_join(scenario_base_median[, c("month_id", "time_diff_base")] %>% dplyr::mutate(percentiles = 0.5)) %>% - dplyr::mutate(time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base) - - # Plotting --------------------------------------------------------------------- - - y_lim <- c(0,350) - - tt_bp_percent <- traveltimes_bp %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, time_diff_percent), - 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::finishAndShowPdf(pdff) - - htmlwidgets::saveWidget( - widget = plotly::ggplotly(tt_bp_total), - file = "traveltimes_all.html" - ) + dplyr::left_join( + scenario_base_median[, c("month_id", "time_diff_base")] %>% + dplyr::mutate(percentiles = 0.5) + ) %>% + dplyr::mutate( + time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base + ) - 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 - - # Save HTML widgets ------------------------------------------------------------ - htmlwidgets::saveWidget( - widget = plotly::ggplotly(tt_bp_total), - title = "Boxplot: median traveltime total", - file = "boxplot_traveltimes-median_total.html" - ) + # Plotting - htmlwidgets::saveWidget( - plotly::ggplotly(tt_bp_quarter), - title = "Boxplot: median traveltime by quarter", - file = "boxplot_traveltimes-median_quarter.html" - ) + y_lim <- c(0,350) + + tt_bp_percent <- traveltimes_bp %>% + ggplot2::ggplot(ggplot2::aes( + x = forcats::fct_reorder(scenario_name, time_diff_percent), + 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::finishAndShowPdf(pdff) + + htmlwidgets::saveWidget( + widget = plotly::ggplotly(tt_bp_total), + file = "traveltimes_all.html" + ) + + 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 + + # 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" + ) } -### check model + ### check model model_dir <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" - atm_files <- fs::dir_ls(model_dir, recurse = TRUE, type = "file", regexp = "tracer.*ATMOSPH.IN") + atm_files <- fs::dir_ls( + model_dir, + recurse = TRUE, + type = "file", + regexp = "tracer.*ATMOSPH.IN" + ) atm_df <- lapply(atm_files, function(file) { atm <- kwb.hydrus1d::read_atmosph(file) - tibble::tibble(path = file, - Prec_sum = sum(atm$data$Prec, na.rm = TRUE)) - }) %>% dplyr::bind_rows() + tibble::tibble( + path = file, + Prec_sum = sum(atm$data$Prec, na.rm = TRUE) + ) + }) %>% + dplyr::bind_rows() - path_split <- stringr::str_split_fixed(atm_df$path,"/", 13) + path_split <- stringr::str_split_fixed(atm_df$path, "/", 13) atm_df %>% - dplyr::mutate(irrig_period = path_split[,9], - duration_extreme = path_split[,10], - scenario= path_split[,12] %>% stringr::str_remove("_soil-column_.*")) %>% + dplyr::mutate( + irrig_period = path_split[, 9], + duration_extreme = path_split[, 10], + scenario = path_split[, 12] %>% + stringr::str_remove("_soil-column_.*") + ) %>% dplyr::group_by(irrig_period, duration_extreme, scenario) %>% dplyr::summarise(Prec_mean = sum(Prec_sum)/dplyr::n()) %>% View() - - } From 5e61a4c1432658b21065a150b1790cd3a3bb3545 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:28:38 +0200 Subject: [PATCH 03/29] Loop over names instead of indices --- R/.scenarios_parallel.R | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index ddf835c..127057b 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -992,7 +992,6 @@ if (FALSE) # Set up parallel plan system.time(expr = { future::plan(future::multisession) - traveltimes_list <- future.apply::future_lapply(model_paths, function(model_path) { setNames(nm = (scenarios), future.apply::future_lapply(scenarios, function(scenario) { try({ @@ -1006,30 +1005,29 @@ if (FALSE) }) })) }) - } - ) + }) # Close the parallel plan future::plan(future::sequential) - traveltimes_list_backup <- traveltimes_list - - sapply(seq_along(traveltimes_list), function(i) { + traveltimes_list_backup <- traveltimes_list + sapply(names(traveltimes_list), function(name) { htmlwidgets::saveWidget( - flextreat.hydrus1d::plot_traveltimes( - traveltimes_list[[i]] %>% - dplyr::bind_rows(), - title = sprintf("%s", extrahiere_letzte_drei_teile(names(traveltimes_list)[i])), - ylim = c(0, 650) - ), - file = sprintf("traveltimes_%s.html", names(traveltimes_list)[i])) + traveltimes_list[[name]] %>% + dplyr::bind_rows() %>% + flextreat.hydrus1d::plot_traveltimes( + title = sprintf("%s", extrahiere_letzte_drei_teile(name)), + ylim = c(0, 650) + ), + file = sprintf("traveltimes_%s.html", name)) }) # extrahiere_letzte_drei_teile ------------------------------------------------- extrahiere_letzte_drei_teile <- function(pfad) { sapply(pfad, function(pf) { - # Teile den Pfad anhand des Schrägstrichs auf + + # Teile den Pfad anhand des Schraegstrichs auf teile <- unlist(strsplit(pf, "/")) # Waehle die letzten drei Teile aus @@ -1104,7 +1102,6 @@ if (FALSE) print(tt_bp_total) - kwb.utils::finishAndShowPdf(pdff) pdff <- "traveltimes_all.pdf" From 9983cb0556944b93947610080d7a753ec3948951 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:32:31 +0200 Subject: [PATCH 04/29] Move function definition out of if (FALSE) --- R/.scenarios_parallel.R | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 127057b..3983c20 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -1022,21 +1022,6 @@ if (FALSE) file = sprintf("traveltimes_%s.html", name)) }) - # extrahiere_letzte_drei_teile ------------------------------------------------- - extrahiere_letzte_drei_teile <- function(pfad) - { - sapply(pfad, function(pf) { - - # Teile den Pfad anhand des Schraegstrichs auf - teile <- unlist(strsplit(pf, "/")) - - # Waehle die letzten drei Teile aus - letzte_drei_teile <- tail(teile, 3) - - paste0(letzte_drei_teile, collapse = "_") - }) - } - # Plotting if (FALSE) { @@ -1345,3 +1330,18 @@ if (FALSE) dplyr::summarise(Prec_mean = sum(Prec_sum)/dplyr::n()) %>% View() } + +# extrahiere_letzte_drei_teile ------------------------------------------------- +extrahiere_letzte_drei_teile <- function(pfad) +{ + sapply(pfad, function(pf) { + + # Teile den Pfad anhand des Schraegstrichs auf + teile <- unlist(strsplit(pf, "/")) + + # Waehle die letzten drei Teile aus + letzte_drei_teile <- tail(teile, 3) + + paste0(letzte_drei_teile, collapse = "_") + }) +} From b3c59c09fd3fee02d747184eb50b88e155fe88da Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:39:26 +0200 Subject: [PATCH 05/29] Simplify and generalise extrahiere_letzte_[drei_]teile() --- R/.scenarios_parallel.R | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 3983c20..fd91dda 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -1016,7 +1016,7 @@ if (FALSE) traveltimes_list[[name]] %>% dplyr::bind_rows() %>% flextreat.hydrus1d::plot_traveltimes( - title = sprintf("%s", extrahiere_letzte_drei_teile(name)), + title = sprintf("%s", extrahiere_letzte_teile(name)), ylim = c(0, 650) ), file = sprintf("traveltimes_%s.html", name)) @@ -1031,7 +1031,7 @@ if (FALSE) # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows() # - # label <- extrahiere_letzte_drei_teile(path) + # label <- extrahiere_letzte_teile(path) # traveltime_bp <- traveltimes_sel %>% # dplyr::bind_rows() %>% @@ -1067,7 +1067,7 @@ if (FALSE) y_lim <- c(0, 350) tt_bp_total <- traveltime_bp %>% - dplyr::mutate(scenario_short = extrahiere_letzte_drei_teile(scenario)) %>% + dplyr::mutate(scenario_short = extrahiere_letzte_teile(scenario)) %>% ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, median), y = time_diff)) + ggplot2::geom_boxplot(outliers = FALSE) + ggplot2::geom_jitter( @@ -1098,7 +1098,7 @@ if (FALSE) 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(), + extrahiere_letzte_teile(), quarter = lubridate::quarter(date) %>% as.factor(), soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% @@ -1152,7 +1152,7 @@ if (FALSE) 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(), + extrahiere_letzte_teile(), quarter = lubridate::quarter(date) %>% as.factor(), soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% @@ -1331,17 +1331,11 @@ if (FALSE) View() } -# extrahiere_letzte_drei_teile ------------------------------------------------- -extrahiere_letzte_drei_teile <- function(pfad) +# extrahiere_letzte_teile ------------------------------------------------------ +extrahiere_letzte_teile <- function(pfad, n = 3L) { - sapply(pfad, function(pf) { - - # Teile den Pfad anhand des Schraegstrichs auf - teile <- unlist(strsplit(pf, "/")) - - # Waehle die letzten drei Teile aus - letzte_drei_teile <- tail(teile, 3) - - paste0(letzte_drei_teile, collapse = "_") - }) + sapply( + strsplit(pfad, "/"), + function(x) paste0(tail(x, n), collapse = "_") + ) } From 729431039d364117386c12a81a9ef58079d98cca Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:46:39 +0200 Subject: [PATCH 06/29] Do not next "if (FALSE)" --- R/.scenarios_parallel.R | 530 ++++++++++++++++++++-------------------- 1 file changed, 265 insertions(+), 265 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index fd91dda..c721f2b 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -905,7 +905,7 @@ if (FALSE) }) scenario_dirs <- fs::dir_ls( - path = root_path, + path = root_path, recurse = TRUE, regexp = "retardation_no/.*hydrus_scenarios.xlsx$", type = "file" @@ -1021,288 +1021,288 @@ if (FALSE) ), file = sprintf("traveltimes_%s.html", name)) }) +} - # Plotting - if (FALSE) - { - pdff <- "traveltimes_per-scenario.pdf" - kwb.utils::preparePdf(pdff) - # sapply(names(traveltimes_list), function(path) { - - # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows() - # - # label <- extrahiere_letzte_teile(path) - - # traveltime_bp <- traveltimes_sel %>% - # dplyr::bind_rows() %>% - # 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()) - - traveltime_bp <- lapply(traveltimes_list, function(x) { - x %>% - dplyr::bind_rows() %>% - 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() - ) +# Plotting --------------------------------------------------------------------- +if (FALSE) +{ + pdff <- "traveltimes_per-scenario.pdf" + kwb.utils::preparePdf(pdff) + # sapply(names(traveltimes_list), function(path) { + + # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows() + # + # label <- extrahiere_letzte_teile(path) + + # traveltime_bp <- traveltimes_sel %>% + # dplyr::bind_rows() %>% + # 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()) + + traveltime_bp <- lapply(traveltimes_list, function(x) { + x %>% + dplyr::bind_rows() %>% + 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 %>% - dplyr::mutate(scenario_short = extrahiere_letzte_teile(scenario)) %>% - ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, 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 = "Scenario", - title = "Boxplot: median traveltime total" - ) + - ggplot2::theme_bw() + - ggplot2::theme( - axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1) - ) + scenario_by_median_traveltime <- traveltime_bp %>% + dplyr::group_by(scenario) %>% + dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>% + dplyr::arrange(median) - print(tt_bp_total) - - kwb.utils::finishAndShowPdf(pdff) - - pdff <- "traveltimes_all.pdf" - kwb.utils::preparePdf(pdff) - - traveltimes_df <- lapply(traveltimes_list, 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_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) + traveltime_bp <- traveltime_bp %>% + dplyr::left_join(scenario_by_median_traveltime) - htmlwidgets::saveWidget( - widget = plotly::ggplotly(tt_bp_total), - file = "traveltimes_all.html" + y_lim <- c(0, 350) + + tt_bp_total <- traveltime_bp %>% + dplyr::mutate(scenario_short = extrahiere_letzte_teile(scenario)) %>% + ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, 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 = "Scenario", + title = "Boxplot: median traveltime total" + ) + + ggplot2::theme_bw() + + ggplot2::theme( + axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1) ) - pdff <- "traveltimes_all_percent.pdf" - kwb.utils::preparePdf(pdff) - - traveltimes_df <- lapply(traveltimes_list, 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_teile(), - quarter = lubridate::quarter(date) %>% - as.factor(), - soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% - stringr::str_remove_all("soil-|m") %>% as.factor() - ) + print(tt_bp_total) - 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) - - traveltimes_bp <- traveltimes_df %>% - dplyr::filter(percentiles == 0.5) %>% - dplyr::left_join( - scenario_base_median[, c("month_id", "time_diff_base")] %>% - dplyr::mutate(percentiles = 0.5) - ) %>% - dplyr::mutate( - time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base - ) + kwb.utils::finishAndShowPdf(pdff) - # Plotting - - y_lim <- c(0,350) - - tt_bp_percent <- traveltimes_bp %>% - ggplot2::ggplot(ggplot2::aes( - x = forcats::fct_reorder(scenario_name, time_diff_percent), - 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::finishAndShowPdf(pdff) + pdff <- "traveltimes_all.pdf" + kwb.utils::preparePdf(pdff) - htmlwidgets::saveWidget( - widget = plotly::ggplotly(tt_bp_total), - file = "traveltimes_all.html" - ) + traveltimes_df <- lapply(traveltimes_list, 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_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") - 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 - - # Save HTML widgets - htmlwidgets::saveWidget( - widget = plotly::ggplotly(tt_bp_total), - title = "Boxplot: median traveltime total", - file = "boxplot_traveltimes-median_total.html" + print(tt_bp_total) + + kwb.utils::finishAndShowPdf(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, 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_teile(), + quarter = lubridate::quarter(date) %>% + as.factor(), + soil_depth = stringr::str_extract(scenario_name, "soil-.*m") %>% + stringr::str_remove_all("soil-|m") %>% as.factor() ) - htmlwidgets::saveWidget( - plotly::ggplotly(tt_bp_quarter), - title = "Boxplot: median traveltime by quarter", - file = "boxplot_traveltimes-median_quarter.html" + 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) + + traveltimes_bp <- traveltimes_df %>% + dplyr::filter(percentiles == 0.5) %>% + dplyr::left_join( + scenario_base_median[, c("month_id", "time_diff_base")] %>% + dplyr::mutate(percentiles = 0.5) + ) %>% + dplyr::mutate( + time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base ) - } + # Plotting + + y_lim <- c(0,350) + + tt_bp_percent <- traveltimes_bp %>% + ggplot2::ggplot(ggplot2::aes( + x = forcats::fct_reorder(scenario_name, time_diff_percent), + 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") - ### check model + print(tt_bp_percent) - model_dir <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" + kwb.utils::finishAndShowPdf(pdff) + htmlwidgets::saveWidget( + widget = plotly::ggplotly(tt_bp_total), + file = "traveltimes_all.html" + ) + + 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 + + # 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" + ) + +} + +# Check model ------------------------------------------------------------------ +if (FALSE) +{ atm_files <- fs::dir_ls( - model_dir, + path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed", recurse = TRUE, type = "file", regexp = "tracer.*ATMOSPH.IN" From cb60bc016ad584366a86a0ecd4b382a5163cec72 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 07:57:21 +0200 Subject: [PATCH 07/29] Clean section "Read and plot solute" --- R/.scenarios_parallel.R | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index c721f2b..dbb6160 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -687,25 +687,32 @@ if (FALSE) { parallel::stopCluster(cl) } -# After main loop -------------------------------------------------------------- +# Read and plot solute --------------------------------------------------------- if (FALSE) { - solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>% - dplyr::mutate(difftime = c(0,diff(time))) + read_solute_with_difftime <- function(file) { + kwb.hydrus1d::read_solute(file) %>% + dplyr::mutate(difftime = c(0, diff(time))) + } - max(solute$sum_cv_top) + min(solute$sum_cv_bot) + min(solute$cv_ch1) + plot_solute <- function(solute) { + plot(solute$time, solute$c_top) + points(solute$c_bot, col = "red") + } - plot(solute$time, solute$c_top) - points(solute$c_bot, col = "red") + solute <- read_solute_with_difftime(paths$solute_vs) + with(solute, max(sum_cv_top) + min(sum_cv_bot) + min(cv_ch1)) + plot_solute(solute) (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))) - + solute <- read_solute_with_difftime(paths$solute) (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 +} +# Rest ------------------------------------------------------------------------- +if (FALSE) +{ sum(atmos$data$Prec) max(tlevel$sum_infil) max(tlevel$sum_evap) From 85b03bce0d5aa144ced1170f0378e24012ec9fe5 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 08:10:32 +0200 Subject: [PATCH 08/29] Add and clean section "Read and plot balance" --- R/.scenarios_parallel.R | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index dbb6160..99e6518 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -708,36 +708,28 @@ if (FALSE) paths$solute_vs2 <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out" solute <- read_solute_with_difftime(paths$solute) (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 -} -# Rest ------------------------------------------------------------------------- -if (FALSE) -{ sum(atmos$data$Prec) max(tlevel$sum_infil) max(tlevel$sum_evap) +} - balance <- kwb.hydrus1d::read_balance(paths$balance) - - 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() - - balance %>% - dplyr::filter( - subdomain_id == 0, - time < 400, - parameter == "CncBalR" - ) %>% - ggplot2::ggplot(mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))) + +# Read and plot balance -------------------------------------------------------- +if (FALSE) +{ + kwb.hydrus1d::read_balance(paths$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() +} +# Rest ------------------------------------------------------------------------- +if (FALSE) +{ sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) solute_day <- flextreat.hydrus1d::aggregate_solute(solute) From 0f7c0a6cac8777756310e46eeb89adde3fe0ce53 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 08:23:52 +0200 Subject: [PATCH 09/29] Move things to where they belong --- R/.scenarios_parallel.R | 41 ++++++++++++++++++++--------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 99e6518..879b1f3 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -712,24 +712,7 @@ if (FALSE) sum(atmos$data$Prec) max(tlevel$sum_infil) max(tlevel$sum_evap) -} -# Read and plot balance -------------------------------------------------------- -if (FALSE) -{ - kwb.hydrus1d::read_balance(paths$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() -} - -# Rest ------------------------------------------------------------------------- -if (FALSE) -{ sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) solute_day <- flextreat.hydrus1d::aggregate_solute(solute) @@ -755,7 +738,26 @@ if (FALSE) sum(solute$cv_top[condition]) solute_aggr_date <- flextreat.hydrus1d::aggregate_solute(solute) + solute_aggr_date + View(solute_aggr_date) +} +# Read and plot balance -------------------------------------------------------- +if (FALSE) +{ + kwb.hydrus1d::read_balance(paths$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() +} + +# Rest ------------------------------------------------------------------------- +if (FALSE) +{ obsnode <- kwb.hydrus1d::read_obsnode(paths$obs_node) obsnode %>% @@ -778,15 +780,12 @@ if (FALSE) plotly::ggplotly(p) - solute_aggr_date - View(tlevel_aggr_date) - View(solute_aggr_date) - t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) t_level ## t_level aggregate tlevel_aggr_date <- flextreat.hydrus1d::aggregate_tlevel(t_level) + View(tlevel_aggr_date) tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel( t_level, From 3beca4512abe9337364b7258d3d3e7c2b8f39576 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 08:31:21 +0200 Subject: [PATCH 10/29] Add/clean more sections --- R/.scenarios_parallel.R | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 879b1f3..0c7db34 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -755,7 +755,7 @@ if (FALSE) ggplot2::geom_point() } -# Rest ------------------------------------------------------------------------- +# Read and plot obsnode -------------------------------------------------------- if (FALSE) { obsnode <- kwb.hydrus1d::read_obsnode(paths$obs_node) @@ -765,10 +765,10 @@ if (FALSE) dplyr::group_by(node_id) %>% dplyr::summarise(sum = sum(value)) - profile <- kwb.hydrus1d::read_profile(paths$profile) - p <- obsnode %>% - dplyr::left_join(profile[,c("node_id", "x")]) %>% + dplyr::left_join( + kwb.hydrus1d::read_profile(paths$profile)[, c("node_id", "x")] + ) %>% dplyr::filter(variable == "conc1") %>% ggplot2::ggplot(mapping = ggplot2::aes( x = time, @@ -779,7 +779,11 @@ if (FALSE) ggplot2::theme_bw() plotly::ggplotly(p) +} +# Read and aggregate tlevel ---------------------------------------------------- +if (FALSE) +{ t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) t_level @@ -797,15 +801,24 @@ if (FALSE) col_aggr = "year_hydrologic" ) %>% dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october +} - 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) +# Plot water balance ----------------------------------------------------------- +if (FALSE) +{ + plots <- list( + tlevel_aggr_date = tlevel_aggr_date, + tlevel_aggr_yearmonth = tlevel_aggr_yearmonth, + tlevel_aggr_year_hydrologic = tlevel_aggr_year_hydrologic + ) %>% + flextreat.hydrus1d::plot_waterbalance - wb_date_plot - wb_yearmonth_plot - wb_yearhydrologic_plot + print(plots) +} +# Rest ------------------------------------------------------------------------- +if (FALSE) +{ solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))] paths$model_dir_vs From dc35e525f41c5dac2488538f0490f7f29ddb31ce Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 08:33:28 +0200 Subject: [PATCH 11/29] Move tlevel expressions after tlevel assignment --- R/.scenarios_parallel.R | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 0c7db34..9eeb011 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -781,7 +781,7 @@ if (FALSE) plotly::ggplotly(p) } -# Read and aggregate tlevel ---------------------------------------------------- +# Read and aggregate tlevel, plot water balance -------------------------------- if (FALSE) { t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) @@ -801,11 +801,8 @@ if (FALSE) col_aggr = "year_hydrologic" ) %>% dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october -} -# Plot water balance ----------------------------------------------------------- -if (FALSE) -{ + # Plot water balance plots <- list( tlevel_aggr_date = tlevel_aggr_date, tlevel_aggr_yearmonth = tlevel_aggr_yearmonth, From 2d035d935fe46a81535a392b69e884cf1cd349e2 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sat, 12 Oct 2024 08:44:37 +0200 Subject: [PATCH 12/29] Move expressions to where they seem to belong --- R/.scenarios_parallel.R | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 9eeb011..fd8859c 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -593,6 +593,9 @@ inner_function <- function(config, atm_data, soil_columns, helper) atm_default <- atmos tlevel <- kwb.hydrus1d::read_tlevel(paths$t_level) + max(tlevel$sum_infil) + max(tlevel$sum_evap) + sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) vs_atm <- flextreat.hydrus1d::recalculate_ctop_with_virtualstorage( atm = atm_default$data, @@ -709,11 +712,9 @@ if (FALSE) solute <- read_solute_with_difftime(paths$solute) (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 - sum(atmos$data$Prec) - max(tlevel$sum_infil) - max(tlevel$sum_evap) + solute$time[min(which.max(solute$sum_cv_top))] - sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) + sum(atmos$data$Prec) solute_day <- flextreat.hydrus1d::aggregate_solute(solute) @@ -816,8 +817,6 @@ if (FALSE) # Rest ------------------------------------------------------------------------- if (FALSE) { - solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))] - paths$model_dir_vs #scenarios_solutes <- paste0(scenarios, "_soil-column") From 6d26794153afd6aac9445be165a357e9dcde22d8 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:02:59 +0200 Subject: [PATCH 13/29] Get rid of duplication in "dry" and "wet" case --- R/.scenarios_parallel.R | 53 +++++++++++------------------------------ 1 file changed, 14 insertions(+), 39 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index fd8859c..8d1870e 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -35,42 +35,19 @@ get_atm <- function(atm, extreme_rain = NULL) 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)) - - 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) - - return(atm_wet_input) - } + atm_dry_wet <- atm_stats %>% + dplyr::filter(rain_mm == (if (extreme_rain == "dry") min else max)(rain_mm)) + + atm_dry_wet_sel <- atm %>% + dplyr::filter(year == atm_dry_wet$year) %>% + dplyr::select(day_of_year, rain_mm) + + atm %>% + dplyr::select(- rain_mm, - hydrologic_year, - year) %>% + dplyr::left_join(atm_dry_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) } # prepare_solute_input --------------------------------------------------------- @@ -817,8 +794,6 @@ if (FALSE) # Rest ------------------------------------------------------------------------- if (FALSE) { - paths$model_dir_vs - #scenarios_solutes <- paste0(scenarios, "_soil-column") scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median") @@ -826,7 +801,7 @@ if (FALSE) root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" scenario_dirs <- fs::dir_ls( - path = root_path, + path = root_path, recurse = TRUE, regexp = "retardation.*vs$", type = "directory" From c459abc659f946dc13200771a479983a06bc524e Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:19:03 +0200 Subject: [PATCH 14/29] Integrate kd() into prepare_solute_input() --- R/.scenarios_parallel.R | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 8d1870e..11ff65e 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -56,22 +56,35 @@ prepare_solute_input <- function( selector, Ks = NULL, SnkL1 = NULL, - diff_w = 0, diff_g = 0, - kd = NULL, + diff_w = 0, + diff_g = 0, halftime_to_firstorderrate = NULL ) { `%>%` <- magrittr::`%>%` + # helper function to calculate kd + # https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html + kd <- function(porosity, retardation, bulk_density) { + (retardation - 1) * porosity / bulk_density + } + stopifnot(nrow(dat) <= 10) if (selector$solute$No.Solutes != nrow(dat)) { selector$solute$No.Solutes <- nrow(dat) } - solute_names <- sprintf("solute_%d", seq_len(nrow(dat))) + indices <- seq_len(nrow(dat)) + solute_names <- sprintf("solute_%d", indices) + named_indices <- setNames(indices, solute_names) - solutes_new <- setNames(nm = solute_names, lapply(seq_len(nrow(dat)), function(i) { + cols <- c( + "Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'", + "SnkS1'", "SnkG1'", "SnkL0", "SnkS0", "SnkG0", "Alfa" + ) + + solutes_new <- lapply(named_indices, function(i) { dat_sel <- dat[i,] @@ -81,11 +94,6 @@ prepare_solute_input <- function( bulk_density = selector$solute$transport$Bulk.d. ) - 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() @@ -110,7 +118,7 @@ prepare_solute_input <- function( reaction = reaction ) - })) + }) sel_tmp <- selector$solute[!names(selector$solute) %in% solute_names] @@ -134,13 +142,6 @@ get_mean <- function(col) round(rowMeans(x), digits = 2) } -# 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 --------------------------------------------------- halftime_to_firstorderrate <- function(half_time) { @@ -545,7 +546,6 @@ inner_function <- function(config, atm_data, soil_columns, helper) selector = selector, diff_w = 0, diff_g = 0, - kd = helper("kd"), halftime_to_firstorderrate = helper("halftime_to_firstorderrate") ) @@ -618,8 +618,7 @@ helper <- kwb.utils::createAccessor(list( prepare_solute_input = prepare_solute_input, halftime_to_firstorderrate = halftime_to_firstorderrate, get_valid_exe_path = get_valid_exe_path, - generate_periods = generate_periods, - kd = kd + generate_periods = generate_periods )) # Main loop -------------------------------------------------------------------- From e581ea98122e48c1a44c60ac2237bfdbac54d6f8 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:21:54 +0200 Subject: [PATCH 15/29] Integrate get_mean() into provide_soil_columns() --- R/.scenarios_parallel.R | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 11ff65e..5da2a03 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -63,7 +63,6 @@ prepare_solute_input <- function( { `%>%` <- magrittr::`%>%` - # helper function to calculate kd # https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html kd <- function(porosity, retardation, bulk_density) { (retardation - 1) * porosity / bulk_density @@ -134,14 +133,6 @@ prepare_solute_input <- function( ) } -# get_mean --------------------------------------------------------------------- -get_mean <- function(col) -{ - x <- stringr::str_split_fixed(col, "-", n = 2) - mode(x) <- "numeric" - round(rowMeans(x), digits = 2) -} - # halftime_to_firstorderrate --------------------------------------------------- halftime_to_firstorderrate <- function(half_time) { @@ -158,6 +149,12 @@ provide_soil_columns <- function(path) { `%>%` <- magrittr::`%>%` + get_mean <- function(col) { + x <- stringr::str_split_fixed(col, "-", n = 2) + mode(x) <- "numeric" + round(rowMeans(x), digits = 2) + } + kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE) %>% janitor::clean_names() %>% dplyr::mutate(half_life_days = dplyr::case_when( From ee72f27363f9470597b3fa9f6e5384b702f40ba8 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:27:42 +0200 Subject: [PATCH 16/29] Integrate halftime_to_firstorderrate() into prepare_solute_input() --- R/.scenarios_parallel.R | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 5da2a03..f94a45e 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -57,8 +57,7 @@ prepare_solute_input <- function( Ks = NULL, SnkL1 = NULL, diff_w = 0, - diff_g = 0, - halftime_to_firstorderrate = NULL + diff_g = 0 ) { `%>%` <- magrittr::`%>%` @@ -68,6 +67,15 @@ prepare_solute_input <- function( (retardation - 1) * porosity / bulk_density } + # 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 + halftime_to_firstorderrate <- function(half_time) { + if (half_time != 0) { + 0.693 / half_time + } else { + 0 + } + } + stopifnot(nrow(dat) <= 10) if (selector$solute$No.Solutes != nrow(dat)) { @@ -133,17 +141,6 @@ prepare_solute_input <- function( ) } -# 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 - } else { - 0 - } -} - # provide_soil_columns --------------------------------------------------------- provide_soil_columns <- function(path) { @@ -542,8 +539,7 @@ inner_function <- function(config, atm_data, soil_columns, helper) SnkL1 = if (tracer) 0, # else NULL selector = selector, diff_w = 0, - diff_g = 0, - halftime_to_firstorderrate = helper("halftime_to_firstorderrate") + diff_g = 0 ) kwb.hydrus1d::write_selector(solutes_new, paths$selector) @@ -613,7 +609,6 @@ helper <- kwb.utils::createAccessor(list( 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, generate_periods = generate_periods )) From e1a551b516ac422a25eef9266edaebabbfbaf739 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:34:30 +0200 Subject: [PATCH 17/29] Use helper function index_of() --- R/.scenarios_parallel.R | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index f94a45e..cb35756 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -128,16 +128,15 @@ prepare_solute_input <- function( }) 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)] - ) + index_of <- function(name) which(names(sel_tmp) == name) c( selector[names(selector) != "solute"], - list(solute = solutes_new_list) + list(solute = c( + sel_tmp[1:index_of("transport")], + solutes_new, + sel_tmp[index_of("kTopSolute"):length(sel_tmp)] + )) ) } From 758effddccf1769714df26d3bb22ab89ef45b91e Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:36:46 +0200 Subject: [PATCH 18/29] Move MAIN section to the top of the script and put it in "if (FALSE)" --- R/.scenarios_parallel.R | 53 +++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index cb35756..b8a90af 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -8,6 +8,33 @@ if (FALSE) { library(magrittr) library(flextreat.hydrus1d) +# MAIN ------------------------------------------------------------------------- +if (FALSE) +{ + #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 <- 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 %>% + 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("tracer", "ka", "o3"), # "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 = c(TRUE, FALSE), + duration_string = "long", #c("short", "long"), + retardation_scenario = c("retardation_yes", "retardation_no", "tracer") + ) +} + # get_atm ---------------------------------------------------------------------- get_atm <- function(atm, extreme_rain = NULL) { @@ -175,32 +202,6 @@ provide_soil_columns <- function(path) 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 <- 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 %>% - 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("tracer", "ka", "o3"), # "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 = c(TRUE, FALSE), - duration_string = "long", #c("short", "long"), - retardation_scenario = c("retardation_yes", "retardation_no", "tracer") -) - - # generate_solute_ids ---------------------------------------------------------- generate_solute_ids <- function(n) { From 97e68f03d18ae448d0c5cf20fb6d69fe65d47993 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:49:53 +0200 Subject: [PATCH 19/29] Move provide_soil_columns() up --- R/.scenarios_parallel.R | 70 ++++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index b8a90af..d75c4fd 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -35,6 +35,41 @@ if (FALSE) ) } +# provide_soil_columns --------------------------------------------------------- +provide_soil_columns <- function(path) +{ + `%>%` <- magrittr::`%>%` + + get_mean <- function(col) { + x <- stringr::str_split_fixed(col, "-", n = 2) + mode(x) <- "numeric" + round(rowMeans(x), digits = 2) + } + + 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) +} + # get_atm ---------------------------------------------------------------------- get_atm <- function(atm, extreme_rain = NULL) { @@ -167,41 +202,6 @@ prepare_solute_input <- function( ) } -# provide_soil_columns --------------------------------------------------------- -provide_soil_columns <- function(path) -{ - `%>%` <- magrittr::`%>%` - - get_mean <- function(col) { - x <- stringr::str_split_fixed(col, "-", n = 2) - mode(x) <- "numeric" - round(rowMeans(x), digits = 2) - } - - 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) -} - # generate_solute_ids ---------------------------------------------------------- generate_solute_ids <- function(n) { From d45c3a1b36374c735fc11ddbef31df9a3239bdb7 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 07:50:31 +0200 Subject: [PATCH 20/29] Extract period definition --- R/.scenarios_parallel.R | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index d75c4fd..4065bd9 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -419,13 +419,14 @@ inner_function <- function(config, atm_data, soil_columns, helper) atm[which(!lubridate::month(atm$date) %in% 4:9), IRRIGATION_COLUMNS] <- 0 } - if (config$duration_string == "test") { - atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2018-04-30") - } else 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") - } + periods <- list( + test = c("2017-05-01", "2018-04-30"), + short = c("2017-05-01", "2020-04-30"), + long = c("2017-05-01", "2023-12-31") + ) + + period <- kwb.utils::selectElements(periods, config$duration_string) + atm <- dplyr::filter(atm, date >= period[1L] & date <= period[2L]) days_monthly <- lubridate::days_in_month( seq.Date(from = min(atm$date), to = max(atm$date), by = "month") From f0177904085e78d46f1d25600ff44a81cf1fa15c Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 08:18:27 +0200 Subject: [PATCH 21/29] Use kwb.utils::startsToRanges() and compare with what generate_solute_ids() returns --- R/.scenarios_parallel.R | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 4065bd9..0274eb8 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -217,10 +217,26 @@ generate_solute_ids <- function(n) seq_end <- c(seq_end, n) } - tibble::tibble( + result <- tibble::tibble( start = as.integer(seq_start), end = as.integer(seq_end) ) + + ranges <- kwb.utils::startsToRanges( + starts = seq.int(from = 1L, to = n, by = 10L), + startOffset = 0L, + stopOffset = 1L, + lastStop = as.integer(n) + ) + + result_2 <- tibble::tibble( + start = ranges$from, + end = ranges$to + ) + + stopifnot(identical(result, result_2)) + + result_2 } # generate_periods ------------------------------------------------------------- From 9831b6adc9d7021406001a3f8deb03c4298fe991 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 08:19:35 +0200 Subject: [PATCH 22/29] Keep only the new version --- R/.scenarios_parallel.R | 23 +---------------------- 1 file changed, 1 insertion(+), 22 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 0274eb8..8dc32b4 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -205,23 +205,6 @@ prepare_solute_input <- function( # 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 - } - - if (length(seq_end) < length(seq_start)) { - seq_end <- c(seq_end, n) - } - - result <- tibble::tibble( - start = as.integer(seq_start), - end = as.integer(seq_end) - ) - ranges <- kwb.utils::startsToRanges( starts = seq.int(from = 1L, to = n, by = 10L), startOffset = 0L, @@ -229,14 +212,10 @@ generate_solute_ids <- function(n) lastStop = as.integer(n) ) - result_2 <- tibble::tibble( + tibble::tibble( start = ranges$from, end = ranges$to ) - - stopifnot(identical(result, result_2)) - - result_2 } # generate_periods ------------------------------------------------------------- From 272ad784e30a4e55a21d82b03200490c2e71aa6b Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 08:33:40 +0200 Subject: [PATCH 23/29] Rename generate_solute_ids() to generate_index_ranges() and reused it. Delete generate_periods() (not needed any more) --- R/.scenarios_parallel.R | 34 +++++++++++----------------------- 1 file changed, 11 insertions(+), 23 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 8dc32b4..659dbf9 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -202,8 +202,8 @@ prepare_solute_input <- function( ) } -# generate_solute_ids ---------------------------------------------------------- -generate_solute_ids <- function(n) +# generate_index_ranges -------------------------------------------------------- +generate_index_ranges <- function(n) { ranges <- kwb.utils::startsToRanges( starts = seq.int(from = 1L, to = n, by = 10L), @@ -218,19 +218,6 @@ generate_solute_ids <- function(n) ) } -# 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) - } - ) -} - # prepare_files_for_irrig_int -------------------------------------------------- prepare_files_for_irrig_int <- function(paths) { @@ -427,11 +414,13 @@ inner_function <- function(config, atm_data, soil_columns, helper) seq.Date(from = min(atm$date), to = max(atm$date), by = "month") ) - loop_df <- if (tracer) { - helper("generate_periods")(n = length(days_monthly)) - } else { - helper("generate_solute_ids")(n = nrow(soil_columns)) - } + loop_df <- helper("generate_index_ranges")( + n = if (tracer) { + length(days_monthly) + } else { + nrow(soil_columns) + } + ) } kwb.utils::catAndRun( @@ -600,13 +589,12 @@ inner_function <- function(config, atm_data, soil_columns, helper) # helper ----------------------------------------------------------------------- helper <- kwb.utils::createAccessor(list( get_atm = get_atm, - generate_solute_ids = generate_solute_ids, + generate_index_ranges = generate_index_ranges, 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, - get_valid_exe_path = get_valid_exe_path, - generate_periods = generate_periods + get_valid_exe_path = get_valid_exe_path )) # Main loop -------------------------------------------------------------------- From 49ee82b19513d37d749484e833bfdff47ed57f8f Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 10:55:00 +0200 Subject: [PATCH 24/29] Move main loop up into MAIN section --- R/.scenarios_parallel.R | 85 +++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 45 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 659dbf9..5600944 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -33,6 +33,46 @@ if (FALSE) duration_string = "long", #c("short", "long"), retardation_scenario = c("retardation_yes", "retardation_no", "tracer") ) + + atm_data <- flextreat.hydrus1d::prepare_atmosphere_data() + + + #arg_combis <- arg_combis[arg_combis$scenario == "soil-3m_irrig-10days" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$duration_string == "long" & arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$extreme_rain == "", ] + #arg_combis <- arg_combis[arg_combis$treatment != "tracer" & arg_combis$retardation_scenario != "tracer" | arg_combis$treatment == "tracer" & arg_combis$retardation_scenario == "tracer", ] + #arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario == "soil-1m_irrig-10days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$extreme_rain == "",] + arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario != "soil-1m_irrig-01days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == TRUE & arg_combis$extreme_rain == "",] + + + configs <- lapply(seq_len(nrow(arg_combis)), function(i) { + as.list(arg_combis[i, ]) + }) + + # Sequential loop + for (config in configs) { + #config <- configs[[1L]] + inner_function( + config = config, + atm_data = atm_data, + soil_columns = soil_columns_selected, + helper = helper + ) + } + + # Parallel loop + ncores <- parallel::detectCores() + ncores <- 8 + cl <- parallel::makeCluster(ncores) + + parallel::parLapply( + cl = cl, + X = configs, + fun = inner_function, + atm_data = atm_data, + soil_columns = soil_columns_selected, + helper = helper + ) + + parallel::stopCluster(cl) } # provide_soil_columns --------------------------------------------------------- @@ -597,51 +637,6 @@ helper <- kwb.utils::createAccessor(list( get_valid_exe_path = get_valid_exe_path )) -# Main loop -------------------------------------------------------------------- - -if (FALSE) { - - atm_data <- flextreat.hydrus1d::prepare_atmosphere_data() - - - #arg_combis <- arg_combis[arg_combis$scenario == "soil-3m_irrig-10days" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$duration_string == "long" & arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$extreme_rain == "", ] - #arg_combis <- arg_combis[arg_combis$treatment != "tracer" & arg_combis$retardation_scenario != "tracer" | arg_combis$treatment == "tracer" & arg_combis$retardation_scenario == "tracer", ] - #arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario == "soil-1m_irrig-10days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$extreme_rain == "",] - arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario != "soil-1m_irrig-01days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == TRUE & arg_combis$extreme_rain == "",] - - - configs <- lapply(seq_len(nrow(arg_combis)), function(i) { - as.list(arg_combis[i, ]) - }) - - # Sequential loop - for (config in configs) { - #config <- configs[[1L]] - inner_function( - config = config, - atm_data = atm_data, - soil_columns = soil_columns_selected, - helper = helper - ) - } - - # Parallel loop - ncores <- parallel::detectCores() - ncores <- 8 - cl <- parallel::makeCluster(ncores) - - parallel::parLapply( - cl = cl, - X = configs, - fun = inner_function, - atm_data = atm_data, - soil_columns = soil_columns_selected, - helper = helper - ) - - parallel::stopCluster(cl) -} - # Read and plot solute --------------------------------------------------------- if (FALSE) { From ceabdb9fcddfb20051d7e67cbe4d9e0633c6d805 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 12:33:36 +0200 Subject: [PATCH 25/29] Provide directory structure and files locally --- R/.scenarios_parallel.R | 84 ++++++++++++++++++++++++++++------------- 1 file changed, 57 insertions(+), 27 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 5600944..1813740 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -8,17 +8,38 @@ if (FALSE) { library(magrittr) library(flextreat.hydrus1d) +# MAIN: Provide directory structure and required files locally ----------------- +if (FALSE) +{ + grammar <- provide_paths() + + if (!file.exists(grammar$AP_3_1_4_local)) { + kwb.utils::copyDirectoryStructure( + sourcedir = grammar$AP_3_1_4_server, + targetdir = grammar$AP_3_1_4_local + ) + } + + paths_server <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_server") + paths_local <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_local") + + file.copy(paths_server$file_soil_columns, paths_local$file_soil_columns) + paths$file_substance_classes +} + # MAIN ------------------------------------------------------------------------- if (FALSE) { - #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" + grammar <- provide_paths() - soil_columns <- provide_soil_columns(path) + #paths <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_server") + paths <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_local") - ### 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 <- provide_soil_columns(file = paths$file_soil_columns) + selected_substances <- readr::read_csv(file = paths$file_substance_classes) + + ### Select 1 substance for 5 different half life classes defined in this table soil_columns_selected <- soil_columns %>% dplyr::filter(substanz_nr %in% selected_substances$substance_id) %>% dplyr::mutate(id = 1:dplyr::n()) @@ -36,13 +57,11 @@ if (FALSE) atm_data <- flextreat.hydrus1d::prepare_atmosphere_data() - #arg_combis <- arg_combis[arg_combis$scenario == "soil-3m_irrig-10days" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$duration_string == "long" & arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$extreme_rain == "", ] #arg_combis <- arg_combis[arg_combis$treatment != "tracer" & arg_combis$retardation_scenario != "tracer" | arg_combis$treatment == "tracer" & arg_combis$retardation_scenario == "tracer", ] #arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario == "soil-1m_irrig-10days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == FALSE & arg_combis$extreme_rain == "",] arg_combis <- arg_combis[arg_combis$retardation_scenario == "tracer" & arg_combis$treatment == "tracer" & arg_combis$scenario != "soil-1m_irrig-01days" & arg_combis$duration_string == "long" & arg_combis$irrig_only_growing_season == TRUE & arg_combis$extreme_rain == "",] - configs <- lapply(seq_len(nrow(arg_combis)), function(i) { as.list(arg_combis[i, ]) }) @@ -76,7 +95,7 @@ if (FALSE) } # provide_soil_columns --------------------------------------------------------- -provide_soil_columns <- function(path) +provide_soil_columns <- function(file) { `%>%` <- magrittr::`%>%` @@ -86,7 +105,7 @@ provide_soil_columns <- function(path) round(rowMeans(x), digits = 2) } - kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE) %>% + kwb.db::hsGetTable(file, "my_results2", stringsAsFactors = FALSE) %>% janitor::clean_names() %>% dplyr::mutate(half_life_days = dplyr::case_when( grepl(">", hwz_tage) ~ hwz_tage %>% @@ -351,23 +370,24 @@ get_valid_exe_path <- function(exe_dir) } # provide_paths ---------------------------------------------------------------- -provide_paths <- function(config, start, end) +provide_paths <- function(config = NULL, 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 - tracer <- config$treatment == "tracer" - # Define a path grammar + #\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no + PATH_GRAMMAR <- list( - exe_dir = sprintf("C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed///%s", config$retardation_scenario), + AP_3_1_4_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell", + AP_3_1_4_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell", + USER = Sys.getenv("USERNAME"), + file_soil_columns = "/StofflicheModellrandbedingungen.xlsx", + file_substance_classes = "inst/extdata/input-data/substance_classes.csv", + exe_dir = "/Vivian/Rohdaten/irrig_fixed///", model_name_org = "model_to_copy", model_name = "__soil-column_", ###model_gui_path_org = "/.h1d", - model_gui_path_org = "/.h1d", - + model_gui_path_org = "/.h1d", model_gui_path = "/.h1d", modelvs_gui_path = "/_vs.h1d", - model_dir_org = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", - model_dir = "/", model_dir_vs = "/_vs", atmosphere = "/ATMOSPH.IN", @@ -384,19 +404,20 @@ provide_paths <- function(config, start, end) 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" + soil_data = "/input-data/soil/soil_geolog.csv" ) + if (is.null(config)) { + return(PATH_GRAMMAR) + } + + tracer <- config$treatment == "tracer" + # Define a path grammar + # Resolve the path grammar by replacing the placeholders recursively kwb.utils::resolve( PATH_GRAMMAR, + AP_3_1_4 = "AP_3_1_4_server", irrig_dir_string = ifelse( config$irrig_only_growing_season, "irrig-period_growing-season", @@ -409,7 +430,16 @@ provide_paths <- function(config, start, end) "" }, #final_subdir = ifelse(tracer, "tracer", config$retardation_scenario), - scenario = config$scenario + scenario = config$scenario, + retardation_scenario = config$retardation_scenario, + solute_id_start = sprintf("%02d", start), + solute_id_end = sprintf("%02d", end), + location = if (tracer) { + "tracer" + } else { + # e.g. "ablauf_ka_median" + sprintf("ablauf_%s_median", config$treatment) + } ) } From e9dc4906f6bc574471138dee396f3754cfba99c0 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 12:35:43 +0200 Subject: [PATCH 26/29] Move provide_paths() up --- R/.scenarios_parallel.R | 148 ++++++++++++++++++++-------------------- 1 file changed, 74 insertions(+), 74 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 1813740..4347d13 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -94,6 +94,80 @@ if (FALSE) parallel::stopCluster(cl) } +# provide_paths ---------------------------------------------------------------- +provide_paths <- function(config = NULL, start = "", end = "") +{ + #\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no + + PATH_GRAMMAR <- list( + AP_3_1_4_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell", + AP_3_1_4_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell", + USER = Sys.getenv("USERNAME"), + file_soil_columns = "/StofflicheModellrandbedingungen.xlsx", + file_substance_classes = "inst/extdata/input-data/substance_classes.csv", + exe_dir = "/Vivian/Rohdaten/irrig_fixed///", + 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 = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/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" + ) + + if (is.null(config)) { + return(PATH_GRAMMAR) + } + + tracer <- config$treatment == "tracer" + # Define a path grammar + + # Resolve the path grammar by replacing the placeholders recursively + kwb.utils::resolve( + PATH_GRAMMAR, + AP_3_1_4 = "AP_3_1_4_server", + 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, + retardation_scenario = config$retardation_scenario, + solute_id_start = sprintf("%02d", start), + solute_id_end = sprintf("%02d", end), + location = if (tracer) { + "tracer" + } else { + # e.g. "ablauf_ka_median" + sprintf("ablauf_%s_median", config$treatment) + } + ) +} + # provide_soil_columns --------------------------------------------------------- provide_soil_columns <- function(file) { @@ -369,80 +443,6 @@ get_valid_exe_path <- function(exe_dir) } } -# provide_paths ---------------------------------------------------------------- -provide_paths <- function(config = NULL, start = "", end = "") -{ - #\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no - - PATH_GRAMMAR <- list( - AP_3_1_4_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell", - AP_3_1_4_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell", - USER = Sys.getenv("USERNAME"), - file_soil_columns = "/StofflicheModellrandbedingungen.xlsx", - file_substance_classes = "inst/extdata/input-data/substance_classes.csv", - exe_dir = "/Vivian/Rohdaten/irrig_fixed///", - 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 = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/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" - ) - - if (is.null(config)) { - return(PATH_GRAMMAR) - } - - tracer <- config$treatment == "tracer" - # Define a path grammar - - # Resolve the path grammar by replacing the placeholders recursively - kwb.utils::resolve( - PATH_GRAMMAR, - AP_3_1_4 = "AP_3_1_4_server", - 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, - retardation_scenario = config$retardation_scenario, - solute_id_start = sprintf("%02d", start), - solute_id_end = sprintf("%02d", end), - location = if (tracer) { - "tracer" - } else { - # e.g. "ablauf_ka_median" - sprintf("ablauf_%s_median", config$treatment) - } - ) -} - # inner_function --------------------------------------------------------------- inner_function <- function(config, atm_data, soil_columns, helper) { From 662f10714d54d57fd2b9ab2ba9696399dbd30434 Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 12:41:47 +0200 Subject: [PATCH 27/29] Replace path constant with placeholder TODO: replace the placeholder (back) in a cleaner way --- R/.scenarios_parallel.R | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 4347d13..8c66df8 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -24,7 +24,6 @@ if (FALSE) paths_local <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_local") file.copy(paths_server$file_soil_columns, paths_local$file_soil_columns) - paths$file_substance_classes } # MAIN ------------------------------------------------------------------------- @@ -112,7 +111,7 @@ provide_paths <- function(config = NULL, start = "", end = "") model_gui_path_org = "/.h1d", model_gui_path = "/.h1d", modelvs_gui_path = "/_vs.h1d", - model_dir_org = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/", + model_dir_org = "/Vivian/Rohdaten/irrig_fixed/", model_dir = "/", model_dir_vs = "/_vs", atmosphere = "/ATMOSPH.IN", @@ -685,7 +684,7 @@ if (FALSE) plot_solute(solute) (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" + paths$solute_vs2 <- "/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out" solute <- read_solute_with_difftime(paths$solute) (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100 @@ -798,7 +797,7 @@ if (FALSE) scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median") #root_path <- "D:/hydrus1d/irrig_fixed_01" - root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" + root_path <- "/Vivian/Rohdaten/irrig_fixed" scenario_dirs <- fs::dir_ls( path = root_path, @@ -901,7 +900,7 @@ if (FALSE) ) # load_default <- res_stats$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/retardation_no/ablauf_ka_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx` - name <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed/irrig-period_status-quo/long/retardation_no/ablauf_o3_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx" + name <- "/Vivian/Rohdaten/irrig_fixed/irrig-period_status-quo/long/retardation_no/ablauf_o3_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx" load_default <- res_stats[[name]] %>% # dplyr::mutate( # retardation = basename(dirname(dirname(path))), @@ -931,7 +930,7 @@ if (FALSE) View(res_stats_df) #root_path <- "D:/hydrus1d/irrig_fixed_01" - root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed" + root_path <- "/Vivian/Rohdaten/irrig_fixed" model_paths <- fs::dir_ls( path = root_path, @@ -1284,7 +1283,7 @@ if (FALSE) if (FALSE) { atm_files <- fs::dir_ls( - path = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed", + path = "/Vivian/Rohdaten/irrig_fixed", recurse = TRUE, type = "file", regexp = "tracer.*ATMOSPH.IN" From 1cddc08a89379d54eaf5b7d20c2e667387fd43fd Mon Sep 17 00:00:00 2001 From: hsonne Date: Sun, 13 Oct 2024 13:01:02 +0200 Subject: [PATCH 28/29] Return early, use inline assignment --- R/.scenarios_parallel.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 8c66df8..65e46c9 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -435,11 +435,11 @@ sum_per_interval <- function(data, interval) # 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() + if (file.exists(file <- file.path(exe_dir, "H1D_CALC.exe"))) { + return(file) } + + kwb.hydrus1d::check_hydrus_exe() } # inner_function --------------------------------------------------------------- From e7b2062cd6fc7d6bf4b3e444ab679c3a0081c06f Mon Sep 17 00:00:00 2001 From: hsonne Date: Mon, 14 Oct 2024 01:03:26 +0200 Subject: [PATCH 29/29] WIP --- R/.scenarios_parallel.R | 638 +++++++++++++++++++++------------------- 1 file changed, 329 insertions(+), 309 deletions(-) diff --git a/R/.scenarios_parallel.R b/R/.scenarios_parallel.R index 65e46c9..3de797d 100644 --- a/R/.scenarios_parallel.R +++ b/R/.scenarios_parallel.R @@ -35,11 +35,10 @@ if (FALSE) paths <- kwb.utils::resolve(grammar, AP_3_1_4 = "AP_3_1_4_local") soil_columns <- provide_soil_columns(file = paths$file_soil_columns) - selected_substances <- readr::read_csv(file = paths$file_substance_classes) ### Select 1 substance for 5 different half life classes defined in this table - 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()) @@ -96,8 +95,7 @@ if (FALSE) # provide_paths ---------------------------------------------------------------- provide_paths <- function(config = NULL, start = "", end = "") { - #\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no - + # Define a path grammar PATH_GRAMMAR <- list( AP_3_1_4_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell", AP_3_1_4_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell", @@ -107,7 +105,7 @@ provide_paths <- function(config = NULL, start = "", end = "") exe_dir = "/Vivian/Rohdaten/irrig_fixed///", model_name_org = "model_to_copy", model_name = "__soil-column_", - ###model_gui_path_org = "/.h1d", + ### model_gui_path_org = "/.h1d", model_gui_path_org = "/.h1d", model_gui_path = "/.h1d", modelvs_gui_path = "/_vs.h1d", @@ -136,7 +134,6 @@ provide_paths <- function(config = NULL, start = "", end = "") } tracer <- config$treatment == "tracer" - # Define a path grammar # Resolve the path grammar by replacing the placeholders recursively kwb.utils::resolve( @@ -153,15 +150,14 @@ provide_paths <- function(config = NULL, start = "", end = "") } else { "" }, - #final_subdir = ifelse(tracer, "tracer", config$retardation_scenario), scenario = config$scenario, retardation_scenario = config$retardation_scenario, solute_id_start = sprintf("%02d", start), solute_id_end = sprintf("%02d", end), + # e.g. "ablauf_ka_median" location = if (tracer) { "tracer" } else { - # e.g. "ablauf_ka_median" sprintf("ablauf_%s_median", config$treatment) } ) @@ -244,192 +240,69 @@ get_atm <- function(atm, extreme_rain = NULL) dplyr::relocate(rain_mm, .after = clearwater.mmPerDay) } -# prepare_solute_input --------------------------------------------------------- -prepare_solute_input <- function( - dat, - selector, - Ks = NULL, - SnkL1 = NULL, - diff_w = 0, - diff_g = 0 -) -{ - `%>%` <- magrittr::`%>%` - - # https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html - kd <- function(porosity, retardation, bulk_density) { - (retardation - 1) * porosity / bulk_density - } - - # 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 - halftime_to_firstorderrate <- function(half_time) { - if (half_time != 0) { - 0.693 / half_time - } else { - 0 - } - } - - stopifnot(nrow(dat) <= 10) - - if (selector$solute$No.Solutes != nrow(dat)) { - selector$solute$No.Solutes <- nrow(dat) - } - - indices <- seq_len(nrow(dat)) - solute_names <- sprintf("solute_%d", indices) - named_indices <- setNames(indices, solute_names) - - cols <- c( - "Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'", - "SnkS1'", "SnkG1'", "SnkL0", "SnkS0", "SnkG0", "Alfa" - ) - - solutes_new <- lapply(named_indices, 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. - ) - - reaction <- matrix(data = 0, ncol = length(cols), nrow = length(ks)) %>% - 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)) { - round(halftime_to_firstorderrate(dat_sel$half_life_days), 5) - } else { - SnkL1 - } - - list( - diffusion = tibble::tibble(DifW = diff_w, DifG = diff_g), - reaction = reaction - ) - - }) - - sel_tmp <- selector$solute[!names(selector$solute) %in% solute_names] - index_of <- function(name) which(names(sel_tmp) == name) - - c( - selector[names(selector) != "solute"], - list(solute = c( - sel_tmp[1:index_of("transport")], - solutes_new, - sel_tmp[index_of("kTopSolute"):length(sel_tmp)] - )) - ) -} - -# generate_index_ranges -------------------------------------------------------- -generate_index_ranges <- function(n) +# extract_irrig_interval ------------------------------------------------------- +# Return the string that is used as "irrig_interval" +extract_irrig_interval <- function(model_dir) { - ranges <- kwb.utils::startsToRanges( - starts = seq.int(from = 1L, to = n, by = 10L), - startOffset = 0L, - stopOffset = 1L, - lastStop = as.integer(n) - ) - - tibble::tibble( - start = ranges$from, - end = ranges$to + string <- stringr::str_extract(model_dir, "[0-9][0-9]?days") + paste( + as.integer(stringr::str_extract(string, "\\d+")), + stringr::str_extract(string, "[a-z]+") ) } -# prepare_files_for_irrig_int -------------------------------------------------- -prepare_files_for_irrig_int <- function(paths) +# copy_model_files ------------------------------------------------------------- +copy_model_files <- function( + from_model_dir, + to_model_dir, + from_model_file, + to_model_file, + overwrite = FALSE +) { - `%>%` <- magrittr::`%>%` - - p <- kwb.utils::createAccessor(paths) - - copy <- function(fun, from, to) { + copy <- function(fun, from, to, ...) { kwb.utils::catAndRun( sprintf("Copying from\n %s to\n %s", from, to), - fun(path = from, new_path = to) + fun(path = from, new_path = to, ...) ) } - if (!fs::dir_exists(p("model_dir"))) { - copy(fs::dir_copy, from = p("model_dir_org"), to = p("model_dir")) + if (!fs::dir_exists(to_model_dir) || overwrite) { + copy(fs::dir_copy, from_model_dir, to_model_dir, overwrite = overwrite) } - if (!fs::file_exists(p("model_gui_path"))) { - copy(fun = fs::file_copy, from = p("model_gui_path_org"), to = p("model_gui_path")) + if (!fs::file_exists(to_model_file) || overwrite) { + copy(fs::file_copy, from_model_file, to_model_file, overwrite = overwrite) } +} - model_out_files <- list.files( - p("model_dir"), - pattern = "\\.out$", - full.names = TRUE - ) +# delete_model_output_files ---------------------------------------------------- +delete_model_output_files <- function(model_dir) +{ + files <- list.files(model_dir, pattern = "\\.out$", full.names = TRUE) - if (length(model_out_files) > 0L) { - fs::file_delete(model_out_files) + if (length(files) > 0L) { + fs::file_delete(files) } +} - soil_depth_cm <- 100 * - stringr::str_extract(p("model_name"), "soil-[0-9]+?m") %>% +# get_soil_depth_cm ------------------------------------------------------------ +get_soil_depth_cm <- function(model_name) +{ + `%>%` <- magrittr::`%>%` + 100 * + stringr::str_extract(model_name, "soil-[0-9]+?m") %>% + # @mrustl: Why not "[0-9]+" ? 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 = p("profile")) - } - - string_irrig_int <- stringr::str_extract(p("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]+") - ) } -# sum_per_interval ------------------------------------------------------------- -sum_per_interval <- function(data, interval) +# update_soil_profile ---------------------------------------------------------- +update_soil_profile <- function(file, x_end) { - `%>%` <- magrittr::`%>%` - - 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 + p <- kwb.hydrus1d::read_profile(file) + p$profile <- kwb.hydrus1d::extend_soil_profile(p$profile, x_end = x_end) + kwb.hydrus1d::write_profile(p, file) } # get_valid_exe_path ----------------------------------------------------------- @@ -447,51 +320,70 @@ inner_function <- function(config, atm_data, soil_columns, helper) { `%>%` <- magrittr::`%>%` - { - # Define constants - IRRIGATION_COLUMNS <- c("groundwater.mmPerDay", "clearwater.mmPerDay") + # Define helper function get_c_tops() + get_c_tops <- function(days_monthly, n) { + days_total <- cumsum(days_monthly) + indices <- seq_along(days_total) + dplyr::bind_cols(lapply(indices, function(i) { + x <- rep(0, n) + 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 == indices[1]) { + "cTop" + } else { + sprintf("cTop%d", which(indices == i)) + } + tib + })) + } - # Provide variables from config - extreme_rain <- if (config$extreme_rain == "") { - NULL - } else { - config$extreme_rain - } - tracer <- config$treatment == "tracer" - retardation <- config$retardation == "retardation_yes" + # Provide variables from config + extreme_rain <- if (config$extreme_rain == "") { + NULL + } else { + config$extreme_rain + } - if (retardation) { - soil_columns$retard <- 1 - } + tracer <- config$treatment == "tracer" + retardation <- config$retardation == "retardation_yes" - atm <- helper("get_atm")(atm_data, extreme_rain) + if (retardation) { + soil_columns$retard <- 1 + } - if (config$irrig_only_growing_season) { - atm[which(!lubridate::month(atm$date) %in% 4:9), IRRIGATION_COLUMNS] <- 0 - } + atm <- helper("get_atm")(atm_data, extreme_rain) - periods <- list( - test = c("2017-05-01", "2018-04-30"), - short = c("2017-05-01", "2020-04-30"), - long = c("2017-05-01", "2023-12-31") - ) + if (config$irrig_only_growing_season) { + atm <- reset_irrigation(atm, skip_months = 4:9) + } - period <- kwb.utils::selectElements(periods, config$duration_string) - atm <- dplyr::filter(atm, date >= period[1L] & date <= period[2L]) + periods <- list( + test = c("2017-05-01", "2018-04-30"), + short = c("2017-05-01", "2020-04-30"), + long = c("2017-05-01", "2023-12-31") + ) - days_monthly <- lubridate::days_in_month( - seq.Date(from = min(atm$date), to = max(atm$date), by = "month") - ) + period <- kwb.utils::selectElements(periods, config$duration_string) + atm <- dplyr::filter(atm, date >= period[1L] & date <= period[2L]) - loop_df <- helper("generate_index_ranges")( - n = if (tracer) { - length(days_monthly) - } else { - nrow(soil_columns) - } - ) + days_monthly <- lubridate::days_in_month( + seq.Date(from = min(atm$date), to = max(atm$date), by = "month") + ) + + n <- if (tracer) { + length(days_monthly) + } else { + nrow(soil_columns) } + loop_df <- kwb.utils::startsToRanges( + starts = seq.int(from = 1L, to = n, by = 10L), + startOffset = 0L, + stopOffset = 1L, + lastStop = as.integer(n) + ) + kwb.utils::catAndRun( messageText = sprintf( "Running '%s' period for treatment '%s'", @@ -502,106 +394,61 @@ inner_function <- function(config, atm_data, soil_columns, helper) sapply(seq_len(nrow(loop_df)), function(i) { #i <- 1L - { - paths <- helper("provide_paths")( - config, - start = loop_df$start[i], - end = loop_df$end[i] - ) - - 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_int <- stringr::str_detect(paths$model_dir, "irrig-[0-9][0-9]?days") - } + solute_ids <- loop_df$from[i]:loop_df$to[i] - if (irrig_int) { - irrig_interval <- helper("prepare_files_for_irrig_int")(paths) - } + paths <- helper("provide_paths")( + config, + start = solute_ids[1L], + end = solute_ids[length(solute_ids)] + ) - # no-irrigation - if (no_irrig) { - atm[, IRRIGATION_COLUMNS] <- 0 - } + irrig_int <- stringr::str_detect(paths$model_dir, "irrig-[0-9][0-9]?days") + irrig_interval <- extract_irrig_interval(paths$model_dir) if (irrig_int) { - atm <- helper("sum_per_interval")(data = atm, interval = irrig_interval) - } - - days_total <- cumsum(days_monthly) - indices <- seq_along(days_total) + # @mrustl: Why all these steps only in case of "irrig_int"? - c_tops <- lapply(indices, function(i) { - - x <- rep(0, nrow(atm)) - x_min <- ifelse(i == 1, 1, days_total[i - 1] + 1) - x[x_min:days_total[i]] <- rep(100, days_monthly[i]) + copy_model_files( + from_model_dir = paths$model_dir_org, + to_model_dir = paths$model_dir, + from_model_file = paths$model_gui_path_org, + to_model_file = paths$model_gui_path + ) - tib <- data.frame(x) + delete_model_output_files(paths$model_dir) - colnames(tib) <- if (i == indices[1]) { - "cTop" - } else { - sprintf("cTop%d", which(indices %in% i)) + if (depth <- get_soil_depth_cm(paths$model_name) != 200) { + update_soil_profile(file = paths$profile, x_end = -depth) } - tib - }) %>% dplyr::bind_cols() - - 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 - ) - } 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 - ) } - n_tsteps <- nrow(atm_prep) - - writeLines( - kwb.hydrus1d::write_atmosphere(atm = atm_prep), - paths$atmosphere + # Prepare the model run by updating model files + update_atmosphere( + atm = atm, + no_irrig = stringr::str_detect(paths$model_dir, "no-irrig"), + irrig_int = irrig_int, + irrig_interval = irrig_interval, + tracer = tracer, + c_tops_tracer = get_c_tops(days_monthly, n = nrow(atm))[, solute_ids], + c_tops_no_tracer = soil_columns[solute_ids, paths$location], + file = paths$atmosphere ) - selector <- kwb.hydrus1d::read_selector(path = paths$selector) - - selector$time$tMax <- n_tsteps - selector$time$MPL <- 250 - - selector$time$TPrint <- seq(0, n_tsteps, n_tsteps/selector$time$MPL) - - if (selector$time$TPrint[1] == 0) { - selector$time$TPrint[1] <- 1 - } - - 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 + update_selector( + file = paths$selector, + n_tsteps = n_tsteps, + dat = soil_columns[solute_ids, ], + tracer = tracer ) - 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) + update_hydrus1d( + file = paths$hydrus1d, + n_solutes = length(solute_ids) + ) + # Run the model exe_path <- helper("get_valid_exe_path")(paths$exe_dir) kwb.hydrus1d::run_model( @@ -610,39 +457,35 @@ inner_function <- function(config, atm_data, soil_columns, helper) print_output = FALSE ) + # Prepare next model run atmos <- kwb.hydrus1d::read_atmosph(paths$atmosphere) - - #atmos$data[names(c_tops)] <- c_tops - - atm_default <- atmos - tlevel <- kwb.hydrus1d::read_tlevel(paths$t_level) + max(tlevel$sum_infil) max(tlevel$sum_evap) sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top)) - vs_atm <- flextreat.hydrus1d::recalculate_ctop_with_virtualstorage( - atm = atm_default$data, + vs_atmos_data <- flextreat.hydrus1d::recalculate_ctop_with_virtualstorage( + atm = atmos$data, tlevel = tlevel, crit_v_top = - 0.05 ) - atmos$data[names(vs_atm$df)] <- vs_atm$df - - fs::dir_copy(paths$model_dir, paths$model_dir_vs, overwrite = TRUE) - fs::file_copy(paths$model_gui_path, paths$modelvs_gui_path, overwrite = TRUE) + atmos$data[names(vs_atmos_data)] <- vs_atmos_data - model_vs_out_files <- list.files(paths$model_dir_vs, pattern = "\\.out$", full.names = TRUE) + copy_model_files( + from_model_dir = paths$model_dir, + to_model_dir = paths$model_dir_vs, + from_model_file = paths$model_gui_path, + to_model_file = paths$modelvs_gui_path, + overwrite = TRUE + ) - if (length(model_vs_out_files) > 0) { - fs::file_delete(model_vs_out_files) - } + delete_model_output_files(paths$model_dir_vs) - writeLines( - kwb.hydrus1d::write_atmosphere(atm = atmos$data), - paths$atmosphere_vs - ) + write_atmosphere_file(atm = atmos$data, file = paths$atmosphere_vs) + # Run the model again kwb.hydrus1d::run_model( exe_path = exe_path, model_path = paths$model_dir_vs, @@ -655,6 +498,183 @@ inner_function <- function(config, atm_data, soil_columns, helper) ) } +# reset_irrigation ------------------------------------------------------------- +reset_irrigation <- function(atm, skip_months = integer()) +{ + i <- if (length(skip_months) > 0L) { + which(!lubridate::month(atm$date) %in% skip_months) + } else { + seq_len(nrow(atm)) + } + atm[i, c("groundwater.mmPerDay", "clearwater.mmPerDay")] <- 0 + atm +} + +# update_atmosphere ------------------------------------------------------------ +update_atmosphere <- function( + atm, no_irrig, irrig_int, irrig_interval, + tracer, c_tops_tracer, c_tops_no_tracer, + file +) +{ + # Define helper function + sum_per_interval <- function(data, interval) { + `%>%` <- magrittr::`%>%` + columns <- c("groundwater.mmPerDay", "clearwater.mmPerDay") + data_summed <- data %>% + dplyr::select(tidyselect::all_of(c("date", columns))) + # Konvertiere date in datetime-Format + dplyr::mutate(group = lubridate::floor_date(date, unit = interval)) %>% + # Gruppiere nach Zeitintervallen + dplyr::group_by(group) %>% + # Berechne die Summe fuer jedes Intervall + dplyr::summarise_at(.vars = tidyselect::all_of(columns), .funs = sum) %>% + dplyr::rename(date = group) + data[, columns] <- 0 + data[data$date %in% data_summed$date, columns] <- data_summed[, columns] + data + } + + # no-irrigation + if (no_irrig) { + atm <- reset_irrigation(atm) + } + + if (irrig_int) { + atm <- sum_per_interval(data = atm, interval = irrig_interval) + } + + conc_default <- if (tracer) { + c_tops_tracer + } else { + 0 + } + + atm_prep <- flextreat.hydrus1d::prepare_atmosphere( + atm = atm, + conc_irrig_clearwater = if (tracer) { + c_tops_tracer + } else { + c_tops_no_tracer + }, + conc_irrig_groundwater = conc_default, + conc_rain = conc_default + ) + + n_tsteps <- nrow(atm_prep) + + write_atmosphere_file(atm = atm_prep, file) +} + +# write_atmosphere_file -------------------------------------------------------- +write_atmosphere_file <- function(atm, file) +{ + writeLines(kwb.hydrus1d::write_atmosphere(atm), file) +} + +# update_selector -------------------------------------------------------------- +update_selector <- function(file, n_tsteps, dat, tracer) +{ + # Define helper functions + + # https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html + kd <- function(porosity, retardation, bulk_density) { + (retardation - 1) * porosity / bulk_density + } + + # 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 + halftime_to_firstorderrate <- function(half_time) { + if (half_time != 0) { + 0.693 / half_time + } else { + 0 + } + } + + prepare_solute_input <- function( + dat, selector, Ks = NULL, SnkL1 = NULL, diff_w = 0, diff_g = 0 + ) { + `%>%` <- magrittr::`%>%` + stopifnot(nrow(dat) <= 10) + if (selector$solute$No.Solutes != nrow(dat)) { + selector$solute$No.Solutes <- nrow(dat) + } + indices <- seq_len(nrow(dat)) + solute_names <- sprintf("solute_%d", indices) + named_indices <- setNames(indices, solute_names) + cols <- c( + "Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'", + "SnkS1'", "SnkG1'", "SnkL0", "SnkS0", "SnkG0", "Alfa" + ) + solutes_new <- lapply(named_indices, 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. + ) + reaction <- matrix(data = 0, ncol = length(cols), nrow = length(ks)) %>% + 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)) { + round(halftime_to_firstorderrate(dat_sel$half_life_days), 5) + } else { + SnkL1 + } + list( + diffusion = tibble::tibble(DifW = diff_w, DifG = diff_g), + reaction = reaction + ) + }) + sel_tmp <- selector$solute[!names(selector$solute) %in% solute_names] + index_of <- function(name) which(names(sel_tmp) == name) + c( + selector[names(selector) != "solute"], + list(solute = c( + sel_tmp[1:index_of("transport")], + solutes_new, + sel_tmp[index_of("kTopSolute"):length(sel_tmp)] + )) + ) + } + + selector <- kwb.hydrus1d::read_selector(path = file) + + selector$time$tMax <- n_tsteps + selector$time$MPL <- 250 + selector$time$TPrint <- seq(0, n_tsteps, n_tsteps/selector$time$MPL) + + if (selector$time$TPrint[1] == 0) { + selector$time$TPrint[1] <- 1 + } + + solutes_new <- prepare_solute_input( + dat = dat, + Ks = if (tracer) 0, # else NULL + SnkL1 = if (tracer) 0, # else NULL + selector = selector, + diff_w = 0, + diff_g = 0 + ) + + kwb.hydrus1d::write_selector(solutes_new, file) +} + +# update_hydrus1d -------------------------------------------------------------- +update_hydrus1d <- function(file, n_solutes) +{ + h <- kwb.hydrus1d::read_hydrus1d(file) + h$Main$NumberOfSolutes <- n_solutes + kwb.hydrus1d::write_hydrus1d(h, file) +} + # helper ----------------------------------------------------------------------- helper <- kwb.utils::createAccessor(list( get_atm = get_atm,