Downloading libraries for project:

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.5.2
library(ggplot2)
library(readr)
library(cowplot)
## Warning: package 'cowplot' was built under R version 4.5.2
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:patchwork':
## 
##     align_plots
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(rlang)
## 
## Attaching package: 'rlang'
## 
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
library(knitr)
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.2
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

Introduction - Need to do more Literature Review for the Narrative: Revisit when other sections are done: Methods - Re- do diagrams, explain flux theory, and theory needed for stats, locations, data cleaning, and more revisit when Results is done.

Results Section

Functions for creating figures:

# importing data sets
Flux_Data_With_Covars <- read_csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/Flux_Data_With_Covars.csv")
## Rows: 21090 Columns: 82
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (5): Pile, Chamber_Corrected, Turning_Status, Date, DATE_TIME
## dbl (77): Chamber_Elevation, DOY, air_temperature, RH, VPD, Tdew, wind_speed...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Functions to clean data

Plotting Moisture and Temperature for the experimental cycle.

plot_points_err <- function(df,
                            x, y_mean, y_se,
                            group,
                            cols = NULL,
                            shapes = NULL,
                            point_size = 2.0,
                            point_alpha = 0.85,
                            err_alpha = 0.6,
                            err_lwd = 0.4,
                            err_width = 0,
                            base_size = 12,
                            legend_position = "top",
                            drop_na = TRUE,
                            x_lab = NULL,
                            y_lab = NULL,
                            legend_title = NULL,
                            vline = NULL,
                            vlines_red = NULL,
                            vlines_black = NULL,
                            vline_lwd = 0.4,
                            vline_alpha = 0.8,
                            vline_linetype = "dashed",
                            vline_red_linetype   = "solid",
                            vline_black_linetype = "dashed",
                            vline_red_colour     = "#9b1c31",
                            vline_black_colour   = "#000000"){

  # capture column names
  x      <- enquo(x)
  y_mean <- enquo(y_mean)
  y_se   <- enquo(y_se)
  group  <- enquo(group)

  # optionally drop NA rows for required columns
  if (drop_na) {
    df <- df[complete.cases(
      df[, c(as_name(x), as_name(y_mean), as_name(y_se), as_name(group))]
    ), ]
  }

  p <- ggplot(df, aes(x = !!x, y = !!y_mean)) +

    # error bars first (so points draw on top)
    geom_errorbar(
      aes(ymin = (!!y_mean) - (!!y_se),
          ymax = (!!y_mean) + (!!y_se),
          colour = !!group),
      width = err_width,
      alpha = err_alpha,
      linewidth = err_lwd
    ) 
  
  if (!is.null(vlines_red)) {
  p <- p + geom_vline(
    xintercept = vlines_red,
    colour     = vline_red_colour,
    linetype   = vline_red_linetype,
    linewidth  = vline_lwd,
    alpha      = vline_alpha
  )
}

if (!is.null(vlines_black)) {
  p <- p + geom_vline(
    xintercept = vlines_black,
    colour     = vline_black_colour,
    linetype   = vline_black_linetype,
    linewidth  = vline_lwd,
    alpha      = vline_alpha
  )
  }
  p <- p +

    # points (no lines)
    geom_point(
      aes(colour = !!group, shape = !!group),
      size = point_size,
      alpha = point_alpha
    ) +

    theme_bw(base_size = base_size) +
    theme(
      legend.position = legend_position,
      panel.grid.minor = element_blank()
    )

  # manual scales if provided (keeps consistent styling across plots)
  if (!is.null(cols)) {
    p <- p + scale_colour_manual(values = cols)
  }
  if (!is.null(shapes)) {
    p <- p + scale_shape_manual(values = shapes)
  }
# ── LABELS SECTION ─────────────────────────────
  if (!is.null(x_lab) || !is.null(y_lab) || !is.null(legend_title)) {
    p <- p + labs(
      x = x_lab,
      y = y_lab,
      colour = legend_title,
      shape  = legend_title
    )
  }
  return(p)
}

Pile Temperature Plot using Daily Average Values Pile Moisture Using Daily Average Values

#turning data to add to plots if needed
turns_both_red <- c(159, 186, 208, 234, 262, 290)
turns_control_black <- c(173, 191, 222, 249, 276)
turns_both_red_week <- c(159, 186, 208, 234, 262, 290)/7
turns_control_black_week <- c(173, 191, 222, 249, 276)/7


pile_cols   <- c(C = "#000000", E = "#9b1c31")
pile_shapes <- c(C = 16, E = 17)

