Vignette 4 of 4 for the pep725 R package. The most current version is available on GitHub and CRAN.

Introduction

Phenological timing varies systematically across space due to environmental gradients. Understanding these spatial patterns is fundamental to phenological research because it allows us to:

Scale observations: Extrapolate from point measurements to landscapes
Predict impacts: Forecast how climate change will shift phenological timing
Design networks: Optimize monitoring station placement
Validate models: Compare observed gradients with model predictions

This vignette covers spatial analysis and visualization functions in the pep725 package, organized into four parts:

What You’ll Learn

Section	Key Concepts	Functions
Part 1: Gradients	How phenology changes with altitude & latitude	`pheno_gradient()`
Part 2: Synchrony	Spatial coherence of phenological events	`pheno_synchrony()`
Part 3: Combined	Integrating gradient and synchrony results	`pheno_gradient()`, `pheno_synchrony()`
Part 4: Mapping	Visualizing station networks and patterns	`pheno_leaflet()`, `pheno_map()`

Prerequisites

This vignette assumes you have completed the “Getting Started” vignette and are familiar with:

The basic structure of PEP725 data
BBCH phenological phase codes
Day of Year (DOY) as a timing metric

Setup

library(pep725)
library(data.table)
library(ggplot2)

# Download the synthetic dataset
pep <- pep_download()

# Overview of the data
cat("Stations:", length(unique(pep$s_id)), "\n")
#> Stations: 10927
cat("Altitude range:", min(pep$alt, na.rm = TRUE), "-",
    max(pep$alt, na.rm = TRUE), "m\n")
#> Altitude range: 0 - 1933 m
cat("Latitude range:", round(min(pep$lat, na.rm = TRUE), 2), "-",
    round(max(pep$lat, na.rm = TRUE), 2), "N\n")
#> Latitude range: 14.44 - 68.44 N

Part 1: Phenological Gradients

What are Phenological Gradients?

A phenological gradient describes systematic changes in the timing of biological events along an environmental axis. Such gradients reflect spatial variations in climatic conditions and are fundamental for understanding large-scale patterns in phenology. The two most important phenological gradients are associated with elevation and latitude:

Altitudinal Gradient

With increasing elevation, air temperature generally decreases. As a result, phenological events tend to occur progressively later at higher altitudes. The rate at which phenological timing shifts with elevation is commonly referred to as the phenological lapse rate.

                     Mountain Profile
                          /\
                         /  \  <- Higher elevation = LATER phenology
                        /    \
                       /      \
         ____________/        \____________
              ^                     ^
         Earlier                  Earlier
        (lowland)                (lowland)

Typical values: In temperate regions of Europe, phenological phases are delayed by 2-4 days per 100 meters of elevation gain. In practical terms: - A station at 500m elevation might observe flowering 10-20 days later than a lowland station. - High-elevation stations (2000m+) can exhibit delays of 40-80 days relative to lowland sites

Latitudinal Gradient

Moving northward in the Northern Hemisphere is generally associated with cooler climatic conditions and changes in seasonal radiation regimes. Consequently, phenological events tend to occur later at higher latitudes.

Typical values: Across Europe, phenology is delayed by 2-5 days per degree of latitude toward the north. For example: - The latitudinal difference between southern France (43°N) and northern Germany (54°N) is about 11 degrees - This corresponds to an expected difference of 22-55 days in phenology

Why Study Gradients?

Understanding phenological gradients is valuable for several reasons:

Scaling observations: Gradients allow phenological timing to be extrapolated to locations without direct observations, based on elevation or latitude
Validating models: Phenological and crop models should reproduce observed spatial gradients
Detecting change: Temporal changes in the strength or shape of phenological gradients may signal climate-driven alterations in biological responses
Understanding adaptation: Systematic deviations from expected gradients may reveal local adaptation, genetic variation or the influence of site-specific factors.

Altitudinal Gradient Analysis

Let’s calculate the elevational gradient for apple flowering:

# Calculate elevational gradient for apple flowering
grad_alt <- pheno_gradient(
  pep = pep,
  variable = "alt",
  species = "Malus domestica",
  phase_id = 60,
  method = "robust"
)

print(grad_alt)
#> Phenological Gradient Analysis
#> -------------------------------------------------- 
#> Gradient variable: alt
#> Regression method: robust
#> Species: Malus domestica
#> Phase ID: 60
#> -------------------------------------------------- 
#> 
#>    variable method    slope intercept  r_squared     rmse n_points
#>      <char> <char>    <num>     <num>      <num>    <num>    <int>
#> 1:      alt robust 1.211219  175.2758 0.03001954 26.46126     8307
#>                                                           interpretation
#>                                                                   <char>
#> 1: Phenology is 1.2 days later per 100m elevation increase (R-sq = 0.03)
#> 
#> Interpretation:
#>    Phenology is 1.2 days later per 100m elevation increase (R-sq = 0.03)

Understanding the Function Parameters

Parameter	Purpose	Your Choice
`pep`	Input data	Your phenological dataset
`variable`	Environmental gradient	`"alt"` for elevation, `"lat"` for latitude
`species`	Species to analyze	Must match exactly (e.g., “Malus domestica”)
`phase_id`	BBCH phase code	60 = flowering
`method`	Regression method	`"robust"` recommended for phenological data

Interpreting the Results

The output contains several important metrics:

Metric	What It Tells You	How to Interpret
`slope`	Days delay per 100m elevation	The phenological lapse rate
`intercept`	Predicted DOY at 0m altitude	Theoretical baseline (extrapolation)
`r_squared`	Proportion of variance explained	How well elevation predicts phenology
`n_points`	Number of stations used	Sample size (affects reliability)

Key interpretation points:

Positive slope: Later phenology at higher elevations (expected)
Negative slope: Earlier phenology at higher elevations (unusual - investigate!)
Slope magnitude: Compare to literature values (2-4 days/100m typical)
R-squared: Values > 0.3 indicate meaningful gradient; < 0.1 suggests weak control

Latitudinal Gradient Analysis

The same function works for latitude gradients:

