Applying Robust Scaling to Skewed Spectral Bands

Use scikit-learn RobustScaler on skewed, outlier-heavy spectral bands: median/IQR scaling, when to log-transform first, per-band fitting, and leakage-free persistence.

TL;DR Answer

Spectral bands are heavy-tailed. Cloud edges, sensor saturation, water bodies, and shadow pixels push a handful of values thousands of units away from the mode, and StandardScaler — which subtracts the mean and divides by the standard deviation — lets those few outliers dominate both statistics. Use sklearn.preprocessing.RobustScaler, which centres on the median and scales by the interquartile range (IQR), so the transform is anchored to the bulk of the distribution rather than its tail. Fit the scaler only on your training pixels, persist it with joblib, and apply the identical object to every new tile at inference. This is a core step in Feature Scaling for Geospatial Inputs.

from sklearn.preprocessing import RobustScaler
import numpy as np

# X_train: (n_pixels, n_bands) array of training-fold pixels only
scaler = RobustScaler(quantile_range=(25.0, 75.0)).fit(X_train)
X_train_scaled = scaler.transform(X_train)   # median ~0, IQR ~1 per band

Part of: Feature Scaling for Geospatial Inputs

Why StandardScaler Breaks on Spectral Data

StandardScaler assumes each feature is roughly Gaussian. Spectral bands rarely are. A Sentinel-1 SAR backscatter band in linear power units is strongly right-skewed — most land pixels sit near a low value while double-bounce structures and metallic targets stretch the distribution far to the right. A thermal band contains cold cloud tops and hot bare soil in the same scene. An optical reflectance band that was supposed to be masked still carries residual cloud and shadow pixels at the extremes. In every case the distribution has a heavy tail, and often a physical floor or ceiling from sensor saturation.

The mean and standard deviation are not robust statistics: a single pixel at 10,000 in a band whose bulk sits near 200 drags the mean up and inflates the standard deviation. After StandardScaler, the informative variation in the central 90% of pixels is compressed into a tiny slice of the output range while the outliers occupy everything else. Gradient-based models and distance-based models both suffer — the scaled feature carries almost no usable gradient across the pixels you actually care about.

The damage propagates downstream:

  • Distance and kernel models: k-nearest-neighbours, SVMs, and any RBF kernel weight features by scaled magnitude. An outlier-inflated scale silently down-weights a genuinely predictive band.
  • Neural network convergence: unnormalised heavy tails produce large, erratic activations early in training, forcing tiny learning rates or triggering divergence.
  • Feature comparability: when you concatenate bands with band-math outputs from Raster Band Math and Index Calculation, one poorly scaled band dominates the others in the combined matrix.

RobustScaler fixes this by construction. The median ignores the tail entirely, and the IQR is the spread of the middle 50% of pixels — both are unaffected by how extreme the extremes are. The output is centred near zero with the informative bulk spread across roughly ±1.

Core Principles for Scaling Spectral Bands

  • Transform before scaling when the band is multiplicatively skewed. RobustScaler recentres and rescales but does not reshape — a long right tail stays long. For strictly positive, heavy-tailed bands (SAR in linear power, radiance, aerosol depth) apply np.log1p or a PowerTransformer first, then scale the transformed values.
  • Scale per band, never globally. A single RobustScaler fitted on a (pixels, bands) matrix already computes independent median and IQR per column. Never flatten all bands into one column — a kelvin thermal band and a 0–1 reflectance band do not share a centre.
  • Fit on training pixels only. The median and IQR are learned parameters. Fitting on the full dataset, or refitting per tile at inference, leaks distribution information and inflates validation scores.
  • Mask nodata and saturated pixels before fitting. Fill values (-9999, 0 over nodata) are not real observations. Exclude them from the fit or they corrupt the very quantiles you are relying on.
  • Persist the fitted object. Save the scaler with joblib alongside the model. Inference must apply the exact training-time median and IQR to every incoming tile.
  • Widen the quantile range for very heavy tails. The default IQR is the 25th–75th percentile. For extremely peaked bands, a wider quantile_range=(10.0, 90.0) produces a more stable scale.

Production-Ready Code

The function below fits a per-band RobustScaler on a training array, optionally applies a log1p transform to named skewed bands first, masks nodata, and persists both the transform choice and the scaler to disk. It mirrors the fit-on-train-only discipline used throughout Spatial Cross-Validation Strategies — the scaler is a model parameter and must never see validation data during fitting.

import logging
from pathlib import Path

import joblib
import numpy as np
from sklearn.preprocessing import RobustScaler

logger = logging.getLogger(__name__)


