Summary Ancillary Data for SE-Svb

GRP_TREE GRP_AGB GRP_GAI GRP_LITTERPNT

Top

Data flow

This R-markdown document visualises ancillary data for forest ecosystems within the ICOS network. To come to this report, the data was collected and processed as follows:

1. Data collection following the instruction documents as published on the ICOS-ETC website. Data submission occurs via BADM to ETC.

2. ETC always performs data quality checks (e.g. missing data, outliers, consistency with instructions,… ).

3. Twice per year, ETC aggregates all data (fluxes, meteo, ancillary,…) collected within the ICOS network. For forest ancillary data, the general aggregation rule for all variable groups is to (1) average values within plots and (2) average these averages across plots. Statistics like standard deviation, min, max and percentiles are based on the values from (1) and thus represents the spatial heterogeneity among plots within the target area.
   These aggregations are published on the carbon portal as the L2 archive.

4. End users can unzip the downloaded L2 archive folder and find all .csv data products, as well as .pdf and BIF.csv documents describing the content of these products. The BIF ancillary file “ICOSETC_CC-Xxx_ANCILLARY_L2.csv” contains ALL ancillary data in a machine readable format.

5. ETC made an online BIF to BIFTAB converter available, which results in more digestable .csv files per variable group. Carefully follow the 4 steps by switching back and forth between the colab.research tab and the google.drive tab in your browser. These BIFTAB files are the direct input files for the graphs in this document.

6. Please cite and acknowledge ICOS if you make use of this data in any possible way.

GRP_TREE data at stand level

Species, DBH, tree height and health status is measured for all trees in CPs (every 3 years - 2000m²/CP) and SP-Is (every 10 years - 700m²/CP). Biomass is calculated for each individual tree via a species-specific allometric relation determined on site, from a nearby site or based on the literature. Biomass (tons/ha), basal area (m²wood/ha) and tree density (#/ha) at stand level is here represented for the most common species (contribute >5% to the total stand). Details on less common species can be found in the L2 archive.

library(readxl)
library(dplyr) 
library(openxlsx)
library(stringr)
library(ggplot2)
library(patchwork)
library(writexl)
library(tidyr)
# install.packages("knitr")
# install.packages("kableExtra")
library(knitr)
library(kableExtra)
library(akima)
library(viridis)
library(tinytex)
library(RColorBrewer)
library(lubridate) 
library(scales)
library(rmarkdown)

#### READ in data based on the BIFTAB files

path <- "I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/MOST_RECENT_BIFTAB_FOREST/"

read_latest_file <- function(path, pattern) {
  files <- list.files(path, pattern = pattern, full.names = TRUE)
  if (length(files) == 1) {
    return(read.csv(files[1]))
  } else if (length(files) > 1) {
    stop("More than one file following this pattern.")
  } else {
    stop("No such file found in this folder")
  }
}

add_columns <- function(df) {
  stat_num_col <- grep("STATISTIC_NUMBER$", colnames(df), value = TRUE)
  date_col <- grep("DATE$", colnames(df), value = TRUE)
  date_start_col <- grep("DATE_START$", colnames(df), value = TRUE)
  stat_col <- grep("STATISTIC$", colnames(df), value = TRUE)
  approach_col <- grep("APPROACH$", colnames(df), value = TRUE)
  spp_col <- grep("SPP$", colnames(df), value = TRUE)
  
  if (length(stat_col) > 0) {
    df <- df %>%
      filter(!!sym(stat_col) %in% c("Mean", "Standard Deviation"))
  }
  
  df <- df %>%
    mutate(
      Year = ifelse(!is.na(!!sym(date_col)), substr(!!sym(date_col), 1, 4), substr(!!sym(date_start_col), 1, 4)),
      Plot = ifelse(!!sym(stat_num_col) > 10, "SP", "CP"),
      PlotYear = paste0(Plot, "_", Year)
    )
  
  if (length(spp_col) > 0 && length(approach_col) > 0) {
    df <- df %>%
      mutate(!!sym(spp_col) := case_when(
        str_starts(!!sym(approach_col), "DEAD") ~ "Dead Standing Trees",
        str_starts(!!sym(approach_col), "TOTAL") ~ ".Total",
        str_starts(!!sym(approach_col), "ALL SPECIES") ~ ".Total",
        TRUE ~ !!sym(spp_col)
      ))
  }
  
  return(df)
}

safe_read_add_columns <- function(path, pattern) {
  tryCatch({
    add_columns(read_latest_file(path, pattern))
  }, error = function(e) {
    message(paste("No file for pattern:", pattern))
    return(NULL)
  })
}

SESvb_LAI <- safe_read_add_columns(path, "BIFTAB__GRP_LAI__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_LITTER_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_LITTER__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_BIOMASS_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_BIOMASS__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_DBH <- safe_read_add_columns(path, "BIFTAB__GRP_DBH__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_HEIGHT <- safe_read_add_columns(path, "BIFTAB__GRP_HEIGHTC__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_SPP_O_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_SPP_O__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_TREES_NUM_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_TREES_NUM__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")
SESvb_BASAL_AREA_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_BASAL_AREA__ICOSETC_SE-Svb_ANCILLARY__.*\\.csv")


### Select the most common species (>5% of species contribution)
SPP_Threshold <- SESvb_SPP_O_L2data %>%
  filter(SPP_O_STATISTIC == "Mean") %>%
  filter(!grepl("^DEAD", SPP_O_APPROACH))

SPP_Threshold2 <- SPP_Threshold %>%
  group_by(SPP_O_SPP) %>%
  summarize(mean_SPP_O = mean(SPP_O, na.rm = TRUE))

selected_species2 <- SPP_Threshold2 %>%
  filter(mean_SPP_O > 5) %>%
  pull(SPP_O_SPP)

dominant_species <- SPP_Threshold2 %>%
  filter(SPP_O_SPP != ".Total" & SPP_O_SPP != "") %>%
  slice_max(order_by = mean_SPP_O, n = 1) %>%  
  pull(SPP_O_SPP)

### Load latest version of the TREE database, which can be found at
# https://drive.google.com/drive/folders/1asrxu2TTWLH_QuzR944hvRoK7ZrEyc5l

FULL_TREE_DATABASE_GD <- read.xlsx("I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/DATABASE_exports/DATAFILE_Biomass_ICOS_GoogleDrive.xlsx")


TREEs_Site <- FULL_TREE_DATABASE_GD %>%
  filter(
    Site == "SE-Svb",
    !TREE_STATUS %in% c("Dead", "Removed", "Fallen")
  )


TREEs_Site_CPs <- TREEs_Site %>%
  filter(
    grepl("^CP", TREE_PLOT)
  ) %>%
  mutate(
    Year = as.numeric(Year),  # Zorg dat Year een numerieke kolom is
    Year = case_when(
      Year == 2019 ~ 2017,
      Year == 2018 ~ 2017,
      TRUE ~ Year
    )
  )

    TREEs_Site2 <- TREEs_Site %>%
  mutate(
    Year = as.numeric(Year),  # Zorg dat Year een numerieke kolom is
    Year = case_when(
      Year == 2019 ~ 2017,
      Year == 2018 ~ 2017,
      TRUE ~ Year
    )
  )
    
### Load coordinates of the plot locations
SESvb_Spatial <- read_excel("I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/DATABASE_exports/GRP_PLOT_AllSites.xlsx")
SESvb_Spatial2 <- SESvb_Spatial %>% 
  filter(Site == "SE-Svb")

transformed_data <- SESvb_BIOMASS_L2data %>%
  filter(grepl("^average|^TOTAL", BIOMASS_APPROACH, ignore.case = TRUE)) %>%

  mutate(
    Biomass_TonsPerHa = BIOMASS * 10, #from kg/m² in L2 file to tons/ha
    TREE_SPP = BIOMASS_SPP,
    Species = paste0("<i>", BIOMASS_SPP, "</i>")
  ) %>%
  select(
    TREE_SPP, Year, Plot, PlotYear, Biomass_TonsPerHa,
    BIOMASS_STATISTIC, BIOMASS_SPP, Species
  )

summarized_data <- transformed_data %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, BIOMASS_SPP, Species, BIOMASS_STATISTIC) %>%
  summarize(Biomass_TonsPerHa = mean(Biomass_TonsPerHa, na.rm = TRUE), .groups = 'drop')

final_data <- summarized_data %>%
  filter(BIOMASS_STATISTIC == "Mean") %>%
  select(-BIOMASS_STATISTIC) %>%
  rename(Biomass_TonsPerHa = Biomass_TonsPerHa) %>%
  left_join(
    summarized_data %>%
      filter(BIOMASS_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, BIOMASS_SPP, Species, Biomass_TonsPerHa) %>%
      rename(Biomass_TonsPerHa_SD = Biomass_TonsPerHa),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "BIOMASS_SPP", "Species")
  )


