GRP_TREE GRP_AGB GRP_GAI GRP_LITTERPNT GRP_SOIL GRP_FLSM PHENOCAM

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. Leaf and soil samples are sent to ETC, which are analysed at INRAE labs. All ancillary data are always firmly quality checked (e.g. missing data, outliers, consistency, …), prior to entering the ICOS database.

  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 CP or SP-I plots and (2) average these averages across these 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 acite 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(openxlsx)
library(stringr)
library(ggplot2)
library(patchwork)
library(writexl)
# install.packages("knitr")
# install.packages("kableExtra")
library(knitr)
library(kableExtra)
library(akima)
library(viridis)
library(tinytex)
library(RColorBrewer)
library(lubridate) 
library(scales)
library(rmarkdown)
library(raster)
library(stars)
library(tmap)
library(tidyverse)
library(gridExtra)
library(grid)

#### 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"),
      Plot_Year = 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, Plot_Year, Biomass_TonsPerHa,
    BIOMASS_STATISTIC, BIOMASS_SPP, Species
  )

summarized_data <- transformed_data %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, BIOMASS_SPP, Species, Biomass_TonsPerHa) %>%
      rename(Biomass_TonsPerHa_SD = Biomass_TonsPerHa),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "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 = Plot_Year)) +
  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 = "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"
  )





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, Plot_Year, BasalArea_m2Ha,
    BASAL_AREA_STATISTIC, BASAL_AREA_SPP, Species
  )

summarized_data2 <- transformed_data2 %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, BASAL_AREA_SPP, Species, BasalArea_m2Ha) %>%
      rename(BasalArea_m2Ha_SD = BasalArea_m2Ha),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "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 = Plot_Year)) +
  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 = "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"
  )





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, Plot_Year, Density_TreesPerHa,
    TREES_NUM_STATISTIC, TREES_NUM_SPP, Species
  )

summarized_data3 <- transformed_data3 %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, TREES_NUM_SPP, Species, Density_TreesPerHa) %>%
      rename(Density_TreesPerHa_SD = Density_TreesPerHa),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "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 = Plot_Year)) +
  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 = "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"
  )





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(Plot_Year, desc(SPP_O))

SESvb_SPP_O_L2data_sorted <- SESvb_SPP_O_L2data_sorted %>%
  group_by(Plot_Year) %>%
  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 = Plot_Year, 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),
    Plot_Year = 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, Plot_Year, DBH,
    DBH_STATISTIC, DBH_SPP, Species
  )

summarized_data5 <- transformed_data5 %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, DBH_SPP, Species, DBH) %>%
      rename(DBH_SD = DBH),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "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 = Plot_Year)) +
  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, Plot_Year, HEIGHTC,
    HEIGHTC_STATISTIC, HEIGHTC_SPP, Species
  )

summarized_data4 <- transformed_data4 %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, HEIGHTC_SPP, Species, HEIGHTC) %>%
      rename(HEIGHTC_SD = HEIGHTC),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "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 = Plot_Year)) +
  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, Plot_Year, Biomass_TonsPerHa, BIOMASS_STATISTIC, BIOMASS_SPP, Species)

summarized_data <- SESvb_BIOMASS_DeadTrees %>%
  group_by(TREE_SPP, Year, Plot, Plot_Year, 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, Plot_Year, BIOMASS_SPP, Species, Biomass_TonsPerHa) %>%
      rename(Biomass_TonsPerHa_SD = Biomass_TonsPerHa),
    by = c("TREE_SPP", "Year", "Plot", "Plot_Year", "BIOMASS_SPP", "Species")
  )

# 🔥 **Make graph** 🔥
ggplot(final_data, aes(x = TREE_SPP, y = Biomass_TonsPerHa, fill = Plot_Year)) +
  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 = "Plot_Year"
  ) +
  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 in 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:



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)
  )

GRP_SOIL

Soil Organic Carbon (SOC) and Nitrogen (SON) stocks are assessed across a vertical soil profile from 0.0 to 1.0 meters depth. Sampling is conducted at standardized mineral soil depths: 0–5, 5–15, 15–30, 30–60, and 60–100 cm. Organic horizons are sampled separately above the mineral soil, typically between 0 and 15 cm, based on horizon thickness. The sampling design includes 20 discontinuous plots (SP-I), each containing 3 to 5 subplots (SP-II) for soil collection. This ensures spatial coverage across the target area. Soil sampling is repeated every 10 years at ICOS class1 and class2 stations.
All samples are analysed at INRAE labs, supervised by the French ETC team.
More detailed soil reports can be found on the French ICOS-ETC website.

