Privacy Risk Scoring Frameworks for GIS

A privacy risk scoring framework for GIS is a quantitative methodology that assigns a normalized composite score to each spatial dataset or feature, enabling automated, threshold-driven decisions about masking, sharing, and retention — replacing ad-hoc subjective assessments with auditable, repeatable engineering controls.

When to Use Risk Scoring vs. Alternative Approaches permalink

Not every privacy problem requires a composite scoring pipeline. The diagram below helps you choose the right technique for your situation.

Choosing a spatial privacy approach: risk scoring vs. alternatives A decision flow diagram showing when to apply composite risk scoring versus simpler approaches such as blanket k-anonymity grouping, coordinate jittering, or differential privacy noise injection. Spatial dataset to evaluate for sharing / publication Multiple sensitivity dimensions vary per record? Yes Composite risk scoring (this page — SG·AS·CE·LP·LS) No Point data only? Yes Coordinate jittering or DP noise No No Blanket k-anonymity grouping per layer Score → tier → automated masking / quarantine / review Validated dataset released

If your dataset has uniform geometry type and a single dominant sensitivity attribute, a simpler technique — k-anonymity grouping for location traces or coordinate jittering — may suffice. Composite risk scoring earns its overhead when records differ meaningfully in granularity, linked attributes, or sharing channel.

Algorithmic Specification permalink

The composite spatial risk score aggregates five orthogonal dimensions, each normalized to [0,1][0, 1]:

R=w1SG+w2AS+w3CE+w4LP+w5LS,i=15wi=1R = w_1 \cdot SG + w_2 \cdot AS + w_3 \cdot CE + w_4 \cdot LP + w_5 \cdot LS, \quad \sum_{i=1}^{5} w_i = 1
Dimension Symbol Description Typical range Notes
Spatial Granularity SGSG Precision of coordinate resolution 0.2–1.0 Sub-metre GPS ≈ 0.95; census-tract centroid ≈ 0.25
Attribute Sensitivity ASAS Privacy impact of non-spatial fields 0.0–1.0 Direct identifiers → 1.0; aggregated counts → 0.0–0.2
Contextual Exposure CECE Audience, sharing channel, retention window 0.0–0.9 Public API → 0.8–0.9; air-gapped analytics → 0.1–0.2
Linkage Potential LPLP Feasibility of cross-referencing auxiliary datasets 0.0–1.0 High overlap with property records or social check-ins → 0.7–1.0
Lifecycle Stage LSLS Data maturity and applied transformations 0.0–0.9 Raw ingestion → 0.8–0.9; aggregated archive → 0.1–0.3

SGSG is derived from the feature’s bounding-box area AA (in km²) using a log-scaled inversion so that higher precision (smaller area) maps to a higher risk score:

SG=1clip ⁣(ln(1+A)ln(1+Amax), 0, 1)SG = 1 - \text{clip}\!\left(\frac{\ln(1 + A)}{\ln(1 + A_{\max})},\ 0,\ 1\right)

where Amax=1000 km2A_{\max} = 1000\ \text{km}^2 is a calibration ceiling appropriate for most municipal datasets. Adjust AmaxA_{\max} for national or global datasets.

ASAS is the maximum sensitivity score across all non-spatial columns for a given record, derived from an organisation-maintained sensitivity map that assigns values to attribute strings (direct identifiers → 1.0, quasi-identifiers → 0.3–0.8, public reference data → 0.0–0.1).

CECE, LPLP, and LSLS require metadata that must be injected at pipeline ingestion time from your data-governance registry. The implementation below shows how to pass them as per-record Series once that registry exists.

Prerequisites & Data Requirements permalink

Before deploying the scoring pipeline, satisfy the following:

  1. Standardized CRS. All geometries must share a single CRS. Mixed projections silently distort the bounding-box area calculation that drives SGSG. Validate with gdf.crs.to_epsg() before ingestion; reproject to EPSG:4326 (WGS 84) for the granularity calculation, then optionally back to a local projected CRS for other operations.
  2. Attribute classification schema. Each non-spatial column must carry a sensitivity tier tag (public, internal, confidential, restricted) stored in your data-governance catalogue. This drives ASAS without requiring hard-coded column names.
  3. Baseline threat inventory. Document known adversary capabilities, available auxiliary datasets, and intended data consumers before calibrating LPLP and CECE weights. Organizations that have already completed re-identification risk assessment for geospatial datasets will find the LPLP inputs immediately available.
  4. Python dependencies. geopandas>=0.13, pandas>=2.0, numpy>=1.24. scikit-learn is optional for clustering-based threshold calibration.
  5. Regulatory mapping. Identify which jurisdictions’ rules apply. The compliance mapping for GDPR and CCPA location data page covers the key provisions that scoring thresholds must reflect.

