Utility Preservation Metrics for Masked Maps

Utility preservation metrics quantify how much analytical value survives a privacy-preserving transformation of a geospatial dataset, giving engineers a measurable bridge between theoretical privacy guarantees and operational data fitness.

When spatial datasets undergo masking under the broader framework of Differential Privacy for Location Data, the engineering challenge shifts from mere anonymization to continuous validation. Every noise injection, aggregation, or synthetic trajectory operation introduces distortion. Utility metrics bound that distortion within acceptable thresholds while respecting the privacy budget — the finite store of privacy loss governed by the ε parameter.

Algorithmic Specification # permalink

Each metric category has a formal definition and a recommended operating range. The table below summarises the core metrics, their mathematical meaning, and typical thresholds for production release.

Parameter guidance:

• Thresholds above assume a metric CRS such as EPSG:3857 or a local UTM zone. Distance values are in metres; all calculations must happen in projected space, not on geographic (degree-based) coordinates.
• As the privacy budget ε decreases below 0.5, noise magnitude grows sharply; recalibrate Green/Amber thresholds toward the Amber column for high-noise regimes.
• When Laplace or Gaussian noise is applied to coordinate pairs, mean displacement approximates the noise scale parameter: $\bar{d} \approx b = \Delta / \varepsilon$ for Laplace. This identity lets you sanity-check the pipeline by comparing the observed mean displacement against the theoretical noise scale.

Prerequisites & Data Requirements # permalink

Before running utility validation, verify that each of the following is in place. Skipping any item frequently causes silent metric corruption or false-positive utility scores.

1. Paired datasets — a raw reference dataset and its masked counterpart sharing an identical schema, attribute set, and temporal window, with 1:1 row alignment enforced by a stable join key.
2. CRS alignment in metric space — both datasets must be re-projected to a metric CRS (e.g., EPSG:3857, EPSG:32633) before any distance calculation. Angular coordinates (EPSG:4326) produce invalid Euclidean metrics.
3. Privacy parameter documentation — recorded ε, δ, sensitivity Δ, and the specific noise mechanism (Laplace, Gaussian, planar isotropy) used during masking. These are required for audit trails.
4. Baseline query set — a predefined set of analytical queries representing the real downstream workload: point-in-polygon counts, density heatmaps, nearest-neighbour lookups, hotspot detection. Metrics must be anchored to actual use cases.
5. Python stack — geopandas ≥ 0.13, shapely ≥ 2.0, numpy ≥ 1.24, scipy ≥ 1.10, pyproj ≥ 3.5.
6. Dataset fingerprints — SHA-256 hashes of both datasets recorded before metric computation to guarantee reproducibility and detect tampering in audit reviews.

Step-by-Step Implementation # permalink

Step 1 — Validate and Project # permalink

Enforce CRS parity, re-project to metric coordinates, and hash both datasets before any measurement. Mismatched projections at this stage contaminate every downstream metric.

import hashlib
import geopandas as gpd
import numpy as np

def fingerprint_gdf(gdf: gpd.GeoDataFrame) -> str:
    """Return a SHA-256 fingerprint of the serialised geometry column."""
    raw = gdf.geometry.to_wkb().values.tobytes()
    return hashlib.sha256(raw).hexdigest()

def align_and_project(
    raw_gdf: gpd.GeoDataFrame,
    masked_gdf: gpd.GeoDataFrame,
    metric_crs: str = "EPSG:3857",
) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
    """
    Validate CRS parity, re-project both datasets to a metric CRS,
    and assert row alignment.  Raises on any mismatch.
    """
    if raw_gdf.crs != masked_gdf.crs:
        raise ValueError(
            f"CRS mismatch: raw={raw_gdf.crs!r}, masked={masked_gdf.crs!r}. "
            "Re-project to a common CRS before calling this function."
        )
    if len(raw_gdf) != len(masked_gdf):
        raise ValueError(
            f"Row count mismatch: raw={len(raw_gdf)}, masked={len(masked_gdf)}. "
            "Enforce 1:1 alignment via a stable join key before validation."
        )
    raw_m = raw_gdf.to_crs(metric_crs)
    masked_m = masked_gdf.to_crs(metric_crs)
    return raw_m, masked_m

