Reshape the raster to (n_pixels, n_bands), drop the invalid rows, fit PCA on a random sample, persist the fitted object, and transform the full raster block by block. The mechanics are simple; the parts that go wrong are nodata leaking into the covariance, standardisation applied (or not) without thinking, and a component matrix that is refitted at inference instead of loaded.
This page is the raster-specific recipe. For when reduction helps at all and how to choose a component count, see Dimensionality Reduction for Spatial Data.
Why This Fails in Geospatial ML Pipelines
Nodata destroys a covariance matrix silently. A -9999 in one band puts a row thousands of standard deviations from the mean, and because PCA maximises variance, the first component becomes “is this pixel nodata”. It explains 97% of the variance, looks like a wonderfully compressible dataset, and carries no information at all. Every scene with a different nodata fraction then produces different components.
The second failure is standardisation applied without a reason. StandardScaler before PCA is a habit imported from tabular work. For a set of reflectance bands that already share units, standardising discards genuinely meaningful differences in variance — the shortwave bands vary more than the blue band because the surface really does vary more there. For a mixed stack of reflectance, elevation and slope it is essential, because elevation in metres would otherwise define component one by itself. The rule is about units, not about habit.
Third, refitting at inference. A fitted PCA is a rotation matrix; refit it on a different scene and component one points somewhere else. A model trained on the old rotation then receives inputs that are numerically valid and semantically unrelated. This is the same contract described in scaling features consistently between training and inference, and it fails the same way: quietly.
Core Principles
- Mask before fitting. Drop rows with any invalid band, do not impute them.
- Decide standardisation by units. Same units → centre only; mixed units → standardise.
- Fit on a sample. A few hundred thousand valid pixels is enough.
- Persist the fitted object and load it everywhere else.
- Fix component signs so maps are comparable between runs.
- Keep the explained-variance vector with the artifact; it is how you justify the component count.
Production-Ready Code
from __future__ import annotations
import logging
import joblib
import numpy as np
import rasterio
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
logger = logging.getLogger(__name__)
def raster_to_matrix(path: str, bands: list[int] | None = None
) -> tuple[np.ndarray, np.ndarray, dict]:
"""Read a raster as (n_pixels, n_bands) float32 plus a validity mask."""
with rasterio.open(path) as src:
idx = bands or list(range(1, src.count + 1))
arr = src.read(idx, masked=True).filled(np.nan).astype("float32")
profile = src.profile.copy()
n_bands, h, w = arr.shape
flat = arr.reshape(n_bands, -1).T
valid = np.isfinite(flat).all(axis=1)
logger.info("%s: %d pixels, %.1f%% valid in all %d bands",
path, flat.shape[0], 100 * valid.mean(), n_bands)
return flat, valid, {**profile, "shape2d": (h, w)}
def fit_band_pca(flat: np.ndarray, valid: np.ndarray, n_components: int = 3,
standardise: bool = False, sample: int = 300_000,
seed: int = 0) -> Pipeline:
"""Fit PCA on a random sample of valid pixels.
Args:
standardise: True only when the bands have different units.
"""
rows = np.flatnonzero(valid)
if rows.size < 1_000:
raise ValueError(f"only {rows.size} valid pixels — too few to fit PCA")
take = np.random.default_rng(seed).choice(rows, size=min(sample, rows.size),
replace=False)
steps = ([("scale", StandardScaler())] if standardise else []) + [
("pca", PCA(n_components=n_components, svd_solver="full", random_state=seed))]
pipe = Pipeline(steps).fit(flat[take])
pca = pipe.named_steps["pca"]
# Eigenvectors are sign-ambiguous: pin the largest loading positive.
flip = np.sign(pca.components_[np.arange(n_components),
np.argmax(np.abs(pca.components_), axis=1)])
pca.components_ *= flip[:, None]
logger.info("explained variance ratio: %s",
np.round(pca.explained_variance_ratio_, 4).tolist())
return pipe
def transform_raster(src_path: str, out_path: str, pipe: Pipeline,
bands: list[int] | None = None) -> None:
"""Apply a fitted pipeline block by block and write a component raster."""
n_components = pipe.named_steps["pca"].n_components_
with rasterio.open(src_path) as src:
idx = bands or list(range(1, src.count + 1))
profile = src.profile.copy()
profile.update(count=n_components, dtype="float32", nodata=np.nan,
compress="deflate", tiled=True, blockxsize=512, blockysize=512)
with rasterio.open(out_path, "w", **profile) as dst:
for _, window in src.block_windows(1):
block = src.read(idx, window=window, masked=True).filled(np.nan)
nb, bh, bw = block.shape
flat = block.reshape(nb, -1).T.astype("float32")
ok = np.isfinite(flat).all(axis=1)
out = np.full((flat.shape[0], n_components), np.nan, dtype="float32")
if ok.any():
out[ok] = pipe.transform(flat[ok]).astype("float32")
dst.write(out.T.reshape(n_components, bh, bw), window=window,
indexes=list(range(1, n_components + 1)))
for i in range(n_components):
logger.info("band %d = PC%d", i + 1, i + 1)Step-by-Step Walkthrough
Step 1 — Read, and look at the validity fraction
flat, valid, profile = raster_to_matrix("s2_stack.tif")
print(f"valid pixels: {valid.sum():,} of {valid.size:,}")If validity is low, find out why before fitting. A 40% valid scene usually means cloud masking has been applied to some bands and not others, which will bias the covariance towards whatever the clear pixels happen to be.
Step 2 — Decide standardisation from the band units
import numpy as np
print(np.nanstd(flat[valid], axis=0).round(4))If the standard deviations differ by an order of magnitude or more, the bands are not on a common scale and you should standardise. If they are within a factor of two or three, centring alone preserves the real structure.
Step 3 — Fit and inspect the loadings
pipe = fit_band_pca(flat, valid, n_components=3, standardise=False, seed=0)
pca = pipe.named_steps["pca"]
import pandas as pd
BANDS = ["blue", "green", "red", "nir", "swir1", "swir2"]
loadings = pd.DataFrame(pca.components_, columns=BANDS,
index=[f"PC{i+1}" for i in range(pca.n_components_)])
print(loadings.round(3))
print("cumulative variance:", np.cumsum(pca.explained_variance_ratio_).round(4))The loadings are the interpretation. On a typical optical stack, PC1 loads positively on everything (overall brightness), PC2 contrasts near-infrared against visible (greenness), and PC3 contrasts shortwave against near-infrared (moisture). If PC1 loads almost entirely on one band, that band has an outlier problem — usually unmasked nodata.
Step 4 — Persist, then transform
import joblib
joblib.dump(pipe, "artifacts/band_pca.joblib")
transform_raster("s2_stack.tif", "s2_pca.tif", pipe)At inference, joblib.load this file. Never call fit again outside the training script — a guard is worth adding in the serving module.
Verification
import numpy as np
import pytest
def test_components_are_orthonormal():
pipe = fit_band_pca(flat, valid, n_components=3)
c = pipe.named_steps["pca"].components_
assert np.allclose(c @ c.T, np.eye(3), atol=1e-6)
def test_nodata_rows_never_reach_the_fit():
dirty = flat.copy()
dirty[0] = -9999.0
bad_valid = valid.copy()
bad_valid[0] = False
pipe = fit_band_pca(dirty, bad_valid, n_components=2, sample=5_000)
assert pipe.named_steps["pca"].explained_variance_ratio_[0] < 0.95
def test_transform_is_deterministic_across_runs():
a = fit_band_pca(flat, valid, n_components=3, seed=0)
b = fit_band_pca(flat, valid, n_components=3, seed=0)
assert np.allclose(a.named_steps["pca"].components_,
b.named_steps["pca"].components_)
def test_signs_are_pinned():
pca = fit_band_pca(flat, valid, n_components=3).named_steps["pca"]
dominant = pca.components_[np.arange(3), np.argmax(np.abs(pca.components_), axis=1)]
assert (dominant > 0).all()
def test_output_raster_matches_the_input_grid(tmp_path):
out = tmp_path / "pca.tif"
transform_raster("s2_stack.tif", str(out), fit_band_pca(flat, valid, 3))
with rasterio.open("s2_stack.tif") as a, rasterio.open(str(out)) as b:
assert a.transform.almost_equals(b.transform, precision=1e-6)
assert (a.width, a.height) == (b.width, b.height)Visually, render PC1–PC3 as an RGB composite and compare against the true-colour image. A good decomposition shows recognisable landscape structure with less correlation between the three channels; a bad one shows the cloud mask.
FAQ
Should I standardise bands before PCA?
Standardise when units differ — reflectance next to elevation and slope — because PCA maximises variance and metres would otherwise define component one. For a set of reflectance bands already on the same scale, centre only, which preserves their real relative variances.
Can I fit PCA on the whole scene?
Rarely necessary. A random sample of a few hundred thousand valid pixels estimates the covariance as well as the full scene and fits in memory; transform the rest block by block.
Why do the component signs flip between runs?
Eigenvectors are defined only up to sign. It is harmless for a model but confusing on a map, so pin the largest-magnitude loading positive as the code above does.
Related
- Dimensionality Reduction for Spatial Data — when reduction helps and how many components to keep
- Raster Band Math and Index Calculation — the interpretable alternative to component bands
- Building a scikit-learn Pipeline for Raster Features — packaging this transform with the model
- Model Explainability for Spatial Predictions — explaining a model whose inputs are components
Part of: Dimensionality Reduction for Spatial Data Part of: Training Geospatial Predictive Models in Python