Âdd read_selector()

with support for "time" and "solute" blocks, others still missing. In addition a write_selector() function would be helpful

in addition:
- remove "RootDepth" as required parameter
- improve read_obsnode() to also calculate "mass" (but offset compared to "input flux" (i.e. Prec*cTop)

NAMESPACE

@@ -13,6 +13,7 @@ export(read_meta_general)
 export(read_obsnode)
 export(read_profile)
 export(read_runinf)
+export(read_selector)
 export(read_solute)
 export(read_tlevel)
 export(run_model)

R/prepare_atmosphere_input.R

@@ -8,7 +8,7 @@
 get_atmosphere_headers <- function() {
 
   c("tAtm", "Prec", "rSoil", "rRoot", "hCritA", "rB", "hB", "ht",
-    "tTop", "tBot", "Ampl", "cTop", "cBot", "RootDepth")
+    "tTop", "tBot", "Ampl", "cTop", "cBot")
 }
 
 #' Prepare Atmosphere Input

R/read_obsnode.R

@@ -1,12 +1,17 @@
 #' Read Obs_Node.out
 #'
 #' @param path path to Obs_Node.out
+#' @param calculate_mass should masses for all concentrations be calculated, i.e.
+#' flux*conc(1-n) (default: TRUE)
 #' @param to_longer convert table to longer format (default: TRUE)
+#' @param debug print debug messages? (default: TRUE)
 #' @return tibble with Obs_Node time series data
 #' @importFrom stringr str_trim str_split str_replace_all str_remove
 #' @importFrom readr read_csv
 #' @export
