Are thermal cameras sensitive enough to reflect the sap and dendro patterns

TODO MEETING NOV 13:

• filter out night sap flow for total

• SWC to ∆T (midday)

• ET (flux data) to ∆T (midday)

• Compare 30min sap vs ∆T 30min

• for dendro:plot diurnal and asses MDS after that

• add time series of sap vs ∆T WITH rain fall pulses on the same panel to see if there are responses

• for daily patterns, add points to the lines to make it clearer

• Work In Progress: add good timeseries.

• DONE go to Dendro script and double check the units (microns) and put that on

• DONE add precip

• DONE add raw, growth values to trace how TWD is calculated

  • Tree Physiology Figure

Globals

LAB_JS_30MIN <- expression(J[s] ~ (g ~ m^-2 ~ s^-1))

LAB_JS_DAILY <- expression(J[s] ~ (g ~ m^-2 ~ day^-1))

LAB_JS_30MIN_TXT <- "Js (g m^-2 s^-1)"

LAB_JS_DAILY_TXT <- "Js (g m^-2 day^-1)"

MIDDAY_START <- 10

MIDDAY_END <- 14

MPJ_LAT <- 34.43

MPJ_LON <- -106.13

MPJ_SOLAR_TZ <- "America/Phoenix"

SPECIES_COLS <- c("piñon" = "#40B0A6", "juniper" = "#E1BE6A")

# CHANGED: centralize repeated parsing logic here so later chunks use the same
# tree/species naming rules instead of re-implementing them.
species_from_tree <- function(tree_id) {
  case_when(
    str_starts(tree_id, "p") ~ "piñon",
    str_starts(tree_id, "j") ~ "juniper",
    TRUE ~ NA_character_
  )
}

tree_from_sensor <- function(sensor_id) {
  str_remove(sensor_id, "[ab]_sap$")
}

# CHANGED: pull the twd_norm.Rmd predawn/midday normalization logic into one
# reusable helper so daily MDS plus both predawn- and TWD4-based normalization
# stay in sync.
calc_normalized_dendro_metrics <- function(
  data,
  lat = MPJ_LAT,
  lon = MPJ_LON,
  solar_tz = MPJ_SOLAR_TZ
) {
  required_cols <- c("tree", "species", "ts", "value", "twd")

  if (!all(required_cols %in% names(data))) {
    stop(
      "calc_normalized_dendro_metrics() requires columns: ",
      paste(required_cols, collapse = ", ")
    )
  }

  date_range <- range(as_date(data$ts), na.rm = TRUE)
  ts_tz <- attr(data$ts, "tzone")
  na_tz <- if (length(ts_tz) == 0 || is.na(ts_tz) || ts_tz == "") {
    "UTC"
  } else {
    ts_tz
  }

  sun_times <- suncalc::getSunlightTimes(
    date = seq.Date(date_range[1], date_range[2], by = "day"),
    lat = lat,
    lon = lon,
    keep = c("sunrise", "solarNoon", "sunset"),
    tz = solar_tz
  ) |>
    mutate(date = as_date(date)) |>
    select(date, sunrise, solarNoon, sunset)

  dendro_with_sun <- data |>
    mutate(date = as_date(ts)) |>
    left_join(sun_times, by = "date") |>
    mutate(
      is_predawn = ts < sunrise,
      is_midday = ts >= solarNoon & ts <= sunset
    )

  twd_daily <- dendro_with_sun |>
    group_by(tree, species, date) |>
    summarise(
      twd_predawn = if (any(is_predawn & !is.na(twd))) min(twd[is_predawn], na.rm = TRUE) else NA_real_,
      twd_predawn_ts = if (any(is_predawn & !is.na(twd))) {
        ts[is_predawn & !is.na(twd)][which.min(twd[is_predawn & !is.na(twd)])][1]
      } else {
        as.POSIXct(NA, tz = na_tz)
      },
      twd_4 = if (all(is.na(twd))) NA_real_ else max(twd, na.rm = TRUE),
      .groups = "drop"
    )

  mds_daily <- dendro_with_sun |>
    group_by(tree, species, date) |>
    summarise(
      predawn_max = if (any(is_predawn & !is.na(value))) max(value[is_predawn], na.rm = TRUE) else NA_real_,
      predawn_max_ts = if (any(is_predawn & !is.na(value))) {
        ts[is_predawn & !is.na(value)][which.max(value[is_predawn & !is.na(value)])][1]
      } else {
        as.POSIXct(NA, tz = na_tz)
      },
      midday_min = if (any(is_midday & !is.na(value))) min(value[is_midday], na.rm = TRUE) else NA_real_,
      midday_min_ts = if (any(is_midday & !is.na(value))) {
        ts[is_midday & !is.na(value)][which.min(value[is_midday & !is.na(value)])][1]
      } else {
        as.POSIXct(NA, tz = na_tz)
      },
      mds = predawn_max - midday_min,
      .groups = "drop"
    )

  twd_daily |>
    left_join(mds_daily, by = c("tree", "species", "date")) |>
    group_by(tree) |>
    mutate(
      mds_max = if (all(is.na(mds))) {
        NA_real_
      } else {
        quantile(mds, probs = 0.99, na.rm = TRUE, names = FALSE, type = 7)
      },
      twd_norm_pd = twd_predawn / mds_max,
      twd_norm_4 = twd_4 / mds_max,
      mds_norm = mds / mds_max
    ) |>
    ungroup()
}

# CHANGED: use shared plot builders inside the Rmd so the species panels and
# correlation plots stay in sync without copy-paste edits.
build_daily_species_stack <- function(species_name, title_text, rain_plot) {
  species_tc <- daily_tc_midday |>
    filter(species == species_name)

  species_sap <- daily_sap |>
    filter(species == species_name) |>
    select(date, sensor, total_sap)

  species_twd4 <- daily_twd4 |>
    filter(species == species_name) |>
    select(date, tree, twd_4)

  p_ts_dtc <- plot_ly(species_tc) |>
    add_lines(
      x = ~date,
      y = ~dtc_mean_midday,
      color = ~tree
    ) |>
    layout(
      yaxis = list(title = "midday avg ∆T (°C)"),
      showlegend = TRUE
    )

  p_ts_sap <- plot_ly(species_sap) |>
    add_lines(
      x = ~date,
      y = ~total_sap,
      color = ~sensor
    ) |>
    layout(
      yaxis = list(title = LAB_JS_DAILY_TXT),
      showlegend = FALSE
    )

  p_ts_twd <- plot_ly(species_twd4) |>
    add_lines(
      x = ~date,
      y = ~twd_4,
      color = ~tree
    ) |>
    layout(
      yaxis = list(title = "TWD4 (µm)"),
      showlegend = FALSE
    )

  subplot(
    p_ts_dtc,
    p_ts_sap,
    p_ts_twd,
    rain_plot,
    nrows = 4,
    shareX = TRUE,
    titleY = TRUE
  ) |>
    layout(
      title = title_text,
      xaxis = list(title = "Date")
    )
}