fig_moist_doy <- plot_points_err(
  daily_pile_doy,
  x = DOY,
  y_mean = pile_moist_mean,
  y_se = pile_moist_se,
  group = Pile,
  cols = pile_cols,
  shapes = pile_shapes,
  point_size = 1.8,
  x_lab = "Day of year",
  y_lab = expression("Volumetric water content (m"^3*" m"^-3*")"),
  legend_title = "Pile",

  vlines_red   = turns_both_red,
  vlines_black = turns_control_black,

  vline_red_linetype   = "solid",
  vline_black_linetype = "dashed"
)

fig_moist_doy

fig_moist_weekly <- plot_points_err(
  weekly_moist,
  x = week,
  y_mean = mean,
  y_se =  se,
  group = Pile,
  cols = pile_cols,
  shapes = pile_shapes,
  point_size = 1.8,
  x_lab = "Week of year",
  y_lab = expression("Volumetric water content (m"^3*" m"^-3*")"),
  legend_title = "Pile",
  
  vlines_red   = turns_both_red_week,
  vlines_black = turns_control_black_week,

  vline_red_linetype   = "solid",
  vline_black_linetype = "dashed"
)

fig_moist_weekly

fig_temp_doy <- plot_points_err(
  daily_pile_doy,
  x = DOY,
  y_mean =pile_temp_mean,
  y_se =  pile_temp_se,
  group = Pile,
  cols = pile_cols,
  shapes = pile_shapes,
  point_size = 1.8,
  x_lab = "Day of year",
  y_lab = expression("Pile Temperature Degrees Celsius"),
  legend_title = "Pile",
  
  vlines_red   = turns_both_red,
  vlines_black = turns_control_black,

  vline_red_linetype   = "solid",
  vline_black_linetype = "dashed"
)

fig_temp_doy

fig_temp_weekly <- plot_points_err(
  weekly_temp,
  x = week,
  y_mean =mean,
  y_se =  se,
  group = Pile,
  cols = pile_cols,
  shapes = pile_shapes,
  point_size = 1.8,
  x_lab = "Week of year",
  y_lab = expression("Pile Temperature Degrees Celsius"),
  legend_title = "Pile",
  
  vlines_red   = turns_both_red_week,
  vlines_black = turns_control_black_week,

  vline_red_linetype   = "solid",
  vline_black_linetype = "dashed"
)

fig_temp_weekly

Creating Tables with Pile Variables:

pile_vars_temp <- Cleaned_Data %>%
  filter(!is.na(pile_temp)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_temp = mean(pile_temp),
    sd_temp   = sd(pile_temp),
    se_temp   = sd_temp/sqrt(n_obs)
  )

pile_vars_water <- Cleaned_Data %>%
  filter(!is.na(pile_moist)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_moist = mean(pile_moist),
    sd_moist   = sd(pile_moist),
    se_moist   = sd_moist/sqrt(n_obs)
  )

Displaying Tables of Summary Stats for Moisture and Temperature

Summary statistics for volumetric water content by pile.
Pile n_obs mean_moist sd_moist se_moist
C 9177 0.490 0.327 0.003
E 9403 0.787 0.244 0.003
Summary statistics for Temperature by pile.
Pile n_obs mean_temp sd_temp se_temp
C 9192 54.069 8.992 0.094
E 9420 51.778 12.812 0.132

Examining Measured Fluxes by Pile

library(readxl)
co2_raw <-  read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/CO2_R_2_filtered.csv")

ch4_raw <-  read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/CH4_R_2_filtered.csv")

n2o_raw <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/N2O_R_2_filtered.csv")

Making Chart with Average values by Pile and SD,SE

pile_vars_co2 <- co2_raw %>%
  filter(!is.na(FCO2_DRY.LIN)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_co2 = mean(FCO2_DRY.LIN),
    sd_co2 = sd(FCO2_DRY.LIN),
    se_co2   = sd_co2/sqrt(n_obs)
  )

pile_vars_ch4 <- ch4_raw %>%
  filter(!is.na(FCH4_DRY.LIN)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_ch4 = mean(FCH4_DRY.LIN),
    sd_ch4   = sd(FCH4_DRY.LIN),
    se_ch4   = sd_ch4/sqrt(n_obs)
  )

pile_vars_n2o <- n2o_raw %>%
  filter(!is.na(FN2O_DRY.LIN)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_n2o = mean(FN2O_DRY.LIN),
    sd_n2o   = sd(FN2O_DRY.LIN),
    se_n2o   = sd_n2o/sqrt(n_obs)
  )

