Running MaxEnt from Python with elapid

For two decades the standard way to run MaxEnt was to shell out to a Java program, hand it directories of ASCII grids, and parse the results back. A pure-Python implementation removes that boundary: background sampling, fitting, cross-validation and raster prediction all live in one process, which makes the model something a pipeline can call rather than something a person runs. This guide covers doing that with elapid, inside MaxEnt Model Training & Tuning, part of Species Distribution Modeling with MaxEnt.

When to use each implementation

Use elapid for anything automated. Where a result must be reproducible against published Java output, run both once on the same inputs and record the difference — typically a few thousandths of AUC, which is smaller than the variation between cross-validation folds.

Minimal reproducible example

The API follows scikit-learn conventions, so the fit is three lines once the data is assembled. Assembling it is where the geospatial care goes.

import elapid as ela
import geopandas as gpd
import numpy as np

RASTERS = ["bio1.tif", "bio4.tif", "bio12.tif", "bio15.tif", "slope.tif", "twi.tif"]


def fit_maxent(presence_path: str, study_area_path: str,
               n_background: int = 10_000, seed: int = 20260809):
    """Fit a MaxEnt model from a presence layer and a study-area polygon."""
    presences = gpd.read_file(presence_path)
    area = gpd.read_file(study_area_path)

    # Background sampled uniformly across the study area, not across the raster extent.
    background = ela.sample_from_vector(area.geometry.union_all(),
                                        count=n_background, overestimate=2)
    background = gpd.GeoDataFrame(geometry=background, crs=area.crs)

    merged = ela.stack_geodataframes(presences, background, add_class_label=True)
    annotated = ela.annotate(merged, RASTERS, drop_na=True)

    x = annotated[[c for c in annotated.columns if c.startswith("bio") or
                   c in ("slope", "twi")]].to_numpy()
    y = annotated["class"].to_numpy()

    model = ela.MaxentModel(
        feature_types=["linear", "hinge", "product"],
        beta_multiplier=2.0,
        tau=1e-6,
        clamp=True,
        scorer="roc_auc",
    )
    model.fit(x, y)
    return model, annotated

Two choices in there deserve their own paragraph. The background is sampled from the study-area polygon, not from the raster extent — a rectangular extent includes ocean, other countries, and habitat the species could never occupy, and background drawn from it makes the model’s job artificially easy and its AUC artificially high. And drop_na=True removes records that fall on a NoData pixel in any predictor, which is the right behaviour and worth logging, because a silent loss of a third of the presences is otherwise invisible.

Feature classes and what they cost

MaxEnt builds features from the raw predictors: linear, quadratic, product (pairwise interactions), hinge (piecewise linear), and threshold (step). More feature classes mean a more flexible fit and worse transferability.

The conventional rule ties feature richness to sample size: linear only below about 15 presences, linear plus quadratic to 50, add hinge to a few hundred, and reserve the full set for very large datasets. Applied work with 40 forest tree records has no business fitting product and threshold features, however good the resulting AUC looks.

Predict to a raster

import elapid as ela


def predict_raster(model, raster_paths: list[str], out_path: str,
                   transform: str = "cloglog", windowed: bool = True) -> str:
    """Apply a fitted model across the predictor stack, window by window."""
    ela.apply_model_to_rasters(
        model,
        raster_paths,
        out_path,
        transform=transform,
        template_idx=0,
        windowed=windowed,
        nodata=-9999,
        driver="GTiff",
        compress="deflate",
    )
    return out_path

transform="cloglog" is the modern default and the one to use. The older logistic transform assumes a prevalence of 0.5, which is meaningless for presence-only data; cloglog has a cleaner interpretation as a relative occurrence probability and is what current practice reports. Mixing the two between a model and its published comparison is a common source of apparent disagreement.

windowed=True matters at scale: a continental stack at 1 km is tens of millions of pixels per band, and the windowed path keeps memory bounded.

Parameter reference

