Cleaning GBIF Occurrence Records with Python

import numpy as np
import pandas as pd

BAD_BASIS = {"FOSSIL_SPECIMEN", "LIVING_SPECIMEN", "MATERIAL_CITATION"}
BAD_ESTABLISHMENT = {"MANAGED", "INTRODUCED", "CULTIVATED"}


def clean_gbif(df: pd.DataFrame, min_year: int = 1970,
               max_uncertainty_m: float = 1000.0,
               min_decimals: int = 3) -> tuple[pd.DataFrame, dict]:
    """Filter a GBIF occurrence export into a modelling-ready presence set."""
    n0 = len(df)
    log = {"raw": n0}

    d = df.copy()
    d = d.dropna(subset=["decimalLatitude", "decimalLongitude"])
    log["has_coords"] = len(d)

    lat, lon = d["decimalLatitude"], d["decimalLongitude"]
    d = d[(lat.between(-90, 90)) & (lon.between(-180, 180))
          & ~((lat == 0) & (lon == 0)) & (lat != lon)]
    log["valid_coords"] = len(d)

    if "coordinateUncertaintyInMeters" in d:
        unc = d["coordinateUncertaintyInMeters"]
        d = d[unc.isna() | (unc <= max_uncertainty_m)]
    log["precise_enough"] = len(d)

    # Coordinates rounded to the degree or the tenth of a degree are label transcriptions.
    def decimals(series):
        s = series.astype(str).str.split(".").str[-1]
        return s.str.len().where(series.astype(str).str.contains(r"\."), 0)

    d = d[(decimals(d["decimalLatitude"]) >= min_decimals)
          & (decimals(d["decimalLongitude"]) >= min_decimals)]
    log["enough_decimals"] = len(d)

    if "basisOfRecord" in d:
        d = d[~d["basisOfRecord"].isin(BAD_BASIS)]
    if "establishmentMeans" in d:
        d = d[~d["establishmentMeans"].str.upper().isin(BAD_ESTABLISHMENT).fillna(False)]
    log["wild_observed"] = len(d)

    if "year" in d:
        d = d[d["year"].between(min_year, 2100)]
    log["in_date_window"] = len(d)

    return d.reset_index(drop=True), log

import numpy as np
import pandas as pd


def flag_pile_ups(d: pd.DataFrame, min_count: int = 20,
                  round_to: int = 5) -> pd.DataFrame:
    """Flag coordinates shared by an implausible number of records."""
    key = (d["decimalLatitude"].round(round_to).astype(str) + "," +
           d["decimalLongitude"].round(round_to).astype(str))
    counts = key.map(key.value_counts())
    out = d.copy()
    out["coord_pile_up"] = counts >= min_count
    return out


def flag_near_points(d: pd.DataFrame, reference: pd.DataFrame,
                     tol_m: float = 2000.0) -> pd.DataFrame:
    """Flag records within tol_m of any reference point (centroids, institutions)."""
    from scipy.spatial import cKDTree

    # Approximate metres per degree; adequate for a 2 km tolerance test.
    lat0 = float(d["decimalLatitude"].mean())
    mx = 111_320.0 * np.cos(np.radians(lat0))
    my = 110_540.0

    ref = np.column_stack([reference["lon"] * mx, reference["lat"] * my])
    pts = np.column_stack([d["decimalLongitude"] * mx, d["decimalLatitude"] * my])
    dist, _ = cKDTree(ref).query(pts)
    out = d.copy()
    out["near_reference"] = dist <= tol_m
    return out

def audit_clean(log: dict, cleaned) -> dict:
    """Report retention at each stage and flag suspicious losses."""
    stages = list(log.items())
    losses = [{"stage": b[0], "removed": a[1] - b[1],
               "removed_pct": round(100.0 * (a[1] - b[1]) / max(a[1], 1), 1)}
              for a, b in zip(stages, stages[1:])]
    return {
        "raw": log["raw"],
        "final": len(cleaned),
        "retained_pct": round(100.0 * len(cleaned) / max(log["raw"], 1), 1),
        "by_stage": losses,
    }


rep = audit_clean(log, cleaned)
assert rep["final"] >= 30, f"too few records left to model: {rep}"
worst = max(rep["by_stage"], key=lambda r: r["removed_pct"])
assert worst["removed_pct"] < 80, \
    f"one filter removed {worst['removed_pct']}% at {worst['stage']} — check it"

import numpy as np
import rasterio


def check_against_predictors(cleaned, raster_path: str) -> dict:
    """Every presence must fall on a valid predictor pixel."""
    with rasterio.open(raster_path) as src:
        coords = list(zip(cleaned["decimalLongitude"], cleaned["decimalLatitude"]))
        vals = np.array([v[0] for v in src.sample(coords)], dtype="float64")
        nodata = src.nodata
    off = ~np.isfinite(vals) | (vals == nodata)
    return {"n": len(vals), "off_raster": int(off.sum()),
            "off_pct": round(100.0 * off.mean(), 2)}


chk = check_against_predictors(cleaned, "bio1.tif")
assert chk["off_pct"] < 5.0, f"presences falling outside the predictors: {chk}"