Urban Plant Phenology and Frost Risk (Digitized Herbarium Images)

Category: Botany · Size: 1.1 GB · Format: CSV License: CC-BY-4.0 · Zenodo record · Data sheet on the CSDH

Digitised herbarium image dataset for 200 plant species from the eastern USA, with four reproductive phenological phases plus PRISM climate data and population density.

The data is mounted read-only at /srv/data/urban-plant-phenology/. Save anything you produce in your personal folder (~/).

⚠️ Large dataset (1.1 GB). Your Hub session has 4 GB RAM — do not load the whole file into memory or the kernel will crash. Work like the pros: read only the columns you need, process the file in chunks, or query it straight from disk with DuckDB (no full load). Copy-paste patterns are in "Working with data larger than memory" near the end of this notebook.

What's in the folder

This dataset comes from a crowdsourced herbarium-scoring project: volunteers looked at digitised herbarium specimens of 200 wildflower species from the eastern USA and scored, for each image, whether it showed buds, flowers or fruit. Those votes were then aggregated and joined to PRISM climate data and human population density for the specimen's county.

The folder holds two kinds of file:

• crowdsourcedData_upload.csv (~1 GB) — the raw per-image, per-volunteer votes.
• *_modeling_data (1).csv (4-8 MB) — the analysis-ready tables, one row per specimen, with the day-of-year (doy) a phase was reached plus climate and urbanisation covariates. There is one file per phase: flowering, peakflowering, fruiting, peakfruiting.

We'll do the science on the small modeling tables, and only peek at the big raw file with DuckDB so we never load a gigabyte into memory.

from pathlib import Path

DATA = Path('/srv/data/urban-plant-phenology')

for f in sorted(DATA.glob('*.csv')):
    print(f"{f.name:<42} {f.stat().st_size/1e6:>9,.1f} MB")

crowdsourcedData_upload.csv                  1,044.1 MB
flowering_modeling_data (1).csv                  7.9 MB
fruiting_modeling_data (1).csv                   4.5 MB
peakflowering_modeling_data (1).csv              7.9 MB
peakfruiting_modeling_data (1).csv               4.5 MB

Where the phenology scores come from (peeking at the 1 GB file)

Instead of loading crowdsourcedData_upload.csv, we let DuckDB stream it from disk and return only a one-row summary. This shows the sheer scale of the citizen-science effort behind the tidy tables we analyse below.

import duckdb

provenance = duckdb.sql("""
    SELECT count(*)                 AS n_votes,
           count(DISTINCT coreid)   AS n_images,
           count(DISTINCT workerID) AS n_volunteers,
           count(DISTINCT binomial_species) AS n_species
    FROM read_csv_auto('/srv/data/urban-plant-phenology/crowdsourcedData_upload.csv')
""").df()
provenance

	n_votes	n_images	n_volunteers	n_species
0	3576825	49827	4642	200

Load the analysis-ready flowering table

We focus on first flowering. Each row is one herbarium specimen. Missing values are written as the string "NA", so we tell pandas to treat them as NaN. Key columns:

column	meaning
doy	day of year the specimen was in flower (1-365)
year	collection year (1895-2018)
binomial_species	species (e.g. Arisaema_triphyllum)
tm_spr	mean spring temperature (deg C)
tm_win	mean winter temperature (deg C)
pop_den	human population density (people / km2) - our urbanisation proxy
alt	altitude (m)
county.lat, county.lon	specimen location
tm_spr_anm	spring-temperature anomaly vs the site's long-term average

import pandas as pd
import numpy as np

flr = pd.read_csv(DATA / 'flowering_modeling_data (1).csv', na_values=['NA'])
print('shape:', flr.shape)
print('species:', flr.binomial_species.nunique(), '| years:', flr.year.min(), '-', flr.year.max())
flr[['doy','year','binomial_species','tm_spr','tm_win','pop_den','alt','county.lat','county.lon']].head()

shape: (32135, 22)
species: 200 | years: 1895 - 2018

	doy	year	binomial_species	tm_spr	tm_win	pop_den	alt	county.lat	county.lon
0	235	1895	Datura_stramonium	0.498333	NaN	89.861395	121.216868	41.2735	-73.389
1	192	1895	Lilium_philadelphicum	0.498333	NaN	89.861395	121.216868	41.2735	-73.389
2	127	1895	Uvularia_sessilifolia	0.498333	NaN	89.861395	121.216868	41.2735	-73.389
3	127	1895	Uvularia_sessilifolia	0.498333	NaN	89.861395	121.216868	41.2735	-73.389
4	212	1895	Helianthus_divaricatus	0.498333	NaN	89.861395	121.216868	41.2735	-73.389