MySite <- "SE-Svb"
dir.L2 <- file.path("I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/MOST_RECENT_BIFTAB_FOREST/SOIL", MySite, paste0(MySite, "_BIFTAB"))

soil_chem_file <- list.files(dir.L2, pattern = paste0(MySite, "_GRP_SOIL_CHEM.*\\.csv$"), full.names = TRUE)
soil_stock_file <- list.files(dir.L2, pattern = paste0(MySite, "_GRP_SOIL_STOCK.*\\.csv$"), full.names = TRUE)

GRP_SOIL_CHEM <- readr::read_csv(soil_chem_file)
GRP_SOIL_STOCK <- readr::read_csv(soil_stock_file)



GRP_SOIL_CHEM_BIF <- GRP_SOIL_CHEM %>%
  select(
    SITE_ID,
    SOIL_CHEM_PROFILE_MIN,
    SOIL_CHEM_PROFILE_MAX,
    SOIL_CHEM_C_ORG,
    SOIL_CHEM_C_ORG_STATISTIC,
    SOIL_CHEM_CN_RATIO,
    SOIL_CHEM_CN_RATIO_STATISTIC,
    SOIL_CHEM_BD,
    SOIL_CHEM_BD_STATISTIC,
    SOIL_CHEM_N_TOT,
    SOIL_CHEM_N_TOT_STATISTIC
  ) %>%
  # reshape into long format
  pivot_longer(
    cols = c(SOIL_CHEM_C_ORG, SOIL_CHEM_CN_RATIO, SOIL_CHEM_BD, SOIL_CHEM_N_TOT),
    names_to = "VARIABLE",
    values_to = "VALUE", 
    values_drop_na = TRUE
    
  ) %>%
  # match values with the right statistic column
  mutate(
    STATISTIC = case_when(
      VARIABLE == "SOIL_CHEM_C_ORG"    ~ SOIL_CHEM_C_ORG_STATISTIC,
      VARIABLE == "SOIL_CHEM_CN_RATIO" ~ SOIL_CHEM_CN_RATIO_STATISTIC,
      VARIABLE == "SOIL_CHEM_BD"       ~ SOIL_CHEM_BD_STATISTIC,
      VARIABLE == "SOIL_CHEM_N_TOT"    ~ SOIL_CHEM_N_TOT_STATISTIC
    )) %>%
  select(SITE_ID, VARIABLE, STATISTIC, VALUE, SOIL_CHEM_PROFILE_MIN, SOIL_CHEM_PROFILE_MAX) %>%
  pivot_wider(
    names_from = STATISTIC,
    values_from = VALUE
  ) %>%
  rename(Mean = Mean, SD = `Standard Deviation`) %>% 
  arrange (VARIABLE, SOIL_CHEM_PROFILE_MAX) %>% 
  ungroup()




GRP_SOIL_STOCK_BIF <- GRP_SOIL_STOCK %>%
  filter(SOIL_STOCK_COMMENT != "Aggregation into O and M layers") %>%
  select(
    SITE_ID,
    SOIL_STOCK_PROFILE_MIN,
    SOIL_STOCK_PROFILE_MAX,
    SOIL_STOCK_C_ORG, SOIL_STOCK_C_ORG_STATISTIC,
    SOIL_STOCK_N_TOT, SOIL_STOCK_N_TOT_STATISTIC
  ) %>%
  pivot_longer(
    cols = c(SOIL_STOCK_C_ORG, SOIL_STOCK_N_TOT),
    names_to = "VARIABLE",
    values_to = "VALUE"
  ) %>%
  mutate(
    STATISTIC = if_else(VARIABLE == "SOIL_STOCK_C_ORG",
                        SOIL_STOCK_C_ORG_STATISTIC,
                        SOIL_STOCK_N_TOT_STATISTIC)
  ) %>%
  select(SITE_ID, VARIABLE, STATISTIC, VALUE, SOIL_STOCK_PROFILE_MIN, SOIL_STOCK_PROFILE_MAX) %>%
  pivot_wider(
    names_from = STATISTIC,
    values_from = VALUE
  ) %>%
  rename(Mean = Mean, SD = `Standard Deviation`) %>%
  arrange(VARIABLE, SOIL_STOCK_PROFILE_MAX) %>%
  mutate(
    min_depth   = min(SOIL_STOCK_PROFILE_MIN),
    max_depth   = max(SOIL_STOCK_PROFILE_MAX),
    AGGREGATION = if_else(
      SOIL_STOCK_PROFILE_MIN == min_depth & SOIL_STOCK_PROFILE_MAX == max_depth,
      "Cumulative", "Layer"
    ),
    MinOrg = case_when(
      AGGREGATION == "Cumulative" ~ "Whole",
      SOIL_STOCK_PROFILE_MIN < 0  ~ "Organic",
      TRUE                        ~ "Mineral"
    )
  ) %>%
  ungroup()