build_correlation_plot_list <- function(
  data,
  group_var,
  x_var,
  y_var,
  facet_formula,
  x_lab,
  y_lab,
  title_prefix,
  add_smooth = FALSE
) {
  group_values <- unique(data[[group_var]])
  plot_list <- setNames(vector("list", length(group_values)), group_values)

  for (group_value in group_values) {
    df_group <- data |>
      filter(.data[[group_var]] == group_value) |>
      filter(!is.na(.data[[x_var]]), !is.na(.data[[y_var]]))

    p <- ggplot(df_group, aes(x = .data[[x_var]], y = .data[[y_var]])) +
      geom_point(alpha = 0.7)

    if (isTRUE(add_smooth)) {
      p <- p + geom_smooth(method = "lm", se = TRUE)
    }

    plot_list[[as.character(group_value)]] <- p +
      stat_cor(
        aes(label = paste(..r.label.., ..p.label.., sep = "~")),
        method = "pearson",
        size = 3
      ) +
      facet_grid(facet_formula, scales = "free") +
      labs(
        title = paste(title_prefix, group_value),
        x = x_lab,
        y = y_lab
      ) +
      theme_bw() +
      theme(strip.text = element_text(size = 9))
  }

  plot_list
}

dir.create(here("output/research_days"), recursive = TRUE, showWarnings = FALSE)

Load Data

Time series

Looking at:

• ∆T

• stem radial growth (dendro) and max growth line to show how TWD is derived

• Total daily sap flow

Sap trees:

• P2

• J4 - ems & don’t have yet

Prep data

Syntax	Description
Header	Title
Paragraph	Text

# CHANGED: reuse a single P2 dendro object for the derivation and later
# time-series panels so the example stays internally consistent.
p_dendro_raw <- ggplot(dendro_dat, aes(x = ts, y = value, color = tree)) +
  geom_line() +
  labs(
    title = "Raw dendrometer values",
    y = "Stem radius displacement (µm)",
    x = "Timestamp"
  )

# to speed things up since plottly can be glitchy
# p_dendro_raw

# ggplotly(p_dendro_raw)

p2_dendro <- dendro_dat |>
  filter(tree == "p2") |>
  select(ts, value) |>
  arrange(ts) |>
  mutate(
    value = value + abs(first(value)),
    max = cummax(value),
    twd = max - value
  )

p_twd <- ggplot(p2_dendro, aes(x = ts)) +
  geom_line(aes(y = value)) +
  geom_line(aes(y = max), color = "red") +
  geom_line(aes(y = twd), color = "blue") +
  labs(
    title = "Radial growth and TWD components",
    y = "Stem radius (µm)",
    x = "Timestamp"
  )

p2_dt <- tc_dat |>
  transmute(ts, p2 = p2_dt)

p_dt <- ggplot(p2_dt, aes(ts, p2)) +
  geom_line()

# p2 <- p2_dt |>
#   left_join(
#     p2_dat_test |> select(ts, twd),
#     by = "ts"
#   ) |>
#   fill(twd, .direction = "downup")

# raw ∆T values
# p_a <- plot_ly(p2_dt, x = ~ts, y = ~p2, type = 'scatter', mode = 'lines',name = "p2 raw")

# CHANGED: reshape sap data once here so later chunks can reuse the same long
# and daily tables instead of rebuilding them with slightly different rules.
sap_long <- sap_dat |>
  pivot_longer(
    cols = -ts,
    names_to = "sensor",
    values_to = "sap_flow"
  ) |>
  mutate(
    date = as_date(ts),
    tree = tree_from_sensor(sensor),
    species = species_from_tree(tree)
  ) |>
  filter(!is.na(species))

sap_long |> 
  ggplot(aes(x = ts, y = sap_flow, colour = sensor))+
  geom_line() +
  facet_wrap(~species, ncol = 1)

## Warning: Removed 74994 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplotly(sap_long |> 
  ggplot(aes(x = ts, y = sap_flow, colour = sensor))+
  geom_line())

# Daily max & mean sap flow
daily_sap <- sap_long |>
  group_by(sensor, tree, species, date) |>
  summarise(
    total_sap = sum(sap_flow, na.rm = TRUE),
    .groups = "drop"
  )

total_sap_p <- ggplot(daily_sap, aes(x = date, y = total_sap, color = sensor)) +
  geom_line() +
  facet_wrap(~species, ncol = 1) +
  labs(
    title = "Total Daily Sap",
    y = LAB_JS_DAILY,
    x = "Date"
  )

# ggplotly(total_sap_p)

# to smooth data...
p2_dt <- p2_dt |>
  arrange(ts) |>
  mutate(
    p2_smooth = rollmean(p2, k = 5, fill = NA, align = "center")
  )

Using P2 for test run

# CHANGED: keep both 30-min and daily precip panels so sub-daily and daily
# figures share a matching time scale.
daily_precip <- flux |>
  mutate(date = as_date(ts)) |>
  group_by(date) |>
  summarise(
    precip = sum(precip, na.rm = TRUE),
    .groups = "drop"
  )

p_rain <- plot_ly(flux) |>
  add_lines(
    x = ~ts,
    y = ~precip,
    name = "precip",
    line = list(color = "blue")
  ) |>
  layout(
    yaxis = list(title = "precip (mm)")
  )

p_rain_daily <- plot_ly(daily_precip) |>
  add_lines(
    x = ~date,
    y = ~precip,
    name = "daily precip",
    line = list(color = "blue")
  ) |>
  layout(
    yaxis = list(title = "precip (mm)")
  )

