Saturo Calculations & Data Validation

Notes to self: hide code? Units for pressure, flux, etc Add units to Kfs_table

1. Introduction

Output excel files from Saturo runs contain a calculated field-saturated hydraulic conductivity (Kfs) and it’s associated error (Kfs error), however, these values are based on the last cycle of any run. Paying attention to the live data outputs on the Saturo control unit in the field may be most important during the beginning of the run to avoid wasting time and water. However, in some cases, the last run can still contain noisy or irregular data due to unavoidable issues. Kfs and Kfs error from a previous cycle may be able to be substituted for the last cycle.

This document provides an annotated explanation of the R script used to manually calculate Kfs and Kfs error for every cycle during a Saturo run. The code below is broken up into sections to make clear what each portion of the code is doing.

# Yellow highlighted box means you need to manually change something within that section of the code yourself.

2. Load libraries

Ensure that you have these packages installed prior to running the code.

library(tidyverse)
library(readxl)
library(shiny)

3. Bring in your data

Change the path directory to where the file is stored (file_path)
Change the name of the file to the one of interest (file_name)

file_path = "C:/Users/Amanda.Pennino/OneDrive - USDA/Documents/methods/Saturo/data/"
file_name = "Run1.xlsx"

# Change directory to the file you want to analyze
dat <- read_excel(paste0(file_path, file_name), sheet = "Raw Data")


# Create summary of run settings
# Change directory to the file you want to analyze
summary <- read_excel(paste0(file_path, file_name), 
                      sheet = "Summary",
                      col_names = c("Settings", "Value"),
                      skip = 2,
                      col_types = c("guess", "numeric"))

4. Extract variables for calculations

# Reformat summary to wide format
summary <- summary %>% 
  drop_na(Settings) %>%
  spread(key = Settings, value = Value)

# Extract out units of measurement
## This is set on the control unit prior to the run
units <- data.frame(names = names(summary)) %>% 
  mutate(units = str_extract(.[[1]], "(?<=\\()([^()]*?)(?=\\)[^()]*$)")) %>%
  mutate(names = str_replace(.[[1]], " \\s*\\([^\\)]+\\)", "")) 

units2 <- data.frame(names = names(dat)) %>%
  mutate(units = str_extract(.[[1]], "(?<=\\()([^()]*?)(?=\\)[^()]*$)")) %>%
  mutate(names = str_replace(.[[1]], " \\s*\\([^\\)]+\\)", "")) 

units <- rbind(units, units2)

#Change column names in the data sheet to take out units of measurement
names(dat) <- units2[[1]]
rm(units2) #remove data frame


# Extract variables for calculations
presoak <- summary[[1,16]]
hold <- summary[[1,5]]
depth <- summary[[1,7]] # Ring depth
cycles <- summary[[1,10]] 
radius <- 7.5  # Ring radius
delta <- round((0.993 * depth) + (0.578 * radius), digits = 1)  # Calculate delta (ring geometry/shape factor)

5. Determine high/low pressure periods

# Determine high/low pressure periods
dat <- dat %>%
  mutate(setting = ifelse(Time <= presoak, "presoak", 
                          ifelse(Pressure > 6, "high", 
                                 ifelse(Pressure < 6, "low", NA) ))) 

nhigh <- as.numeric(nrow(dat[dat$setting == "high",])) 
nlow <- as.numeric(nrow(dat[dat$setting == "low",]))


# Add cycle numbers
dat <- dat %>%
  group_by(setting) %>%
  mutate(cycle = ifelse(setting == "presoak", 0,
                        ifelse(setting == "high", rep(c(1:cycles), each = nhigh/cycles, length.out = n()), 
                               ifelse(setting == "low", rep(c(1:cycles), each = nlow/cycles, length.out = n()), NA)))) %>%
  ungroup() %>%
  group_by(cycle, setting) %>%
  slice(3:n()) #takes out the first 2 minutes of each pressure setting for every cycle

5. Calculate the Kfs & Kfs error for each cycle