-read_obsnode <- function(path, to_longer = TRUE) {
+read_obsnode <- function(path, to_longer = TRUE,
+                         calculate_mass = TRUE,
+                         debug = TRUE) {
   # Lese die Datei ein
   lines <- readLines(path)
 
@@ -22,7 +27,7 @@ read_obsnode <- function(path, to_longer = TRUE) {
     as.vector() %>%
     stringr::str_replace_all(pattern = "\\(\\s?", replacement = "") %>%
     stringr::str_remove("\\)") %>%
-    tolower
+    tolower()
 
   timeseries_idx <- grep("time", lines)
 
@@ -44,7 +49,7 @@ read_obsnode <- function(path, to_longer = TRUE) {
     is_conc <- grepl("conc", headers_sel)
 
     if(sum(is_conc) > 0) {
-      headers_sel[is_conc] <- sprintf("%s_%d",
+      headers_sel[is_conc] <- sprintf("%s%d",
                                       headers_sel[is_conc],
                                       seq_len(sum(is_conc)))
     }
@@ -62,15 +67,57 @@ read_obsnode <- function(path, to_longer = TRUE) {
 
   writeLines(dat_csv, path_csv)
 
-  dat <- readr::read_csv(path_csv)
+  dat <- readr::read_csv(path_csv,
+                         show_col_types = if(!debug) FALSE)
+
+
+  calc_mass <- function(dat) {
+
+    # Berechne die Massen dynamisch
+    conc_cols <- names(dat)[grepl("^conc", names(dat))]
+    mass_cols <- paste0("mass", seq_along(conc_cols))
+
+    dat_mass <- dat %>%
+      dplyr::mutate(dplyr::across(tidyselect::all_of(conc_cols), ~ flux * ., .names = "mass{col}"))
+
+    names(dat_mass) <- gsub("massconc", "mass", names(dat_mass))
+
+    dat_mass
+  }
 
-  if(to_longer) {
-    dat %>%
+  dat_long <- dat %>%
       tidyr::pivot_longer( - time) %>%
       tidyr::separate(col = "name", into = c("node_id", "variable"), sep = "_") %>%
-      dplyr::mutate(node = stringr::str_remove(node, "node") %>% as.integer())
-  } else {
-    dat
+      dplyr::mutate(node_id = stringr::str_remove(node_id, "node") %>% as.integer())
+
+
+    if(calculate_mass) {
+      n_conc <- dat_long %>%
+        dplyr::filter(grepl("^conc", variable)) %>%
+        dplyr::pull(variable) %>%
+        unique() %>%
+        length()
+
+      dat_long <- kwb.utils::catAndRun(sprintf("Calculating 'mass' for %d substance concentrations (flux*conc[%s]",
+                                               n_conc,
+                                               paste0(seq_len(n_conc), collapse = ",")),
+                           expr = {
+
+      dat_long %>%
+      tidyr::pivot_wider(names_from = "variable") %>%
+      calc_mass() %>%
+      tidyr::pivot_longer(names_to = "variable",
+                          - tidyselect::all_of(c("time","node_id")))
+
+    },
+    dbg = debug)
+    }
+
+
+  if(!to_longer) {
+    dat_long <- dat_long %>%
+      tidyr::pivot_wider(names_from = "variable")
   }
+  dat_long
 
 }

R/read_profile.R

@@ -7,7 +7,7 @@
 #' @importFrom stringr str_replace
 read_profile <- function(path) {
 
-  lines <- readLines(paths$profile)
+  lines <- readLines(path)
 
   header_idx <- grep("x", lines)
 
@@ -27,9 +27,9 @@ read_profile <- function(path) {
                            as.vector() %>% tolower())
 
   header_clean <- if(median(ncols) > length(header_names_file)) {
-    string_conc <- sprintf("conc_%d", seq_len(median(ncols) - length(header_names_file))+1)
+    string_conc <- sprintf("conc%d", seq_len(median(ncols) - length(header_names_file))+1)
 
-    c(stringr::str_replace(header_names_file, "conc", "conc_1"),
+    c(stringr::str_replace(header_names_file, "conc", "conc1"),
       string_conc)
 
   } else {

R/read_selector.R

@@ -0,0 +1,126 @@
+#' Read SELECTOR.in
+#'
+#' @param path path to SELECTOR.in
+#' @export
+
+read_selector <- function(path) {
+
+  lines <- readLines(path)
+
+  clean_line <- function(line, pattern = "\\s+") {
+    line %>%
+      stringr::str_trim() %>%
+      stringr::str_split(pattern, simplify = TRUE) %>%
+      as.vector()
+  }
+
+  header_values_to_list <- function (headers, values) {
+
+    lapply(values, function(value) {
+
+      is_num_val <- !is.na(suppressWarnings(as.numeric(value)))
+      if(is_num_val) {
+        as.numeric(value)
+      } else {
+        value
+      }
+    }) %>%
+      stats::setNames(headers)
+  }
+
+  blocks_idx_start <- grep("BLOCK", lines)
+  end_idx <- grep("END OF INPUT FILE 'SELECTOR.IN'", lines)
+
+  blocks_idx_end <- c(blocks_idx_start[seq_len(length(blocks_idx_start)-1)+1] - 1,
+                      end_idx - 1)
+
+
+  blocks_title_start <- lines[blocks_idx_start]
+  blocks_title_start_clean <- blocks_title_start %>%
+    stringr::str_remove_all("\\*|BLOCK") %>%
+    stringr::str_trim() %>%
+    stringr::str_replace_all(" ", "") %>%
+    stringr::str_replace_all(":", "_") %>%
+    stringr::str_remove_all("INFORMATION")
+
+  blocks <- tibble::tibble(name_clean = blocks_title_start_clean,
+                           name_org = blocks_title_start,
+                           start_idx = blocks_idx_start + 1,
+                           end_ix = blocks_idx_end)
+
+
+  "Pcp_File_Version=4"
+
+  block_time <- blocks[blocks$name_clean == "C_TIME",]
+  block_time_txt <- lines[block_time$start_idx:block_time$end_ix]
+
+  time <- c(
+    general = list(lapply(c(1,3,5), function(i) {
+    header_values_to_list(headers = clean_line(block_time_txt[i]),
+                          values = clean_line(block_time_txt[i + 1]))
+    })),
+    "TPrint" = list(lapply((grep("TPrint", block_time_txt)+1):length(block_time_txt),
+           function(i) {
+           clean_line(block_time_txt[i])
+           }
+    ) %>% unlist() %>% as.double()))
+
+
+  block_solute <- blocks[blocks$name_clean == "F_SOLUTETRANSPORT",]
+  block_solute_txt <- lines[block_solute$start_idx: block_solute $end_ix]
+
+  header_val_idx <- grep("Epsi|iNonEqul|kTopSolute|tPulse", block_solute_txt)
+  solute_transport_idx <- grep("Bulk.d.", block_solute_txt)
+  solute_reaction_idx <- grep("DifW", block_solute_txt)
+
+
+solute_transport <- list(transport =
+  lapply((solute_transport_idx+1):(min(solute_reaction_idx)-1), function(i) {
+  vec <- clean_line(block_solute_txt[i],  pattern = "\\s{2,}") %>% as.numeric()
+  names(vec) <- clean_line(block_solute_txt[solute_transport_idx],  pattern = "\\s{2,}")
+  vec
+}) %>% dplyr::bind_rows())
+
+
+solute_reaction <- list(reaction = stats::setNames(lapply(solute_reaction_idx, function(reac_idx) {
+  reac_max_idx <- if(reac_idx == max(solute_reaction_idx)) {
+    grep("kTopSolute", block_solute_txt) - 1
+  } else {
+    solute_reaction_idx[which(solute_reaction_idx == reac_idx)+1]-1
+  }
+
+   list(header_values_to_list(headers =  clean_line(block_solute_txt[reac_idx])[1:2],
+                                 values = clean_line(block_solute_txt[reac_idx+1])),
+           lapply((reac_idx+3):reac_max_idx, function(i) {
+             vec <- clean_line(block_solute_txt[i],
+                               pattern = "\\s{2,}") %>%
+               as.numeric()
+             names(vec) <- clean_line(block_solute_txt[reac_idx+2],
+                                      pattern = "\\s{2,}")
+             vec
+             }) %>%
+             dplyr::bind_rows()
+    )}),
+ nm = sprintf("solute_%d", seq_along(solute_reaction_idx))
+ )
+)
+
+
+
+
+solute <- c(general_1 = list(lapply(header_val_idx[1:2], function(i) {
+header_values_to_list(headers = clean_line(block_solute_txt[i]),
+                      values = clean_line(block_solute_txt[i + 1]))
+})),
+solute_transport,
+solute_reaction,
+general_2 = list(lapply(header_val_idx[3:4], function(i) {
+  header_values_to_list(headers = clean_line(block_solute_txt[i]),
+                        values = clean_line(block_solute_txt[i + 1]))
+})))
+
+list(time = time,
+     solute = solute)
+
+}
+

R/write_atmosphere.R

@@ -1,7 +1,7 @@
 #' Write "ATMOSPH.IN" input file
 #'
 #' @param atm tibble of input data as defined in \code{prepare_atmospherice_input}
-#' @param MaxAL Number of meteorological records
+#' @param MaxAL Number of meteorological records (default: nrow(atm))
 #' @param DailyVar TRUE if HYDRUS-1D is to generate daily variations in evaporation
 #' and transpiration (see section 2.7.2.)., otherwise: FALSE (default: FALSE)
 #' @param SinusVar TRUE if HYDRUS-1D is to generate sinusoidal variations in
@@ -31,7 +31,7 @@
 
 write_atmosphere <- function (
   atm,
-  MaxAL = 365,
+  MaxAL = nrow(atm),
   DailyVar = FALSE,
   SinusVar = FALSE,
   lLai = FALSE,