Computing Focal Window Statistics on Rasters

Compute focal mean, standard deviation, range and entropy on rasters with SciPy: window sizing, nodata-safe kernels, multi-scale stacks and correct tile-edge halos.

A focal statistic summarises a moving window around every pixel — mean, standard deviation, range, entropy. It is the raster counterpart of the spatial lag computed from a weights matrix, and it is often the single cheapest way to give a pixel-wise model the spatial context it otherwise lacks: texture, local contrast, heterogeneity.

The mechanics are scipy.ndimage filters. What makes them non-trivial in production is nodata (which must not be averaged in as zero), window sizing (a modelling decision, not a default), and tile edges (which need a halo). For the vector-side equivalent, see choosing a spatial weights matrix with libpysal; the wider context is in Spatial Lag and Neighborhood Statistics.

Two convolutions give a mean over valid cells only The source array with nodata replaced by zero is convolved with a uniform kernel to give a sum of valid values. A binary mask of the same array is convolved with the same kernel to give a count of valid cells. Dividing the first by the second yields the focal mean over valid cells, unaffected by the number of nodata cells in the window. A single convolution treats nodata as a real zero — two do not values, NaN → 0 convolve with kernel validity mask 0/1 convolve with kernel sum of valid count of valid divide where count > 0 focal mean Windows straddling a gap are averaged over the cells that actually exist, instead of being dragged towards zero by the gap.

Why This Fails in Geospatial ML Pipelines

ndimage.uniform_filter does exactly what it says: it averages every cell in the window, including the ones you meant to exclude. If nodata was left as 0, a window overlapping a cloud gap reports a mean pulled towards zero, and the size of the error depends on how much of the window the gap covers — so the artefact has structure, and the model learns it. If nodata was left as NaN, the filter returns NaN for every window touching the gap, which converts a small hole into a hole the width of the window.

Window size is the second problem, and it is a modelling decision disguised as a parameter. The window has no meaning in pixels; it has meaning in metres. A 5×5 window is 50 m on a Sentinel-2 band, 150 m on Landsat, and 10 m on aerial imagery, so the “same” feature computed in two pipelines is measuring completely different processes. Any window size that appears in code should be derived from a distance in configuration, not typed as an integer.

Third, tile edges. Focal statistics are exactly the kind of neighbourhood operation that breaks when a raster is processed in blocks. Without a halo, the kernel near a tile boundary sees zero-padding or edge replication, so every boundary gets a band of wrong values — a grid of faint lines across the output that is invisible at national zoom and obvious in any tile-level inspection. The halo mechanics are the same as for terrain derivatives; see DEM and terrain derivative features.

Core Principles

  • Express windows in ground units. Derive the pixel window from a metre radius and the cell size.
  • Never average nodata. Use the sum-over-count pattern with a validity mask.
  • Require a minimum valid fraction. A window that is 90% gap should return NaN, not a confident mean over two pixels.
  • Use separable filters where possible. uniform_filter is O(n) in window size; a generic filter is O(n²).
  • Compute the anomaly, not just the mean. value − focal_mean is usually the stronger predictor.
  • Read a halo when tiling. At least window // 2 pixels, discarded before writing.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
from scipy import ndimage

logger = logging.getLogger(__name__)


def window_from_metres(radius_m: float, cell_size_m: float) -> int:
    """Odd pixel window covering a given ground radius."""
    if cell_size_m <= 0:
        raise ValueError("cell_size_m must be positive")
    px = int(round(2 * radius_m / cell_size_m)) + 1
    if px < 3:
        raise ValueError(f"radius {radius_m} m is smaller than one cell ({cell_size_m} m)")
    return px if px % 2 == 1 else px + 1


