Phenology in R: day of year and event timing

phenology

ggplot2

ecology tutorial

Turn field dates into day of year and measure phenological event timing in R: first, mean, median, peak and duration, and which summary to report.

Author

Tidy Ecology

Published

2026-07-22

Phenology is the study of when biological events happen: when a plant first flowers, when a butterfly first appears, when birds start to breed. In R the raw material is almost always a calendar date, and the first job is to turn that date into a number you can model. That number is the day of year, and this post is about getting it right and then deciding what “when did it happen” even means, because a flowering season is not one date but a whole distribution of dates.

From a date to a day of year

R stores a date as a Date object, and format(x, "%j") returns the day of year as a character string. Wrap it in as.integer() and you have a number from 1 to 365 (or 366) that a regression can use.

d   <- as.Date(c("2023-03-01", "2024-03-01"))
doy <- as.integer(format(d, "%j"))
doy

[1] 60 61

The same calendar date can map to a different day of year. The first of March is day 60 in 2023 but day 61 in 2024, because 2024 is a leap year and the extra 29 February pushes everything after it forward by one. If you are comparing early spring dates across many years, decide up front whether you want the raw day of year (which carries this one-day wobble) or a day count from a fixed anchor. For most spring phenology the wobble is small and harmless, but it is worth knowing it is there.

A season is a distribution, not a date

Imagine monitoring a population of a spring wildflower and recording the day each plant opens its first flower. The result is not a single number but a spread of dates: a few early plants, a peak in the middle, and a long tail of stragglers. Real flowering seasons are usually skewed this way, with the early edge sharper than the late edge.

set.seed(101)
onset <- 105
flow  <- round(onset + rgamma(400, shape = 4.5, scale = 5.5))   # first-flower day per plant

first_d  <- min(flow)
last_d   <- max(flow)
mean_d   <- mean(flow)
median_d <- median(flow)
tab      <- table(flow)
peak_d   <- as.integer(names(tab)[which.max(tab)])              # busiest single day
span     <- last_d - first_d
iqr80    <- as.numeric(quantile(flow, 0.9) - quantile(flow, 0.1))
c(first = first_d, peak = peak_d, median = median_d,
  mean = round(mean_d, 1), last = last_d)

 first   peak median   mean   last 
 110.0  120.0  128.0  130.6  185.0

Five people could summarise this same season with five different numbers. The first date is day 110, the busiest day is day 120, the median (the day by which half the plants have flowered) is day 128, the mean is day 130.6, and the last plant flowers on day 185. Two more numbers describe the spread rather than the centre: the full span from first to last is 75 days, and the interval holding the middle 80 percent of plants is 30.1 days.

marks <- data.frame(
  day = c(first_d, peak_d, median_d, mean_d, last_d),
  lab = c("first", "peak", "median", "mean", "last"),
  col = c(terra, gold, forest, "#16241d", terra))
ggplot(data.frame(flow), aes(flow)) +
  geom_histogram(binwidth = 2, fill = "#c9d3c2", colour = "white") +
  geom_vline(data = marks, aes(xintercept = day, colour = col), linewidth = 0.8) +
  geom_text(data = marks, aes(x = day, y = 34, label = lab, colour = col),
            angle = 90, vjust = -0.3, hjust = 1, size = 3.4) +
  scale_colour_identity() +
  labs(x = "day of year", y = "plants flowering") +
  coord_cartesian(clip = "off") +
  theme(plot.margin = margin(12, 14, 6, 6))

A histogram of flowering dates rising sharply near day 110, peaking near day 120, and trailing off to day 185. Vertical lines mark the first, peak, median, mean and last dates, with mean to the right of median to the right of peak. — Figure 1: One simulated flowering season of 400 plants, with five common phenophase summaries marked. The season is right skewed, so the mean sits to the right of the median, which sits to the right of the peak.

Because the season is right skewed, these summaries fall in a fixed order: the peak (day 120) comes first, then the median (day 128), then the mean (day 130.6). That ordering is not a quirk of this dataset; it is what skewness does, and it means the choice of summary can move your reported date by 10.6 days without anything in the biology changing. Moussus and colleagues compared ten such estimators on simulated seasons and found their accuracy differs a lot, so the summary is a real modelling decision, not a formality.

