wbspatialstats

Pure-Rust spatial statistics and geostatistics library for kriging, variography, spatial autocorrelation, spatial regression, and point-process analysis, designed to serve as the statistical engine for Whitebox.

Mission
The Whitebox Project
Is wbspatialstats Only for Whitebox?
What wbspatialstats Is Not
Design Goals
Current Capabilities
Installation
Quick Start
Examples
Known Limitations
Architecture & Dependencies
Development
License

Mission

Provide robust spatial statistics and geostatistics capabilities for Whitebox applications and workflows.
Implement kriging variants, variogram estimation/fitting, spatial autocorrelation measures, spatial regression, and point-process analysis in pure Rust.
Support production-grade spatial inference with proper diagnostics, cross-validation, and uncertainty quantification.
Prioritize numerical stability and performance while maintaining minimal external dependencies.

The Whitebox Project

Whitebox is a suite of open-source geospatial data analysis software with roots at the University of Guelph, Canada, where Dr. John Lindsay began the project in 2009. Over more than fifteen years it has grown into a widely used platform for geomorphometry, spatial hydrology, LiDAR processing, and remote sensing research. In 2021 Dr. Lindsay and Anthony Francioni founded Whitebox Geospatial Inc. to ensure the project's long-term, sustainable development. Whitebox Next Gen is the current major iteration of that work, and this crate is part of that larger effort.

Whitebox Next Gen is a ground-up redesign that improves on its predecessor in nearly every dimension:

CRS & reprojection — Full read/write of coordinate reference system metadata across raster, vector, and LiDAR data, with multiple resampling methods for raster reprojection.
Raster I/O — More robust GeoTIFF handling (including Cloud-Optimized GeoTIFFs), plus newly supported formats such as GeoPackage Raster and JPEG2000.
Vector I/O — Expanded from Esri Shapefile-only to 11 formats, including GeoPackage, FlatGeobuf, GeoParquet, and other modern interchange formats.
Vector topology — A new, dedicated topology engine (wbtopology) enabling robust overlay, buffering, and related operations.
LiDAR I/O — Full support for LAS 1.0–1.5, LAZ, COPC, E57, and PLY via wblidar, a high-performance, modern LiDAR I/O engine.
Spatial Statistics — A comprehensive spatial statistics engine (wbspatialstats) with kriging, variography, autocorrelation, regression, and point-process tools.
Frontends — Whitebox Workflows for Python (WbW-Python), Whitebox Workflows for R (WbW-R), and a QGIS 4-compliant plugin are in active development.

Is wbspatialstats Only for Whitebox?

No. wbspatialstats is developed primarily to support Whitebox, but it is not restricted to Whitebox projects.

Whitebox-first: API and roadmap decisions prioritize Whitebox geospatial analysis needs.
General-purpose: the crate is usable as a standalone spatial statistics library in other Rust geospatial applications.
Production-ready: proper error handling, numerical stability checks, and comprehensive test coverage make it suitable for production workflows.

What wbspatialstats Is Not

wbspatialstats is a spatial statistics computation engine. It is not a full statistics framework.