Pop_TREE_data_reduced2 <- final_data %>%
  filter(TREE_SPP %in% selected_species2)


ggplot(Pop_TREE_data_reduced2, aes(x = TREE_SPP, y = Biomass_TonsPerHa, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  
  scale_fill_brewer(palette = "Accent") +  
  geom_text(aes(label = sprintf("%.1f", Biomass_TonsPerHa)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, 
            size = 2.5,  
            color = "black") +
  geom_errorbar(aes(ymin = Biomass_TonsPerHa - Biomass_TonsPerHa_SD, 
                    ymax = Biomass_TonsPerHa + Biomass_TonsPerHa_SD), 
                position = position_dodge(width = 0.9), 
                width = 0.25, 
                color = "black") +
  labs(
    title = "Standing (dry) biomass at SE-Svb",
    subtitle = "only includes species that contribute >5% to the total tree count",  
    x = "Species",
    y = "Biomass (tons/ha)",
    fill = "PlotYear"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),  
    legend.position = "right"
  )

transformed_data2 <- SESvb_BASAL_AREA_L2data %>%
  filter(grepl("^average|^TOTAL", BASAL_AREA_APPROACH, ignore.case = TRUE)) %>%
  mutate(
    BasalArea_m2Ha = BASAL_AREA,
    TREE_SPP = BASAL_AREA_SPP,
    Species = paste0("<i>", BASAL_AREA_SPP, "</i>")
  ) %>%
  select(
    TREE_SPP, Year, Plot, PlotYear, BasalArea_m2Ha,
    BASAL_AREA_STATISTIC, BASAL_AREA_SPP, Species
  )

summarized_data2 <- transformed_data2 %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, BASAL_AREA_SPP, Species, BASAL_AREA_STATISTIC) %>%
  summarize(BasalArea_m2Ha = mean(BasalArea_m2Ha, na.rm = TRUE), .groups = 'drop')

final_data2 <- summarized_data2 %>%
  filter(BASAL_AREA_STATISTIC == "Mean") %>%
  select(-BASAL_AREA_STATISTIC) %>%
  left_join(
    summarized_data2 %>%
      filter(BASAL_AREA_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, BASAL_AREA_SPP, Species, BasalArea_m2Ha) %>%
      rename(BasalArea_m2Ha_SD = BasalArea_m2Ha),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "BASAL_AREA_SPP", "Species")
  )


