Building a Competitor Proximity Index in PostGIS

This page solves one exact task: computing, for every candidate site, the distance to its nearest competitor and the count of competitors within a set of radii, then materializing those into an indexed proximity table that the scoring pipeline can join in milliseconds.

Doing this in Python with a nested loop is $O(n \times m)$ and collapses at portfolio scale. PostGIS answers the same questions with a GiST spatial index and two operators — ST_DWithin for radius counts and the <-> KNN operator for nearest-neighbor distance — turning an all-pairs scan into an index-assisted lookup. The output is a per-candidate proximity index table that feeds the competitive saturation term of the suitability score without any further geometry work at scoring time.

Prerequisites #

Before running this task you need:

• PostGIS installed and enabled (CREATE EXTENSION postgis;) on a PostgreSQL 13+ instance, set up per Setting Up PostGIS for Retail Analytics.
• Two point tables: candidate_sites and competitors, each with a geometry(Point, 5070) column. Storing geometry in the projected EPSG:5070 (Conus Albers) means ST_Distance and ST_DWithin return metres directly, avoiding the geography-vs-geometry pitfalls discussed below.
• Python packages for the driver: psycopg[binary] and, optionally, geopandas to read the result back. Install with pip install "psycopg[binary]" geopandas.
• The parent context. This page assumes familiarity with Competitor Mapping & Cannibalization Analysis, which defines the saturation index this proximity table feeds.

Parameters #

Storing points in a projected SRID rather than geography is a deliberate choice for a bounded national portfolio: planar distance in EPSG:5070 is accurate to well under a percent over retail trade-area scales and lets the <-> operator use the GiST index directly. geography distances are geodesically exact but the KNN ordering semantics differ, and mixing the two silently changes results.

Building the index and the proximity table #

The SQL below creates GiST indexes on both point layers, then materializes a proximity table. Nearest-competitor distance uses a LATERAL join with the <-> operator so the index returns candidates in distance order and LIMIT 1 stops at the first; the radius counts use ST_DWithin, which is index-assisted and far cheaper than computing ST_Distance for every pair and filtering.

-- 1. GiST indexes are what make ST_DWithin and <-> index-assisted.
CREATE INDEX IF NOT EXISTS idx_candidate_geom
    ON candidate_sites USING GIST (geom);
CREATE INDEX IF NOT EXISTS idx_competitor_geom
    ON competitors USING GIST (geom);
ANALYZE candidate_sites;
ANALYZE competitors;

-- 2. Materialize per-candidate proximity metrics in one pass.
DROP TABLE IF EXISTS competitor_proximity;
CREATE TABLE competitor_proximity AS
SELECT
    c.site_id,
    nn.dist_m               AS nearest_competitor_m,
    nn.brand                AS nearest_brand,
    (SELECT count(*) FROM competitors k
       WHERE ST_DWithin(k.geom, c.geom, 1000)) AS n_within_1km,
    (SELECT count(*) FROM competitors k
       WHERE ST_DWithin(k.geom, c.geom, 3000)) AS n_within_3km,
    (SELECT count(*) FROM competitors k
       WHERE ST_DWithin(k.geom, c.geom, 5000)) AS n_within_5km
FROM candidate_sites AS c
CROSS JOIN LATERAL (
    -- KNN: <-> uses the GiST index to return the closest competitor first.
    SELECT k.brand, k.geom <-> c.geom AS dist_m
    FROM competitors AS k
    ORDER BY k.geom <-> c.geom
    LIMIT 1
) AS nn;

-- 3. Index the output so the scoring join is fast too.
CREATE UNIQUE INDEX ON competitor_proximity (site_id);

The <-> operator inside an ORDER BY ... LIMIT is the idiomatic PostGIS KNN pattern: the planner walks the GiST index in ascending distance order and stops after one row, so nearest-neighbor is a bounded index traversal rather than a full scan. Both point layers carry SRID 5070, so dist_m and the ST_DWithin thresholds are all in metres with no unit conversion.

Driving the build and reading it back from Python #

