Merge pull request #53 from Plant-Functional-Trait-Course/next_version

Next version

R/flux_calc.R

@@ -7,6 +7,11 @@
 #' @param datetime_col column containing the datetime of each gas concentration
 #' measurements in slopes_df. The first one after cutting will be kept as
 #' datetime of each flux in the output.
+#' @param conc_unit unit in which the concentration of gas was measured
+#' ppm or ppb
+#' @param flux_unit unit in which the calculated flux will be
+#' `mmol` outputs fluxes in mmol*m^(-2)*h^(-1);
+#' `micromol` outputs fluxes in micromol*m^(-2)*h^(-1)
 #' @param cut_col column containing cutting information
 #' @param keep_arg name in cut_col of data to keep
 #' @param chamber_volume volume of the flux chamber in L,
@@ -16,7 +21,8 @@
 #' can also be a column in case it is a variable
 #' @param atm_pressure atmospheric pressure, assumed 1 atm,
 #' can be a constant (numerical) or a variable (column name)
-#' @param plot_area area of the plot in m^2, default for Three-D
+#' @param plot_area area of the plot in m^2, default for Three-D,
+#' can also be a column in case it is a variable
 #' @param cols_keep columns to keep from the input to the output.
 #' Those columns need to have unique values for each flux,
 #' as distinct() is applied.
@@ -30,7 +36,8 @@
 #' @param fit_type (optional) model used in flux_fitting, exponential,
 #' quadratic or linear.
 #' Will be automatically filled if slopes_df was produced using flux_quality().
-#' @return a dataframe containing fluxID, fluxes (in mmol*m^(-2)*h^(-1)),
+#' @return a dataframe containing fluxID, fluxes (in mmol*m^(-2)*h^(-1)
+#' or micromol*m^(-2)*h^(-1), depending on the value of `flux_unit`),
 #' temperature average for each flux,
 #' slope used for each flux calculation,
 #' the model used in `flux_fitting`,
@@ -40,14 +47,19 @@
 #' ungroup mutate case_when distinct left_join across everything
 #' @examples
 #' data(slopes0)