Pop_TREE_data_reduced2ba <- final_data2 %>%
  filter(TREE_SPP %in% selected_species2)


ggplot(Pop_TREE_data_reduced2ba, aes(x = TREE_SPP, y = BasalArea_m2Ha, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  
  scale_fill_brewer(palette = "Accent") +  
  geom_text(aes(label = sprintf("%.1f", BasalArea_m2Ha)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, 
            size = 2.5,  
            color = "black") +
  geom_errorbar(aes(ymin = BasalArea_m2Ha - BasalArea_m2Ha_SD, 
                    ymax = BasalArea_m2Ha + BasalArea_m2Ha_SD), 
                position = position_dodge(width = 0.9), 
                width = 0.25, 
                color = "black") +
  labs(
    title = "Basal Area (m²/ha) at SE-Svb ",
    subtitle = "only includes species that contribute >5% to the total tree count",  
    x = "Species",
    y = "Basal Area (m²wood/ha)",
    fill = "PlotYear"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),  
    legend.position = "right"
  )

transformed_data3 <- SESvb_TREES_NUM_L2data %>%
  filter(grepl("^average|^TOTAL", TREES_NUM_APPROACH, ignore.case = TRUE)) %>%
  mutate(
    Density_TreesPerHa = TREES_NUM,
    TREE_SPP = TREES_NUM_SPP,
    Species = paste0("<i>", TREES_NUM_SPP, "</i>")
  ) %>%
  select(
    TREE_SPP, Year, Plot, PlotYear, Density_TreesPerHa,
    TREES_NUM_STATISTIC, TREES_NUM_SPP, Species
  )

summarized_data3 <- transformed_data3 %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, TREES_NUM_SPP, Species, TREES_NUM_STATISTIC) %>%
  summarize(Density_TreesPerHa = mean(Density_TreesPerHa, na.rm = TRUE), .groups = 'drop')


final_data3 <- summarized_data3 %>%
  filter(TREES_NUM_STATISTIC == "Mean") %>%
  select(-TREES_NUM_STATISTIC) %>%
  left_join(
    summarized_data3 %>%
      filter(TREES_NUM_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, TREES_NUM_SPP, Species, Density_TreesPerHa) %>%
      rename(Density_TreesPerHa_SD = Density_TreesPerHa),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "TREES_NUM_SPP", "Species")
  )

Pop_TREE_data_reduced3 <- final_data3 %>%
  filter(TREE_SPP %in% selected_species2)