def fit_spectral_scaler(
    X_train: np.ndarray,
    log_bands: tuple[int, ...] = (),
    nodata: float | None = None,
    quantile_range: tuple[float, float] = (25.0, 75.0),
    out_path: str | Path = "spectral_scaler.joblib",
) -> dict:
    """Fit a per-band RobustScaler on training pixels and persist it.

    Parameters
    ----------
    X_train : np.ndarray of shape (n_pixels, n_bands)
        Training-fold pixels ONLY. Never pass validation or inference tiles.
    log_bands : tuple of int
        Column indices of strictly positive, right-skewed bands to log1p first.
    nodata : float or None
        Sentinel value to mask out before computing quantiles.
    quantile_range : tuple of float
        Lower/upper percentile bounds for the IQR. Widen for heavy tails.
    out_path : str or Path
        Destination for the persisted transform bundle.

    Returns
    -------
    dict
        The persisted bundle: {"scaler", "log_bands", "nodata", "quantile_range"}.
    """
    if X_train.ndim != 2:
        raise ValueError("X_train must be 2-D: (n_pixels, n_bands).")
    if any(b < 0 or b >= X_train.shape[1] for b in log_bands):
        raise ValueError("log_bands contains an out-of-range column index.")

    X = X_train.astype(np.float64, copy=True)

    # Mask nodata by turning sentinels into NaN so we can drop them per column.
    if nodata is not None:
        X[X == nodata] = np.nan

    # Compress skewed bands before scaling. log1p requires values > -1.
    for b in log_bands:
        col = X[:, b]
        if np.nanmin(col) <= -1:
            raise ValueError(
                f"Band {b} has values <= -1; log1p is undefined. "
                "Offset the band or drop it from log_bands."
            )
        X[:, b] = np.log1p(col)

    # Drop rows with any NaN so the scaler fits only on valid observations.
    valid = ~np.isnan(X).any(axis=1)
    n_dropped = int((~valid).sum())
    if n_dropped:
        logger.info("Dropping %d pixel rows with nodata before fit.", n_dropped)
    X_valid = X[valid]
    if X_valid.shape[0] == 0:
        raise ValueError("No valid pixels remain after nodata masking.")

    scaler = RobustScaler(quantile_range=quantile_range).fit(X_valid)

    bundle = {
        "scaler": scaler,
        "log_bands": log_bands,
        "nodata": nodata,
        "quantile_range": quantile_range,
    }
    joblib.dump(bundle, out_path)
    logger.info("Persisted spectral scaler to %s", out_path)
    return bundle


def apply_spectral_scaler(X: np.ndarray, bundle: dict) -> np.ndarray:
    """Apply a persisted spectral scaler to a new tile.

    Reproduces the training-time transform exactly: same log bands,
    same nodata handling, same median/IQR. Nodata pixels are returned as NaN.
    """
    scaler = bundle["scaler"]
    log_bands = bundle["log_bands"]
    nodata = bundle["nodata"]

    Xt = X.astype(np.float64, copy=True)
    if nodata is not None:
        Xt[Xt == nodata] = np.nan
    for b in log_bands:
        Xt[:, b] = np.log1p(Xt[:, b])

    return scaler.transform(Xt)

Step-by-Step Walkthrough

1. Read the training raster and flatten to a pixel matrix

import numpy as np
import rioxarray

# Multi-band training scene: SAR, thermal, and reflectance stacked
da = rioxarray.open_rasterio("train_scene.tif")   # (band, y, x)
bands, height, width = da.shape

# Reshape to (n_pixels, n_bands); RobustScaler fits per column
X_train = da.values.reshape(bands, -1).T.astype(np.float64)
print(f"Pixel matrix: {X_train.shape}")   # (height*width, bands)

2. Fit the scaler on training pixels only

Band 0 here is SAR backscatter in linear power units — strongly right-skewed — so it is log-transformed before scaling. The nodata sentinel -9999 is masked out of the fit.

bundle = fit_spectral_scaler(
    X_train,
    log_bands=(0,),
    nodata=-9999.0,
    quantile_range=(25.0, 75.0),
    out_path="models/spectral_scaler.joblib",
)

3. Confirm the per-band statistics

scaler = bundle["scaler"]
for b, (center, scale) in enumerate(zip(scaler.center_, scaler.scale_)):
    print(f"Band {b}: median={center:.3f}, IQR={scale:.3f}")

Each band gets its own centre and spread — proof that scaling is independent per column and no single band dominates.

4. Apply the persisted scaler at inference

At inference you load the bundle and transform new tiles without ever refitting. This is what keeps production consistent with training and stops leakage.

