In bmcnellis/SDEF.analysis: Data Analysis for 2016 San Dimas Experimental Forest Project

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This is the 'final-product' processing vignette for the 2016 Adenostoma sapflux data.

Data product is saved in data/Adenostoma_2016_mphillips_process.rda.

Pre-processing is in vignettes/Adenostoma_2016_mphillips_preprocess.Rmd.

Used sapflux version 1.0.0.

Import and initial view:

library(SDEF.analysis)
library(sapflux)
data(Adenostoma_2016_mphillips_preprocessed)
flux_data <- Adenostoma_2016_mphillips_preprocessed
#sapflux::plot(flux_data, fastplot = TRUE)

Scale up:

flux_data <- sapflux::ScaleFlux(flux_data, to = "stem flux")
#print(flux_data@metadata[, c("stem.tag", "diameter")])
flux_data <- sapflux::ScaleFlux(flux_data, to = "plant flux")
#print(flux_data@metadata[, c("plant.tag", "diameter")])
#flux_data <- sapflux::ScaleFlux(flux_data, to = "site flux")
#sapflux::plot(flux_data, fastplot = TRUE)

The site flux should be compared to previous literature values in order to assess how well the analysis performed. Values may match more or less - TDP sapflux on a shrub hasn't been done much before and the methods are somewhat unproven, so this data should be presented as relative values:

time <- flux_data@time
time <- as.Date(time)
time <- as.character(time)
flux <- flux_data@data
suppressWarnings(max_daily <- aggregate(
    x = flux, by = list(time), FUN = function(x) { max(x, na.rm = TRUE) }
    ))
max_daily[, 1] <- as.POSIXct(strptime(max_daily[, 1], format = "%Y-%m-%d"))
max_daily <- as.data.frame(lapply(max_daily, function(x) {
  ifelse(is.infinite(x), NA, x)
}))
#max_daily[, 2] <- ifelse(is.infinite(max_daily[, 2]), NA, max_daily[, #2])
o0 <- lubridate::origin
max_daily[, 1] <- as.POSIXct(max_daily[, 1], origin = o0)
#colnames(max_daily) <- c("day", "max_flux")
colnames(max_daily)[-1] <- [e-mail]
save(max_daily, file = 'max_daily_adenostoma_unaggregated.Rdata')
plot(max_daily, xlab = "Time", ylab = "Maximum Daily Flux, Q (kg hr-1 m-2)")

Now smooth it out a bit. The NDVI measurements are from Landsat 8, so every 2 weeks.

max_daily_avg <- max_daily
max_daily_avg[, 2] <- zoo::rollapply(data = max_daily_avg[, 2], width = 14, fill = NA,
                                 FUN = function(x) {
                                   mean(x, na.rm = TRUE)
                                 })
plot(max_daily_avg, xlab = "Time", ylab = "Max Q, 2-wk running mean (kg hr-1 m-2)")

Now express as units relative to 'yearly peak activity'. Two objects are made here - one without the running average in the smoothing chunk, and one with it.

# No avg:
max_daily[, 2] <- max_daily[, 2] / (max(max_daily[, 2], na.rm = TRUE))
plot(max_daily, xlab = "Time", ylab = "Q relative to peak, unitless")
Adenostoma_2016_mphillips_processed_noavg <- max_daily
#use_data(Adenostoma_2016_mphillips_processed_noavg)
# With avg:
max_daily_avg[, 2] <- max_daily_avg[, 2] / (max(max_daily_avg[, 2], na.rm = TRUE))
plot(max_daily_avg, xlab = "Time", ylab = "2-wk avg. Q relative to peak, unitless")
Adenostoma_2016_mphillips_processed_withavg <- max_daily_avg
#use_data(Adenostoma_2016_mphillips_processed_withavg)

This final product is what goes into the final paper figure.

bmcnellis/SDEF.analysis documentation built on June 4, 2019, 10 a.m.