Not a format I/O library (file reading/writing belongs in wbraster and wbvector).
Not a rendering or visualization engine.
Not a coordinate reference system engine (CRS and projection belong in wbprojection).
Not a full statistical modeling framework (Bayesian inference, advanced model selection, simulation methods are beyond current scope).
Not a spatial database (point storage and management is the caller's responsibility).

Design Goals

Pure Rust only with no unsafe code in core algorithms.
Keep dependencies minimal; only use nalgebra for linear algebra and thiserror for error handling.
Prioritize numerical stability with appropriate scaling, regularization, and stability checks.
Support both synchronous batch-mode computation and single-prediction workflows.
Interoperability with wbraster and wbvector via coordinate/value tuples and standard Rust collections.

Current Capabilities

Kriging

Deterministic spatial interpolation with uncertainty quantification (kriging variance):

Ordinary Kriging: Local prediction with mean estimation from data
Local Kriging: Neighborhood-based kriging for large datasets
Simple Kriging: Prediction with known, fixed mean
Universal Kriging: Linear trend-surface integration (drift polynomial: constant, linear, quadratic)
Space-Time Kriging: Spatio-temporal prediction with optional temporal autocorrelation weighting
Ordinary CoKriging: Multivariate kriging leveraging auxiliary variables via cross-variograms for improved predictions

Features:

Kriging variance estimation for uncertainty quantification
Cross-covariance matrix computation for error analysis
Neighborhood-size control for scalability
Support for Gaussian and exponential variogram models
Anisotropic kriging support via variogram anisotropy parameters
CoKriging block-structured system matrix for multivariate prediction

Variography

Empirical and theoretical variogram computation for spatial dependence characterization, including anisotropy analysis:

Empirical Variogram: Lag-based semi-variance estimation with configurable bin sizes
Directional Variogram: Azimuthal analysis across multiple directions (0-180°) with tolerance control for detecting and quantifying spatial anisotropy
Anisotropy Modeling: Automatic detection and fitting of directional variation in spatial continuity with anisotropy ratio and principal direction estimation
Cross-Variogram: Spatial dependence between primary and auxiliary variables (required for CoKriging)
Robust Variogram Fitting: Multiple robust loss functions (Cressie, Dowd, Genton)
Variogram Model Families:
- Exponential (isotropic and anisotropic)
- Gaussian (isotropic and anisotropic)
- Linear
- Power
- Spherical

Features:

Cloud-based semi-variance diagnostics (pairwise analysis for outlier detection)
Nested variogram support for complex spatial structures
Automatic nugget effect estimation
Cross-variogram fitting for multivariate kriging workflows
Directional rose diagrams for visualization
Anisotropic distance transformation for directional kriging applications

Spatial Autocorrelation

Measures of spatial structure and clustering for exploratory analysis:

Global Moran's I: Global spatial autocorrelation statistic with asymptotic significance testing
Local Moran's I (LISA): Local indicators of spatial association with significance maps
Getis-Ord Gi*: Local clustering/hotspot identification with asymptotic inference
Spatial Weights Construction:
- Inverse-distance weights (customizable decay exponent)
- Queen/Rook contiguity (vector polygon neighbors)
- K-nearest neighbor weights
- Distance-band (fixed radius) weights

Features:

Multiple testing correction (Bonferroni, FDR-Benjamini-Hochberg)
Island policy configuration (how to handle isolated features)
Standardized/row-normalized weight matrices
Comprehensive diagnostics (weight matrix sparsity, neighbor counts)

Spatial Regression

Spatial inference combining regression with spatial structure:

Spatial Lag (SAR): Include neighboring dependent variable values as predictor
Spatial Error (SEM): Model spatial autocorrelation in residuals
Geographically Weighted Regression (GWR): Local regression with distance-decay kernel weighting
- Kernel functions: Gaussian, Biweight, Tricube
- Adaptive (nearest-k) and fixed-radius bandwidth options
- Local coefficient estimation with uncertainty quantification

Features:

Ordinary least squares (OLS) baseline for comparison
Lagged variable/residual spatial dependence
Local goodness-of-fit diagnostics per regression window
Bandwidth optimization via cross-validation

Point Process Analysis

Statistical tools for point pattern analysis and hypothesis testing:

Ripley's K Function: Second-order intensity analysis (clustering vs. dispersion)
Ripley's L Transform: Centered version of K for easier interpretation
Envelope Testing: Monte Carlo confidence envelopes for point patterns
- Configurable null model (homogeneous Poisson or random shifting)
- Adjustable simulation count for trade-off between precision and speed
Inhomogeneous Intensity: Kernel density estimation for spatially varying rate processes
Quadrat Count Test: Goodness-of-fit test for spatial randomness (Poisson vs. observed counts)

Features:

Distance-based (isotropic) and directional K variants
Bootstrap resampling for residual uncertainty
Support for marked point patterns (optional attribute weighting)

Cross-Validation & Diagnostics

Model validation and diagnostic tools for spatial model assessment:

Leave-One-Out Cross-Validation (LOOCV): Remove each observation, predict it, and assess error
Spatial k-Fold Cross-Validation: Spatially stratified train/test splitting
Prediction Interval Calculation: Confidence bounds from kriging variance
Residual Analysis: Residuals, standardized residuals, and regression diagnostics
Prediction Intervals: Gaussian-based confidence and prediction intervals
Bootstrap Confidence Intervals: Resampling-based uncertainty quantification for kriging

Features:

Root-mean-square error (RMSE) and mean-absolute-error (MAE) computation
Quantile-quantile (Q-Q) plots for normality checking
Spatial residual autocorrelation diagnostics

Installation

Crates.io dependency:

[dependencies]
wbspatialstats = "0.1"

Local workspace/path dependency:

[dependencies]
wbspatialstats = { path = "../wbspatialstats" }

Quick Start

Ordinary Kriging Example

use wbspatialstats::kriging::OrdinaryKriging;
use wbspatialstats::variogram::{VariogramModel, VariogramModelFamily};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Sample data: (x, y) coordinates with z-values
    let primary_coords = vec![
        (0.0, 0.0),
        (1.0, 0.0),
        (0.0, 1.0),
        (1.0, 1.0),
    ];
    let primary_values = vec![10.0, 12.0, 11.0, 13.0];

    // Create variogram model
    let vgm = VariogramModel::exponential(
        1.0,    // range parameter
        1.0,    // sill (total variance)
        0.1,    // nugget
        "z",    // variable name
    );

    // Initialize kriging predictor
    let kriging = OrdinaryKriging::new(
        &vgm,
        &primary_coords,
        &primary_values,
        16,  // neighborhood size
    )?;

    // Predict at a new location
    let prediction = kriging.predict((0.5, 0.5))?;
    println!("Predicted value: {}", prediction.value);
    println!("Kriging variance: {}", prediction.variance);
    println!("Kriging std error: {}", prediction.variance.sqrt());

    Ok(())
}

