Set Up
Average Shrub Cover
calculate_uncovered_cm <- function(data, transect_length = 3000) {
# Filter out rows where start or end is NA or both are zero
data <- data %>%
filter(!is.na(Line_Start), !is.na(Line_End)) %>%
filter(!(Line_Start == 0 & Line_End == 0))
# If no valid shrubs, entire transect is uncovered
if (nrow(data) == 0) {
return(transect_length)
}
# Sort by Line_Start
data <- data[order(data$Line_Start, data$Line_End), ]
# Create merged list
merged <- list()
current <- data[1, ]
for (i in 2:nrow(data)) {
next_row <- data[i, ]
# If current and next overlap, merge them
if (!is.na(next_row$Line_Start) && !is.na(current$Line_End) &&
next_row$Line_Start <= current$Line_End) {
current$Line_End <- max(current$Line_End, next_row$Line_End, na.rm = TRUE)
} else {
merged <- append(merged, list(current))
current <- next_row
}
}
merged <- append(merged, list(current))
# Total covered length
total_covered <- sum(sapply(merged, function(x) x$Line_End - x$Line_Start), na.rm = TRUE)
# Uncovered length
uncovered <- transect_length - total_covered
return(uncovered)
}
no_shrub_cover <- shrub_data %>%
group_by(Plot) %>%
summarise(No_Shrub_Cover_cm = calculate_uncovered_cm(cur_data())) %>%
ungroup()
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `No_Shrub_Cover_cm = calculate_uncovered_cm(cur_data())`.
## ℹ In group 1: `Plot = "BI200"`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
Herbaceous Height
herbaceous_height <- fuel_data %>%
mutate(Height = as.numeric(Height)) %>% # Convert to numeric
filter(!is.na(Height)) %>% # Remove NA heights
group_by(Plot) %>%
summarise(Avg_Height = mean(Height, na.rm = TRUE)) %>%
ungroup()
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Height = as.numeric(Height)`.
## Caused by warning:
## ! NAs introduced by coercion
Herbacous Height by Status
# 25th, 50th, and 75th percentiles of height by invasion status
Herb_height_status <- herbaceous_height %>%
left_join(fuel_data %>% select(Plot, Status) %>% distinct(), by = "Plot") %>%
group_by(Status) %>%
summarise(
Height_25th = quantile(Avg_Height, 0.25, na.rm = TRUE),
Height_50th = quantile(Avg_Height, 0.50, na.rm = TRUE),
Height_75th = quantile(Avg_Height, 0.75, na.rm = TRUE)
) %>%
ungroup()
print(Herb_height_status)
## # A tibble: 2 × 4
## Status Height_25th Height_50th Height_75th
## <chr> <dbl> <dbl> <dbl>
## 1 Invaded 68.2 83.8 103.
## 2 Non_Invaded 15 21.3 33.7
Shrub Percentile
no_shrub_cover <- no_shrub_cover %>%
mutate(Percent_Shrub_Cover = 100 * (1 - No_Shrub_Cover_cm / 3000))
fuel_bed_data <- no_shrub_cover %>%
left_join(fuel_data %>% select(Plot, Status) %>% distinct(), by = "Plot")
# Shrub percentile by status
shrub_percentile <- fuel_bed_data %>%
group_by(Status) %>%
summarise(
Shrub_25th = quantile(Percent_Shrub_Cover, 0.25, na.rm = TRUE),
Shrub_50th = quantile(Percent_Shrub_Cover, 0.50, na.rm = TRUE),
Shrub_75th = quantile(Percent_Shrub_Cover, 0.75, na.rm = TRUE)
) %>%
ungroup()
print(shrub_percentile)
## # A tibble: 2 × 4
## Status Shrub_25th Shrub_50th Shrub_75th
## <chr> <dbl> <dbl> <dbl>
## 1 Invaded 25.7 47.2 77.6
## 2 Non_Invaded 39.3 72.4 92.1
Option 2: Binned Quantiles
herb_status_joined <- herbaceous_height %>%
left_join(fuel_data %>% select(Plot, Status) %>% distinct(), by = "Plot")
quantile_thresholds <- herb_status_joined %>%
group_by(Status) %>%
summarise(
Q25 = quantile(Avg_Height, 0.25, na.rm = TRUE),
Q50 = quantile(Avg_Height, 0.50, na.rm = TRUE),
Q75 = quantile(Avg_Height, 0.75, na.rm = TRUE),
.groups = "drop"
)
# Bin to nearest quantile
herbaceous_height_bins <- herb_status_joined %>%
left_join(quantile_thresholds, by = "Status") %>%
rowwise() %>%
mutate(
Height_Class = c("Q25", "Q50", "Q75")[which.min(
c(abs(Avg_Height - Q25), abs(Avg_Height - Q50), abs(Avg_Height - Q75))
)]
) %>%
ungroup() %>%
select(Plot, Status, Avg_Height, Height_Class)
print(herbaceous_height_bins)
## # A tibble: 159 × 4
## Plot Status Avg_Height Height_Class
## <chr> <chr> <dbl> <chr>
## 1 BI200 Invaded 64.3 Q25
## 2 BI201 Invaded 65.3 Q25
## 3 BI202 Invaded 80.3 Q50
## 4 BI97 Invaded 90 Q50
## 5 BI99 Invaded 68 Q25
## 6 BN210 Non_Invaded 13.3 Q25
## 7 BN211 Non_Invaded 15.3 Q25
## 8 BN212 Non_Invaded 27.3 Q50
## 9 BN96 Non_Invaded 33.7 Q75
## 10 BN98 Non_Invaded 8.33 Q25
## # ℹ 149 more rows
# Average height in each bin
average_height_bins <- herbaceous_height_bins %>%
group_by(Status, Height_Class) %>%
summarise(
Avg_Height = mean(Avg_Height, na.rm = TRUE),
.groups = "drop"
)
Shrub Cover by Height Bin
fuel_bed_Q <- no_shrub_cover %>%
left_join(herbaceous_height_bins %>% select(Plot, Status, Height_Class, Avg_Height) %>% distinct(), by = "Plot")
# Calculate average shrub cover per height bin (and status)
shrub_cover_by_height_bin <- fuel_bed_Q %>%
group_by(Status, Height_Class) %>%
summarise(
Avg_Height = mean(Avg_Height, na.rm = TRUE),
Avg_Shrub_Cover = mean(Percent_Shrub_Cover, na.rm = TRUE),
Count = n(),
.groups = "drop"
)
print(shrub_cover_by_height_bin)
## # A tibble: 6 × 5
## Status Height_Class Avg_Height Avg_Shrub_Cover Count
## <chr> <chr> <dbl> <dbl> <int>
## 1 Invaded Q25 59.8 66.9 27
## 2 Invaded Q50 84.1 49.8 21
## 3 Invaded Q75 110. 33.8 30
## 4 Non_Invaded Q25 11.6 65.7 30
## 5 Non_Invaded Q50 21.8 64.1 17
## 6 Non_Invaded Q75 40.9 64.2 34
kable(shrub_cover_by_height_bin, caption = "Average Shrub Cover by Height Bin and Status", format = "html", table.attr = "class='table table-striped'")
Average Shrub Cover by Height Bin and Status
|
Status
|
Height_Class
|
Avg_Height
|
Avg_Shrub_Cover
|
Count
|
|
Invaded
|
Q25
|
59.84691
|
66.91852
|
27
|
|
Invaded
|
Q50
|
84.12698
|
49.75397
|
21
|
|
Invaded
|
Q75
|
109.51111
|
33.79778
|
30
|
|
Non_Invaded
|
Q25
|
11.55000
|
65.73778
|
30
|
|
Non_Invaded
|
Q50
|
21.82353
|
64.05686
|
17
|
|
Non_Invaded
|
Q75
|
40.88235
|
64.15196
|
34
|