In FoRTExperiment/fortedata: Data from the Forest Resilience and Threshold Experiment

knitr::opts_chunk$set(
  collapse = TRUE,
  require(fortedata),
  require(ggplot2),
  require(viridis),
  require(magrittr),
  require(dplyr))

Forest Inventory Background and Introduction

Each tree in FoRTE's 32 subplots above 8 cm in diameter-at-breast-height (DBH) is tagged, measured, and monitored over time. These data are compiled in the fd_inventory data set, described in detail below. In short, these data are traditional forest inventory data. They include:

• Tree DBH in cm (dbh_cm) measured as the average of two caliper measurements taken with a Haglof Postex PDII Inventory Unit and Caliper Set. Initial measurements were taken in 2018, with intention of remeasurement at 4 to 5 year intervals with the first remeasurmennt slated for 2022. This time interval is chosen based on previous experience at the University of Michigan Biological Station. Given the slow tree growth rates at the station's latitude, a 4 to 5 year measurement interval insures that growth will outpace measurment error.

• Canopy status (canopy_status) is recorded in this dataset where OD = overstory dominant, UN = understory, OS = overstory submissive, SA = sapling, and NA is a blank or missing record.

• Tree health status (health_status) where D = dead, M = moribund, and L = live

• Species (species) using the USDA Taxon system (e.g. FAGR is Fagus grandfolia, QURU is Quercus rubra, etc.).

• Tree individual identifier (tag) as the 4 digit tree tag number.

These data were originally collected in 2018, pre-disturbance. For stem-girdled tree information, i.e. which trees were targeted for mortality, see fd_mortatlity().

Data Availability

fortedata is an evolving, open-science data package with data updated in near-real time. The current inventory observations available as of r Sys.Date() are detailed in Figure 1. Remeasurement is slated to occur in the summer of 2022.

no_of_records.df <- fd_observations()

no_of_records <- subset(no_of_records.df, table == 'fd_inventory')


ggplot2::ggplot(no_of_records, ggplot2::aes(x = as.factor(month), y = as.integer(year), fill= no_of_obs)) +
  ggplot2::geom_tile(ggplot2::aes(fill = no_of_obs), color = "black") +
  ggplot2::geom_text(ggplot2::aes(label = no_of_obs), color = "white") +
  ggplot2::coord_equal()+
  ggplot2::scale_fill_gradient(low = "#450d54", high = "#450d54", na.value = 'white')+
  ggplot2::scale_y_reverse()+
  ggplot2::theme_minimal()+
  ggplot2::theme(legend.position = "none")+
  ggplot2::ylab("Year")+
  ggplot2::xlab("Month")+
  ggplot2::ggtitle(paste("Figure 1: No.forest inventory observations as of:", Sys.Date()))+
  ggplot2::facet_grid(table ~ .,  space = "free")+
  ggplot2::theme(strip.text.y = element_text(size = 9), strip.background = element_rect(
    color="black", fill="white", size= 0.5, linetype="solid"))

Forest Inventory Data

The fd_inventory() script within fortedata currently includes one function:

• fd_inventory() returns a single dataset of the forest inventory data, including diameter-at-breast height (DBH), latitude, longitude, and biomass for each measured stem, as well as information on vitality and canopy position. There are 3165 observations in the dataset, all measured in 2018.

fd_inventory()

inv <- data.frame(fd_inventory())


inv <- subset(inv, species != "????")
inv$species <- as.factor(inv$species)

# function to make limits
stat_box_data <- function(y, upper_limit = max(df$dbh_cm, na.rm = TRUE) * 1.15) {
  return(
    data.frame(
      y = 0.95 * upper_limit,
      label = paste('n =', length(y), sep = " ")
    )
  )
}

# good looking box plot with jitter plot on top
ggplot(inv, aes(x = species, y = dbh_cm, fill = species)) +
  geom_boxplot(outlier.size = 0) +
  scale_fill_viridis(discrete = TRUE, alpha=0.6) +
  #geom_jitter(color = "black", size = 0.5, alpha=0.2) +
  stat_summary(
    fun.data = stat_box_data,
    geom = "text",
    hjust = 0.5,
    vjust = 0.9,
    size = 3
  )+
  theme_bw()+
  theme(
    legend.position = "none",
    plot.title = element_text(size=11)
  )+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  xlab("")+
  ylab("Diameter at Breast Height [cm]")+
  facet_grid(replicate ~ species, scales = "free", space = "free")+
  theme(
    strip.background.x = element_blank(),
    strip.text.x = element_blank(),
    panel.spacing.x = unit(0, "lines"))+
  ggplot2::ggtitle(paste("Figure 2:  DBH by species, by replicate"))

