2 Summarize raw fuel data

Code

library(dplyr)
library(tidyr)
library(purrr)
library(stringr)
library(forcats)
source("./scripts/process_datasheets.r")
source("./scripts/test_funs.r")

Before I convert our data to fuel loads, I would like to have at least a summary of the actual data that we collected in the field.

2.1 Load data from datasheets

Data in quasi-csv format, described in Section 1.1. I’m ommitting the medium density dispersed plot (MD), which only exists at one site.

Counts data is straight forward

Sampling cylinder data uses a mix of NA and 0 if there was no vegetation. I’m making the height 0 if there was no vegetation. Also, if there is no vegetation, species is “-”.

I’m converting to duff/litter depth and percent litter to duff depth and litter depth.

Code

transectid <- c("phase", "site", "treatment", "corner", "azi")

d <- combine_fuels_datasheets("../data") |>
  map(\(x) filter(x, treatment != "md"))

warn_duplicates(d$transects, phase, site, treatment, corner, azi)

counts <- d$transect |>
  select(all_of(transectid), slope, matches("(one|ten|hund)hr_count")) |>
  replace_na(list(slope = 0))

stations <- d$transect |>
  select(all_of(transectid), matches("\\w+[12]$")) |>
  pivot_longer(
    !all_of(transectid),
    names_to = c(".value", "station"),
    names_pattern = "(\\w+)([12])"
  ) |>
  # there was some inconsistency in whether heights were zero or blank if no veg
  # was present, here I sort that out.
  mutate(
    avg_w_ht = if_else(live_woody == 0 & dead_woody == 0, 0, avg_w_ht),
    avg_h_ht = if_else(live_herb == 0 & dead_herb == 0, 0, avg_h_ht),
    species = tolower(str_replace(species, "hb", "eh")),
    species = if_else(avg_w_ht == 0, "-", species),
    litter_depth = duff_litter * pct_litter / 100,
    duff_depth = duff_litter - litter_depth
  ) |>
  select(-c(metermark, duff_litter, pct_litter)) |>
  relocate(fbd, .after = last_col())

most_freq_string <- function(x) names(which.max(table(x)))

cwd <- d$coarse_woody

Heights and decay class observations are conditional on percent cover and diameter respectively, it might make more sense to report weighted averages here, but I’m not doing that.

For coarse woody debris, counts are averaged differently than diameter and decay class. I’ll report average diameter of observed particles only, and average of counts at the transect level (including zeros). I need to ensure I’m making zero counts explicit, as they are not in the data.

Code

sum_funcs <- list(
  mean = function(x) mean(x, na.rm = TRUE),
  se = function(x) mean(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))
)

stations_sum <- stations |>
  group_by(phase, treatment) |>
  summarise(across(c(where(is.double), -azi), sum_funcs)) |>
  pivot_longer(
    -c(treatment, phase),
    names_to = c("class", ".value"),
    names_pattern = "(.*)_(mean|se)"
  )

counts_sum <- counts |>
  group_by(phase, treatment) |>
  summarise(across(c(where(is.double), -azi), sum_funcs)) |>
  pivot_longer(
    -c(treatment, phase),
    names_to = c("class", ".value"),
    names_pattern = "(.*)_(mean|se)"
  )

# I need to know all possible transects to know where there are zeros for coarse
# woody
all_transects <- d$transects |>
  expand(nesting(phase, site, treatment, corner, azi))

# Average decay and size across all measured logs, average counts for all
# transects
cwd_sum <- cwd |>
  group_by(phase, treatment) |>
  mutate(across(c(dia, decay), sum_funcs)) |>
  group_by(phase, treatment, site, corner, azi) |>
  summarize(count = n(), across(matches("_(se|mean)$"), first)) |>
  right_join(all_transects) |>
  replace_na(list(count = 0)) |>
  group_by(phase, treatment) |>
  summarize(across(count, sum_funcs), across(matches("_(se|mean)$"), first)) |>
  pivot_longer(
    -c(treatment, phase),
    names_to = c("class", ".value"),
    names_pattern = "(.*)_(mean|se)"
  )

