Reverse-Geocoding POI Nodes to Administrative Boundaries

Analysts ask the graph for “all charging stations in Bavaria” or “delivery hubs inside the Paris city limits,” and the query has nothing to filter on — the POI nodes carry a location point but no notion of the administrative area that contains them. The naive fix, a point.distance() scan against a boundary centroid, is both wrong (a point can be near a centroid yet outside the polygon) and slow (it reads every boundary for every POI). The correct operation is point-in-polygon containment: which admin polygons actually enclose this coordinate, across every level from country down to postal code. This page resolves each POI against nested boundaries with a shapely STRtree prefilter, then writes the result as a :WITHIN hierarchy in the graph so area filters become a cheap relationship traversal instead of a geometry test at query time. It builds on the POI enrichment workflows that established these nodes.

Prerequisites & Versions

The STRtree.query(point, predicate="contains") call — an R-tree bounding-box prefilter fused with the exact containment test — requires shapely 2.x. The write side is the official async driver.

pip install "shapely>=2.0" "neo4j>=5.18"

This guide assumes the POI nodes already exist with a stable id and a native location point, produced upstream by the OSM data ingestion pipelines stage and following the node and edge spatial mapping conventions — coordinates as a real point, not detached lat/lon strings. Create the anchoring constraints so every MERGE seeks an index:

CREATE CONSTRAINT poi_id_unique IF NOT EXISTS
FOR (p:POI) REQUIRE p.id IS UNIQUE;

CREATE CONSTRAINT admin_area_id IF NOT EXISTS
FOR (a:AdminArea) REQUIRE a.id IS UNIQUE;

Implementation

The boundaries load once into an STRtree; each POI is a point query that returns the containment chain coarsest-first. The graph write upserts every area on that chain, stitches the hierarchy from each area’s parent_id, and attaches the POI to its most specific enclosing area.

import asyncio
import json
from dataclasses import dataclass

from neo4j import AsyncGraphDatabase
from shapely import STRtree, Point
from shapely.geometry import shape


@dataclass(frozen=True)
class AdminArea:
    admin_id: str
    name: str
    level: int           # OSM admin_level: 2 country, 4 region, 8 city, 10 postal
    parent_id: str | None


def load_boundaries(geojson_path: str) -> list[tuple[AdminArea, object]]:
    """Load admin polygons (WGS84 GeoJSON) as (AdminArea, shapely geometry) pairs."""
    with open(geojson_path) as fh:
        collection = json.load(fh)
    pairs = []
    for feature in collection["features"]:
        props = feature["properties"]
        area = AdminArea(
            admin_id=props["admin_id"],
            name=props["name"],
            level=int(props["admin_level"]),
            parent_id=props.get("parent_id"),
        )
        pairs.append((area, shape(feature["geometry"])))
    return pairs


class ReverseGeocoder:
    """Resolves a WGS84 coordinate to its nested admin containment chain."""

    def __init__(self, boundaries: list[tuple[AdminArea, object]]):
        self.areas = [area for area, _ in boundaries]
        self.geoms = [geom for _, geom in boundaries]
        self.tree = STRtree(self.geoms)   # R-tree over polygon bounding boxes

    def locate(self, lat: float, lon: float) -> list[AdminArea]:
        # shapely is planar and axis order is (x=lon, y=lat).
        point = Point(lon, lat)
        # query() prefilters by bounding box; predicate="contains" then applies the
        # exact point-in-polygon test, so only truly enclosing polygons return.
        idx = self.tree.query(point, predicate="contains")
        matched = [self.areas[i] for i in idx]
        matched.sort(key=lambda a: a.level)   # coarsest (country) first
        return matched


WRITE_HIERARCHY = """
UNWIND $rows AS row
// 1. upsert every admin area on the containment chain
UNWIND row.areas AS area
MERGE (a:AdminArea {id: area.admin_id})
  ON CREATE SET a.name = area.name, a.level = area.level
// 2. materialise the hierarchy: each area WITHIN its parent
WITH row, area
WHERE area.parent_id IS NOT NULL
MERGE (child:AdminArea {id: area.admin_id})
MERGE (parent:AdminArea {id: area.parent_id})
MERGE (child)-[:WITHIN]->(parent)
// 3. attach the POI to its most specific enclosing area
WITH DISTINCT row
MATCH (p:POI {id: row.poi_id})
MATCH (finest:AdminArea {id: row.finest_id})
MERGE (p)-[:WITHIN]->(finest)
"""


async def _write(tx, rows):
    await tx.run(WRITE_HIERARCHY, rows=rows)


def _area_dict(area: AdminArea) -> dict:
    return {"admin_id": area.admin_id, "name": area.name,
            "level": area.level, "parent_id": area.parent_id}


async def stamp_pois(driver, geocoder: ReverseGeocoder, pois: list[tuple[str, float, float]]):
    rows = []
    for poi_id, lat, lon in pois:
        chain = geocoder.locate(lat, lon)
        if not chain:
            # A point that lands in no boundary (offshore, border sliver, or a
            # coverage gap) is skipped rather than mis-assigned. See failures below.
            continue
        rows.append({
            "poi_id": poi_id,
            "finest_id": chain[-1].admin_id,
            "areas": [_area_dict(a) for a in chain],
        })
    if rows:
        async with driver.session() as session:
            await session.execute_write(_write, rows)
    return len(rows)