# Calculate latitude gradient
grad_lat <- pheno_gradient(
  pep = pep,
  variable = "lat",
  species = "Malus domestica",
  phase_id = 60,
  method = "robust"
)

print(grad_lat)
#> Phenological Gradient Analysis
#> -------------------------------------------------- 
#> Gradient variable: lat
#> Regression method: robust
#> Species: Malus domestica
#> Phase ID: 60
#> -------------------------------------------------- 
#> 
#>    variable method     slope intercept    r_squared     rmse n_points
#>      <char> <char>     <num>     <num>        <num>    <num>    <int>
#> 1:      lat robust 0.1380235  171.1081 0.0003640754 26.49627     8403
#>                                                            interpretation
#>                                                                    <char>
#> 1: Phenology is 0.1 days later per degree latitude increase (R-sq = 0.00)
#> 
#> Interpretation:
#>    Phenology is 0.1 days later per degree latitude increase (R-sq = 0.00)

Comparing Regression Methods

The pheno_gradient() function offers three regression methods. Choosing the right method matters because phenological data often contains outliers.

# Robust regression (default) - handles outliers
grad_robust <- pheno_gradient(pep, variable = "alt",
                              species = "Malus domestica",
                              phase_id = 60, method = "robust")

# Ordinary least squares
grad_ols <- pheno_gradient(pep, variable = "alt",
                           species = "Malus domestica",
                           phase_id = 60, method = "ols")

cat("Robust slope:", round(grad_robust$summary$slope, 2), "days/100m\n")
#> Robust slope: 1.21 days/100m
cat("OLS slope:", round(grad_ols$summary$slope, 2), "days/100m\n")
#> OLS slope: 0.22 days/100m

Which Method Should You Use?

Method	Strengths	Weaknesses	When to Use
`robust`	Resistant to outliers and leverage points	Slightly lower efficiency under ideal conditions (clean data)	Recommended default choice for phenological data
`ols`	Efficient when assumptions are met	Highly sensitive to outliers	Use after careful data screening
`quantile`	Models conditional medians or other quantiles	Less precise estimates (wider confidence intervals)	When focusing on median responses

Recommendation: Start with method = "robust" for exploratory and applied analyses. Substantial differences between robust and OLS estimates often indicate influential observations and should prompt further inspection.

Comparing Species: Grapevine Gradient

Different species may show different gradient strengths. Let’s compare apple with grapevine:

# Calculate gradient for grapevine
vine_data <- pep[species == "Vitis vinifera"]

# Verify data exists before analysis
if (nrow(vine_data) > 0) {
  grad_vine <- pheno_gradient(
    pep = pep,
    variable = "alt",
    species = "Vitis vinifera",
    phase_id = 65,  # Full flowering
    method = "robust"
  )

  cat("Comparison of elevation gradients:\n")
  cat("Apple:     ", round(grad_robust$summary$slope, 2), "days/100m",
      "(R² =", round(grad_robust$summary$r_squared, 3), ")\n")
  cat("Grapevine: ", round(grad_vine$summary$slope, 2), "days/100m",
      "(R² =", round(grad_vine$summary$r_squared, 3), ")\n")
} else {
  cat("No grapevine data available for gradient analysis.\n")
}
#> Comparison of elevation gradients:
#> Apple:      1.21 days/100m (R² = 0.03 )
#> Grapevine:  0.05 days/100m (R² = 0 )

Gradients by Region

Phenological gradients are not uniform across space. Their magnitude and shape can vary between regions due to: - Different climate zones - Local adaptation of species - Data quality differences

Estimating gradients separately for different regions:

# Calculate gradient by country
grad_by_country <- pheno_gradient(
  pep = pep,
  variable = "alt",
  species = "Malus domestica",
  phase_id = 60,
  by = "country",
  method = "robust"
)

