Skip to content

Commit

Permalink
Progress on conservative substance travel times
Browse files Browse the repository at this point in the history
  • Loading branch information
@mrustl
mrustl committed Jul 3, 2024
1 parent 32c0130 commit 231b98a
Showing 1 changed file with 182 additions and 7 deletions.
189 changes: 182 additions & 7 deletions R/.virtual_storage.R
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
remotes::install_github("kwb-r/kwb.hydrus1d@dev")
remotes::install_github("kwb-r/flextreat.hydrus1d@dev")

paths_list <- list(
#extdata = system.file("extdata", package = "flextreat.hydrus1d"),
Expand All @@ -7,7 +8,7 @@ paths_list <- list(
#root_local = "C:/kwb/projects/flextreat/hydrus/Szenarien_10day",
#root_local = system.file("extdata/model", package = "flextreat.hydrus1d"),
exe_dir = "<root_local>",
model_name = "test_fracht1", #"1a2a_BTA_korr_test_40d",
model_name = "1a2a_tracer_3140", #"1a2a_BTA_korr_test_40d",
model_gui_path = "<exe_dir>/<model_name>.h1d",
modelvs_gui_path = "<exe_dir>/<model_name>_vs.h1d",
model_dir = "<exe_dir>/<model_name>",
Expand All @@ -30,15 +31,73 @@ paths_list <- list(


paths <- kwb.utils::resolve(paths_list)
paths$solute

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)


library(flextreat.hydrus1d)
atm <- flextreat.hydrus1d::prepare_atmosphere_data()
atm_selected <- flextreat.hydrus1d::select_hydrologic_years(atm)
# atm_prep <- flextreat.hydrus1d::prepare_atmosphere(atm = atm_selected,
# conc_irrig_clearwater = c(6738,
# 875,
# 4291,
# 2884,
# 1062),
# conc_irrig_groundwater = 0,
# conc_rain = 0
# )
atm <- atm_selected
days_monthy <- lubridate::days_in_month(seq.Date(from = min(atm$date),
to = max(atm$date),
by = "month"))

days_total <- cumsum(days_monthy)

indeces <- 31:40

c_tops <- lapply(indeces, function(i) {

x <- rep(0, nrow(atm))
if(i == 1) {
x_min = 1
} else {
x_min = days_total[i - 1] + 1
}
x[x_min:days_total[i]] <- rep(100, days_monthy[i])


tib <- data.frame(x)
colnames(tib) <- if(i == indeces[1]) {
"cTop"} else {
sprintf("cTop%d", which(indeces %in% i))
}

tib
}) %>% dplyr::bind_cols()




atm_prep <- flextreat.hydrus1d::prepare_atmosphere(atm = atm_selected,
conc_irrig_clearwater = c_tops,
conc_irrig_groundwater = 0,
conc_rain = 0
)


writeLines(kwb.hydrus1d::write_atmosphere(atm = atm_prep),
paths$atmosphere)

kwb.hydrus1d::run_model(model_path = paths$model_dir)

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)
Expand All @@ -49,24 +108,28 @@ vs_atm <- flextreat.hydrus1d::recalculate_ctop_with_virtualstorage(
crit_v_top = - 0.05
)

atmos$data$cTop <- vs_atm$cTop

atmos$data[names(vs_atm$df)] <- vs_atm$df

writeLines(kwb.hydrus1d::write_atmosphere(atm = atmos$data),
paths$atmosphere_vs)


kwb.hydrus1d::run_model(model_path = paths$model_dir_vs)

solute <- kwb.hydrus1d::read_solute(paths$solute) %>%
solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>%
dplyr::mutate(difftime = c(0,diff(time)))

plot(solute$time, solute$c_top)
points(solute$c_bot, col = "red")
(1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop)) * 100

solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>%

paths$solute_vs2 <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out"

solute <- kwb.hydrus1d::read_solute(paths$solute) %>%
dplyr::mutate(difftime = c(0,diff(time)))

(1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop)) * 100
(1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100


sum(atmos$data$Prec)
Expand Down Expand Up @@ -139,6 +202,118 @@ p <- obsnode %>%

plotly::ggplotly(p)

ssolute_aggr_date
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)
tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel(t_level,
col_aggr = "yearmonth")
tlevel_aggr_year_hydrologic <- flextreat.hydrus1d::aggregate_tlevel(t_level,
col_aggr = "year_hydrologic") %>%
dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october


wb_date_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_date)
wb_yearmonth_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_yearmonth)
wb_yearhydrologic_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_year_hydrologic)

