---
title: "Piekmoment stuifmeelproductie bij berk en grassen"
date: 2026-03-25T10:00:00+01:00
bibliography: ../references.bib
link-citations: TRUE
thema:
  - Klimaat
keywords:
  - stuifmeel
  - pollen
  - gras
  - berk
lang: nl
tab: indicator
verantwoordelijke:
  - Lieve Vriens <lieve.vriens@inbo.be>
output: html_document
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)
library(git2rdata)
library(INLA)
library(plotly)
library(tidyverse)
library(effectclass)
library(INBOtheme)
library(htmlwidgets)
set.seed(20200615)
setWidgetIdSeed(20200615)
options(htmlwidgets.TOJSON_ARGS = list(pretty = TRUE))
```

```{r inlezen-gegevens}
read_vc("stuifmeel") %>%
  mutate(
    cjaar = jaar - min(jaar)
  ) -> stuifmeel
```

```{r analyse-berk}
model_berk_lineair <- inla(
  berk ~ cjaar,
  data = stuifmeel,
  control.compute = list(waic = TRUE, return.marginals.predictor = TRUE)
)
model_berk_niet_lineair <- inla(
  berk ~ f(
    cjaar,
    model = "rw2",
    hyper = list(theta = list(prior = "pc.prec", param = c(0.08, 0.01)))
  ),
  data = stuifmeel,
  control.compute = list(waic = TRUE, return.marginals.predictor = TRUE)
)

if (model_berk_lineair$waic$waic - model_berk_niet_lineair$waic$waic < 2) {
  warning(
    "Het lineaire model heeft een WAIC die minder dan 2 punten slechter is dan
    het niet-lineaire model."
  )
}

stuifmeel %>%
  bind_cols(
    sapply(
      model_berk_lineair$marginals.fitted.values,
      inla.qmarginal,
      p = c(0.5, 0.05, 0.2, 0.35, 0.65, 0.8, 0.95)
    ) %>%
      round(2) %>%
      t() %>%
      `colnames<-`(paste0(
        "berk_",
        c("mediaan", "lcl90", "lcl60", "lcl30", "ucl30", "ucl60", "ucl90")
      )) %>%
      as.data.frame()
  ) -> stuifmeel
piek_berk_10 <- format_ci(
  estimate = -10 * model_berk_lineair$summary.fixed["cjaar", "mean"],
  se = sqrt(10) * model_berk_lineair$summary.fixed["cjaar", "sd"],
  interval = 0.9
)
```

```{r analyse-grassen, include = FALSE}
model_grassen_lineair <- inla(
  grassen ~ cjaar,
  data = stuifmeel,
  control.compute = list(waic = TRUE, return.marginals.predictor = TRUE)
)
model_grassen_niet_lineair <- inla(
  grassen ~ f(
    cjaar,
    model = "rw2",
    hyper = list(theta = list(prior = "pc.prec", param = c(0.08, 0.01)))
  ),
  data = stuifmeel,
  control.compute = list(waic = TRUE, return.marginals.predictor = TRUE)
)
if (
  model_grassen_lineair$waic$waic - model_grassen_niet_lineair$waic$waic < 2
) {
  warning(
    "Het lineaire model heeft een WAIC die minder dan 2 punten slechter is dan
    het niet-lineaire model."
  )
}
stuifmeel %>%
  bind_cols(
    sapply(
      model_grassen_lineair$marginals.fitted.values,
      inla.qmarginal,
      p = c(0.5, 0.05, 0.2, 0.35, 0.65, 0.8, 0.95)
    ) %>%
      round(2) %>%
      t() %>%
      `colnames<-`(paste0(
        "grassen_",
        c("mediaan", "lcl90", "lcl60", "lcl30", "ucl30", "ucl60", "ucl90")
      )) %>%
      as.data.frame()
  ) -> stuifmeel
piek_grassen_10 <- format_ci(
  estimate = -10 * model_grassen_lineair$summary.fixed["cjaar", "mean"],
  se = sqrt(10) * model_grassen_lineair$summary.fixed["cjaar", "sd"],
  interval = 0.9
)
```

```{r bewaar-grafiekdata, results = "hide"}
stuifmeel %>%
  select(jaar, starts_with("berk"), starts_with("grassen")) %>%
  write_vc(
    "stuifmeel_grafiek", sorting = "jaar", optimize = FALSE, stage = TRUE
  )