print(grad_by_country$summary)
#>              country      slope intercept    r_squared      rmse n_points
#>               <char>      <num>     <num>        <num>     <num>    <int>
#>  1:           France -0.5098732 125.16238 0.0008288219 77.910433      109
#>  2:          Austria  4.3428535 160.84490 0.0795285751 39.639291     1234
#>  3:           Norway         NA        NA           NA        NA        2
#>  4:    Germany-South  1.9631768 168.49517 0.1352181144 15.117422     1924
#>  5:          Czechia -8.8919479 242.31970 0.2560345827 21.314550       12
#>  6:          Hungary 65.0277647  75.32681 0.9348290289 45.149454        6
#>  7:            Spain -0.2565161 156.88915 0.0006775534 35.110126       36
#>  8:         Slovakia  2.4919661 144.08313 0.5354808286  8.420740       74
#>  9:    Germany-North  0.8295158 177.84818 0.0122156400 19.037512     4568
#> 10:      Switzerland  3.1563711  99.53822 0.6686215332 12.101512      168
#> 11:             <NA> 55.6580206 180.04261 0.2556121362 23.787673       36
#> 12:           Poland -3.9848006 185.69537 0.3800259994  8.128931       11
#> 13:       Montenegro  2.5621524 179.25280 0.5421907136  8.021958        5
#> 14:        Lithuania  1.5653666 166.35464 0.0046741701 11.311140       28
#> 15:      Netherlands         NA        NA           NA        NA        1
#> 16:       Luxembourg         NA        NA           NA        NA        1
#> 17:          Croatia  0.7422652 192.98568 0.6199968934 34.405313        5
#> 18:         Slovenia  3.0580002 111.01714 0.8373262029  2.822048       13
#> 19:           Latvia         NA        NA           NA        NA        1
#> 20: Bosnia and Herz.  0.4356856 152.91546 0.0110994625  7.600161        6
#> 21:    Liechtenstein         NA        NA           NA        NA        1
#> 22:           Sweden  4.4369767 135.09253 0.1836587725 27.457113       58
#> 23:            Italy  2.3636444  96.69996 0.5716057177  3.974141        7
#> 24:           Russia         NA        NA           NA        NA        2
#> 25:            Yemen         NA        NA           NA        NA        1
#>              country      slope intercept    r_squared      rmse n_points
#>               <char>      <num>     <num>        <num>     <num>    <int>
#>                                                              interpretation
#>                                                                      <char>
#>  1: Phenology is 0.5 days earlier per 100m elevation increase (R-sq = 0.00)
#>  2:   Phenology is 4.3 days later per 100m elevation increase (R-sq = 0.08)
#>  3:                                                                    <NA>
#>  4:   Phenology is 2.0 days later per 100m elevation increase (R-sq = 0.14)
#>  5: Phenology is 8.9 days earlier per 100m elevation increase (R-sq = 0.26)
#>  6:  Phenology is 65.0 days later per 100m elevation increase (R-sq = 0.93)
#>  7: Phenology is 0.3 days earlier per 100m elevation increase (R-sq = 0.00)
#>  8:   Phenology is 2.5 days later per 100m elevation increase (R-sq = 0.54)
#>  9:   Phenology is 0.8 days later per 100m elevation increase (R-sq = 0.01)
#> 10:   Phenology is 3.2 days later per 100m elevation increase (R-sq = 0.67)
#> 11:  Phenology is 55.7 days later per 100m elevation increase (R-sq = 0.26)
#> 12: Phenology is 4.0 days earlier per 100m elevation increase (R-sq = 0.38)
#> 13:   Phenology is 2.6 days later per 100m elevation increase (R-sq = 0.54)
#> 14:   Phenology is 1.6 days later per 100m elevation increase (R-sq = 0.00)
#> 15:                                                                    <NA>
#> 16:                                                                    <NA>
#> 17:   Phenology is 0.7 days later per 100m elevation increase (R-sq = 0.62)
#> 18:   Phenology is 3.1 days later per 100m elevation increase (R-sq = 0.84)
#> 19:                                                                    <NA>
#> 20:   Phenology is 0.4 days later per 100m elevation increase (R-sq = 0.01)
#> 21:                                                                    <NA>
#> 22:   Phenology is 4.4 days later per 100m elevation increase (R-sq = 0.18)
#> 23:   Phenology is 2.4 days later per 100m elevation increase (R-sq = 0.57)
#> 24:                                                                    <NA>
#> 25:                                                                    <NA>
#>                                                              interpretation
#>                                                                      <char>
#>     variable method
#>       <char> <char>
#>  1:      alt robust
#>  2:      alt robust
#>  3:      alt robust
#>  4:      alt robust
#>  5:      alt robust
#>  6:      alt robust
#>  7:      alt robust
#>  8:      alt robust
#>  9:      alt robust
#> 10:      alt robust
#> 11:      alt robust
#> 12:      alt robust
#> 13:      alt robust
#> 14:      alt robust
#> 15:      alt robust
#> 16:      alt robust
#> 17:      alt robust
#> 18:      alt robust
#> 19:      alt robust
#> 20:      alt robust
#> 21:      alt robust
#> 22:      alt robust
#> 23:      alt robust
#> 24:      alt robust
#> 25:      alt robust
#>     variable method
#>       <char> <char>

What to look for: - Consistency of slopes: Are slopes consistent across regions? Similar slope estimates across regions suggest a robust and generalizable gradient pattern. - Explained variance: Are R-squared values similar? (Comparable R-squared values indicate that the gradient explains phenological variability consistently across regions.) - Sample size adequacy: Do sample sizes support reliable estimates? (Reliable gradient estimates require sufficient data coverage; small sample sizes (e.g. n < 20) should be interpreted with caution.)

Visualizing Gradients

Visualization plays a key role in interpreting and communicating phenological gradients, helping to assess linearity, uncertainty, and regional differences at a glance:

# Plot the gradient relationship
if (!is.null(grad_alt$data) && nrow(grad_alt$data) > 2) {
  p <- plot(grad_alt)
  print(p)
}
#> `geom_smooth()` using formula = 'y ~ x'

Reading the plot: - Each point represents a station-level mean phenophase (averaged across years). - The fitted line shows the estimated relationship between phenophase and elevation. - Points deviating strongly from the line indicate locations where elevation alone explains phenology poorly. - Systematic clustering of residuals may point to regional influences or additional controlling factors.

Expected Values and Troubleshooting

Reference Values from Literature

The following table summarises phenological gradient values reported in the literature. These can serve as benchmarks when interpreting your own results:

Gradient	Typical Range	Region	Sources
Altitude	2–4 days/100m	Europe (temperate)	Pellerin et al. (2012); Vitasse et al. (2009); Ziello et al. (2009)
Altitude	1–3 days/100m	Europe (warm / low elevation)	Vitasse et al. (2009); Ziello et al. (2009)
Latitude	2–4 days/degree	Europe	Rötzer & Chmielewski (2001)
Latitude	3–5 days/degree	North America	Hopkins (1918); Richardson et al. (2019)

Notes:

Altitude gradients are reported as delay in days per 100 m of elevation gain. Values vary by species: deciduous broadleaves (e.g. oak, hazel) tend toward the upper end, while conifers (e.g. spruce) show weaker responses (Ziello et al. 2009).
Latitude gradients are reported as delay in days per degree northward. Hopkins’
1. bioclimatic law predicts approximately 4 days per degree of latitude, which remains a useful first approximation.
Gradient strength may change over time: Chen et al. (2018) showed that altitude gradients in Europe have been declining since the 1980s, indicating more uniform phenology across elevations under warming.

Key references:

Pellerin, M. et al. (2012). Spring tree phenology in the Alps. European Journal of Forest Research, 131, 1957–1965. doi:10.1007/s10342-012-0646-1
Vitasse, Y. et al. (2009). Leaf phenology sensitivity to temperature in European trees. Agricultural and Forest Meteorology, 149, 735–744. doi:10.1016/j.agrformet.2008.10.019
Ziello, C. et al. (2009). Influence of altitude on phenology of selected plant species in the Alpine region. Climate Research, 39, 227–234. doi:10.3354/cr00822
Rötzer, T. & Chmielewski, F.-M. (2001). Phenological maps of Europe. Climate Research, 18, 249–257. doi:10.3354/cr018249
Richardson, A.D. et al. (2019). Testing Hopkins’ Bioclimatic Law with PhenoCam data. Applications in Plant Sciences, 7, e01228. doi:10.1002/aps3.1228
Chen, L. et al. (2018). Spring phenology at different altitudes is becoming more uniform under global warming in Europe. Global Change Biology, 24, 3969–3975. doi:10.1111/gcb.14288

