all use cases

Forest Health Assessment

Objective

Evaluate graphically to what extent the worst drought event in recorded meteorological history in Germany during the year of 2018 (Schuldt et al. 2020) can be recognized in meteorological and forest condition data from permanent ICP-Forest level II study sites in the state of Brandenburg. To limit the scope of work and exclude data unavailability periods, the study period is set to include ten years from 2014 to 2024 and only a subset of all available plots (1201, 1202, 1203, 1205, 1206, 1208) is assessed.

Project Setup

Importing libraries and setting a default theme for later plotting using ggplot2.

# Load necessary libraries
require(httr2)
require(dotenv)
require(data.table)
require(ggplot2)
require(RColorBrewer)
# Default theme for all plots
theme_set(
  theme_minimal() +
  theme(
    plot.title = element_text(
      color = "black", size = 12, hjust=0),
    plot.title.position = "plot",
    axis.title = element_text(size=8),
    axis.text = element_text(size=7),
    plot.caption.position = "plot",
    legend.title = element_text(size=8),
    legend.text = element_text(size=7),
    legend.position = "right",
  )
)
# Define custom colours
temp.red <- "#B22222"
prec.blue <- "#1874CD"

The home directory, desired database schemas and tables, as well as study period and survey plots to include are defined.

# Specify script directory
script_dir <- "/home/linus/Documents/FIT/semester_1/Forest_Inventory_and_Tree_Monitoring/llissel_code"
# Load environment variables from .env file
dotenv::load_dot_env(file = file.path(script_dir, ".env"))
# Supabase connection parameters
host <- "https://db.forstliche-umweltkontrolle.de/"
apikey <- Sys.getenv("ANON_KEY") # Get apikey from .env
# Check if apikey is available
if (apikey == "") {
    stop("API_KEY is not set in the .env file")
}
# Database configuration
data_schema <- "icp_download"
dictionaries_schema <- "icp_dictionaries"
# Tables to read
mm.mem <- "mm_mem"
cc.trc <- "cc_trc"
gr.ipm <- "gr_ipm"
tree.spec <- "d_tree_spec_fo"
# Time period
start_date = "2014-01-01"
end_date = "2024-12-31"
code_plots = "1201,1202,1203,1205,1206,1208"

Two functions for data fetching and descriptive statistics are created to limit recurring code later.

# Function to fetch data
fetch_data <- function(apikey, host, schema, table, date_col="date_observation", start=start_date, end=end_date, plot=code_plots, variable=NULL, select ){
  if(!is.null(variable)){
  url <- paste0(host, "rest/v1/", table,
                "?", date_col, "=gte.", start_date,
                "&", date_col, "=lte.", end_date,
                "&code_plot=in.(", plot, ")",
                "&code_variable=eq.", variable,
                paste0("&select=", select))
  } else{
   url <- paste0(host, "rest/v1/", table,
                "?", date_col, "=gte.", start_date,
                "&", date_col, "=lte.", end_date,
                "&code_plot=in.(", plot, ")",
                paste0("&select=", select)
                )
   }
  # Create and execute GET request
  response <- httr2::request(url) |>
    httr2::req_headers(
      "apikey" = apikey,
      "Authorization" = paste("Bearer", apikey),
      "Accept-Profile" = schema,
      "Accept" = "text/csv"
    ) |>
    httr2::req_perform()
  # Catch GET errors
  if (httr2::resp_status(response) != 200) {
    stop("Failed to fetch data: ", httr2::resp_status(response), " - ", httr2::resp_body_string(response))
  }
  data <- httr2::resp_body_string(response) |>
    read.csv(text = _)
  return(data)
  }

# Function to retrieve main statistical parameters
summary_short <- function(d){
  results <- c(
    sprintf("%.2f", min(d, na.rm=TRUE)),
    sprintf("%.2f", max(d, na.rm=TRUE)),
    sprintf("%.2f", mean(d, na.rm=TRUE)),
    sprintf("%.2f", median(d, na.rm=TRUE)),
    sprintf("%.2f", sd(d, na.rm=TRUE)),
    sprintf("%d", length(d))
  )
  labels<-c("MIN", "MAX", "MEAN", "MED", "STDEV", "COUNT")
  return(data.frame(row.names = labels, Param = results))
}

In the following, the previously defined data fetching function is used to retrieve meteorological and crown condition data for the chosen survey plots.