Tuning beta_multiplier is not optional and is covered in Tuning MaxEnt Regularization with betamultiplier Grid Search. The default of 1.0 routinely overfits small occurrence sets.

Expected output and verification

Before believing any map, check the fit against a spatially blocked split rather than a random one.

import numpy as np
from sklearn.metrics import roc_auc_score


def evaluate_blocked(model_factory, x, y, block_id) -> dict:
    """Leave-one-block-out evaluation; block_id assigns each record to a spatial block."""
    aucs = []
    for b in np.unique(block_id):
        train, test = block_id != b, block_id == b
        if y[test].sum() < 3 or y[train].sum() < 10:
            continue
        m = model_factory()
        m.fit(x[train], y[train])
        aucs.append(roc_auc_score(y[test], m.predict(x[test])))
    return {
        "n_folds": len(aucs),
        "auc_mean": round(float(np.mean(aucs)), 3),
        "auc_sd": round(float(np.std(aucs)), 3),
        "auc_min": round(float(np.min(aucs)), 3),
    }


res = evaluate_blocked(lambda: ela.MaxentModel(feature_types=["linear", "hinge"],
                                               beta_multiplier=2.0), x, y, blocks)
assert res["auc_mean"] > 0.7, f"blocked AUC too low to use: {res}"
assert res["auc_sd"] < 0.15, f"unstable across blocks — check for a spatial outlier: {res}"

A large fold-to-fold spread is more informative than a low mean. It usually means one block contains a distinct habitat that the rest of the data does not represent, which is exactly the situation in which a projected map will be confidently wrong in that region.

The prediction raster itself needs two checks that catch the plumbing errors:

import numpy as np
import rasterio


def check_prediction(out_path: str, presences, background) -> dict:
    """Suitability must be bounded, and higher at presences than at background."""
    with rasterio.open(out_path) as src:
        arr = src.read(1, masked=True)
        p = np.array([v[0] for v in src.sample(
            zip(presences.geometry.x, presences.geometry.y))])
        b = np.array([v[0] for v in src.sample(
            zip(background.geometry.x, background.geometry.y))])
    return {
        "min": round(float(arr.min()), 4),
        "max": round(float(arr.max()), 4),
        "median_at_presence": round(float(np.nanmedian(p)), 4),
        "median_at_background": round(float(np.nanmedian(b)), 4),
    }


chk = check_prediction("suitability.tif", presences, background)
assert 0.0 <= chk["min"] and chk["max"] <= 1.0, f"cloglog output outside 0-1: {chk}"
assert chk["median_at_presence"] > chk["median_at_background"], \
    f"presences score no higher than background — the fit or the sampling is broken: {chk}"

The second assertion is the one worth keeping permanently. It feeds the finished raster the simplest possible test — do the places the species was actually recorded score higher than random ground — and a model that fails it has a wiring error, not a subtle statistical problem.

Common pitfalls

• Background from the raster extent. Including ocean and unreachable regions inflates AUC and distorts every response curve.
• The default regularisation. beta_multiplier=1.0 overfits small occurrence sets reliably; tune it against spatially blocked folds.
• Feature richness beyond the sample size. Product and threshold features on 40 records fit sampling noise and transfer badly.
• Mixing cloglog and logistic outputs. They are different transforms of the same model; thresholds and comparisons do not carry between them.
• An unrecorded background seed. Background sampling is random, so an unseeded run is not reproducible even with identical inputs.
• Ignoring drop_na losses. Presences falling on NoData disappear silently; log the count before and after annotation.

Frequently Asked Questions

Related

• MaxEnt Model Training & Tuning — the wider fitting workflow
• Tuning MaxEnt Regularization with betamultiplier Grid Search — choosing the regularisation this fit uses
• Choosing MaxEnt Feature Classes for Small Samples — matching feature richness to sample size
• Spatial Block Cross-Validation for SDM in Python — the evaluation this page assumes

Up: MaxEnt Model Training & Tuning · Species Distribution Modeling with MaxEnt