```


De trend van het **piekmoment van de berk en grassen toont een significante vervroeging** over de jaren heen.
Het piekmoment van berkenstuifmeel vervroegt elke 10 jaar met `r piek_berk_10` dagen.
Het piekmoment van grassenstuifmeel vervroegt elke 10 jaar met `r piek_grassen_10` dagen.
Ook in 2025 was er een vervroeging in lijn met de piekmomenten van de laatste decennia.

```{r berk, fig.cap = "Figuur 1: evolutie van het piekmoment van berkenstuifmeel. De lijn toont de gemodelleerde trend met daarrond het 90% geloofwaardigheidsinterval (CI)."}
breaks <- as.Date(
  c(
    "2020-03-01", "2020-03-15", "2020-04-01", "2020-04-15", "2020-05-01",
    "2020-05-15"
  )
)
labels <- as.character(format(breaks, paste("%d", "%m", sep = "/")))
breaks <- as.integer(breaks - as.Date("2020-01-01")) + 1
p <- ggplot(stuifmeel, aes(x = jaar, y = berk_mediaan)) +
  geom_ribbon(aes(ymin = berk_lcl90, ymax = berk_ucl90,
                  fill = "90% CI"), alpha = 0.2) +
  # geom_ribbon(aes(ymin = berk_lcl60, ymax = berk_ucl60), alpha = 0.2) +
  # geom_ribbon(aes(ymin = berk_lcl30, ymax = berk_ucl30), alpha = 0.2) +
  geom_line(aes(color = "Trend berk")) +
  geom_point(aes(y = berk, color = "Piekmoment")) +
  scale_y_continuous("Dag van het jaar", breaks = breaks, labels = labels) +
  scale_color_manual(name = "Legende",
                     labels = c("Trend berk", "90% CI", "Piekmoment"),
                     values = c(inbo_palette(1), inbo_palette(1),
                                inbo_palette(1))) +
  guides(color = guide_legend("Legende"),
         fill = guide_legend("")) +
  ggtitle("") +
  theme(
    axis.title.x = element_blank()
  )

if (interactive() || opts_knit$get("rmarkdown.pandoc.to") == "html") {
  pp <- ggplotly(p, dynamicTicks = TRUE)
  for (i in seq_along(pp$x$data)) {
    pp$x$data[[i]]$name <- str_replace_all(pp$x$data[[i]]$name,
                                             c(",1" = "",
                                               "\\)" = "",
                                               "\\(" = ""))
  }

  pp <- pp %>% plotly::config(
      modeBarButtonsToRemove = list(
        "lasso2d", "select2d", "hoverClosestCartesian",
        "hoverCompareCartesian", "toggleSpikelines"
      ),
      displaylogo = FALSE
    )
  pp
} else {
  p
}
```

```{r grassen, fig.cap = "Figuur 2: evolutie van het piekmoment van grassenstuifmeel.De lijn toont de gemodelleerde trend met daarrond het 90% geloofwaardigheidsinterval (CI)."}
breaks <- as.Date(
  c(
    "2020-05-15", "2020-06-01", "2020-06-15", "2020-07-01", "2020-07-15"
  )
)
labels <- as.character(format(breaks, paste("%d", "%m", sep = "/")))
breaks <- as.integer(breaks - as.Date("2020-01-01")) + 1
p <- ggplot(stuifmeel, aes(x = jaar, y = grassen_mediaan)) +
  geom_ribbon(aes(ymin = grassen_lcl90, ymax = grassen_ucl90,
                  fill = "90% CI"), alpha = 0.2) +
  # geom_ribbon(aes(ymin = grassen_lcl60, ymax = grassen_ucl60), alpha = 0.2) +
  # geom_ribbon(aes(ymin = grassen_lcl30, ymax = grassen_ucl30), alpha = 0.2) +
  geom_line(aes(color = "Trend grassen")) +
  geom_point(aes(y = grassen, color = "Piekmoment")) +
  scale_y_continuous("Dag van het jaar", breaks = breaks, labels = labels) +
  scale_color_manual(name = "Legende",
                     labels = c("Trend grassen", "90% CI", "Piekmoment"),
                     values = c(inbo_palette(1), inbo_palette(1),
                                inbo_palette(1))) +
  guides(color = guide_legend("Legende"),
         fill = guide_legend("")) +
  theme(
    axis.title.x = element_blank()
  )