ggplot(Pop_TREE_data_reduced3, aes(x = TREE_SPP, y = Density_TreesPerHa, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  
  scale_fill_brewer(palette = "Accent") +  
  geom_text(aes(label = sprintf("%.1f", Density_TreesPerHa)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, 
            size = 2.5,  
            color = "black") +
  geom_errorbar(aes(ymin = Density_TreesPerHa - Density_TreesPerHa_SD, 
                    ymax = Density_TreesPerHa + Density_TreesPerHa_SD), 
                position = position_dodge(width = 0.9), 
                width = 0.25, 
                color = "black") +
  labs(
      title = "Number of trees per hectare at SE-Svb ",
    subtitle = "only includes species that contribute >5% to the total tree count",  
    x = "Species",
    y = "Tree density (#/ha)",
    fill = "PlotYear"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),  
    legend.position = "right"
  )

SESvb_SPP_O_L2data_filtered <- SESvb_SPP_O_L2data %>%
  filter(
    SPP_O_SPP != "Dead Standing Trees" & 
      SPP_O_SPP != ".Total" & 
      SPP_O_STATISTIC != "Standard Deviation"  
    # SPP_O != 0
  )

SESvb_SPP_O_L2data_sorted <- SESvb_SPP_O_L2data_filtered %>%
  arrange(PlotYear, desc(SPP_O))

SESvb_SPP_O_L2data_sorted <- SESvb_SPP_O_L2data_sorted %>%
  group_by(PlotYear) %>%
  mutate(SPP_O_SPP = factor(SPP_O_SPP, levels = unique(SPP_O_SPP[order(SPP_O)])))  

unique_species <- unique(SESvb_SPP_O_L2data_sorted$SPP_O_SPP)
top_8_species <- unique_species[1:8]
remaining_species <- unique_species[9:length(unique_species)]

top_8_colors <- brewer.pal(n = 8, name = "Accent")
names(top_8_colors) <- top_8_species

remaining_colors <- viridis(length(remaining_species), option = "D")
names(remaining_colors) <- remaining_species

species_colors <- c(top_8_colors, remaining_colors)

ggplot(SESvb_SPP_O_L2data_sorted, aes(x = PlotYear, y = SPP_O, fill = SPP_O_SPP)) +
  geom_bar(stat = "identity", position = "stack", color = "black", size = 0.1) +  
  geom_text(aes(label = ifelse(SPP_O > 5, round(SPP_O, 1), "")), 
            position = position_stack(vjust = 0), 
            size = 3, 
            color = "black", 
            vjust = -0.5) +  
  scale_fill_manual(values = species_colors, name = "Species") +  
  labs(
    title = "Average tree species ratio per plot/campaign at SE-Svb",
    # x = "Plot Year",
    y = "mean tree species percentage (%)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), 
    axis.title.x = element_blank(), 
    axis.title.y = element_text(face = "bold"), 
    axis.text.x = element_text(face = "bold"),  
    axis.text.y = element_text(face = "bold"),  
    legend.title = element_text(face = "bold", hjust = 0.5) 
  )

############Frequency distribution for dominant species

Pop_TREE_data_dominant <- TREEs_Site_CPs %>%
  filter(TREE_SPP %in% dominant_species)

Pop_TREE_data_dominant <- Pop_TREE_data_dominant %>%
  mutate(Year = as.factor(Year))

ggplot(Pop_TREE_data_dominant, aes(x = TREE_DBH_CM, fill = Year, group = Year)) +
  geom_histogram(binwidth = 5, color = "black", position = "dodge", boundary = 0, alpha = 0.7) +
  geom_density(aes(y = after_stat(count) * 5, color = Year, fill = Year), linewidth = 0.5, alpha = 0, adjust = 1.5, position = "identity") +
  labs(
    title = "DBH distribution per observation year within CPs at SE-Svb",
    subtitle = bquote("for the most common species: " ~ italic(.(dominant_species))),
    x = "DBH (cm)", y = "Tree count"
  ) +
  scale_x_continuous(breaks = seq(0, max(Pop_TREE_data_dominant$TREE_DBH_CM, na.rm = TRUE), by = 5)) +
  scale_fill_brewer(palette = "Accent") +
  scale_color_brewer(palette = "Accent") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), 
    plot.subtitle = element_text(hjust = 0.5, size = 10),  
    axis.title.x = element_text(face = "bold"),  
    axis.title.y = element_text(face = "bold"), 
    axis.text.x = element_text(face = "bold"), 
    axis.text.y = element_text(face = "bold"),  
    legend.title = element_text(face = "bold"),  
    legend.position = "right", 
    strip.text = element_text(face = "bold"),
    panel.grid.minor.x = element_blank()   
  )

##### Spatial component Biomass ####

# Average Biomass_KG per TREE_PLOT and Year
Biomass_PLOT <- TREEs_Site2 %>%
  group_by(TREE_PLOT, Year) %>%
  summarize(Sum_Biomass_Ton = round(sum(Biomass_KG, na.rm = TRUE) / 1000, 2)) %>%
  mutate(
    Area_ha = ifelse(startsWith(TREE_PLOT, "CP"), 0.2, 0.07),
    Biomass_TonsPerHa = round(Sum_Biomass_Ton / Area_ha, 2),
    PlotYear = paste(TREE_PLOT, Year, sep = "_")
  )

SESvb_Spatial_filtered <- SESvb_Spatial2 %>%
  select(PLOT_ID, PLOT_LOCATION_LAT, PLOT_LOCATION_LONG)
Biomass_PLOT <- Biomass_PLOT %>%
  left_join(SESvb_Spatial_filtered, by = c("TREE_PLOT" = "PLOT_ID"))

#Double check: plot SP-I_02 turned into CP_06 during labelling. so it can be ignored here
#Biomass_PLOT2 <- Biomass_PLOT %>%
# filter(TREE_PLOT != "SP-I_02")

# Double check: Filter de data voor Year == year of labelling - this is very site specific
Biomass_PLOT_labelling <- Biomass_PLOT %>%
  filter(Year == "2017")

interp_data <- interp(
  x = Biomass_PLOT_labelling$PLOT_LOCATION_LONG,  
  y = Biomass_PLOT_labelling$PLOT_LOCATION_LAT,   
  z = Biomass_PLOT_labelling$Biomass_TonsPerHa,
  xo = seq(min(Biomass_PLOT_labelling$PLOT_LOCATION_LONG), max(Biomass_PLOT_labelling$PLOT_LOCATION_LONG), length = 100),
  yo = seq(min(Biomass_PLOT_labelling$PLOT_LOCATION_LAT), max(Biomass_PLOT_labelling$PLOT_LOCATION_LAT), length = 100),
  linear = TRUE,
  duplicate = "mean"  
)

# Turn this into a dataframe
interp_df <- data.frame(
  x = rep(interp_data$x, times = length(interp_data$y)),
  y = rep(interp_data$y, each = length(interp_data$x)),
  z = as.vector(interp_data$z)
)

Biomass_PLOT_labelling_clean <- Biomass_PLOT_labelling %>%
  mutate(
    plot_type = factor(ifelse(startsWith(TREE_PLOT, "CP"), "CP plot", "SP-I plot"))
  )