# SOIL CHEMICAL VARIABLES




chem_labels <- c(
  "SOIL_CHEM_C_ORG"    = "SOC content (gC/kg)",
  "SOIL_CHEM_CN_RATIO" = "C:N ratio",
  "SOIL_CHEM_BD"       = "Bulk density (g/cm³)",
  "SOIL_CHEM_N_TOT"    = "Total N (gN/kg)"
)

depth_breaks <- sort(unique(GRP_SOIL_CHEM_BIF$SOIL_CHEM_PROFILE_MAX))

library(purrr)
GRP_SOIL_CHEM_BIF <- GRP_SOIL_CHEM_BIF %>%
  mutate(
    Mean = map_dbl(Mean, ~ .x[1]),
    SD   = map_dbl(SD, ~ .x[1])
  )


# Visualisatie
ggplot(
  GRP_SOIL_CHEM_BIF %>%
    filter(VARIABLE != "SOIL_CHEM_CN_RATIO") %>%
    mutate(
      Depth_mid = (SOIL_CHEM_PROFILE_MIN + SOIL_CHEM_PROFILE_MAX) / 2
    ),
  aes(x = SOIL_CHEM_PROFILE_MIN, y = Mean, fill = VARIABLE)
) +
  # Verticale referentielijnen
  geom_vline(xintercept = c(0, 5, 15, 30, 60, 100), linetype = "dashed", color = "gray", size = 0.5) +
  # Mean balken
  geom_rect(aes(xmin = SOIL_CHEM_PROFILE_MAX, xmax = SOIL_CHEM_PROFILE_MIN,
                ymin = Mean, ymax = Mean, color = VARIABLE), alpha = 1, size = 1) +
  # SD balken
  geom_rect(aes(xmin = SOIL_CHEM_PROFILE_MAX, xmax = SOIL_CHEM_PROFILE_MIN,
                ymin = Mean - SD, ymax = Mean + SD), alpha = 0.2) +
  # Diepte-as omkeren
  scale_x_reverse(breaks = c(0, 5, 15, 30, 60, 100)) +
  coord_flip() +
  xlab("Soil depth (cm)") +
  facet_wrap(~VARIABLE, scales = "free_x", labeller = as_labeller(chem_labels)) +
  theme_bw() +
  
  theme(
    strip.placement = "outside",
    panel.background = element_rect(fill = "white", colour = "white"),
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold", size = 12),  
    legend.position = "none",
    axis.text = element_text(face = "bold", size = 12),         
    axis.title = element_text(face = "bold", size = 14),      
    plot.title = element_text(hjust = 0.5, face = "bold", size = 16), 
  ) +
  labs(
    title = paste0("Soil characteristics (Mean ± SD) across soil depth at ", MySite,"\n")
    # subtitle = "Solid vertical lines represent reference depths\n Shaded ribbons indicate standard deviation",
    #  caption = "Data source: ICOS Ancillary Data"
  )

stock_labels <- c(
  "SOIL_STOCK_C_ORG" = "bold(SOC~stock~(g*C/m^{-2}))",
  "SOIL_STOCK_N_TOT" = "bold(Total~N~stock~(g*N/m^{-2}))"
)


depth_breaks <- sort(unique(GRP_SOIL_STOCK_BIF$SOIL_STOCK_PROFILE_MAX))