Step-by-Step Implementation permalink

Step 1: Prepare the GeoDataFrame and sensitivity map permalink

import geopandas as gpd
import pandas as pd
import numpy as np
from typing import Dict

# EPSG:4326 (WGS 84) is required for the bounding-box area calculation.
# Reproject here to avoid silent distortion in the SG dimension.
def load_and_normalize(path: str, crs_target: str = "EPSG:4326") -> gpd.GeoDataFrame:
    gdf = gpd.read_file(path)
    if gdf.crs is None:
        raise ValueError("GeoDataFrame has no CRS — set it before scoring.")
    return gdf.to_crs(crs_target)

The CRS reprojection must happen before any spatial calculation. A dataset projected in EPSG:3857 (Web Mercator) uses metre units; feeding those bounding-box values into the log-scaled formula without conversion produces SGSG scores that are meaningless.

Step 2: Compute the five scoring dimensions permalink

def score_spatial_granularity(gdf: gpd.GeoDataFrame, area_max_km2: float = 1000.0) -> pd.Series:
    """
    Compute SG dimension from bounding-box area.
    Assumes gdf is already in EPSG:4326 (degrees).
    1 degree ≈ 111.32 km on both axes at low latitudes; adequate for most scoring use cases.
    """
    b = gdf.geometry.bounds
    # Convert degree-unit bounding box to approximate square kilometres.
    deg_to_km = 111.32
    area_km2 = (b["maxx"] - b["minx"]) * (b["maxy"] - b["miny"]) * (deg_to_km ** 2)
    log_norm = np.log1p(area_km2) / np.log1p(area_max_km2)
    return 1.0 - np.clip(log_norm, 0.0, 1.0)


def score_attribute_sensitivity(
    gdf: gpd.GeoDataFrame,
    sensitivity_map: Dict[str, float],
) -> pd.Series:
    """
    Compute AS dimension: maximum sensitivity score across non-spatial columns.
    sensitivity_map maps column-value strings to floats in [0, 1].
    Direct identifiers (name, NI number, passport ID) → 1.0.
    Quasi-identifiers (occupation, borough, age band) → 0.3–0.7.
    Public reference data → 0.0–0.1.
    """
    geom_col = gdf.geometry.name
    attr_cols = [c for c in gdf.columns if c != geom_col]
    # Replace attribute values with their sensitivity scores; non-mapped values → 0.0.
    scores = (
        gdf[attr_cols]
        .replace(sensitivity_map)
        .apply(pd.to_numeric, errors="coerce")
        .fillna(0.0)
    )
    return scores.max(axis=1).clip(0.0, 1.0)

Privacy implication of Step 2: scores.max(axis=1) deliberately takes the worst-case attribute rather than an average. A single direct identifier in a row elevates the whole record to high sensitivity regardless of how many low-risk attributes accompany it. This is intentional: an average would mask the presence of the identifier.

Step 3: Assemble the composite score permalink