ggplot() +
  geom_raster(data = interp_df, aes(x = x, y = y, fill = z)) +  
  geom_contour(data = interp_df, aes(x = x, y = y, z = z), bins = 10, color = "white") + 
  geom_point(data = Biomass_PLOT_labelling_clean, aes(x = PLOT_LOCATION_LONG, y = PLOT_LOCATION_LAT, shape = plot_type, color = plot_type), size = 3) +  
  geom_text(data = Biomass_PLOT_labelling_clean, aes(x = PLOT_LOCATION_LONG, y = PLOT_LOCATION_LAT, label = TREE_PLOT), vjust = -0.5, hjust = 1.5, size = 3, check_overlap = TRUE) +  
  scale_fill_viridis_c(name = "Biomass (tons/ha)") +  
  scale_color_manual(values = c("CP plot" = "red", "SP-I plot" = "pink"), name = "Plot Type") +  
  scale_shape_manual(values = c("CP plot" = 15, "SP-I plot" = 16), name = "Plot Type") +  
  labs(
    title = "Spatial variation in Biomass (tons/ha) at SE-Svb",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold"),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(face = "bold"),
    legend.position = "right"  
  )

GRP_TREE data at individual level

DBH and tree height at individual level is here represented for the most common species (contribute >5% to the total stand).

transformed_data5 <- SESvb_DBH %>%
  filter(grepl("^average", DBH_APPROACH, ignore.case = TRUE)) %>%
  mutate(
    TREE_SPP = DBH_SPP,
    Species = paste0("<i>", DBH_SPP, "</i>")
  ) %>%
  select(
    TREE_SPP, Year, Plot, PlotYear, DBH,
    DBH_STATISTIC, DBH_SPP, Species
  )

summarized_data5 <- transformed_data5 %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, DBH_SPP, Species, DBH_STATISTIC) %>%
  summarize(DBH = mean(DBH, na.rm = TRUE), .groups = 'drop')

final_data5 <- summarized_data5 %>%
  filter(DBH_STATISTIC == "Mean") %>%
  select(-DBH_STATISTIC) %>%
  left_join(
    summarized_data5 %>%
      filter(DBH_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, DBH_SPP, Species, DBH) %>%
      rename(DBH_SD = DBH),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "DBH_SPP", "Species")
  )


Pop_TREE_data_reduced5 <- final_data5 %>%
  filter(TREE_SPP %in% selected_species2)

ggplot(Pop_TREE_data_reduced5, aes(x = TREE_SPP, y = DBH, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  
  scale_fill_brewer(palette = "Accent") +  
  geom_text(aes(label = sprintf("%.1f", DBH)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, 
            size = 2.5,  
            color = "black") +
  geom_errorbar(aes(ymin = DBH - DBH_SD, 
                    ymax = DBH + DBH_SD), 
                position = position_dodge(width = 0.9), 
                width = 0.25, 
                color = "black") +
  labs(
    title = "Average tree DBH per species and plot/campaign at SE-Svb",
    subtitle = "only includes species that contribute >5% to the total tree count",  
    x = "Species",
    y = "mean +/- SD tree DBH (cm)",
    fill = "Plot Year"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),  
    legend.position = "right"
  )

#####Height

transformed_data4 <- SESvb_HEIGHT %>%
  filter(grepl("^average", HEIGHTC_APPROACH, ignore.case = TRUE)) %>%
  mutate(
    TREE_SPP = HEIGHTC_SPP,
    Species = paste0("<i>", HEIGHTC_SPP, "</i>")
  ) %>%
  select(
    TREE_SPP, Year, Plot, PlotYear, HEIGHTC,
    HEIGHTC_STATISTIC, HEIGHTC_SPP, Species
  )

summarized_data4 <- transformed_data4 %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, HEIGHTC_SPP, Species, HEIGHTC_STATISTIC) %>%
  summarize(HEIGHTC = mean(HEIGHTC, na.rm = TRUE), .groups = 'drop')

final_data4 <- summarized_data4 %>%
  filter(HEIGHTC_STATISTIC == "Mean") %>%
  select(-HEIGHTC_STATISTIC) %>%
  left_join(
    summarized_data4 %>%
      filter(HEIGHTC_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, HEIGHTC_SPP, Species, HEIGHTC) %>%
      rename(HEIGHTC_SD = HEIGHTC),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "HEIGHTC_SPP", "Species")
  )

Pop_TREE_data_reduced4 <- final_data4 %>%
  filter(TREE_SPP %in% selected_species2)

