Private Set Intersection for Spatial Joins

Two ride-hail companies want to know how many census tracts are visited by both their rider populations — a competitive-overlap figure that would let each size the other’s footprint without a data-sharing agreement. The catch is that the only way to answer it naively is to trade rider location histories, which neither legal team will ever sign off on and which would expose precise home and work coordinates the instant the files change hands. Private set intersection (PSI) is the cryptographic answer: each party discretizes its riders’ locations into cell tokens, blinds every token under a secret key, and exchanges only the blinded values, so the two sides can compute the size of their shared tract set while the raw coordinate lists — and even the plaintext cell lists — never cross the boundary between them.

This guide implements the Diffie–Hellman / oblivious-PRF construction of PSI adapted to spatial joins, the canonical “private set intersection, joint proximity” workload named in the Secure Multi-Party Computation in Spatial Analytics section. It covers reducing continuous coordinates to comparable geohash tokens, the double-exponentiation trick that makes equal cells collide only after both keys are applied, and returning intersection cardinality rather than the intersection itself. The equality-matching variant here answers “same cell”; the “near, not equal” question — matching tracts within a radius rather than by exact cell — is the fuzzy/threshold extension covered by the child guide on spatial PSI for proximity queries.

Prerequisites & Cryptographic Environment

The reference implementation runs on the Python standard library alone: hashlib for the full-domain hash that maps a cell token into the group, secrets for cryptographically secure key generation and shuffling, and the three-argument pow() built-in for modular exponentiation. Production deployments add the cryptography library for the mutually authenticated TLS 1.3 transport that carries blinded elements between parties, and a KMS or HSM binding for the per-party secret exponents. A geohash or H3 encoder is the only spatial dependency; the code below implements geohash inline so nothing external is required to reproduce it.

Three invariants must hold before the first token is blinded:

• A canonical, agreed cell grid. Both parties must discretize with identical parameters — same scheme (geohash vs H3), same precision, same CRS (EPSG:4326 decimal degrees). A one-level mismatch in geohash precision means no two cells ever collide and the intersection is silently, wrongly zero. This is the front-end contract that the coordinate masking protocols normalise before geometry reaches this layer.
• Shared group parameters. Both parties fix the same prime modulus and the same hash-to-group mapping out of band. Correctness of the match relies only on a shared modulus and consistent exponentiation; security additionally relies on the Decisional Diffie–Hellman (DDH) assumption in the subgroup the tokens are hashed into.
• A declared leakage budget. PSI-cardinality still reveals the size of the overlap and each party’s set size. Decide up front what that discloses about your population and bound it — the privacy model comparison frames when PSI’s exact-count leakage is acceptable versus when the count itself must be differentially-private noised before release.

For the modulus we use the Mersenne prime $p = 2^{521} - 1$: a large, well-documented value that keeps every group element far below the field boundary. Where PSI feeds a downstream computation on the matched cells rather than a bare count, homomorphic encryption basics and secret sharing for coordinates are the neighbouring tools in the same pipeline.

Step 1: Discretize Coordinates to Cell Tokens

PSI is a set operation over equal elements, and two GPS fixes are essentially never bit-for-bit equal. The reduction that makes spatial joins tractable is discretization: snap every coordinate to the cell of a fixed grid, and treat the cell identifier as the set element. Two riders in the same census tract now produce the same token, so “same cell” becomes the equality the protocol tests. Geohash precision is the privacy–utility dial — geohash-6 is roughly a $1.2\,\text{km} \times 0.6\,\text{km}$ cell (census-tract scale), geohash-7 tightens to about $153\,\text{m} \times 153\,\text{m}$. Coarser cells raise the population per cell (stronger k-anonymity, blurrier join); finer cells sharpen the join but shrink the token domain toward brute-forceable size.

from __future__ import annotations

from typing import Iterable, List, Set

_BASE32: str = "0123456789bcdefghjkmnpqrstuvwxyz"