wb_date_plot
wb_yearmonth_plot
wb_yearhydrologic_plot


solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))]

paths$model_dir_vs

solute_files <- fs::dir_ls(paths$exe_dir,
regexp = "1a2a_tracer.*_vs/solute\\d\\d?.out",
recurse = TRUE)



sol_travel <- lapply(solute_files, function(path) {

solute <- kwb.hydrus1d::read_solute(path)

tibble::tibble(
model_name = basename(dirname(path)),
solute_name = kwb.utils::removeExtension(basename(path)),
p01_top = mean(solute$time[which(max(solute$sum_cv_top)*0.005 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.015)], na.rm = TRUE),
p01_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.015 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.005)], na.rm = TRUE),
p01_diff = p01_bot - p01_top,
p05_top = mean(solute$time[which(max(solute$sum_cv_top)*0.045 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.055)], na.rm = TRUE),
p05_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.055 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.045)], na.rm = TRUE),
p05_diff = p05_bot - p05_top,
p10_top = mean(solute$time[which(max(solute$sum_cv_top)*0.09 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.11)], na.rm = TRUE),
p10_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.11 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.09)], na.rm = TRUE),
p10_diff = p10_bot - p10_top,
p25_top = mean(solute$time[which(max(solute$sum_cv_top)*0.24 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.26)], na.rm = TRUE),
p25_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.26 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.24)], na.rm = TRUE),
p25_diff = p25_bot - p25_top,
p50_top = mean(solute$time[which(max(solute$sum_cv_top)*0.48 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.52)], na.rm = TRUE),
p50_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.52 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.48)], na.rm = TRUE),
p50_diff = p50_bot - p50_top,
p75_top = mean(solute$time[which(max(solute$sum_cv_top)*0.73 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.77)], na.rm = TRUE),
p75_bot =mean(solute$time[which(- max(solute$sum_cv_top)*0.77 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.73)], na.rm = TRUE),
p75_diff = p75_bot - p75_top,
p90_top = mean(solute$time[which(max(solute$sum_cv_top)*0.88 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.92)], na.rm = TRUE),
p90_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.92 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.88)], na.rm = TRUE),
p90_diff = p90_bot - p90_top,
p95_top = mean(solute$time[which(max(solute$sum_cv_top)*0.935 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*0.965)], na.rm = TRUE),
p95_bot = mean(solute$time[which(- max(solute$sum_cv_top)*0.965 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.935)], na.rm = TRUE),
p95_diff = p95_bot - p95_top,
p99_top = mean(solute$time[which(max(solute$sum_cv_top)*0.98 <= solute$sum_cv_top & solute$sum_cv_top <= max(solute$sum_cv_top)*1)], na.rm = TRUE),
p99_bot = mean(solute$time[which(- max(solute$sum_cv_top)*1 >= solute$sum_cv_bot & solute$sum_cv_bot <= - max(solute$sum_cv_top)*0.98)], na.rm = TRUE),
p99_diff = p99_bot - p99_top
)
}) %>% dplyr::bind_rows()


sol_travel <- sol_travel %>%
dplyr::mutate(solute_name = stringr::str_remove(solute_name, "solute") %>% as.integer()) %>%
dplyr::rename(solute_id = solute_name) %>%
dplyr::arrange(model_name, solute_id)

lookup_model <- sol_travel %>%
dplyr::group_by(model_name) %>%
dplyr::summarise(n = dplyr::n()) %>%
dplyr::mutate(n_cum = cumsum(n),
n_cum_fix = n_cum - dplyr::first(n))


library(lubridate)

get_last_day_of_months <- function(ids) {
sapply(ids, function(id) {
start_date <- lubridate::ymd("2017-05-01")
month_date <- start_date %m+% months(id - 1)
last_day <- lubridate::ceiling_date(month_date, "month") - days(1)

return(last_day)
}) %>% as.Date()
}


sol_travel_tot <- lookup_model %>%
dplyr::left_join(sol_travel) %>%
dplyr::mutate(month_id = solute_id + n_cum_fix,
date = get_last_day_of_months(month_id))


p1 <- sol_travel_tot %>%
dplyr::select(date, tidyselect::ends_with("diff")) %>%
tidyr::pivot_longer(- date) %>%
dplyr::mutate(name = stringr::str_remove(name, "_diff")) %>%
ggplot2::ggplot(ggplot2::aes(x = date, y = value, col = name)) +
ggplot2::geom_point() +
ggplot2::geom_line() +
ggplot2::theme_bw()

plotly::ggplotly(p1)

0 comments on commit 231b98a

Please sign in to comment.