How to extract climate covariates at occurrence points

On this page

Problem statement

A distribution model needs the environment at each record: temperature, precipitation, seasonality, elevation. The layers are global rasters, the records are points, and the extraction is a one-liner — which is why the mistakes in it are quiet ones.

Sampling at the wrong resolution, sampling a point whose coordinate uncertainty is larger than the cell, keeping records that fall in a nodata cell, and feeding nineteen highly correlated bioclimatic variables into a model that assumes they are independent: each of those produces a model that fits and misleads.

Quick answer

Sample every layer at once, and check what came back:

import glob
import numpy as np, pandas as pd, rasterio
from gbif import search

layers = sorted(glob.glob("bio/wc/wc2.1_10m_bio_*.tif"),
                key=lambda p: int(p.split("_")[-1].split(".")[0]))
print(f"{len(layers)} WorldClim bioclimatic layers")

lynx = search(taxonKey=2435263, hasCoordinate="true")
lynx = lynx.dropna(subset=["decimalLongitude", "decimalLatitude"])
coords = list(zip(lynx.decimalLongitude, lynx.decimalLatitude))

values = {}
for path in layers:
    name = "bio" + path.split("_")[-1].split(".")[0]
    with rasterio.open(path) as src:
        sampled = np.array([v[0] for v in src.sample(coords)], dtype="float32")
        sampled[sampled == src.nodata] = np.nan
    values[name] = sampled
env = pd.DataFrame(values)

print(f"{len(env):,} records sampled")
print(f"records with any missing covariate: {env.isna().any(axis=1).sum():,}")
print(env[["bio1", "bio5", "bio6", "bio12"]].describe().loc[
    ["mean", "min", "max"]].round(1).to_string())

19 WorldClim bioclimatic layers
6,353 records sampled
records with any missing covariate: 10
      bio1  bio5  bio6   bio12
mean  -1.2  21.3 -23.0   551.1
min  -11.0   7.6 -34.9   138.0
max   12.1  33.9   1.8  3751.0

src.sample takes an iterable of coordinates and walks the raster once. It is the fast path, and the only thing to remember is that it returns the nodata value rather than a NaN.

Four steps, and three of them are checks.

Raster CRS, nodata as NaN, resolution against uncertainty

1. Match the raster's CRS, not the other way round

Sample in the raster's own coordinate system. Reprojecting a 19-band global raster to sample a few thousand points is the expensive way round.

2. Convert nodata to NaN immediately

rasterio's sample returns the raw value, including the nodata sentinel. In WorldClim that is a large negative float that will quietly become a very cold place.

3. Compare the resolution with the coordinate uncertainty

A 10-minute cell is about 18 km. A record with a 20 km uncertainty samples an arbitrary neighbour of the correct cell.

4. Screen for collinearity before modelling

The nineteen bioclim variables are derived from the same twelve monthly temperature and precipitation series. Several pairs correlate above 0.95.

5. Drop or impute missing values deliberately

Coastal records frequently fall in a sea cell. Dropping them loses coastal species; nearest-land sampling keeps them at the cost of a small position change.

6. Keep the extraction with the model

Layer version, resolution, sampling date. WorldClim 2.1 differs from 1.4 enough to change a model's coefficients.

At 10 arc-minutes the cell is 18 km across; a quarter of these records claim a radius bigger than that.

Nodata costs and collinearity among bioclim variables

The data used on this page

Records come from the GBIF API via gbif.search, the paging function from the download guide plus eventDate and countryCode. Climate is WorldClim 2.1 bioclim at 10′ in bio/wc/; pa.protected_areas() is British Columbia's protected-area layer (EPSG:3005). GBIF changes daily, so the setup answers each request with the reply it got on 9 October 2026.

import hashlib
import io
import os
import 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