When Your Results Differ from Expected

Observation	Possible Causes	What to Do
Slope too steep (>5 days/100m)	Outliers, microclimate effects	Check for data errors, use robust method
Slope too shallow (<1 day/100m)	Limited elevation range, species adaptation	Expand geographic scope
Negative slope	Data errors, southern aspect effects	Investigate unusual observations
Very low R²	Other factors dominate (precipitation, soil)	Consider multiple regression

Part 2: Phenological Synchrony

What is Phenological Synchrony?

Phenological synchrony describes the degree to which phenological events occur at similar times across different locations within a region during a given period. It addresses the question: Do stations experience the same phenological phase at roughly the same time?

Visualizing Synchrony Concepts

High Synchrony (SD = 5 days)        Low Synchrony (SD = 25 days)

Station A: ●────────●               Station A: ●────────●
Station B:  ●────────●              Station B:          ●────────●
Station C: ●────────●               Station C:   ●────────●
Station D:  ●────────●              Station D:              ●────────●
           ▲                                    ▲
           All stations flower                  Stations flower at
           within ~10 day window                very different times

Why Study Synchrony?

Ecosystem coherence: High synchrony suggests strong regional climate control, whereas low synchrony indicates a greater influence of local factors.
Monitoring network design: When synchrony is high, fewer stations may be sufficient to characterize a region; low synchrony implies the need for denser observation networks.
Detection of change: Temporal changes in synchrony can indicate shifts in climate variability, increasing spatial heterogeneity, or emerging local adaptations.
Ecological consequences: Synchrony influences ecological interactions such as pollination, pest dynamics, and trophic mismatch.

Calculating Synchrony

The pheno_synchrony() function calculates several metrics to quantify spatial variability:

# Calculate synchrony by country and year
sync <- pheno_synchrony(
  pep = pep,
  species = "Malus domestica",
  phase_id = 60,
  by = c("country", "year"),
  min_stations = 3
)

print(sync)
#> Phenological Synchrony Analysis
#> -------------------------------------------------- 
#> Species: Malus domestica
#> Phase ID: 60
#> Grouping: country, year
#> Min stations required: 3
#> -------------------------------------------------- 
#> 
#> Overall Statistics:
#>   Total groups: 1175 (valid: 709)
#>   Mean SD across stations: 9.6 days
#>   Mean CV: 8.1%
#>   Mean stations per group: 241.5
#> 
#> Synchrony Data:
#>     country  year n_stations mean_doy sd_doy cv_pct range_doy iqr_doy min_doy
#>      <char> <int>      <int>    <num>  <num>  <num>     <int>   <num>   <int>
#>  1:    <NA>  1951          9    142.8   8.03   5.63        27     6.5     129
#>  2:    <NA>  1952          7    136.6   9.02   6.60        23    10.8     122
#>  3:    <NA>  1953          5    136.0   6.28   4.62        13    12.0     130
#>  4:    <NA>  1954          9    140.4   7.17   5.10        23     9.5     131
#>  5:    <NA>  1955         12    141.5   9.59   6.78        28    15.5     126
#>  6:    <NA>  1956          9    141.6  10.97   7.75        34     9.0     121
#>  7:    <NA>  1957         11    136.6   8.45   6.18        25    13.8     120
#>  8:    <NA>  1958         12    133.7   9.57   7.16        31    15.2     115
#>  9:    <NA>  1959         14    129.9  13.46  10.37        54    14.0     106
#> 10:    <NA>  1960         13    135.2   9.45   6.99        27    14.0     122
#>     max_doy
#>       <int>
#>  1:     156
#>  2:     145
#>  3:     143
#>  4:     154
#>  5:     154
#>  6:     155
#>  7:     145
#>  8:     146
#>  9:     160
#> 10:     149
#> 
#> ... and 1165 more rows
#> 
#> Trend Analysis (change in SD over time):
#>              country n_years year_min year_max   slope intercept r_squared
#>               <char>   <int>    <int>    <int>   <num>     <num>     <num>
#>  1:             <NA>      73     1951     2023 -0.0201     48.06     0.036
#>  2:          Austria      87     1926     2025  0.0097     -7.34     0.034
#>  3: Bosnia and Herz.      30     1971     2024 -0.0708    153.54     0.069
#>  4:          Croatia      52     1967     2020 -0.0287     66.69     0.010
#>  5:          Czechia      25     1951     1980  0.1452   -276.73     0.059
#>  6:    Germany-North      73     1951     2023 -0.0289     66.53     0.463
#>  7:    Germany-South      73     1951     2023 -0.0457    100.76     0.633
#>  8:          Hungary       6     1969     2000      NA        NA        NA
#>  9:            Italy      15     1989     2003 -0.3111    630.23     0.411
#> 10:        Lithuania      45     1960     2004  0.0263    -45.71     0.051
#> 11:       Montenegro      36     1952     2021  0.0198    -28.97     0.009
#> 12:           Poland      42     1954     2007  0.0411    -75.34     0.045
#> 13:         Slovakia      25     2000     2024  0.1096   -211.88     0.569
#> 14:         Slovenia      60     1961     2020  0.0317    -53.15     0.060
#> 15:            Spain      22     1987     2020 -0.2292    471.56     0.440
#> 16:           Sweden       3     2016     2019      NA        NA        NA
#> 17:      Switzerland      42     1959     2024 -0.0629    138.43     0.286
#>     p_value            direction significant
#>       <num>               <char>      <lgcl>
#>  1:  0.1758 increasing synchrony       FALSE
#>  2:  0.2056 decreasing synchrony       FALSE
#>  3:  0.0984 increasing synchrony       FALSE
#>  4:  0.4160 increasing synchrony       FALSE
#>  5:  0.4501 decreasing synchrony       FALSE
#>  6:  0.0000 increasing synchrony        TRUE
#>  7:  0.0000 increasing synchrony        TRUE
#>  8:      NA                 <NA>          NA
#>  9:  0.3575 increasing synchrony       FALSE
#> 10:  0.1792 decreasing synchrony       FALSE
#> 11:  0.6400 decreasing synchrony       FALSE
#> 12:  0.2425 decreasing synchrony       FALSE
#> 13:  0.0000 decreasing synchrony        TRUE
#> 14:  0.0343 decreasing synchrony        TRUE
#> 15:  0.0005 increasing synchrony        TRUE
#> 16:      NA                 <NA>          NA
#> 17:  0.0260 increasing synchrony        TRUE
#> 
#> Significant trends detected:
#>   Germany-North: increasing synchrony (p = 0.0000)
#>   Germany-South: increasing synchrony (p = 0.0000)
#>   Slovakia: decreasing synchrony (p = 0.0000)
#>   Slovenia: decreasing synchrony (p = 0.0343)
#>   Spain: increasing synchrony (p = 0.0005)
#>   Switzerland: increasing synchrony (p = 0.0260)