Here is a summary of all the data we collected.

Code

# this pivots wider and formats mean and se to print
stations_sum |>
  mutate(across(c(mean, se), \(x) format(round(x, 1)))) |>
  mutate(mean = glue::glue("{mean} ({se})", .trim = FALSE), .keep = "unused") |>
  pivot_wider(names_from = "class", values_from = "mean") |>
  mutate(
    treatment = fct_relevel(treatment, c("gs", "ld", "ha", "hd")),
    phase = factor(phase, levels = c("prepct", "postpct"))
  ) |>
  arrange(phase, treatment) |>
  rename(
    "Phase" = phase,
    "Treatment" = treatment,
    "Live Woody Cover (%)" = live_woody,
    "Dead Woody Cover (%)" = dead_woody,
    "Average Woody Height (m)" = avg_w_ht,
    "Live Herb Cover (%)" = live_herb,
    "Dead Herb Cover (%)" = dead_herb,
    "Average Herb Height (m)" = avg_h_ht,
    "Litter Depth (cm)" = litter_depth,
    "Duff Depth (cm)" = duff_depth,
    "Fuel Bed Depth (cm)" = fbd
  ) |>
  knitr::kable(align = c("l", "l", rep("r", 9)))

Table 2.1: Raw raw summary of data collected within sampling cylinders for each treatment (n = 64). Average heights are for live and dead vegetation particles, and don’t represent an average height across an area. They are meant to be multiplied by percent cover, but are presented indepenently here to document our actual measurements.

Phase	Treatment	Live Woody Cover (%)	Dead Woody Cover (%)	Average Woody Height (m)	Live Herb Cover (%)	Dead Herb Cover (%)	Average Herb Height (m)	Litter Depth (cm)	Duff Depth (cm)	Fuel Bed Depth (cm)
prepct	gs	38.4 (4.8)	6.2 (0.8)	3.5 (0.4)	6.2 (0.8)	1.4 (0.2)	0.6 (0.1)	2.8 (0.3)	3.5 (0.4)	8.9 (1.3)
prepct	ld	29.8 (3.7)	1.9 (0.2)	2.9 (0.4)	8.2 (1.0)	1.5 (0.2)	0.3 (0.0)	2.0 (0.2)	3.8 (0.5)	6.7 (1.0)
prepct	ha	23.2 (2.9)	1.7 (0.2)	2.0 (0.3)	5.7 (0.7)	2.1 (0.3)	0.6 (0.1)	2.6 (0.3)	2.6 (0.3)	12.5 (1.8)
prepct	hd	30.0 (3.8)	1.5 (0.2)	2.5 (0.3)	4.5 (0.6)	1.7 (0.2)	0.3 (0.0)	2.4 (0.3)	4.7 (0.6)	9.1 (1.3)
postpct	gs	8.3 (1.0)	0.7 (0.1)	1.0 (0.1)	5.5 (0.7)	0.9 (0.1)	0.2 (0.0)	20.0 (2.5)	2.2 (0.3)	35.7 (4.5)
postpct	ld	10.0 (1.3)	0.5 (0.1)	0.7 (0.1)	8.0 (1.0)	1.3 (0.2)	0.3 (0.0)	14.8 (1.8)	4.3 (0.5)	38.8 (4.9)
postpct	ha	10.4 (1.3)	0.8 (0.1)	1.0 (0.1)	7.8 (1.0)	0.6 (0.1)	0.3 (0.0)	9.1 (1.1)	1.5 (0.2)	22.0 (2.8)
postpct	hd	11.3 (1.4)	0.4 (0.1)	0.7 (0.1)	4.3 (0.5)	0.3 (0.0)	0.2 (0.0)	11.4 (1.4)	2.2 (0.3)	23.0 (2.9)

Code

