Derived variables

.github/workflows/check-workflow.yaml

@@ -22,6 +22,7 @@ jobs:
           packages:
             any::quarto
             any::tidyr
+            any::dplyr
             any::readr
             any::lubridate
             any::ggplot2

pipeline/L0.qmd

@@ -34,7 +34,7 @@ This script
 
 -   Reshapes from wide to long; only `Logger`, `Table`, and `TIMESTAMP` don't get reshaped
 
--   Adds a unique observation ID (using `digest::digest()` )
+-   Adds a source file ID (using `digest::digest()` )
 
 -   Writes as CSV files with raw file hash info in filename
 
@@ -99,7 +99,7 @@ f <- function(fn, new_dir) {
         note <- "Used edited version"
     }
 
-    short_hash <- substr(digest::digest(file = fn, algo = "md5"), 1, 6)
+    short_hash <- substr(digest::digest(file = fn, algo = "md5"), 1, 8)
 
     # Check whether this file is on the removal list
     if(any(grepl(basefn, removal_list))) {
@@ -171,13 +171,7 @@ f <- function(fn, new_dir) {
             dat_long <- dat_long[!dupes,]
         }
 
-        # Add a unique ID column with the "md5" algorithm, using 16 of 32
-        # characters; "crc32" algorithm quickly generated duplicate IDs
-        # This is not *guaranteed* to be unique, but collisions are highly unlikely
-        dat_long$ID <- sapply(apply(dat_long, 1, paste, collapse = ""),
-                              FUN = function(x) {
-                                  substr(digest::digest(x, algo = "md5"), 1, 16)
-                              })
+        dat_long$Source_file <- short_hash
 
         # Write the new file
         message("\tWriting ", new_fqfn)

pipeline/L2-utils.R

@@ -4,6 +4,8 @@
 
 library(testthat)
 
+# Gap-filling functions -----------------------------------
+
 # Fill a gap in a vector x, based on the values immediately before
 # and after the gap and the mean annual cycle (mac)
 # Returns a numeric vector exactly as long as the gap
@@ -55,6 +57,8 @@ expect_error(fill_gap(1:2, -1, 2), regexp = "gap_begin")
 # Find the gaps (NAs) in a vector and return a data frame listing their
 # start and end positions
 find_gaps <- function(x) {
+    stopifnot(length(x) > 0)
+
     rle_x <- rle(is.na(x))
     # Compute endpoints of run
     end <- cumsum(rle_x$lengths)
@@ -72,7 +76,7 @@ expect_equivalent(find_gaps(NA), data.frame(start = 1, end = 1))
 
 expect_identical(nrow(find_gaps(c(1, 2))), 0L) # 0-row d.f. if no gaps
 expect_identical(nrow(find_gaps(c(1))), 0L) # 0-row d.f. if no gaps
-
+expect_error(find_gaps())
 
 # Fill all the gaps in a vector based on the mean annual cycle
 fill_all_gaps <- function(x, mac) {
@@ -91,5 +95,83 @@ fill_all_gaps <- function(x, mac) {
         if(gaps$end[i] < length(x))
             gapfilled[gaps$end[i] + 1] <- x[gaps$end[i] + 1]
     }
-    return(gapfilled)
+    return(round(gapfilled, 4))
+}
+
+# Derived variable functions -----------------------------------
+# Name of the function must start with "CALC_DERIVED" followed by
+# name of new research name (from derived variables table)
+
+# Generate full range of possible VWC values and corresponding roots
+# This code is outside a function because we only need to do it once
+VWCs <- seq(0, 1, length.out = 100)
+get_roots <- function(VWC) {
+    # Polynomial: -10.848-VWC + 1.302e-2*x - 5.105e-6*x^2 + 6.771e-10*x^3 = 0
+    coeffs <- c(-10.848 - VWC, 1.302e-2, -5.105e-6, 6.771e-10)
+    roots <- polyroot(coeffs)
+    # Return real part of real roots
+    rrt <- Re(roots[abs(Im(roots)) < 1e-8])
+    if(length(rrt)) {
+        return(rrt)
+    } else {
+        return(NA)
+    }
+}
+raws <- sapply(VWCs, get_roots)
+raws <- zoo::na.approx(raws) # fill in five missing values (when no real roots)
+
+compute_salinity <- function(T, VWC, EC) {
+    # the three data frames passed in should *normally* have the
+    # exact same dimensions and ordering, since they're linked
+    # TEROS data from the same site and plot
+    if(length(unique(sapply(list(T, VWC, EC), length))) != 1) {
+        stop("The T,VWC,EC data frames are not the same size!")
+    }
+
+    # The following steps are from Fausto Machado-Silva
+    # See also Hilhorst 2000, A pore water conductivity sensor,
+    # Soil Science Society of America Journal 64:6 1922–1925.
+
+    # 1. Convert TEROS EC to porewater EC
+    # METER recommends that σp not be calculated in soils with
+    # VWC < 0.10 m3/m3 using this method
+    VWC[VWC < 0.1] <- NA
+
+    # Convert conductivity to porewater conductivity
+    ep <- 80.3 - 0.37 * (T - 20)
+    esb0 <- 4.1
+    # Interpolate raw values for given VWCs
+    raw <- approx(VWCs, raws, xout = VWC)$y
+
+    eb <- (2.887e-9 * raw ** 3 - 2.080e-5 * raw ** 2 + 5.276e-2 * raw - 43.39) ** 2
+    cond_porewater <- (ep * EC) / (eb - esb0)
+
+    # 2. Convert porewater EC to salinity
+    -2.060e-1 + 5.781e-4 * cond_porewater
+}
+
+`CALC_DERIVED_soil-salinity-10cm` <- function(x) {
+
+    T <- x$`soil-temp-10cm`$Value
+    VWC <- x$`soil-vwc-10cm`$Value
+    EC <- x$`soil-EC-10cm`$Value
+
+    # Pick one of our input data frames, overwrite its Value column,
+    # and return it. The calling code in L2.qmd will take care of changing
+    # the research_name and Type columns
+    x$`soil-EC-10cm`$Value <- compute_salinity(T, VWC, EC)
+    return(x$`soil-EC-10cm`)
+}
+
+`CALC_DERIVED_soil-salinity-30cm` <- function(x) {
+
+    T <- x$`soil-temp-30cm`$Value
+    VWC <- x$`soil-vwc-30cm`$Value
+    EC <- x$`soil-EC-30cm`$Value
+
+    # Pick one of our input data frames, overwrite its Value column,
+    # and return it. The calling code in L2.qmd will take care of changing
+    # the research_name and Type columns
+    x$`soil-EC-30cm`$Value <- compute_salinity(T, VWC, EC)
+    return(x$`soil-EC-30cm`)
 }

pipeline/L2.qmd

@@ -15,6 +15,7 @@ params:
   METADATA_COLUMNS_TABLE: "L2_metadata_columns.csv"
   RELEASE_README_FILES: "readme_files/"
   L2_METADATA: "L2_metadata/"
+  DERIVED_VARS_TABLE: "derived_variables.csv"
   L2_RELEASE_DATE: "???"
   L2_DATA_TIMEZONE: "Etc/GMT+5"
   L2_VERSION: "???"
@@ -49,53 +50,120 @@ library(tidyr)
 library(readr)
 library(lubridate)
 library(ggplot2)
+theme_set(theme_bw())
 
 L2_QAQC <- file.path(params$DATA_ROOT, params$L2_QAQC)
 L2 <- file.path(params$DATA_ROOT, params$L2)
 
-# Find all L1_qaqc files and build a table of their
+# Find all L2_qaqc files and build a table of their
 # site/plot/research_name combinations
 files_to_process <- list.files(L2_QAQC, 
-                               pattern = "^[A-Z]{3}_[A-Z]{2,}_[0-9]{4}_.+csv$", 
+                               pattern = "^[A-Z]{3}_[A-Z]{1,}_[0-9]{4}_.+csv$", 
                                full.names = TRUE, recursive = TRUE)
 ftp_table <- tibble(file = files_to_process, filename = basename(files_to_process))
 ftp_table <- separate(ftp_table, filename, 
                       into = c("site", "plot", "year", "research_name", 
                                "level", "qaqc", "version"), sep = "_", 
                       remove = FALSE)
 ftp_table$level <- ftp_table$qaqc <- ftp_table$version <- NULL
-site_rn_table <- unique(ftp_table[c("site", "plot", "research_name")])
 
 QAQC_TABLE <- file.path(params$METADATA_ROOT, params$QAQC_TABLE)
 qaqct <- read_csv(QAQC_TABLE, col_types = "cccccc")
 
+# Get the variables metadata
+L2_VAR_MD <- read.csv(file.path(params$METADATA_ROOT,
+                                params$METADATA_VARS_TABLE))
+
 # Get the column metadata file
 column_md <- read.csv(file.path(params$METADATA_ROOT,
                                 params$L2_METADATA,
                                 params$METADATA_COLUMNS_TABLE),
                       comment.char = "#",
                       stringsAsFactors = FALSE)
 
+# Read the derived variables files
+derived_vars <- read_csv(file.path(params$METADATA_ROOT,
+                                   params$L2_METADATA,
+                                   params$DERIVED_VARS_TABLE))
+
 source("helpers.R")
 source("L2-utils.R")
 source("metadata-utils.R")
 ```
 
-I see `r nrow(site_rn_table)` site-plot-research_name combinations to process 
-in `r L2_QAQC`, representing `r length(files_to_process)` files.
+The derived variables table has `r nrow(derived_vars)` entries.
 
-QAQC table `r QAQC_TABLE` and has `r nrow(qaqct)` entries.
+For gap-filling, I see `r length(files_to_process)` files.
+
+QAQC table has `r nrow(qaqct)` entries.
 
 L2 column metadata has `r nrow(column_md)` entries.
 
 Output directory is `r L2`.
 
 HTML outfile is "`r params$html_outfile`".
 
-## MAC gap filling
+## Derived variables
 
-```{r processing}
+```{r derived-vars}
 errors <- 0
+
+# Process the derived variables table row by row
+for(i in seq_len(nrow(derived_vars))) {
+    dv <- derived_vars$research_name[i]
+    # Parse the 'needs_research_names' column and figure out which files are needed
+    needs <- trimws(strsplit(derived_vars$needs_research_names[i], ",")[[1]])
+    message("Calculating ", dv, " - need ", length(needs), 
+            " variables: ", derived_vars$needs_research_names[i])
+    files_needed <- ftp_table[ftp_table$research_name %in% needs,]
+
+    # Split by site and plot and year, and then for each split...
+    files_needed_split <- split(files_needed, ~ site + plot + year, drop = TRUE)
+    for(grp in files_needed_split) {
+        if(nrow(grp) == 0) stop("This shouldn't happen!") 
+        if(nrow(grp) < length(needs)) {
+            message("\tRequires ", length(needs), " variables but only ", nrow(grp), "are available")
+            next
+        }
+        si <- grp$site[1]
+        pl <- grp$plot[1]
+        yr <- grp$year[1]
+        message("\tLoading data for ", si, " ", pl, " ", yr)
+        # Read the files needed and put into a list, using research names as names
+        dat_needed <- lapply(grp$file, read_csv, col_types = "Tcccccccdii")
+
+        names(dat_needed) <- grp$research_name
+        # Call CALC_DERIVED_{dv} (the calculation functions are in L2-utils.R)
+        newdat <- do.call(paste0("CALC_DERIVED_", dv), list(dat_needed))
+        # Set columns to show that this is derived data
+        newdat$research_name <- dv
+        newdat$N_avg <- newdat$N_drop <- NA
+        fqfn <- write_to_folders(newdat, 
+                                 root_dir = L2, 
+                                 data_level = "L2",
+                                 site = si, 
+                                 plot = pl, 
+                                 derived_tempfile = TRUE)
+        # Add this file to the files-to-process table
+        # This is a one-by-one addition, not efficient, maybe TODO improve later
+        new_record <- tibble(file = fqfn, filename = basename(fqfn),
+                             site = si, plot = pl, year = yr,
+                             research_name = dv)
+        ftp_table <- rbind(ftp_table, new_record)
+    }
+    # add dv info to variables metadata table
+}
+
+```
+
+## MAC gap filling
+
+```{r mac-processing}
+
+# The site-research name table holds the unique combinations of those attributes
+# We will process all files at once that match the attributes for a given row
+site_rn_table <- unique(ftp_table[c("site", "plot", "research_name")])
+
 smry <- site_rn_table
 smry$Files <- 0
 smry$Dates <- ""
@@ -127,8 +195,13 @@ for(i in seq_len(nrow(site_rn_table))) {
     })
 
     # Compute the mean annual cycle (MAC) across all sensors and files (years)
-    message("\tComputing MAC...")
     all_data <- do.call(rbind, file_data)
+    if(all(is.na(all_data$Value))) {
+        message("\tNo non-NA data in this file! Moving on")
+        next
+    }
+
+    message("\tComputing MAC...")
     smry$N_vals[i] <- sum(!is.na(all_data$Value))
     yrs <- year(all_data$TIMESTAMP)
     smry$Dates[i] <- paste(format(min(all_data$TIMESTAMP), format = "%Y%m"),
@@ -151,25 +224,28 @@ for(i in seq_len(nrow(site_rn_table))) {
         # split doesn't work with NA values (???) so replace those...
         dat <- replace_na(dat, list(Instrument_ID = "", Sensor_ID = ""))
         # ...then split and apply
-        dat_split <- split(dat, ~ Instrument_ID + Sensor_ID)
+        dat_split <- split(dat, ~ Instrument_ID + Sensor_ID, drop = TRUE)
         dat_gf_list <- lapply(dat_split, function(x) {
+            if(nrow(x) == 0) {
+                stop("This isn't good; in ", paste(i, j, si, pl, rn))
+            }
             x <- x[order(x$TIMESTAMP),]
-            x$GF_MAC <- fill_all_gaps(x$Value, x$mac)
+            x$Value_GF_MAC <- fill_all_gaps(x$Value, x$mac)
             x
         })
         dat_gf <- do.call(rbind, dat_gf_list)
         rownames(dat_gf) <- NULL
         # Transfer good values into gap-fill column
-        nagf <- is.na(dat_gf$GF_MAC)
-        dat_gf$GF_MAC[nagf] <- dat_gf$Value[nagf]
+        nagf <- is.na(dat_gf$Value_GF_MAC)
+        dat_gf$Value_GF_MAC[nagf] <- dat_gf$Value[nagf]
         message("\tGap filled ", 
-                sum(is.na(dat_gf$Value)) - sum(is.na(dat_gf$GF_MAC)),
+                sum(is.na(dat_gf$Value)) - sum(is.na(dat_gf$Value_GF_MAC)),
                 " values in ", these_files$filename[j])
         # Update the summary and files-to-process tables
         f <- which(these_files$filename[j] == ftp_table$filename)
         ftp_table$N_vals[f] <- sum(!is.na(dat_gf$Value))
-        ftp_table$N_vals_gf[f] <- sum(!is.na(dat_gf$GF_MAC))
-        smry$N_vals_gf[i] <- smry$N_vals_gf[i] + sum(!is.na(dat_gf$GF_MAC))
+        ftp_table$N_vals_gf[f] <- sum(!is.na(dat_gf$Value_GF_MAC))
+        smry$N_vals_gf[i] <- smry$N_vals_gf[i] + sum(!is.na(dat_gf$Value_GF_MAC))
 
         # Clean up and check that the data frame has correct columns
         dat_gf <- check_drop_sort_columns(dat_gf, column_md, params$METADATA_COLUMNS_TABLE)
@@ -179,6 +255,7 @@ for(i in seq_len(nrow(site_rn_table))) {
                          data_level = "L2",
                          site = si,
                          plot = pl,
+                         variable_metadata = L2_VAR_MD,
                          version = params$L2_VERSION,
                          write_plots = params$write_plots) 
     }
@@ -239,7 +316,8 @@ col_md_for_insert <- paste(sprintf("%-15s", column_md$Column), column_md$Descrip
 # Get the variable metadata
 var_md_for_insert <- md_variable_info(file.path(params$METADATA_ROOT,
                                                 params$METADATA_VARS_TABLE),
-                                      only_these = qaqct$variable)
+                                      only_these = qaqct$variable,
+                                      derived_vars = derived_vars)
 
 message("Main template has ", length(L2_metadata_template), " lines")
 message("Column metadata info has ", length(col_md_for_insert), " lines")