How to correct richness for sampling effort

On this page

Problem statement

A cell with 26,000 records holds more species than a cell with 200, and almost all of that difference is effort. Comparing raw richness between cells compares how hard people looked.

Rarefaction is the standard correction: ask how many species each cell would have shown if the same number of records had been collected in all of them. It needs the abundance of each species in each cell, which a faceted GBIF query returns for free — and it does not fix the deeper problem, which is that a cell with 200 records cannot be rarefied to 1,000.

Quick answer

Rarefy each cell to a common sample size with the analytic formula:

import numpy as np, requests, time
from scipy.special import gammaln

API = "https://api.gbif.org/v1/occurrence/search"

def cell_species(west, south, east, north, **params):
    """Species counts within one cell, straight from GBIF's facet."""
    wkt = (f"POLYGON(({west} {south},{east} {south},{east} {north},"
           f"{west} {north},{west} {south}))")
    for attempt in range(6):
        r = requests.get(API, params={**params, "geometry": wkt, "hasCoordinate": "true",
                                      "facet": "speciesKey", "facetLimit": 1000,
                                      "limit": 0}, timeout=180)
        if r.status_code == 429:
            time.sleep(2 ** attempt)
            continue
        r.raise_for_status()
        body = r.json()
        time.sleep(0.2)
        return np.array([c["count"] for c in body["facets"][0]["counts"]], dtype=float)
    raise RuntimeError("rate limited")

def rarefy(counts, m):
    """Expected species in a random subsample of m records (Hurlbert 1971)."""
    n = counts.sum()
    if m > n:
        return np.nan
    log_missing = (gammaln(n - counts + 1) - gammaln(n - counts - m + 1)
                   - gammaln(n + 1) + gammaln(n - m + 1))
    return float(len(counts) - np.exp(log_missing).sum())

vancouver = cell_species(-123.5, 49.0, -122.5, 49.5, taxonKey=359, country="CA")
interior = cell_species(-126.0, 54.0, -125.0, 55.0, taxonKey=359, country="CA")
for name, c in [("near Vancouver", vancouver), ("central interior", interior)]:
    print(f"{name:18} {int(c.sum()):>7,} records, {len(c):>3} species observed, "
          f"{rarefy(c, 100):>5.1f} expected in 100, {rarefy(c, 200):>5.1f} in 200")

near Vancouver      26,390 records,  94 species observed,  25.4 expected in 100,  31.6 in 200
central interior       150 records,  33 species observed,  28.5 expected in 100,   nan in 200

Observed, the Vancouver cell has nearly three times the species of the interior one. Rarefied to 100 records the ranking reverses: the interior cell would have shown 28.5 species to Vancouver's 25.4. The whole of the apparent difference was effort, and a little more. Rarefied to 200 the sparse cell returns nan, because it only has 150 records — a refusal that is the correct answer, since rarefaction interpolates downwards and cannot invent records a cell never had.

Comparing at the same sample size is the whole idea; comparing at the observed sizes is the mistake.

Abundance distributions, a common sample, the analytic formula

1. Get the abundance distribution, not just the richness

Rarefaction needs how many records each species has in the cell. The facet gives exactly that.

2. Choose the common sample size

Usually the smallest cell you are willing to keep. Cells with fewer records than that cannot be rarefied and must be dropped, not extrapolated.

3. Use the analytic formula, not a simulation

Hurlbert's formula gives the expected richness exactly, in one vectorised line. Resampling gives the same answer with noise and a thousand times the work.

4. Work in log space

The factorials overflow quickly; gammaln keeps the arithmetic finite for cells with tens of thousands of records.

5. Report how many cells you lost

Rarefying to 200 records may drop half the map. That is the real cost of the correction and belongs in the caption.

6. Remember what it does not fix

Rarefaction equalises sample size. It does not fix which places were sampled, or that one platform contributed most of the records.

A curve that is still climbing means the cell has not been surveyed enough to know its richness.

Curves and a corrected grid

The data used on this page