Species Abbreviations for Figure 2:

ACPE - Acer pensylvanicum; ACRU - Acer rubrum; ACSA3 - Acer saccharum; AMELA - Amelanchier Medik.; BEAL2 - Betula alleghaniensis; BEPA - Betula papyrifera; FAGR - Fagus grandifolia; PIRE - Pinus resinosa; PIST - Pinus strobus; POGR4 - Populus grandidentata; POTR - Populus tremuloides; QURU - Quercus rubra; TSCA - Tsuga canadensis. Abbreviations are from the USDA PLANTS list.

Inventory Data

Biomass can be calculated from fd_inventory() by using the function calc_biomass().

# import the biomass data
x <- fortedata::calc_biomass()

# # return density plot of replicate biomass
# ggplot2::ggplot(x, ggplot2::aes(x = biomass, fill = as.factor(replicate)))+
#   ggplot2::geom_density(alpha = 0.4)+
#   ggplot2::theme(legend.position = "none")+
#   ggplot2::facet_grid(.~replicate)+
#     ggplot2::ggtitle(paste("Figure 3:  Aboveground Woody Biomass by Replicate"))
# 



# bring in metadata via the plot_metadata() function
df <- fortedata::fd_plot_metadata()

# now we convert the tibble to a data frame
df <- data.frame(df)

# First we want to concatenate our replicate, plot and subplot data to make a subplot_id column 
df$subplot_id <- paste(df$replicate, 0, df$plot, df$subplot, sep = "")
df$subplot_id <- as.factor(df$subplot_id)

# Now that we have our data in the form for this analysis, let's filter our metadata to the subplot level.
df %>%
  select(subplot_id, disturbance_severity, treatment) %>%
  distinct() %>%
  data.frame() -> dis.meta.data

# this filters the metadata down to the subplot_id level
dis.meta.data <- dis.meta.data[c(1:32), ]

# Then we merge with the metadata from above
x <- merge(x, dis.meta.data)


# For this analysis we want to code both disturbance severity and treatment as factors
x$disturbance_severity <- as.factor(x$disturbance_severity)
x$treatment <- as.factor(x$treatment)

# forte color palette
forte_pal <- forte_colors()

# first let's make some new, more informative labels for our facets
facet.labs <- c("B" = "Bottom-Up", "T" = "Top-Down")



x %>%
  group_by(subplot_id, disturbance_severity, treatment) %>%
  summarize(biomass_plot = sum(biomass)) -> y

y$biomass_mg_ha <- y$biomass_plot * 0.001 * 10

ggplot2::ggplot(y, aes(y = biomass_mg_ha, x = disturbance_severity, fill = disturbance_severity))+
  geom_boxplot(color = "black")+
  geom_jitter(position = position_jitter(0.2), shape = 21, alpha = 0.4, size = 4)+
  xlab("Disturbance Severity")+
  ylab(expression("Above-ground Woody Biomass [Mg per "~Ha^-1*"]"))+
  theme_minimal()+
  scale_color_manual(values = forte_pal, guide = FALSE)+
  scale_fill_manual(values = forte_pal,
                    name = "Disturbance Severity",
                    labels = c("0%", "45%", "65%", "85%"))+
  theme(legend.position = "bottom")+
  ggplot2::ggtitle(paste("Figure 3:  Aboveground Woody Biomass \n by Treatment, Disturbance Severity "))+
  facet_grid(. ~ treatment, labeller = labeller(treatment = facet.labs)) 

In figure 3 we show above-ground woody biomass at the plot level. This can be done by summing the biomass estimates to the plot level and then multiplying by 0.001 to convert kg to Mg and then multiplying that value by 10 to scale from 0.1 ha to 1 ha as follows:

# group the data by your factor of interest, here we show how the data are summed in Figure 3.
x %>%
  group_by(subplot_id, disturbance_severity, treatment) %>%
  summarize(biomass_plot = sum(biomass)) -> y

# Above, y$biomass is in units of kg per 0.1 ha. We want to convert to Mg per ha
y$biomass_mg_ha <- y$biomass_plot * 0.001 * 10

FoRTExperiment/fortedata documentation built on May 24, 2023, 4:51 a.m.