A cleaner way to see the percentile-based summaries is the cumulative curve: the fraction of plants that have flowered by each day. The median is simply where this curve crosses one half, and the 10th and 90th percentiles bracket the bulk of the season.

qs   <- quantile(flow, c(0.1, 0.5, 0.9))
ecdf_df <- data.frame(day = sort(unique(flow)))
ecdf_df$cum <- ecdf(flow)(ecdf_df$day)
ggplot(ecdf_df, aes(day, cum)) +
  geom_step(colour = forest, linewidth = 0.9) +
  geom_hline(yintercept = c(0.1, 0.5, 0.9), colour = "#b9c2b0", linetype = 2) +
  geom_vline(xintercept = qs, colour = gold, linetype = 3) +
  annotate("text", x = qs, y = c(0.1, 0.5, 0.9) + 0.05,
           label = c("10%", "50% (median)", "90%"), hjust = -0.05, size = 3.3,
           colour = "#46604a") +
  labs(x = "day of year", y = "fraction flowered") +
  scale_y_continuous(limits = c(0, 1))

An S-shaped cumulative curve rising from 0 near day 108 to 1 near day 185. Horizontal dashed lines at 0.1, 0.5 and 0.9 meet the curve at the 10th percentile, median and 90th percentile dates. — Figure 2: The cumulative flowering curve for the same season. The median date is where the curve crosses 0.5; the 10th and 90th percentile dates bracket the middle 80 percent of plants.

When the linear day of year breaks

The day of year works because a spring season lives comfortably in the middle of the calendar, far from the December-January boundary. That is not always true. Consider events clustered around midwinter, some in late December and some in early January.

winter   <- c(350, 355, 358, 360, 3, 6, 10, 12)   # days of year, spanning the boundary
lin_mean <- mean(winter)
ang       <- winter / 365 * 2 * pi
circ_mean <- (atan2(mean(sin(ang)), mean(cos(ang))) / (2 * pi) * 365) %% 365
c(linear_mean = round(lin_mean, 1), circular_mean = round(circ_mean, 1))

  linear_mean circular_mean 
        181.8         364.3

The arithmetic mean of these dates is day 181.8, which lands in the middle of summer, the opposite side of the year from every single observation. The problem is that day 360 and day 5 are five days apart on the calendar but 355 apart as plain numbers. The fix is to treat the dates as angles on a circle: the circular mean is day 364.3, back in late December where it belongs. Any phenological question that wraps around the year, or that asks about time of day, needs this circular treatment rather than a linear day of year. The circular data and the von Mises distribution tutorial covers exactly that machinery.

For a normal spring season none of this bites, and the linear day of year is the right, simple choice. But it is worth carrying the distinction in your head, because the moment your events cross the boundary the arithmetic mean stops meaning anything.

Where to go next

You now have a season represented as a distribution and a menu of summaries for it. The next tutorial takes the most popular summary of all, the first date, and shows that it is the least trustworthy: it moves with how hard you look, not just with the biology. That single fact reshapes how phenological trends should be read.

References

Moussus J-P, Julliard R, Jiguet F 2010. Methods in Ecology and Evolution 1(2):140-150 (10.1111/j.2041-210X.2010.00020.x)

Miller-Rushing AJ, Inouye DW, Primack RB 2008. Journal of Ecology 96(6):1289-1296 (10.1111/j.1365-2745.2008.01436.x)

CaraDonna PJ, Iler AM, Inouye DW 2014. Proceedings of the National Academy of Sciences 111(13):4916-4921 (10.1073/pnas.1323073111)