ggplot(Pop_TREE_data_reduced4, aes(x = TREE_SPP, y = HEIGHTC, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  
  scale_fill_brewer(palette = "Accent") +  
  geom_text(aes(label = sprintf("%.1f", HEIGHTC)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, 
            size = 2.5,  
            color = "black") +
  geom_errorbar(aes(ymin = HEIGHTC - HEIGHTC_SD, 
                    ymax = HEIGHTC + HEIGHTC_SD), 
                position = position_dodge(width = 0.9), 
                width = 0.25, 
                color = "black") +
  labs(
    title = "Average tree height per species and plot/campaign at SE-Svb",
    subtitle = "only includes species that contribute >5% to the total tree count",  
    x = "Species",
    y = "mean +/- SD tree height (m)",
    fill = "Plot Year"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),  
    legend.position = "right"
  )

Extra - Dead standing trees

Please note that standing dead trees were not consistently monitored during the labelling phase of the site. Dead standing trees were consistently measured in remeasurements campaigns from 2023 onwards.

SESvb_BIOMASS_DeadTrees <- SESvb_BIOMASS_L2data %>% 
  filter(startsWith(BIOMASS_APPROACH, "DEAD")) %>%
  mutate(
    Biomass_TonsPerHa = BIOMASS * 10,  # from kg/m² to ton/ha
    TREE_SPP = BIOMASS_SPP,
    Species = paste0("<i>", BIOMASS_SPP, "</i>")
  ) %>%
  select(TREE_SPP, Year, Plot, PlotYear, Biomass_TonsPerHa, BIOMASS_STATISTIC, BIOMASS_SPP, Species)

summarized_data <- SESvb_BIOMASS_DeadTrees %>%
  group_by(TREE_SPP, Year, Plot, PlotYear, BIOMASS_SPP, Species, BIOMASS_STATISTIC) %>%
  summarize(Biomass_TonsPerHa = mean(Biomass_TonsPerHa, na.rm = TRUE), .groups = 'drop')

final_data <- summarized_data %>%
  filter(BIOMASS_STATISTIC == "Mean") %>%
  select(-BIOMASS_STATISTIC) %>%
  rename(Biomass_TonsPerHa = Biomass_TonsPerHa) %>%
  left_join(
    summarized_data %>%
      filter(BIOMASS_STATISTIC == "Standard Deviation") %>%
      select(TREE_SPP, Year, Plot, PlotYear, BIOMASS_SPP, Species, Biomass_TonsPerHa) %>%
      rename(Biomass_TonsPerHa_SD = Biomass_TonsPerHa),
    by = c("TREE_SPP", "Year", "Plot", "PlotYear", "BIOMASS_SPP", "Species")
  )

# 🔥 **Make graph** 🔥
ggplot(final_data, aes(x = TREE_SPP, y = Biomass_TonsPerHa, fill = PlotYear)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  scale_fill_brewer(palette = "Accent") +
  geom_text(
    aes(label = sprintf("%.1f", Biomass_TonsPerHa)),
    position = position_dodge(width = 0.9),
    vjust = -0.5,
    size = 2.5,
    color = "black"
  ) +
  geom_errorbar(
    aes(ymin = Biomass_TonsPerHa - Biomass_TonsPerHa_SD, ymax = Biomass_TonsPerHa + Biomass_TonsPerHa_SD),
    position = position_dodge(width = 0.9),
    width = 0.25,
    color = "black"
  ) +
  labs(
    title = "Biomass standing dead trees at SE-Svb",
    x = "", 
    y = "Biomass (tons/ha)",
    fill = "PlotYear"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, face = "plain"), 
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 10, hjust = 0.5),
    legend.position = "right"
  )

GRP_AGB - Understory

In addition to the TREE data (DBH >5cm), above ground biomass is determined for 5 understory plant types (moss, herbs, ferns, shrubs and saplings (DBH <5cm)). As detailed in the instruction document, this is done at 5 specific locations per CP.

filtered_data <- SESvb_BIOMASS_L2data %>%
  filter(grepl("^AGB understory", BIOMASS_APPROACH) & BIOMASS_STATISTIC == "Mean") %>%
  mutate(
    AGB_VEGTYPE = case_when(
      BIOMASS_VEGTYPE == "Herb" ~ "Ferns",
      BIOMASS_VEGTYPE == "Annual Herb" ~ "Herb",
      BIOMASS_VEGTYPE == "Non-vascular" ~ "Moss",
      BIOMASS_VEGTYPE == "Tree" ~ "Sapling",
      BIOMASS_VEGTYPE == "Understory" ~ NA_character_,
      TRUE ~ BIOMASS_VEGTYPE
    )
  ) %>%
  filter(!is.na(AGB_VEGTYPE))

date_n_mapping <- filtered_data %>%
  distinct(Year, BIOMASS_STATISTIC_NUMBER)

plot <- ggplot(filtered_data, aes(x = as.factor(Year), y = BIOMASS, fill = AGB_VEGTYPE)) +
  geom_bar(stat = "identity", color = "black") +
  scale_fill_brewer(palette = "Accent") +  
  labs(
    x = "Year",
    y = "AGB understory (kg/m²)",
    title = "Annual variation in AGB understory at SE-Svb",
    subtitle = "(n = number of CPs)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(labels = function(x) {
    n_value <- date_n_mapping$BIOMASS_STATISTIC_NUMBER[match(x, date_n_mapping$Year)]
    paste0(x, "\n(n = ", n_value, ") ")
  })

print(plot)

GRP_GAI

Green Area Index is estimated either with digital hemispherical photography (DHP) at 9 standard locations per CP or with a ceptometer in 25 locations per CP. Values are first averaged within plots and then averaged across plots to come to a single value per observation date.

SESvb_LAI$DATE2 <- ifelse(is.na(SESvb_LAI$LAI_DATE), SESvb_LAI$LAI_DATE_START, SESvb_LAI$LAI_DATE)