Understanding the Parameters

Parameter	Purpose	Typical Values
`pep`	Input dataset	Phenological dataset prepared with `pep725`
`species`	Species to analyze	Single species for clear interpretation
`phase_id`	Phenophase coded by the BBCH-scale (Meier, 2018)	e.g. 60 (flowering), 65 (full flowering)
`by`	Grouping variables	`c("country", "year")` for regional, year-specific analysis
`min_stations`	Minimum number of stations per group	3-5 as a lower bound; 10+ recommended for robust estimates
`compute_trend`	Estimate temporal trends in synchrony	TRUE (default) when assessing changes over time

Understanding Synchrony Metrics

Different synchrony metrics capture complementary aspects of spatial variability. Selecting an appropriate metric depends on the analysis goal and data characteristics:

Metric	Formula	Interpretation	When to Use
`sd_doy`	Standard deviation of day-of-year	Absolute variability in days	Comparing synchrony across phases or regions
`cv_pct`	(SD/mean) × 100	Relative variability	Comparing phases with different mean timing
`range_doy`	Maximum minus minimum doy	Total observed spread	Assessing extreme differences
`iqr_doy`	75th - 25th percentile	Robust measure of spread	Preferred when outliers are present

Interpretation guide:

The following table provides a rule-of-thumb interpretation of synchrony based on the standard deviation of phenological timing across stations. Thresholds are indicative and should be interpreted in relation to species, phenophase, and spatial extent:

SD Value	Interpretation
< 5 days	Very high synchrony; stations behave very similarly
5-10 days	Moderate synchrony; regional pattern dominates
10-20 days	Low synchrony; local variability is substantial
> 20 days	Very low synchrony; local factors dominate

summary(sync)
#> Phenological Synchrony Summary
#> ================================================== 
#> 
#> Species: Malus domestica
#> Phase ID: 60
#> 
#> Groups analyzed: 1175
#> Groups with sufficient stations: 709 (60.3%)
#> 
#> Synchrony Statistics (across all valid groups):
#>   Mean SD: 9.6 days
#>   Median SD: 9.4 days
#>   Mean CV: 8.1%
#> 
#> SD Distribution:
#>   Min: 1.7 days
#>   25th percentile: 7.8 days
#>   Median: 9.4 days
#>   75th percentile: 11.5 days
#>   Max: 26.9 days
#> 
#> Trend Summary:
#>   Regions with significant trends: 6 of 17
#>   Increasing synchrony: 4
#>   Decreasing synchrony: 2

Temporal Trends in Synchrony

An important question in phenological analysis is whether spatial synchrony changes over time.

Decreasing synchrony: means growing spatial heterogeneity, for example due to more variable weather or diverging local conditions.
Increasing synchrony: suggests stronger regional coherence, potentially driven by a common climate signal overriding local factors

# Check if trend analysis was performed
if (!is.null(sync$trend) && nrow(sync$trend) > 0) {
  cat("Trend Analysis Results:\n")
  print(sync$trend)

  # Interpret significant trends
  sig_trends <- sync$trend[!is.na(significant) & significant == TRUE]
  if (nrow(sig_trends) > 0) {
    cat("\nSignificant trends detected:\n")
    print(sig_trends[, .(country, slope, direction, p_value)])
  }
}
#> Trend Analysis Results:
#> Index: <significant>
#>              country n_years year_min year_max   slope intercept r_squared
#>               <char>   <int>    <int>    <int>   <num>     <num>     <num>
#>  1:             <NA>      73     1951     2023 -0.0201     48.06     0.036
#>  2:          Austria      87     1926     2025  0.0097     -7.34     0.034
#>  3: Bosnia and Herz.      30     1971     2024 -0.0708    153.54     0.069
#>  4:          Croatia      52     1967     2020 -0.0287     66.69     0.010
#>  5:          Czechia      25     1951     1980  0.1452   -276.73     0.059
#>  6:    Germany-North      73     1951     2023 -0.0289     66.53     0.463
#>  7:    Germany-South      73     1951     2023 -0.0457    100.76     0.633
#>  8:          Hungary       6     1969     2000      NA        NA        NA
#>  9:            Italy      15     1989     2003 -0.3111    630.23     0.411
#> 10:        Lithuania      45     1960     2004  0.0263    -45.71     0.051
#> 11:       Montenegro      36     1952     2021  0.0198    -28.97     0.009
#> 12:           Poland      42     1954     2007  0.0411    -75.34     0.045
#> 13:         Slovakia      25     2000     2024  0.1096   -211.88     0.569
#> 14:         Slovenia      60     1961     2020  0.0317    -53.15     0.060
#> 15:            Spain      22     1987     2020 -0.2292    471.56     0.440
#> 16:           Sweden       3     2016     2019      NA        NA        NA
#> 17:      Switzerland      42     1959     2024 -0.0629    138.43     0.286
#>     p_value            direction significant
#>       <num>               <char>      <lgcl>
#>  1:  0.1758 increasing synchrony       FALSE
#>  2:  0.2056 decreasing synchrony       FALSE
#>  3:  0.0984 increasing synchrony       FALSE
#>  4:  0.4160 increasing synchrony       FALSE
#>  5:  0.4501 decreasing synchrony       FALSE
#>  6:  0.0000 increasing synchrony        TRUE
#>  7:  0.0000 increasing synchrony        TRUE
#>  8:      NA                 <NA>          NA
#>  9:  0.3575 increasing synchrony       FALSE
#> 10:  0.1792 decreasing synchrony       FALSE
#> 11:  0.6400 decreasing synchrony       FALSE
#> 12:  0.2425 decreasing synchrony       FALSE
#> 13:  0.0000 decreasing synchrony        TRUE
#> 14:  0.0343 decreasing synchrony        TRUE
#> 15:  0.0005 increasing synchrony        TRUE
#> 16:      NA                 <NA>          NA
#> 17:  0.0260 increasing synchrony        TRUE
#> 
#> Significant trends detected:
#>          country   slope            direction p_value
#>           <char>   <num>               <char>   <num>
#> 1: Germany-North -0.0289 increasing synchrony  0.0000
#> 2: Germany-South -0.0457 increasing synchrony  0.0000
#> 3:      Slovakia  0.1096 decreasing synchrony  0.0000
#> 4:      Slovenia  0.0317 decreasing synchrony  0.0343
#> 5:         Spain -0.2292 increasing synchrony  0.0005
#> 6:   Switzerland -0.0629 increasing synchrony  0.0260