def focal_stats(arr: np.ndarray, window: int, min_valid_frac: float = 0.5) -> dict:
    """Nodata-safe focal mean, standard deviation and range.

    Args:
        arr: float array with nodata as NaN.
        window: odd pixel window size.
        min_valid_frac: windows with less valid data than this return NaN.

    Returns:
        dict of "mean", "std", "range", "count" arrays.
    """
    if window % 2 == 0:
        raise ValueError("window must be odd so the target cell sits at the centre")

    valid = np.isfinite(arr)
    filled = np.where(valid, arr, 0.0).astype("float64")

    count = ndimage.uniform_filter(valid.astype("float64"), size=window,
                                   mode="nearest") * window ** 2
    total = ndimage.uniform_filter(filled, size=window, mode="nearest") * window ** 2
    total_sq = ndimage.uniform_filter(filled ** 2, size=window, mode="nearest") * window ** 2

    enough = count >= max(1.0, min_valid_frac * window ** 2)
    mean = np.full(arr.shape, np.nan, dtype="float32")
    std = np.full(arr.shape, np.nan, dtype="float32")

    np.divide(total, count, out=mean, where=enough, casting="unsafe")
    variance = np.zeros_like(mean, dtype="float64")
    np.divide(total_sq, count, out=variance, where=enough)
    variance -= np.where(enough, mean.astype("float64") ** 2, 0.0)
    np.sqrt(np.maximum(variance, 0.0), out=variance)
    std[enough] = variance[enough].astype("float32")

    # Range needs order statistics; NaN must not participate.
    big = np.where(valid, arr, -np.inf)
    small = np.where(valid, arr, np.inf)
    hi = ndimage.maximum_filter(big, size=window, mode="nearest")
    lo = ndimage.minimum_filter(small, size=window, mode="nearest")
    rng = np.where(enough & np.isfinite(hi) & np.isfinite(lo), hi - lo, np.nan)

    logger.info("focal window %d px: %.1f%% of cells had enough valid data",
                window, 100 * enough.mean())
    return {"mean": mean, "std": std, "range": rng.astype("float32"),
            "count": count.astype("float32")}


def focal_anomaly(arr: np.ndarray, window: int) -> np.ndarray:
    """value minus its neighbourhood mean — local contrast, scale-free."""
    return (arr - focal_stats(arr, window)["mean"]).astype("float32")

Step-by-Step Walkthrough

Step 1 — Derive the window from a distance, not a habit

import rasterio

with rasterio.open("ndvi.tif") as src:
    cell = abs(src.transform.a)
    ndvi = src.read(1, masked=True).filled(np.nan).astype("float32")
    profile = src.profile.copy()

W_LOCAL = window_from_metres(radius_m=50, cell_size_m=cell)     # field-scale texture
W_REGION = window_from_metres(radius_m=500, cell_size_m=cell)   # landscape context
print(f"cell {cell} m -> windows {W_LOCAL} px and {W_REGION} px")

Step 2 — Compute a multi-scale stack

local = focal_stats(ndvi, W_LOCAL)
region = focal_stats(ndvi, W_REGION)

features = {
    "ndvi": ndvi,
    "ndvi_mean_50m": local["mean"],
    "ndvi_std_50m": local["std"],
    "ndvi_anom_50m": ndvi - local["mean"],
    "ndvi_mean_500m": region["mean"],
    "ndvi_anom_500m": ndvi - region["mean"],
}
for name, arr in features.items():
    print(f"{name:>16}: {np.nanmin(arr):+.3f} to {np.nanmax(arr):+.3f}")

The anomaly bands are usually the most predictive of the set: they encode how a pixel differs from its surroundings, which is scale-free and robust to sensor-level offsets in a way the raw value is not.

Step 3 — Check that the statistics behave near gaps

