PHENTHAUproc

PHENTHAUproc is an early warning system, modelling the life cycle of oak processionary moth (OPM, Thaumetopoea processionea) with temperature data. It was created by Halbig et al. 2024¹ at FVA – Forest Research Institute Baden-Wuerttemberg, Germany and Institute Baden-Wuerttemberg, Germany and BOKU - University of Natural Resources and Life Sciences, Vienna, Austria.

The aim is to support practitioners in risk assessment and planning of measures. It consists of temperature-sum based models to predict the phenological development of OPM stages, bud swelling and leaf unfolding of oak (Quercus robur), the starvation-dependent mortality of OPM L1 larvae as well as the application period of possible plant protection agents (ppa) and biocides.

This package implements the phenological models of PHENTHAUproc in R for local weather station data and spatial data, using the terra package created by Robert J. Hijmans.

The models were developed, adjusted and tested using temperature data and OPM observation data from Brandenburg and Baden-Wuerttemberg, Germany. Consider this in calculating OPM development in other regions.

Installation & setup

PHENTHAUproc is available on CRAN.

library(PHENTHAUproc)

# Packages used for demonstration:
library(terra)
library(tidyterra)
library(ggplot2)
library(grDevices)

Required data structure

PHENTHAUproc works with temperature-sum based models. The calculation requires daily mean, min and max air temperature in Celsius. For local prognosis, a data frame with the columns date, hour and tmean is required. For regional calculations, a named list of SpatRasters with equal time attribute is needed. The names of the list objects have to be tmean, tmin and tmax. The last available year will be used as year of prediction. Be aware that some models require temperature data from the previous year (Default: 1. Sept from previous year).

Weather station data and raster data for Germany can be found at Climate Data Center of DWD (German Meteorological Service), the basis of our sample record. In our regional example, we use a 4*4 pixel cutout centered at FVA from the spatial data of the DWD dataset “HYRAS”.

Spatial resolution: 5 km x 5 km
Projection: ETRS89 / LCC Europe (EPSG:3034)
Parameter: mean min and max air temperature at 2 m
time: from 2019-09-01 to 2020-09-30

For local data, we use hourly temperature data in Celsius from the weather station of Freiburg im Breisgau(Stations_id: 01443) from 2019 to 2022. (You can use ?load_test() to load example datasets.)

# local data for PHENTHAUproc 
freiburg <- load_test("hour")

head(freiburg)
#>       date hour tmean
#> 1 20190901    0  21.4
#> 2 20190901    1  21.9
#> 3 20190901    2  20.7
#> 4 20190901    3  20.8
#> 5 20190901    4  21.2
#> 6 20190901    5  20.3

# spatial data for PHENTHAUproc
fva <-  load_test("SpatRaster")
fva
#> $tmean
#> class       : SpatRaster 
#> dimensions  : 4, 4, 396  (nrow, ncol, nlyr)
#> resolution  : 5000, 5000  (x, y)
#> extent      : 3835000, 3855000, 2360000, 2380000  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89-extended / LCC Europe 
#> source      : fva_tmean.nc 
#> names       : fva_tmean_1, fva_tmean_2, fva_tmean_3, fva_tmean_4, fva_tmean_5, fva_tmean_6, ... 
#> time (days) : 2019-09-01 to 2020-09-30 
#> 
#> $tmax
#> class       : SpatRaster 
#> dimensions  : 4, 4, 396  (nrow, ncol, nlyr)
#> resolution  : 5000, 5000  (x, y)
#> extent      : 3835000, 3855000, 2360000, 2380000  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89-extended / LCC Europe 
#> source      : fva_tmax.nc 
#> names       : fva_tmax_1, fva_tmax_2, fva_tmax_3, fva_tmax_4, fva_tmax_5, fva_tmax_6, ... 
#> time (days) : 2019-09-01 to 2020-09-30 
#> 
#> $tmin
#> class       : SpatRaster 
#> dimensions  : 4, 4, 396  (nrow, ncol, nlyr)
#> resolution  : 5000, 5000  (x, y)
#> extent      : 3835000, 3855000, 2360000, 2380000  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89-extended / LCC Europe 
#> source      : fva_tmin.nc 
#> names       : fva_tmin_1, fva_tmin_2, fva_tmin_3, fva_tmin_4, fva_tmin_5, fva_tmin_6, ... 
#> time (days) : 2019-09-01 to 2020-09-30

PHENTHAUproc

Calculation

With the phenthau function, all model predictions are calculated. For a data frame input, it will return a data frame with the computed start dates of the phenological events (like bud swelling, different larval stages, etc.) or % for OPM L1 larval mortality. For a SpatRasterlist input, phenthau returns a SpatRasterlist. The output will always be limited from 1. February to 30. September. For more output information see: ?phenthau

# The parameter default method is "dailymeanminmax". If you want to calculate with hourly data, you have to change the parameter manually.
regional <- phenthau(fva)
local <- phenthau(freiburg, params = parameter("hour", year = 2020))

# we can create a parameter list
params <- parameter("dailymeanminmax", year = 2020)

# and change single parameter
params$budswelling$ldt <- 5 # change lower development threshold for budswelling from Default to 5