def calculate_spatial_risk_score(
    gdf: gpd.GeoDataFrame,
    weights: Dict[str, float],
    sensitivity_map: Dict[str, float],
    ce_scores: pd.Series,
    lp_scores: pd.Series,
    ls_scores: pd.Series,
    area_max_km2: float = 1000.0,
) -> pd.Series:
    """
    Compute normalized composite privacy risk scores for a GeoDataFrame.

    Args:
        gdf: GeoDataFrame in EPSG:4326 with point or polygon geometries.
        weights: Keys must be 'SG', 'AS', 'CE', 'LP', 'LS'; values must sum to 1.0.
        sensitivity_map: Maps attribute value strings to sensitivity scores (0–1).
        ce_scores: Per-record contextual exposure scores from metadata registry.
        lp_scores: Per-record linkage potential scores from auxiliary-overlap analysis.
        ls_scores: Per-record lifecycle stage scores from ETL pipeline tags.
        area_max_km2: Calibration ceiling for bounding-box normalization.

    Returns:
        pandas Series of composite scores in [0.0, 1.0] aligned with gdf.index.
    """
    required_keys = {"SG", "AS", "CE", "LP", "LS"}
    if set(weights.keys()) != required_keys:
        raise ValueError(f"weights must have exactly these keys: {required_keys}")
    if not np.isclose(sum(weights.values()), 1.0, atol=1e-4):
        raise ValueError("Weights must sum to 1.0")

    sg = score_spatial_granularity(gdf, area_max_km2)
    as_ = score_attribute_sensitivity(gdf, sensitivity_map)

    composite = (
        weights["SG"] * sg
        + weights["AS"] * as_
        + weights["CE"] * ce_scores.clip(0.0, 1.0)
        + weights["LP"] * lp_scores.clip(0.0, 1.0)
        + weights["LS"] * ls_scores.clip(0.0, 1.0)
    )
    return composite.clip(0.0, 1.0)

Step 4: Threshold and route records permalink

def apply_risk_tiers(scores: pd.Series) -> pd.Series:
    """
    Map composite scores to actionable risk tiers.
    Thresholds are a starting point — calibrate against your baseline distribution.
    """
    bins = [0.0, 0.30, 0.60, 0.80, 1.01]
    labels = ["low", "medium", "high", "critical"]
    return pd.cut(scores, bins=bins, labels=labels, right=True, include_lowest=True)


# Example usage
if __name__ == "__main__":
    gdf = load_and_normalize("incidents.gpkg")

    weights = {"SG": 0.25, "AS": 0.30, "CE": 0.20, "LP": 0.15, "LS": 0.10}
    sensitivity_map = {
        "full_name": 1.0, "passport_no": 1.0,
        "postcode": 0.65, "occupation": 0.45,
        "borough": 0.20, "year": 0.05,
    }

    # In production these come from your metadata registry and auxiliary-overlap engine.
    ce = pd.Series(0.6, index=gdf.index)  # Replace with registry lookup
    lp = pd.Series(0.5, index=gdf.index)  # Replace with spatial overlap analysis
    ls = pd.Series(0.8, index=gdf.index)  # Replace with ETL pipeline stage tag

    scores = calculate_spatial_risk_score(gdf, weights, sensitivity_map, ce, lp, ls)
    gdf["risk_score"] = scores
    gdf["risk_tier"] = apply_risk_tiers(scores)

    # Route: quarantine high/critical; publish low/medium after masking
    quarantine = gdf[gdf["risk_tier"].isin(["high", "critical"])]
    publish_candidates = gdf[gdf["risk_tier"].isin(["low", "medium"])]
    print(gdf[["risk_score", "risk_tier"]].describe())

Wrap the call to calculate_spatial_risk_score in a CI/CD validation step that asserts the score distribution against a rolling 90-day baseline. A sudden population shift toward higher scores indicates schema drift, a new high-sensitivity attribute, or a CRS mismatch.

Validation & Re-identification Testing permalink

A composite score is a model — validate it before trusting it to govern data releases.

Entropy checks. For each scoring tier, compute the Shannon entropy of the spatial distribution. Records scoring critical should cluster in a small area of the attribute space; records scoring low should span it broadly. If low-scored records show low entropy, the ASAS weight is understating identifier risk.

from scipy.stats import entropy

def tier_entropy_audit(gdf: gpd.GeoDataFrame, score_col: str = "risk_score") -> pd.Series:
    """
    Compute per-tier coordinate entropy as a validation signal.
    Low-scored records should show high spatial entropy (spread out);
    high-scored records may show low entropy (dense, re-identifiable clusters).
    """
    results = {}
    for tier in ["low", "medium", "high", "critical"]:
        subset = gdf[gdf["risk_tier"] == tier]
        if subset.empty:
            results[tier] = float("nan")
            continue
        # Round to 2 decimal places to create a discrete distribution.
        coords = subset.geometry.centroid
        lat_bins = np.round(coords.y, 2)
        counts = lat_bins.value_counts(normalize=True).values
        results[tier] = float(entropy(counts, base=2))
    return pd.Series(results)