def geohash_encode(lat: float, lon: float, precision: int = 6) -> str:
    """Encode an EPSG:4326 (lat, lon) pair to a geohash cell token.

    The cell token is the ONLY thing derived from a raw coordinate that ever
    enters the protocol. Precision fixes the cell size and therefore both the
    join resolution and the k-anonymity floor of every matched cell.
    """
    if not (-90.0 <= lat <= 90.0) or not (-180.0 <= lon <= 180.0):
        raise ValueError(f"Coordinate out of WGS84 range: ({lat}, {lon})")
    lat_lo, lat_hi = -90.0, 90.0
    lon_lo, lon_hi = -180.0, 180.0
    bits = [16, 8, 4, 2, 1]
    out: List[str] = []
    bit, ch, even = 0, 0, True
    while len(out) < precision:
        if even:                              # longitude bit
            mid = (lon_lo + lon_hi) / 2
            if lon > mid:
                ch |= bits[bit]; lon_lo = mid
            else:
                lon_hi = mid
        else:                                 # latitude bit
            mid = (lat_lo + lat_hi) / 2
            if lat > mid:
                ch |= bits[bit]; lat_lo = mid
            else:
                lat_hi = mid
        even = not even
        if bit < 4:
            bit += 1
        else:
            out.append(_BASE32[ch]); bit, ch = 0, 0
    return "".join(out)


def discretize(points: Iterable[tuple[float, float]], precision: int = 6) -> Set[str]:
    """Reduce a party's raw location list to a SET of cell tokens.

    Deduplication is intrinsic: PSI is a set operation, so a tract a rider
    visits ten times contributes one token, not ten — and the raw visit
    coordinates are discarded the moment the token set is built.
    """
    return {geohash_encode(lat, lon, precision) for lat, lon in points}

Building a set here is the first privacy decision, not just a convenience: collapsing repeat visits to one token per cell prevents the overlap count from being inflated by frequency, and it discards the raw coordinate the instant the token exists. Everything downstream operates on cell tokens only.

Step 2: Hash Cell Tokens into the Group

A cell token like "dr5ru7" cannot be exponentiated directly; it must first become a group element. A full-domain hash maps the token into $\mathbb{F}_p$, and squaring the result lands it in the subgroup of quadratic residues, which strips the low-order component a small-subgroup attack would exploit. This hash-to-group step is what turns DH-PSI into an oblivious pseudo-random function (OPRF): the blinded value $H(\text{cell})^k$ is a PRF of the cell under key $k$, computable by the key-holder without ever learning the input.

import hashlib

# 2^521 - 1 is a Mersenne prime — a large, well-documented modulus. Match
# correctness needs only a shared modulus and consistent exponentiation;
# security relies on DDH in the quadratic-residue subgroup we hash into.
GROUP_PRIME: int = (1 << 521) - 1


def hash_to_group(token: str, prime: int = GROUP_PRIME) -> int:
    """Map a cell token to a quadratic residue in F_p via a full-domain hash.

    Squaring forces the element into the QR subgroup, removing the low-order
    part that would otherwise leak a bit under a small-subgroup attack. The
    mapping is deterministic, so identical cells hash identically across parties
    — the property the whole intersection depends on.
    """
    digest = hashlib.sha512(token.encode("utf-8")).digest()
    h = int.from_bytes(digest, "big") % prime
    if h == 0:                                # avoid the degenerate zero element
        h = 1
    return pow(h, 2, prime)                   # into the quadratic-residue subgroup

Determinism is the load-bearing property: both parties must derive the same integer from the same cell, so the hash and the modulus are part of the out-of-band agreement, never chosen per-party.

Step 3: Party A Blinds with Key a

Each party is one instance of SpatialPSI holding a secret exponent and a token set. Party A raises every hashed cell to its private key $a$, producing $H(x)^a$ for each cell $x$, then shuffles the results before sending them. Shuffling severs the positional link between an input cell and its blinded output, so Party B — who will operate on these values next — cannot map a returned element back to an ordering or a rider. The secret key never leaves the party.

import secrets