async def main():
    geocoder = ReverseGeocoder(load_boundaries("admin_boundaries.geojson"))
    pois = [
        ("poi:cafe/1", 48.8566, 2.3522),    # central Paris
        ("poi:depot/2", 43.2965, 5.3698),   # Marseille port
    ]
    driver = AsyncGraphDatabase.driver(
        "neo4j+s://your-cluster.databases.neo4j.io",
        auth=("neo4j", "secure-password"),
        max_connection_pool_size=20,
    )
    try:
        stamped = await stamp_pois(driver, geocoder, pois)
        print(f"Stamped {stamped} POIs with their admin hierarchy.")
    finally:
        await driver.close()


if __name__ == "__main__":
    asyncio.run(main())

How It Works

Three decisions carry the correctness and the speed:

• STRtree prefilter fused with the exact test. STRtree.query(point, predicate="contains") first descends the R-tree to the handful of polygons whose bounding box covers the point, then runs the real contains predicate on only those candidates. A metropolitan boundary set can hold thousands of postal polygons; the tree turns a linear scan into a logarithmic descent plus a few exact tests, which is what keeps per-POI cost near constant.
• Hierarchy from parent_id, not from geometry. Rather than infer nesting by testing polygon-inside-polygon (fragile on shared borders), each area declares its parent, and the write stitches (:AdminArea)-[:WITHIN]->(:AdminArea) directly. The POI links only to its finest enclosing area; the chain up to country is reached by traversing :WITHIN. That normalization means “everything in this region” is a single variable-length traversal — MATCH (p:POI)-[:WITHIN*]->(r:AdminArea {id: $region}) — with no geometry evaluated at query time.
• Coarsest-first ordering. Sorting the match by level makes chain[-1] the most specific area deterministically, so finest_id is unambiguous even when a coordinate legitimately sits inside four levels at once.

Storing containment as edges is what lets downstream analytics correlate POIs, demographics, and routes by area without a spatial join on every read — the geometric join is paid once here. When you do need to intersect POIs against arbitrary external polygons at query time rather than a fixed hierarchy, that is the province of spatial join techniques, and demographic attributes attach to the same nodes through enriching POI data with real-time demographics.

Common Failure Patterns

1. Points on borders and boundary precision. A POI exactly on a shared edge, or a simplified boundary whose vertices drift a few meters, can match zero polygons or two adjacent ones. contains treats the boundary as open, so a point on the line is not contained. When a coordinate must resolve, fall back from contains to intersects and pick the smallest matching area, or buffer the query point by a sub-meter tolerance:

idx = geocoder.tree.query(point, predicate="intersects")   # includes the boundary
if len(idx) == 0:
    idx = geocoder.tree.query(point.buffer(1e-6), predicate="intersects")

2. Skipping the R-tree — O(N) per point. Iterating for geom in geoms: if geom.contains(point) looks equivalent and is catastrophically slower: every POI tests every polygon’s full vertex ring. On a national boundary set this turns a minutes job into hours. Always query the STRtree so the bounding-box prefilter runs first; the exact test only touches candidates.

3. CRS and axis-order mismatch. shapely is a planar library that trusts whatever coordinates you hand it. Two traps compound here: constructing Point(lat, lon) instead of Point(lon, lat) silently mirrors every POI across the diagonal, and testing WGS84 points against boundaries stored in a projected CRS (say EPSG:3857) produces containment that is wrong everywhere. Normalize both point and polygons to the same CRS at load time, and keep the (lon, lat) axis order explicit:

point = Point(lon, lat)   # x = longitude, y = latitude — never (lat, lon)

Performance Notes

The dominant cost is the containment pass, not the graph write. For $Q$ POIs tested against $N$ boundary polygons where the tree returns $k$ bounding-box candidates per point and each candidate polygon has $\bar{v}$ vertices, total work is

$$C \approx Q\bigl(\log N + k \cdot \bar{v}\bigr)$$

The $\log N$ term is the tree descent; $k \cdot \bar{v}$ is the exact test over the few survivors. Without the tree the cost degrades to $Q \cdot N \cdot \bar{v}$ — the difference between log N and N candidates is the entire justification for the STRtree. Build the tree once and reuse it across the whole POI batch; rebuilding per point discards the amortization. On the write side, batch the resolved rows into a single UNWIND per transaction (the code sends all POIs in one execute_write) so the admin-area MERGE set amortizes across the batch, and keep boundary geometries pre-simplified to a tolerance matched to your smallest meaningful area — over-detailed polygons inflate $\bar{v}$ without changing which POIs they contain. These area edges then serve the same downstream role as the rest of the POI enrichment workflows.

Related

• Enriching POI data with real-time demographics — attaching mobility and census attributes to the same POI nodes.
• Spatial join techniques — intersecting nodes against external polygons at query time rather than a fixed hierarchy.
• Building automated OSM-to-graph ETL pipelines — where the POI id and location anchors are first established.
• Node and edge spatial mapping — the point-property conventions the containment test depends on.

This guide is part of POI Enrichment Workflows, within the broader Spatial Graph Construction & OSM Ingestion guide.