# Kfs
Kfs <- dat %>%
  filter(setting != "presoak") %>%
  group_by(setting, cycle) %>%
  summarise(mFlux = mean(Flux),
            mPress = mean(Pressure)) %>%
  pivot_wider(names_from = setting, values_from = c(mFlux, mPress)) %>%
  mutate(Kfs = (delta * (mFlux_high - mFlux_low)) / (mPress_high - mPress_low)) %>%
  select(cycle, Kfs)


Kfs <- dat %>%
  filter(setting != "presoak") %>%
  group_by(setting, cycle) %>%
  summarise(mFlux = mean(Flux),
            mPress = mean(Pressure)) %>%
  pivot_wider(names_from = setting, values_from = c(mFlux, mPress)) %>%
  mutate(Kfs = (delta * (mFlux_high - mFlux_low)) / (mPress_high - mPress_low)) %>%
  select(cycle, Kfs)


# Kfs error
high <- dat %>%
  filter(setting == "high") %>%
  group_by(setting) %>%
  rename(highP = Pressure,
         highQ = Flux) %>%
  ungroup()%>%
  select(highP, highQ, cycle)
  
low <- dat %>%
  filter(setting == "low") %>%
  group_by(setting) %>%
  rename(lowP = Pressure,
         lowQ = Flux) %>%
  ungroup()%>%
  select(lowP, lowQ)


dat2 <- cbind(high, low)
dat2 <- left_join(dat2, Kfs, by = "cycle")

Kfs_error <- dat2 %>%
  mutate(row_Kfs = (delta * (highQ - lowQ)) / (highP-lowP),
         row_error = (row_Kfs - Kfs)^2 ) %>%
  group_by(cycle) %>%
  summarise(Kfs_error = sqrt(sum(row_error)) / (hold - 2) )


Kfs_data <- left_join(Kfs, Kfs_error, by = "cycle")

6. View the calculated results by cycle

# Change view of data to scientific notation with 3 significant digits
Kfs_data$Kfs <- formatC(Kfs_data$Kfs, format = "e", 2)
Kfs_data$Kfs_error <- formatC(Kfs_data$Kfs_error, format = "e", 2)


names(Kfs_data) <- c("Cycle", 
                     paste0("Kfs (", units$units[units$names == "Kfs"], ")"), 
                     paste0("Kfs Error (", units$units[units$names == "Kfs Error"], ")"))

# Show output table
as.matrix(Kfs_data)

##      Cycle Kfs (cm/s) Kfs Error (cm/s)
## [1,] "1"   "1.85e-03" "7.34e-05"      
## [2,] "2"   "1.49e-03" "3.93e-05"      
## [3,] "3"   "1.45e-03" "5.73e-05"

# Shiny app for data viewing
## Define UI
ui <- fluidPage(
  titlePanel("Explore the data"),
  sidebarLayout(
    sidebarPanel(
      selectInput(
        inputId = "y_var",
        label = "Choose variable:",
        choices = c("Water Level", "Pressure", "Flux", "Volume"),
        selected = "Water Level"
      )
    ),
    mainPanel(
      plotOutput("line_plot", height = "400px"), #ggplot2 line plot
      hr(), #horizontal rule for separation
      tableOutput("static_table") 
    )
  )
)

# Define server logic
server <- function(input, output) {
  # Render the ggplot2 line plot
  output$line_plot <- renderPlot({
    ggplot(data = dat, aes(x = Time, y = .data[[input$y_var]])) +
      geom_line(color = "blue", size = 1) +
      labs(
        title = paste(input$y_var, "vs Time"),
        x = "Time",
        y = input$y_var
      ) +
      theme_minimal(base_size = 15)
  })
  
  # Render the static table
  output$static_table <- renderTable({
    Kfs_data
  })
}

# Run the application
shinyApp(ui = ui, server = server)

Shiny applications not supported in static R Markdown documents

7. Read out your results

This will be saved in the same location as where the original datasheet was.

#Export data as a .csv file
write.csv(Kfs_data, paste0(file_path, "Results_", file_name))