gap = ~np.isfinite(ndvi)
edge = ndimage.binary_dilation(gap, iterations=W_LOCAL // 2) & ~gap
print(f"cells adjacent to a gap: {edge.sum():,}")
print(f"  focal mean is NaN for {np.isnan(local['mean'][edge]).mean():.1%} of them")
print(f"  median valid count in their windows: {np.median(local['count'][edge]):.0f} "
      f"of {W_LOCAL ** 2}")

If almost none of them are NaN, min_valid_frac is too permissive and windows dominated by gaps are reporting confident means over a handful of pixels.

Step 4 — Tile with a halo

HALO = max(W_LOCAL, W_REGION) // 2 + 1

with rasterio.open("ndvi.tif") as src:
    out_profile = src.profile.copy()
    out_profile.update(count=len(features), dtype="float32", nodata=np.nan,
                       compress="deflate", tiled=True)
    with rasterio.open("focal_stack.tif", "w", **out_profile) as dst:
        for _, win in src.block_windows(1):
            padded = win_with_halo(win, HALO, src.width, src.height)
            tile = src.read(1, window=padded, masked=True).filled(np.nan).astype("float32")
            stats = focal_stats(tile, W_LOCAL)
            trimmed = crop_halo(stats["mean"], win, padded)
            dst.write(trimmed, 2, window=win)
Window size decides which process the feature measures A single target pixel is shown twice over the same landscape. With a small window the statistic responds to variation inside one field. With a large window it spans several fields and a woodland, so it responds to landscape composition instead. The two are different features, not two settings of one feature. A 5×5 window is 50 m on Sentinel-2 and 150 m on Landsat radius 50 m — within-field texture stays inside one parcel → measures crop uniformity and stress patches radius 500 m — landscape context wood spans parcels and woodland → measures landscape composition

Verification

import numpy as np


def test_focal_mean_on_a_constant_field():
    arr = np.full((40, 40), 7.5, dtype="float32")
    assert np.allclose(focal_stats(arr, 5)["mean"], 7.5)
    assert np.allclose(focal_stats(arr, 5)["std"], 0.0, atol=1e-5)


def test_nodata_is_not_averaged_as_zero():
    """A window that is half nodata must return the mean of the valid half."""
    arr = np.full((9, 9), 4.0, dtype="float32")
    arr[:, :4] = np.nan
    out = focal_stats(arr, 3, min_valid_frac=0.1)["mean"]
    assert np.allclose(out[4, 5:], 4.0), "valid cells were diluted by nodata"


def test_sparse_windows_return_nan():
    arr = np.full((11, 11), 2.0, dtype="float32")
    arr[1:, :] = np.nan                       # only the top row is valid
    out = focal_stats(arr, 5, min_valid_frac=0.5)["mean"]
    assert np.isnan(out[5, 5])


def test_window_from_metres_is_odd_and_covers_the_radius():
    w = window_from_metres(radius_m=50, cell_size_m=10)
    assert w % 2 == 1 and w >= 11


def test_halo_removes_tile_seams():
    """Tiled computation with a halo must equal the whole-array computation."""
    rng = np.random.default_rng(0)
    arr = rng.normal(size=(256, 256)).astype("float32")
    whole = focal_stats(arr, 9)["mean"]
    tiled = tiled_focal(arr, window=9, tile=64, halo=5)
    assert np.allclose(whole[10:-10, 10:-10], tiled[10:-10, 10:-10], atol=1e-5)

The last test is the one worth keeping in CI. Halo bugs are easy to introduce during a refactor and produce artefacts that look like real spatial structure.

Without a halo, the tiling shows up in the feature Two focal output rasters are shown. The first, computed tile by tile without a halo, has visible lines of incorrect values along every internal tile boundary. The second, computed with a halo, is continuous with no trace of the tiling. no halo halo = window // 2 + 1 a grid of wrong values, one per boundary tiling is invisible in the output The artefact is faint enough to survive review and strong enough for a model to learn.

FAQ

How do I compute a focal mean that ignores nodata?

Convolve twice: once over the array with nodata set to zero, once over a binary validity mask. Dividing the first by the second gives the mean over valid cells. A single convolution treats nodata as a genuine zero and drags every nearby window down.

What window size should I use?

Whichever matches the process. Convert to ground units first — a 5×5 window is 50 m at 10 m resolution — and pick from domain knowledge. When the scale is unknown, compute two or three windows as separate bands and let importance decide, as in model explainability for spatial predictions.

Why is my focal output wrong at tile edges?

Because the kernel sees padding rather than the neighbouring tile. Read each tile with a halo of at least half the window width, compute over the padded array, and discard the halo before writing.

Part of: Spatial Lag and Neighborhood Statistics Part of: Spatial Feature Engineering for Machine Learning