counts_sum |>
  mutate(across(c(mean, se), \(x) format(round(x, 1)))) |>
  mutate(mean = glue::glue("{mean} ({se})", .trim = FALSE), .keep = "unused") |>
  pivot_wider(names_from = "class", values_from = "mean") |>
  mutate(
    treatment = fct_relevel(treatment, c("gs", "ld", "ha", "hd")),
    phase = factor(phase, levels = c("prepct", "postpct"))
  ) |>
  arrange(phase, treatment) |>
  rename(
    "Phase" = phase,
    "Treatment" = treatment,
    "Slope" = slope,
    "1-hr Count" = onehr_count,
    "10-hr Count" = tenhr_count,
    "100-hr Count" = hundhr_count
  ) |>
  knitr::kable(align = c("l", "l", rep("r", 4)))

Table 2.2: Raw raw summary of data collected along fuel transects for each treatment (n = 32). Transect lengths for 1-hr, 10-hr, and 100-hr fuels were 1, 2, and 4 m, respectively. Redwood leaf sprays were only counted if they were ⪆2 mm.

Phase	Treatment	Slope	1-hr Count	10-hr Count	100-hr Count
prepct	gs	29.3 (5.2)	15.5 (2.7)	7.8 (1.4)	4.4 (0.8)
prepct	ld	22.2 (3.9)	18.2 (3.2)	7.3 (1.3)	3.8 (0.7)
prepct	ha	31.2 (5.5)	31.5 (5.6)	6.6 (1.2)	3.8 (0.7)
prepct	hd	29.3 (5.2)	27.8 (4.9)	6.2 (1.1)	3.7 (0.6)
postpct	gs	25.9 (4.6)	49.2 (8.7)	19.6 (3.5)	8.5 (1.5)
postpct	ld	27.5 (4.9)	48.2 (8.5)	13.7 (2.4)	5.9 (1.0)
postpct	ha	29.9 (5.3)	24.2 (4.3)	8.2 (1.4)	4.3 (0.8)
postpct	hd	28.6 (5.1)	40.5 (7.2)	7.1 (1.2)	2.9 (0.5)

Code

cwd_sum |>
  mutate(across(c(mean, se), \(x) format(round(x, 1)))) |>
  mutate(mean = glue::glue("{mean} ({se})", .trim = FALSE), .keep = "unused") |>
  pivot_wider(names_from = "class", values_from = "mean") |>
  mutate(
    treatment = fct_relevel(treatment, c("gs", "ld", "ha", "hd")),
    phase = factor(phase, levels = c("prepct", "postpct"))
  ) |>
  arrange(phase, treatment) |>
  rename(
    "Phase" = phase,
    "Treatment" = treatment,
    "Count" = count,
    "Diameter (cm)" = dia,
    "Decay Class" = decay
  ) |>
  knitr::kable(align = c("l", "l", rep("r", 3)))

Table 2.3: Raw raw summary of data collected for coarse woody debris (> 7.62 cm) along 10-m transects for each treatment (n = 32). Diameter and decay class are for observed particles averaged at the treatment level (48 ≤ n ≤ 84). Counts are made at the transect level and averaged at the treatment level (n = 32).

Phase	Treatment	Count	Diameter (cm)	Decay Class
prepct	gs	2.4 (0.4)	17.2 (2.0)	3.7 (0.4)
prepct	ld	2.3 (0.4)	15.8 (1.8)	3.6 (0.4)
prepct	ha	1.7 (0.3)	17.3 (2.3)	3.8 (0.5)
prepct	hd	2.0 (0.3)	17.9 (2.3)	3.7 (0.5)
postpct	gs	2.6 (0.5)	16.8 (1.8)	2.7 (0.3)
postpct	ld	2.5 (0.4)	15.4 (1.7)	3.1 (0.4)
postpct	ha	1.5 (0.3)	17.3 (2.5)	3.7 (0.5)
postpct	hd	1.5 (0.3)	21.4 (3.1)	3.7 (0.5)

Code

save(d, file = "summarize_collected_data.rda")