-#' flux_calc(slopes0, slope_col = "f_slope")
+#' flux_calc(slopes0,
+#' slope_col = "f_slope",
+#' conc_unit = "ppm",
+#' flux_unit = "mmol")
 #' @export
 
 
 
 flux_calc <- function(slopes_df,
                       slope_col,
                       datetime_col = "f_datetime",
+                      conc_unit,
+                      flux_unit,
                       cut_col = c(),
                       keep_arg = c(),
                       chamber_volume = 24.5,
@@ -60,19 +72,6 @@ flux_calc <- function(slopes_df,
                       temp_air_col = "temp_air",
                       temp_air_unit = "celsius",
                       fit_type = c()) {
-  fit_type <- flux_fit_type(
-    slopes_df,
-    fit_type = ((fit_type))
-  )
-
-  temp_air_unit <- match.arg(
-    ((temp_air_unit)),
-    c("celsius", "fahrenheit", "kelvin")
-  )
-
-  if (!is.double(((atm_pressure)))) stop("atm_pressure has to be a double")
-  if (!is.double(((plot_area)))) stop("plot_area has to be a double")
-
   colnames <- colnames(slopes_df)
   if (!(((slope_col)) %in% ((colnames)))) {
     stop("could not find slope_col in slopes_df")
@@ -83,15 +82,60 @@ flux_calc <- function(slopes_df,
   if (!(((temp_air_col)) %in% ((colnames)))) {
     stop("could not find temp_air_col in slopes_df")
   }
-
-
   if (length(setdiff(((cols_keep)), ((colnames)))) > 0) {
     stop("some names in cols_keep cannot be found in slopes_df")
   }
   if (length(setdiff(((cols_ave)), ((colnames)))) > 0) {
     stop("some names in cols_ave cannot be found in slopes_df")
   }
 
+
+  slopes_df_check <- slopes_df |>
+    select(
+      all_of(((slope_col))),
+      all_of(((temp_air_col))),
+      all_of(((datetime_col)))
+    )
+
+  df_ok <- flux_fun_check(slopes_df_check,
+                          fn = list(
+                            is.numeric,
+                            is.numeric,
+                            is.POSIXct
+                          ),
+                          msg = rep(c(
+                            "has to be numeric",
+                            "has to be POSIXct"
+                          ),
+                          c(2, 1)
+                          ),
+                          origdf = slopes_df)
+
+
+  if (any(!df_ok))
+    stop("Please correct the arguments", call. = FALSE)
+
+
+  fit_type <- flux_fit_type(
+    slopes_df,
+    fit_type = ((fit_type))
+  )
+
+  temp_air_unit <- match.arg(
+    ((temp_air_unit)),
+    c("celsius", "fahrenheit", "kelvin")
+  )
+
+  conc_unit <- match.arg(
+    ((conc_unit)),
+    c("ppm", "ppb")
+  )
+
+  flux_unit <- match.arg(
+    ((flux_unit)),
+    c("micromol", "mmol")
+  )
+
   if (is.double((chamber_volume))) {
     slopes_df <- slopes_df |>
       mutate(
@@ -134,6 +178,19 @@ flux_calc <- function(slopes_df,
       )
   }
 
+  if (is.double((plot_area))) {
+    slopes_df <- slopes_df |>
+      mutate(
+        plot_area = ((plot_area))
+      )
+  }
+
+  if (is.character(((plot_area)))) {
+    slopes_df <- slopes_df |>
+      rename(
+        plot_area = all_of(((plot_area)))
+      )
+  }
 
 
   slopes_df <- slopes_df |>
@@ -159,11 +216,11 @@ flux_calc <- function(slopes_df,
   slope_temp <- slopes_df |>
     select(
       "f_slope_calc", "f_fluxID", "air_temp", "chamber_volume",
-      "tube_volume", "atm_pressure", "f_datetime"
+      "tube_volume", "atm_pressure", "f_datetime", "plot_area"
     ) |>
     group_by(
       .data$f_fluxID, .data$f_slope_calc, .data$chamber_volume,
-      .data$tube_volume, .data$atm_pressure
+      .data$tube_volume, .data$atm_pressure, .data$plot_area
     ) |>
     summarise(
       temp_air_ave = mean(.data$air_temp, na.rm = TRUE),
@@ -210,6 +267,20 @@ flux_calc <- function(slopes_df,
   message("Calculating fluxes...")
 
   r_const <- 0.082057
+  message("R constant set to 0.082057")
+
+
+  # putting slope in ppm/s
+  if (((conc_unit)) == "ppm") {
+    message("Concentration was measured in ppm")
+  }
+  if (((conc_unit)) == "ppb") {
+    message("Concentration was measured in ppb")
+    slope_ave <- slope_ave |>
+      mutate(
+        f_slope_calc = .data$f_slope_calc * 0.001 # now the slope is in ppm/s
+      )
+  }
 
   fluxes <- slope_ave |>
     mutate(
@@ -218,9 +289,8 @@ flux_calc <- function(slopes_df,
         (.data$f_slope_calc * .data$atm_pressure * .data$volume_setup)
         / (((r_const)) *
            .data$temp_air_ave
-           * ((plot_area))) # flux in micromol/s/m^2
-        * 3600 # secs to hours
-        / 1000, # micromol to mmol flux is now in mmol/m^2/h
+           * .data$plot_area) # flux in micromol/s/m^2
+        * 3600, # secs to hours, flux is now in micromol/m^2/h
       temp_air_ave = case_when(
         ((temp_air_unit)) == "celsius" ~ .data$temp_air_ave - 273.15,
         ((temp_air_unit)) == "fahrenheit"
@@ -229,5 +299,18 @@ flux_calc <- function(slopes_df,
       ),
       model = ((fit_type))
     )
+
+  # output unit
+  if (((flux_unit)) == "micromol") {
+    message("Fluxes are in micromol/m2/h")
+  }
+  if (((flux_unit)) == "mmol") {
+    fluxes <- fluxes |>
+      mutate(
+        flux = .data$flux / 1000
+      )
+    message("Fluxes are in mmol/m2/h")
+  }
+
   fluxes
 }

R/flux_check_item.R

@@ -0,0 +1,25 @@
+#' check the items inside `flux_fun_check`
+#' @param arg argument to be checked by `fn`
+#' @param fn function to check `arg`
+#' @param msg message to display in case `arg` is the wrong class
+#' @param narg name of `arg`
+#' @param df_name name of `arg` in case it is a data frame
+#' @author Adam Klimes
+
+
+flux_check_item <- function(arg,
+                            fn,
+                            msg,
+                            narg,
+                            df_name = NA) {
+  isok <- fn(arg)
+  if (!isok) {
+    msg_parts <- if (!is.na(df_name))
+      c("Column ", paste0(" of data frame ", df_name, " "))
+    else c("Argument ", " ")
+    message(paste0(msg_parts[1], narg, msg_parts[2], msg))
+    FALSE
+  } else {
+    TRUE
+  }
+}

R/flux_fitting.R

@@ -51,6 +51,40 @@ flux_fitting <- function(conc_df,
                          roll_width = 15,
                          t_zero = 0,
                          fit_type) {
+  args_ok <- flux_fun_check(list(
+    start_cut = ((start_cut)),
+    end_cut = ((end_cut))
+  ),
+  fn = list(is.numeric, is.numeric),
+  msg = rep("has to be numeric", 2))
+
+  conc_df_check <- conc_df |>
+    select(
+      all_of(((conc_col))),
+      all_of(((start_col))),
+      all_of(((end_col))),
+      all_of(((datetime_col)))
+    )
+
+  conc_df_ok <- flux_fun_check(conc_df_check,
+                               fn = list(
+                                 is.numeric,
+                                 is.POSIXct,
+                                 is.POSIXct,
+                                 is.POSIXct
+                               ),
+                               msg = rep(c(
+                                 "has to be numeric",
+                                 "has to be POSIXct"
+                               ),
+                               c(1, 3)
+                               ),
+                               origdf = conc_df)
+
+
+  if (any(!c(args_ok, conc_df_ok)))
+    stop("Please correct the arguments", call. = FALSE)
+
   fit_type <- flux_fit_type(
     ((conc_df)),
     fit_type = ((fit_type))

R/flux_fitting_exp.R

@@ -48,6 +48,31 @@ flux_fitting_exp <- function(conc_df,
                              datetime_col = "f_datetime",
                              conc_col = "f_conc",
                              fluxid_col = "f_fluxID") {
+
+  args_ok <- flux_fun_check(list(
+    t_window = ((t_window)),
+    cz_window = ((cz_window)),
+    b_window = ((b_window)),
+    a_window = ((a_window)),
+    roll_width = ((roll_width)),
+    start_cut = ((start_cut)),
+    end_cut = ((end_cut))
+  ),
+  fn = list(
+    is.numeric,
+    is.numeric,
+    is.numeric,
+    is.numeric,
+    is.numeric,
+    is.numeric,
+    is.numeric
+  ),
+  msg = rep("has to be numeric", 7))
+
+  if (any(!args_ok))
+    stop("Please correct the arguments", call. = FALSE)
+
+
   conc_df <- conc_df |>
     rename(
       f_start = all_of((start_col)),
@@ -57,19 +82,6 @@ flux_fitting_exp <- function(conc_df,
       f_fluxID = all_of((fluxid_col))
     )
 
-
-
-  if (!is.double(t_window)) stop("t_window has to be a double")
-  if (!is.double(cz_window)) stop("cz_window has to be a double")
-  if (!is.double(b_window)) stop("b_window has to be a double")
-  if (!is.double(a_window)) stop("a_window has to be a double")
-  if (!is.double(roll_width)) stop("roll_width has to be a double")
-  if (!is.double(start_cut)) stop("start_cut has to be a double")
-  if (!is.double(end_cut)) stop("end_cut has to be a double")
-
-
-
-
   length_flux_max <- conc_df |>
     mutate(
       length_flux = difftime(.data$f_end, .data$f_start, units = "sec"),
@@ -222,7 +234,11 @@ flux_fitting_exp <- function(conc_df,
     left_join(tz_df, by = "f_fluxID") |>
     group_by(.data$f_fluxID) |>
     mutate(
-      f_Cb = .data$f_conc[.data$f_time_cut == .data$tz_est - ((b_window))]
+      diff = .data$f_time_cut - .data$tz_est + ((b_window))
+    ) |>
+    distinct(.data$diff, .keep_all = TRUE) |>
+    mutate(
+      f_Cb = .data$f_conc[which.min(abs(.data$diff))]
     ) |>
     ungroup() |>
     select("f_fluxID", "f_Cb") |>
@@ -232,7 +248,11 @@ flux_fitting_exp <- function(conc_df,
     group_by(.data$f_fluxID) |>
     mutate(
       ta = .data$length_window - ((a_window)),
-      Ca = .data$f_conc[.data$f_time_cut == .data$ta]
+      ta_diff = .data$f_time_cut - .data$ta
+    ) |>
+    distinct(.data$ta_diff, .keep_all = TRUE) |>
+    mutate(
+      Ca = .data$f_conc[which.min(abs(.data$ta_diff))]
     ) |>
     ungroup() |>
     select("f_fluxID", "ta", "Ca") |>

R/flux_fitting_lin.R

@@ -27,6 +27,16 @@ flux_fitting_lin <- function(conc_df,
                              datetime_col = "f_datetime",
                              conc_col = "f_conc",
                              fluxid_col = "f_fluxID") {
+  args_ok <- flux_fun_check(list(
+    start_cut = ((start_cut)),
+    end_cut = ((end_cut))
+  ),
+  fn = list(is.numeric, is.numeric),
+  msg = rep("has to be numeric", 2))
+
+  if (any(!args_ok))
+    stop("Please correct the arguments", call. = FALSE)
+
   conc_df <- conc_df |>
     rename(
       f_start = all_of(((start_col))),
@@ -36,8 +46,6 @@ flux_fitting_lin <- function(conc_df,
       f_fluxID = all_of(((fluxid_col)))
     )
 
-  if (!is.double(((start_cut)))) stop("start_cut has to be a double")
-  if (!is.double(((end_cut)))) stop("end_cut has to be a double")
 
   length_flux_max <- conc_df |>
     mutate(