Calculating Spatial Coverage Gaps in Raster Datasets

A raster that is missing 1.5% of its pixels rarely announces itself. The file opens, the overviews render, and every cardinality check passes because the tile count is unchanged — yet a band of nodata where a sensor swath dropped out, or a reprojection that quietly clipped an edge, will skew every zonal statistic computed downstream. This page is a focused procedure for detecting and quantifying those coverage gaps deterministically, from a chunked valid-pixel mask through topology-clean gap polygons to threshold-driven alerting. It sits under Geometry Validity & Topology Checks, the structural-integrity gate within the broader Spatial Data Freshness & Quality Metrics program, and is written for the data engineers, GIS platform administrators, and SREs who own raster ingestion integrity.

Coverage-gap detection flow from reference extent to pass, or to gap extraction, topology validation, and alerting A reference extent drawn from the master tile index feeds a chunked coverage calculation that builds a valid-pixel mask. A decision gate asks whether coverage is at least 98.5 percent. If yes, the raster passes. If no, gap polygons are extracted, validated with ST_MakeValid, and routed to an alert and incident step that re-fetches the affected tiles. Reference extent master tile index Chunked coverage calc valid-pixel mask Extract gap polygons connected components Validate topology ST_MakeValid Alert & incident re-fetch tiles Pass emit coverage_ratio Coverage ≥ 98.5%? no yes

Problem Framing

A coverage gap is any region inside the expected extent of a raster that carries no valid observation. It arises from three distinct pipeline stages, and conflating them is the first mistake teams make. At ingestion, a partial write or an interrupted tile-publish job leaves contiguous blocks of nodata. During reprojection, a gdalwarp pass with a mismatched target extent silently trims rows or columns at the margins. At the sensor, cloud cover, swath dropouts, or instrument downtime produce gaps that are real absences of data rather than processing faults — and those must be classified, not “repaired.”

The signal that distinguishes a gap from healthy data is the ratio of valid pixels to expected pixels, evaluated against a known reference boundary rather than against the raster’s own self-reported extent. A file can be internally consistent and still cover the wrong area, which is exactly why coordinate reference system drift must be ruled out before any pixel is counted. Running Coordinate Reference System Validation ahead of the coverage calculation prevents a projection fault from being misread as missing data — an origin offset of a fraction of a degree shifts the whole grid and manufactures a phantom gap along one edge. This procedure is the raster-domain complement to vector Spatial Coverage & Extent Monitoring, which answers the same “does this dataset cover what it should” question for feature collections.

Anchor the calculation to a reference polygon — typically the master tile index or a contractual service boundary — stored in a spatially indexed table so bounding-box intersection stays fast. That baseline is the ground truth for every difference operation that follows.

Implementation

The core measurement is a windowed read that never loads the whole raster into memory, a binary valid-pixel mask, and a morphological erosion that removes resampling artifacts along compression and tile boundaries before the ratio is computed. The function below uses rasterio windowed reads so a multi-terabyte cloud-optimized GeoTIFF is processed in cache-friendly chunks, and returns both the pass/fail verdict and the measured ratio for emission as a metric.

import rasterio
import numpy as np
from scipy import ndimage
from rasterio.windows import Window

def calculate_coverage_gaps(
    raster_path: str,
    coverage_threshold: float = 0.985,
    erosion_radius: int = 3
) -> tuple[bool, float]:
    with rasterio.open(raster_path) as src:
        total_pixels = src.width * src.height
        valid_pixels = 0
        chunk_size = 2048  # windowed read; bounds peak RSS regardless of raster size

        for i in range(0, src.height, chunk_size):
            for j in range(0, src.width, chunk_size):
                width = min(chunk_size, src.width - j)
                height = min(chunk_size, src.height - i)
                window = Window(j, i, width, height)

                # Band 1 -> binary mask: 1 = valid observation, 0 = nodata/missing
                band = src.read(1, window=window)
                nodata = src.nodata
                mask = (band != nodata).astype(np.uint8) if nodata is not None else np.ones_like(band, dtype=np.uint8)

                # Erode edges so resampling/compression fringes are not counted as coverage
                kernel = np.ones(
                    (erosion_radius * 2 + 1, erosion_radius * 2 + 1),
                    dtype=np.uint8
                )
                eroded = ndimage.binary_erosion(mask, structure=kernel)
                valid_pixels += int(np.sum(eroded))

        coverage_ratio = valid_pixels / total_pixels if total_pixels > 0 else 0.0
        return coverage_ratio >= coverage_threshold, coverage_ratio