# smoothed 2.5hr ∆T data
p_dt <- plot_ly(p2_dt) |>
  add_markers(
    x = ~ts,
    y = ~p2_smooth,
    name = "p2_smoothed ∆T",
    marker = list(size = 4, color = "black")
  ) |>
  layout(
    yaxis = list(title = "∆T (°C)")
  )

p_rad_growth <- plot_ly(
  p2_dendro,
  x = ~ts, y = ~value,
  type = "scatter", mode = "lines",
  name = "Stem radius"
) |>
  add_lines(
    x = ~ts,
    y = ~max,
    name = "Max radius",
    line = list(color = "red", width = 3)
  ) |>
  layout(
    yaxis = list(title = "Stem Radius (µm)")
  )

p_twd <- plot_ly(p2_dendro, x = ~ts, y = ~twd, type = "scatter", mode = "lines", name = "TWD", line = list(color = "blue")) |>
  layout(
    yaxis = list(title = "TWD (µm)")
  )

p_sap_tot <- plot_ly(daily_sap |> filter(tree == "p2"), x = ~date, y = ~total_sap, type = "scatter", mode = "lines", name = "Daily total Js", line = list(color = "blue")) |>
  layout(
    yaxis = list(title = LAB_JS_DAILY_TXT, size = 0.1)
  )

p_flow <- plot_ly(sap_long |> filter(tree == "p2"), x = ~ts, y = ~sap_flow, type = "scatter", mode = "lines", name = "Sap flux Js") |>
  layout(
    yaxis = list(title = LAB_JS_30MIN_TXT)
  )

fig <- subplot( p_dt, p_rad_growth, p_twd, p_sap_tot, p_flow, nrows = 5, shareX = TRUE, titleY = TRUE) |>
  layout(
    xaxis = list(title = "Time")
  )

## Warning: Ignoring 1907 observations

# fig

Dendro

Getting Dendrometer Stats

MDS — Maximum Daily Shrinkage (Blanco & Kalcsits, 2024)

We determined MDS of absolute SDV by calculating the differences between the daily pre-dawn maximum and the daily afternoon minimum in stem radius

• MDS = dailymax(radius) − daily min(radius)

• Hypothesis: Higher MDS -> more transpiration -> lower ∆T

TWD4 — Daily Tree Water Deficit (Dietrich et al.,2018)

• TWD4 = daily maximum TWD (the TWD value at the minimum stem radius) This version relates best to midday water potential (Ψ).

• Hypothesis: Higher MDS -> more transpiration -> lower ∆T

  • TWD4 increases with water stress and should be positively correlated with ∆T (Dietrich et al.,2018).

  • Because TWD related more strongly to midday Ψ than to MDS, midday ∆T is likely a better comparison

    Further, daily TWD was a better indicator for daily midday Ψ, particularly under dry conditions, than maximum daily shrinkage… (Dietrich et al., 2018)

TODO

• make sure the max is limited to the morning and the min is after, just so the max isn’t taken from the the end of the day (midnight)

• IDEA: maybe use how TWD differs from the previous day to make it daily instead of cumulative

How are TWD and TWD4 derived?

# CHANGED: reuse the same P2 derivation object created above instead of
# recalculating it here with slightly different inputs.
twd4 <- p2_dendro |> 
  mutate(date = as_date(ts)) |>
  group_by(date) |>
  filter(twd == max(twd, na.rm = TRUE)) |>
  slice(1) |>                 # in case of ties
  ungroup() |>
  select(ts, date, twd4 = twd)

p_rad_growth <- plot_ly(
  p2_dendro,
  x = ~ts, y = ~value,
  type = "scatter", mode = "lines",
  name = "Stem radius",
  line = list(color = "green", width = 3)
) |>
  add_trace(
    x = ~ts,
    y = ~max,
    name = "Max radius envelope",
    mode = "lines",
    fill = "tonexty",
    fillcolor = "lightblue",
    line = list(color = "red", width = 3)
  ) |>
  layout(
    yaxis = list(title = "Stem Radius (µm)")
  )

p_twd <- plot_ly(p2_dendro) |> 
  add_lines(x = ~ts, y = ~twd, name = "TWD", line = list(color = "lightblue")) |>
  add_markers(
    data = twd4,
    x = ~ts,
    y = ~twd4,
    name = "TWD4 (daily max)",
    marker = list(color = "black", size = 8, symbol = "circle")
  ) |> 
  layout(
    yaxis = list(title = "TWD (µm)")
  )

p_twd4 <- plot_ly(
  twd4,
  x = ~date,
  y = ~twd4,
  type = "scatter",
  mode = "lines",
  line = list(color = "black", width = 3),
  name = "TWD4"
) |>
  layout(
    yaxis = list(title = "TWD4 (µm)"),
    xaxis = list(title = "Date")
  )

fig <- subplot(p_rad_growth, p_twd,p_twd4, p_rain, nrows = 4, shareX = TRUE, titleY = TRUE) |>
  layout(
    xaxis = list(title = "Time")
  )

fig

How is MSD derived?

# CHANGED: switch the daily MDS derivation to the same predawn/midday windows
# used in twd_norm.Rmd so `mds`, `twd_norm_pd`, `twd_norm_4`, and `mds_norm`
# share one definition.
daily_mds <- calc_normalized_dendro_metrics(dendro_dat)

# Filtered datasets
df_p2 <- dendro_dat |>
  filter(tree == "p2") |>
  arrange(ts)

dm_p2 <- daily_mds |> filter(tree == "p2")

# Build plot
p_p2 <- plot_ly(df_p2) |>
  add_lines(
    x = ~ts,
    y = ~value,
    name = "Stem Radius",
    line = list(color = "green", width = 3)
  ) |> 
  # predawn max (red)
  add_markers(
    data = dm_p2,
    x = ~predawn_max_ts,
    y = ~predawn_max,
    marker = list(color = "red", size = 7),
    name = "Predawn max"
  ) |>
  # midday minimum (blue)
  add_markers(
    data = dm_p2,
    x = ~midday_min_ts,
    y = ~midday_min,
    marker = list(color = "blue", size = 7),
    name = "Midday min"
  ) |> 
  layout(
    title = "Diurnal Radius + MDS Points",
    yaxis = list(title = "Stem Radius (µm)")
  )

  # CHANGED: anchor the daily MDS trace on the midday minimum timestamp so it
  # lines up with the normalized predawn-to-midday window logic above.