Displaying Average Values for Each Trace Gas:

kable(
  pile_vars_co2,
  digits = 3,
  caption = "Summary statistics for co2 flux by pile."
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("bordered", "striped")
  )
Summary statistics for co2 flux by pile.
Pile n_obs mean_co2 sd_co2 se_co2
C 5064 190.076 202.288 2.843
E 4215 191.757 193.987 2.988 
kable(
  pile_vars_ch4,
  digits = 3,
  caption = "Summary statistics for ch4 flux by pile."
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("bordered", "striped")
  )
Summary statistics for ch4 flux by pile.
Pile n_obs mean_ch4 sd_ch4 se_ch4
C 5355 7055.016 11741.80 160.456
E 4937 6594.771 10324.23 146.935 
kable(
  pile_vars_n2o,
  digits = 3,
  caption = "Summary statistics for n2o flux by pile."
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("bordered", "striped")
  )
Summary statistics for n2o flux by pile.
Pile n_obs mean_n2o sd_n2o se_n2o
C 6906 157.089 295.036 3.550
E 6911 98.481 245.689 2.955

Creating Daily and Weekly Plots for all Fluxes

# N2O
daily_n2o <- n2o_raw %>%
  mutate(DOY = floor(DOY)) %>%        # collapse fractional DOY
  group_by(Pile, DOY) %>%
  summarise(
    n_obs       = sum(!is.na(FN2O_DRY.LIN)),
    n2o_mean    = mean(FN2O_DRY.LIN, na.rm = TRUE),
    n2o_se      = se(FN2O_DRY.LIN),
    .groups     = "drop"
  ) %>%
  filter(n_obs >= 3)                  # tweak threshold if needed

# CO2
daily_co2 <- co2_raw %>%
  mutate(DOY = floor(DOY)) %>%
  group_by(Pile, DOY) %>%
  summarise(
    n_obs       = sum(!is.na(FCO2_DRY.LIN)),
    co2_mean    = mean(FCO2_DRY.LIN, na.rm = TRUE),
    co2_se      = se(FCO2_DRY.LIN),
    .groups     = "drop"
  ) %>%
  filter(n_obs >= 3)

# CH4
daily_ch4 <- ch4_raw %>%
  mutate(DOY = floor(DOY)) %>%
  group_by(Pile, DOY) %>%
  summarise(
    n_obs       = sum(!is.na(FCH4_DRY.LIN)),
    ch4_mean    = mean(FCH4_DRY.LIN, na.rm = TRUE),
    ch4_se      = se(FCH4_DRY.LIN),
    .groups     = "drop"
  ) %>%
  filter(n_obs >= 3)

