How to structure a geospatial database for multi-state retail chains

CREATE TABLE retail_sites (
    site_id       UUID PRIMARY KEY DEFAULT gen_random_uuid(),
    state_code    CHAR(2) NOT NULL,
    store_name    VARCHAR(150),
    lease_status  VARCHAR(20) CHECK (lease_status IN ('active', 'pending', 'closed', 'under_review')),
    geom          GEOMETRY(Point, 4326) NOT NULL,
    geom_valid    BOOLEAN GENERATED ALWAYS AS (ST_IsValid(geom)) STORED,
    created_at    TIMESTAMPTZ DEFAULT NOW(),
    updated_at    TIMESTAMPTZ DEFAULT NOW(),
    CONSTRAINT enforce_srid       CHECK (ST_SRID(geom) = 4326),
    CONSTRAINT enforce_valid_geom CHECK (ST_IsValid(geom))
);

-- Partitioned table (PostgreSQL 11+)
CREATE TABLE retail_sites (
    site_id      UUID          NOT NULL,
    state_code   CHAR(2)       NOT NULL,
    store_name   VARCHAR(150),
    lease_status VARCHAR(20),
    geom         GEOMETRY(Point, 4326) NOT NULL,
    created_at   TIMESTAMPTZ DEFAULT NOW(),
    updated_at   TIMESTAMPTZ DEFAULT NOW(),
    PRIMARY KEY (site_id, state_code)   -- partition key must be part of PK
) PARTITION BY LIST (state_code);

CREATE TABLE retail_sites_ca      PARTITION OF retail_sites FOR VALUES IN ('CA');
CREATE TABLE retail_sites_tx      PARTITION OF retail_sites FOR VALUES IN ('TX');
CREATE TABLE retail_sites_ny      PARTITION OF retail_sites FOR VALUES IN ('NY');
CREATE TABLE retail_sites_default PARTITION OF retail_sites DEFAULT;

-- Spatial index per partition
CREATE INDEX CONCURRENTLY idx_retail_sites_ca_geom
    ON retail_sites_ca USING GIST (geom);
CREATE INDEX CONCURRENTLY idx_retail_sites_tx_geom
    ON retail_sites_tx USING GIST (geom);
CREATE INDEX CONCURRENTLY idx_retail_sites_ny_geom
    ON retail_sites_ny USING GIST (geom);

-- Partial index for active-lease proximity queries
CREATE INDEX CONCURRENTLY idx_retail_sites_tx_active_geom
    ON retail_sites_tx USING GIST (geom)
    WHERE lease_status = 'active';

import logging
import geopandas as gpd
import psycopg2
from psycopg2.extras import execute_values
from shapely.geometry import Point
from shapely.validation import make_valid

logging.basicConfig(level=logging.INFO)


def validate_and_ingest(
    gdf: gpd.GeoDataFrame,
    states: gpd.GeoDataFrame,   # boundary layer with a 'state_code' column
    conn_str: str,
) -> int:
    # 1. Assert the CRS explicitly; never assume an inbound layer is WGS84.
    if gdf.crs is None:
        gdf = gdf.set_crs(epsg=4326)
    gdf = gdf.to_crs(epsg=4326)
    states = states.to_crs(epsg=4326)

    # 2. Drop null geometries and repair any invalid topology.
    gdf = gdf.dropna(subset=["geometry"]).copy()
    gdf["geometry"] = gdf["geometry"].apply(make_valid)

    # 3. Keep only valid points (make_valid can return non-point types).
    valid_mask = gdf["geometry"].apply(lambda g: isinstance(g, Point) and g.is_valid)
    gdf = gdf[valid_mask].copy()

    # 4. Derive state_code by spatial containment, not a hard-coded bbox.
    gdf = gpd.sjoin(gdf, states[["state_code", "geometry"]],
                    how="left", predicate="within")
    gdf["state_code"] = gdf["state_code"].fillna("ZZ")  # routes to DEFAULT partition

    # 5. Bulk insert using EWKT geometry literals.
    insert_sql = """
        INSERT INTO retail_sites (site_id, state_code, store_name, lease_status, geom)
        VALUES %s
        ON CONFLICT (site_id, state_code) DO NOTHING
    """
    records = [
        (
            str(row["site_id"]),
            row["state_code"],
            row["store_name"],
            row.get("lease_status", "pending"),
            f"SRID=4326;POINT({row['geometry'].x} {row['geometry'].y})",
        )
        for _, row in gdf.iterrows()
    ]

    with psycopg2.connect(conn_str) as conn:
        with conn.cursor() as cur:
            execute_values(cur, insert_sql, records, page_size=5000)
    logging.info("Ingested %d validated records.", len(records))
    return len(records)

-- 1. Row count and how many landed in DEFAULT (should be near zero).
SELECT tableoid::regclass AS partition, count(*)
FROM retail_sites
GROUP BY 1 ORDER BY 2 DESC;

-- 2. No invalid geometries survived the constraints.
SELECT count(*) FROM retail_sites WHERE NOT ST_IsValid(geom);  -- expect 0

-- 3. Bounding-box sanity: every point inside the continental envelope.
SELECT count(*) FROM retail_sites
WHERE NOT ST_Within(geom, ST_MakeEnvelope(-125, 24, -66, 50, 4326));  -- expect 0

-- 4. Pruning works: this should scan only the CA partition.
EXPLAIN SELECT site_id FROM retail_sites WHERE state_code = 'CA';

-- 5. Refresh planner statistics after the bulk load.
ANALYZE retail_sites;