GRP_SOIL_STOCK_BIF %>%
  filter(AGGREGATION == "Layer") %>%
  mutate(Depth = as.numeric(SOIL_STOCK_PROFILE_MAX)) %>%
  arrange(VARIABLE, SITE_ID, Depth) %>%
  ggplot(aes(x = SOIL_STOCK_PROFILE_MIN, y = Mean, fill = VARIABLE, group = VARIABLE)) +
  geom_vline(xintercept = c(0, 5, 15, 30, 60, 100), size = 0.5, linetype = "dashed", color = "gray") +
  geom_rect(aes(xmin = SOIL_STOCK_PROFILE_MAX, xmax = SOIL_STOCK_PROFILE_MIN, 
                ymin = Mean, ymax = Mean, color = VARIABLE), alpha = 1, size = 1) +
  geom_rect(aes(xmin = SOIL_STOCK_PROFILE_MAX, xmax = SOIL_STOCK_PROFILE_MIN,
                ymin = Mean - SD, ymax = Mean + SD, fill = VARIABLE), alpha = 0.2) +
  scale_x_reverse(breaks = c(0, 5, 15, 30, 60, 100)) +
  coord_flip() +
  facet_wrap(
    ~VARIABLE, scales = "free_x",
    labeller = as_labeller(stock_labels, label_parsed)
  ) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    panel.background = element_rect(fill = "white", colour = "white"),
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold", size = 14),         
    axis.text = element_text(face = "bold", size = 12),          
    axis.title = element_text(face = "bold", size = 14),        
    plot.title = element_text(hjust = 0.5, face = "bold", size = 16),  
    plot.subtitle = element_text(hjust = 0.5, size = 12),       
    plot.caption = element_text(size = 10, face = "italic"),
    plot.margin = margin(t = 30, r = 10, b = 10, l = 10),        
    legend.position = "none"
  ) +
  labs(
    title = paste0("SOC and Total N Stocks (Mean ± SD) across soil depth at ", MySite),
    #    subtitle = "Solid vertical lines represent reference depths\nShaded ribbons indicate standard deviation",
    #   caption  = "Data source: ICOS Ancillary Data",
    x = "Soil depth (cm)",
    y = "Stock (g/m²)"
  )

###cummulative plot

GRP_SOIL_STOCK_BIF %>%
  filter(AGGREGATION == "Cumulative") %>%
  ggplot(aes(x = 0, y = Mean, fill = VARIABLE)) +
  
  geom_bar(stat = 'identity', position = position_dodge(0.9)) +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), 
                width = 0.2,
                position = position_dodge(0.9)) +
  
  geom_text(aes(y = Mean + SD + ((Mean + SD) * 0.05), 
                label = paste(round(Mean, 0))),
            position = position_dodge(0.9), 
            check_overlap = TRUE,
            fontface = "bold",
            color = "darkgreen",
            vjust = 0.5, 
            size = 5) +
  
  # Facetten met parsed labels
  facet_wrap(
    ~VARIABLE, scales = "free_y",
    labeller = as_labeller(stock_labels, label_parsed)
  ) +
  
  # Titels en labels
  labs(
    title = paste0("Cumulative SOC and Total N Stocks (Mean ± SD) at ", MySite),
    # subtitle = "Mean ± SD values shown with error bars",
    #  caption  = "Data source: ICOS Ancillary Data",
    x = NULL,
    y = "Stock (g/m²)"
  ) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    panel.background = element_rect(fill = "white", colour = "white"),
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold", size = 14),     
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_text(face = "bold", size = 12),
    axis.title.y = element_text(face = "bold", size = 14, margin = margin(r = 10)),
    axis.title.x = element_text(face = "bold", size = 14, margin = margin(t = 10)),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 16), 
    plot.subtitle = element_text(hjust = 0.5, size = 12),          
    plot.caption = element_text(size = 10, face = "italic"),
    plot.margin = margin(t = 30, r = 10, b = 30, l = 10),           
    legend.position = "none"
  )

########## spatial plot



dir.strata <- "I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/SOIL"

strata_file <- file.path(dir.strata, paste0(MySite, "__TA.strata.csv"))

TA_strata <- read_csv(strata_file) %>%
  rename(
    SOSM_PLOT_ID = sp1,
    Site = site
  )


sosm_file <- list.files(dir.L2, pattern = paste0("^", MySite, "__.*__GRP_SOSM_SPI\\.csv$"), full.names = TRUE)
GRP_SOSM <- read_csv(sosm_file)



SOC_stock_cum <- GRP_SOSM %>%
  group_by(SOSM_PLOT_ID) %>%
  reframe(SOSM_STOCK_C = round(sum(SOSM_STOCK_C, na.rm = TRUE), 0))


TA_strata<- left_join(TA_strata, SOC_stock_cum, by = c ("SOSM_PLOT_ID"))