Privacy note: re-projection does not add or remove privacy guarantees; it only ensures that distance-based metrics are computed in a geometrically consistent space.

Step 2 — Geometric Fidelity Metrics # permalink

Compute point-level displacement and Hausdorff distance to quantify positional distortion.

from scipy.spatial.distance import directed_hausdorff

def geometric_fidelity(
    raw_m: gpd.GeoDataFrame,
    masked_m: gpd.GeoDataFrame,
) -> dict[str, float]:
    """
    Returns displacement statistics and the directed Hausdorff distance
    (both in the units of the input CRS, typically metres).
    """
    raw_coords = np.column_stack([raw_m.geometry.x, raw_m.geometry.y])
    masked_coords = np.column_stack([masked_m.geometry.x, masked_m.geometry.y])

    # Vectorised L2 displacement per point
    displacements = np.linalg.norm(raw_coords - masked_coords, axis=1)

    # Directed Hausdorff: worst-case shape-level distortion
    h_fwd = directed_hausdorff(raw_coords, masked_coords)[0]
    h_rev = directed_hausdorff(masked_coords, raw_coords)[0]

    return {
        "mean_displacement_m": float(np.mean(displacements)),
        "median_displacement_m": float(np.median(displacements)),
        "p95_displacement_m": float(np.percentile(displacements, 95)),
        "hausdorff_m": float(max(h_fwd, h_rev)),
        "n": len(displacements),
    }

For large datasets (> 1 M rows), replace directed_hausdorff with an approximate nearest-neighbour approach using scipy.spatial.cKDTree to avoid O(n²) memory usage.

Step 3 — Distributional Consistency # permalink

Spatial density patterns must survive the masking step even when individual coordinates shift. Kernel density estimation (KDE) correlation and Moran’s I delta are the canonical checks.

from scipy.stats import pearsonr
from scipy.ndimage import gaussian_filter
import numpy as np