#only data from CP - thus excluding SP-I DHPs from labelling
SESvb_LAI <- SESvb_LAI %>%
  filter(LAI_APPROACH != "Up to 7 DHP pictures were taken in each sparse plot (SP). LAI was calculated for each DHP passing the quality check and then averaged per plot if at least 70% of DHP was available. Values reported are the average and SD across the SP-Is used (number reported in the STATISTIC_NUMBER variable). For more info visit www.icos-etc.eu.")

SESvb_LAI_mean <- SESvb_LAI %>% filter(LAI_STATISTIC == "Mean")
SESvb_LAI_sd <- SESvb_LAI %>% filter(LAI_STATISTIC == "Standard Deviation")

combined_data <- SESvb_LAI_mean %>%
  rename(Mean_LAI = LAI) %>%
  select(DATE2, Mean_LAI) %>%
  inner_join(SESvb_LAI_sd %>% rename(SD_LAI = LAI) %>% select(DATE2, SD_LAI), by = "DATE2")

combined_data$DATE2 <- as.Date(as.character(combined_data$DATE2), format = "%Y%m%d")
min_date <- min(combined_data$DATE2, na.rm = TRUE)
start_year <- as.numeric(format(min_date, "%Y"))
start_date <- as.Date(paste0(start_year, "-01-01"), format = "%Y-%m-%d")

plot <- ggplot(combined_data, aes(x = DATE2)) +
  geom_point(aes(y = Mean_LAI), color = "blue", size = 2.5, position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = Mean_LAI - SD_LAI, ymax = Mean_LAI + SD_LAI), width = 0.2, color = "blue", position = position_dodge(width = 0.2), color = "black") +
  #geom_line(aes(y = Mean_LAI), color = "grey", group = 1) + 
  labs(x = "Date", y = "LAI +/- SD (m²/m²)", title = "Seasonal variation LAI at SE-Svb") +
  scale_x_date(
    date_labels = "%d-%m-%Y",
    breaks = seq(start_date, max(combined_data$DATE2), by = "1 year"), # Ticks every year
    minor_breaks = seq(start_date, max(combined_data$DATE2), by = "3 months") # Minor ticks every 3 months
  ) +
  ylim(0, NA) +  # Set the lower limit of the y-axis to 0
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, face = "bold"), 
    axis.text.y = element_text(face = "bold"), 
    axis.title = element_text(face = "bold"), 
    plot.title = element_text(face = "bold", hjust = 0.5),  # Center the plot title
    legend.position = "bottom" 
  )

print(plot)

GRP_LITTERPNT

Three types of litter are collected in class1 sites: coarse woody debris (diameter >10cm that cross one of the three transects per CP - 1x per year), fine woody debris (diameter >2cm and <10cm in 5x1m² subplots per CP - 1x per year) and non-woody litter, sorted per fraction (min. 5 litter traps per CP - up to 16x per year).

SESvb_LITTER <- SESvb_LITTER_L2data %>%
  mutate(LITTER_DATE = as.Date(as.character(LITTER_DATE), format = "%Y%m%d"))
#Double check: Site specific: exclude Year because it is not a complete year.
#SESvb_LITTER_filtered <- SESvb_LITTER %>%
#  filter(format(LITTER_DATE, "%Y") != "2022")



filtered_data <- SESvb_LITTER %>% ####Double check: adjust if needed
  filter(LITTER_STATISTIC %in% c("Mean", "Standard Deviation") &
           LITTER_ORGAN %in% c("Total", "Branches", "Trunks")) %>%
  mutate(
    Year = ifelse(!is.na(LITTER_DATE), substr(LITTER_DATE, 1, 4), substr(LITTER_DATE_START, 1, 4))
  ) %>%
  pivot_wider(
    names_from = LITTER_STATISTIC,
    values_from = LITTER,
    names_prefix = "LITTER_"
  )


aggregated_data <- filtered_data %>%
  group_by(Year) %>%
  summarise(
    NonWoody = sum(ifelse(LITTER_ORGAN == "Total", LITTER_Mean, 0), na.rm = TRUE),
    NonWoody_SD = sum(ifelse(LITTER_ORGAN == "Total", `LITTER_Standard Deviation`, 0), na.rm = TRUE),
    FineWoody = sum(ifelse(LITTER_ORGAN == "Branches", LITTER_Mean, 0), na.rm = TRUE),
    FineWoody_SD = sum(ifelse(LITTER_ORGAN == "Branches", `LITTER_Standard Deviation`, 0), na.rm = TRUE),
    CoarseWoody = sum(ifelse(LITTER_ORGAN == "Trunks", LITTER_Mean, 0), na.rm = TRUE),
    CoarseWoody_SD = sum(ifelse(LITTER_ORGAN == "Trunks", `LITTER_Standard Deviation`, 0), na.rm = TRUE)
  ) %>%
  filter(!is.na(Year))



plot_data <- aggregated_data %>%
  pivot_longer(
    cols = c(NonWoody, FineWoody, CoarseWoody),
    names_to = "LITTER_TYPE",
    values_to = "LITTER"
  ) %>%
  mutate(
    LITTER_SD = case_when(
      LITTER_TYPE == "FineWoody" ~ FineWoody_SD,
      LITTER_TYPE == "CoarseWoody" ~ CoarseWoody_SD,
      LITTER_TYPE == "NonWoody" ~ NonWoody_SD,
    )
  )