r_soc <- rasterFromXYZ(TA_strata[, c("lon", "lat", "SOSM_STOCK_C")], crs=4326)
r_soc <- st_as_stars(r_soc)
names(r_soc) <- "SOSM_STOCK_C"

val_range_soc <- range(r_soc[[1]], na.rm = TRUE)
breaks_soc <- seq(min(val_range_soc), max(val_range_soc), length.out = 5)

suppressMessages(tmap_mode("plot"))

# https://github.com/cols4all/cols4all-R?tab=readme-ov-file

centroids = TA_strata %>% 
  mutate (idstratum = as.factor (idstratum)) %>% 
  group_by(idstratum) %>%
  summarize(lon = mean(lon), 
            lat = mean(lat),
            SOSM_STOCK_C = mean(SOSM_STOCK_C),
            .groups = 'drop') %>% 
  mutate(SOSM_STOCK_C = round(SOSM_STOCK_C, 0))

centroid.sf <- st_as_sf(centroids, coords = c("lon", "lat"), crs = 4326)

soc_map <- tm_shape(r_soc) +
  tm_raster(
    col = "SOSM_STOCK_C",
    col.scale = tm_scale_intervals(
      breaks = breaks_soc,
      values = "area9"
    ),
    col.legend = tm_legend(title = "SOC stock (g C/m²)")
  ) +
  tm_shape(centroid.sf) +
  tm_text(
    text = "SOSM_STOCK_C", 
    size = 0.6, 
    col = "black",
    bg.color = "white",
    bg.alpha = 0.7
  ) +
  tm_compass() + 
  tm_scalebar(position = c("right", "bottom")) +
  tm_layout(
    frame = TRUE,
    frame.col = "grey50",      
    frame.lwd = 1.1,
    legend.frame = TRUE,       
    legend.frame.col = "grey50", 
    legend.frame.lwd = 1,     
    legend.outside = TRUE,
    legend.outside.position = "bottom",
    asp = 1.4,
    outer.margins = c(0.01, 0.01, 0.01, 0.01),
    inner.margins = c(0.05, 0.01, 0.05, 0.01),
    title.size = 1.4,
    title.fontface = "bold",
    title.align = "center"
  ) +
  tm_credits(
    "SOC stock",
    position = c("CENTER", "TOP"),   
    size = 1.1,
    fontface = "bold"
  )


N_stock_cum <- GRP_SOSM %>% 
  group_by(SOSM_PLOT_ID) %>% 
  reframe (SOSM_STOCK_N = sum (SOSM_STOCK_N, na.rm = TRUE))


TA_strata<- left_join(TA_strata, N_stock_cum, by = c ("SOSM_PLOT_ID"))


r_n <- rasterFromXYZ(TA_strata[, c("lon", "lat", "SOSM_STOCK_N")], crs=4326)
r_n <- st_as_stars(r_n)
names(r_n) <- "SOSM_STOCK_N"

val_range_n <- range(r_n[[1]], na.rm = TRUE)
breaks_n <- seq(min(val_range_n), max(val_range_n), length.out = 5)

centroids = TA_strata %>% 
  mutate (idstratum = as.factor (idstratum)) %>% 
  group_by(idstratum) %>%
  summarize(lon = mean(lon), 
            lat = mean(lat),
            SOSM_STOCK_N = mean(SOSM_STOCK_N),
            .groups = 'drop') %>% 
  mutate(SOSM_STOCK_N = round(SOSM_STOCK_N, 0))

centroid.sf <- st_as_sf(centroids, coords = c("lon", "lat"), crs = 4326)

#tmap_mode("view")
tmap_mode("plot")
tmap_options(component.autoscale = FALSE)

# https://github.com/cols4all/cols4all-R?tab=readme-ov-file