class SpatialPSI:
    """One party in a Diffie–Hellman PSI over discretized spatial cells.

    A party holds a set of cell tokens and a secret blinding key. It exchanges
    only blinded group elements; the raw cell set — and therefore the raw
    coordinates behind it — never crosses the trust boundary. Equal cells
    collide only after BOTH parties' keys have been applied.
    """

    def __init__(self, cells: Iterable[str], prime: int = GROUP_PRIME) -> None:
        self.prime: int = prime
        # Set semantics: duplicates would inflate the cardinality the peer learns.
        self.cells: List[str] = sorted(set(cells))
        # Secret blinding key in [2, p-2]; transmitted to no one, ever.
        self.key: int = secrets.randbelow(prime - 3) + 2

    def blind_own(self) -> List[int]:
        """Round 1: H(cell)^key for each owned cell, returned SHUFFLED.

        Shuffling breaks the positional correlation between input cells and
        emitted elements, so the peer cannot re-link a value to a rider.
        """
        blinded = [pow(hash_to_group(c, self.prime), self.key, self.prime)
                   for c in self.cells]
        secrets.SystemRandom().shuffle(blinded)   # unlink order from identity
        return blinded

Because the key is drawn independently on each side and never exchanged, no single party ever holds both exponents — the structural property that keeps either side from unblinding the other’s set alone.

Step 4: Party B Double-Blinds with Key b

Party B receives Party A’s set of $H(x)^a$ values and raises each to its own key $b$, yielding $H(x)^{ab}$ — the double-blinded token. B performs the identical blind_own step on its own cells to emit $H(y)^b$, which A then re-blinds to $H(y)^{ab}$. The commutativity of exponentiation is the whole trick: $\left(H(c)^a\right)^b = \left(H(c)^b\right)^a = H(c)^{ab}$, so a cell present in both sets lands on the same integer regardless of which party applied which key first. A cell in only one set has no counterpart to collide with.

    def reblind(self, elements: Iterable[int]) -> List[int]:
        """Round 2: raise the PEER's already-blinded elements to our key.

        Applied to H(peer_cell)^peer_key this yields H(peer_cell)^{ab}: the
        double-blinded token that is directly comparable across parties. We
        learn nothing about the peer's cells — the inputs are already blinded
        under a key we do not hold.
        """
        return [pow(e, self.key, self.prime) for e in elements]

Note that reblind operates on integers that are already blinded under the peer’s secret key, so this step reveals nothing about the peer’s cells to the party running it. Each side contributes exactly one exponentiation of the other’s data and never sees an unblinded token.

Step 5: Compare Double-Blinded Tokens

Once both sides hold their two double-blinded sets, the intersection is a plain set intersection over integers. Equal underlying cells produced equal $H(c)^{ab}$ values; unequal cells produced values that are, under DDH, computationally indistinguishable from random and will not collide except with negligible probability. Whichever party is designated to compute the result takes the two double-blinded sets and intersects them.

def matched_tokens(double_blinded_a: Iterable[int],
                   double_blinded_b: Iterable[int]) -> Set[int]:
    """Set of double-blinded integers common to both parties.

    Membership here corresponds exactly to a shared cell, but the returned
    integers are H(cell)^{ab} — opaque tokens, not cells. Recovering which
    tracts matched would require inverting the OPRF, which neither party can do
    without the other's key.
    """
    return set(double_blinded_a) & set(double_blinded_b)

The matched tokens are opaque: they identify that a cell overlaps, not which cell. To turn them back into tracts, a party would need to invert the OPRF, which is infeasible without the counterpart key — so even the party computing the intersection learns only anonymous collision points unless it deliberately reveals its own blinded-to-cell map.

Step 6: Return Cardinality Only

For the mobility-overlap question, the honest minimum disclosure is the count of shared tracts — PSI-cardinality — not the set of tracts. Returning the intersection itself would tell each party exactly which census tracts its competitor’s riders also visit; returning only $|A \cap B|$ answers “how much overlap” while withholding “where.” The cardinality function is a one-line specialization of the compare step, and it is the value you release (optionally after adding differentially-private noise, per the leakage budget from the prerequisites).

def intersection_cardinality(double_blinded_a: Iterable[int],
                             double_blinded_b: Iterable[int]) -> int:
    """Return ONLY |A ∩ B| — the count of shared cells, never their identities.

    This is the PSI-cardinality variant: it answers 'how many tracts overlap'
    while withholding 'which tracts', the minimal disclosure for a competitive
    footprint estimate. Noise this count before release if the exact size is
    itself sensitive.
    """
    return len(matched_tokens(double_blinded_a, double_blinded_b))

Choosing cardinality over the full intersection is a design decision, not a default. If a downstream workflow genuinely needs the matched cells — say, to run a further encrypted aggregation over them — keep the tokens opaque and feed them into that computation rather than decoding them to tracts on either side.