p_mds <- plot_ly(data = dm_p2) |>
  add_lines(    
    x = ~midday_min_ts,
    y = ~mds,
    line = list(color = "purple", width = 3),
    name = "MDS") |> 
  layout(
    yaxis = list(title = "MDS (µm)")
  )

fig <- subplot(p_p2,p_mds, p_rain, nrows = 3, shareX = TRUE, titleY = TRUE) |>
  layout(
    xaxis = list(title = "Time")
  )

## Warning: Ignoring 1 observations

fig

TWD and MDS together

fig <- subplot(p_rad_growth, p_twd, p_twd4, p_mds, p_rain, nrows = 5, shareX = TRUE, titleY = TRUE) |>
  layout(
    xaxis = list(title = "Time")
  )

fig

# CHANGED: keep the normalized dendro outputs next to TWD4 here so later joins
# can use `twd_norm_pd`, `twd_norm_4`, and `mds_norm` without rebuilding the
# daily metrics.
dendro_daily <- dendro_dat |> 
  mutate(
    date = as_date(ts)
  )

# ---- DAILY TWD4 ----
# TWD4 = TWD at the DAILY MINIMUM radius (== daily maximum TWD)
# Because:
#   twd = max_envelope - radius
#   smallest radius -> largest TWD
daily_twd4 <- dendro_daily |>
  group_by(tree, species, date) |>
  summarise(
    twd_4 = if (all(is.na(twd))) NA_real_ else max(twd, na.rm = TRUE), # TWD4 described in https://doi.org/10.1093/treephys/tpy023
    .groups = "drop"
  )

# merge
daily_mds_twd <- daily_mds |>
  select(
    tree,
    species,
    date,
    daily_max_value = predawn_max,
    daily_min_value = midday_min,
    mds,
    twd_4,
    twd_norm_pd,
    twd_norm_4,
    mds_norm
  )

base_p <- ggplot(daily_mds_twd, aes(x = date))

twd_p <- base_p + geom_line(aes(y = twd_4, colour = tree)) +
  facet_wrap(~species, ncol = 1) +
  labs(
    title = "Max Daily TWD (TWD4)",
    y = "TWD (µm)",
    x = "Time"
  )

ggplotly(twd_p)

mds_p <- base_p + geom_line(aes(y = mds, colour = tree)) +
  facet_wrap(~species, ncol = 1) +
  labs(
    title = "Max Daily Shrinkage",
    y = "MDS (µm)",
    x = "Time"
  )

ggplotly(mds_p)

Sap Flow Metrics

Total Daily Sapflux (Snyder et al., 2024)

Note: I tested daily max and mean sap flow, but after discussing with Marcy and Will, I’m using total daily sap flow going forward.

TODO:

• Look into the Sap Flow Index (SFI) and how it compares to absolute sap flow measurements: https://doi.org/10.1093/treephys/19.13.885 - Decide whether max, mean, or SFI is most appropriate for comparison with ∆T and dendrometer metrics.
• Use Sap Total Daily Sap Flow (Done 11/10/2025)

# CHANGED: reuse the sap tables built in the prep chunk so this section only
# handles plotting and interpretation.

ggplotly(sap_long |> 
  ggplot(aes(x = ts, y = sap_flow, color = sensor)) +
  geom_line())

# ggplot code for total
p_total_sap <- ggplot(daily_sap, aes(x = date, y = total_sap, color = sensor)) +
  geom_line() +
  facet_wrap(~species, ncol = 1) +
  ggtitle(" sap") +
  labs(
    title = "Total Daily Sap",
    y = LAB_JS_DAILY,
    x = "Date"
  )

ggplotly(p_total_sap)

plots <- lapply(split(sap_long, sap_long$species), function(df_sp) {
  plot_ly(df_sp,
          x = ~ts,
          y = ~sap_flow,
          color = ~sensor,
          type = "scatter",
          mode = "lines") |>
    layout(
      yaxis = list(title = paste(unique(df_sp$species), LAB_JS_30MIN_TXT)),
      xaxis = list(title = "Date")
    )
})

# Combine into a column of subplots
p_total_sap <- subplot(plots,
                       nrows = length(plots),
                       shareX = TRUE,
                       titleY = TRUE)

fig <- subplot(
  p_total_sap,
  p_rain_daily,
  nrows = 2,
  shareX = TRUE,
  titleY = TRUE,
  heights = c(0.7, 0.3)   # larger height for sap, smaller for rain
)

fig

# smoothed 2.5hr ∆T data
p_a <- plot_ly(p2_dt) |>
  add_markers(
    x = ~ts,
    y = ~p2_smooth,
    name = "p2_smoothed",
    marker = list(size = 4, color = "black")
  ) |>
  layout(
    yaxis = list(title = "∆T (°C)")
  )

p_b <- plot_ly(sap_long |> filter(tree == "p2")) |>
  add_lines(
    x = ~ts,
    y = ~sap_flow,
    name = "Js"
  ) |>
  layout(
    yaxis = list(title = LAB_JS_30MIN_TXT)
  )

fig <- subplot(p_a, p_b, p_rain, nrows = 3, shareX = TRUE, titleY = TRUE) |>
  layout(
    xaxis = list(title = "Time")
  )

## Warning: Ignoring 1907 observations

fig

Canopy Temperature (∆T) Metrics

Midday window (10:00–14:00): Midday was chosen by visually inspecting the thermal time series and identifying when canopy temperatures consistently peaked and when shadow effects were minimal. This window best captures the period when evaporation cooling is most limited and water stress signals in ∆T are strongest.

Max midday ∆T: Represents the highest canopy temperature anomaly (Tc–Ta) during the midday period. This metric may reflect peak water stress or reduced transpiration when trees are most thermally exposed.

Mean midday ∆T: Calculated as the average Tc–Ta between 10:00 and 14:00. This smooths out short-term fluctuations and better represents sustained midday canopy water status.

# CHANGED: exclude non-tree panel ROIs up front so later joins stay focused on
# tree-level thermal signals.
tc <- tc_dat |>
  mutate(
    date = as_date(ts),
    hour = hour(ts)
  )

roi_cols <- tc |>
  select(ends_with("_dt"), -contains("panel")) |>
  colnames()

