Writing GDAL-Based Audit Scripts for Quality Policies

You have a written data quality policy — a document that says parcel layers must be in EPSG:27700, must contain no invalid geometries, and must have a non-null parcel_id on every feature — and you need a script that proves whether a given dataset actually satisfies it. This page shows how to turn those clauses into runnable audit scripts with GDAL (Geospatial Data Abstraction Library) and its vector component OGR, using ogrinfo for fast metadata probes, the osgeo.ogr Python API for per-feature checks, and a machine-readable JSON report as the output. It is the executable counterpart to defining spatial data quality policies: the policy defines what good looks like, and the GDAL script decides whether this file passes.

Prerequisites

• GDAL 3.6+ with Python bindings (pip install gdal==$(gdal-config --version) or, more reliably, conda install gdal). Confirm with python -c "from osgeo import gdal; print(gdal.__version__)".
• ogrinfo and ogr2ogr on the PATH — these ship with the GDAL binaries. The -json output flag on ogrinfo requires GDAL 3.7+; on older builds, parse the text output or upgrade.
• PROJ 9.0+ — GDAL delegates all coordinate reference system (CRS) handling to PROJ. A stale PROJ install produces wrong authority codes, so keep the two in step.
• Enable exceptions explicitly. From GDAL 3.7 the Python API warns that silent-error mode is deprecated; call gdal.UseExceptions() at the top of every script so a failed Open raises instead of returning None.
• Gotcha — layer name is not the file name. A GeoPackage can hold many layers; always pass the layer name explicitly to ogrinfo and GetLayerByName(), never assume it matches the file stem.

Step-by-Step Procedure

Step 1 — Encode the written policy as thresholds

The audit script cannot read prose. Translate each policy clause into a structured record with a metric name, a comparison operator, a threshold, and a severity. This file is the contract between the policy author and the script.

# policy_thresholds.py — the machine-readable form of the written policy
POLICY = {
    "policy_id": "PARCELS-QC-2026",
    "expected_epsg": 27700,
    "expected_geometry_type": "Polygon",
    "checks": [
        {"id": "CRS-001",  "metric": "epsg_code",        "op": "==", "value": 27700, "severity": "critical"},
        {"id": "GEOM-001", "metric": "invalid_geom_ratio","op": "<=", "value": 0.0,   "severity": "critical"},
        {"id": "COMP-001", "metric": "null_ratio.parcel_id", "op": "<=", "value": 0.0, "severity": "critical"},
        {"id": "COMP-002", "metric": "null_ratio.land_use",  "op": "<=", "value": 0.02, "severity": "warning"},
    ],
}

Verification: Every clause in the written policy document should map to exactly one entry in checks. Read the policy aloud and confirm no sentence beginning “must” or “shall” is left without a corresponding row — this one-to-one mapping is what makes the audit defensible.

Step 2 — Probe the layer metadata with ogrinfo

Before iterating features, capture the cheap facts — feature count, geometry type, CRS, and field list — in one summary call. The -so flag means “summary only” and skips geometry reads.

import json
import subprocess

def probe_layer(source: str, layer: str) -> dict:
    """Return ogrinfo summary metadata as a parsed dict (GDAL 3.7+ -json)."""
    result = subprocess.run(
        ["ogrinfo", "-so", "-json", source, layer],
        capture_output=True, text=True, check=True,
    )
    info = json.loads(result.stdout)
    layer_info = info["layers"][0]
    crs = layer_info.get("geometryFields", [{}])[0].get("coordinateSystem", {})
    return {
        "feature_count": layer_info["featureCount"],
        "geometry_type": layer_info["geometryFields"][0]["type"],
        "crs_wkt": crs.get("wkt", ""),
        "fields": [f["name"] for f in layer_info["fields"]],
    }

Verification: Run probe_layer("parcels.gpkg", "parcels") and confirm feature_count is greater than zero and fields contains every attribute your null_ratio checks reference. A missing field here is itself a policy failure — a check that silently skips an absent column is worse than no check.

Step 3 — Open the source with the osgeo.ogr API

ogrinfo is enough for metadata, but per-feature geometry and attribute tests need programmatic access. Open the datasource once and reuse the layer handle.

from osgeo import ogr, osr
ogr.UseExceptions()

def open_layer(source: str, layer_name: str):
    ds = ogr.Open(source, 0)  # 0 = read-only
    if ds is None:
        raise FileNotFoundError(f"OGR could not open {source}")
    layer = ds.GetLayerByName(layer_name)
    if layer is None:
        raise KeyError(f"Layer {layer_name!r} not found in {source}")
    return ds, layer  # keep ds alive; layer is invalid once ds is GC'd

Verification: Confirm layer.GetFeatureCount() matches the feature_count from Step 2. A mismatch usually means an attribute or spatial filter is still attached to the layer from an earlier call — reset it with layer.SetAttributeFilter(None) and layer.SetSpatialFilter(None).

Step 4 — Run the geometry, CRS, and attribute checks

Walk the layer once, accumulating the raw counts every metric needs. A single pass keeps the audit fast even on large files because OGR streams features rather than materializing them all.

def measure_layer(layer, tracked_fields: list[str]) -> dict:
    """Single-pass scan collecting the metrics the policy references."""
    srs = layer.GetSpatialRef()
    epsg = int(srs.GetAuthorityCode(None)) if srs and srs.GetAuthorityCode(None) else None

    total = 0
    invalid_geom = 0
    null_counts = {f: 0 for f in tracked_fields}

    layer.ResetReading()
    for feature in layer:
        total += 1
        geom = feature.GetGeometryRef()
        if geom is None or not geom.IsValid():
            invalid_geom += 1
        for field in tracked_fields:
            value = feature.GetField(field) if feature.GetFieldIndex(field) >= 0 else None
            if value is None or (isinstance(value, str) and value.strip() == ""):
                null_counts[field] += 1

    denom = max(total, 1)
    return {
        "epsg_code": epsg,
        "invalid_geom_ratio": invalid_geom / denom,
        **{f"null_ratio.{f}": null_counts[f] / denom for f in tracked_fields},
    }