Step 7: Validate Correctness and Non-Crossing of Raw Sets

Two properties must hold and both must be tested in CI: the protocol computes the true cardinality (no false matches, no missed matches), and no raw cell token ever appears in anything a party transmits. The harness below drives a full two-party run, checks the computed cardinality against the ground-truth set intersection over known geohash cells, confirms that key order does not change the result, and asserts that every transmitted value is a blinded integer distinct from any plaintext-derived hash.

def _validate() -> None:
    # Two providers with overlapping and disjoint coverage (lat, lon points).
    a_points = [(40.7484, -73.9857), (40.7580, -73.9855),   # Midtown, shared
                (40.6892, -74.0445)]                          # Liberty Is., A-only
    b_points = [(40.7484, -73.9857), (40.7580, -73.9855),   # Midtown, shared
                (34.0522, -118.2437)]                         # LA, B-only

    a_cells = discretize(a_points, precision=6)
    b_cells = discretize(b_points, precision=6)
    ground_truth = len(a_cells & b_cells)   # plaintext oracle for the test only

    party_a = SpatialPSI(a_cells)
    party_b = SpatialPSI(b_cells)

    # Round 1: each party blinds and ships only shuffled group elements.
    a_blinded = party_a.blind_own()          # H(x)^a  — crosses the boundary
    b_blinded = party_b.blind_own()          # H(y)^b  — crosses the boundary

    # Round 2: each re-blinds the OTHER's set to reach H(.)^{ab}.
    ab_from_a = party_b.reblind(a_blinded)   # H(x)^{ab}
    ab_from_b = party_a.reblind(b_blinded)   # H(y)^{ab}

    card = intersection_cardinality(ab_from_a, ab_from_b)

    # 1. Correctness: PSI cardinality equals the plaintext intersection size.
    assert card == ground_truth == 2, f"cardinality mismatch: {card} != {ground_truth}"

    # 2. Commutativity: swapping which party re-blinds first changes nothing.
    swapped = intersection_cardinality(party_a.reblind(b_blinded),
                                       party_b.reblind(a_blinded))
    assert swapped == card

    # 3. Raw sets never cross: no transmitted value equals a plaintext-cell hash.
    plain_hashes = {hash_to_group(c) for c in a_cells | b_cells}
    on_the_wire = set(a_blinded) | set(b_blinded)
    assert plain_hashes.isdisjoint(on_the_wire), "unblinded token leaked onto the wire"

    # 4. Blinding is non-trivial: keys actually transformed the elements.
    assert all(v not in plain_hashes for v in a_blinded)
    print("all spatial-PSI invariants hold — cardinality correct, no raw leakage")


if __name__ == "__main__":
    _validate()

Invariant 3 is the one auditors care about most: it mechanically proves that the values leaving each party are blinded under a secret key and cannot be matched against a plaintext-cell dictionary — the on-wire form is disjoint from every hash of a known cell.

Threat Model Considerations

DH-PSI changes the adversary’s job from “read the other party’s locations” to “break DDH or exploit what the protocol legitimately leaks.” The surfaces specific to a spatial join:

• Semi-honest vs malicious parties. The construction above is secure against semi-honest (honest-but-curious) adversaries who follow the protocol but analyse what they see. A malicious party can deviate — most dangerously by padding its input set with a fine-grained sweep of candidate cells to probe whether specific tracts are in the peer’s set. Defences are set-size commitments, a proof that inputs were honestly derived, or moving to a maliciously-secure PSI variant; assume semi-honest only if a contract or audit binds both sides to it.
• Set-size leakage. Even ideal PSI reveals each party’s set cardinality and the intersection size. For a spatial join, the number of distinct tracts a provider covers is itself competitively sensitive. Pad both sets to a common size with dummy cells drawn from outside the real domain, or release only a differentially-private noised cardinality.
• Small-domain brute force. The token domain is the number of possible cells, which for high geohash precision over a metro area can be small enough to enumerate. An adversary who knows the OPRF key (or colludes to obtain a doubly-keyed oracle) can hash every candidate cell and match it against the blinded set, recovering the peer’s tracts. Keep keys strictly single-party, coarsen the cell grid so each cell holds a defensible k-anonymity population, and never expose a blind-then-return oracle that lets a peer blind arbitrary probe cells.
• Cell-alignment side channel. If one party can choose the grid origin or precision after seeing hints about the other’s data, it can craft cells that isolate a target. Fix the grid parameters out of band and identically for both sides, before any exchange.
• Ordering and timing correlation. Failing to shuffle blinded outputs, or emitting them in input order, lets a peer correlate positions across rounds. The SystemRandom().shuffle in blind_own closes this; batch the exponentiations so wall-clock timing does not correlate with set content.

