Auditing CRS and nodata Drift Across a Collection

To find the moment a collection’s conventions changed, build a header-only inventory and count distinct values per field, ordered by date:

import pandas as pd

inv = pd.read_parquet("inventory.parquet")          # one row per object, headers only
drift = (inv.groupby("collection")[["dtype", "nodata", "crs", "compress", "blockxsize"]]
            .nunique()
            .sort_values("nodata", ascending=False))
print(drift[drift.max(axis=1) > 1])                  # collections with more than one convention

This is the analysis that makes the sweep in Extracting and Parsing Raster Metadata worth keeping rather than discarding.

Why This Arises in Remote Sensing Workflows

Archives change underneath their consumers. A provider reprocesses a mission to a new baseline and the scale factor or offset moves. An internal pipeline is refactored and starts declaring nodata=0 where it previously declared nothing. A cost review switches a collection from DEFLATE to ZSTD. A backfill re-ingests three years of data with a different tile size.

None of these break a read. Every affected file opens, returns pixels and passes any test that checks values. What changes is an assumption the downstream code made silently — that the fill value is always the same, that every scene shares a dtype, that a windowed read is cheap — and the symptom appears somewhere else entirely, as a biased statistic or a job that suddenly runs four times slower.

An audit is the only practical detector, because the evidence is distributed across thousands of headers and visible only in aggregate. And it is cheap: the fields that drift are exactly the fields a header read already returns, so the audit is a grouping over data you should be collecting anyway.

Drift, drawn against the ingest timeline A collection ingested over three years changes nodata at one point, dtype at another, and compression at a third. Each change is invisible to a reader but breaks a downstream assumption: masking, statistics, and reader compatibility respectively. One collection, three silent changes 2023-01 2024-06 2025-12 nodata = None nodata = 0 masking starts dropping real zeros dtype = uint16 dtype = int16 stacks promote to float64, memory doubles compress = DEFLATE ZSTD an older consumer can no longer open the newest scenes Every boundary here is a date the audit can name — and nobody would find it by reading one file.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 Header, tag and structure access
pandas ≥2.0 Grouping, diffing, parquet I/O
pyarrow ≥14 Parquet engine for the inventory 
pip install "rasterio>=1.3.0" "pandas>=2.0" "pyarrow>=14"

Complete Working Example

The audit has three parts: collect, summarise, and compare against the previous run. Collection is a header-only sweep; the rest is dataframe work.

Drift classes and what each one breaks Nodata drift breaks masking, dtype drift breaks stacking and doubles memory, compression drift breaks older readers, tiling drift breaks windowed read performance, and scale drift changes the values themselves — the only class that is invisible to every structural check. Five things that drift, and what each one costs field symptom downstream severity nodata fill values enter statistics high — biases every number scale / offset a step change in a time series highest — looks like real change dtype stacks promote to float64 medium — memory and speed block size windowed reads become strip reads medium — cost and latency compression an older consumer cannot open it low, until it is not The scale row is the one worth alerting on: everything else is slow or expensive, not wrong. 
from concurrent.futures import ThreadPoolExecutor, as_completed

import pandas as pd
import rasterio

STRUCTURAL = ["dtype", "nodata", "crs", "compress", "blockxsize", "blockysize", "count"]


def header_row(uri: str) -> dict:
    """One inventory row. Reads the header only — never src.read()."""
    try:
        with rasterio.open(uri) as src:
            block = src.block_shapes[0]
            return {
                "uri": uri,
                "dtype": src.dtypes[0],
                "nodata": src.nodata,
                "crs": str(src.crs),
                "compress": (src.profile.get("compress") or "none").lower(),
                "blockxsize": block[1],
                "blockysize": block[0],
                "count": src.count,
                "overviews": len(src.overviews(1)),
                "scale": (src.scales or [1.0])[0],
                "offset": (src.offsets or [0.0])[0],
                "error": None,
            }
    except Exception as exc:                      # unreadable files are findings, not crashes
        return {"uri": uri, "error": repr(exc)}


def collect(uris: list[str], workers: int = 16) -> pd.DataFrame:
    """Header sweep. Threads, because every worker is waiting on the network."""
    rows = []
    with ThreadPoolExecutor(max_workers=workers) as pool:
        futures = [pool.submit(header_row, u) for u in uris]
        for fut in as_completed(futures):
            rows.append(fut.result())
    return pd.DataFrame(rows)


def summarise(inv: pd.DataFrame, group: str = "collection") -> pd.DataFrame:
    """Distinct values per structural field, per group."""
    present = [c for c in STRUCTURAL if c in inv.columns]
    counts = inv.groupby(group)[present].nunique()
    counts["files"] = inv.groupby(group).size()
    return counts