n_map <- tm_shape(r_n) +
  tm_raster(
    col = "SOSM_STOCK_N",
    col.scale = tm_scale_intervals(breaks = breaks_n, values = "area9"),
    col.legend = tm_legend(title = "Total N stock (g N/m²)")
  ) +
  tm_shape(centroid.sf) +
  tm_text(text = "SOSM_STOCK_N", size = 0.6, col = "black", bg.color = "white", bg.alpha = 0.7) +
  tm_compass() + 
  tm_scalebar() +
  tm_layout(
    frame = TRUE,
    frame.col = "grey50",      
    frame.lwd = 1.1,
    legend.frame = TRUE,       
    legend.frame.col = "grey50", 
    legend.frame.lwd = 1,     
    legend.outside = TRUE,
    legend.outside.position = "bottom",
    asp = 1.4,
    outer.margins = c(0.01, 0.01, 0.01, 0.01),
    inner.margins = c(0.05, 0.01, 0.05, 0.01),
    title.size = 1.4,
    title.fontface = "bold",
    title.align = "center"
  ) +
  tm_credits(
    "Total N stock",
    position = c("CENTER", "TOP"),   
    size = 1.1,
    fontface = "bold"
  )


combined_map <- tmap_arrange(
  soc_map, n_map,
  ncol = 2,
  outer.margins = c(0.02, 0.02, 0.14, 0.02) # extra ruimte bovenaan (boven, rechts, onder, links)
)

print(combined_map)

grid.text("\n Spatial variation in cumulative soil stock at SE-Svb",
          y = 0.99, gp = gpar(fontsize = 16, fontface = "bold"))
grid.text("\n\n The values denote the stocks within each stratum",
          y = 0.955, gp = gpar(fontsize = 12))



GRP_FLSM

Leaf sampling focusses on both, leaf mass to area ratio (LMA) and on nutrient analysis (NA) of multiple macronutrients. Within each continuous measurement plot (CP), between 12 and 30 dominant or co-dominant healthy trees must be selected, ensuring a comparable number of trees per CP. LMA is measured once per year at both class1 and class2 sites. NA samples are also collected once per year, with class1 stations permitted to submit up to 90 samples annually, thus allowing multiple campaigns per year. Across all campaigns, the same trees are consistently resampled.
All samples are analysed at INRAE labs, supervised by the French ETC team. 

read_latest_biftab <- function(site_code, dir_path) {
  pattern <- paste0("BIFTAB__GRP_VEG_CHEM__ICOSETC_", site_code, "_ANCILLARY__.*\\.csv$")
  files <- list.files(dir_path, pattern = pattern, full.names = TRUE)
  
  if (length(files) == 0) {
    stop("❌  No file for site: ", site_code)
  }
  
  latest <- files[which.max(file.info(files)$mtime)]
  df <- read_csv(latest, show_col_types = FALSE, na = "")
  
  date_col <- grep("DATE$", names(df), value = TRUE)
  stat_col <- grep("STATISTIC$", names(df), value = TRUE)
  
  if (length(stat_col) > 0) {
    df <- df %>% filter(.data[[stat_col]] %in% c("Mean", "Standard Deviation"))
  }
  
  if (length(date_col) > 0) {
    df <- df %>% mutate(Year = substr(.data[[date_col]], 1, 4))
  } else {
    df <- df %>% mutate(Year = NA_character_)
  }
  
  df
}

data_dir <- "I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/MOST_RECENT_BIFTAB_FOREST/"
SESvb_VEG_CHEM_L2data <- read_latest_biftab("SE-Svb", data_dir)
SESvb_VEG_CHEM_L2data <- SESvb_VEG_CHEM_L2data %>%
  mutate(
    across(
      c(VEG_CHEM_LMA, VEG_CHEM_CONC_C, VEG_CHEM_CONC_N),
      ~ suppressWarnings(as.numeric(.))
    )
  )


