Geostatistical Modeling Spatial Statistics with Python

Reducing Memory Bottlenecks in Geospatial Workflows

TL;DR: Cast DataFrame columns to float32 with df[col].astype(np.float32), read rasters via rasterio.windows.from_bounds, replace dense distance arrays with scipy.sparse CSR matrices using libpysal’s sparse constructors, and wrap Dask pipelines so .compute() fires only at final aggregation. These four changes together cut peak RAM usage by 60–80% on typical nationwide vector or 1-m DEM workflows.

Why This Matters

High-resolution DEMs, dense environmental sensor networks, and nationwide parcel layers routinely exceed available RAM when processed with default library configurations. MemoryError exceptions during variogram estimation, spatial lag regression, or Monte Carlo uncertainty propagation are not hardware problems — they are pipeline design problems. This page drills into the concrete techniques covered in the Memory-Efficient Processing cluster, translating each principle into copy-pasteable code you can drop into any geostatistical workflow built on the Python Workflows for Spatial Modeling & Regression stack.

The bottlenecks follow a predictable pattern: naive read() calls load entire files into contiguous memory; default float64 dtypes double every numeric column’s footprint; dense spatial weight matrices scale quadratically with observation count; and eager evaluation materialises every intermediate DataFrame before any reduction has occurred. Fixing all four in sequence usually brings a crashing pipeline back within 8–16 GB on commodity hardware.

Environment and Version Pinning

bash
pip install \
  geopandas==0.14.4 \
  pyogrio==0.7.2 \
  rasterio==1.3.10 \
  rioxarray==0.15.5 \
  dask==2024.3.1 \
  dask-geopandas==0.3.1 \
  xarray==2024.2.0 \
  libpysal==4.10.0 \
  scipy==1.13.0 \
  numpy==1.26.4 \
  memory-profiler==0.61.0 \
  psutil==5.9.8
python
import numpy as np
import geopandas as gpd
import rasterio
from rasterio.windows import from_bounds
import dask_geopandas as dgpd
from libpysal.weights import KNN
from scipy.sparse import csr_matrix
import tracemalloc
import psutil, os

Enable the pyogrio engine globally to reduce serialisation overhead on every vector read:

python
gpd.options.io_engine = "pyogrio"

Step-by-Step Implementation

Step 1 — Baseline Memory Profiling

Before optimising anything, measure. gdf.memory_usage(deep=True) reports per-column byte counts including object-column serialisation costs. tracemalloc tracks Python-heap allocations; psutil tracks the OS-level RSS delta that includes C extensions and GDAL buffers.

python
def report_memory(gdf: gpd.GeoDataFrame, label: str = "") -> None:
    """Print per-column and total memory for a GeoDataFrame."""
    usage = gdf.memory_usage(deep=True)
    total_mb = usage.sum() / 1024**2
    print(f"[{label}] Total: {total_mb:.1f} MB")
    for col, bytes_ in usage.items():
        print(f"  {col}: {bytes_ / 1024**2:.2f} MB  ({gdf[col].dtype if col != 'Index' else 'index'})")

Run this before and after every transformation stage. If the delta does not match your expectation, you have a hidden copy or implicit type promotion.

Step 2 — Numeric Downcasting

Spatial DataFrames default to float64 and int64 even when attribute ranges and coordinate precision do not require it. Downcasting numeric columns to float32 or int32 typically halves memory usage without degrading geostatistical accuracy.