weekly_n2o <- n2o_raw %>%
  mutate(
    week = floor((DOY - 1) / 7) + 1
  ) %>%
  group_by(Pile, week) %>%
  summarise(
    n_obs    = sum(!is.na(FN2O_DRY.LIN)),
    n2o_mean = mean(FN2O_DRY.LIN, na.rm = TRUE),
    n2o_se   = sd(FN2O_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
    DOY      = mean(DOY, na.rm = TRUE),   # x-position for plotting (optional but useful)
    .groups  = "drop"
  ) %>%
  filter(n_obs >= 3)

# CO2
weekly_co2 <- co2_raw %>%
  mutate(
    week = floor((DOY - 1) / 7) + 1
  ) %>%
  group_by(Pile, week) %>%
  summarise(
    n_obs    = sum(!is.na(FCO2_DRY.LIN)),
    co2_mean = mean(FCO2_DRY.LIN, na.rm = TRUE),
    co2_se   = sd(FCO2_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
    DOY      = mean(DOY, na.rm = TRUE),
    .groups  = "drop"
  ) %>%
  filter(n_obs >= 3)

# CH4
weekly_ch4 <- ch4_raw %>%
  mutate(
    week = floor((DOY - 1) / 7) + 1
  ) %>%
  group_by(Pile, week) %>%
  summarise(
    n_obs    = sum(!is.na(FCH4_DRY.LIN)),
    ch4_mean = mean(FCH4_DRY.LIN, na.rm = TRUE),
    ch4_se   = sd(FCH4_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
    DOY      = mean(DOY, na.rm = TRUE),
    .groups  = "drop"
  ) %>%
  filter(n_obs >= 3)

Plotting Values

fig_co2_daily <- plot_points_err(
  daily_co2,
  x = DOY, y_mean = co2_mean, y_se = co2_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red,
  vlines_black = turns_control_black,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_co2_daily

fig_co2_weekly <- plot_points_err(
  weekly_co2,
  x = week, y_mean = co2_mean, y_se = co2_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red_week,
  vlines_black = turns_control_black_week,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_co2_weekly

fig_ch4_daily <- plot_points_err(
  daily_ch4,
  x = DOY, y_mean = ch4_mean, y_se = ch4_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red,
  vlines_black = turns_control_black,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_ch4_daily

fig_ch4_weekly <- plot_points_err(
  weekly_ch4,
  x = week, y_mean = ch4_mean, y_se = ch4_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red_week,
  vlines_black = turns_control_black_week,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_ch4_weekly

fig_n2o_daily <- plot_points_err(
  daily_n2o,
  x = DOY, y_mean = n2o_mean, y_se = n2o_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red,
  vlines_black = turns_control_black,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_n2o_daily

fig_n2o_weekly <- plot_points_err(
  weekly_n2o,
  x = week, y_mean = n2o_mean, y_se = n2o_se,
  group = Pile,
  cols = pile_cols, shapes = pile_shapes,
  point_size = 1.4,
  vlines_red = turns_both_red_week,
  vlines_black = turns_control_black_week,
  vline_red_linetype = "solid",
  vline_black_linetype = "dashed"
)

fig_n2o_weekly

Looking at A Combination of all three gases using established 100 year conversion, these can easily be changed:

Co2 1 Ch4 28 N2o 265

To combine Trace Gases as cleaning was done on a per measurement basis:

Method 1: Looking at trace gases where all there trace gases where taken at once Method 2: Averaging Per Day and Comparing Daily Average GHG Potential Method 3: Weekly Averaging and Comparing Weekly GHG Potential

Method to keep only rows with all three trace gases

Flux_All_Gases_QC <- Flux_Data_With_Covars %>%
  filter(
    # CO2
    FCO2_DRY.LIN > 0,
    FCO2_DRY.LIN_R2 > 0.9,

    # CH4
    FCH4_DRY.LIN > 0,
    FCH4_DRY.LIN_R2 > 0.9,

    # N2O
    FN2O_DRY.LIN > 0,
    FN2O_DRY.LIN_R2 > 0.9
  )


Flux_Data_With_Covars %>%
  summarise(
    total_obs = n(),
    pass_all  = sum(
      FCO2_DRY.LIN > 0 & FCO2_DRY.LIN_R2 > 0.9 &
      FCH4_DRY.LIN > 0 & FCH4_DRY.LIN_R2 > 0.9 &
      FN2O_DRY.LIN > 0 & FN2O_DRY.LIN_R2 > 0.9,
      na.rm = TRUE
    ),
    percent_retained = 100 * pass_all / total_obs
  )
## # A tibble: 1 × 3
##   total_obs pass_all percent_retained
##       <int>    <int>            <dbl>
## 1     21090     5558             26.4
Flux_All_Gases_QC <- Flux_All_Gases_QC %>%
  mutate(
    # keep units explicit
    CO2_umol = FCO2_DRY.LIN,                  # µmol m-2 s-1
    N2O_umol = FN2O_DRY.LIN,                  # µmol m-2 s-1
    CH4_umol = FCH4_DRY.LIN * 0.001,          # nmol -> µmol

    # GWP100 CO2-equivalent (µmol CO2-eq m-2 s-1)
    GHG_CO2e_umol = CO2_umol + 28 * CH4_umol + 265 * N2O_umol
  )

Summary Stats of Co2e by pile when looking at measurements where all three gases passed qa/qc

pile_vars_co2e <- Flux_All_Gases_QC %>%
  filter(!is.na(GHG_CO2e_umol)) %>%
  group_by(Pile) %>%
  summarise(
    n_obs = n(),
    mean_CO2e = mean(GHG_CO2e_umol),
    sd_CO2e   = sd(GHG_CO2e_umol),
    se_CO2e = sd_CO2e/sqrt(n_obs)
  )


kable(
  pile_vars_co2e,
  digits = 3,
  caption = "Summary statistics for co2e flux by pile."
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("bordered", "striped")
  )
Summary statistics for co2e flux by pile.
Pile n_obs mean_CO2e sd_CO2e se_CO2e
C 3269 59570.66 101672.24 1778.259
E 2289 40046.13 96185.23 2010.414

We are going to write the code to plot this but we will plot all three together - as this is only the first method of plotting the values.