def first_seen(inv: pd.DataFrame, field: str, date_col: str = "date") -> pd.DataFrame:
    """The earliest date each distinct value of `field` appears — i.e. when it changed."""
    ordered = inv.sort_values(date_col)
    return (ordered.groupby([field])[date_col].agg(["min", "max", "count"])
            .rename(columns={"min": "first_seen", "max": "last_seen", "count": "files"})
            .sort_values("first_seen"))


if __name__ == "__main__":
    inv = collect(uris)                                # uris from a bucket listing or a STAC search
    inv["collection"] = inv["uri"].str.split("/").str[3]
    inv["date"] = inv["uri"].str.extract(r"(\d{8})")[0]

    print(summarise(inv))
    print(first_seen(inv[inv.collection == "s2_l2a_ndvi"], "nodata"))

    inv.to_parquet("inventory_2026-08-06.parquet", index=False)

first_seen is the function that turns “this collection has two nodata values” into “the second value appeared on 2024-06-14”, which is what makes the finding actionable — you can go and look at what shipped that day.

Variant Patterns

1. Diffing two sweeps

The most useful alert is not “there are two values” but “there is a value now that was not there last week”.

def new_values(previous: pd.DataFrame, current: pd.DataFrame, field: str) -> set:
    """Values of `field` present in the current sweep but absent from the previous one."""
    return set(current[field].dropna().unique()) - set(previous[field].dropna().unique())


for field in ["dtype", "nodata", "crs", "compress", "scale"]:
    added = new_values(prev_inv, inv, field)
    if added:
        print(f"NEW {field}: {sorted(added)}")

Keeping every sweep in a dated parquet file makes this a two-line query and gives you a history to look back through when something is finally noticed.

2. Auditing scale factors specifically

Scaling drift is the most damaging kind, because it changes values rather than performance, and the value change is a clean multiplicative factor that looks like a real trend.

What a scaling change looks like downstream The median reflectance of a stable target follows a seasonal cycle until a processing baseline change introduces an offset, after which the whole series shifts by a constant amount. Without an audit this reads as an abrupt environmental change; with one, the date matches a metadata change exactly. Median reflectance over a stable target 2024-06-14 baseline change: offset −1000 introduced seasonal variation, stable level 2023 2025 The step is exactly the offset, and it appears on the day the audit records a new scale/offset pair. Harmonising the archive is easy once the date is known; finding the date without an inventory is not. 
scaling = inv.assign(pair=inv["scale"].astype(str) + "/" + inv["offset"].astype(str))
print(first_seen(scaling, "pair"))

The consequences of getting scaling wrong, and the fix, are covered in Handling nodata and Scale Factors in Band Math.

3. Auditing structure, not just semantics

The same sweep answers “why did this collection get slow?” — usually because tiling or overviews changed.

slow = inv[(inv.blockxsize < 256) | (inv.overviews < 3)]
print(slow.groupby("collection").size().sort_values(ascending=False).head())

Files without adequate tiling or overviews are the ones that make windowed reads expensive, as quantified in Writing and Validating Cloud-Optimized GeoTIFFs.

Reading the Audit Without Overreacting

Not all drift is a defect, and an audit that cries wolf gets muted. Three distinctions keep it useful.

Expected multiplicity is not drift. A collection spanning UTM zones has many CRS values by design, and a collection of derived products may legitimately hold several dtypes. Record the expected set per collection and alert only on values outside it.

A change with a clean date boundary is a processing change; a change scattered through time is an ingest bug. The first needs a harmonisation decision, the second needs a fix in the writer. first_seen distinguishes them at a glance, because a processing change shows one date and a bug shows overlapping ranges.

Unreadable files are the most urgent finding and the easiest to overlook, because they appear as an error column rather than as a value change. Count them every run: a collection that gains unreadable objects is losing data, and no amount of downstream defensiveness recovers it.

Common Errors

The sweep is slow and the bucket returns SlowDown

Too many concurrent header reads. Sixteen threads is a reasonable ceiling for a shared bucket; past thirty-two most providers start throttling, and the retries cost more than the parallelism gained.

nodata comes back as nan for some files and None for others

Both mean “no usable declaration” in practice but compare differently. Normalise to None in the inventory before grouping, or the same convention will appear as two distinct values.

The audit reports drift that does not exist

The grouping key is wrong — often a path-derived collection name that changed when the prefix layout changed. Derive the collection from the STAC item where possible, not from a path fragment.

Frequently Asked Questions

Q: Is more than one CRS in a collection always a problem? No. Sentinel-2 spans UTM zones by design, so several EPSG codes are expected. The problem is an unexpected code, or a code that appears only after a certain date, which indicates a change in processing rather than in geography.

Q: How often should the audit run? As often as the collection grows. A daily incremental sweep over new objects costs minutes; a full sweep weekly or monthly catches retrospective reprocessing that rewrites existing objects.

Q: What should the audit do when it finds drift? Report, not block. Drift is often legitimate — a new baseline, a deliberate compression change. The value of the audit is that the change is noticed on the day it happens rather than during an investigation months later.