Framework for final L2 processing

pipeline/L1.qmd

@@ -54,6 +54,7 @@ options(warn = 2)
 
 source("helpers.R")
 source("L1-utils.R")
+source("metadata-utils.R")
 
 # Get the column metadata file
 column_md <- read.csv(file.path(params$METADATA_ROOT,
@@ -112,15 +113,7 @@ f <- function(dir_name, dirs_to_process, out_dir) {
         site <- dat$Site[1]
         plot <- dat$Plot[1]
 
-        # Check for metadata columns that are missing...
-        if(!all(column_md$Column %in% colnames(dat))) {
-            stop("Column metadata file ", params$METADATA_COLUMNS_TABLE,
-                 " has entries not in data: ", setdiff(column_md$Column, colnames(dat)))
-        }
-        # ...order and remove columns not in the metadata...
-        dat <- dat[column_md$Column]
-        # ...and sort rows
-        dat <- dat[order(dat$TIMESTAMP),]
+        dat <- check_drop_sort_columns(dat, column_md, params$METADATA_COLUMNS_TABLE)
 
         write_to_folders(dat, 
                          root_dir = out_dir, 
@@ -178,20 +171,17 @@ options(warn = oldwarn)
 L1 metadata template directory is `r params$L1_METADATA`.
 
 ```{r metadata}
-
-source("metadata-utils.R")
-
 # Write the overall README
 readme_fn <- file.path(params$METADATA_ROOT,
+                       params$L1_METADATA,
                        params$RELEASE_README_FILES, 
                        paste0("README_v", params$L1_VERSION, ".txt"))
-if(!file.exists(readme_fn)) stop("Couldn't find ", readme_fn)
-readme <- readLines(readme_fn)
-readme <- gsub("[VERSION]", params$L1_VERSION, readme, fixed = TRUE)
-readme <- gsub("[DATESTAMP]", params$L1_RELEASE_DATE, readme, fixed = TRUE)
-readme <- gsub("[OBSERVATIONS]", n_obs, readme, fixed = TRUE)
-readme <- gsub("[GIT_COMMIT]", GIT_COMMIT, readme, fixed = TRUE)
-readme <- gsub("[TIMEZONE]", params$L1_DATA_TIMEZONE, readme, fixed = TRUE)
+readme <- md_readme_substitutions(readme_fn,
+                                  params$L1_VERSION, 
+                                  params$L1_RELEASE_DATE, 
+                                  n_obs,
+                                  GIT_COMMIT,
+                                  params$L1_DATA_TIMEZONE)
 readme_outfn <- file.path(L1, basename(readme_fn))
 message("Writing overall README ", readme_outfn, "...")
 writeLines(readme, readme_outfn)
@@ -215,7 +205,7 @@ message("Main template has ", length(L1_metadata_template), " lines")
 message("Column metadata info has ", length(col_md_for_insert), " lines")
 message("Variable metadata info has ", length(var_md_for_insert), " lines")
 
-# Identify the main data directories in L1/{version}/, which are <site>_<year>
+# Identify the main data directories in L1, which are <site>_<year>
 data_dirs <- list.files(L1, pattern = "^[a-zA-Z]+_[0-9]{4}$")
 site_files_folder <- file.path(params$METADATA_ROOT, 
                               params$METADATA_SITE_FILES)
@@ -229,7 +219,7 @@ for(dd in data_dirs) {
     md <- gsub("[FOLDER_NAME]", dd, md, fixed = TRUE)
 
     # Insert info on data files into metadata
-    md <- md_insert_fileinfo(dd_full, md)
+    md <- md_insert_fileinfo(dd_full, md, filename_spacing = 55)
 
     # Insert column metadata
     col_info_pos <- grep("[COLUMN_INFO]", md, fixed = TRUE)

pipeline/L2-utils.R

@@ -0,0 +1,95 @@
+# L2-utils.R
+# Utility functions, and test code, used exclusively by L1.qmd
+# BBL June 2025
+
+library(testthat)
+
+# 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
+fill_gap <- function(x, gap_begin, gap_end, mac) {
+    if(!gap_end >= gap_begin) {
+        stop("I am seeing ", gap_begin, " ", gap_end, " in ", length(x))
+    }
+    stopifnot(gap_end >= gap_begin)
+    stopifnot(gap_begin > 0 && gap_begin <= length(x))
+    stopifnot(gap_end > 0 && gap_end <= length(x))
+    stopifnot(length(x) == length(mac))
+
+    # x and mac are numeric vectors of equal length
+    # gap_begin and gap_end are the indices of the first and last NAs of the gap
+
+    if(gap_begin > 1) {
+        left_adjust <- x[gap_begin - 1] - mac[gap_begin - 1]
+    } else {
+        left_adjust <- 0
+    }
+    #message("left_adjust ", left_adjust)
+    if(gap_end < length(x)) {
+        right_adjust <- x[gap_end + 1] - mac[gap_end + 1]
+    } else {
+        right_adjust <- 0
+    }
+    #message("right_adjust ", right_adjust)
+
+    # Calculate adjustments for the gap points plus the pre- and post-gap
+    # good data points that make up the calculation of left_adjust and right_adjust
+    adj <- seq(left_adjust, right_adjust, length.out = gap_end - gap_begin + 3)
+    # ...but then drop those out-of-gap adjustment values
+    adj <- adj[c(-1, -length(adj))]
+    # Return the adjusted mean annual cycle values to fill the gap
+    mac[gap_begin:gap_end] + adj
+}
+
+# Test code
+expect_identical(fill_gap(c(1, 2, NA), 3, 3, 1:3), 3)
+expect_identical(fill_gap(c(1, 2, NA), 3, 3, c(1, 2, NA)), NA_real_)
+expect_equivalent(fill_gap(c(1, NA, NA, 4), 2, 3, 1:4), 2:3)
+expect_equivalent(fill_gap(c(2, NA, NA, 5), 2, 3, 1:4), 3:4)
+
+expect_error(fill_gap(1:2, 3, 2), regexp = "I am seeing")
+expect_error(fill_gap(1:2, 1, 1, 1), regexp = "length")
+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) {
+    rle_x <- rle(is.na(x))
+    # Compute endpoints of run
+    end <- cumsum(rle_x$lengths)
+    start <- c(1, head(end, -1) + 1)
+    data.frame(start, end)[rle_x$values,] # return only TRUE value gaps
+}
+
+# Test code
+expect_equivalent(find_gaps(c(1, NA, 3)), data.frame(start = 2, end = 2))
+expect_equivalent(find_gaps(c(1, NA, NA, 3)), data.frame(start = 2, end = 3))
+expect_equivalent(find_gaps(c(1, NA)), data.frame(start = 2, end = 2))
+expect_equivalent(find_gaps(c(1, NA, 3, NA, 5)),
+                  data.frame(start = c(2, 4), end = c(2, 4)))
+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
+
+
+# Fill all the gaps in a vector based on the mean annual cycle
+fill_all_gaps <- function(x, mac) {
+    stopifnot(length(x) == length(mac))
+
+    gaps <- find_gaps(x)
+    gapfilled <- rep(NA_real_, length(x))
+
+    for(i in seq_len(nrow(gaps))) {
+        thisgap <- gaps$start[i]:gaps$end[i]
+        gapfilled[thisgap] <- fill_gap(x, gaps$start[i], gaps$end[i], mac)
+        # Just to make the graphs look good, return the boundary points too
+        # This code (these two ifs) should not be in production code
+        if(gaps$start[i] > 1)
+            gapfilled[gaps$start[i] - 1] <- x[gaps$start[i] - 1]
+        if(gaps$end[i] < length(x))
+            gapfilled[gaps$end[i] + 1] <- x[gaps$end[i] + 1]
+    }
+    return(gapfilled)
+}