python
def downcast_numeric(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """
    Downcast float64 → float32 and int64 → int32 for all numeric columns.
    Coordinate precision rarely exceeds 10⁻⁵ degrees (~1 m), so float32
    is sufficient for regional analyses.
    """
    gdf = gdf.copy()
    for col in gdf.select_dtypes(include=["float64"]).columns:
        gdf[col] = gdf[col].astype(np.float32)
    for col in gdf.select_dtypes(include=["int64"]).columns:
        gdf[col] = gdf[col].astype(np.int32)
    return gdf

Verify statistical parity after downcasting:

python
original_mean = gdf["elevation"].mean()
gdf_dc = downcast_numeric(gdf)
assert np.allclose(original_mean, gdf_dc["elevation"].mean(), rtol=1e-4), \
    "Downcast changed summary statistics beyond tolerance"

For geometry vertices, shapely.set_precision(gdf.geometry, grid_size=1e-5) snaps coordinates to a fixed grid before any spatial operation, reducing both memory and topology edge cases during subsequent joins.

Step 3 — Windowed Raster Reading

Full-tile raster loads are unsustainable for modern high-resolution datasets. Fetch only the bounding box that intersects your analysis points using rasterio’s windowed reading API.

python
def sample_raster_windowed(
    raster_path: str,
    points_gdf: gpd.GeoDataFrame,
    chunk_size: int = 5000
) -> gpd.GeoDataFrame:
    """
    Memory-efficient point sampling using windowed raster reads.
    Avoids OOM by processing points in spatial chunks and reading
    only pixels that intersect each chunk's bounding box.

    Parameters
    ----------
    raster_path  : Path to a GeoTIFF (any resolution or size).
    points_gdf   : GeoDataFrame of sampling locations in the raster CRS.
    chunk_size   : Points per iteration; 5,000–20,000 keeps peak RAM < 500 MB
                   for typical 1-m DEMs at 10,000 km² extent.

    Returns
    -------
    GeoDataFrame with an 'extracted_value' column appended.
    """
    # Pre-allocate results array — avoids dynamic resizing overhead
    extracted = np.full(len(points_gdf), np.nan, dtype=np.float32)

    with rasterio.open(raster_path) as src:
        raster_bounds = src.bounds

        for i in range(0, len(points_gdf), chunk_size):
            chunk = points_gdf.iloc[i : i + chunk_size]
            bounds = chunk.total_bounds  # (minx, miny, maxx, maxy)

            # Compute tight window for the current spatial chunk
            window = from_bounds(*bounds, transform=src.transform)

            # Clip to raster extent to guard against edge-of-study-area chunks
            col_off = max(int(window.col_off), 0)
            row_off = max(int(window.row_off), 0)
            col_end = min(int(window.col_off + window.width), src.width)
            row_end = min(int(window.row_off + window.height), src.height)

            if col_end <= col_off or row_end <= row_off:
                continue  # Chunk lies entirely outside raster extent

            clipped = rasterio.windows.Window(
                col_off=col_off, row_off=row_off,
                width=col_end - col_off,
                height=row_end - row_off,
            )

            # Read only necessary pixels; masked=True preserves NaN for NoData
            tile = src.read(1, window=clipped, masked=True)

            # Map geographic coordinates → array row/col indices
            rows, cols = rasterio.transform.rowcol(
                src.window_transform(clipped),
                chunk.geometry.x.values,
                chunk.geometry.y.values,
            )

            # Discard any indices that fall outside the clipped tile
            valid = (
                (rows >= 0) & (rows < clipped.height) &
                (cols >= 0) & (cols < clipped.width)
            )
            extracted[i : i + chunk_size][valid] = tile[rows[valid], cols[valid]]

    result = points_gdf.copy()
    result["extracted_value"] = extracted
    return result

Chunk size is tunable based on available RAM. The table below gives typical working-memory figures for a 0.5-m DEM:

Chunk size Pixels read per iteration Approx. peak RAM
1,000 pts ~50,000 px ~40 MB
5,000 pts ~250,000 px ~180 MB
20,000 pts ~1,000,000 px ~700 MB

Step 4 — Sparse Spatial Weights

Geostatistical models require spatial weight matrices or distance structures that scale quadratically with observation count. Building dense adjacency arrays for datasets exceeding 50,000 points is not feasible. Replace them with scipy.sparse CSR matrices and libpysal’s sparse weight constructors.

python
from libpysal.weights import KNN, DistanceBand
from scipy.sparse import csr_matrix

def build_sparse_knn_weights(
    gdf: gpd.GeoDataFrame,
    k: int = 8,
) -> csr_matrix:
    """
    Construct K-nearest-neighbour spatial weights as a CSR sparse matrix.
    Memory scales with k * n, not n², making this viable for 500k+ observations.
    """
    # libpysal's KNN uses a spatial index internally — no O(n²) pairwise step
    w = KNN.from_dataframe(gdf, k=k)
    return csr_matrix(w.sparse)

Avoid calling .toarray() unless a legacy algorithm explicitly requires a dense matrix. Operations like matrix-vector multiplication (w_sparse @ y) work directly on CSR format and are equally fast via BLAS routines.

Step 5 — Lazy Evaluation with Dask

Defer heavy computation until aggregation or model fitting. Dask builds a task graph instead of materialising intermediate arrays. This pattern aligns with the deferred-execution model described in cross-validation strategies, where fold assignment and feature extraction must not trigger premature materialisation.

python
import dask_geopandas as dgpd

def lazy_spatial_pipeline(path: str, target_gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """
    Build a Dask-backed pipeline that reads, filters, and spatially joins a large
    vector dataset without materialising intermediate partitions.
    """
    # Read into Dask partitions — no data loaded yet
    ddf = dgpd.read_file(path, npartitions=16)

    # .assign() and boolean indexing build the task graph lazily
    ddf = ddf.assign(area_km2=ddf.geometry.area / 1e6)
    ddf = ddf[ddf["area_km2"] > 0.01]

    # .compute() fires only here, after all reductions are defined
    result = ddf.sjoin(target_gdf, how="inner").compute()
    return result

Use .persist() only when the same intermediate result is reused across multiple downstream operations. Calling .compute() prematurely materialises the full graph and negates the memory benefit.

Interpreting the Output

After running sample_raster_windowed, inspect the extracted_value column for NaN patterns:

python
nan_rate = result["extracted_value"].isna().mean()
print(f"NaN rate: {nan_rate:.1%}")
# Expected: < 1% for points fully within the raster extent
# High NaN rate: check CRS alignment — points and raster must share the same CRS

A NaN rate above 5% usually indicates a CRS mismatch between the point GeoDataFrame and the raster. Always reproject points to match the raster’s native CRS before calling sample_raster_windowed — never reproject the raster, which would trigger a full re-read and resample.

After downcasting, verify the column dtype and memory reduction:

python
report_memory(gdf, label="before")
gdf = downcast_numeric(gdf)
report_memory(gdf, label="after")
# Expect 40–55% reduction in total DataFrame memory for typical attribute tables

Critical Best Practices

Avoid Hidden Copies from Chained Assignments

Each .copy() duplicates geometry and attribute arrays in memory. Use .loc[] or boolean masking for in-place filtering. When a copy is unavoidable — such as when returning a modified GeoDataFrame from a function — make it explicit and document the intent.

Enforce CRS Alignment Before Any Spatial Operation

CRS mismatches are the most common source of unexpected NaN floods and inflated memory use from failed spatial indexes. Always call assert gdf.crs == target_crs at pipeline entry points. The GeoPandas data preparation workflow covers systematic CRS harmonisation before chunking begins.

Bound Spatial Joins with Explicit Distance Thresholds

gpd.sjoin() on large datasets creates Cartesian products when geometries overlap extensively. Always filter bounding boxes first or use sjoin_nearest() with an explicit max_distance parameter to prevent unbounded join expansion.

Control Dask’s .compute() Boundaries

Calling .compute() inside a loop re-materialises the task graph on every iteration. Build the complete transformation chain outside the loop, call .compute() once, then iterate over the result. For iterative algorithms like EM-based spatial regression, use .persist() to keep intermediate results in distributed memory across iterations.

Monitor Implicit Type Promotion

Integer arithmetic in NumPy promotes silently: int32 + float64 → float64. This can undo all your downcasting gains in a single arithmetic expression. Use np.result_type(a.dtype, b.dtype) to preview the promotion before combining columns, and explicitly cast the result back with .astype(np.float32).

Troubleshooting

Symptom Likely cause Fix
MemoryError during variogram estimation All observation pairs loaded into dense lag matrix Sample to ≤ 10,000 observations or use scikit-gstat’s n_lags + maxlag parameters to reduce pair count
NaN rate > 10% after windowed sampling CRS mismatch between points and raster Reproject points_gdf to src.crs before calling sample_raster_windowed
Memory usage identical before and after downcasting Column contains Python object dtype, not float64 Convert with pd.to_numeric(df[col], errors='coerce') first, then downcast
dask.compute() produces MemoryError Partition size too large for available RAM Increase npartitions so each partition fits in ~512 MB; check with ddf.memory_usage(deep=True).compute()
Sparse weight construction hangs for > 500k observations KNN fallback to brute-force distance Ensure libpysal ≥ 4.9 uses scikit-learn’s ball-tree; set silence_warnings=True and check w.n_components
Peak RSS spikes despite low Python heap GDAL/rasterio C-level buffer not freed Call src.close() explicitly or use the context manager; avoid keeping rasterio.DatasetReader objects in long-lived scopes

Next Steps

For the broader context of fitting spatial regression models under memory constraints — including sparse lag matrices for SAR and SEM models and chunked GMM estimation — see the Memory-Efficient Processing cluster. For workflows that feed the sampled raster values into spatial cross-validation folds without leakage, see Spatial K-Fold Cross-Validation Setup.

Memory-efficient geospatial pipeline Four sequential pipeline stages: (1) Raw Data In with downcast types, (2) Windowed Raster Read limited to intersecting pixels, (3) Sparse Weight Matrix using CSR format, (4) Lazy Model Graph deferred until compute(). Arrows connect each stage. Peak RAM decreases from left to right. Raw Data In Windowed Read Sparse Weights Lazy Graph float64 / int64 dense columns full-tile loads → downcast float32 / int32 from_bounds() chunk_size=5000 pixels only reads only intersecting px CSR format KNN k=8 O(n·k) memory vs O(n²) dense array task graph no materialise until .compute() deferred aggregation Peak RAM reduction across pipeline stages baseline –25% –50% –75%

FAQ

Why does GeoPandas default to float64 even for elevation data? NumPy’s default dtype promotion rules upcast mixed-type inputs to float64. GeoPandas inherits this behaviour. Explicitly cast with astype(np.float32) after reading to reclaim half the column memory.

When should I use windowed rasterio reads versus rioxarray with Dask chunks? Use windowed rasterio reads for point-sampling workflows where you need exact pixel values at arbitrary locations. Prefer rioxarray with Dask chunks when you need whole-grid computations (NDVI, convolutions, zonal statistics) that benefit from parallelism across bands or tiles.

Does downcasting to float32 affect kriging accuracy? For the vast majority of regional analyses, no. float32 provides roughly 7 decimal digits of precision, corresponding to sub-centimetre coordinate accuracy at global scales. Verify with np.allclose(original, downcast, rtol=1e-4) on summary statistics after conversion.

How large can sparse spatial weight matrices get before scipy.sparse runs out of memory? A CSR matrix for 500,000 observations with k=8 neighbours stores roughly 4 million non-zero entries — about 32 MB for float32 indices and data arrays combined, well within 16 GB RAM. Dense equivalents of the same size would require 2 TB.

Related:

← Back to Memory-Efficient Processing for Spatial Statistics