PhenoVision: Machine Annotations of Phenology for iNaturalist Plant Photos¶

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

Global flower and fruit presence data produced by automated labelling of iNaturalist images up to March 2024, with per-observation detection metrics.

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

⚠️ Large dataset (12.9 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 actually in the file¶

This dataset is a single ~12 GB CSV (annotations_all_w_headers_9cf8ad8.csv) with about 27 million rows — one row per iNaturalist photo in which a deep-learning model detected a reproductive structure. Every row is a positive detection (prediction_class = Detected, certainty = High), so the table tells us where and when plants were photographed flowering or fruiting.

We never load it whole. We use DuckDB to read straight from disk and only pull small aggregated tables into memory. First, let's look at the columns and a few sample rows with a LIMIT query (no full scan).

In [1]:
import duckdb
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style="whitegrid")
CSV = "/srv/data/phenovision/annotations_all_w_headers_9cf8ad8.csv"

con = duckdb.connect()
con.execute("PRAGMA threads=4")

# Column names + types, read from the header only (cheap)
schema = con.execute(f"DESCRIBE SELECT * FROM read_csv_auto('{CSV}') LIMIT 1").df()
print(schema[["column_name", "column_type"]].to_string(index=False))

                            column_name column_type
         machine_learning_annotation_id     VARCHAR
                             datasource     VARCHAR
                          verbatim_date        DATE
                            day_of_year      BIGINT
                                   year      BIGINT
                               latitude      DOUBLE
                              longitude      DOUBLE
          coordinate_uncertainty_meters     VARCHAR
                                 family     VARCHAR
                           count_family     VARCHAR
                                  genus     VARCHAR
                        scientific_name     VARCHAR
                             taxon_rank     VARCHAR
                        basis_of_record     VARCHAR
                                  trait     VARCHAR
                    observed_image_guid     VARCHAR
                     observed_image_url     VARCHAR
                  observed_metadata_url     VARCHAR
                              certainty     VARCHAR
                              model_uri     VARCHAR
                 prediction_probability      DOUBLE
                       prediction_class     VARCHAR
        proportion_low_certainty_family     VARCHAR
accuracy_excluding_low_certainty_family     VARCHAR
                        accuracy_family     VARCHAR

A quick peek at a handful of rows so we can see what the values look like:

In [2]:
sample = con.execute(f"""
    SELECT scientific_name, family, trait, year, day_of_year,
           latitude, longitude, prediction_probability
    FROM read_csv_auto('{CSV}')
    LIMIT 5
""".format(CSV=CSV)).df()
sample
Out[2]:
scientific_name family trait year day_of_year latitude longitude prediction_probability
0 Leucadendron grandiflorum Proteaceae flower 1806 365 -34.007000 18.450000 0.965820
1 Dalbergia purpurascens Fabaceae flower 1893 181 -12.292183 49.209821 0.895508
2 Dalbergia baronii Fabaceae fruit 1893 181 -17.626906 49.570417 0.949219
3 Dalbergia pervillei Fabaceae fruit 1893 199 -17.213798 44.611626 0.943359
4 Dalbergia pervillei Fabaceae fruit 1893 199 -17.213798 44.611626 0.805176

The scientific question: phenology¶

Phenology is the study of the timing of recurring biological events — here, when plants flower and when they fruit. Two questions drive our analysis:

  1. Across the calendar (day_of_year), when do flowers vs. fruits peak — and does the pattern flip between the Northern and Southern hemispheres?
  2. Does the flowering "wave" shift with latitude (tropics vs. temperate zones)?

We answer both from three columns — day_of_year, latitude, trait — aggregated on disk with DuckDB. The query below groups 27 M rows into a compact table in ~20 seconds and returns only a few thousand rows.

In [3]:
phen = con.execute(f"""
    SELECT (CASE WHEN latitude >= 0 THEN 'Northern' ELSE 'Southern' END) AS hemisphere,
           CAST(floor(latitude/10)*10 AS INTEGER)                        AS lat_band,
           trait,
           day_of_year,
           count(*)                                                      AS n
    FROM read_csv_auto('{CSV}')
    WHERE latitude BETWEEN -90 AND 90
      AND day_of_year BETWEEN 1 AND 366
    GROUP BY 1, 2, 3, 4
""".format(CSV=CSV)).df()

