Designing Spatial Indexes for Geocoded Data in PostGIS

A geocoding pipeline only earns its keep once the coordinates it produces are queryable at speed. This topic sits inside the caching, deduplication and spatial indexing section and focuses on one job: designing spatial indexes so that proximity and bounding-box queries over geocoded point data stay fast as the table grows from thousands to hundreds of millions of rows. A spatial index is the data structure Postgres consults to avoid scanning every point when you ask “which addresses fall inside this box” or “what are the ten nearest stores to this customer.”

Get the index type, column type, and query shape right and a nearest-neighbour lookup returns in single-digit milliseconds against a table of 100M points. Get them wrong and the same query does a sequential scan and takes minutes. The sections below walk from column design through index selection, query patterns, edge cases, and maintenance, with runnable SQL and Python throughout.

Prerequisites

Before building indexes, confirm the foundations are in place:

If your points still live as raw latitude/longitude floats, promote them to a typed column first — every index and operator below depends on it.

Step 1 — Choose geometry or geography

PostGIS offers two spatial column types and the choice shapes everything downstream.

geometry treats coordinates as points on a flat Cartesian plane. Distance is computed in the units of the SRID (degrees for 4326), so a raw ST_Distance between two geometry points in 4326 is a meaningless “degrees” number — but bounding-box operations, index support, and function coverage are complete and fast. This is the right default for geocoded address points.

geography treats coordinates as points on a spheroid and returns distances in metres from ST_Distance and ST_DWithin. It is convenient for global metric queries but every comparison is more expensive, and fewer functions accept it.

The pragmatic pattern most pipelines use: store geometry(Point, 4326), and when you need a true metric radius, cast to geography inside the query (ST_DWithin(geom::geography, ...)) or pre-project to a metre-based SRID. Keep one type per column and never mix them.

Step 2 — Create the table and load points

Define the table with an explicit point type and SRID so PostGIS can validate every insert:

CREATE TABLE geocoded_address (
    id           bigserial PRIMARY KEY,
    raw_input    text        NOT NULL,
    postal_code  text,
    geom         geometry(Point, 4326) NOT NULL,
    accuracy     text,        -- 'rooftop' | 'interpolated' | 'centroid'
    created_at   timestamptz DEFAULT now()
);

-- Insert a geocoded point (note ST_MakePoint takes lon, lat — X, Y order)
INSERT INTO geocoded_address (raw_input, postal_code, geom, accuracy)
VALUES (
    '1600 Amphitheatre Pkwy, Mountain View, CA',
    '94043',
    ST_SetSRID(ST_MakePoint(-122.084, 37.4220), 4326),
    'rooftop'
);

The single most common loading bug is coordinate-order inversion: ST_MakePoint takes (x, y) which is (longitude, latitude), the opposite of how humans and most APIs write “lat, lon.” The deduplication guidance for fuzzy canonical keys pairs well here: normalize and deduplicate before you index, so you index each distinct location once rather than storing thousands of near-identical rows.

Step 3 — Build the matching index type

PostGIS supports three access methods for spatial columns, each with a distinct trade-off.

GiST (Generalized Search Tree) is a balanced tree of bounding boxes. It is the general-purpose default: it handles arbitrarily distributed points, supports ST_DWithin, bounding-box overlap, and KNN ordering, and its performance is predictable regardless of insert order.

CREATE INDEX idx_geocoded_geom_gist
    ON geocoded_address USING gist (geom);

SP-GiST (Space-Partitioned GiST) uses a quadtree/k-d-tree partitioning. For point data with highly skewed clustering it can be smaller and slightly faster to probe than GiST, but it lacks some operator support and is a niche optimization.

CREATE INDEX idx_geocoded_geom_spgist
    ON geocoded_address USING spgist (geom);

BRIN (Block Range Index) stores only the bounding box of each range of physical table blocks. It is astonishingly small and fast to build, but it only helps when physical row order correlates with spatial location — for example when you bulk-load points sorted by geohash or postal code.

CREATE INDEX idx_geocoded_geom_brin
    ON geocoded_address USING brin (geom) WITH (pages_per_range = 32);