Verification: On a dataset you know is clean, invalid_geom_ratio must be 0.0 and epsg_code must equal POLICY["expected_epsg"]. Because OGR calls GEOS for IsValid(), this ratio matches what PostGIS ST_IsValid and Shapely would report on the same features — see implementing Shapely geometry checks in Python for the equivalent vectorized approach.

Step 5 — Emit a machine-readable report and set the exit code

Compare each measured metric to its threshold, serialize the outcome to JSON, and exit non-zero if any critical check failed so the script can gate a pipeline.

import operator, sys, json
from datetime import datetime, timezone

OPS = {"==": operator.eq, "<=": operator.le, ">=": operator.ge, "<": operator.lt, ">": operator.gt}

def audit(source: str, layer_name: str, policy: dict) -> dict:
    tracked = [c["metric"].split(".", 1)[1] for c in policy["checks"] if c["metric"].startswith("null_ratio.")]
    ds, layer = open_layer(source, layer_name)
    measured = measure_layer(layer, tracked)

    results, critical_failed = [], False
    for check in policy["checks"]:
        actual = measured.get(check["metric"])
        passed = actual is not None and OPS[check["op"]](actual, check["value"])
        if not passed and check["severity"] == "critical":
            critical_failed = True
        results.append({
            "id": check["id"], "metric": check["metric"], "severity": check["severity"],
            "measured": actual, "threshold": f'{check["op"]} {check["value"]}', "passed": passed,
        })

    report = {
        "policy_id": policy["policy_id"], "source": source, "layer": layer_name,
        "audited_at": datetime.now(timezone.utc).isoformat(),
        "passed": not critical_failed, "results": results,
    }
    return report

if __name__ == "__main__":
    from policy_thresholds import POLICY
    report = audit(sys.argv[1], sys.argv[2], POLICY)
    print(json.dumps(report, indent=2))
    sys.exit(0 if report["passed"] else 1)

Verification: Run the script against a known-bad file and confirm the process exits 1 and the JSON results array names the failing check id. In a shell, python audit.py bad.gpkg parcels; echo $? should print 1.

Interpreting Results

The report is designed to be consumed by both machines and auditors. Each entry pairs the measured value with the threshold it was tested against, so a reviewer never has to re-run the tool to understand a verdict. The top-level passed flag reflects only critical checks — warnings are recorded but do not fail the gate.

When the same layer is audited on every ingestion, store each JSON report keyed by dataset version. The trend of invalid_geom_ratio over time is an early-warning signal that an upstream source has regressed, and the archived reports are the evidence artifacts your compliance framework alignment process needs at audit time.

Gotchas & Edge Cases

GetAuthorityCode(None) returns None for custom projections. A layer defined with a bare WKT string but no EPSG registration has no authority code, so the CRS check reads null and fails. This is usually correct behaviour — an unidentifiable CRS should fail policy — but if you legitimately use a custom grid, match on the PROJ string or a datum name instead.

IsValid() is XY-only. OGR evaluates validity in the horizontal plane and ignores Z. A dataset with corrupt elevations passes the geometry check. If elevation integrity is in policy, read geom.GetZ() per vertex and add a separate bounds metric.

Field-level null semantics vary by driver. In a Shapefile an empty string and a true null are indistinguishable for text fields, whereas GeoPackage preserves SQL NULL. The measure_layer function above treats blank strings as null to normalize this, but document the choice — a policy that counts "" as present will disagree with one that does not.

ogr2ogr can pre-filter but also silently transform. It is tempting to use ogr2ogr -where to extract failing features for a remediation queue, but note that ogr2ogr reprojects to the output CRS by default when -t_srs is given. Never run the audit against an ogr2ogr output unless you intend to audit the transformed data, not the source.

Large GeoPackages need ResetReading() between passes. If you scan the layer twice — once for a count, once for checks — the second loop yields nothing unless you call layer.ResetReading(). The single-pass measure_layer above avoids this, but any refactor that adds a second loop will silently under-count.

When to Escalate

A single-pass OGR script is the right tool for per-feature, per-layer policy gates. Move to a heavier method when:

• Cross-feature topology is in policy. Rules like “no two parcels overlap” require pairwise spatial predicates, not per-feature validity. Push these into PostGIS with a GiST index, or build them in a rule engine with GeoPandas that can hold a spatial index across features.
• The catalogue spans dozens of layers with distinct policies. A flat script per layer becomes unmanageable; adopt a declarative registry that maps each layer to its threshold set, as described in defining spatial data quality policies.
• You need regulator-facing conformance statements. Pass/fail JSON is enough for a gate but not for a formal quality report. Encode results against a recognized standard as covered in compliance framework alignment.
• Feature counts exceed a few million per file. OGR streaming stays memory-bounded, but a pure-Python per-feature loop becomes CPU-bound; move geometry checks into PostGIS ST_IsValid over indexed storage.

Frequently Asked Questions

Related:

• Defining Spatial Data Quality Policies — the written policies these scripts enforce
• Building Rule Engines with GeoPandas — when audits grow into a full declarative rule engine
• Compliance Framework Alignment — turning pass/fail output into standards-based conformance evidence

Back to Defining Spatial Data Quality Policies