# Attach a calendar month to each day-of-year (2000 is a leap year, so day 366 is valid)
base = pd.Timestamp("2000-01-01")
phen["month"] = (base + pd.to_timedelta(phen["day_of_year"] - 1, unit="D")).dt.month

print(f"aggregated rows: {len(phen):,}   detections covered: {phen['n'].sum():,}")
phen.head()

aggregated rows: 9,567   detections covered: 26,968,229
Out[3]:
hemisphere lat_band trait day_of_year n month
0 Northern 10 flower 18 1363 1
1 Northern 60 fruit 20 7 1
2 Northern 40 flower 23 839 1
3 Northern 40 flower 24 639 1
4 Southern -10 flower 24 471 1

1. Flowering vs. fruiting through the year¶

We sum detections by day-of-year for each trait, separately per hemisphere, and smooth with a 15-day rolling mean so the seasonal shape is easy to read. Because the two hemispheres have very different sample sizes, we plot each curve as a share of that trait's annual total — this compares timing, not volume.

In [4]:
fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), sharey=True)

colors = {"flower": "#c94f7c", "fruit": "#e08a2b"}
for ax, hemi in zip(axes, ["Northern", "Southern"]):
    sub = phen[phen["hemisphere"] == hemi]
    for trait in ["flower", "fruit"]:
        s = (sub[sub["trait"] == trait]
             .groupby("day_of_year")["n"].sum()
             .reindex(range(1, 367), fill_value=0))
        s = s / s.sum()                              # share of the year
        s = s.rolling(15, center=True, min_periods=1).mean()
        ax.plot(s.index, s.values, label=trait, color=colors[trait], lw=2)
    ax.set_title(f"{hemi} hemisphere")
    ax.set_xlabel("Day of year")
    ax.set_xlim(1, 366)
axes[0].set_ylabel("Share of that trait's detections")
axes[0].legend(title="trait")
fig.suptitle("Seasonal timing of flowering and fruiting", fontsize=13)
plt.tight_layout()
plt.show()
No description has been provided for this image

Two classic phenology signals jump out:

  • The Northern peak sits around northern spring/summer, while the Southern curve is offset by roughly six months — the seasons are mirror images across the equator.
  • Within each hemisphere, fruiting trails flowering: the orange curve peaks after the pink one, exactly as you'd expect since fruit develops from flowers.

2. The flowering wave across latitude¶

Now we look at flowers only and build a latitude × month heatmap. Each row is a 10-degree latitude band; we normalise each row so it shows when in the year that band flowers most. Reading from top (far north) to bottom (far south), the bright diagonal shows the flowering peak marching through the calendar as we cross the equator.

In [5]:
flowers = phen[phen["trait"] == "flower"]
mat = (flowers.groupby(["lat_band", "month"])["n"].sum()
       .unstack("month").reindex(columns=range(1, 13)).fillna(0))

# Keep well-sampled bands and normalise each row to its own maximum
mat = mat[mat.sum(axis=1) > 5000]
norm = mat.div(mat.max(axis=1), axis=0)
norm = norm.sort_index(ascending=False)            # north at the top

fig, ax = plt.subplots(figsize=(9, 6))
sns.heatmap(norm, cmap="viridis", ax=ax,
            cbar_kws={"label": "Flowering intensity (row-normalised)"})
ax.set_xticklabels(["Jan","Feb","Mar","Apr","May","Jun",
                    "Jul","Aug","Sep","Oct","Nov","Dec"], rotation=0)
ax.set_ylabel("Latitude band (deg)")
ax.set_xlabel("Month")
ax.set_title("Flowering season shifts with latitude")
plt.tight_layout()
plt.show()
No description has been provided for this image

The tropical bands (near 0deg) flower fairly evenly year-round, while temperate bands concentrate their flowering into their local warm season — and the Northern and Southern temperate bands peak half a year apart. This latitudinal phase shift is the fingerprint of solar-driven seasonality captured entirely by citizen photographs.