tc_long <- tc |>
  pivot_longer(
    cols = all_of(roi_cols),
    names_to = "roi_raw",
    values_to = "tc_ta"
  )

tc_long <- tc_long |>
  mutate(
    tree = str_remove(roi_raw, "_dt$"),
    species = species_from_tree(tree)
  ) |>
  filter(!is.na(species))

tc_midday <- tc_long |>
  filter(hour >= MIDDAY_START, hour <= MIDDAY_END)

daily_tc_midday <- tc_midday |>
  group_by(date, tree, species) |>
  summarise(
    dtc_max = max(tc_ta, na.rm = TRUE),
    dtc_mean_midday = mean(tc_ta, na.rm = TRUE),
    .groups = "drop"
  )

## Warning: There were 872 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `dtc_max = max(tc_ta, na.rm = TRUE)`.
## ℹ In group 1: `date = 2025-04-08`, `tree = "j20"`, `species = "juniper"`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 871 remaining warnings.

max_dtc_p <- ggplot(daily_tc_midday, aes(x = date, y = dtc_max, color = tree)) +
  geom_point(size = 0.7) +
  facet_wrap(~species, ncol = 1)

# ggplotly(max_dtc_p)

mean_dtc_p <- ggplot(daily_tc_midday, aes(x = date, y = dtc_mean_midday, color = tree)) +
  geom_point() +
  facet_wrap(~species, ncol = 1) +
  labs(
    title = "Average Midday ∆T",
    y = "∆T (°C)",
    x = "Date"
  )

ggplotly(mean_dtc_p)

Piñon: ∆T, Sap, dendro

# CHANGED: call the shared species plot builder here so the piñon/juniper
# sections stay identical except for the species filter.
build_daily_species_stack(
  species_name = "piñon",
  title_text = "Piñon: ∆T, Sap Flow, and TWD4",
  rain_plot = p_rain_daily
)

## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

Juniper: ∆T, Sap, dendro

# CHANGED: reuse the same plot builder as the piñon section so formatting only
# has to be updated in one place.
build_daily_species_stack(
  species_name = "juniper",
  title_text = "Juniper: ∆T, Sap Flow, and TWD4",
  rain_plot = p_rain_daily
)

## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

Pairwise response plots

Combining Daily Metrics

• daily_tc_dendro has midday ∆T values with daily MDS and TWD.
• daily_tc_sap has midday ∆T values with daily max and mean sap flow.

# CHANGED: keep missing normalized dendro values from dropping otherwise valid
# MDS/TWD4 days; the downstream plots filter only the columns they need.
daily_tc_dendro <- daily_tc_midday |>
  left_join(
    daily_mds_twd |> select(-c(daily_max_value, daily_min_value)),
    by = c("tree", "species", "date")
  )

daily_tc_sap <- daily_tc_midday |>
  left_join(
    daily_sap,
    by = c("tree", "species", "date")
  )|>
  na.omit()

# NEED TO WORK ON THIS: doesnt work
# daily_tc_dendro_sap <- daily_tc_midday |>
#   left_join(
#     daily_mds_twd |> select(-c(daily_max_value,daily_min_value)), by = c("tree", "species", "date")
#     ) |>
#     left_join(
#     daily_sap, by = c("tree", "species", "date")
#     )

∆T–Dendrometer Relationships

Midday canopy temperature metrics (∆T max and mean) were paired with daily dendrometer metrics (MDS and TWD) to assess how thermal stress relates to stem - level water status.

For every tree: -∆T variables were plotted against MDS and TWD with p and R values

TODO:

-separate by month

# CHANGED: use the shared plotting helper here so tree- and species-level
# dendro comparisons stay parallel and only use the cleaned midday ∆T metric.
pair_data <- daily_tc_dendro |>
  transmute(
    tree,
    species,
    date,
    dtc_var = "dtc_mean_midday",
    dtc = dtc_mean_midday,
    mds,
    twd_4,
    mds_norm,
    twd_norm_pd,
    twd_norm_4
  ) |>
  pivot_longer(
    cols = c(mds, twd_4, mds_norm, twd_norm_pd, twd_norm_4),
    names_to = "dendro_var",
    values_to = "dendro"
  ) |> 
  na.omit()

tree_plot_list <- build_correlation_plot_list(
  data = pair_data,
  group_var = "tree",
  x_var = "dtc",
  y_var = "dendro",
  facet_formula = dtc_var ~ dendro_var,
  x_lab = "°C",
  y_lab = "stem radius (µm)",
  title_prefix = "Tree:"
)

tree_plot_list

## Warning: The dot-dot notation (`..r.label..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(r.label)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

species_plot_list <- build_correlation_plot_list(
  data = pair_data,
  group_var = "species",
  x_var = "dtc",
  y_var = "dendro",
  facet_formula = dtc_var ~ dendro_var,
  x_lab = "°C",
  y_lab = "(µm)",
  title_prefix = "Species:"
)

species_plot_list

daily_tc_dendro <- daily_tc_dendro |> 
  filter(dtc_mean_midday>0)

p_daily_tc_twd4 <- daily_tc_dendro |>
  filter(!is.na(dtc_mean_midday), !is.na(twd_4)) |>
  ggplot(aes(x = dtc_mean_midday, y = twd_4, color = species)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) + # Add regression line
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson", size = 3
  ) +
  labs(
    title = "Relationship between Midday Canopy Temperature (ΔT) and Max Daily TWD.",
    x = "Midday Average ∆T (°C)",
    y = "Max Daily TWD (µm)",
    # caption = "Max Daily TWD (µm) shows a significant positive correlation with midday canooy temperature in the selected piñon and juniper. Shaded regions indicate 95% confidence intervals. Midday ΔT values were averaged from 10:00 to 16:00."
  ) +
  scale_color_manual(values = SPECIES_COLS) +
  theme(plot.caption = element_textbox_simple())

# CHANGED: split the normalized TWD view into predawn- and TWD4-based versions
# so both normalization choices are visible downstream.
p_daily_tc_norm_pd <- daily_tc_dendro |>
  filter(!is.na(dtc_mean_midday), !is.na(twd_norm_pd)) |>
  ggplot(aes(x = dtc_mean_midday, y = twd_norm_pd, color = species)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) + # Add regression line
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson", size = 3
  ) +
  labs(
    x = "Midday Average ∆T (°C)",
    y = "Normalized TWD_pd (-)") +
  scale_color_manual(values = SPECIES_COLS) +
  theme(plot.caption = element_textbox_simple())

