Presence-Only Data Preparation for MaxEnt in Python

Presence-only data forms the foundational input layer for robust ecological modeling. Unlike stratified presence-absence surveys that rely on randomized field plots, opportunistic records from herbarium archives, timber stand inventories, and biodiversity aggregators lack verified non-detection points. This structural absence introduces distinct spatial and statistical artifacts — spurious clusters along road networks, county-centroid coordinates, marine spillover, and duplicate geometries — that silently inflate model performance if they reach the algorithm unfiltered. A rigorously executed presence-only preparation workflow ensures downstream algorithms receive spatially unbiased, topologically valid coordinates aligned with rasterized ecological predictors. This page is one stage of the broader Species Distribution Modeling with MaxEnt pipeline, which orchestrates the upstream acquisition and downstream training steps this preparation feeds.

A worked scenario anchors the rest of the page: you have ~18,000 occurrence records for a target tree species pulled from GBIF and a regional forest inventory, exported as a Darwin Core CSV, and a 1 km bioclimatic predictor stack in an equal-area projection. The goal is a deduplicated, thinned, raster-masked occurrence set in which every retained point resolves to one complete environmental feature vector.

Clean, deduplicated set — every retained point maps to one complete environmental feature vector, ready for MaxEnt.

Prerequisites

Before running the pipeline, confirm each of the following. A failure in any one item produces a silently corrupted occurrence set rather than an explicit error.

• Python 3.10+ with pandas ≥ 2.0, geopandas ≥ 0.14, shapely ≥ 2.0, rasterio ≥ 1.3, scipy ≥ 1.11, and numpy ≥ 1.24 installed in a clean environment.
• GDAL/PROJ supplied through the rasterio and geopandas binary wheels — avoid mixing a system GDAL with pip-installed bindings, which produces subtle reprojection offsets.
• Occurrence records exported in Darwin Core terms (decimalLatitude, decimalLongitude, coordinateUncertaintyInMeters, eventDate) or mappable to them.
• A terrestrial land mask (e.g. Natural Earth land polygons) loaded as a GeoDataFrame for marine-spillover removal.
• The environmental predictor stack already harmonized via Environmental Predictor Stacking, with one authoritative CRS, a known cell size, and an explicit nodata value (not None).
• The predictor CRS is equal-area with square pixels, so spatial thinning distances and effort surfaces are measured in true metres rather than degrees.

Concept: why presence-only records need spatial correction

MaxEnt estimates the ratio of the conditional density of environmental covariates at presence sites to their density across the landscape background. When occurrences are collected opportunistically, observed point density is the product of two terms — true habitat suitability and observer effort:

$${\lambda }_{\text{obs}}(s)={\lambda }_{\text{true}}(s)\cdot b(s)$$

where ${\lambda }_{\text{obs}}(s)$ is the observed intensity at location $s$, ${\lambda }_{\text{true}}(s)$ is the ecologically meaningful suitability, and $b(s)$ is the sampling-effort bias field. Because $b(s)$ is typically high near roads, trailheads, and research stations, an uncorrected model fits the gradient of access as if it were a gradient of suitability. Two corrections recover ${\lambda }_{\text{true}}$: spatial thinning, which enforces a minimum inter-point distance to break the autocorrelation that makes dense clusters dominate the likelihood; and effort-matched background sampling, which draws background points from the same $b(s)$ surface so the bias term cancels in the density ratio. Both are developed in depth in handling sampling bias in presence-only data.

Step-by-step Python pipeline

Step 1 — Data ingestion and coordinate standardization

Raw occurrence datasets typically arrive as tabular files (CSV, Excel, or Darwin Core archives) containing latitude, longitude, collection dates, and metadata fields such as coordinate uncertainty and observer identifiers. The first step converts these records into a consistent geographic coordinate reference system (CRS), typically EPSG:4326, and removes records with missing coordinates, zero-valued latitudes or longitudes, or coordinates falling outside the continental landmass. Cross-referencing coordinate pairs against a land mask using shapely geometry operations eliminates marine or offshore artifacts that frequently contaminate terrestrial forestry datasets.

import pandas as pd
import geopandas as gpd