regional_manipulated <- phenthau(fva, params = params)

rm(params, regional_manipulated)

Results

Oak phenology

For oak phenology, bud swelling and leaf unfolding are calculated based on the work of Menzel 1997². Bud swelling determines the feeding start of hatched OPM L1 larvae. Leaf unfolding sets the start of the possible ppa/biocide usage. Leaf unfolding can be calculated for Quercus robur with different parametrisations related to the oak clones from which phenology observation data for model creation were obtained.

L1 - hatch

Hatch of OPM defines the beginning of the first larval stage L1. PHENTHAUproc includes three hatch models.

Custers 2003³
Meurisse et al. 2012⁴
Wagenhoff et al. 2014⁵

OPM - stages

Hatch and bud swelling are starting the feeding period of OPM and the calculation of the further development stages:

L2-L6 larval stages
Pp Pupa
Ad Adult

Ppa/biocide application period

The application period of appropriate ppa/biocides starts with leaf unfolding of oak and ends with the beginning of larval stage (L4) of OPM.

Mortality

The starvation-related mortality is calculated. It depends on the temperature sum during the period from OPM L1 hatching to oak bud swelling. Lacking coincidence of these two phenological events can lead to mortality of the young larvae (Wagenhoff et al. 2013⁶, Halbig et al. 2024).

Presentation of local PHENTHAUproc

One graph (currently only in English language) is available for local results and requires ggplot2.

plot_station_step(local)

Presentation of regional PHENTHAUproc

Stages

Stages is a SpatRaster with one layer per day and values from 0 to 8 for each development stage of OPM rom the egg stage to the adult stage (see phenthau_legend() for assignment). The plot_stages function is a wrapper around terra::plot to preset categories and colors and subset the required date (Default: max(date)).

stages <- regional$stages

plot_stages(stages,
            time = "2020-06-05",
            main = "OPM 5. June 2020",
            axes = F,
            box = T)

OPM larval stages 5. June 2020

Mortality

To assign colors and categories directly to the SpatRaster or to use it elsewhere, you can use the get_legend function.

# return mortality from list
mort <- regional$mortality

# get_legend returns a dataframe with ID, category and colors for a spatial PHENTHAUproc output
legend <- get_legend("mortality")

# set levels and colors for mortality
levels(mort) <- legend[,c("ID", "category")]
terra::coltab(mort) <- legend[,c("ID", "colors")]

# plot mortality FVA 2020
terra::plot(mort,
            plg = list(title = "%"),
            main = "Mortality 2020",
            all_levels = T,
            axes = F,
            box = T)

Starvation related Mortality 2020

Examples

Example: Phenology of Quercus robur at Freiburg (FVA) in 2020

Single phenological models can be calculated with the phenology function using different parametrisations (see parameter() for all preset parameters). Single parameter can directly be manipulated to compare phenological models.

# show possible models and parametrisation: parameter()
leafunfolding <- phenology(fva,
                           model = "leafunfolding",
                           parametrisation = "quercus_robur_clone256_type1")

Example: Hatch at Freiburg (FVA) in 2020

custers <- phenology(fva, model = "hatch", parametrisation = "custers")
meurisse <- phenology(fva, model = "hatch", parametrisation = "meurisse")
wagenhoff <- phenology(fva, model = "hatch", parametrisation = "wagenhoff")

Example: Implementation with leaflet

Implementation with leaflet

Halbig, P., Stelzer, A., Baier, P., Pennerstorfer, J., Delb, H., Schopf, A. 2023 PHENTHAUproc – An early warning and decision support system for hazard assessment and control of oak processionary moth (Thaumetopoea processionea) doi ↩︎
Menzel, A. 1997 Phänologie von Waldbäumen unter sich ändernden Klimabedingungen - Auswertung der Beobachtungen in den Internationalen Phänologischen Gärten und Möglichkeiten der Modellierung von Phänodaten. Forstliche Forschungsberichte München, 164, 1-147.↩︎
Custers, C.J.L. 2003 Climate change and trophic synchronisation - A casestudy of the Oak Processionary Caterpillar. Master’s Thesis, Wageningen University.↩︎
Meurisse, N., Hoch, G., Schopf, A., Battisti, A. and Grégoire, J.C. 2012 Low temperature tolerance and starvation ability of the oak processionary moth: implications in a context of increasing epidemics. Agricultural and Forest Entomology, 14, 239-250. doi ↩︎
Wagenhoff, E., Wagenhoff, A., Blum, R., Veit, H., Zapf, D. and Delb, H. 2014 Does the prediction of the time of egg hatch of Thaumetopoea processionea (Lepidoptera: Notodontidae) using a frost day/temperature sum model provide evidence of an increasing temporal mismatch between the time of egg hatch and that of budburst of Quercus robur due to recent global warming? European Journal of Entomology, 111, 207-215. doi ↩︎
Wagenhoff, E., Blum, R., Engel, K., Veit, H., Delb, H., 2013. Temporal synchrony of Thaumetopoea processionea egg hatch and Quercus robur budburst. J. Pest. Sci. 86, 193–202. doi ↩︎