Interpreting Trend Results

Trend Direction	Interpretation	Ecological meaning
Decreasing SD	Increasing synchrony	Phenological timing becoming more spatially coherent
Increasing SD	Decreasing synchrony	Growing spatial variability among stations
No trend	Stable synchrony	Persistent spatial patterns over time

Cautions with trend interpretation: - Changes in the observation network (e.g. stations added or removed) can introduce artificial trends. - Reliable trend detection typically requires at least 15–20 years of data. - Statistical significance should be interpreted alongside ecological relevance and effect size.

Visualizing Synchrony Over Time

# Plot synchrony time series
if (nrow(sync$data[!is.na(sd_doy)]) > 5) {
  p <- plot(sync)
  print(p)
}
#> `geom_smooth()` using formula = 'y ~ x'

Reading the plot: - Y-axis shows the synchrony metric (typically the SD of timing) - The fitted trend line indicates long-term change - Scatter around trend reflects year-to-year variability in synchrony

Synchrony Without Trend Analysis

For exploratory analyses or descriptive summaries, synchrony can be calculated without estimating temporal trends:

sync_simple <- pheno_synchrony(
  pep = pep,
  species = "Malus domestica",
  phase_id = 60,
  by = c("country", "year"),
  min_stations = 3,
  compute_trend = FALSE
)

# Just the synchrony data
head(sync_simple$data)
#>    country  year n_stations mean_doy sd_doy cv_pct range_doy iqr_doy min_doy
#>     <char> <int>      <int>    <num>  <num>  <num>     <int>   <num>   <int>
#> 1:    <NA>  1951          9    142.8   8.03   5.63        27     6.5     129
#> 2:    <NA>  1952          7    136.6   9.02   6.60        23    10.8     122
#> 3:    <NA>  1953          5    136.0   6.28   4.62        13    12.0     130
#> 4:    <NA>  1954          9    140.4   7.17   5.10        23     9.5     131
#> 5:    <NA>  1955         12    141.5   9.59   6.78        28    15.5     126
#> 6:    <NA>  1956          9    141.6  10.97   7.75        34     9.0     121
#>    max_doy
#>      <int>
#> 1:     156
#> 2:     145
#> 3:     143
#> 4:     154
#> 5:     154
#> 6:     155

This approach is useful for rapid assessment or when time series length is insufficient for robust trend estimation.

Comparing Synchrony: Grapevine vs. Apple

Different species may show different synchrony levels due to their biology:

# Calculate synchrony for grapevine
vine_check <- pep[species == "Vitis vinifera" & phase_id == 65]

if (nrow(vine_check) > 0) {
  sync_vine <- pheno_synchrony(
    pep = pep,
    species = "Vitis vinifera",
    phase_id = 65,
    by = c("country", "year"),
    min_stations = 3,
    compute_trend = FALSE
  )

  cat("Synchrony comparison (mean SD across years):\n")
  cat("Apple (flowering):    ", round(sync$overall$mean_sd_doy, 1), "days\n")
  cat("Grapevine (flowering):", round(sync_vine$overall$mean_sd_doy, 1), "days\n")
} else {
  cat("Insufficient grapevine data for synchrony comparison.\n")
}
#> Synchrony comparison (mean SD across years):
#> Apple (flowering):     9.6 days
#> Grapevine (flowering): 11.5 days

Part 3: Combining Gradient and Synchrony Analysis

A comprehensive spatial phenological analysis typically combines gradient and synchrony approaches. While each addresses a different aspect of spatial structure, together they provide a more complete picture of how phenology varies across landscapes.

Analysis	Question Answered
Gradient	How does MEAN phenological timing change across space?
Synchrony	How VARIABLE is phenology within regions?

Gradients describe systematic spatial shifts (e.g. with elevation or latitude), whereas synchrony captures the degree of spatial coherence around these shifts.

# 1. Understand elevation gradient
gradient <- pheno_gradient(
  pep, variable = "alt",
  species = "Malus domestica",
  phase_id = 60
)

cat("Elevation Gradient Analysis:\n")
#> Elevation Gradient Analysis:
cat("  Slope:", round(gradient$summary$slope, 2), "days/100m\n")
#>   Slope: 1.21 days/100m
cat("  R-squared:", round(gradient$summary$r_squared, 3), "\n\n")
#>   R-squared: 0.03

# 2. Assess spatial synchrony
synchrony <- pheno_synchrony(
  pep,
  species = "Malus domestica",
  phase_id = 60,
  min_stations = 3
)

