From 90b009fc906f4938c9d457d90205e81ec5d167ad Mon Sep 17 00:00:00 2001
From: Edward Betts <edward@4angle.com>
Date: Thu, 12 Feb 2026 18:19:17 +0000
Subject: [PATCH] Return nearest parish/polygon for points in the sea near the
 coast

When a point falls in the sea within a country boundary (e.g. Scotland
Q22 or England Q21) but not within any specific parish, find the nearest
civil parish or high-admin-level polygon within ~0.15 degrees (~10-17km).

- Add get_nearest_scotland_code() for Scottish civil parishes
- Add Polygon.nearest() for general OSM polygon lookup
- Use .limit(1).scalar() in scotland.py for cleaner returns
---
 geocode/model.py    | 26 ++++++++++++++++++++++++++
 geocode/scotland.py | 21 ++++++++++++++++++---
 lookup.py           | 30 ++++++++++++++++++++++++++++++
 3 files changed, 74 insertions(+), 3 deletions(-)

diff --git a/geocode/model.py b/geocode/model.py
index e077688..50cc1e8 100644
--- a/geocode/model.py
+++ b/geocode/model.py
@@ -65,6 +65,32 @@ class Polygon(Base):
         return q  # type: ignore
 
 
+    @classmethod
+    def nearest(
+        cls,
+        lat: float,
+        lon: float,
+        min_admin_level: int = 7,
+        max_degrees: float = 0.15,
+    ) -> "Polygon | None":
+        """Find nearest polygon with high admin_level within distance.
+
+        Uses degree-based filter (~0.15° ≈ 10-17km at UK latitudes) to
+        leverage the spatial index on the geometry column.
+        """
+        point = func.ST_SetSRID(func.ST_MakePoint(lon, lat), 4326)
+        return (
+            cls.query.filter(
+                cls.admin_level.isnot(None),
+                cls.admin_level.regexp_match(r"^\d+$"),
+                cast(cls.admin_level, Integer) >= min_admin_level,
+                func.ST_DWithin(cls.way, point, max_degrees),
+            )
+            .order_by(func.ST_Distance(cls.way, point))
+            .first()
+        )
+
+
 class Scotland(Base):
     """Civil parishes in Scotland."""
 
diff --git a/geocode/scotland.py b/geocode/scotland.py
index f377c84..73753cf 100644
--- a/geocode/scotland.py
+++ b/geocode/scotland.py
@@ -9,9 +9,24 @@ from geocode.model import Scotland
 def get_scotland_code(lat: float, lon: float) -> str | None:
     """Find civil parish in Scotland for given lat/lon."""
     point = func.ST_Transform(func.ST_SetSRID(func.ST_MakePoint(lon, lat), 4326), 27700)
-    result = (
+    return (
         session.query(Scotland.code)
         .filter(func.ST_Contains(Scotland.geom, point))
-        .first()
+        .limit(1)
+        .scalar()
+    )
+
+
+def get_nearest_scotland_code(
+    lat: float, lon: float, max_degrees: float = 0.15
+) -> str | None:
+    """Find nearest civil parish in Scotland within distance."""
+    point = func.ST_SetSRID(func.ST_MakePoint(lon, lat), 4326)
+    geom_4326 = func.ST_Transform(Scotland.geom, 4326)
+    return (
+        session.query(Scotland.code)
+        .filter(func.ST_DWithin(geom_4326, point, max_degrees))
+        .order_by(func.ST_Distance(geom_4326, point))
+        .limit(1)
+        .scalar()
     )
-    return result[0] if result else None
diff --git a/lookup.py b/lookup.py
index fa2ff45..70d618b 100755
--- a/lookup.py
+++ b/lookup.py
@@ -136,6 +136,36 @@ def lat_lon_to_wikidata(lat: float, lon: float) -> dict[str, typing.Any]:
     if admin_level >= 7:
         return {"elements": elements, "result": result}
 
+    # Point is in Scotland but not in a specific parish — try nearest parish
+    if result.get("wikidata") == "Q22":
+        nearby_code = scotland.get_nearest_scotland_code(lat, lon)
+        if nearby_code:
+            rows = wikidata.lookup_scottish_parish_in_wikidata(nearby_code)
+            add_missing_commons_cat(rows)
+            hit = wikidata.commons_from_rows(rows)
+            if hit:
+                nearby_result = wikidata.build_dict(hit, lat, lon)
+                if not nearby_result.get("missing"):
+                    return {"elements": elements, "result": nearby_result}
+
+    # Point is in a broad area (e.g. country) — try nearest specific polygon
+    nearby = model.Polygon.nearest(lat, lon)
+    if nearby and nearby.tags:
+        tags: typing.Mapping[str, str] = nearby.tags
+        al = get_admin_level(tags)
+        hit = (
+            hit_from_wikidata_tag(tags)
+            or hit_from_ref_gss_tag(tags)
+            or hit_from_name(tags, lat, lon)
+        )
+        if hit:
+            hit["admin_level"] = al
+            hit["element"] = nearby.osm_id
+            hit["geojson"] = typing.cast(str, nearby.geojson_str)
+            nearby_result = wikidata.build_dict(hit, lat, lon)
+            if not nearby_result.get("missing"):
+                return {"elements": elements, "result": nearby_result}
+
     row = wikidata.geosearch(lat, lon)
     if row:
         hit = wikidata.commons_from_rows([row])