# Fetch meteorological data
for (variable in c("AT", "PR")){
  data <- as.data.table(fetch_data(
    a=apikey, h=host, sc=data_schema, table=mm.mem, v=variable,
    se="code_plot,survey_year,date_observation,daily_mean,daily_min,daily_max"))
  # Convert column type
  data[,code_plot:=as.factor(code_plot)]
  # Assign data to a variable named <variable>.data>
  assign(paste0(variable, ".data"), data)
}
rm(data)
# Drop unnecessary columns
PR.data[, `:=`(daily_min = NULL, daily_max = NULL)]
# Convert column types
AT.data[,date_observation:=as.Date(date_observation)]
PR.data[,date_observation:=as.Date(date_observation)]
# Aggregate to remove duplicate date-plot combinations
PR.data <- PR.data[, .(
  daily_mean = mean(daily_mean, na.rm = T)),
  by = .(code_plot, date_observation)]
# Remove PR values above threshold
PR.pre <- nrow(PR.data)
PR.data <- PR.data[daily_mean <= 130]
cat(sprintf("%i rows were removed from PR", PR.pre - nrow(PR.data)))
# Calculate main statistical parameters for the meteorological datasets
AT.stats <- summary_short(AT.data[,daily_mean])
PR.stats <- summary_short(PR.data[,daily_mean])

# Fetch crown condition data
CC.data <- as.data.table(fetch_data(
    a=apikey, h=host, sc=data_schema, table=cc.trc,
    date_col="date_survey", se="code_plot,survey_year,date_survey,code_tree_species,code_defoliation,code_social_class,code_tree_age"))
# Convert column types
CC.data[,code_plot:=as.factor(code_plot)]
CC.data[,date_survey:=as.Date(date_survey)]
# Fetch dictionary data to get tree species names
response <- httr2::request(paste0(host, "rest/v1/", tree.spec)) |>
  httr2::req_headers(
    "apikey" = apikey,
    "Authorization" = paste("Bearer", apikey),
    "Accept-Profile" = dictionaries_schema,
    "Accept" = "text/csv"
  ) |>
  httr2::req_perform()
spec.data <- httr2::resp_body_string(response) |>
  read.csv(text = _)
spec.data <- spec.data[, c("code", "description")]
# Clean CC.data
CC.data <- CC.data[code_defoliation != -1]
# Compute defoliation levels
CC.data[, defol_class := ifelse(code_defoliation <= 10, "No defoliation",
  ifelse(code_defoliation > 10 & code_defoliation <= 25, "Alert stage", "Severe defoliation"))]
CC.data[, defol_class := factor(defol_class, levels = c("Severe defoliation", "Alert stage", "No defoliation"))]
# Merge species name into CC.data
CC.data <- merge.data.table(CC.data, spec.data, by.x = "code_tree_species", by.y = "code")

To evaluate the validity of precipitation measurements, measurements from a German weather service (DWD) climate station roughly in the center of the assessed survey plots are compared to exclude potentially unrealistic values. The used data can be found at the opendata section of DWD.

# Load DWD precipitation data
dwd.cl.hist <- as.data.table(read.csv2(paste0(script_dir, "/data/tageswerte_KL_00433_hist/produkt_klima_tag_19480101_20241231_00433.txt"), dec = "."))
# Limit data to study period
dwd.cl.hist <- dwd.cl.hist[as.POSIXct(as.character(
  MESS_DATUM), format="%Y%m%d") >= start_date &
    as.POSIXct(as.character(
      MESS_DATUM), format="%Y%m%d") <= end_date]
# Retrieve maximum precipitation value
dwd.prec.max <- max(dwd.cl.hist[,RSK])

Daily and long term average values for precipitation and air temperature are calculated to use them later in the plotting stage.

# Calculation of daily average precipitation values over all survey plots
PR.daily.avg <- PR.data[, .(
  daily_avg = mean(daily_mean, na.rm = T)), by = date_observation]
# Calculation of monthly air temperature values by plot
AT.data[, `:=`(year = year(date_observation),month = month(date_observation))]
AT.monthly <- AT.data[, .(
  monthly_temp = mean(daily_mean, na.rm = TRUE), n_days = .N), by = .(
    code_plot, year, month)]
AT.monthly[, date := as.Date(paste(year, month, "15", sep = "-"))]
# Calculation of monthly precipitation values by plot
PR.data[, `:=`(year = year(date_observation),month = month(date_observation))]
PR.monthly <- PR.data[, .(
  monthly_prec = sum(daily_mean), n_days = .N), by = .(
    code_plot, year, month)]