cat("Synchrony Analysis:\n")
#> Synchrony Analysis:
cat("  Mean SD across stations:", round(synchrony$overall$mean_sd_doy, 1), "days\n")
#>   Mean SD across stations: 9.6 days
cat("  Mean CV:", round(synchrony$overall$mean_cv_pct, 1), "%\n")
#>   Mean CV: 8.1 %

The two analyses answer complementary questions: the gradient quantifies the direction and strength of spatial change, while synchrony indicates how tightly phenological timing aligns around that spatial pattern.

Interpreting Combined Results

Interpreting gradients and synchrony together helps to identify the dominant controls on phenology within a region:

Gradient	Synchrony	Interpretation	Example
Strong	High	Clear environmental control with uniform response	Continental climate regions
Strong	Low	Environmental control with strong local variation	Mountain regions with microclimates
Weak	High	Broad regional climate signal dominates	Coastal or maritime regions
Weak	Low	Phenology shaped by complex local factors	Highly fragmented landscapes

This combined perspective is particularly useful for comparing regions, diagnosing model performance, and identifying where local processes may override large-scale climatic drivers.

Part 4: Mapping Phenological Patterns

Spatial visualization is a critical component of phenological analysis. Before formal modeling, maps help to assess station coverage, identify spatial gaps, and explore geographic patterns in phenological timing. The pep725 package provides two mapping functions:

Function	Type	Best suited for
`pheno_leaflet()`	Interactive	Data exploration and station selection
`pheno_map()`	Static	Analysis summaries and publication-quality figures

Interactive Maps with pheno_leaflet()

The pheno_leaflet() function launches an interactive Shiny-based map for exploring phenological station networks. It is particularly useful during the early stages of analysis for:

exploring spatial coverage and station density,
interactively selecting subsets of stations,
understanding geographic context before aggregation or modeling.

# Launch interactive map (opens in viewer)
# Draw rectangles or polygons to select stations
selected_stations <- pheno_leaflet(pep, label_col = "species")

# The function returns a data.frame of selected stations
print(selected_stations)

#> *[Screenshot: Interactive map showing PEP725 stations across Europe with drawing toolbar for rectangle and polygon selection]*

Features of pheno_leaflet()

The interactive map supports multiple selection and inspection tools:

Feature	Purpose
Click selection	Select individual stations
Rectangle selection	Select rectangular geographic areas using the rectangle tool from the toolbar
Polygon selection	Define custom shapes to select irregular regions
Hover labels	Move mouse over points to inspect station metadata
Done button	Click “Done” to finalize and return selected stations

The function returns a data.frame containing all selected stations with their coordinates and metadata, which can then be directly reused in subsequent analysis.

Best Practices for Interactive Mapping

Tip: For large datasets, filtering before mapping substantially improves performance.

# Filter to a specific species for faster loading (recommended)
apple <- pep[species == "Malus domestica"]
selected <- pheno_leaflet(apple, label_col = "s_id")

# Use selected stations for focused analysis
apple_subset <- apple[s_id %in% selected$s_id]

Why filter first? The full PEP725 dataset contains millions of observations. Rendering all these points in an interactive map can take long. Restricting the dataset to a species or region of interest results in faster loading and smoother interaction.

Static Maps with pheno_map()

The pheno_map() function creates static maps suitable for reporting and publication. Two background options are supported:

`background`	Description	API Key Required?
`"none"`	Country borders from Natural Earth	No
`"google"`	Google Maps satellite/terrain imagery	Yes

For most applications, the background = "none" option is recommended because it: - requires no external services (API key setup), - works reliably in package vignettes, - produces reproducible, publication-ready maps

Basic Station Maps

# Basic station map with country borders (no API key needed)
pheno_map(pep, background = "none", color_by = "none", point_size = 1.5)

This view is useful for assessing overall network coverage and geographic gaps.

Coloring by Data Attributes

# Color stations by number of observations
pheno_map(pep, background = "none", color_by = "n_obs", point_size = 1.5)

# Color by species diversity at each station
pheno_map(pep, background = "none", color_by = "n_species", point_size = 1.5)

Using Google Maps Background

For satellite imagery (requires API key from Google Cloud Console):

# Register your API key first
ggmap::register_google(key = "your_api_key_here")

# Then use background = "google"
pheno_map(pep, background = "google", color_by = "n_obs", zoom = 5)

# Regional focus with Google Maps
pheno_map(
  pep,
  background = "google",
  location = c(lon = 8.2, lat = 46.8),
  zoom = 7,
  color_by = "n_obs"
)

These maps help identify well-monitored stations and multi-species observation sites.

Mapping Phenological Patterns

Beyond station metadata, pheno_map() can visualize derived phenological quantities, allowing spatial interpretation of biological signals.

Mean Phenological Timing

# Map mean phenological timing for flowering (phase 60)
# Earlier timing = darker colors, later = lighter (plasma scale)
pheno_map(pep, background = "none", color_by = "mean_doy",
        phase_id = 60, point_size = 1.5)

This representation reveals spatial gradients in average phenological timing.

Phenological Trends

# Map trends per station
# Blue = phenology getting earlier, Red = getting later
pheno_map(pep, background = "none", color_by = "trend",
        phase_id = 60, period = 1990:2020, min_years = 10, point_size = 1.5)

Trend maps highlight regions where phenology is advancing or delaying over time.

Species-level Variation

# Map species variation at each station
# Higher CV = more variation in timing among species
pheno_map(pep, background = "none", color_by = "species_cv",
        phase_id = 60, min_species = 3, point_size = 1.5)

High inter-species variability may indicate ecological heterogeneity or mixed habitat conditions.

Understanding `color_by` Options

`color_by`	Meaning	Typical application
`"none"`	Station locations	Observation network overview
`"n_obs"`	Observation density	Identification of well-monitored stations
`"n_species"`	Species richness	Finding multi-species stations
`"mean_doy"`	Mean timing of day-of-year	Reveal spatial phenology gradients
`"trend"`	Temporal change	Climate impact assessment
`"species_cv"`	Inter-species variability	Find stations with consistent timing

Parameter Reference for pheno_map()