--- title: "Phenology in R: day of year and event timing" description: "Turn field dates into day of year and measure phenological event timing in R: first, mean, median, peak and duration, and which summary to report." date: "2026-07-22 11:00" categories: [R, phenology, ggplot2, ecology tutorial] image: thumbnail.png image-alt: "A skewed flowering season with first, peak, median and mean dates marked as vertical lines." --- Phenology is the study of when biological events happen: when a plant first flowers, when a butterfly first appears, when birds start to breed. In R the raw material is almost always a calendar date, and the first job is to turn that date into a number you can model. That number is the day of year, and this post is about getting it right and then deciding what "when did it happen" even means, because a flowering season is not one date but a whole distribution of dates. ```{r} #| label: setup #| include: false library(ggplot2) forest <- "#275139"; gold <- "#cda23f"; terra <- "#b5534e"; mown <- "#2f8f63" theme_te <- function(base = 13) { theme_minimal(base_size = base) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6da", linewidth = 0.3), axis.title = element_text(colour = "#46604a"), axis.text = element_text(colour = "#2c3a31"), plot.title = element_text(colour = "#16241d", face = "bold"), plot.subtitle = element_text(colour = "#5d6b61"), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA), legend.position = "none") } theme_set(theme_te()) ``` ## From a date to a day of year R stores a date as a `Date` object, and `format(x, "%j")` returns the day of year as a character string. Wrap it in `as.integer()` and you have a number from 1 to 365 (or 366) that a regression can use. ```{r} #| label: doy d <- as.Date(c("2023-03-01", "2024-03-01")) doy <- as.integer(format(d, "%j")) doy ``` ```{r} #| label: doy-store #| include: false doy23 <- doy[1]; doy24 <- doy[2] ``` The same calendar date can map to a different day of year. The first of March is day `r doy23` in 2023 but day `r doy24` in 2024, because 2024 is a leap year and the extra 29 February pushes everything after it forward by one. If you are comparing early spring dates across many years, decide up front whether you want the raw day of year (which carries this one-day wobble) or a day count from a fixed anchor. For most spring phenology the wobble is small and harmless, but it is worth knowing it is there. ## A season is a distribution, not a date Imagine monitoring a population of a spring wildflower and recording the day each plant opens its first flower. The result is not a single number but a spread of dates: a few early plants, a peak in the middle, and a long tail of stragglers. Real flowering seasons are usually skewed this way, with the early edge sharper than the late edge. ```{r} #| label: season set.seed(101) onset <- 105 flow <- round(onset + rgamma(400, shape = 4.5, scale = 5.5)) # first-flower day per plant first_d <- min(flow) last_d <- max(flow) mean_d <- mean(flow) median_d <- median(flow) tab <- table(flow) peak_d <- as.integer(names(tab)[which.max(tab)]) # busiest single day span <- last_d - first_d iqr80 <- as.numeric(quantile(flow, 0.9) - quantile(flow, 0.1)) c(first = first_d, peak = peak_d, median = median_d, mean = round(mean_d, 1), last = last_d) ``` Five people could summarise this same season with five different numbers. The first date is day `r first_d`, the busiest day is day `r peak_d`, the median (the day by which half the plants have flowered) is day `r median_d`, the mean is day `r round(mean_d, 1)`, and the last plant flowers on day `r last_d`. Two more numbers describe the spread rather than the centre: the full span from first to last is `r span` days, and the interval holding the middle 80 percent of plants is `r round(iqr80, 1)` days. ```{r} #| label: fig-season #| fig-cap: "One simulated flowering season of 400 plants, with five common phenophase summaries marked. The season is right skewed, so the mean sits to the right of the median, which sits to the right of the peak." #| fig-alt: "A histogram of flowering dates rising sharply near day 110, peaking near day 120, and trailing off to day 185. Vertical lines mark the first, peak, median, mean and last dates, with mean to the right of median to the right of peak." #| fig-width: 7 #| fig-height: 4.2 marks <- data.frame( day = c(first_d, peak_d, median_d, mean_d, last_d), lab = c("first", "peak", "median", "mean", "last"), col = c(terra, gold, forest, "#16241d", terra)) ggplot(data.frame(flow), aes(flow)) + geom_histogram(binwidth = 2, fill = "#c9d3c2", colour = "white") + geom_vline(data = marks, aes(xintercept = day, colour = col), linewidth = 0.8) + geom_text(data = marks, aes(x = day, y = 34, label = lab, colour = col), angle = 90, vjust = -0.3, hjust = 1, size = 3.4) + scale_colour_identity() + labs(x = "day of year", y = "plants flowering") + coord_cartesian(clip = "off") + theme(plot.margin = margin(12, 14, 6, 6)) ``` Because the season is right skewed, these summaries fall in a fixed order: the peak (day `r peak_d`) comes first, then the median (day `r median_d`), then the mean (day `r round(mean_d, 1)`). That ordering is not a quirk of this dataset; it is what skewness does, and it means the choice of summary can move your reported date by `r round(mean_d - peak_d, 1)` days without anything in the biology changing. Moussus and colleagues compared ten such estimators on simulated seasons and found their accuracy differs a lot, so the summary is a real modelling decision, not a formality. A cleaner way to see the percentile-based summaries is the cumulative curve: the fraction of plants that have flowered by each day. The median is simply where this curve crosses one half, and the 10th and 90th percentiles bracket the bulk of the season. ```{r} #| label: fig-cumulative #| fig-cap: "The cumulative flowering curve for the same season. The median date is where the curve crosses 0.5; the 10th and 90th percentile dates bracket the middle 80 percent of plants." #| fig-alt: "An S-shaped cumulative curve rising from 0 near day 108 to 1 near day 185. Horizontal dashed lines at 0.1, 0.5 and 0.9 meet the curve at the 10th percentile, median and 90th percentile dates." #| fig-width: 7 #| fig-height: 4.2 qs <- quantile(flow, c(0.1, 0.5, 0.9)) ecdf_df <- data.frame(day = sort(unique(flow))) ecdf_df$cum <- ecdf(flow)(ecdf_df$day) ggplot(ecdf_df, aes(day, cum)) + geom_step(colour = forest, linewidth = 0.9) + geom_hline(yintercept = c(0.1, 0.5, 0.9), colour = "#b9c2b0", linetype = 2) + geom_vline(xintercept = qs, colour = gold, linetype = 3) + annotate("text", x = qs, y = c(0.1, 0.5, 0.9) + 0.05, label = c("10%", "50% (median)", "90%"), hjust = -0.05, size = 3.3, colour = "#46604a") + labs(x = "day of year", y = "fraction flowered") + scale_y_continuous(limits = c(0, 1)) ``` ## When the linear day of year breaks The day of year works because a spring season lives comfortably in the middle of the calendar, far from the December-January boundary. That is not always true. Consider events clustered around midwinter, some in late December and some in early January. ```{r} #| label: winter winter <- c(350, 355, 358, 360, 3, 6, 10, 12) # days of year, spanning the boundary lin_mean <- mean(winter) ang <- winter / 365 * 2 * pi circ_mean <- (atan2(mean(sin(ang)), mean(cos(ang))) / (2 * pi) * 365) %% 365 c(linear_mean = round(lin_mean, 1), circular_mean = round(circ_mean, 1)) ``` The arithmetic mean of these dates is day `r round(lin_mean, 1)`, which lands in the middle of summer, the opposite side of the year from every single observation. The problem is that day 360 and day 5 are five days apart on the calendar but 355 apart as plain numbers. The fix is to treat the dates as angles on a circle: the circular mean is day `r round(circ_mean, 1)`, back in late December where it belongs. Any phenological question that wraps around the year, or that asks about time of day, needs this circular treatment rather than a linear day of year. The [circular data and the von Mises distribution](../circular-data-and-von-mises/) tutorial covers exactly that machinery. For a normal spring season none of this bites, and the linear day of year is the right, simple choice. But it is worth carrying the distinction in your head, because the moment your events cross the boundary the arithmetic mean stops meaning anything. ## Where to go next You now have a season represented as a distribution and a menu of summaries for it. The next tutorial takes the most popular summary of all, the first date, and shows that it is the least trustworthy: it moves with how hard you look, not just with the biology. That single fact reshapes how phenological trends should be read. ## References Moussus J-P, Julliard R, Jiguet F 2010. Methods in Ecology and Evolution 1(2):140-150 (10.1111/j.2041-210X.2010.00020.x) Miller-Rushing AJ, Inouye DW, Primack RB 2008. Journal of Ecology 96(6):1289-1296 (10.1111/j.1365-2745.2008.01436.x) CaraDonna PJ, Iler AM, Inouye DW 2014. Proceedings of the National Academy of Sciences 111(13):4916-4921 (10.1073/pnas.1323073111) ## Related tutorials - [First flowering dates and sampling effort](../first-flowering-date-and-sampling/) - [Phenological trends and temperature](../phenological-trends-and-temperature/) - [Checking a phenology analysis](../checking-a-phenology-analysis/) - [Circular data and the von Mises distribution](../circular-data-and-von-mises/)