pygeostats

Geostatistics for Python with a Rust-accelerated core. Variograms, kriging, point-pattern analysis and spatial autocorrelation, behind a familiar fit / predict API that accepts NumPy arrays, pandas DataFrames and GeoPandas GeoDataFrames.

Documentation · Quickstart · API reference · Known limitations · Changelog · Issues

Alpha release. 0.1.0a2 is on PyPI as a pre-release. The API may still change, and some features are unfinished or unreliable: read the known limitations before relying on anisotropy analysis.

Highlights

• A compiled core. Distances, empirical variograms, model fitting and kriging run in Rust, built with PyO3 and maturin.
• The full workflow. Estimate a variogram, fit a model, krige with prediction variance, and cross-validate the result, all in one package.
• Honest fits. When the data does not constrain a variogram model, fit() says so through converged_ and warnings_, instead of returning a range as though it were reliable.
• Beyond kriging. Point-pattern statistics, spatial clustering and spatial autocorrelation share the same inputs.
• Wheels for every major platform. One stable-ABI wheel per platform covers Python 3.11 to 3.14, with no Rust toolchain needed to install.
• Tested documentation. Every example in the documentation and in this README runs as part of the test suite.

Installation

pip install pygeostats

pygeostats needs Python 3.11 or newer. While no stable release exists, pip installs the pre-release; after one does, use pip install --pre pygeostats to get pre-releases.

Optional extras:

pip install "pygeostats[plotting]"

Wheels are published for:

Elsewhere, pip builds from the source distribution, which needs a Rust toolchain.

Quick start

Estimate a variogram from scattered samples, fit a model, and krige onto a grid:

import numpy as np
from pygeostats.kriging import OrdinaryKriging
from pygeostats.variogram import EmpiricalVariogram, Variogram

rng = np.random.default_rng(0)
coords = rng.uniform(0, 10, size=(100, 2))
values = np.sin(coords[:, 0]) + np.cos(coords[:, 1]) + rng.normal(0, 0.1, 100)

# 1. Empirical variogram: semivariance by distance bin
empirical = EmpiricalVariogram(coords, values, n_bins=12).compute()

# 2. Fit a model, and check that the fit can be trusted
model = Variogram(model="exponential")
model.fit(empirical.distances_, empirical.gamma_, weights=empirical.counts_)
print(model.converged_, model.nugget_, model.sill_, model.range_)

# 3. Krige onto a grid, with the prediction variance
xs = np.linspace(0, 10, 40)
grid_x, grid_y = np.meshgrid(xs, xs)
targets = np.column_stack([grid_x.ravel(), grid_y.ravel()])

kriging = OrdinaryKriging(model).fit(coords, values)
predictions, variance = kriging.predict(targets, return_variance=True)

Check how well the workflow predicts, holding out whole regions at a time:

from pygeostats.validation import (
    default_kriging_builder,
    default_variogram_builder,
    spatial_kfold_cross_validation,
)

cv = spatial_kfold_cross_validation(
    coords,
    values,
    default_variogram_builder("exponential"),
    default_kriging_builder(),
    n_splits=5,
    random_state=0,
)
print(cv.summary())  # {"rmse": ..., "r2": ...}

Point patterns and spatial autocorrelation work directly on coordinates and values:

from pygeostats import (
    morans_i,
    ripley_l_function,
    simulate_poisson_process,
    spatial_weights_knn,
)

# Ripley's L for a random pattern: close to r at every radius
points = simulate_poisson_process(200, bounds=(0.0, 1.0, 0.0, 1.0), random_state=0)
radii = np.linspace(0.01, 0.15, 15)
l_values = ripley_l_function(points, radii, area=1.0)

# Moran's I for the samples above, with a permutation test
weights = spatial_weights_knn(coords, k=8)
print(morans_i(values, weights, permutations=999, random_state=0))

The Quickstart walks through a complete workflow, including plots, and examples/basic_kriging.py, examples/variogram_fitting.py and examples/anisotropic_kriging_example.py are complete scripts.

Features

Status and known limitations

pygeostats is alpha software. The variogram, kriging, point-pattern, autocorrelation and validation workflows are tested and documented, but:

• Fitting an anisotropic model from directional variograms is not implemented. AnisotropicKriging works with parameters set by hand, and measures its rotation angle clockwise, unlike DirectionalVariogram.
• Anisotropy estimates are rough. The ratio from detect_anisotropy() runs low, and on a few hundred samples the estimated axis can be tens of degrees off.
• Matérn models cannot be fitted, only exponential, spherical and Gaussian ones.
• Type annotations are incomplete, and mypy runs as an advisory CI step.

The known limitations page has the details and workarounds.

Development

A Rust toolchain is needed to build from source. The compiler version is pinned in rust-toolchain.toml, and rustup fetches it automatically.

git clone https://github.com/DiogoRibeiro7/pygeostats.git
cd pygeostats
pip install -e ".[dev,test]"

pytest tests/
black --check src/python/ tests/
ruff check src/python/ tests/
cargo test --no-default-features --features parallel

The development guide explains the checks and how to build the documentation.

Contributing

Bug reports, questions and pull requests are welcome in the issue tracker. See CONTRIBUTING.md for the workflow, and the Code of Conduct.

Name

This project was previously called pyspatialstats. It was renamed because that name belongs to an existing, actively maintained package on PyPI by Jasper Roebroek. The two are unrelated.

License

MIT. See LICENSE.