df <- data.frame(
  Name = c('C/N','LMA','Ca','Cu','Fe','Mg','Mn','C','N','P','K','Zn','Dry Ratio',''),
  Unit = c('','kg/m2','g/kg','mg/kg','mg/kg','g/kg','mg/kg','g/kg','g/kg','g/kg','g/kg','mg/kg','%',''),
  Description = c('Carbon-to-nitrogen ratio','Leaf area per leaf dry mass','Total calcium (ICP-Radial) at 65 degC','Total copper (ICP-Radial) at 65 degC','Total iron (ICP-Radial) at 65 degC',
  'Total manganese (ICP-Radial) at 65 degC','Total magnesium (ICP-Radial) at 65 degC','Carbon (Dumas method) at 65 degC','Nitrogen (Dumas method) at 65 degC',
  'Total phosphorus (ICP-Radial) at 65 degC','Total potassium (ICP-Radial) at 65 degC','Total zinc (ICP-Radial) at 65 degC','Dry matter at 65 degC','')
)
kable(df, format = "markdown", caption = "Table of monitored variables ")
Table of monitored variables
Name Unit Description
C/N Carbon-to-nitrogen ratio
LMA kg/m2 Leaf area per leaf dry mass
Ca g/kg Total calcium (ICP-Radial) at 65 degC
Cu mg/kg Total copper (ICP-Radial) at 65 degC
Fe mg/kg Total iron (ICP-Radial) at 65 degC
Mg g/kg Total manganese (ICP-Radial) at 65 degC
Mn mg/kg Total magnesium (ICP-Radial) at 65 degC
C g/kg Carbon (Dumas method) at 65 degC
N g/kg Nitrogen (Dumas method) at 65 degC
P g/kg Total phosphorus (ICP-Radial) at 65 degC
K g/kg Total potassium (ICP-Radial) at 65 degC
Zn mg/kg Total zinc (ICP-Radial) at 65 degC
Dry Ratio % Dry matter at 65 degC 
plot_LMA <- function(df) {
  
  # Filter lege waarden eruit
  df <- df %>%
    filter(!is.na(VEG_CHEM_LMA))
  
  # Splits Mean en SD
  mean_data <- df %>%
    filter(VEG_CHEM_STATISTIC == "Mean") %>%
    rename(Mean_val = VEG_CHEM_LMA) %>%
    select(VEG_CHEM_DATE, VEG_CHEM_SPP, Mean_val)
  
  sd_data <- df %>%
    filter(VEG_CHEM_STATISTIC == "Standard Deviation") %>%
    rename(SD_val = VEG_CHEM_LMA) %>%
    select(VEG_CHEM_DATE, VEG_CHEM_SPP, SD_val)
  
  # Combineer en schaal
  combined_data <- mean_data %>%
    left_join(sd_data, by = c("VEG_CHEM_DATE", "VEG_CHEM_SPP")) %>%
    mutate(
      Mean_val = Mean_val * 10,
      SD_val   = SD_val * 10,
      DATE2    = as.Date(as.character(VEG_CHEM_DATE), format = "%Y%m%d")
    )
  
  # Tijdslijn instellen
  min_date <- min(combined_data$DATE2, na.rm = TRUE)
  max_date <- max(combined_data$DATE2, na.rm = TRUE)
  start_year <- as.numeric(format(min_date, "%Y"))
  start_date <- as.Date(paste0(start_year, "-01-01"))
  
  major_breaks <- seq(start_date, max_date, by = "1 year")
  minor_breaks <- seq(start_date, max_date, by = "6 months")
  
  # Plot
  ggplot(combined_data, aes(x = DATE2, color = VEG_CHEM_SPP)) +
    geom_line(aes(y = Mean_val, group = VEG_CHEM_SPP), linewidth = 0.8) +
    geom_point(aes(y = Mean_val), size = 2.5,
               position = position_jitter(width = 5, height = 0)) +
    geom_errorbar(aes(ymin = Mean_val - SD_val, ymax = Mean_val + SD_val),
                  width = 0.2, na.rm = TRUE) +
    labs(
      title = "Leaf Mass Area at SE-Svb",
      x = "Date",
      y = "LMA ± SD (kg/m²)",
      color = "Species"
    ) +
    scale_x_date(
      breaks = major_breaks,
      labels = date_format("%d-%m-%Y"),
      minor_breaks = minor_breaks
    ) +
    expand_limits(y = 0) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
      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.margin = margin(t = 30, b = 0, l = 0, r = 0),
      panel.grid.major.x = element_line(color = "grey85", linewidth = 0.5),
      panel.grid.minor.x = element_line(color = "grey85", linewidth = 0.3)
    )
}
  plot_LMA(SESvb_VEG_CHEM_L2data)

