A Geospatial ML Species Distribution Model for Mushrooms

A species distribution model is only as good as the environment it can see. SporeCast's model sees sixty features per cell, drawn from sixteen public sources whose native resolutions span 10 metres to 20 kilometres, and turns them into a habitat map for twenty mushroom species with a seven-learner ensemble, a stacked blend, and isotonic calibration. This is the modelling layer, from the feature stack up.

The parts of the system that are about data honesty, observer bias in the sightings and spatial leakage in the evaluation, have their own articles. Here I assume the training set is already clean and ask what the model is built from.

Sixty features, twelve families

Three design decisions shape that table more than any individual feature.

Reconcile resolution without pretending

The sources disagree about resolution by three orders of magnitude, and the model has to be honest about that rather than treat a 20-kilometre climate value as if it varied at 10 metres. Two techniques do most of the work. Climate normals are lapse-rate downscaled: each point's monthly temperatures are shifted by the elevation difference between the point and the centre of its climate cell at the standard 6.5 degrees per kilometre, and annual mean, warm-month mean, and growing degree days are recomputed at terrain resolution. Seasonality is deliberately not recomputed, because a uniform shift does not change a standard deviation, and precipitation is deliberately not downscaled, because orographic rainfall needs a windward-leeward model and naive scaling would be worse than the coarse value.

The second technique fixes the single largest resolution loss the pipeline ever had. The serve grid is about seven kilometres per cell, and the early pipeline sampled one point per cell from sources that are natively 10 to 30 metres. A cell contains roughly 328,000 embedding pixels, and the model read one. The fix is a server-side reduction in Earth Engine over each cell's footprint at 30 metres, returning the mean and the per-band spread of the 64-dimensional embedding and, per species, the fraction of the cell that clears the species' similarity threshold, computed on the similarity itself rather than reconstructed from moments, because the mean of a maximum is not the maximum of a mean.

Satellite embeddings as habitat texture

The AlphaEarth annual embeddings are the most important family in the model. Each 10-metre pixel carries a 64-dimensional unit vector that encodes vegetation texture, moisture signature, and seasonality at a resolution no land-cover class can express. Sixteen principal components go into the model as features, and a separate similarity feature is built per species: the embedding's dot product with the centroids of that species' finds, clustered into a few modes with the background centroid subtracted. For the golden chanterelle that similarity is the top feature by permutation importance, ahead of canopy height and soil pH.

Hosts as biology, not correlation

Most of the twenty species are mycorrhizal, and their habitat is a tree before it is a climate. Each species carries an affinity profile over forest community types from the LANDFIRE vegetation layer ▓▓▓ and a distance to its host class, and at serve time the suitability of mycorrhizal species is gated on a host-affinity surface. Two species go further: the lobster mushroom, a parasite, is gated on the modelled suitability of the russulas it lives on, and matsutake on a candystick indicator plant that shares its mycorrhizal network. The gates respect the biology instead of hoping the correlations find it.

Seven learners and a blend

Seven learners with deliberately different inductive biases are trained per species under spatially blocked cross-validation: random forest and extra trees (bagged), histogram gradient boosting, LightGBM and XGBoost (boosted), a MaxEnt-style formulation (quadratic features into a regularised logistic regression, which is the feature-class expansion MaxEnt performs, kept dependency-free), and a plain logistic regression as a smooth statewide anchor. Gradient boosting's hyperparameters are tuned per species with a short Optuna search whose objective is spatially blocked AUC, cached to disk.

The blend is a logistic meta-learner trained on the out-of-fold predictions, with coefficients clipped at zero and normalised, so a member that adds nothing gets exactly zero weight. That happens: for the golden chanterelle the histogram boosting member has an AUC of 0.975 and a stacking weight of zero, because the other boosters carry the same signal. When a species has fewer than forty usable out-of-fold rows the stacker is not trusted and the blend falls back to weights of AUC minus one half, clipped at zero. In the current build all twenty species have enough data to stack.

meta = LogisticRegression(max_iter=1000)
meta.fit(oof[mask], y[mask])
coef = np.clip(meta.coef_[0], 0.0, None)
wnorm = coef / coef.sum() if coef.sum() > 0 else np.ones(len(names)) / len(names)

Calibration is a first-class output

Suitability scores have to be comparable across species, because the map shows the best of the selected species per cell. Each species' blended score is passed through isotonic regression fitted on the out-of-fold predictions with balanced class weights, deliberately against balanced classes rather than true prevalence, so the result is a well-spread relative habitat score rather than the tiny, peaky probability that presence-only data would otherwise imply. The Brier score is tracked before and after: for the median species it falls from 0.144 to 0.101, and for the hedgehog mushroom by nearly half. A calibration curve requires ten rows per bin before it is reported, and both axes use the same weighting, a fix for a real bug where the reported gap scaled with the class weight ratio rather than with any miscalibration.

Uncertainty and importance that respect the folds

• Per-cell uncertainty is the standard deviation of the seven members' predictions, max-normalised, and its per-species mean feeds the confidence badge on the map.
• Confidence intervals on AUC and the Boyce index come from a block bootstrap, resampling whole spatial blocks two hundred times.
• Permutation importance is computed on held-out folds, because in-sample importance on a tree ensemble mostly measures what the trees used to memorise their training points.
• Response curves per feature are stored with the model so the map can show, for a species, how suitability moves with elevation or soil moisture rather than presenting the surface as an oracle.

The grid and the regions

1. The model trains across the wider Pacific Northwest, from northern California to southern British Columbia, and predicts only onto Washington. State lines are political, not ecological, and truncating training data at one would multiply the sample size problem for rare species.
2. The serve grid is 13,103 cells at 0.0625 degrees, about seven kilometres, adaptively densified in two passes: steep terrain gets extra cells, and so do forest edges where a cell's cover disagrees with its neighbour's.
3. Observations and cells are described on identical terms, sampled over the same footprint at the same scale. The pipeline learned that lesson from a feature that once encoded a fingerprint of which points were finds: the asymmetry was the leak.

None of this is exotic on its own. What makes it a model rather than a classifier is the discipline about resolution, leakage, calibration, and biology, applied at every layer. The full write-up, including the pipeline that assembles all of it from keyless public sources, is the SporeCast case study, and the served map is described in the forecasting article. The place to start reading the code is theFEATURE_COLUMNS assembly in the configuration, which is where every decision above becomes a column.