Every number below comes from GBIF's API. GBIF adds and revises records daily, so the setup answers each GBIF request this page makes from the reply it got on 9 October 2026, and lets any other request through live.

import hashlib, io, zipfile

import requests

# GBIF adds and revises records daily. So that the numbers below reproduce, every GBIF request on
# this page is answered from the reply it got on 9 October 2026 (one file per request, named by a
# SHA-1 of its URL); any other request, or all of them once these lines are deleted, goes live.
SAVED = zipfile.ZipFile(io.BytesIO(requests.get(
    "https://data.spatialworkflow.io/gbif-api-2026-10-09.zip", timeout=300).content))
live_get = requests.get


def saved_get(url, params=None, **kwargs):
    prepared = requests.Request("GET", url, params=params).prepare()
    try:
        body = SAVED.read(hashlib.sha1(prepared.url.encode()).hexdigest()[:16] + ".json")
    except KeyError:
        return live_get(url, params=params, **kwargs)
    reply = requests.models.Response()
    reply.status_code, reply._content, reply.url = 200, body, prepared.url
    reply.headers["Content-Type"] = "application/json"
    return reply


requests.get = saved_get

Example 1 — the curve, and what its shape says

for name, c in [("near Vancouver", vancouver), ("central interior", interior)]:
    total = int(c.sum())
    points = [m for m in (10, 50, 100, 200, 500, 1_000, 5_000) if m <= total]
    curve = [rarefy(c, m) for m in points]
    print(f"{name} ({total:,} records, {len(c)} species)")
    for m, s in zip(points, curve):
        bar = "█" * int(30 * s / len(c))
        print(f"  {m:>6,} records -> {s:>5.1f} species  {bar}")
    slope = (curve[-1] - curve[-2]) / (points[-1] - points[-2])
    print(f"  slope at the end: {slope * 1000:.2f} species per 1,000 more records\n")

near Vancouver (26,390 records, 94 species)
      10 records ->   7.6 species  ██
      50 records ->  19.5 species  ██████
     100 records ->  25.4 species  ████████
     200 records ->  31.6 species  ██████████
     500 records ->  40.7 species  ████████████
   1,000 records ->  48.2 species  ███████████████
   5,000 records ->  67.9 species  █████████████████████
  slope at the end: 4.91 species per 1,000 more records

central interior (150 records, 33 species)
      10 records ->   7.7 species  ███████
      50 records ->  21.1 species  ███████████████████
     100 records ->  28.5 species  █████████████████████████
  slope at the end: 146.57 species per 1,000 more records

A curve that has flattened means more records would add few species; one still climbing means the cell's richness is unknown, whatever the observed number says.

Example 2 — a grid, corrected

import pandas as pd

cells = []
for west in np.arange(-128.0, -118.0, 1.0):
    for south in np.arange(49.0, 55.0, 1.0):
        counts = cell_species(round(west, 3), round(south, 3),
                              round(west + 1, 3), round(south + 1, 3),
                              taxonKey=359, country="CA")
        if len(counts):
            cells.append({"west": west, "south": south, "records": int(counts.sum()),
                          "observed": len(counts),
                          "rarefied_200": rarefy(counts, 200),
                          "rarefied_1000": rarefy(counts, 1_000)})
grid = pd.DataFrame(cells)

print(f"{len(grid)} occupied cells")
for col, label in [("observed", "observed richness"),
                   ("rarefied_200", "rarefied to 200"),
                   ("rarefied_1000", "rarefied to 1,000")]:
    ok = grid[col].notna()
    print(f"  {label:<20} {ok.sum():>3} cells, median {grid.loc[ok, col].median():>5.1f}, "
          f"range {grid.loc[ok, col].min():>5.1f}–{grid.loc[ok, col].max():>5.1f}")

60 occupied cells
  observed richness     60 cells, median  41.0, range  12.0– 94.0
  rarefied to 200       46 cells, median  33.4, range  20.9– 45.1
  rarefied to 1,000     19 cells, median  49.3, range  35.1– 69.8

