Complete first scenarios 2017-2020 and prepare

prolonged simulation period 2017-05-01 - 2023-12-31

R/.ctops_conservative_tracer.R

@@ -0,0 +1,27 @@
+days_monthy <- lubridate::days_in_month(seq.Date(from = min(atm$date),
+                                                 to = max(atm$date),
+                                                 by = "month"))
+
+days_total <- cumsum(days_monthy)
+
+indeces <- 11:20
+
+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()

R/.scenarios_add-trace-organics_vs.R

@@ -60,8 +60,12 @@ kd <- function(porosity, retardation, bulk_density) {
 
 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
+  }
 }
 
 
@@ -82,8 +86,9 @@ soil_columns <- kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE
   dplyr::filter(!is.na(half_life_days)) %>%
   dplyr::mutate(id = 1:dplyr::n()) %>%
   dplyr::relocate(id) %>%
-  dplyr::mutate(retard = 1,
-                half_life_days = 0)
+  dplyr::mutate(retard = 1#,
+                #half_life_days = 0
+                )
 
 
 
@@ -97,7 +102,7 @@ if(length(seq_end) < length(seq_start)) seq_end <- c(seq_end, nrow(soil_columns)
 solute_ids <- tibble::tibble(start = seq_start,
                               end = seq_end)
 treatments <- c("ka", "o3") #"ka" #c("ka", "o3") #"ka"
-
+#treatments <- c("ka")
 sapply(treatments, function(treatment) {
 sapply(scenarios, function(scenario) {
 
@@ -106,7 +111,7 @@ sapply(seq_len(nrow(solute_ids)), function(i) {
   paths_list <- list(
     #extdata = system.file("extdata", package = "flextreat.hydrus1d"),
     #root_server = "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no",
-    root_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/conservative-transport",
+    root_local = "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no",
     #root_local = "C:/kwb/projects/flextreat/hydrus/Szenarien_10day",
     #root_local =  system.file("extdata/model", package = "flextreat.hydrus1d"),
     exe_dir = "<root_local>",
@@ -452,7 +457,7 @@ paths$model_dir_vs
 scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median")
 
 solutes_list <- setNames(lapply(scenarios_solutes, function(scenario) {
-  solute_files <- fs::dir_ls("C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/conservative-transport/",
+  solute_files <- fs::dir_ls("C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no/",
                            regexp = sprintf(".*%s.*_vs/solute\\d\\d?.out",
                                             scenario),
                            recurse = TRUE)
@@ -485,12 +490,13 @@ solutes_df <- solutes_list %>%
   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::mutate(mass_balance_error_percent = 100*(sum_cv_top + cv_ch1 + sum_cv_bot)/sum_cv_top,
+                soil = sum_cv_top - (cv_ch1 + sum_cv_bot)) %>%
   dplyr::arrange(mass_balance_error_percent)
 
-exp_dir <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/conservative-transport/"
-saveRDS(solutes_list, file = file.path(exp_dir, "solutes-list_conservative-transport.rds"))
+exp_dir <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/retardation_no/"
+saveRDS(solutes_list, file = file.path(exp_dir, "solutes-list_retardation_no.rds"))
 
-openxlsx::write.xlsx(solutes_df, file = file.path(exp_dir, "hydrus_scenarios_conservative-transport.xlsx"))
+openxlsx::write.xlsx(solutes_df, file = file.path(exp_dir, "hydrus_scenarios_retardation-no.xlsx"))
 
 View(solutes_list$`soil-1m_irrig-01days_soil-column`)

data-raw/climate_data.R

@@ -17,32 +17,41 @@ url_bs_rain <- rdwd::selectDWD(name = "Braunschweig",
 
 bs_rain <- rdwd::dataDWD(url_bs_rain)
 
-precipitation_hourly <- rdwd::dataDWD(url_bs_rain) %>%
-  dplyr::select(.data$MESS_DATUM, .data$R1) %>%
+precipitation_hourly <- bs_rain %>%
+  dplyr::select(MESS_DATUM, R1) %>%
   dplyr::rename("datetime" = "MESS_DATUM",
                 "precipitation_mm" = "R1")
 
-usethis::use_data(precipitation_hourly)
+usethis::use_data(precipitation_hourly, overwrite = TRUE)
 
 
 precipitation_daily <- precipitation_hourly %>%
   dplyr::mutate("date" = as.Date(datetime)) %>%
   dplyr::group_by(date) %>%
   dplyr::summarise(rain_mm = sum(precipitation_mm))
 
-usethis::use_data(precipitation_daily)
+usethis::use_data(precipitation_daily, overwrite = TRUE)
 
 
 ## DWD datasets
-remotes::install_github("kwb-r/kwb.dwd")
+remotes::install_github("kwb-r/kwb.dwd@get-rid-of-rgdal")
 
 
 shape_file <- system.file("extdata/input-data/gis/Abwasserverregnungsgebiet.shp",
                           package = "flextreat.hydrus1d")
 
 
-yearmonth_start <- "201611"
-yearmonth_end <- "202204"
+# Only data of full months can currently be read!
+evapo_p <- kwb.dwd::read_daily_data_over_shape(
+  file = shape_file,
+  variable = "evapo_p",
+  from = "201705",
+  to = "202312"
+)
+usethis::use_data(evapo_p, overwrite = TRUE)
+
+yearmonth_start <- "201705"
+yearmonth_end <- "202312"
 
 kwb.dwd:::list_daily_grids_germany_tgz("x")
 
@@ -64,7 +73,6 @@ dwd_daily_list <- stats::setNames(lapply(dwd_daily_vars, function(dwd_var) {
 
 
 dwd_daily <- dplyr::bind_rows(dwd_daily_list, .id = "parameter")
-"soil moisture" =
 
 
 kwb.dwd:::list_monthly_grids_germany_asc_gz("x")

data-raw/update_irrigation.R

@@ -0,0 +1,36 @@
+file_irrig <- "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/Bewässerungsmenge_Braunschweig_2020_2023.xlsx"
+
+irrig <- kwb.db::hsGetTable(file_irrig, tbl = "my_data") %>%
+  dplyr::select(! tidyselect::contains("mm"))
+
+names(irrig) <- stringr::str_remove(names(irrig), pattern = "\\s.*")
+
+
+irrig_update <- irrig %>%
+  dplyr::rename("Monat_num" = "Monat") %>%
+  dplyr::mutate("Monat" = withr::with_locale(new = c("LC_TIME" = "de_DE"), {
+    format(as.Date(sprintf("2000-%02d-01", Monat_num)), "%b")
+  })) %>%
+  dplyr::relocate(Monat) %>%
+  dplyr::relocate(Monat_num, .after = Monat) %>%
+  tidyr::gather(key = "Typ", value = "Menge_m3", - Monat, - Monat_num, - Jahr) %>%
+  dplyr::mutate(Typ = stringr::str_remove(Typ, "menge")) %>%
+  dplyr::arrange(Typ, Jahr, Monat_num)
+
+irrig_old <- readr::read_csv2("inst/extdata/input-data/Beregnungsmengen_AVB.csv")
+
+
+
+label_to_remove <- sprintf("%d-%02d_%s", irrig_old$Jahr, irrig_old$Monat_num, irrig_old$Typ) %>%
+  unique()
+
+irrig_update_tmp <- irrig_update %>%
+  dplyr::mutate(label = sprintf("%d-%02d_%s",Jahr,  Monat_num, Typ)) %>%
+  dplyr::filter(!label %in% label_to_remove) %>%
+  dplyr::select(-label)
+
+
+irrig_new <- dplyr::bind_rows(irrig_old, irrig_update_tmp) %>%
+  dplyr::arrange(Typ, Jahr, Monat_num)
+
+readr::write_csv2(irrig_new, "inst/extdata/input-data/Beregnungsmengen_AVB.csv")

inst/extdata/input-data/Beregnungsmengen_AVB.csv

@@ -47,6 +47,42 @@ Sep;9;2020;Grundwasser;0
 Okt;10;2020;Grundwasser;0
 Nov;11;2020;Grundwasser;0
 Dez;12;2020;Grundwasser;0
+Jan;1;2021;Grundwasser;0
+Feb;2;2021;Grundwasser;0
+Mrz;3;2021;Grundwasser;0
+Apr;4;2021;Grundwasser;0
+Mai;5;2021;Grundwasser;0
+Jun;6;2021;Grundwasser;145000
+Jul;7;2021;Grundwasser;38300
+Aug;8;2021;Grundwasser;153500
+Sep;9;2021;Grundwasser;2100
+Okt;10;2021;Grundwasser;0
+Nov;11;2021;Grundwasser;0
+Dez;12;2021;Grundwasser;0
+Jan;1;2022;Grundwasser;0
+Feb;2;2022;Grundwasser;0
+Mrz;3;2022;Grundwasser;0
+Apr;4;2022;Grundwasser;0
+Mai;5;2022;Grundwasser;80100
+Jun;6;2022;Grundwasser;469200
+Jul;7;2022;Grundwasser;387300
+Aug;8;2022;Grundwasser;278600
+Sep;9;2022;Grundwasser;41900
+Okt;10;2022;Grundwasser;0
+Nov;11;2022;Grundwasser;0
+Dez;12;2022;Grundwasser;0
+Jan;1;2023;Grundwasser;0
+Feb;2;2023;Grundwasser;0
+Mrz;3;2023;Grundwasser;0
+Apr;4;2023;Grundwasser;0
+Mai;5;2023;Grundwasser;9300
+Jun;6;2023;Grundwasser;462700
+Jul;7;2023;Grundwasser;243200
+Aug;8;2023;Grundwasser;0
+Sep;9;2023;Grundwasser;0
+Okt;10;2023;Grundwasser;0
+Nov;11;2023;Grundwasser;0
+Dez;12;2023;Grundwasser;0
 Jan;1;2017;Klarwasser;0
 Feb;2;2017;Klarwasser;748645
 Mrz;3;2017;Klarwasser;845499
@@ -95,3 +131,39 @@ Sep;9;2020;Klarwasser;1091727
 Okt;10;2020;Klarwasser;887857
 Nov;11;2020;Klarwasser;571506
 Dez;12;2020;Klarwasser;0
+Jan;1;2021;Klarwasser;0
+Feb;2;2021;Klarwasser;426326
+Mrz;3;2021;Klarwasser;582016
+Apr;4;2021;Klarwasser;756272
+Mai;5;2021;Klarwasser;863896
+Jun;6;2021;Klarwasser;1167427
+Jul;7;2021;Klarwasser;1267537
+Aug;8;2021;Klarwasser;1094828
+Sep;9;2021;Klarwasser;870819
+Okt;10;2021;Klarwasser;766444
+Nov;11;2021;Klarwasser;625077
+Dez;12;2021;Klarwasser;0
+Jan;1;2022;Klarwasser;0
+Feb;2;2022;Klarwasser;444913
+Mrz;3;2022;Klarwasser;610414
+Apr;4;2022;Klarwasser;939634
+Mai;5;2022;Klarwasser;1303181
+Jun;6;2022;Klarwasser;1303918
+Jul;7;2022;Klarwasser;1315510
+Aug;8;2022;Klarwasser;1220701
+Sep;9;2022;Klarwasser;1152751
+Okt;10;2022;Klarwasser;906681
+Nov;11;2022;Klarwasser;739971
+Dez;12;2022;Klarwasser;0
+Jan;1;2023;Klarwasser;0
+Feb;2;2023;Klarwasser;581089
+Mrz;3;2023;Klarwasser;618414
+Apr;4;2023;Klarwasser;754968
+Mai;5;2023;Klarwasser;1043533
+Jun;6;2023;Klarwasser;1237659
+Jul;7;2023;Klarwasser;1309563
+Aug;8;2023;Klarwasser;1120943
+Sep;9;2023;Klarwasser;1008015
+Okt;10;2023;Klarwasser;701111
+Nov;11;2023;Klarwasser;456938
+Dez;12;2023;Klarwasser;0