Empirical Variogram Example

use wbspatialstats::variogram::EmpiricalVariogramBuilder;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Sample data
    let coords = vec![
        (0.0, 0.0),
        (1.0, 0.0),
        (0.0, 1.0),
        (1.0, 1.0),
    ];
    let values = vec![10.0, 12.0, 11.0, 13.0];

    // Build empirical variogram
    let builder = EmpiricalVariogramBuilder::new(&coords, &values);
    let vgm = builder
        .with_max_distance(2.0)
        .with_bin_size(0.5)
        .build()?;

    println!("Empirical variogram: {} lags", vgm.lags.len());
    for lag in &vgm.lags {
        println!("  Distance: {:.2}, Semi-variance: {:.2}", lag.distance, lag.semi_variance);
    }

    Ok(())
}

Spatial Autocorrelation Example

use wbspatialstats::autocorrelation::GlobalMoransI;
use wbspatialstats::weights::SpatialWeightsGraph;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Sample data
    let values = vec![10.0, 12.0, 11.0, 13.0];
    
    // Create spatial weights (inverse distance)
    let weights = SpatialWeightsGraph::inverse_distance(
        &[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
        1.0,  // exponent
    )?;

    // Compute Moran's I
    let morans_i = GlobalMoransI::new(&weights, &values)?;
    println!("Moran's I: {:.4}", morans_i.statistic);
    println!("Expected I (under null): {:.6}", morans_i.expected);
    println!("Variance: {:.8}", morans_i.variance);
    println!("p-value (one-tailed): {:.4}", morans_i.p_value);

    Ok(())
}

Examples

See examples/ directory for additional worked examples:

kriging_workflow.rs — End-to-end kriging pipeline with variogram fitting
spatial_autocorrelation.rs — Computing and interpreting Moran's I and LISA
gwr_local_regression.rs — Geographically weighted regression with local coefficients
point_process.rs — Ripley's K analysis and envelope testing

Known Limitations

Kriging:
- Ordinary kriging assumes stationarity; non-stationary processes may need trend removal first.
- Large datasets (>5000 points) may require local kriging or neighborhood selection for performance.
- Anisotropy support is integrated via directional variogram analysis and anisotropy transformation; arbitrary geometric anisotropy tensors are not yet supported.
- CoKriging supports arbitrary number of auxiliary variables but limited to pairwise cross-variogram relationships.
Variography:
- Nested (multi-structure) variograms are partially supported; some complex nesting patterns may require post-hoc model fitting.
- Cross-variograms are supported but currently limited to primary-auxiliary variable pairs.
- Directional variogram analysis uses azimuthal binning; tolerance-based directional windows are configurable.
Spatial Regression:
- Bayesian spatial models (MCMC, Gibbs sampling) are not yet implemented.
- Spatial quantile regression is not yet implemented.
- Robust M-estimation variants are not yet integrated into regression solvers.
Point Process:
- Marked point patterns are supported but multitype/multivariate patterns have limited statistical support.
- 3D/4D point patterns are not yet implemented.
Numerical:
- Matrix condition number warnings are logged but not enforced; ill-conditioned systems may silently produce unreliable results.
- No preconditioning support for very large matrices.

Architecture & Dependencies

Dependencies

nalgebra — Linear algebra (LU decomposition, matrices, vectors) for kriging solvers and regression
thiserror — Error type derivation and Display formatting
Standard library only for core types and algorithms

Module Structure

variogram/ — Empirical variogram, directional variogram, anisotropy modeling, model families, fitting algorithms, cloud diagnostics
kriging/ — OrdinaryKriging, LocalOrdinaryKriging, SimpleKriging, UniversalKriging, SpaceTimeKriging, OrdinaryCoKriging, and result types
cv/ — Leave-one-out and spatial k-fold cross-validation
weights/ — Spatial weights matrices with multiple construction modes (inverse distance, queen, k-nearest)
autocorrelation/ — Global/local Moran's I, Getis-Ord Gi*, significance testing, multiple testing correction
regression/ — Spatial lag, spatial error, and geographically weighted regression solvers
density_estimation/ — Kernel density estimation for inhomogeneous intensity functions
point_process/ — Ripley's K, envelope testing, quadrat count tests

Error Handling

wbspatialstats::GeostatError is the primary error type, covering:

Invalid variograms (singular matrices, degenerate models)
Kriging solve failures (numerical instability, rank deficiency)
Invalid parameters (negative ranges, inappropriate coordinates)
Numerical instability warnings

Development

Testing

Run the test suite:

cargo test -p wbspatialstats

Run with backtrace on failures:

RUST_BACKTRACE=1 cargo test -p wbspatialstats -- --nocapture

Building Documentation

cargo doc -p wbspatialstats --open

Benchmarking

Benchmarks (if present) can be run with:

cargo bench -p wbspatialstats

License

wbspatialstats is licensed under the MIT license. See LICENSE for details.

wbspatialstats 0.1.1

wbspatialstats

Table of Contents

Mission

The Whitebox Project

Is wbspatialstats Only for Whitebox?

What wbspatialstats Is Not

Design Goals

Current Capabilities

Kriging

Variography

Spatial Autocorrelation

Spatial Regression

Point Process Analysis

Cross-Validation & Diagnostics

Installation

Quick Start

Ordinary Kriging Example

Empirical Variogram Example

Spatial Autocorrelation Example

Examples

Known Limitations

Architecture & Dependencies

Dependencies

Module Structure

Error Handling

Development

Testing

Building Documentation

Benchmarking

License