PR.monthly[, date := as.Date(paste(year, month, "15", sep = "-"))]
# Calculation of long-term average values
AT.lt.avg <- AT.monthly[, .(
  lt_mean_temp = mean(monthly_temp, na.rm = TRUE)), by = .(code_plot, month)]
PR.lt.avg <- PR.monthly[, .(
  lt_mean_precip = mean(monthly_prec, na.rm = TRUE)), by = .(code_plot, month)]
# Calculation of anomaly values
AT.monthly[AT.lt.avg, lt_avg := i.lt_mean_temp, on = "month"]
AT.monthly[, temp_anomaly := monthly_temp - lt_avg]
PR.monthly[PR.lt.avg, lt_avg := i.lt_mean_precip, on = "month"]
PR.monthly[, precip_anomaly := monthly_prec - lt_avg]

@tbl-stats shows the main descriptive statistical parameters for the air temperature and precipitation during the study period, averaged over the evaluated survey plots. These are calculated using the previously defined function summary_short.

	Air Temperature	Precipitation
Minimum	-13.94	0.00
Maximum	29.54	176.10
Mean	9.86	1.70
Median	9.53	0.10
Standard Deviation	7.48	4.46
Count	43649	25569

: Main descriptive statistical parameters for the air temperature and precipitation data sets.

Visualization of Results

In this section first the source code and then its ggplot2 output are presented.

# Air temperature line plot
ggplot() +
  geom_line(AT.data, mapping = aes(
    x = date_observation, y = daily_mean),
    colour = temp.red, linewidth = 0.2) +
  theme(legend.position = "none",
    axis.title.x = element_blank()) +
  scale_x_date(
    date_breaks = "2 years",
    date_labels = "%Y") +
  ylab("Air Temperature [°C]")

# Precipitation average bar plot
ggplot() +
  geom_bar(PR.daily.avg, mapping = aes(
    x = date_observation, y = daily_avg),
    colour = prec.blue, linewidth = 0.3, stat = "identity") +
  theme(
    legend.position = "none",
    axis.title.x = element_blank()) +
  scale_y_continuous(
    breaks = seq(0, 50, by = 10),
    labels = seq(0, 50, by = 10)) +
  scale_x_date(
    date_breaks = "2 years",
    date_labels = "%Y") +
  ylab("Precipitation [mm]")

# Weather anomaly plot
ggplot() +
  geom_hline(
    yintercept = 0, linetype = "solid", color = "black", linewidth = 0.3) +
  geom_col(data = PR.monthly, aes(
    x = date, y = precip_anomaly/20, fill = "Precipitation")) +
  geom_line(data = AT.monthly, aes(
    x = date, y = temp_anomaly, colour = "Temperature"), linewidth = 0.5) +
  scale_y_continuous(
    name = "Temperature Anomaly [°C]",
    sec.axis = sec_axis(~ . * 20, name = "Precipitation Anomaly [mm]")) +
  scale_x_date(breaks = "2 years",
               date_labels = "%Y") +
  scale_fill_manual(values = c("Precipitation" = prec.blue)) +
  scale_colour_manual(values = c("Temperature" = temp.red)) +
  labs(fill = NULL, colour = NULL) +
  theme(
    axis.text.x = element_text(),
    axis.title.x = element_blank(),
    legend.position = "bottom")

# Crown condition stacked bar chart
cc.palette <- c("#D73027", "#FEE08B", "#1A9850")
ggplot() +
  geom_bar(CC.data, mapping = aes(
    y = survey_year, fill = defol_class),
    position = "fill", stat = "count") +
  scale_fill_discrete(palette = cc.palette) +
  guides(fill = guide_legend(reverse = T)) +
  labs(fill = "") +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "bottom") +
   scale_y_continuous(breaks = seq(min(CC.data[,survey_year]),
                                   max(CC.data[,survey_year]), by = 1))