3. Which plant families dominate, and their flower-to-fruit balance¶

Finally, one more disk scan to rank the most-photographed families and compare how often each is caught flowering vs. fruiting. Families we mostly see in flower vs. in fruit hint at both real biology and photographer/detector bias.

In [6]:
fam = con.execute(f"""
    SELECT family, trait, count(*) AS n
    FROM read_csv_auto('{CSV}')
    WHERE family IS NOT NULL
    GROUP BY 1, 2
""".format(CSV=CSV)).df()

pivot = (fam.pivot_table(index="family", columns="trait",
                         values="n", aggfunc="sum", fill_value=0))
pivot["total"] = pivot.sum(axis=1)
top = pivot.sort_values("total", ascending=False).head(15).iloc[::-1]

fig, ax = plt.subplots(figsize=(9, 6))
ax.barh(top.index, top["flower"], color="#c94f7c", label="flower")
ax.barh(top.index, top["fruit"], left=top["flower"], color="#e08a2b", label="fruit")
ax.set_xlabel("Number of detections")
ax.set_title("Top 15 plant families by detection count")
ax.legend()
ax.xaxis.set_major_formatter(lambda x, _: f"{x/1e6:.1f}M")
plt.tight_layout()
plt.show()

top_ratio = (top["fruit"] / top["total"]).sort_values(ascending=False)
print("Highest fruit share:\n", top_ratio.head(3).round(2).to_string())
print("\nLowest fruit share:\n", top_ratio.tail(3).round(2).to_string())
No description has been provided for this image 
Highest fruit share:
 family
Poaceae       0.66
Solanaceae    0.33
Rosaceae      0.30

Lowest fruit share:
 family
Boraginaceae       0.05
Orchidaceae        0.04
Caryophyllaceae    0.03

Asteraceae (daisies/sunflowers) and Fabaceae (legumes) lead by a wide margin — they are species-rich, widespread, and showy, so people photograph them a lot. The fruit-share numbers reveal ecology too: families known for conspicuous fruits show a higher orange share than those photographed almost exclusively in bloom.

Working with data larger than memory¶

This whole notebook touched a 12 GB file on a 4 GB machine without a single crash, because we never called pd.read_csv on the whole thing. The pattern to remember:

  • Let DuckDB do the heavy lifting on disk and pull back only small aggregates:
    import duckdb
duckdb.sql("SELECT trait, count(*) FROM 'file.csv' GROUP BY trait").df()
  • Read only the columns you need if you must use pandas: pd.read_csv(f, usecols=["trait", "day_of_year"]).
  • Or process in chunks and keep just the running result:
    total = 0
for chunk in pd.read_csv(CSV, usecols=["trait"], chunksize=1_000_000):
    total += len(chunk)

2. The flowering wave across latitude¶

Now we look at flowers only and build a latitude × month heatmap. Each row is a 10-degree latitude band; we normalise each row so it shows when in the year that band flowers most. Reading from top (far north) to bottom (far south), the bright diagonal shows the flowering peak marching through the calendar as we cross the equator.

In [7]:
flowers = phen[phen["trait"] == "flower"]
mat = (flowers.groupby(["lat_band", "month"])["n"].sum()
       .unstack("month").reindex(columns=range(1, 13)).fillna(0))

# Keep well-sampled bands and normalise each row to its own maximum
mat = mat[mat.sum(axis=1) > 5000]
norm = mat.div(mat.max(axis=1), axis=0)
norm = norm.sort_index(ascending=False)            # north at the top

fig, ax = plt.subplots(figsize=(9, 6))
sns.heatmap(norm, cmap="viridis", ax=ax,
            cbar_kws={"label": "Flowering intensity (row-normalised)"})
ax.set_xticklabels(["Jan","Feb","Mar","Apr","May","Jun",
                    "Jul","Aug","Sep","Oct","Nov","Dec"], rotation=0)
ax.set_ylabel("Latitude band (deg)")
ax.set_xlabel("Month")
ax.set_title("Flowering season shifts with latitude")
plt.tight_layout()
plt.show()
No description has been provided for this image