Creating Distance Matrices for Spatial Features

Create spatial distance matrices for machine learning features in Python: compute pairwise proximity metrics from vector geometries using scipy and GeoPandas in production pipelines.

Creating distance matrices for spatial features transforms raw vector geometries into structured proximity metrics that machine learning models consume directly. Instead of feeding absolute coordinates—which introduce scale bias and ignore spatial autocorrelation—pairwise distance matrices encode relative neighborhood relationships. This operation sits at the foundation of Spatial Feature Engineering for Machine Learning, where proximity features consistently outperform raw lat/lon inputs in gradient boosting and neural architectures.

The production-standard workflow extracts planar coordinates from validated geometries, enforces a metric coordinate reference system (CRS), and computes pairwise distances using optimized numerical routines. For MLOps environments, the computation must be deterministic, memory-aware, and pipeline-compatible.

Prerequisites & CRS Validation

Accurate distance computation requires planar coordinates. Geographic systems (e.g., EPSG:4326) measure angles, not meters, making Euclidean or Manhattan distances mathematically invalid. Always project to a metric CRS before calculation. Regional UTM zones minimize distortion for localized datasets, while EPSG:3857 (Web Mercator) serves as a global fallback despite high-latitude area inflation. Consult the GeoP CRS transformation guide for projection selection and distortion trade-offs.

Before computing distances, filter invalid geometries to prevent silent NaN propagation or topology errors. Self-intersecting polygons, empty shapes, and null geometries will break downstream matrix operations. A strict validation mask ensures reproducible outputs across training and inference environments.

Production Implementation

The following function generates a symmetric distance matrix from a GeoDataFrame, enforces metric CRS, filters invalid geometries, and supports optional sparsification for memory-constrained workloads.

import numpy as np
import geopandas as gpd
from scipy.spatial.distance import pdist, squareform
from scipy.sparse import csr_matrix
from sklearn.metrics import pairwise_distances
import warnings
from typing import Union

def create_spatial_distance_matrix(
    gdf: gpd.GeoDataFrame,
    metric: str = "euclidean",
    sparse_threshold: float = 0.0,
    target_crs: str = "EPSG:3857"
) -> Union[np.ndarray, csr_matrix]:
    """
    Compute pairwise distance matrix for spatial features.
    Optimized for MLOps pipeline integration, drift detection, and memory efficiency.

    Parameters
    ----------
    gdf : GeoDataFrame
        Input vector data with valid geometries.
    metric : str
        Distance metric passed to scipy/sklearn (e.g., 'euclidean', 'manhattan', 'cosine').
    sparse_threshold : float
        Distances > threshold are zeroed and converted to sparse CSR format.
        Use 0.0 to return a dense numpy array.
    target_crs : str
        Metric CRS for projection (default: EPSG:3857).

    Returns
    -------
    np.ndarray or scipy.sparse.csr_matrix
        Symmetric distance matrix aligned with input GDF index.
    """
    if gdf.empty:
        raise ValueError("Input GeoDataFrame is empty.")

    # Enforce valid geometries
    valid_mask = gdf.geometry.is_valid & gdf.geometry.notna()
    if not valid_mask.all():
        invalid_count = (~valid_mask).sum()
        warnings.warn(f"Dropping {invalid_count} invalid or null geometries.")
        gdf = gdf[valid_mask].copy()

    # CRS validation and transformation
    if gdf.crs is None:
        raise RuntimeError("GeoDataFrame lacks CRS. Assign before distance computation.")
    if gdf.crs.is_geographic:
        gdf = gdf.to_crs(target_crs)
        warnings.warn(f"Transformed from geographic CRS to {target_crs} for metric distances.")

    # Extract planar coordinates (centroids for polygons/lines)
    coords = np.column_stack((gdf.geometry.centroid.x, gdf.geometry.centroid.y))
    n_samples = coords.shape[0]

    # Compute distances
    if metric == "euclidean" and n_samples > 1:
        # scipy.pdist is faster for pure Euclidean
        condensed = pdist(coords, metric="euclidean")
        dist_matrix = squareform(condensed)
    else:
        # sklearn handles arbitrary metrics and single-sample edge cases
        dist_matrix = pairwise_distances(coords, metric=metric)

    # Optional sparsification for memory-constrained environments
    if sparse_threshold > 0.0:
        mask = dist_matrix > sparse_threshold
        dist_matrix[mask] = 0.0
        dist_matrix = csr_matrix(dist_matrix)

    return dist_matrix

MLOps Pipeline Integration

In automated geospatial workflows, distance matrices must behave deterministically across training, validation, and inference stages. Wrap the function in a custom scikit-learn transformer to guarantee consistent index alignment and prevent data leakage during cross-validation.

For drift detection, log baseline distance distributions (mean, median, 95th percentile) to MLflow or a feature store. Production shifts in spatial sampling density—such as new sensor deployments or regional boundary expansions—will skew distance distributions. Monitoring these shifts prevents silent model degradation when spatial autocorrelation patterns change.

When paired with Vector Proximity and Buffer Generation, distance matrices enable both global neighborhood aggregation and localized feature scaling. Proximity matrices feed directly into graph convolutional networks, spatial lag features, and kernel-based similarity metrics.

Memory Optimization & Sparse Conversion

Full distance matrices scale at $O(N^2)$ memory. A dataset of 50,000 points requires ~20 GB as a dense float64 array. For tree-based models or large-scale inference, apply sparse_threshold to retain only meaningful neighborhood connections. The resulting scipy.sparse.csr_matrix reduces memory footprint by 60–90% while preserving model performance, as most algorithms ignore zero-weight edges.

If your pipeline requires strict memory caps, consider:

  • Chunked computation: Process subsets and merge via block-sparse formats.
  • Nearest-neighbor pruning: Use scipy.spatial.KDTree to limit matrix generation to $k$-nearest neighbors.
  • Float32 casting: Halve memory usage with negligible accuracy loss for most spatial ML tasks.

Performance & Scaling Strategies

The implementation above defaults to scipy.spatial.distance.pdist for Euclidean calculations, which leverages vectorized C routines and outperforms generic pairwise loops. For non-Euclidean metrics, sklearn.metrics.pairwise_distances provides a unified interface with optimized BLAS backends. Refer to the SciPy spatial.distance documentation for supported metrics and numerical stability notes.

When scaling to millions of features:

  1. Pre-filter spatial bounds using bounding box intersections before matrix generation.
  2. Parallelize chunking with joblib or Dask to distribute $O(N^2)$ workloads across cores.
  3. Cache transformed coordinates to avoid repeated CRS projections during hyperparameter tuning.

Conclusion

Creating distance matrices for spatial features bridges raw vector data and machine learning-ready proximity representations. By enforcing metric CRS, validating geometries, and selecting appropriate dense or sparse outputs, teams can deploy spatial proximity features reliably across training and inference pipelines. Properly integrated distance matrices capture spatial autocorrelation, reduce coordinate-system bias, and provide a scalable foundation for advanced geospatial modeling.