p_daily_tc_norm_4 <- daily_tc_dendro |>
  filter(!is.na(dtc_mean_midday), !is.na(twd_norm_4)) |>
  ggplot(aes(x = dtc_mean_midday, y = twd_norm_4, color = species)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) + # Add regression line
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson", size = 3
  ) +
  labs(
    x = "Midday Average ∆T (°C)",
    y = "Normalized TWD_4 (-)") +
  scale_color_manual(values = SPECIES_COLS) +
  theme(plot.caption = element_textbox_simple())

p_daily_tc_mds <- daily_tc_dendro |> 
  filter(!is.na(dtc_mean_midday), !is.na(mds)) |>
  ggplot(aes(x = dtc_mean_midday, y = mds, color = species))+
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  # Add regression line
  stat_cor(
      aes(label = paste(..r.label.., ..p.label.., sep = "~")),
      method = "pearson", size = 3
    ) +
  scale_color_manual(values = SPECIES_COLS)+
  labs(
    x = "Midday Average ∆T (°C)",
    y = "MDS (µm)"
  )

p_daily_tc_mds_norm <- daily_tc_dendro |> 
  filter(!is.na(dtc_mean_midday), !is.na(mds)) |>
  ggplot(aes(x = dtc_mean_midday, y = mds_norm, color = species))+
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  # Add regression line
  stat_cor(
      aes(label = paste(..r.label.., ..p.label.., sep = "~")),
      method = "pearson", size = 3
    ) +
  scale_color_manual(values = SPECIES_COLS)+ 
  labs(
    x = "Midday Average ∆T (°C)",
    y = "Normalized MDS_pd (-)"
  )

ggsave(
  filename = here("output/research_days/p_daily_tc_twd4.png"),
  plot = p_daily_tc_twd4,
  width = 7,
  height = 5,
  units = "in"
)

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_twd4

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_norm_pd

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_norm_4

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_mds

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_mds_norm

## `geom_smooth()` using formula = 'y ~ x'

∆T–Sap Flow relationships

∆T max and mean vs daily sap flow metrics (daily max and mean sap flow) to examine how canopy thermal stress relates to tree - level water use. Data were reshaped into long format to generate all combinations of ∆T × sap flow variables.

For each **sap flow sensor** :-∆T values were plotted against sap flow with with p and R values

Note: tried separating by month to get more info on outliers. showed clusters, but got rid of correlations

# CHANGED: mirror the dendro workflow above so sap correlations use the same
# plotting pattern and cleaner daily inputs.
pair_data <- daily_tc_sap |>
  transmute(
    sensor,
    tree,
    species,
    date,
    dtc_var = "dtc_mean_midday",
    dtc = dtc_mean_midday,
    total_sap
  ) |>
  mutate(
    sap = total_sap, # directly assign total_sap
    sap_var = "total_sap", # single label for facet
    month = month(date)
  ) |> 
  filter(sap >0)

plot_list <- build_correlation_plot_list(
  data = pair_data,
  group_var = "sensor",
  x_var = "dtc",
  y_var = "sap",
  facet_formula = dtc_var ~ sap_var,
  x_lab = "Thermal signal (°C)",
  y_lab = LAB_JS_DAILY_TXT,
  title_prefix = "Sensor:",
  add_smooth = TRUE
)

plot_list

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

TREE_COLS <- c("p2" = "#40B0A6", "j4" = "#E1BE6A")


# CHANGED: keep the saved summary figure explicit so it also renders in the
# knitted notebook.
p_daily_tc_sap <- daily_tc_sap |> 
  filter(total_sap>0) |> 
  ggplot(aes(x = dtc_mean_midday, y = total_sap, color = tree)) +
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  stat_cor(
      aes(label = paste(..r.label.., ..p.label.., sep = "~")),
      method = "pearson", size = 3
    ) +
  scale_color_manual(values = TREE_COLS) +
  labs(
      #title = "Relationship between Midday Canopy Temperature (ΔT) and Daily Sap Flux.",
      x = "Midday Average ∆T (°C)",
      y = LAB_JS_DAILY,
      color = "Tree ID"
      # caption = " Daily total sap flux (Js) shows a significant positive correlation with midday thermal gradients in the selected Pinus edulis (P2) but not in selected Juniper monosperma (J4). Shaded regions indicate 95% confidence intervals. Sap Flux values were summed for the entire day. Midday ΔT values were averaged from 10:00 to 16:00."
    ) +
    theme(plot.caption = element_textbox_simple())

p_daily_tc_sap

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_sap_greaterThan0dtc_lm <- daily_tc_sap |> 
  filter(total_sap>0, dtc_mean_midday>0 ) |> 
  ggplot(aes(x = dtc_mean_midday, y = total_sap, color = tree)) +
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  stat_cor(
      aes(label = paste(..r.label.., ..p.label.., sep = "~")),
      method = "pearson", size = 3
    ) +
  scale_color_manual(values = TREE_COLS) +
  labs(
      #title = "Relationship between Midday Canopy Temperature (ΔT) and Daily Sap Flux.",
      x = "Midday Average ∆T (°C)",
      y = LAB_JS_DAILY,
      color = "Tree ID") +
    theme(plot.caption = element_textbox_simple())

p_daily_tc_sap_greaterThan0dtc_lm

## `geom_smooth()` using formula = 'y ~ x'

p_daily_tc_sap_greaterThan0dtc_gam <- daily_tc_sap |> 
  filter(total_sap>0, dtc_mean_midday>0 ) |> 
  ggplot(aes(x = dtc_mean_midday, y = total_sap, color = tree)) +
  geom_point()+
  geom_smooth(method = "gam", se = TRUE)+
  stat_cor(
      aes(label = paste(..r.label.., ..p.label.., sep = "~")),
      method = "pearson", size = 3
    ) +
  scale_color_manual(values = TREE_COLS) +
  labs(
      #title = "Relationship between Midday Canopy Temperature (ΔT) and Daily Sap Flux.",
      x = "Midday Average ∆T (°C)",
      y = LAB_JS_DAILY,
      color = "Tree ID") +
    theme(plot.caption = element_textbox_simple())



p_daily_tc_sap

## `geom_smooth()` using formula = 'y ~ x'

