Summary Ancillary Data for SE-Nor Authors: ETC - Belgium: Arne Iserbyt, Emmy Jacobs, Maarten Op de Beeck and Bert Gielen ETC - France: Bruna Winck, Sebastien Lafont, Christelle Aluome, Nicolas Saby, Jean-Philippe Chenu, Nicolas Proix, Denis Loustau, Benjamin Loubet ETC - Italy: Carlo Trotta, Simone Sabbatini, Giacomo Labbri, Giacomo Nicolini, Adriana Mariotti, Eleonora Canfora, Alessio Ribeca, Dario Papale Latest version: 2025-10-29

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:

Data collection following the instruction documents as published on the ICOS-ETC website. Data submission occurs via BADM to ETC.
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.
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.
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.
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.
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)
  })
}

SENor_LAI <- safe_read_add_columns(path, "BIFTAB__GRP_LAI__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_LITTER_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_LITTER__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_BIOMASS_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_BIOMASS__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_DBH <- safe_read_add_columns(path, "BIFTAB__GRP_DBH__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_HEIGHT <- safe_read_add_columns(path, "BIFTAB__GRP_HEIGHTC__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_SPP_O_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_SPP_O__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_TREES_NUM_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_TREES_NUM__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")
SENor_BASAL_AREA_L2data <- safe_read_add_columns(path, "BIFTAB__GRP_BASAL_AREA__ICOSETC_SE-Nor_ANCILLARY__.*\\.csv")


### Select the most common species (>5% of species contribution)
SPP_Threshold <- SENor_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-Nor",
    !TREE_STATUS %in% c("Dead", "Removed", "Fallen")
  )


TREEs_Site_CPs <- TREEs_Site %>%
  filter(
    grepl("^CP", TREE_PLOT) ) %>%
  mutate(Year = ifelse(Year == "2017", "2018", Year)) #unhashtag if needed for this site

    TREEs_Site2 <- TREEs_Site %>%
  mutate(Year = ifelse(Year == "2017", "2018", Year)) #unhashtag if needed for this site
    
 
### Load coordinates of the plot locations
SENor_Spatial <- read_excel("I:/ETC_ancillary data/Labelled/Aggregated data/Aggregations_L2_Forests/DATABASE_exports/GRP_PLOT_AllSites.xlsx")
SENor_Spatial2 <- SENor_Spatial %>% 
  filter(Site == "SE-Nor")

transformed_data <- SENor_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-Nor",
    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 <- SENor_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-Nor ",
    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 <- SENor_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-Nor ",
    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"
  )

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

SENor_SPP_O_L2data_sorted <- SENor_SPP_O_L2data_filtered %>%
  arrange(Plot_Year, desc(SPP_O))

SENor_SPP_O_L2data_sorted <- SENor_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(SENor_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(SENor_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-Nor",
    # 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-Nor",
    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 = "_")
  )

SENor_Spatial_filtered <- SENor_Spatial2 %>%
  select(PLOT_ID, PLOT_LOCATION_LAT, PLOT_LOCATION_LONG)
Biomass_PLOT <- Biomass_PLOT %>%
  left_join(SENor_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 == "2018")

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-Nor",
    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 <- SENor_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-Nor",
    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 <- SENor_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-Nor",
    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.

SENor_BIOMASS_DeadTrees <- SENor_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 <- SENor_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-Nor",
    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 <- SENor_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))%>%
  filter(Year != 2020)

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-Nor",
    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.

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

#only data from CP - thus excluding SP-I DHPs from labelling
SENor_LAI <- SENor_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.")

SENor_LAI_mean <- SENor_LAI %>% filter(LAI_STATISTIC == "Mean")
SENor_LAI_sd <- SENor_LAI %>% filter(LAI_STATISTIC == "Standard Deviation")

combined_data <- SENor_LAI_mean %>%
  rename(Mean_LAI = LAI) %>%
  select(DATE2, Mean_LAI) %>%
  inner_join(SENor_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-Nor") +
  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

SE-Nor is a Class2 site, hence litter is not collected.

GRP_SOIL

Soil data collection is currently in progress for SE-Nor

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/"
SENor_VEG_CHEM_L2data <- read_latest_biftab("SE-Nor", data_dir)
SENor_VEG_CHEM_L2data <- SENor_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-Nor",
      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(SENor_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(SENor_VEG_CHEM_L2data, "VEG_CHEM_CONC_C", "C concentration ± SD (gC/kg)")
plot_N <- plot_conc(SENor_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-Nor",
    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

ETC did not yet receive phenocam images for SE-Nor.

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