Example 3 — how much of the pattern was effort

from scipy import stats

usable = grid[grid.rarefied_200.notna()]
raw = stats.spearmanr(np.log10(usable.records), usable.observed)
adj = stats.spearmanr(np.log10(usable.records), usable.rarefied_200)
print(f"{len(usable)} cells with at least 200 records")
print(f"correlation with effort, observed richness: {raw.statistic:>5.2f}")
print(f"correlation with effort, rarefied richness: {adj.statistic:>5.2f}")

spread = usable.observed.max() / usable.observed.min()
spread_r = usable.rarefied_200.max() / usable.rarefied_200.min()
print(f"\nrichest/poorest cell, observed: {spread:.1f}×")
print(f"richest/poorest cell, rarefied: {spread_r:.1f}×")

46 cells with at least 200 records
correlation with effort, observed richness:  0.82
correlation with effort, rarefied richness:  0.04

richest/poorest cell, observed: 3.6×
richest/poorest cell, rarefied: 2.2×

The correlation with effort falls from 0.82 to 0.04. That does not show effort stopped mattering: cells that were sampled more are often cells with more habitat and more people, and both of those relate to real richness, so even perfect data would keep some of it. What rarefaction removes is the part that is arithmetic — more draws from the same urn find more species.

Why the analytic formula is worth it

Why the analytic formula is worth the gamma functions

Hurlbert's formula computes, for each species, the probability that a random subsample of m records misses it entirely, and subtracts the sum of those probabilities from the species count. That is exactly what a resampling simulation estimates, without the Monte Carlo error, in one pass. The gammaln form is needed because the binomial coefficients overflow float64 at a few hundred records.

Why cells below the threshold cannot be rescued

Rarefaction interpolates downwards. Asking what a 50-record cell would have shown with 200 records is extrapolation, and the estimators that attempt it — Chao1, Chao2, the coverage-based family — make assumptions about the undetected species that presence-only data cannot check. Dropping the cell is honest; extrapolating it is a model.

Why this is still not comparable richness

Two cells rarefied to the same sample size are comparable in the number of draws, not in what was drawn. One cell's 200 records are a single ten-year camera study; the other's are two hundred people photographing squirrels in a park. The correction handles sample size and nothing else, and saying so is part of using it.

Each rarefaction threshold is a trade between comparability and coverage.

Coverage-based rarefaction, singletons and Chao1

Coverage-based rarefaction (Chao & Jost) equalises sample completeness rather than sample size, and is arguably the better comparison.
Individuals vs records: occurrence records are not individuals, so this is record-based rarefaction; say so.
Duplicates inflate abundances and therefore the curve; clean before rarefying.
Singletons dominate the curve's tail — a cell whose species are nearly all singletons has not been surveyed.
Chao1 estimates asymptotic richness from singletons and doubletons; with presence-only data it is optimistic.
scipy.special.gammaln keeps the formula stable; the factorial form fails above a few hundred records.

Internal links

How to build a species richness grid — the map this corrects
Biodiversity metrics explained — what else richness misses
A richness map is a map of sampling effort — the fix guide
Sampling bias in occurrence data explained — the bias rarefaction cannot touch
How to clean species occurrence records — duplicates distort the curve
Sample design explained: where to measure — what comparable sampling looks like

FAQ

What is rarefaction?

Estimating how many species a sample would have shown at a smaller, common sample size, so that cells with different effort can be compared.

Which sample size should I rarefy to?

The smallest cell you are willing to keep. Every cell below it is dropped.

Can I extrapolate upwards instead?

Estimators exist, and they assume things presence-only data cannot verify. Prefer interpolation and report the dropped cells.

Does rarefaction remove sampling bias?

No. It equalises sample size only. Where people looked is a separate problem.

Do I need to simulate it?

No — Hurlbert's analytic formula gives the exact expectation in one line.

What does a still-rising curve mean?

That the cell's species list is incomplete; its observed richness is a lower bound.