max_values <- filtered_data %>%
  group_by(Year) %>%
  summarise(MaxValue = max(LITTER_STATISTIC_NUMBER, na.rm = TRUE))

plot_data <- plot_data %>%
  left_join(max_values, by = "Year")
# Correctly setting the color override
plot_data <- plot_data %>%
  mutate(LITTER_TYPE_ORDER = factor(LITTER_TYPE, levels = c("CoarseWoody", "FineWoodyPool", "FineWoody", "NonWoody")),
         color_override = ifelse(LITTER_TYPE == "FineWoody" & Year == min(Year), "FineWoodyPool", LITTER_TYPE))

plot_data_coarse <- plot_data %>%
  filter(LITTER_TYPE == "CoarseWoody")
plot_data_fine_non <- plot_data %>%
  filter(LITTER_TYPE %in% c("FineWoody", "FineWoodyPool", "NonWoody"))


# CoarseWoody plot
plot_coarse <- ggplot(plot_data_coarse, aes(x = Year, y = LITTER, fill = color_override)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black") +
  geom_errorbar(aes(ymin = LITTER - LITTER_SD, ymax = LITTER + LITTER_SD),
                position = position_dodge(width = 0.9), width = 0.2, color = "black") +
  scale_fill_manual(values = c("CoarseWoody" = "#BEAED4"), name = "LITTER_FRACTION") + 
  labs(
    x = "Year",
    y = "Coarse woody litter mass (kg/m²)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    legend.position = "bottom"
  ) +
  scale_x_discrete(labels = function(x) {
    year_labels <- max_values$MaxValue[match(x, max_values$Year)]
    paste0(x, " \n(n = ", year_labels, ")")
  })

# FineWoody, FineWoodyPool, NonWoody plot
plot_fine_non <- ggplot(plot_data_fine_non, aes(x = Year, y = LITTER, fill = color_override)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black") +
  geom_errorbar(aes(ymin = LITTER - LITTER_SD, ymax = LITTER + LITTER_SD),
                position = position_dodge(width = 0.9), width = 0.2, color = "black") +
  scale_fill_manual(values = c("FineWoodyPool" = "#BF5B17", "FineWoody" = "#FDC086", "NonWoody" = "#7FC97F"),
                    breaks = c("FineWoodyPool", "FineWoody", "NonWoody"),
                    labels = c("FineWoodyPool", "FineWoody", "NonWoody")) +
  labs(
    x = "Year",
    y = "Fine- and non-woody litter flux (kg/m²/year)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    legend.position = "bottom",
    legend.title = element_blank()  
  ) +
  scale_x_discrete(labels = function(x) {
    year_labels <- max_values$MaxValue[match(x, max_values$Year)]
    paste0(x, " \n(n = ", year_labels, ")")
  })

# Combineer plots
combined_plot <- (plot_coarse + plot_fine_non) +
  plot_layout(ncol = 2) + 
  plot_annotation(
    title = "Mean (+/- SD) annual variation in dry litter mass at SE-Svb",  
    subtitle = "(n = number of CPs)",  
    theme = theme(
      plot.title = element_text(face = "bold", hjust = 0.5, size = 14),  
      plot.subtitle = element_text(hjust = 0.5, size = 10, margin = margin(b = 25))   
    )
  )

combined_plot

Extra comments about the output:

• Mass of coarse woody debris is calculated as descirbed in detail by Woodal & Monleon (2008). This type of litter is left in situ.

• Fine- and non-woody litter is collected from the field on a regular basis, hence this type of litter can be expressed as a flux (kg/m²/year).

• The first fine-woody litter can be seen as the pool of this type of litter.

• Non-woody litter is here presented as the summed litter fractions. Details on litter per campaign, fraction and species can be downloaded from the L2 archive.

filtered_data2 <- SESvb_LITTER %>% #### double check: delete _filtered
  filter(LITTER_APPROACH == "Non-Woody Litter and twigs with diameter <2cm collected in minimum 5 traps (minimum 0.5m-2) per plot (number of plots reported in LITTER_NUMBER). For more info visit www.icos-etc.eu") %>%
  filter(LITTER_STATISTIC == "Mean") %>%
  filter(!is.na(LITTER_DATE)) %>%
  mutate(Quarter = paste0(year(LITTER_DATE), "-Q", quarter(LITTER_DATE))) %>%
  select(Quarter, LITTER_ORGAN, LITTER) %>%
  group_by(Quarter, LITTER_ORGAN) %>%
  summarise(LITTER = sum(LITTER, na.rm = TRUE), .groups = 'drop')

ggplot(filtered_data2, aes(x = Quarter, y = LITTER, fill = LITTER_ORGAN)) +
  geom_bar(stat = "identity", color = "black") +
  scale_fill_brewer(palette = "Accent") +
  labs(x = "Quarter", y = "Summed dry litter mass (kg/m²)", fill = "Litter Type", title = "Stacked non-woody litter per quarter at SE-Svb") +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

If you have questions or suggestions about the above output and/or graphs, please contact icos-etc@unitus.it, arne.iserbyt@uantwerpen.be, maarten.opdebeeck@uantwerpen.be, bert.gielen@uantwerpen.be and trottacarlo@unitus.it