Validation & Compliance Checklist

Wire each control to a measurable pass/fail gate rather than a code review:

1. Cardinality correctness — PASS if the protocol’s intersection_cardinality equals the plaintext set-intersection size for a randomized $\ge 10^4$-cell property test across overlapping and disjoint coverage. Fail the build on any mismatch.
2. Commutativity — PASS if swapping which party applies its key first leaves the cardinality unchanged for every test case; a difference means a key or modulus is not shared correctly.
3. No raw leakage — PASS if the set of on-the-wire values is provably disjoint from the hash of every known cell token (invariant 3 above), asserted on every run.
4. Grid-agreement — PASS if both parties’ discretization is verified identical (scheme, precision, CRS) via a shared config hash; a mismatch must abort before any exchange, not silently return zero overlap.
5. Set-size disclosure bound — PASS if transmitted set sizes are padded to a common ceiling or the released cardinality is DP-noised, with the chosen $\varepsilon$ tied to a control in the compliance framework mapping. Record the disclosure per run.
6. Key isolation — PASS if each secret exponent is generated and held on exactly one party (KMS/HSM audit), and no endpoint offers a blind-arbitrary-input oracle that would enable small-domain brute force.

The _validate() harness in Step 7 exercises criteria 1, 2, and 3 directly and is safe to run in CI.

Failure Modes & Remediation

Spatial PSI fails quietly — usually by returning a plausible wrong count rather than an error.

• Precision mismatch → silent zero overlap. One party discretizes at geohash-6, the other at geohash-7; no tokens ever collide and the intersection is always zero. Detection: compare a config hash of the grid parameters before exchange. Recovery: fix precision, CRS, and scheme in a shared, versioned config and abort on any drift.
• Modulus or hash drift → garbage matches. A build where one side uses a different prime or hash function makes double-blinded tokens incomparable, producing a near-zero or nonsensical count. Detection: the commutativity check (criterion 2) fails. Recovery: pin group parameters in the shared agreement and assert them at handshake.
• Duplicate-inflated cardinality. Feeding a list instead of a set lets repeat visits multiply matches and overstate the overlap. Detection: assert set semantics on ingest. Recovery: route all input through discretize, which deduplicates by construction.
• Set-size fingerprinting. A party with a very small coverage set reveals nearly its whole footprint through the size and intersection alone. Detection: audit set sizes against a minimum-padding policy. Recovery: pad to a common ceiling with out-of-domain dummy cells, and DP-noise the released count.
• Malicious over-fine probing. A counterparty submits a dense sweep of candidate cells to test membership of specific tracts. Detection: input-size commitments and grid-agreement checks. Recovery: enforce a maximum set size, require honest-derivation proofs, or escalate to a maliciously-secure PSI protocol; move the exchange onto a partition-tolerant transport via async routing for MPC so a stalled round cannot be exploited to replay probes.

Frequently Asked Questions

Related

This guide is part of the Secure Multi-Party Computation in Spatial Analytics reference — see it for how PSI, masking, sharing, and encryption compose into one pipeline.

• Spatial PSI for Proximity Queries — the fuzzy/threshold extension that matches “near” rather than “equal”, for radius-based joins.
• Coordinate Masking Protocols — the deterministic front end that normalises geometry and CRS before discretization.
• Secret Sharing for Coordinates — the distribution layer for when matched cells feed a joint computation rather than a bare count.
• Homomorphic Encryption Basics — the compute-on-ciphertext alternative when you must aggregate over the intersection.
• Async Routing for MPC — partition-tolerant transport for the blinded-element exchange rounds.
• Privacy Model Comparison — when PSI’s exact-count leakage is acceptable versus when the count must be DP-noised.

Up: Secure Multi-Party Computation in Spatial Analytics