Neighbor-count audits. For point data, verify that low-scored records have at least k=5k = 5 spatial neighbors within a configurable radius, confirming that the SGSG score correctly reflects sparse-area exposure. Records that score low but have fewer than kk neighbors suggest that AmaxA_{\max} is set too high for your dataset’s geographic extent.

Adversarial auxiliary-join simulation. Periodically execute a mock spatial linkage attack using publicly available datasets (property records, transit stop locations, business registers). Records that join successfully should correspond to high LPLP scores. Any successful join on a low-scored record indicates that LPLP weights or inputs need recalibration.

Score drift monitoring. Track the mean and standard deviation of risk_score across weekly releases. A shift of more than one standard deviation in the mean signals that a schema change or new attribute has entered the pipeline without a corresponding update to the sensitivity map.

Common Failure Modes & Gotchas permalink

  • Mixed-CRS ingestion. If two feature layers are joined before CRS normalization, bounding-box areas computed in EPSG:3857 metres will be orders of magnitude larger than their EPSG:4326 equivalents. The log normalization compresses them all to near-zero, producing uniformly low SGSG scores for a dataset that actually contains sub-metre GPS traces.

  • Static placeholders left in production. The implementation above assigns fixed values to ce_scores, lp_scores, and ls_scores for illustration. In production, all three must be populated from live metadata; leaving the defaults causes the framework to systematically underestimate exposure for data shared on public APIs (CE fixed at 0.6 when 0.9 is appropriate) and raw ingestion records (LS fixed at 0.8, which is already high, but the real value may be 1.0 for unprocessed personal records).

  • Polygon features with large bounding boxes. A municipality boundary polygon may have a bounding box of several hundred km², driving SGSG toward zero even though the feature’s actual area is tiny. For polygon layers, compute SGSG from the feature’s actual area rather than its bounding box when geometries are not axis-aligned or are highly irregular.

  • Sparse-data edge cases. In low-density rural datasets, even census-tract centroids may represent only one or two individuals. A low SGSG score (large area) can mask extreme individual exposure. Combine the scoring framework with the k-anonymity grouping for location traces check to flag geometries that span a large area but represent fewer than k=5k = 5 individuals.

  • Utility collapse from over-weighting. Setting wASw_{AS} above 0.5 in a dataset of medical records can push nearly all features into the critical tier regardless of spatial granularity, blocking legitimate aggregated-area releases. If more than 30% of records hit critical, audit the weight calibration — utility collapse often indicates a calibration error rather than genuinely critical data.

  • Weight drift between jurisdictions. When a pipeline crosses regional boundaries, CECE weights must adjust to reflect local regulatory baselines. A weight calibrated for an internal analytics environment becomes dangerously low if the same dataset is later shared with a third-party API consumer. Encode jurisdiction-specific weight presets in your data-governance configuration rather than hard-coding a single weight set.

Compliance Alignment permalink

Regulatory clause Satisfied by
GDPR Art. 35 (DPIA) Documented composite score serves as the quantitative privacy impact measure; threshold crossings are the trigger for formal review
GDPR Art. 25 (Data Protection by Design) Ingestion-gate integration ensures scoring runs before any disclosure decision
CCPA § 1798.100 (data minimization) critical-tier block prevents unnecessary processing of high-sensitivity spatial records
NIST Privacy Framework PR.DS-1 Sensitivity map and weight versioning satisfies data-at-rest classification controls
NIST Privacy Framework GV.PO-P2 Quarterly threshold review and changelog satisfy governance policy documentation requirements

Automated logs of threshold crossings, applied transformations, and human-review approvals constitute the audit trail. Store these logs alongside the version-controlled scoring configuration so that auditors can trace every disclosure decision back to the exact weight set and sensitivity map in effect at the time. The compliance mapping for GDPR and CCPA location data page maps these clauses to specific GIS pipeline controls in greater detail.

For the privacy risk matrix for municipal GIS use case, align scoring tiers with open-data publication schedules and communicate threshold rationale to policy makers in non-technical language backed by the quantitative score distributions.

Back to Spatial Privacy Fundamentals & Threat Modeling