Parameter	Description	Default
`background`	Map background: “none” or “google”	“google”
`location`	Map center as c(lon, lat)	European center
`zoom`	Zoom level (4=wide, 7=tight)	4
`color_by`	Coloring option	“none”
`phase_id`	BBCH phase for phenological colors	NULL
`period`	Years for trend calculation	All years
`min_years`	Minimum years for trend	10
`min_species`	Minimum species for CV	3
`point_size`	Marker size	0.8
`output_file`	Path to save map	NULL (display only)

Combining Mapping with Analysis

A typical workflow integrates mapping with spatial analysis:

# 1. Explore stations interactively
selected <- pheno_leaflet(pep[genus == "Malus"])

# 2. Subset data to selected region
pep_region <- pep[s_id %in% selected$s_id]

# 3. Analyze gradients in the selected region
gradient <- pheno_gradient(pep_region, variable = "alt",
                           species = "Malus domestica", phase_id = 60)

# 4. Assess synchrony
synchrony <- pheno_synchrony(pep_region, species = "Malus domestica",
                              phase_id = 60)

# 5. Create publication map of the region (no API needed)
pheno_map(pep_region, background = "none", color_by = "mean_doy",
        phase_id = 60, output_file = "regional_phenology.pdf")

Best Practices Summary

For Gradient Analysis

Check sample size: Need sufficient stations across the gradient range (minimum 20 stations, ideally 50+)
Use robust methods: Phenological data often has outliers from observation errors or microclimate effects - always start with method = "robust"
Consider regional differences: Gradients may vary by location - use the by parameter to check consistency
Validate against literature: Compare your results to published lapse rates (2-4 days/100m for altitude is typical in temperate Europe)
Report uncertainty: Always include R-squared and sample size when reporting gradient estimates

For Synchrony Analysis

Set appropriate minimum stations: Too few stations gives unreliable estimates (minimum 3, ideally 10+)
Consider temporal scale: Annual synchrony captures different patterns than decadal averages
Account for network changes: Station additions/removals can create artificial trends - check if network composition changed
Use appropriate metrics: SD for absolute variability, CV for relative comparisons across phases
Interpret trends cautiously: Distinguish statistically significant trends from ecologically meaningful changes

For Mapping

Start with exploration: Use pheno_leaflet() to understand your data before formal analysis
Choose appropriate zoom: Match zoom level to your study area (4=Europe, 6=country, 8=region)
Use meaningful colors: Choose color_by option that matches your research question
Document selections: Save selected station lists for reproducibility
Export for publication: Use output_file to create high-quality figures

Summary

This vignette demonstrated how pep725 supports spatial phenological analysis through complementary methods:

Topic	Key Function	Main Output
Elevation gradients	`pheno_gradient()`	Days delay per 100m
Latitude gradients	`pheno_gradient()`	Days delay per degree
Spatial synchrony	`pheno_synchrony()`	SD/CV of timing across stations
Interactive mapping	`pheno_leaflet()`	Selected station data.frame
Static mapping	`pheno_map()`	Publication-quality maps

Key Take-Home Messages

Phenology delays systematically with increasing altitude and latitude - quantifying these gradients helps scale point observations to landscapes
Synchrony tells you about spatial variability - high synchrony suggests strong regional climate control, low synchrony indicates local factors matter
Combine approaches for complete understanding - gradients show mean patterns, synchrony shows variability around those patterns
Mapping is essential for exploring data, selecting regions, and communicating results

Next Steps

Explore the other vignettes for complementary analyses:

Getting Started: Data access, the pep class, and basic exploration

vignette("getting-started", package = "pep725")

Phenological Analysis: Normals, anomalies, quality assessment, and trends

vignette("phenological-analysis", package = "pep725")

Data Quality Assessment: Outlier detection and data validation

vignette("data-quality", package = "pep725")

Session Info

sessionInfo()
#> R version 4.5.2 (2025-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.0.1
#> 
#> Matrix products: default
#> BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Zurich
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] pep725_1.0.0        ggplot2_4.0.1       data.table_1.18.2.1
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10             generics_0.1.4          tidyr_1.3.2            
#>  [4] class_7.3-23            robustbase_0.99-7       KernSmooth_2.23-26     
#>  [7] lattice_0.22-7          digest_0.6.39           magrittr_2.0.4         
#> [10] evaluate_1.0.5          grid_4.5.2              RColorBrewer_1.1-3     
#> [13] rnaturalearth_1.2.0     fastmap_1.2.0           Matrix_1.7-4           
#> [16] jsonlite_2.0.0          e1071_1.7-17            DBI_1.2.3              
#> [19] mgcv_1.9-3              purrr_1.2.1             viridisLite_0.4.2      
#> [22] scales_1.4.0            jquerylib_0.1.4         cli_3.6.5              
#> [25] rlang_1.1.7             units_1.0-0             splines_4.5.2          
#> [28] withr_3.0.2             cachem_1.1.0            yaml_2.3.12            
#> [31] otel_0.2.0              tools_4.5.2             dplyr_1.2.0            
#> [34] vctrs_0.7.1             R6_2.6.1                proxy_0.4-29           
#> [37] lifecycle_1.0.5         classInt_0.4-11         pkgconfig_2.0.3        
#> [40] pillar_1.11.1           bslib_0.10.0            gtable_0.3.6           
#> [43] glue_1.8.0              Rcpp_1.1.1              sf_1.0-24              
#> [46] rnaturalearthdata_1.0.0 DEoptimR_1.1-4          xfun_0.56              
#> [49] tibble_3.3.1            tidyselect_1.2.1        rstudioapi_0.18.0      
#> [52] knitr_1.51              farver_2.1.2            nlme_3.1-168           
#> [55] htmltools_0.5.9         patchwork_1.3.2         rmarkdown_2.30         
#> [58] labeling_0.4.3          compiler_4.5.2          S7_0.2.1

University of Applied Sciences Northwestern Switzerland (FHNW)↩︎
Swiss Federal Institute for Forest, Snow and Landscape Research (WSL)↩︎

Spatial Phenological Patterns

Matthias Templ1

Barbara Templ2

2026-03-04