A little cleaning: we keep spring/early-summer flowering records (doy between 1 and 250, which drops a handful of odd late-season entries) and rows with the climate fields we need.

flr = flr.dropna(subset=['doy', 'tm_spr', 'tm_win', 'pop_den'])
flr = flr[(flr.doy >= 1) & (flr.doy <= 250)].copy()
print('records kept:', len(flr))

# The best-sampled species carry the clearest signal
counts = flr.binomial_species.value_counts()
counts.head(12)

records kept: 28880

binomial_species
Arisaema_triphyllum       1974
Rudbeckia_hirta           1317
Trillium_erectum           747
Houstonia_caerulea         673
Anemone_quinquefolia       626
Erythronium_americanum     621
Uvularia_sessilifolia      560
Anemone_hepatica           523
Trillium_undulatum         474
Hexastylis_arifolia        417
Uvularia_perfoliata        410
Ruellia_caroliniensis      397
Name: count, dtype: int64

1. When and where do these plants flower?

First, the raw phenology signal. The left panel shows the spread of flowering dates across all species. The right panel maps every specimen by location, coloured by flowering day - a preview of the north-south gradient (plants in the warm south flower earlier).

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style='whitegrid')

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

ax1.hist(flr.doy, bins=60, color='#4c9f70', edgecolor='white')
ax1.set_xlabel('Flowering day of year'); ax1.set_ylabel('Number of specimens')
ax1.set_title('Distribution of flowering dates (200 species)')

sc = ax2.scatter(flr['county.lon'], flr['county.lat'], c=flr.doy, s=6,
                 cmap='viridis_r', alpha=0.5)
ax2.set_xlabel('Longitude'); ax2.set_ylabel('Latitude')
ax2.set_title('Specimen locations, coloured by flowering day')
fig.colorbar(sc, ax=ax2, label='Flowering day of year')
plt.tight_layout(); plt.show()

2. Do plants flower earlier where springs are warmer?

This is the heart of the challenge. We plot flowering day against spring temperature (tm_spr). A hexbin handles the ~29,000 points without turning into a blob, and a fitted line gives the overall temperature sensitivity in days per degree C.

from scipy import stats

x, y = flr.tm_spr.values, flr.doy.values
slope, intercept, r, p, se = stats.linregress(x, y)

fig, ax = plt.subplots(figsize=(9, 6))
hb = ax.hexbin(x, y, gridsize=45, cmap='YlGnBu', mincnt=1)
xs = np.linspace(x.min(), x.max(), 100)
ax.plot(xs, intercept + slope * xs, color='crimson', lw=2.5,
        label=f'{slope:.1f} days per +1 deg C  (r = {r:.2f})')
ax.set_xlabel('Mean spring temperature (deg C)')
ax.set_ylabel('Flowering day of year')
ax.set_title('Warmer springs -> earlier flowering')
ax.legend(loc='upper right')
fig.colorbar(hb, ax=ax, label='specimens per cell')
plt.tight_layout(); plt.show()

print(f'Overall: flowering advances {-slope:.1f} days for every +1 deg C of spring warming '
      f'(p = {p:.1e}, n = {len(flr):,}).')

Overall: flowering advances 2.1 days for every +1 deg C of spring warming (p = 0.0e+00, n = 28,880).

Controlling for species

That overall line mixes 200 species that naturally differ in range and timing. The cleaner test is within each species: fit doy ~ tm_spr separately for every well-sampled species and look at the slopes. If warming advances flowering, almost all slopes should be negative.

well_sampled = counts[counts >= 150].index
rows = []
for sp in well_sampled:
    d = flr[flr.binomial_species == sp]
    s, _, rr, pp, _ = stats.linregress(d.tm_spr, d.doy)
    rows.append({'species': sp, 'n': len(d), 'days_per_degC': s, 'r': rr})
sens = pd.DataFrame(rows).sort_values('days_per_degC')

print(f'{(sens.days_per_degC < 0).mean()*100:.0f}% of {len(sens)} well-sampled species '
      f'flower earlier in warmer springs.')
print(f'Median sensitivity: {sens.days_per_degC.median():.1f} days per deg C.')