# # p_daily_tc_sap_greater has dtc less than zero linear model
ggsave(
  filename = here("output/research_days/p_daily_tc_sap.png"),
  plot = p_daily_tc_sap,
  width = 7,
  height = 5,
  units = "in"
)

## `geom_smooth()` using formula = 'y ~ x'

# 
# # p_daily_tc_sap_greaterThan0dtc linear model
# ggsave(
#   filename = here("output/research_days/p_daily_tc_sap_greaterThan0dtc_lm.png"),
#   plot = p_daily_tc_sap_greaterThan0dtc_lm,
#   width = 10,
#   height = 7,
#   units = "in",
#   dpi = 300
# )
# 
# # p_daily_tc_sap_greaterThan0dtc GAM
# ggsave(
#   filename = here("output/research_days/p_daily_tc_sap_greaterThan0dtc_gam.png"),
#   plot = p_daily_tc_sap_greaterThan0dtc_gam,
#   width = 10,
#   height = 7,
#   units = "in",
#   dpi = 300
# )

p_sap_twdNorm4 <- (p_daily_tc_sap_greaterThan0dtc_lm / p_daily_tc_norm_4) +
  plot_layout(axes = "collect") + 
  plot_annotation(tag_levels = 'a') & 
  labs(title = NULL)

# ggsave(
#   filename = here("output/research_days/p_sap_twdNorm4.png"),
#   plot = p_sap_twdNorm4,
#   width = 10,
#   height = 7,
#   units = "in",
#   dpi = 300
# )

Dirnal

Dayly sap and dt

# CHANGED: reuse the cleaned long sap table so this exploratory chunk still
# runs after the sap import was standardized higher up in the notebook.
test <- sap_long |>
  mutate(
    sap = sap_flow,
    hour = as_hms(ts)
  )

summary(test$sap)

p2_sap <- test |> 
  filter(tree == "p2") 

p2_dt_diurnal <- p2_dt |>
  mutate(hour = as_hms(ts))

p2_dt_diurnal |>
  ggplot(aes(x = hour, y = p2)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "loess", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 1372 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 1372 rows containing missing values or values outside the scale range
## (`geom_point()`).

p2_sap |> 
  ggplot(aes(x = hour, y = sap)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "loess", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 12625 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 12625 rows containing missing values or values outside the scale range
## (`geom_point()`).

test |> 
  ggplot(aes(x = hour, y = sap, color = tree)) +
  geom_point(alpha = 0.3) + 
  geom_smooth(method = "loess", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 87398 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 87398 rows containing missing values or values outside the scale range
## (`geom_point()`).

Species-level sap aggregation

The goal here is to see whether averaging sap and ∆T at the species level changes the relationship compared with looking at individual trees or sensors.

Visual comparison of piñon and juniper sap

# CHANGED: aggregate sap and midday ∆T by species using means rather than sums
# so the comparison is not driven by one species having more sensors/trees.
species_daily_sap <- daily_sap |>
  group_by(species, date) |>
  summarise(
    sap_species_mean = mean(total_sap, na.rm = TRUE),
    sd_sap = sd(total_sap, na.rm = TRUE),
    n_sensors = n_distinct(sensor), # Using n_distinct for the number of replicates
    se_sap = sd_sap / sqrt(n_sensors),
    .groups = "drop"
  )

species_daily_dtc <- daily_tc_midday |>
  filter(month(date)>3) |> 
  group_by(species, date) |>
  summarise(
    dtc_species_midday_mean = mean(dtc_mean_midday, na.rm = TRUE),
    sd_dtc = sd(dtc_mean_midday, na.rm = TRUE),
    n_trees = n_distinct(tree),
    se_dtc = sd_dtc / sqrt(n_trees),
    .groups = "drop"
  ) 


p_species_sap_compare <- species_daily_sap |>
  ggplot(aes(x = date, y = sap_species_mean, color = species, fill = species)) + # Added fill
  geom_ribbon(aes(ymin = sap_species_mean - se_sap, 
                  ymax = sap_species_mean + se_sap), 
              alpha = 0.2, color = NA) + # Shaded Error Area
  geom_line() +
  geom_point(size = 0.9) +
  scale_color_manual(values = SPECIES_COLS) +
  scale_fill_manual(values = SPECIES_COLS) + # Make ribbon match line color
  labs(title = "Species mean daily sap flow", x = "Date", y = LAB_JS_DAILY) +
  theme(legend.position = "right")

p_species_dtc_midday <- species_daily_dtc |>
  ggplot(aes(x = date, y = dtc_species_midday_mean, color = species, fill = species)) +
  geom_ribbon(aes(ymin = dtc_species_midday_mean - se_dtc, 
                  ymax = dtc_species_midday_mean + se_dtc), 
              alpha = 0.2, color = NA) +
  geom_line() +
  geom_point(size = 0.9) +
  scale_color_manual(values = SPECIES_COLS) +
  scale_fill_manual(values = SPECIES_COLS) +
  labs(title = "Species mean midday ∆T", x = "Date", y = "Midday average ∆T (°C)") +
  theme(legend.position = "right")

# Keep your patchwork layout as is
(p_species_sap_compare / p_species_dtc_midday) +
  plot_layout(guides = 'collect', axes = "collect")

## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).

Species diurnal sap and mean midday ∆T

# CHANGED: summarize both sap and dTc diurnally by species so the species-level
# comparison uses the same time-of-day frame and shows uncertainty on both.
species_sap_diurnal <- sap_long |>
  mutate(time_of_day = as_hms(format(ts, "%H:%M:%S"))) |>
  group_by(species, time_of_day) |>
  summarise(
    mean_sap = mean(sap_flow, na.rm = TRUE),
    sd_sap = sd(sap_flow, na.rm = TRUE),
    n = n(),
    se_sap = sd_sap / sqrt(n),
    .groups = "drop"
  )

species_dtc_diurnal <- tc_long |>
  mutate(time_of_day = as_hms(format(ts, "%H:%M:%S"))) |>
  group_by(species, time_of_day) |>
  summarise(
    mean_dtc = mean(tc_ta, na.rm = TRUE),
    sd_dtc = sd(tc_ta, na.rm = TRUE),
    n = n(),
    se_dtc = sd_dtc / sqrt(n),
    .groups = "drop"
  )

p_species_sap_diurnal <- species_sap_diurnal |>
  ggplot(aes(x = time_of_day, y = mean_sap, color = species, fill = species)) +
  geom_ribbon(
    aes(ymin = mean_sap - se_sap, ymax = mean_sap + se_sap),
    alpha = 0.2,
    color = NA
  ) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = SPECIES_COLS) +
  scale_fill_manual(values = SPECIES_COLS) +
  scale_x_time(
    breaks = hms::as_hms(sprintf("%02d:00:00", seq(0, 23, by = 2))),
    labels = scales::label_time("%H:%M")
  ) +
  labs(
    title = "Species mean diurnal sap flow",
    x = "Time of day",
    y = LAB_JS_30MIN
  ) 

p_species_dtc_diurnal <- species_dtc_diurnal |>
  ggplot(aes(x = time_of_day, y = mean_dtc, color = species, fill = species)) +
  geom_ribbon(
    aes(ymin = mean_dtc - se_dtc, ymax = mean_dtc + se_dtc),
    alpha = 0.2,
    color = NA
  ) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = SPECIES_COLS) +
  scale_fill_manual(values = SPECIES_COLS) +
  scale_x_time(
    breaks = hms::as_hms(sprintf("%02d:00:00", seq(0, 23, by = 2))),
    labels = scales::label_time("%H:%M")
  ) +
  labs(
    title = "Species mean diurnal ∆T",
    x = "Time of day",
    y = "∆T (°C)"
  ) 

# Make p_species_sap_diurnal share time of day range of p_species_dtc_diurnal

shared_limits <- c(min(species_dtc_diurnal$time_of_day), max(as_hms("19:00:00")))

# Combine and force the limits
(p_species_sap_diurnal / p_species_dtc_diurnal) +
  plot_layout(guides = 'collect', axes = "collect") & 
  coord_cartesian(xlim = shared_limits) & theme(legend.position = 'bottom')

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).