os.makedirs("bio/wc", exist_ok=True)
with open("gbif.py", "w") as f:
    f.write('''
import pandas as pd, requests

API = "https://api.gbif.org/v1"
MAX_OFFSET = 100_000
FIELDS = ["key", "species", "speciesKey", "decimalLongitude", "decimalLatitude",
          "coordinateUncertaintyInMeters", "year", "month", "basisOfRecord",
          "datasetKey", "institutionCode", "issues", "eventDate", "countryCode"]


def count(**params):
    return requests.get(f"{API}/occurrence/search", params={**params, "limit": 0}, timeout=120).json()["count"]


def search(limit=None, page=300, **params):
    total = count(**params)
    if limit is None and total > MAX_OFFSET:
        raise ValueError(f"{total:,} records exceeds the {MAX_OFFSET:,} offset limit")
    rows, offset = [], 0
    target = total if limit is None else min(limit, total)
    while offset < target:
        body = requests.get(f"{API}/occurrence/search",
                            params={**params, "limit": min(page, target - offset), "offset": offset},
                            timeout=180).json()
        rows += body["results"]
        offset += len(body["results"])
        if body["endOfRecords"] or not body["results"]:
            break
    return pd.DataFrame([{f: r.get(f) for f in FIELDS} for r in rows])
''')
with open("pa.py", "w") as f:
    f.write('''
import os
import geopandas as gpd

URL = ("https://services6.arcgis.com/ubm4tcTYICKBpist/arcgis/rest/services/"
       "British_Columbia_Parks_Ecological_Reserves_and_Protected_Areas/FeatureServer/0/query"
       "?where=1%3D1&outFields=PROTECTED_LANDS_NAME&outSR=3005&f=geojson")


def protected_areas():
    if not os.path.exists("bio/bc_protected.gpkg"):
        gpd.read_file(URL).to_file("bio/bc_protected.gpkg")
    return gpd.read_file("bio/bc_protected.gpkg")
''')
if not os.path.exists("bio/wc/wc2.1_10m_bio_1.tif"):
    blob = requests.get("https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_bio.zip", timeout=600).content
    zipfile.ZipFile(io.BytesIO(blob)).extractall("bio/wc")
print(len([n for n in os.listdir("bio/wc") if n.endswith(".tif")]), "WorldClim layers")

19 WorldClim layers

Example 1 — what the nodata handling costs

with rasterio.open(layers[0]) as src:
    raw = np.array([v[0] for v in src.sample(coords)], dtype="float32")
    print(f"raster nodata value: {src.nodata}")
    print(f"resolution: {src.res[0] * 60:.0f} arc-minutes "
          f"(~{src.res[0] * 111:.0f} km at the equator)")

missing = raw == -3.4e+38
print(f"records landing on nodata: {missing.sum():,} ({missing.mean():.1%})")
print(f"mean with nodata treated as data: {raw.mean():>12,.1f}")
print(f"mean with nodata masked:          {raw[~missing].mean():>12,.1f} °C")

raster nodata value: -3.3999999521443642e+38
resolution: 10 arc-minutes (~18 km at the equator)
records landing on nodata: 10 (0.2%)
mean with nodata treated as data:         -inf
mean with nodata masked:                  -1.2 °C

Example 2 — collinearity among the nineteen

clean = env.dropna()
corr = clean.corr().abs()
pairs = (corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))
         .stack().sort_values(ascending=False))
print(f"{len(clean):,} complete records")
print("most correlated pairs:")
for (a, b), r in pairs.head(6).items():
    print(f"  {a:>6} / {b:<6} r = {r:.3f}")
print(f"\npairs above 0.9: {(pairs > 0.9).sum()} of {len(pairs)}")

keep, dropped = [], []
for col in clean.columns:
    if all(corr.loc[col, k] < 0.8 for k in keep):
        keep.append(col)
    else:
        dropped.append(col)
print(f"greedy selection at r < 0.8 keeps {len(keep)}: {keep}")

6,343 complete records
most correlated pairs:
   bio14 / bio17  r = 0.995
   bio13 / bio16  r = 0.992
    bio6 / bio11  r = 0.982
    bio4 / bio7   r = 0.966
   bio12 / bio17  r = 0.953
    bio4 / bio6   r = 0.953

pairs above 0.9: 14 of 361
greedy selection at r < 0.8 keeps 9: ['bio1', 'bio2', 'bio3', 'bio4', 'bio5', 'bio8', 'bio9', 'bio12', 'bio15']

Nineteen variables reduce to a handful once correlated pairs are removed. Which handful depends on the order they are considered in, which is an argument for choosing them on ecological grounds rather than by a greedy loop.

Example 3 — uncertainty against cell size

u = lynx.coordinateUncertaintyInMeters.to_numpy()
cell_km = 18.5
print(f"cell size at 10 arc-minutes: about {cell_km:.1f} km")
for label, mask in [("uncertainty <= half a cell", u <= cell_km * 500),
                    ("uncertainty > one cell", u > cell_km * 1000),
                    ("no uncertainty recorded", np.isnan(u))]:
    print(f"  {label:<28} {np.sum(mask):>6,} records ({np.mean(mask):>5.1%})")