Line by line: chunk_size = 2048 caps peak resident memory so the same code runs on a 200 MB tile and a 2 TB archive. The mask collapses each band to 1 for a real observation and 0 for nodata; when the raster declares no nodata value, every pixel is treated as valid rather than guessing. The erosion step with a three-pixel radius is the part teams most often omit — without it, the soft anti-aliased fringe of a resampled tile reads as partial coverage and pushes the ratio just under threshold, generating false incidents on otherwise healthy data. The returned coverage_ratio is what you publish as gis.spatial.raster.coverage_ratio.

The coverage ratio itself is the simple proportion the threshold is checked against:

C=PvalidPexpected,alert when C<τ    (τ=0.985)C = \frac{P_{\text{valid}}}{P_{\text{expected}}}, \qquad \text{alert when } C < \tau \;\;(\tau = 0.985)

where PvalidP_{\text{valid}} is the eroded valid-pixel count and PexpectedP_{\text{expected}} is the total pixel count across the reference extent. To locate where the gap is rather than only that one exists, run connected-component labeling (scipy.ndimage.label) on the inverse mask, take the bounding box of each component above a minimum-pixel floor, and export those regions as a GeoJSON or Parquet gap layer. Those polygons must clear the same ST_IsValid / ST_MakeValid discipline as any other geometry before they reach the operational datastore — vectorized gap outlines from a sensor sweep pattern frequently self-intersect, and a sliver polygon silently corrupts the spatial join that attributes the gap to a tile.

Emit the result with the attribute namespace defined in the Geospatial Metric Taxonomy for ETL so a ratio from this raster job and a validity count from a PostGIS trigger land in the same metric space. The Prometheus rule that consumes it:

groups:
  - name: raster_coverage_alerts
    rules:
      - alert: RasterCoverageGapDetected
        expr: gis_spatial_raster_coverage_ratio < 0.985
        for: 5m
        labels:
          severity: critical
          team: spatial-platform-ops
          routing_key: raster-incident-channel
        annotations:
          summary: "Spatial coverage gap exceeds threshold in {{ $labels.dataset_id }}"
          description: "Coverage ratio dropped to {{ $value }}. Expected > 98.5%. Initiate gap extraction and pipeline validation."
      - alert: RasterCRSDeviation
        expr: gis_spatial_raster_crs_offset_degrees > 0.0001
        for: 1m
        labels:
          severity: warning
          team: data-ingestion-ops
        annotations:
          summary: "Coordinate reference system drift detected"
          description: "Origin offset exceeds 0.0001 degrees. Pipeline halted pending CRS reconciliation."

The RasterCRSDeviation rule fires before the coverage rule deliberately: a CRS warning halts the pipeline so a projection fault never reaches the coverage stage and presents as a phantom gap.

Verification & Testing

Confirm the implementation works by injecting a synthetic gap into a known-good raster and asserting the ratio responds. Burn a rectangular block of nodata into a copy of a passing tile, then assert that calculate_coverage_gaps returns False and a ratio close to the area you removed:

import numpy as np, rasterio

with rasterio.open("known_good.tif") as src:
    profile = src.profile
    data = src.read(1)

# Punch a ~2% nodata hole, larger than the erosion fringe so it is unambiguous
data[1000:1300, 1000:1300] = src.nodata
with rasterio.open("synthetic_gap.tif", "w", **profile) as dst:
    dst.write(data, 1)

ok, ratio = calculate_coverage_gaps("synthetic_gap.tif")
assert ok is False, "threshold should fail on an injected gap"
assert ratio < 0.985, f"expected sub-threshold coverage, got {ratio:.4f}"

Two further checks complete the verification. First, run the unmodified known_good.tif through the function and assert the ratio sits above threshold — this catches an erosion radius set so aggressively that healthy edges fail. Second, validate the exported gap polygons in PostGIS: SELECT ST_IsValid(geom), ST_Area(geom::geography) FROM raster_gaps; should return all-valid geometries whose summed area approximates the injected hole. Aligning the gap timestamp to ingestion windows defined by Temporal Baseline Alignment for Time-Series GIS confirms the calculation is scoring the batch you think it is, not bleeding two acquisition cycles into one comparison.