The trade-off table below summarizes when each access method wins. For a deeper decision walkthrough with a benchmark harness, see choosing GiST vs BRIN in PostGIS.

Index type reference

Access method Index size Build time Best for Update-friendly KNN (nearest-neighbour) support
GiST Large (~ table size) Slow–moderate Arbitrary point distributions, mixed read/write Yes Yes
SP-GiST Moderate Moderate Heavily skewed / clustered points Yes Yes
BRIN Tiny (kilobytes) Very fast Bulk-loaded data physically sorted by location Poor (correlation decays on random insert) No

The mental model: GiST indexes the values, BRIN indexes the physical layout. GiST is the safe default for a geocoding table that receives ongoing inserts; BRIN is a specialist tool for append-only, spatially sorted archives where index size matters.

GiST bounding-box tree over geocoded points A three-level GiST structure. The root node is a bounding box covering the whole extent. It contains two internal bounding boxes, each of which contains several leaf point entries. Arrows connect parent boxes to their children, illustrating how a query descends only into boxes that overlap the search region. Root bounding box covers full extent Internal box A west cluster Internal box B east cluster leaf points (geocoded addresses)

Step 4 — Write proximity and KNN queries

Two query shapes dominate geocoding workloads: radius filters and nearest-neighbour ranking. Both must be written so the planner uses the spatial index.

Radius / “within distance” — ST_DWithin. This is the correct operator for “all points within N metres.” It is index-assisted because PostGIS first tests the bounding-box expansion against the index, then refines. Cast to geography for a true metre radius:

-- All addresses within 500 metres of a query point
SELECT id, raw_input, ST_Distance(geom::geography, q.g) AS metres
FROM geocoded_address,
     (SELECT ST_SetSRID(ST_MakePoint(-122.084, 37.4220), 4326)::geography AS g) q
WHERE ST_DWithin(geom::geography, q.g, 500)
ORDER BY metres;

Nearest neighbour — the <-> KNN operator. To get the k closest points, order by the index-backed distance operator. The GiST index can return rows in distance order without scanning the whole table:

-- The 10 nearest addresses to a query point (KNN, index-ordered)
SELECT id, raw_input,
       geom <-> ST_SetSRID(ST_MakePoint(-122.084, 37.4220), 4326) AS dist_deg
FROM geocoded_address
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(-122.084, 37.4220), 4326)
LIMIT 10;

The <-> ordering distance is in SRID units (degrees for 4326), which is fine for ranking but not for reporting real distances — recompute the metric distance with a geography cast on just the returned rows. Confirm the plan with EXPLAIN (ANALYZE, BUFFERS) and look for Index Scan using idx_geocoded_geom_gist with an Order By: (geom <-> ...) line.

Step 5 — Python loader with index creation and KNN

The loader below uses psycopg 3 to create the table, build the GiST index, bulk-insert points, and run a nearest-neighbour query. It compiles the SQL once and parameterizes coordinates to avoid injection.

from __future__ import annotations

import logging
from typing import Iterable, Sequence

import psycopg

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger(__name__)

DDL = """
CREATE TABLE IF NOT EXISTS geocoded_address (
    id          bigserial PRIMARY KEY,
    raw_input   text NOT NULL,
    geom        geometry(Point, 4326) NOT NULL,
    accuracy    text
);
"""

CREATE_INDEX = """
CREATE INDEX IF NOT EXISTS idx_geocoded_geom_gist
    ON geocoded_address USING gist (geom);
"""

INSERT_POINT = """
INSERT INTO geocoded_address (raw_input, geom, accuracy)
VALUES (%(raw)s, ST_SetSRID(ST_MakePoint(%(lon)s, %(lat)s), 4326), %(acc)s);
"""

KNN_QUERY = """
SELECT id, raw_input,
       geom <-> ST_SetSRID(ST_MakePoint(%(lon)s, %(lat)s), 4326) AS dist_deg
FROM geocoded_address
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(%(lon)s, %(lat)s), 4326)
LIMIT %(k)s;
"""


def setup_schema(conn: psycopg.Connection) -> None:
    """Create the table and GiST index if they do not already exist."""
    with conn.cursor() as cur:
        cur.execute(DDL)
        cur.execute(CREATE_INDEX)
    conn.commit()
    logger.info("Schema and GiST index ready")