import joblib
import rioxarray

bundle = joblib.load("models/spectral_scaler.joblib")

tile = rioxarray.open_rasterio("new_tile.tif")
bands, h, w = tile.shape
X_new = tile.values.reshape(bands, -1).T

X_new_scaled = apply_spectral_scaler(X_new, bundle)

# Reshape back to (band, y, x) for the model or for writing out
scaled_cube = X_new_scaled.T.reshape(bands, h, w)

5. Wire the fit into a cross-validation loop

Fit inside each training fold, never before the split, so the held-out fold never influences the median or IQR.

from sklearn.model_selection import KFold

for train_idx, val_idx in KFold(n_splits=5).split(X):
    bundle = fit_spectral_scaler(X[train_idx], log_bands=(0,))
    X_tr = bundle["scaler"].transform(_prep(X[train_idx], bundle))
    X_va = apply_spectral_scaler(X[val_idx], bundle)   # transform, do NOT fit
    # ... fit and score the model on X_tr / X_va ...

How Robust Scaling Anchors to the Distribution

The diagram contrasts how StandardScaler and RobustScaler respond to the same skewed band. The mean and standard deviation are dragged toward the outlier tail; the median and IQR stay anchored to the dense core of real pixels.

StandardScaler versus RobustScaler on a skewed spectral band A right-skewed histogram of pixel values. The mean is pulled toward the outlier tail while the median stays on the dense peak. StandardScaler subtracts the mean and divides by standard deviation; RobustScaler subtracts the median and divides by the interquartile range, keeping the informative bulk centred. Right-skewed spectral band (e.g. SAR backscatter) pixel value outliers median mean (pulled right) IQR (25th–75th) StandardScaler (x − mean) / std outliers distort both terms RobustScaler (x − median) / IQR anchored to the dense core

Verification

Scaling is a data transform, so verify it against the data rather than trusting the call succeeded. Confirm each band’s scaled median sits at zero and that the outlier tail no longer dominates the spread.

import numpy as np

X_scaled = apply_spectral_scaler(X_train, bundle)
valid = ~np.isnan(X_scaled).any(axis=1)
Xv = X_scaled[valid]

# 1. Median of every band should be ~0 after robust centering.
medians = np.median(Xv, axis=0)
assert np.allclose(medians, 0.0, atol=1e-6), medians
print("Per-band scaled medians:", np.round(medians, 4))

# 2. The middle 50% of each band should span ~1.0 (one IQR).
q75, q25 = np.percentile(Xv, [75, 25], axis=0)
print("Per-band scaled IQR:", np.round(q75 - q25, 3))   # ~1.0 each

# 3. Outliers are recentred, not amplified: compare tail vs core spread.
core = np.percentile(np.abs(Xv), 90, axis=0)
print("90th-percentile magnitude per band:", np.round(core, 2))

A per-band median that rounds to zero and an IQR near one confirm the scaler is anchored to the training distribution as intended. If a band’s median drifts from zero on a new tile, its distribution has shifted relative to training — a covariate-shift signal worth monitoring rather than a reason to refit. For memory-efficient handling of the pixel matrices these transforms operate on, see Optimizing Memory Usage for Large Vector Datasets.

FAQ

When should I log-transform a band before applying RobustScaler?

Log or power transform first when a band is multiplicatively skewed and strictly positive — SAR backscatter in linear power units, aerosol optical depth, or radiance are classic cases. RobustScaler recentres and rescales but preserves the shape of the distribution, so a long right tail stays long. Apply np.log1p or a Yeo-Johnson PowerTransformer to compress the tail, then fit RobustScaler on the transformed values. If a band is already roughly symmetric, skip the transform and scale directly.

Should I scale each spectral band separately or all bands together?

Scale each band separately. Spectral bands occupy different physical ranges and have different outlier structures — a thermal band in kelvin, a reflectance band in 0–1, and SAR backscatter in decibels share no common median or IQR. Fitting one global scaler across flattened bands lets the widest band dominate and destroys the per-band centering. Passing a 2-D (pixels, bands) array to a single RobustScaler already fits independent statistics per column, which is the behaviour you want.

How do I stop the scaler from leaking test data into training?

Call fit only on the training split, then transform both training and validation tiles with that fitted scaler. The median and IQR are statistics estimated from data, so fitting on the full dataset — or refitting per tile at inference — leaks distribution information across the split boundary and inflates cross-validation scores. Persist the fitted scaler with joblib and load it unchanged for every new tile so production inference uses exactly the training-time statistics.

Part of: Feature Scaling for Geospatial Inputs — Spatial Feature Engineering for Machine Learning