--- /dev/null
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+# This file is part of Nominatim. (https://nominatim.org)
+#
+# Copyright (C) 2023 by the Nominatim developer community.
+# For a full list of authors see the git log.
+"""
+Implementation of reverse geocoding.
+"""
+from typing import Optional
+
+import sqlalchemy as sa
+from geoalchemy2 import WKTElement
+from geoalchemy2.types import Geometry
+
+from nominatim.typing import SaColumn, SaSelect, SaTable, SaLabel, SaClause
+from nominatim.api.connection import SearchConnection
+import nominatim.api.results as nres
+from nominatim.api.logging import log
+from nominatim.api.types import AnyPoint, DataLayer, LookupDetails, GeometryFormat
+
+def _select_from_placex(t: SaTable, wkt: Optional[str] = None) -> SaSelect:
+ """ Create a select statement with the columns relevant for reverse
+ results.
+ """
+ if wkt is None:
+ distance = t.c.distance
+ else:
+ distance = t.c.geometry.ST_Distance(wkt)
+
+ return sa.select(t.c.place_id, t.c.osm_type, t.c.osm_id, t.c.name,
+ t.c.class_, t.c.type,
+ t.c.address, t.c.extratags,
+ t.c.housenumber, t.c.postcode, t.c.country_code,
+ t.c.importance, t.c.wikipedia,
+ t.c.parent_place_id, t.c.rank_address, t.c.rank_search,
+ t.c.centroid,
+ distance.label('distance'),
+ t.c.geometry.ST_Expand(0).label('bbox'))
+
+
+def _interpolated_housenumber(table: SaTable) -> SaLabel:
+ # Entries with startnumber = endnumber are legacy from version < 4.1
+ return sa.cast(table.c.startnumber
+ + sa.func.round(((table.c.endnumber - table.c.startnumber) * table.c.position)
+ / table.c.step) * table.c.step,
+ sa.Integer).label('housenumber')
+
+
+def _is_address_point(table: SaTable) -> SaClause:
+ return sa.and_(table.c.rank_address == 30,
+ sa.or_(table.c.housenumber != None,
+ table.c.name.has_key('housename')))
+
+
+class ReverseGeocoder:
+ """ Class implementing the logic for looking up a place from a
+ coordinate.
+ """
+
+ def __init__(self, conn: SearchConnection, max_rank: int, layer: DataLayer,
+ details: LookupDetails) -> None:
+ self.conn = conn
+ self.max_rank = max_rank
+ self.layer = layer
+ self.details = details
+
+
+ def _add_geometry_columns(self, sql: SaSelect, col: SaColumn) -> SaSelect:
+ if not self.details.geometry_output:
+ return sql
+
+ out = []
+
+ if self.details.geometry_simplification > 0.0:
+ col = col.ST_SimplifyPreserveTopology(self.details.geometry_simplification)
+
+ if self.details.geometry_output & GeometryFormat.GEOJSON:
+ out.append(col.ST_AsGeoJSON().label('geometry_geojson'))
+ if self.details.geometry_output & GeometryFormat.TEXT:
+ out.append(col.ST_AsText().label('geometry_text'))
+ if self.details.geometry_output & GeometryFormat.KML:
+ out.append(col.ST_AsKML().label('geometry_kml'))
+ if self.details.geometry_output & GeometryFormat.SVG:
+ out.append(col.ST_AsSVG().label('geometry_svg'))
+
+ return sql.add_columns(*out)
+
+
+ def _filter_by_layer(self, table: SaTable) -> SaColumn:
+ if self.layer & DataLayer.MANMADE:
+ exclude = []
+ if not (self.layer & DataLayer.RAILWAY):
+ exclude.append('railway')
+ if not (self.layer & DataLayer.NATURAL):
+ exclude.extend(('natural', 'water', 'waterway'))
+ return table.c.class_.not_in(tuple(exclude))
+
+ include = []
+ if self.layer & DataLayer.RAILWAY:
+ include.append('railway')
+ if not (self.layer & DataLayer.NATURAL):
+ include.extend(('natural', 'water', 'waterway'))
+ return table.c.class_.in_(tuple(include))
+
+
+ async def _find_closest_street_or_poi(self, wkt: WKTElement) -> SaRow:
+ """ Look up the clostest rank 26+ place in the database.
+ """
+ t = self.conn.t.placex
+
+ sql = _select_from_placex(t, wkt)\
+ .where(t.c.geometry.ST_DWithin(wkt, distance))\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(sa.or_(t.c.geometry.ST_GeometryType().not_in(('ST_Polygon', 'ST_MultiPolygon')),
+ t.c.centroid.ST_Distance(wkt) < distance))\
+ .order_by('distance')\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, t.c.geometry)
+
+ restrict = []
+
+ if self.layer & DataLayer.ADDRESS:
+ restrict.append(sa.and_(t.c.rank_address >= 26,
+ t.c.rank_address <= self.max_rank))
+ if self.max_rank == 30:
+ restrict.append(_is_address_point(t))
+ if self.layer & DataLayer.POI and max_rank == 30:
+ restrict.append(sa.and_(t.c.rank_search == 30,
+ t.c.class_.not_in(('place', 'building')),
+ t.c.geometry.ST_GeometryType() != 'ST_LineString'))
+ if self.layer & (DataLayer.RAILWAY | DataLayer.MANMADE | DataLayer.NATURAL):
+ restrict.append(sa.and_(t.c.rank_search >= 26,
+ tc.rank_search <= self.max_rank,
+ self._filter_by_layer(t)))
+
+ if restrict:
+ sql = sql.where(sa.or_(*restrict))
+
+ return (await self.conn.execute(sql)).one_or_none()
+
+
+ async def _find_housenumber_for_street(self, parent_place_id: int,
+ wkt: WKTElement) -> Optional[SaRow]:
+ t = conn.t.placex
+
+ sql = _select_from_placex(t, wkt)\
+ .where(t.c.geometry.ST_DWithin(wkt, 0.001))\
+ .where(t.c.parent_place_id == parent_place_id)\
+ .where(_is_address_point(t))\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .order_by('distance')\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, t.c.geometry)
+
+ return (await self.conn.execute(sql)).one_or_none()
+
+
+ async def _find_interpolation_for_street(self, parent_place_id: Optional[int],
+ wkt: WKTElement) -> Optional[SaRow]:
+ t = self.conn.t.osmline
+
+ inner = sa.select(t,
+ t.c.linegeo.ST_Distance(wkt).label('distance'),
+ t.c.linegeo.ST_LineLocatePoint(wkt).label('position'))\
+ .where(t.c.linegeo.ST_DWithin(wkt, distance))\
+ .order_by('distance')\
+ .limit(1)
+
+ if parent_place_id is not None:
+ inner = inner.where(t.c.parent_place_id == parent_place_id)
+
+ inner = inner.subquery()
+
+ sql = sa.select(inner.c.place_id, inner.c.osm_id,
+ inner.c.parent_place_id, inner.c.address,
+ _interpolated_housenumber(inner),
+ inner.c.postcode, inner.c.country_code,
+ inner.c.linegeo.ST_LineInterpolatePoint(inner.c.position).label('centroid'),
+ inner.c.distance)
+
+ if self.details.geometry_output:
+ sub = sql.subquery()
+ sql = self._add_geometry_columns(sql, sub.c.centroid)
+
+ return (await self.conn.execute(sql)).one_or_none()
+
+
+ async def _find_tiger_number_for_street(self, parent_place_id: int,
+ wkt: WKTElement) -> Optional[SaRow]:
+ t = self.conn.t.tiger
+
+ inner = sa.select(t,
+ t.c.linegeo.ST_Distance(wkt).label('distance'),
+ sa.func.ST_LineLocatePoint(t.c.linegeo, wkt).label('position'))\
+ .where(t.c.linegeo.ST_DWithin(wkt, 0.001))\
+ .where(t.c.parent_place_id == parent_place_id)\
+ .order_by('distance')\
+ .limit(1)\
+ .subquery()
+
+ sql = sa.select(inner.c.place_id,
+ inner.c.parent_place_id,
+ _interpolated_housenumber(inner),
+ inner.c.postcode,
+ inner.c.linegeo.ST_LineInterpolatePoint(inner.c.position).label('centroid'),
+ inner.c.distance)
+
+ if self.details.geometry_output:
+ sub = sql.subquery()
+ sql = self._add_geometry_columns(sql, sub.c.centroid)
+
+ return (await conn.execute(sql)).one_or_none()
+
+
+ async def lookup_street_poi(self, wkt: WKTElement) -> Optional[nres.ReverseResult]:
+ """ Find a street or POI/address for the given WKT point.
+ """
+ log().section('Reverse lookup on street/address level')
+ result = None
+ distance = 0.006
+ parent_place_id = None
+
+ row = await self._find_closest_street_or_poi(wkt)
+ log().var_dump('Result (street/building)', row)
+
+ # If the closest result was a street, but an address was requested,
+ # check for a housenumber nearby which is part of the street.
+ if row is not None:
+ if self.max_rank > 27 \
+ and self.layer & DataLayer.ADDRESS \
+ and row.rank_address <= 27:
+ distance = 0.001
+ parent_place_id = row.place_id
+ log().comment('Find housenumber for street')
+ addr_row = await self._find_housenumber_for_street(parent_place_id, wkt)
+ log().var_dump('Result (street housenumber)', addr_row)
+
+ if addr_row is not None:
+ row = addr_row
+ distance = addr_row.distance
+ elif row.country_code == 'us' and parent_place_id is not None:
+ log().comment('Find TIGER housenumber for street')
+ addr_row = await self._find_tiger_number_for_street(parent_place_id, wkt)
+ log().var_dump('Result (street Tiger housenumber)', addr_row)
+
+ if addr_row is not None:
+ result = nres.create_from_tiger_row(addr_row)
+ else:
+ distance = row.distance
+
+ # Check for an interpolation that is either closer than our result
+ # or belongs to a close street found.
+ if self.max_rank > 27 and self.layer & DataLayer.ADDRESS:
+ log().comment('Find interpolation for street')
+ addr_row = await self._find_interpolation_for_street(parent_place_id, wkt)
+ log().var_dump('Result (street interpolation)', addr_row)
+ if addr_row is not None:
+ result = nres.create_from_osmline_row(addr_row)
+
+ return result or nres.create_from_placex_row(row)
+
+
+ async def _lookup_area_address(self, wkt: WKTElement) -> Optional[SaRow]:
+ """ Lookup large addressable areas for the given WKT point.
+ """
+ log().comment('Reverse lookup by larger address area features')
+ t = self.conn.t.placex
+
+ # The inner SQL brings results in the right order, so that
+ # later only a minimum of results needs to be checked with ST_Contains.
+ inner = sa.select(t, sa.literal(0.0).label('distance'))\
+ .where(t.c.rank_search.between(5, self.max_rank))\
+ .where(t.c.rank_address.between(5, 25))\
+ .where(t.c.geometry.ST_GeometryType().in_(('ST_Polygon', 'ST_MultiPolygon')))\
+ .where(t.c.geometry.intersects(wkt))\
+ .where(t.c.name != None)\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(t.c.type != 'postcode')\
+ .order_by(sa.desc(t.c.rank_search))\
+ .limit(50)\
+ .subquery()
+
+ sql = _select_from_placex(inner)\
+ .where(inner.c.geometry.ST_Contains(wkt))\
+ .order_by(sa.desc(inner.c.rank_search))\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ address_row = (await self.conn.execute(sql)).one_or_none()
+ log().var_dump('Result (area)', address_row)
+
+ if address_row is not None and address_row.rank_search < max_rank:
+ log().comment('Search for better matching place nodes inside the area')
+ inner = sa.select(t,
+ t.c.geometry.ST_Distance(wkt).label('distance'))\
+ .where(t.c.osm_type == 'N')\
+ .where(t.c.rank_search > address_row.rank_search)\
+ .where(t.c.rank_search <= max_rank)\
+ .where(t.c.rank_address.between(5, 25))\
+ .where(t.c.name != None)\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(t.c.type != 'postcode')\
+ .where(t.c.geometry
+ .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
+ .intersects(wkt))\
+ .order_by(sa.desc(t.c.rank_search))\
+ .limit(50)\
+ .subquery()
+
+ touter = conn.t.placex.alias('outer')
+ sql = _select_from_placex(inner)\
+ .where(touter.c.place_id == address_row.place_id)\
+ .where(touter.c.geometry.ST_Contains(inner.c.geometry))\
+ .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
+ .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ place_address_row = (await self.conn.execute(sql)).one_or_none()
+ log().var_dump('Result (place node)', place_address_row)
+
+ if place_address_row is not None:
+ return place_address_row
+
+ return address_row
+
+
+ async def _lookup_area_others(self, wkt: WKTElement) -> Optional[SaRow]:
+ t = conn.t.placex
+
+ inner = sa.select(t, t.c.geometry.ST_Distance(wkt).label('distance'))\
+ .where(t.c.rank_address == 0)\
+ .where(t.c.rank_search.between(5, self.max_rank))\
+ .where(t.c.name != None)\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(self._filter_by_layer(t))\
+ .where(sa.func.reverse_buffered_extent(t.c.geometry, type_=Geometry)
+ .intersects(wkt))\
+ .order_by(sa.desc(t.c.rank_search))\
+ .limit(50)
+
+ sql = _select_from_placex(inner)\
+ .where(sa._or(inner.c.geometry.ST_GeometryType().not_in(('ST_Polygon', 'ST_MultiPolygon')),
+ inner.c.geometry.ST_Contains(wkt)))\
+ .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ row = (await self.conn.execute(sql)).one_or_none()
+ log().var_dump('Result (non-address feature)', row)
+
+ return row
+
+
+ async def lookup_area(self, wkt: WKTElement) -> Optional[nres.ReverseResult]:
+ """ Lookup large areas for the given WKT point.
+ """
+ log().section('Reverse lookup by larger area features')
+ t = self.conn.t.placex
+
+ if self.layer & DataLayer.ADDRESS:
+ address_row = await self._lookup_area_address(wkt)
+ address_distance = address_row.distance
+ else:
+ address_row = None
+ address_distance = 1000
+
+ if self.layer & (~DataLayer.ADDRESS & ~DataLayer.POI):
+ other_row = await self._lookup_area_others(wkt)
+ other_distance = other_row.distance
+ else:
+ other_row = None
+ other_distance = 1000
+
+ result = address_row if address_distance <= other_distance else other_row
+
+ return nres.create_from_placex_row(result)
+
+
+ async def lookup_country(self, wkt: WKTElement) -> Optional[nres.ReverseResult]:
+ """ Lookup the country for the given WKT point.
+ """
+ log().section('Reverse lookup by country code')
+ t = self.conn.t.country_grid
+ sql = sa.select(t.c.country_code).distinct()\
+ .where(t.c.geometry.ST_Contains(wkt))
+
+ ccodes = tuple((r[0] for r in await self.conn.execute(sql)))
+ log().var_dump('Country codes', ccodes)
+
+ if not ccodes:
+ return None
+
+ if self.layer & DataLayer.ADDRESS and self.max_rank > 4:
+ log().comment('Search for place nodes in country')
+
+ t = conn.t.placex
+ inner = sa.select(t,
+ t.c.geometry.ST_Distance(wkt).label('distance'))\
+ .where(t.c.osm_type == 'N')\
+ .where(t.c.rank_search > 4)\
+ .where(t.c.rank_search <= self.max_rank)\
+ .where(t.c.rank_address.between(5, 25))\
+ .where(t.c.name != None)\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(t.c.type != 'postcode')\
+ .where(t.c.country_code.in_(ccodes))\
+ .where(t.c.geometry
+ .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
+ .intersects(wkt))\
+ .order_by(sa.desc(t.c.rank_search))\
+ .limit(50)\
+ .subquery()
+
+ sql = _select_from_placex(inner)\
+ .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
+ .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ address_row = (await self.conn.execute(sql)).one_or_none()
+ log().var_dump('Result (addressable place node)', address_row)
+ else:
+ address_row = None
+
+ if layer & (~DataLayer.ADDRESS & ~DataLayer.POI) and self.max_rank > 4:
+ log().comment('Search for non-address features inside country')
+
+ t = conn.t.placex
+ inner = sa.select(t, t.c.geometry.ST_Distance(wkt).label('distance'))\
+ .where(t.c.rank_address == 0)\
+ .where(t.c.rank_search.between(5, self.max_rank))\
+ .where(t.c.name != None)\
+ .where(t.c.indexed_status == 0)\
+ .where(t.c.linked_place_id == None)\
+ .where(self._filter_by_layer(t))\
+ .where(t.c.country_code.in_(ccode))\
+ .where(sa.func.reverse_buffered_extent(t.c.geometry, type_=Geometry)
+ .intersects(wkt))\
+ .order_by(sa.desc(t.c.rank_search))\
+ .limit(50)\
+ .subquery()
+
+ sql = _select_from_placex(inner)\
+ .where(sa._or(inner.c.geometry.ST_GeometryType().not_in(('ST_Polygon', 'ST_MultiPolygon')),
+ inner.c.geometry.ST_Contains(wkt)))\
+ .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
+ .limit(1)
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ other_row = (await self.conn.execute(sql)).one_or_none()
+ log().var_dump('Result (non-address feature)', other_row)
+ else:
+ other_row = None
+
+ if layer & DataLayer.ADDRESS and address_row is None and other_row is None:
+ # Still nothing, then return a country with the appropriate country code.
+ t = conn.t.placex
+ sql = _select_from_placex(t, wkt)\
+ .where(t.c.country_code.in_(ccodes))\
+ .where(t.c.rank_address == 4)\
+ .where(t.c.rank_search == 4)\
+ .where(t.c.linked_place_id == None)\
+ .order_by('distance')
+
+ sql = self._add_geometry_columns(sql, inner.c.geometry)
+
+ address_row = (await self.conn.execute(sql)).one_or_none()
+
+ return nres.create_from_placex_row(_get_closest_row(address_row, other_row))
+
+
+ async def lookup(self, coord: AnyPoint) -> Optional[nres.ReverseResult]:
+ """ Look up a single coordinate. Returns the place information,
+ if a place was found near the coordinates or None otherwise.
+ """
+ log().function('reverse_lookup',
+ coord=coord, max_rank=self.max_rank,
+ layer=self.layer, details=self.details)
+
+
+ wkt = WKTElement(f'POINT({coord[0]} {coord[1]})', srid=4326)
+
+ result: Optional[ReverseResult] = None
+
+ if max_rank >= 26:
+ result = await self.lookup_street_poi(wkt)
+ if result is None and max_rank > 4:
+ result = await self.lookup_area(wkt)
+ if result is None:
+ result = await self.lookup_country(wkt)
+ if result is not None:
+ await nres.add_result_details(self.conn, result, self.details)
+
+ return result
projects / nominatim.git / commitdiff