plot_conc <- function(df, value_col, y_label) {
  
  mean_data <- df %>%
    filter(VEG_CHEM_STATISTIC == "Mean") %>%
    rename(Mean_val = !!sym(value_col)) %>%
    select(VEG_CHEM_DATE, VEG_CHEM_SPP, Mean_val)
  
  sd_data <- df %>%
    filter(VEG_CHEM_STATISTIC == "Standard Deviation") %>%
    rename(SD_val = !!sym(value_col)) %>%
    select(VEG_CHEM_DATE, VEG_CHEM_SPP, SD_val)
  
  combined_data <- mean_data %>%
    inner_join(sd_data, by = c("VEG_CHEM_DATE", "VEG_CHEM_SPP")) %>%
    mutate(DATE2 = as.Date(as.character(VEG_CHEM_DATE), format = "%Y%m%d"))
  
  min_date <- min(combined_data$DATE2, na.rm = TRUE)
  max_date <- max(combined_data$DATE2, na.rm = TRUE)
  start_year <- as.numeric(format(min_date, "%Y"))
  start_date <- as.Date(paste0(start_year, "-01-01"))
  
  major_breaks <- seq(start_date, max_date, by = "1 year")
  minor_breaks <- seq(start_date, max_date, by = "6 months")
  
  ggplot(combined_data, aes(x = DATE2, color = VEG_CHEM_SPP)) +
    geom_line(aes(y = Mean_val, group = VEG_CHEM_SPP), size = 0.8) +
    geom_point(aes(y = Mean_val), size = 2.5,
               position = position_jitter(width = 5, height = 0)) +
    geom_errorbar(aes(ymin = Mean_val - SD_val, ymax = Mean_val + SD_val),
                  width = 0.2) +
    labs(
      x = "Date",
      y = y_label,
      color = "Species",
      title = NULL
    ) +
    scale_x_date(
      breaks = major_breaks,                      
      labels = scales::date_format("%d-%m-%Y"),
      minor_breaks = minor_breaks                 
    ) +
    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")
    ) +
    guides(color = guide_legend(title.theme = element_text(face = "bold")))
}


# Maak de plots
plot_C <- plot_conc(SESvb_VEG_CHEM_L2data, "VEG_CHEM_CONC_C", "C concentration ± SD (gC/kg)")
plot_N <- plot_conc(SESvb_VEG_CHEM_L2data, "VEG_CHEM_CONC_N", "N concentration ± SD (gN/kg)")


combined_plot <- (plot_C + plot_spacer() + plot_N) +
  plot_layout(ncol = 3, widths = c(1, 0.05, 1), guides = "collect") +
  plot_annotation(
    title = "Leaf nutrient analysis for C and N at SE-Svb",
    theme = theme(
      plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
      legend.position = "none",
      legend.margin = margin(t = 30, b = 0, l = 0, r = 0)
    )
  )
print(combined_plot)



PHENOCAM

The phenocam images are processed following the method developped and described in detail by:

library(dplyr)
 library(phenocamr)
 #get the correct path
 df <- read_phenocam("I:/Phenocam/SE-Svb/ROI/SE-Svb_EN_1008_3day.csv")
 
 df <- expand_phenocam(df)
 df <- detect_outliers(df)
 df <- smooth_ts(df)
 transition_dates <- transition_dates(
   df,
   percentile = 0.8 )
 phenology_dates <- phenophases(df, internal = TRUE)
 
 plot(as.Date(df$data$date),
      df$data$smooth_gcc_90,
      type = "l",
      xlab = "",
      ylab = "")
 
 #adjust title
 title(main = "Phenocam based phenology at SE-Svb", font.main = 2)
 points(as.Date(df$data$date),df$data$gcc_90)
 mtext("Date", side = 1, line = 2.5, font = 2)  
 mtext("Gcc", side = 2, line = 2.5, font = 2) 
 
 # rising "spring" greenup dates
 abline(v = phenology_dates$rising$transition_50,
        col = "green")
 # falling "autumn" senescence dates
 abline(v = phenology_dates$falling$transition_50,
        col = "brown")

 rising_dates <- phenology_dates$rising %>%
   filter(gcc_value == "gcc_90") %>%
   pull(transition_50)  
 falling_dates <- phenology_dates$falling %>%
   filter(gcc_value == "gcc_90") %>%
   pull(transition_50) 
 # Converteer to julian day
 rising_dates <- as.Date(rising_dates, origin = "1970-01-01")
 falling_dates <- as.Date(falling_dates, origin = "1970-01-01")
 all_dates <- c(rising_dates, falling_dates)
 all_labels <- format(all_dates, "%y-%m-%d")  
 text(x = all_dates,
      y = par("usr")[3] - 0.005,  
      labels = all_labels,
      srt = 45,            
      adj = 1,           
      xpd = TRUE,          
      cex = 0.7,          
      col = "blue")  

 phenology_dates <- phenophases(df, out_dir = "I:/Phenocam/SE-Svb/ROI/", internal=FALSE)
#Plot transition dates on x-axis

Notes for SE-Svb:



If you have questions or suggestions about the above output and/or graphs, please contact , , , and