Species-level sap–∆T pairwise

# CHANGED: pair species-aggregated daily sap with species-aggregated mean
# midday ∆T to test whether piñon and juniper show different species-level
# sap-temperature relationships.
species_daily_sap_dtc <- species_daily_sap |>
  left_join(species_daily_dtc, by = c("species", "date")) |>
  filter(sap_species_mean > 0)

p_species_sap_dtc_pair <- species_daily_sap_dtc |>
  filter(dtc_species_midday_mean>=0) |> # i need to look into if i should do this. 
  ggplot(aes(
    x = dtc_species_midday_mean,
    y = sap_species_mean,
    color = species
  )) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "gam", se = TRUE) +
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson",
    size = 3
  ) +
  scale_color_manual(values = SPECIES_COLS) +
  labs(
    title = "Species-aggregated sap flow vs species mean midday ∆T",
    x = "Species mean midday ∆T (°C)",
    y = LAB_JS_DAILY
  ) +
  theme(legend.position = "none")

p_species_sap_dtc_pair

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

# ggsave(
#   filename = here("output/research_days/p_species_sap_dtc_gam.png"),
#   plot = p_species_sap_dtc_pair,
#   width = 7,
#   height = 5,
#   units = "in",
#   dpi = 300
# )

p_species_sap_dtc_pair <- species_daily_sap_dtc |>
    filter(dtc_species_midday_mean>=0) |> # i need to look into if i should do this. 
  ggplot(aes(
    y = dtc_species_midday_mean,
    x = sap_species_mean,
    color = species
  )) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "gam", se = TRUE) +
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson",
    size = 3
  ) +
  scale_color_manual(values = SPECIES_COLS) +
  labs(
    title = "Species-aggregated sap flow vs species mean midday ∆T",
    y = "Species mean midday ∆T (°C)",
    x = LAB_JS_DAILY
  ) +
  theme(legend.position = "none")

p_species_sap_dtc_pair

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

Species-level midday sap–∆T pairwise

# CHANGED: repeat the species-level sap–∆T comparison using sap summed only
# from 10:00 to 14:00 so the sap window matches the midday ∆T metric.
species_midday_sap <- sap_long |>
  filter(hour(ts) >= MIDDAY_START, hour(ts) <= MIDDAY_END) |>
  group_by(sensor, tree, species, date) |>
  summarise(
    midday_total_sap = sum(sap_flow, na.rm = TRUE),
    .groups = "drop"
  ) |>
  group_by(species, date) |>
  summarise(
    midday_sap_species_mean = mean(midday_total_sap, na.rm = TRUE),
    sd_midday_sap = sd(midday_total_sap, na.rm = TRUE),
    n_sensors = n_distinct(sensor),
    se_midday_sap = sd_midday_sap / sqrt(n_sensors),
    .groups = "drop"
  )

species_midday_sap_dtc <- species_midday_sap |>
  left_join(species_daily_dtc, by = c("species", "date")) |>
  filter(midday_sap_species_mean > 0)

p_species_midday_sap_dtc_pair <- species_midday_sap_dtc |>
  filter(dtc_species_midday_mean >= 0) |>
  ggplot(aes(
    x = dtc_species_midday_mean,
    y = midday_sap_species_mean,
    color = species
  )) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "gam", se = TRUE) +
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson",
    size = 3
  ) +
  scale_color_manual(values = SPECIES_COLS) +
  labs(
    title = "Species-aggregated midday sap flow vs species mean midday ∆T",
    x = "Species mean midday ∆T (°C)",
    y = "Midday total sap (10:00-14:00)"
  ) +
  theme(legend.position = "none")

p_species_midday_sap_dtc_pair

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

p_species_midday_sap_dtc_pair_flipped <- species_midday_sap_dtc |>
  filter(dtc_species_midday_mean >= 0) |>
  ggplot(aes(
    y = dtc_species_midday_mean,
    x = midday_sap_species_mean,
    color = species
  )) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "gam", se = TRUE) +
  stat_cor(
    aes(label = paste(..r.label.., ..p.label.., sep = "~")),
    method = "pearson",
    size = 3
  ) +
  scale_color_manual(values = SPECIES_COLS) +
  labs(
    title = "Species-aggregated midday sap flow vs species mean midday ∆T",
    y = "Species mean midday ∆T (°C)",
    x = "Midday total sap (10:00-14:00)"
  ) +
  theme(legend.position = "none")

p_species_midday_sap_dtc_pair_flipped

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'