def standardize_occurrences(csv_path: str, land_mask: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Ingest tabular occurrence data, validate coordinates, and filter to land."""
    df = pd.read_csv(csv_path)

    # Coerce to numeric and drop rows with missing or invalid coordinates.
    df["decimalLatitude"] = pd.to_numeric(df["decimalLatitude"], errors="coerce")
    df["decimalLongitude"] = pd.to_numeric(df["decimalLongitude"], errors="coerce")
    df = df.dropna(subset=["decimalLatitude", "decimalLongitude"]).copy()

    # Remove exact-zero artifacts (a common null-island placeholder) and out-of-range values.
    df = df[(df["decimalLatitude"] != 0) | (df["decimalLongitude"] != 0)]
    df = df[df["decimalLatitude"].between(-90, 90) & df["decimalLongitude"].between(-180, 180)]

    gdf = gpd.GeoDataFrame(
        df,
        geometry=gpd.points_from_xy(df["decimalLongitude"], df["decimalLatitude"]),
        crs="EPSG:4326",
    )

    # Spatial join against the land mask removes marine / offshore spillover.
    land_mask = land_mask.to_crs(gdf.crs)
    terrestrial = gdf.sjoin(land_mask, how="inner", predicate="intersects")
    return terrestrial.drop(columns=[c for c in terrestrial.columns if c.startswith("index_right")])

For comprehensive spatial data manipulation standards, consult the official GeoPandas documentation.

Step 2 — Spatial accuracy and uncertainty filtering

Many legacy forestry records and early biodiversity databases report coordinates at coarse resolutions, sometimes aggregated to county centroids or ten-kilometre grid cells. Records whose reported uncertainty (coordinateUncertaintyInMeters) exceeds the native raster cell size of the predictor stack cannot be reliably attributed to a specific environmental pixel and introduce irreducible positional error. Retain only observations that meet a threshold compatible with the raster resolution — for a 1 km stack, keep records with uncertainty ≤ 1000 m and treat missing uncertainty conservatively.

def filter_by_uncertainty(
    gdf: gpd.GeoDataFrame,
    cell_size_m: float,
    keep_missing: bool = False,
) -> gpd.GeoDataFrame:
    """Drop records whose positional uncertainty exceeds one predictor cell."""
    unc = pd.to_numeric(gdf.get("coordinateUncertaintyInMeters"), errors="coerce")
    within = unc <= cell_size_m
    if keep_missing:
        within = within | unc.isna()
    return gdf[within.fillna(False)].copy()

Step 3 — Spatial thinning and effort-matched background

Roadside surveys, accessible trail networks, and proximity to research stations create dense spatial clusters that artificially inflate model performance. Spatial thinning enforces a minimum inter-point distance, breaking the autocorrelation described in the concept section above. A scipy.spatial.cKDTree over projected coordinates gives a fast, greedy thinning pass.

import numpy as np
from scipy.spatial import cKDTree


def spatial_thin(gdf: gpd.GeoDataFrame, min_dist_m: float, seed: int = 0) -> gpd.GeoDataFrame:
    """Greedily retain points so no two survivors are closer than min_dist_m."""
    proj = gdf.to_crs(gdf.estimate_utm_crs())
    coords = np.column_stack([proj.geometry.x, proj.geometry.y])

    rng = np.random.default_rng(seed)
    order = rng.permutation(len(coords))  # randomized order avoids row-order bias
    tree = cKDTree(coords)

    keep = np.zeros(len(coords), dtype=bool)
    suppressed = np.zeros(len(coords), dtype=bool)
    for idx in order:
        if suppressed[idx]:
            continue
        keep[idx] = True
        neighbors = tree.query_ball_point(coords[idx], r=min_dist_m)
        for n in neighbors:
            if n != idx:
                suppressed[n] = True
    return gdf.iloc[keep].copy()

Thinning removes redundancy but does not by itself neutralize effort gradients. The companion step, target-group background sampling — drawing background points from all records of the same taxonomic group so observer effort cancels in the density ratio — is covered in full under handling sampling bias in presence-only data. Adhering to established data quality frameworks such as the GBIF Data Quality Guidelines keeps bias correction aligned with global biodiversity informatics standards.

Step 4 — Raster alignment and predictor masking

Once spatially validated and thinned, occurrence coordinates must be projected into the exact CRS of the predictor stack before feature extraction; mismatched projections or extents cause silent failures. Mask the occurrence layer to the valid-data extent of the stack, discarding points that fall on nodata or ocean-masked cells, so every retained occurrence maps to a complete, non-null environmental feature vector.

import numpy as np
import rasterio


def filter_to_valid_stack_extent(occ_gdf: gpd.GeoDataFrame, stack_path: str) -> gpd.GeoDataFrame:
    """Remove occurrences that land on nodata pixels of the predictor stack."""
    with rasterio.open(stack_path) as src:
        occ_proj = occ_gdf.to_crs(src.crs)
        coords = [(geom.x, geom.y) for geom in occ_proj.geometry]
        sampled = list(src.sample(coords))  # one tuple per point, value per band
        nodata = src.nodata
        valid = np.array([
            all(np.isfinite(v) and (nodata is None or v != nodata) for v in vals)
            for vals in sampled
        ])
    return occ_gdf[valid].copy()

Step 5 — Deduplication, logging, and export

The final step exports a clean GeoDataFrame optimized for algorithmic ingestion: drop duplicate geometries (multiple records snapping to the same predictor cell add no information and re-weight the likelihood), standardize column names to Darwin Core where possible, and append a processing-log column documenting which filters each record survived. This curated set feeds directly into MaxEnt model training and tuning.

def deduplicate_and_export(gdf: gpd.GeoDataFrame, out_path: str, stage: str) -> gpd.GeoDataFrame:
    """Snap to unique geometries, log the surviving stage, and write GeoPackage."""
    gdf = gdf.copy()
    gdf["wkt_key"] = gdf.geometry.apply(lambda g: g.wkb_hex)
    gdf = gdf.drop_duplicates(subset="wkt_key").drop(columns="wkt_key")
    gdf["prep_stage"] = stage  # provenance: last filter this record passed
    gdf.to_file(out_path, layer="occurrences", driver="GPKG")
    return gdf

PDAL or library configuration

The pipeline reads its thresholds from a single declarative config so a run is reproducible and auditable. Keeping the parameters out of the code lets a reviewer confirm — and a CI job re-run — exactly which cell size, thinning distance, and uncertainty cap produced a given occurrence set.

# prep_config.yaml — one source of truth for a preparation run
occurrence_csv: data/raw/gbif_export.csv
predictor_stack: data/stack/bioclim_1km.tif      # authoritative CRS + nodata
land_mask: data/masks/ne_10m_land.gpkg

coordinate:
  target_crs: EPSG:4326          # ingestion CRS; reprojected to stack CRS at masking
  drop_null_island: true         # remove exact (0, 0) placeholders

uncertainty:
  cell_size_m: 1000              # match predictor resolution
  keep_missing: false            # treat unknown uncertainty as failing

thinning:
  min_dist_m: 1000               # >= one predictor cell; break autocorrelation
  seed: 42                       # deterministic survivor selection

export:
  output: data/clean/occurrences_prepared.gpkg
  stage_label: prepared_v1

Validation and verification

Treat preparation as a stage with explicit pass criteria, not a one-shot script. After a run, assert that no point landed on nodata, that the minimum pairwise distance respects the thinning radius, and that geometries are unique. These checks catch the failure modes below before they reach training.

import numpy as np
from scipy.spatial import cKDTree


def verify_prepared_set(gdf, stack_path, min_dist_m):
    # 1. Every point resolves to a complete, non-null feature vector.
    assert len(filter_to_valid_stack_extent(gdf, stack_path)) == len(gdf), "nodata point survived"

    # 2. Thinning radius honored (allow tiny float tolerance).
    proj = gdf.to_crs(gdf.estimate_utm_crs())
    coords = np.column_stack([proj.geometry.x, proj.geometry.y])
    if len(coords) > 1:
        dist, _ = cKDTree(coords).query(coords, k=2)
        assert dist[:, 1].min() >= min_dist_m - 1e-6, "points closer than thinning radius"

    # 3. No duplicate geometries remain.
    assert not gdf.geometry.apply(lambda g: g.wkb_hex).duplicated().any(), "duplicate geometry"
    print(f"OK: {len(gdf)} records verified")

A useful sanity check beyond the assertions: plot the prepared points over the predictor extent and confirm the spatial pattern no longer traces the road network. A residual train-versus-test AUC gap during model validation and AUC metrics is the clearest signal that thinning or background sampling was too weak.

Failure modes and gotchas

• CRS mismatch at masking. Sampling the stack with EPSG:4326 coordinates when the raster is in an equal-area CRS places points kilometres off; always reproject occurrences to src.crs immediately before src.sample.
• nodata left as None. If the predictor stack has no declared nodata, ocean or fill cells return real-looking numbers and survive Step 4 — set an explicit nodata before stacking.
• Thinning in degrees. Running cKDTree on raw lat/long treats one degree of longitude as constant length; project to UTM (or the stack’s equal-area CRS) so min_dist_m is true metres.
• Row-order thinning bias. A non-randomized greedy pass systematically keeps whichever record happens to come first; seed and permute the order so survivor selection is reproducible but unbiased.
• Dropping uncertainty rows you needed. Setting keep_missing=False can discard an entire historical archive that simply never recorded coordinateUncertaintyInMeters; inspect the null rate before choosing.
• Deduplicating before thinning. Removing exact-coordinate duplicates is fine, but thin only after uncertainty filtering, or you may keep a coarse low-quality point and suppress a precise neighbour.

Performance and scale notes

• Chunk the ingest. For multi-million-row Darwin Core archives, read the CSV with chunksize and apply Steps 1–2 per chunk before concatenating; coordinate validation is embarrassingly parallel and keeps peak memory bounded.
• Tile the masking. When the predictor stack is larger than RAM, process occurrences grouped by raster window (or use a windowed rasterio read) rather than loading the full array; only the sampled cells are needed.
• Vectorize the land join. Build a spatial index on the land mask (land_mask.sindex) before sjoin; for continental extents this turns an O(n·m) scan into a near-linear lookup.
• Cache the projection. estimate_utm_crs and to_crs are repeated across steps — compute the working CRS once and reuse it so thinning and masking share a single reprojection.

Frequently Asked Questions

Related

• Handling Sampling Bias in Presence-Only Data — quantify effort and build a proportional MaxEnt background.
• Environmental Predictor Stacking — produce the aligned, nodata-defined raster stack this pipeline masks against.
• MaxEnt Model Training and Tuning — regularization and feature-class selection on the prepared occurrence set.
• Model Validation and AUC Metrics — spatially blocked evaluation that confirms bias correction held.

Up: Species Distribution Modeling with MaxEnt