# The 15 most temperature-sensitive species
top15 = sens.head(15)
fig, ax = plt.subplots(figsize=(9, 6))
ax.barh(top15.species.str.replace('_', ' '), -top15.days_per_degC, color='#d1495b')
ax.set_xlabel('Flowering advance (days earlier per +1 deg C)')
ax.set_title('Most temperature-sensitive species')
ax.invert_yaxis()
plt.tight_layout(); plt.show()

95% of 61 well-sampled species flower earlier in warmer springs.
Median sensitivity: -2.5 days per deg C.

3. The frost-risk twist

Flowering earlier isn't free. A plant that opens its flowers in early spring can get caught by a late frost that kills the blooms. The species that advance fastest per degree of warming (the spring ephemerals at the top of the chart above - trout lily, Jack-in-the-pulpit, spring beauties) are exactly the ones most exposed as the climate warms.

We can see the raw ingredient of that risk directly: plants that flower early do so when winter has been colder and the ground has barely warmed - i.e. close to frost conditions.

fig, ax = plt.subplots(figsize=(9, 6))
sc = ax.scatter(flr.doy, flr.tm_win, c=flr['county.lat'], cmap='coolwarm_r',
                s=8, alpha=0.4)
# binned mean to show the trend clearly
bins = np.arange(0, 251, 15)
flr['doy_bin'] = pd.cut(flr.doy, bins)
trend = flr.groupby('doy_bin', observed=True).agg(
    doy=('doy', 'mean'), tm_win=('tm_win', 'mean')).dropna()
ax.plot(trend.doy, trend.tm_win, color='black', lw=2.5, marker='o', label='binned mean')
ax.axhline(0, color='steelblue', ls='--', lw=1.5, label='0 deg C (freezing)')
ax.set_xlabel('Flowering day of year')
ax.set_ylabel('Mean winter temperature (deg C)')
ax.set_title('Early flowerers experience colder winters -> higher late-frost exposure')
ax.legend()
fig.colorbar(sc, ax=ax, label='Latitude')
plt.tight_layout(); plt.show()

A note on urbanisation

The challenge also asks whether more urban areas flower earlier (urban heat islands). In this dataset the effect is real but confounded by geography: the densest counties are in the cool north-east, so a naive correlation of doy with pop_den is weak and even points the 'wrong' way. The mechanistic driver is temperature - urban warming matters only insofar as it raises tm_spr. A good next step is to test pop_den within a single species and latitude band.

corr_raw = flr[['doy', 'pop_den']].assign(log_pop=np.log10(flr.pop_den)).corr().loc['doy','log_pop']
print(f'Raw correlation doy vs log10(pop_den): {corr_raw:+.2f}  (weak & confounded by latitude)')
print(f'Correlation doy vs spring temperature:  {r:+.2f}  (strong - the real driver)')

Raw correlation doy vs log10(pop_den): +0.06  (weak & confounded by latitude)
Correlation doy vs spring temperature:  -0.30  (strong - the real driver)

Working with data larger than memory

Your session has a 4 GB RAM limit, but you can analyse files of 10 GB or more without loading them whole:

• Read only the columns you need: pd.read_csv(f, usecols=[...]) / pd.read_parquet(f, columns=[...]).
• Process in chunks and keep only the result:

  total = 0
  for chunk in pd.read_csv(file, chunksize=1_000_000):
      total += len(chunk)
  

• Query with SQL without loading anything — DuckDB (already installed) reads CSV and Parquet straight from disk and only brings the result into memory:

  import duckdb
  duckdb.sql("SELECT column, count(*) FROM '/srv/data/.../file.parquet' GROUP BY column").df()
  

Your turn

This is just the starting point. Some ideas:

• Check the dataset challenge on its CSDH data sheet.
• Work on a copy: right-click the file → Duplicate (or Save Notebook As…). Your changes only live in your Hub space — they're never pushed to GitHub.
• Edited this notebook and want the original back? Use the Restore cell below (or the restore.ipynb notebook).
• Questions and results: on the platform forum.

Attribution: data from Urban Plant Phenology and Frost Risk (Digitized Herbarium Images), license CC-BY-4.0. Notebook from the Citizen Science Data Hub (CSDH) — Fundación Ibercivis.

# ⚠️ RESTORE: this DISCARDS YOUR CHANGES to this notebook and resets it to the original.
# 1. Uncomment the line below (remove the #)   2. Run this cell
# 3. Then: menu File → Reload Notebook from Disk

# !git -C ~/citizen-science-data fetch -q origin && git -C ~/citizen-science-data checkout origin/main -- urban-plant-phenology.ipynb && echo "Restored. Now: File → Reload Notebook from Disk"