if (interactive() || opts_knit$get("rmarkdown.pandoc.to") == "html") {
pp <- ggplotly(p, dynamicTicks = TRUE)

  for (i in seq_along(pp$x$data)) {
    pp$x$data[[i]]$name <- str_replace_all(pp$x$data[[i]]$name,
                                             c(",1" = "",
                                               "\\)" = "",
                                               "\\(" = ""))
  }

  pp <- pp %>% plotly::config(
      modeBarButtonsToRemove = list(
        "lasso2d", "select2d", "hoverClosestCartesian",
        "hoverCompareCartesian", "toggleSpikelines"
      ),
      displaylogo = FALSE
    )
  pp
} else {
  p
}
```

## Definitie

Deze indicator toont de evolutie van het **piekmoment van de stuifmeelproductie** van berk en diverse grassoorten.

## Bespreking

Sinds 1974 meet het Belgisch aerobiologisch surveillance netwerk ([AirAllergy](https://airallergy.sciensano.be)) in Elsene de concentratie stuifmeel van berk en diverse grassen in de lucht.
De datum waarop de **stuifmeelconcentratie de hoogste waarde bereikt** (= piekmoment), kan sterk variëren van jaar tot jaar.
De trendlijn van het piekmoment toont zowel voor berk als voor grassen een **significante vervroeging** over de jaren heen.
In de periode 1975-1985 situeerde de piek voor berk zich rond 21 april, terwijl die in de periode 2010-2024 meer dan tien dagen vroeger ligt (Figuur 1).
Ook in 2025 zette deze vervroeging zich verder, met een piekdatum van 5 april.
Een gelijkaardig patroon wordt waargenomen bij grassen, waar de trend in de periode 2010-2024 rond 6 juni ligt, of bijna tien dagen vroeger dan in de periode 1975-1985 (Figuur 2).
Met een piekmoment op 31 mei volgt 2025 goed de trend van de voorbije jaren.
Er werd een **stijging** vastgesteld **van de dagelijkse pollenconcentraties** van diverse bomen en grassen.
Deze trend blijkt gecorreleerd te zijn met een graduele verandering van **meteorologische parameters** zoals temperatuur, straling, vochtigheid en regenval.
Hoewel deze taxonspecifieke gevoeligheid voor klimaatverandering mogelijk vermengd is met **andere factoren** zoals luchtvervuiling, biogeochemische aspecten en landgebruik [@bruffaerts.desmedt.ea2018ComparativeLongtermTrend] zijn de vervroeging en de stijging in pollenconcentratie in overeenstemming met de trend op andere locaties in het noordelijk halfrond [@ziska.makra.ea2019TemperaturerelatedChangesAirborne].

De verandering van het tijdstip van bloeien van graslandsoorten heeft ook **impact op de biodiversiteit**.
@sherry.zou.ea2006DivergenceReproductivePhenology rapporteerden dat soorten die bloeien voor de zomerpiek hun bloei- en vruchtperiode naar voor schuiven, terwijl de soorten die na de zomerpiek bloeien hun fenologische activiteiten naar achter schuiven.
Door deze verschuivingen ontstaat er een periode in het midden van het seizoen waar er geen of slechts weinig bloemen en vruchten voorhanden zijn.
Dat kan belangrijke implicaties hebben voor soorten, zoals bloembestuivers, die daarvan afhankelijk zijn [@burkle.alarcon2011FuturePlantPollinator].

Soortspecifieke reacties op klimaatwijziging kunnen eveneens leiden tot een **veranderde soortensamenstelling en structuur** van de graslandgemeenschap [@jones1997ClimateChangeImplications].

## Referenties