rng = np.random.default_rng(0)
usable = np.where(np.isfinite(u) & (u < 50_000))[0]
with rasterio.open(layers[0]) as src:
    here = np.array([v[0] for v in src.sample(
        [coords[i] for i in usable])], dtype="float32")
    angle = rng.uniform(0, 2 * np.pi, len(usable))
    radius = u[usable] * np.sqrt(rng.uniform(0, 1, len(usable)))
    jittered = [(coords[i][0] + (radius[j] * np.cos(angle[j])) / (111_320 * np.cos(
                    np.radians(coords[i][1]))),
                 coords[i][1] + (radius[j] * np.sin(angle[j])) / 111_320)
                for j, i in enumerate(usable)]
    there = np.array([v[0] for v in src.sample(jittered)], dtype="float32")
ok = (here > -1e30) & (there > -1e30)
print(f"\nredrawing each position inside its own uncertainty circle changes "
      f"the sampled annual mean temperature by a median of "
      f"{np.median(np.abs(here[ok] - there[ok])):.2f} °C")
print(f"  90th percentile of the change: "
      f"{np.percentile(np.abs(here[ok] - there[ok]), 90):.2f} °C")

cell size at 10 arc-minutes: about 18.5 km
  uncertainty <= half a cell    2,115 records (33.3%)
  uncertainty > one cell        1,647 records (25.9%)
  no uncertainty recorded       2,217 records (34.9%)

redrawing each position inside its own uncertainty circle changes the sampled annual mean temperature by a median of 0.00 °C
  90th percentile of the change: 0.48 °C

The median change is zero, because most redrawn positions stay inside the same 18 km cell — that is what a step function does. The tail is what matters: a tenth of the records would sample a climate half a degree different, and those are the records with the largest uncertainty, which are also the oldest and the most spatially distinctive.

Sampling is not interpolation

Why sampling is not interpolation

src.sample returns the value of the cell the point falls in — a step function, not a surface. Two records 100 m apart either side of a cell boundary get values that differ by the full between-cell difference. Bilinear interpolation smooths that, at the cost of inventing values the dataset does not contain; for a categorical layer it is simply wrong.

Why the bioclim variables are so correlated

All nineteen are summaries of the same twelve monthly temperature and twelve monthly precipitation grids: annual mean temperature, mean of the warmest quarter, maximum of the warmest month and so on are different reductions of one series. Feeding all nineteen to a model that assumes independent predictors inflates the apparent importance of whatever the correlated block measures.

Why the resolution has to match the question

A 10-minute grid describes climate at about 18 km. It cannot represent a valley, a north-facing slope or an urban heat island, and a model fitted on it predicts at that grain no matter how finely it is rendered. Finer WorldClim grids exist (30 arc-seconds, about 1 km), and they are interpolations from the same weather stations — finer output, not more information.

Nineteen variables, a handful of independent dimensions.

Antimeridian, sea cells and units

Antimeridian and poles: global rasters are in EPSG:4326 and sampling near ±180° needs care.
Sea cells: coastal records commonly land in nodata; decide between dropping and nearest-land sampling.
Elevation from a DEM is a useful extra covariate and is available at far finer resolution than climate.
Future climate layers exist for the same variables, which is how a model is projected forward.
rasterstats handles buffered sampling — the mean within a record's uncertainty radius — if you want to propagate it.
Units: WorldClim 2.1 temperature is in °C, precipitation in mm; version 1.4 used °C × 10.

Internal links

Species distribution models explained — what the covariates feed
How to fit a species distribution model in Python — the model itself
Coordinate uncertainty explained — why the radius matters here
How to extract raster values at point locations — the general operation
A hazard raster's nodata is read as zero — the same sentinel problem
How to thin occurrence records spatially — matching records to the grain
How to count heatwave days — thresholds, baselines and warm nights

FAQ

How do I get climate values at occurrence points?

Open each raster and call src.sample(coords), then replace the nodata value with NaN.

Should I reproject the points or the raster?

The points. Reprojecting a global multi-band raster to sample a few thousand points is much more work.

Why are some values missing?

The point fell on a nodata cell — usually the sea, sometimes outside the layer's extent.

Do I need all nineteen bioclim variables?

No. They are derived from the same monthly series and several pairs correlate above 0.95; choose a subset on ecological grounds.

What resolution should I use?

One that matches the coordinate uncertainty of the records and the grain of the question, not the finest available.

Does the record's uncertainty matter?

Yes — a record with a 20 km radius on an 18 km grid may be sampling a neighbouring cell entirely.