### Temperature time series decomposition plot
# Aggregate values
AT.ts <- AT.data[, .(
  temp_mean = mean(daily_mean, na.rm = TRUE)
), by = .(year = year(date_observation), month = month(date_observation))]
setorder(AT.ts, year, month)
# Create a time series object
AT.ts <- ts(
  data = AT.ts$temp_mean,
  start = c(min(AT.ts$year), min(AT.ts$month)),
  frequency = 12
)
# Decompose and convert to data.table
temp.decomp <- decompose(AT.ts, type = "additive")
temp.decomp.dt <- data.table(
  date = as.numeric(time(temp.decomp$x)),
  observed = as.numeric(temp.decomp$x),
  trend = as.numeric(temp.decomp$trend),
  seasonal = as.numeric(temp.decomp$seasonal),
  random = as.numeric(temp.decomp$random)
)
# Convert to long format
temp.decomp.long <- rbindlist(list(
  temp.decomp.dt[, .(date, component = "Observed values", value = observed)],
  temp.decomp.dt[, .(date, component = "Trend component", value = trend)],
  temp.decomp.dt[, .(date, component = "Seasonal component", value = seasonal)],
  temp.decomp.dt[, .(date, component = "Random component", value = random)]
))
temp.decomp.long[, component := factor(component, levels = c(
  "Observed values", "Trend component", "Seasonal component", "Random component"
))]
# Create custom ggplot
ggplot(temp.decomp.long, aes(x = date, y = value)) +
  geom_line(data = temp.decomp.long[component != "Random component"],
            color = temp.red, linewidth = 0.5) +
  geom_col(data = temp.decomp.long[component == "Random component"],
           fill = temp.red, width = 0.05) +
  facet_wrap(~ component, ncol = 1, scales = "free_y") +
  ylab("Temperature [°C]") +
  scale_x_continuous(
    breaks = seq(2014, 2024, by = 2),
    labels = seq(2014, 2024, by = 2)) +
  theme(
    strip.text = element_text(size = 7, face = "bold"),
    axis.title.x = element_blank())

### Precipitation time series decomposition plot
# Aggregate values
# First calculate monthly sums per plot
PR.ts <- PR.data[, .(
  prec_sum = sum(daily_mean, na.rm = TRUE)
), by = .(code_plot, year = year(date_observation), month = month(date_observation))]
# Then average across all plots
PR.ts <- PR.ts[, .(
  prec_avg = mean(prec_sum, na.rm = TRUE)
), by = .(year, month)]
setorder(PR.ts, year, month)
# Create timeseries object
PR.ts <- ts(
  data = PR.ts$prec_avg,
  start = c(min(PR.ts$year), min(PR.ts$month)),
  frequency = 12
)
# Decompose and convert to data.table
prec.decomp <- decompose(PR.ts, type = "additive")
prec.decomp.dt <- data.table(
  date = as.numeric(time(prec.decomp$x)),
  observed = as.numeric(prec.decomp$x),
  trend = as.numeric(prec.decomp$trend),
  seasonal = as.numeric(prec.decomp$seasonal),
  random = as.numeric(prec.decomp$random)
)
# Reshape to long format
prec.decomp.long <- rbindlist(list(
  prec.decomp.dt[, .(date, component = "Observed values", value = observed)],
  prec.decomp.dt[, .(date, component = "Trend component", value = trend)],
  prec.decomp.dt[, .(date, component = "Seasonal component", value = seasonal)],
  prec.decomp.dt[, .(date, component = "Random component", value = random)]
))
prec.decomp.long[, component := factor(component, levels = c(
  "Observed values", "Trend component", "Seasonal component", "Random component"
))]
# Create timeseries plot
ggplot(prec.decomp.long, aes(x = date, y = value)) +
  geom_line(data = prec.decomp.long[component %in% c("Trend component", "Seasonal component")],
            color = prec.blue, linewidth = 0.5) +
  geom_col(data = prec.decomp.long[component %in% c("Observed values","Random component")],
           fill = prec.blue, width = 0.05) +
  facet_wrap(~ component, ncol = 1, scales = "free_y") +
  # ylab("Precipitation [mm]") +
  scale_x_continuous(
    breaks = seq(2014, 2024, by = 2),
    labels = seq(2014, 2024, by = 2)) +
  theme(
    strip.text = element_text(size = 7, face = "bold"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    )

Conclusion

The recorded and analysed data from the selected six ICP-Forests level II monitoring sites in Brandenburg exhibits the same weather patterns as observed across large parts of (central) Europe during 2018 and the following years, that led to the worst drought conditions in forestry history so far and unprecedented damage, increased by secondary factors such as wildfire and an explosion of tree damaging insect populations.

The full report on the matter can be found at [input report].