Annotated raster valid-pixel mask with a nodata gap block, an eroded tile-boundary fringe, and the connected-component bounding box exported as a gap polygon A raster grid of pixels forms the valid mask where each cell is a valid observation. A contiguous rectangular block of nodata pixels marks a coverage gap. Along an internal tile boundary a hatched band shows the three-pixel erosion fringe that is stripped before counting so resampling artifacts are not scored as coverage. A dashed bounding box wraps the gap block as a connected component, and an arrow exports it to a validated GeoJSON gap polygon on the right. Valid-pixel mask → connected-component gap polygon mask = 1 · valid observation nodata = 0 3 px erosion fringe stripped at tile boundary component bbox export GeoJSON gap polygon ST_MakeValid · ST_IsValid valid nodata fringe

Gotchas & Failure Modes

Erosion radius masking small real gaps. The same three-pixel erosion that removes resampling fringe will also erase genuine gaps narrower than the kernel — a one-pixel sensor scan-line dropout disappears entirely. If your sensor produces thin linear voids, drop the radius to one pixel for that dataset and compensate by raising the minimum connected-component floor instead, so you keep noise suppression without dissolving real structure.

nodata not declared, so every gap reads as full coverage. When a raster carries no nodata value in its header, the mask defaults every pixel to valid and the ratio is a constant 1.0 no matter how much data is missing — the check passes on a raster that is half empty. Assert src.nodata is not None as a precondition and fail loud when it is absent rather than silently trusting the header.

Coverage passing while the geometry covers the wrong area. A raster can score 100% valid pixels and still be shifted off its reference extent by a SRID mismatch, so the gap is outside the counted window entirely. This is the raster analogue of a validity check passing on a wrongly projected vector; gate CRS validation ahead of coverage and treat a green coverage ratio as meaningful only once the datum is confirmed. When the source feed degrades rather than fails outright, route through the tiers in Fallback Chains for Spatial API Failures so a partial tile set is held back instead of being scored as a gap and re-fetched in a loop.

FAQ

Why erode the valid mask before counting pixels?

Resampled and compressed tiles carry a soft fringe at their boundaries where interpolation produces near-nodata values. Counted raw, that fringe reads as partial coverage and pushes the ratio fractionally below threshold, firing false incidents on healthy data. A three-pixel erosion strips the fringe so the ratio reflects real observations. The trade-off is that gaps narrower than the kernel vanish, so tune the radius per sensor.

Should a coverage gap automatically trigger a tile re-fetch?

No. Distinguish the cause first. A gap from a partial write or interrupted publish should re-fetch the missing tiles; a gap from cloud cover or a sensor swath dropout is a real absence of observation and must be classified, not re-fetched. Re-fetching a sensor gap loops forever against a source that has nothing to return. Cross-reference the gap timestamp against the ingestion cadence before any repair.

What coverage threshold should I set?

98.5% is a defensible default for tiled imagery archives, but the right value is workload-specific. Datasets feeding zonal statistics or regulatory reporting warrant a tighter ceiling; opportunistic feeds with known cloud-cover seasons tolerate more. Set the threshold from the rolling distribution of healthy runs rather than a guessed constant, and pair it alongside — not inside — the freshness contract in Tracking Spatial Data Freshness SLAs.

Why validate the extracted gap polygons with ST_IsValid?

Connected-component outlines traced from a noisy mask — especially sensor sweep patterns — frequently self-intersect or produce zero-width slivers. An invalid gap polygon silently corrupts the spatial join that attributes the gap to a tile or maintenance window, so a real gap gets mis-located or dropped. Run ST_MakeValid before persisting, exactly as you would for any vector geometry entering the curated store.

How do I keep the calculation memory-safe on terabyte-scale rasters?

Read in windows rather than loading the full array. The chunk_size = 2048 window bounds peak resident memory independent of raster dimensions, and for cloud-optimized GeoTIFFs you can compute against a coarser overview level when an exact pixel count is not required. Parallelize windows across workers for very large archives, accumulating the valid-pixel counter per chunk.

For authoritative implementation details on cloud-optimized raster formats and spatial metadata standards, consult the GDAL Cloud Optimized GeoTIFF Driver and the OGC GeoTIFF Standard.