def load_points(
    conn: psycopg.Connection,
    rows: Iterable[dict[str, object]],
) -> int:
    """Bulk-insert geocoded points.

    Each row must contain 'raw', 'lon', 'lat', and 'acc' keys. Note the
    lon/lat order — ST_MakePoint expects (x, y) = (longitude, latitude).
    """
    count = 0
    with conn.cursor() as cur:
        for row in rows:
            cur.execute(INSERT_POINT, row)
            count += 1
    conn.commit()
    logger.info("Inserted %d points", count)
    return count


def nearest(
    conn: psycopg.Connection,
    lon: float,
    lat: float,
    k: int = 10,
) -> Sequence[tuple]:
    """Return the k nearest addresses to (lon, lat) using the KNN operator."""
    if not (-180.0 <= lon <= 180.0 and -90.0 <= lat <= 90.0):
        raise ValueError(f"coordinate out of range: lon={lon}, lat={lat}")
    with conn.cursor() as cur:
        cur.execute(KNN_QUERY, {"lon": lon, "lat": lat, "k": k})
        return cur.fetchall()


if __name__ == "__main__":
    with psycopg.connect("postgresql:///geocoding") as conn:
        setup_schema(conn)
        load_points(conn, [
            {"raw": "1600 Amphitheatre Pkwy", "lon": -122.084, "lat": 37.4220, "acc": "rooftop"},
            {"raw": "1 Hacker Way, Menlo Park", "lon": -122.148, "lat": 37.4848, "acc": "rooftop"},
        ])
        for row in nearest(conn, -122.10, 37.44, k=5):
            logger.info("match: %s", row)

Vectorized geopandas variant

When points arrive as a DataFrame, geopandas writes them to PostGIS in one call and lets you build the index in the same session. This is dramatically faster than row-by-row inserts for large loads.

from __future__ import annotations

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
from sqlalchemy import create_engine, text


def load_dataframe_to_postgis(
    df: pd.DataFrame,
    table: str = "geocoded_address",
    dsn: str = "postgresql:///geocoding",
) -> None:
    """Vectorized load of a lon/lat DataFrame into PostGIS with a GiST index.

    Expects columns 'lon', 'lat', and 'raw_input'. Builds geometry from the
    coordinate columns, writes in bulk, then creates the spatial index.
    """
    geometry = gpd.points_from_xy(df["lon"], df["lat"], crs="EPSG:4326")
    gdf = gpd.GeoDataFrame(df.drop(columns=["lon", "lat"]), geometry=geometry)

    engine = create_engine(dsn)
    gdf.to_postgis(table, engine, if_exists="append", index=False)

    with engine.begin() as conn:
        conn.execute(text(
            f"CREATE INDEX IF NOT EXISTS idx_{table}_geom_gist "
            f"ON {table} USING gist (geometry);"
        ))
        conn.execute(text(f"VACUUM ANALYZE {table};"))

Edge cases and failure modes

Antimeridian wrap-around

A bounding box that crosses the 180°/−180° meridian (for example a query spanning the Bering Strait or Fiji) breaks naive ST_MakeEnvelope boxes, which assume min_lon < max_lon. The box silently inverts to cover the entire globe minus your intended region. Split the query into two boxes on either side of the antimeridian and union the results, or store and query in a projected SRID that does not straddle the seam for that region.

Poles and geography distortion

Near the poles, one degree of longitude collapses to almost zero metres. A planar geometry distance in degrees becomes wildly misleading, and even geography spheroid math loses precision in a tiny neighbourhood around the exact pole. For polar datasets prefer a polar stereographic projection and index the projected geometry.

Mixed SRIDs in one column

If some rows are 4326 and others are 3857 (Web Mercator) in the same column, operators either error out or, worse, compare incompatible coordinates. Enforce the SRID in the column type (geometry(Point, 4326)) so Postgres rejects mismatched inserts, and ST_Transform at load time rather than query time so the index stays usable.

Duplicate stacked points

Geocoders frequently return the exact same centroid for many distinct inputs (a building, a ZIP centroid, a “no match” fallback). Thousands of identical points inflate the index and skew KNN results. Deduplicate to a canonical location first — the fuzzy canonical key workflow collapses these before they reach the spatial table.