The Python driver runs the DDL, then reads the proximity table into a GeoDataFrame for the scoring stage. It asserts the SRID it expects before trusting any distance, mirroring the CRS discipline enforced everywhere else in the pipeline.

import psycopg
import geopandas as gpd

EXPECTED_SRID = 5070

def build_proximity_index(dsn: str, ddl_sql: str) -> None:
    """Run the GiST + proximity DDL as a single transaction."""
    with psycopg.connect(dsn) as conn, conn.cursor() as cur:
        # Guard against a silent SRID mismatch before we compute distances.
        cur.execute("SELECT DISTINCT ST_SRID(geom) FROM candidate_sites;")
        srids = {row[0] for row in cur.fetchall()}
        assert srids == {EXPECTED_SRID}, f"candidate SRID {srids} != {EXPECTED_SRID}"
        cur.execute("SELECT DISTINCT ST_SRID(geom) FROM competitors;")
        assert {row[0] for row in cur.fetchall()} == {EXPECTED_SRID}, "competitor SRID mismatch"
        cur.execute(ddl_sql)     # the CREATE INDEX / CREATE TABLE block above
        conn.commit()

def read_proximity(dsn: str) -> gpd.GeoDataFrame:
    """Load the materialized proximity metrics joined back to site geometry."""
    sql = """
        SELECT c.site_id, c.geom,
               p.nearest_competitor_m, p.nearest_brand,
               p.n_within_1km, p.n_within_3km, p.n_within_5km
        FROM competitor_proximity p
        JOIN candidate_sites c USING (site_id)
    """
    gdf = gpd.read_postgis(sql, dsn, geom_col="geom")
    assert gdf.crs is not None and gdf.crs.to_epsg() == EXPECTED_SRID, "unexpected CRS on read"
    return gdf

Reading through read_postgis preserves the SRID from the geometry column metadata, and the assertion turns a missing or mismatched CRS into an immediate failure rather than a distorted saturation score three stages later.

Failure modes and debugging #

The geography-versus-geometry decision is the subtlest trap. geography gives geodesically correct distances on the sphere, but for a bounded regional portfolio a projected geometry in EPSG:5070 is both accurate enough and materially faster, and it makes the <-> KNN operator behave predictably. Pick one representation and enforce it with the SRID assertions above; do not mix them within a single proximity calculation.

Verification #

Confirm the index is actually used and the counts are sane before wiring the table into scoring:

1. Index usage: run EXPLAIN (ANALYZE, BUFFERS) on the KNN subquery and confirm an Index Scan using idx_competitor_geom rather than a Seq Scan. A sequential scan means the GiST index is missing or the planner’s statistics are stale.
2. Row completeness: the proximity table row count must equal the candidate count — SELECT count(*) FROM competitor_proximity should match candidate_sites. A shortfall means the KNN LATERAL found no competitor for some sites, usually an empty competitor table or a SRID mismatch.
3. Monotone radii: for every row, n_within_1km <= n_within_3km <= n_within_5km; a violation signals a units error in one of the ST_DWithin calls.
4. Nearest bounds the counts: if nearest_competitor_m > 5000 then n_within_5km must be 0. This cross-check catches a KNN result computed against the wrong table.

EXPLAIN (ANALYZE, BUFFERS)
SELECT k.brand, k.geom <-> (SELECT geom FROM candidate_sites LIMIT 1) AS d
FROM competitors k
ORDER BY k.geom <-> (SELECT geom FROM candidate_sites LIMIT 1)
LIMIT 1;

Seeing an Index Scan in the plan is the definitive confirmation that the GiST index is doing its job — the same tuning discipline that governs any large-scale spatial join in the platform. Once verified, the competitor_proximity table becomes a plain indexed join at scoring time, and the nearest-distance and radius-count columns drop straight into the competitive saturation index.

Related #

• Setting Up PostGIS for Retail Analytics — the database, SRID, and index conventions this task depends on.
• Location Intelligence Architecture & Data Foundations — where the spatial database sits in the wider platform.

← Back to Competitor Mapping & Cannibalization Analysis