def kde_correlation(
    raw_coords: np.ndarray,
    masked_coords: np.ndarray,
    grid_cells: int = 256,
) -> float:
    """
    Estimate KDE on a square grid and return the Pearson r between
    the raw and masked density surfaces.  Higher r = better distributional
    preservation.  grid_cells controls spatial resolution; 256 is a
    practical default for city-scale datasets.
    """
    def _rasterise(coords: np.ndarray, xmin: float, xmax: float,
                   ymin: float, ymax: float) -> np.ndarray:
        ix = np.clip(
            ((coords[:, 0] - xmin) / (xmax - xmin) * grid_cells).astype(int),
            0, grid_cells - 1,
        )
        iy = np.clip(
            ((coords[:, 1] - ymin) / (ymax - ymin) * grid_cells).astype(int),
            0, grid_cells - 1,
        )
        grid = np.zeros((grid_cells, grid_cells), dtype=float)
        np.add.at(grid, (iy, ix), 1)
        return gaussian_filter(grid, sigma=max(1, grid_cells // 32))

    all_coords = np.vstack([raw_coords, masked_coords])
    xmin, ymin = all_coords.min(axis=0)
    xmax, ymax = all_coords.max(axis=0)

    raw_surf = _rasterise(raw_coords, xmin, xmax, ymin, ymax)
    masked_surf = _rasterise(masked_coords, xmin, xmax, ymin, ymax)

    r, _ = pearsonr(raw_surf.ravel(), masked_surf.ravel())
    return float(r)

Privacy note: KDE correlation is computed on already-masked outputs; it does not expose raw coordinates. The raw surface must be generated in a secure enclave and discarded after the metric is recorded.

Step 4 — Query Accuracy Benchmarking # permalink

Relative error across a representative query set is the most operationally meaningful metric. Define queries against your actual downstream workload.

from shapely.geometry import box

def point_in_polygon_relative_error(
    raw_m: gpd.GeoDataFrame,
    masked_m: gpd.GeoDataFrame,
    query_bounds: list[tuple[float, float, float, float]],
) -> dict[str, float]:
    """
    For each (minx, miny, maxx, maxy) bounding box, count points in the
    raw and masked datasets, then return mean and max relative error.

    query_bounds: list of axis-aligned rectangles in the metric CRS.
    Privacy implication: counts below a minimum suppression threshold
    should not be published; set count=0 for cells where raw_count < k_min.
    """
    errors = []
    for bounds in query_bounds:
        poly = box(*bounds)
        q_raw = raw_m.geometry.within(poly).sum()
        q_masked = masked_m.geometry.within(poly).sum()
        if q_raw == 0:
            continue  # skip empty cells; relative error is undefined
        errors.append(abs(q_raw - q_masked) / q_raw)

    arr = np.array(errors)
    return {
        "mean_relative_error": float(np.mean(arr)),
        "max_relative_error": float(np.max(arr)),
        "n_queries": len(arr),
    }

Step 5 — Classify Tier and Emit Audit Record # permalink

Aggregate all metrics into a single audit record with a tier classification.

import json
from datetime import datetime, timezone

def classify_tier(metrics: dict[str, float]) -> str:
    """
    Map metric values to Green / Amber / Red operational tier.
    A single Red indicator promotes the whole dataset to Red.
    """
    checks = [
        metrics.get("mean_displacement_m", 0) <= 25,
        metrics.get("p95_displacement_m", 0) <= 75,
        metrics.get("hausdorff_m", 0) <= 150,
        metrics.get("kde_pearson_r", 0) >= 0.95,
        metrics.get("moran_i_delta", 1) <= 0.05,
        metrics.get("mean_relative_error", 1) <= 0.05,
    ]
    amber_checks = [
        metrics.get("mean_displacement_m", 0) <= 100,
        metrics.get("p95_displacement_m", 0) <= 300,
        metrics.get("hausdorff_m", 0) <= 500,
        metrics.get("kde_pearson_r", 0) >= 0.80,
        metrics.get("moran_i_delta", 1) <= 0.15,
        metrics.get("mean_relative_error", 1) <= 0.15,
    ]
    if all(checks):
        return "Green"
    if all(amber_checks):
        return "Amber"
    return "Red"

def emit_audit_record(
    metrics: dict[str, float],
    tier: str,
    epsilon: float,
    delta: float,
    mechanism: str,
    raw_hash: str,
    masked_hash: str,
) -> str:
    """Return a JSON audit record suitable for a compliance registry."""
    record = {
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "utility_tier": tier,
        "epsilon": epsilon,
        "delta": delta,
        "mechanism": mechanism,
        "raw_dataset_sha256": raw_hash,
        "masked_dataset_sha256": masked_hash,
        "metrics": metrics,
    }
    return json.dumps(record, indent=2)

Validation and Re-identification Testing # permalink

Utility validation and re-identification testing are two sides of the same coin. High utility scores must not come at the expense of residual privacy risk.

Neighbour-count audit. For each masked point, count how many raw points fall within the noise radius $b = \Delta / \varepsilon$. If any masked point is the sole near-neighbour of a raw point, it may be individually traceable. Flag records where the neighbourhood count falls below the k-anonymity grouping threshold.

from scipy.spatial import cKDTree

def neighbour_count_audit(
    raw_coords: np.ndarray,
    masked_coords: np.ndarray,
    radius_m: float,
    k_min: int = 5,
) -> dict[str, object]:
    """
    For each masked point, count raw points within radius_m.
    Returns the fraction of masked points with fewer than k_min neighbours —
    a proxy for re-identification exposure in sparse areas.
    """
    tree = cKDTree(raw_coords)
    counts = tree.query_ball_point(masked_coords, r=radius_m, return_length=True)
    below_k = np.array(counts) < k_min
    return {
        "fraction_below_k": float(below_k.mean()),
        "k_min": k_min,
        "radius_m": radius_m,
        "n_flagged": int(below_k.sum()),
    }

Auxiliary-join simulation. Merge the masked dataset with a synthetic auxiliary table containing public geographic attributes (census block identifiers, building footprint centroids). Measure the join hit rate: if more than 5% of masked points uniquely match a public record, treat the masking radius as insufficient and increase ε-compensation or apply additional spatial fuzzing.

Entropy check on coordinate pairs. Discretise masked coordinates into a fine grid (e.g., 10 m cells) and compute the Shannon entropy of the frequency distribution. A drop greater than 15% relative to the raw entropy signals density collapse, a failure mode where noise has concentrated records into a small number of grid cells and paradoxically made them more identifiable.

Common Failure Modes and Gotchas # permalink

Compliance Alignment # permalink

Utility preservation metrics satisfy several regulatory and standards requirements:

• GDPR Article 5(1)(f) and WP29 Opinion 05/2014 — the anonymisation test requires that data “cannot be used to single out, link, or infer” individuals. Structured utility reports with documented ε values and neighbour-count audits provide concrete evidence in a GDPR data protection impact assessment (DPIA).
• CCPA / CPRA — “deidentified data” must have technical safeguards preventing re-linkage. Query accuracy benchmarks and auxiliary-join simulations directly address the re-linkage prevention requirement.
• NIST SP 800-188 (De-Identification of Government Datasets) — recommends documenting both the privacy technique applied and the resulting analytical utility. The audit record schema above is directly mappable to SP 800-188 Annex B.
• ISO/IEC 20889:2018 — the standard for privacy-enhancing data de-identification defines utility and risk metrics. The metric categories in this guide align with ISO 20889 Section 8 (utility) and Section 9 (re-identification risk).

Document all metric outputs alongside dataset hashes and ε/δ values in a versioned compliance registry. Automate audit record emission in CI/CD pipelines to maintain an unbroken evidence chain for GDPR and CCPA compliance.

Validation Checklist # permalink

Before publishing any masked spatial dataset to a downstream consumer:

• Raw and masked datasets share identical CRS, schema, and temporal scope, verified by hash comparison
• Both datasets re-projected to a metric CRS (EPSG:3857 Both datasets re-projected to a metric CRS (`EPSG:3857` or local UTM) before all metric calculations
• Geometric fidelity: mean displacement, P95 displacement, and Hausdorff distance within operational tier thresholds
• Distributional consistency: KDE Pearson r ≥ 0.80 (Amber) or ≥ 0.95 (Green); Moran's I delta within threshold
• Query accuracy: relative error across baseline query set within operational tier thresholds
• Neighbour-count audit: fraction below k_min Neighbour-count audit: fraction below `k_min` does not exceed 5%
• Auxiliary-join simulation: join hit rate on public data below 5%
• Privacy parameters (ε, δ, Δ, mechanism) documented alongside metric values
• Audit record emitted with SHA-256 hashes and timestamp, stored in compliance registry
• Masked outputs validated with shapely.validation.make_valid() Masked outputs validated with `shapely.validation.make_valid()` before publication
• Utility metrics integrated into CI/CD pipeline with regression alerts on tier degradation

Treat utility preservation as a continuous validation discipline. As accuracy vs. utility tradeoffs shift with changing ε budgets or new analytical workloads, recalibrate thresholds, expand the baseline query set, and re-run the full metric pipeline against each new masked release.

Related # permalink

• Privacy Budget Allocation for Spatial Queries
• Accuracy vs. Utility Tradeoffs in Geospatial DP
• Laplace & Gaussian Noise for Coordinate Data
• Re-identification Risk Assessment for Geospatial Datasets
• Spatial Fuzzing & Buffer Zone Implementation

← Back to Differential Privacy for Location Data