Index bloat after bulk deletes

Deleting a large fraction of rows leaves dead tuples that GiST does not immediately reclaim. Subsequent scans read empty pages. Run VACUUM to reclaim space and REINDEX INDEX CONCURRENTLY if the index has grown well past its healthy size.

Performance and scaling

Scenario Technique Effect
Point table > 10M rows GiST + ST_DWithin with geography cast Sub-10 ms radius filters with correct index usage
Append-only archive sorted by geohash BRIN index Index shrinks from gigabytes to kilobytes; builds in seconds
Frequent “k nearest” lookups <-> KNN operator + LIMIT Index-ordered scan avoids sorting the whole table
Existence-only bounding-box checks Index-only scan on GiST Skips heap fetch when visibility map is current
Bulk load of tens of millions of points geopandas.to_postgis then build index Index-after-load is far faster than maintaining it per insert

Always build the index after a large bulk load, not before — maintaining a GiST index during millions of inserts is much slower than one index build over the finished table. After loading, run VACUUM ANALYZE so the planner has fresh statistics and row estimates.

Troubleshooting

Sequential scan despite an existing index

Run EXPLAIN (ANALYZE, BUFFERS). If the table is small (a few thousand rows) a sequential scan is genuinely cheaper and Postgres is correct. On a large table, check for an SRID mismatch or a function wrapping the indexed column — ST_DWithin(ST_Transform(geom, 3857), ...) cannot use an index built on the untransformed geom.

KNN returns wrong-order results

The <-> operator orders by bounding-box distance for some geometry types. For points it is exact, but if you accidentally index a different column or mix geography and geometry, ordering degrades. Confirm the ORDER BY expression references the same column and literal SRID as the index.

ST_DWithin radius looks wrong

If your radius is in metres but you forgot the geography cast, ST_DWithin(geom, point, 500) on a 4326 geometry column interprets 500 as degrees — an enormous distance. Cast both operands to geography, or reproject to a metre-based SRID.

REINDEX locks the table

Plain REINDEX INDEX foo takes an exclusive lock. On a live geocoding table use REINDEX INDEX CONCURRENTLY foo so reads and writes continue during the rebuild, at the cost of a slower rebuild.

FAQ

Should I store geocoded points as geometry or geography?

Use geometry(Point, 4326) for the common case: it is faster, supports every index type, and ST_DWithin with a geography cast handles radius queries. Switch a column to geography only when you routinely compute true great-circle distances over wide areas and cannot tolerate planar distortion. Geography indexes exist but geography operators are more expensive per comparison.

When does BRIN beat GiST for spatial data?

BRIN wins when rows are physically stored in an order that correlates with location — for example a table bulk-loaded by geohash, postal code, or a spatial sort. In that case a BRIN index is thousands of times smaller than GiST and builds in seconds, while still pruning most blocks. If inserts arrive in random order the correlation breaks and BRIN scans far more blocks than GiST, so GiST is the safe default for point queries.

Why is my ST_DWithin query not using the index?

The most common causes are an SRID mismatch between the column and the query literal, wrapping the indexed column in a function (such as ST_Transform) so the planner cannot use the index, or a table small enough that a sequential scan is genuinely cheaper. Run EXPLAIN ANALYZE, confirm both sides share one SRID, and keep the raw indexed column on the left of the operator.

How do I get an index-only scan on a spatial index?

GiST supports index-only scans when the query only needs the indexed geometry and the table is well vacuumed so the visibility map is current. Add the geometry to the index, keep autovacuum aggressive, and select only the covered column. Note that KNN ordering with the distance operator still reads heap rows for the returned tuples, so index-only scans mostly help bounding-box existence checks.

How often should I REINDEX a spatial index?

Rarely on a schedule. Let autovacuum handle bloat and only REINDEX when EXPLAIN shows the index has grown far past its healthy size or after a large bulk delete. Use REINDEX INDEX CONCURRENTLY on live tables to avoid a long exclusive lock, and prefer a nightly VACUUM ANALYZE to keep planner statistics and the visibility map fresh.