1 # SPDX-License-Identifier: GPL-3.0-or-later
3 # This file is part of Nominatim. (https://nominatim.org)
5 # Copyright (C) 2023 by the Nominatim developer community.
6 # For a full list of authors see the git log.
8 Implementation of reverse geocoding.
10 from typing import Optional, List, Callable, Type, Tuple, Dict, Any, cast, Union
13 import sqlalchemy as sa
15 from nominatim.typing import SaColumn, SaSelect, SaFromClause, SaLabel, SaRow,\
16 SaBind, SaLambdaSelect
17 from nominatim.api.connection import SearchConnection
18 import nominatim.api.results as nres
19 from nominatim.api.logging import log
20 from nominatim.api.types import AnyPoint, DataLayer, ReverseDetails, GeometryFormat, Bbox
21 from nominatim.db.sqlalchemy_types import Geometry
22 import nominatim.db.sqlalchemy_functions as snfn
24 # In SQLAlchemy expression which compare with NULL need to be expressed with
26 # pylint: disable=singleton-comparison
28 RowFunc = Callable[[Optional[SaRow], Type[nres.ReverseResult]], Optional[nres.ReverseResult]]
30 WKT_PARAM: SaBind = sa.bindparam('wkt', type_=Geometry)
31 MAX_RANK_PARAM: SaBind = sa.bindparam('max_rank')
33 def no_index(expr: SaColumn) -> SaColumn:
34 """ Wrap the given expression, so that the query planner will
35 refrain from using the expression for index lookup.
37 return sa.func.coalesce(sa.null(), expr) # pylint: disable=not-callable
40 def _select_from_placex(t: SaFromClause, use_wkt: bool = True) -> SaSelect:
41 """ Create a select statement with the columns relevant for reverse
45 distance = t.c.distance
46 centroid = t.c.centroid
48 distance = t.c.geometry.ST_Distance(WKT_PARAM)
49 centroid = sa.case((t.c.geometry.is_line_like(), t.c.geometry.ST_ClosestPoint(WKT_PARAM)),
50 else_=t.c.centroid).label('centroid')
53 return sa.select(t.c.place_id, t.c.osm_type, t.c.osm_id, t.c.name,
55 t.c.address, t.c.extratags,
56 t.c.housenumber, t.c.postcode, t.c.country_code,
57 t.c.importance, t.c.wikipedia,
58 t.c.parent_place_id, t.c.rank_address, t.c.rank_search,
60 distance.label('distance'),
61 t.c.geometry.ST_Expand(0).label('bbox'))
64 def _interpolated_housenumber(table: SaFromClause) -> SaLabel:
65 return sa.cast(table.c.startnumber
66 + sa.func.round(((table.c.endnumber - table.c.startnumber) * table.c.position)
67 / table.c.step) * table.c.step,
68 sa.Integer).label('housenumber')
71 def _interpolated_position(table: SaFromClause) -> SaLabel:
72 fac = sa.cast(table.c.step, sa.Float) / (table.c.endnumber - table.c.startnumber)
73 rounded_pos = sa.func.round(table.c.position / fac) * fac
75 (table.c.endnumber == table.c.startnumber, table.c.linegeo.ST_Centroid()),
76 else_=table.c.linegeo.ST_LineInterpolatePoint(rounded_pos)).label('centroid')
79 def _locate_interpolation(table: SaFromClause) -> SaLabel:
80 """ Given a position, locate the closest point on the line.
82 return sa.case((table.c.linegeo.is_line_like(),
83 table.c.linegeo.ST_LineLocatePoint(WKT_PARAM)),
84 else_=0).label('position')
87 def _is_address_point(table: SaFromClause) -> SaColumn:
88 return sa.and_(table.c.rank_address == 30,
89 sa.or_(table.c.housenumber != None,
90 table.c.name.has_key('addr:housename')))
93 def _get_closest(*rows: Optional[SaRow]) -> Optional[SaRow]:
94 return min(rows, key=lambda row: 1000 if row is None else row.distance)
97 class ReverseGeocoder:
98 """ Class implementing the logic for looking up a place from a
102 def __init__(self, conn: SearchConnection, params: ReverseDetails) -> None:
106 self.bind_params: Dict[str, Any] = {'max_rank': params.max_rank}
110 def max_rank(self) -> int:
111 """ Return the maximum configured rank.
113 return self.params.max_rank
116 def has_geometries(self) -> bool:
117 """ Check if any geometries are requested.
119 return bool(self.params.geometry_output)
122 def layer_enabled(self, *layer: DataLayer) -> bool:
123 """ Return true when any of the given layer types are requested.
125 return any(self.params.layers & l for l in layer)
128 def layer_disabled(self, *layer: DataLayer) -> bool:
129 """ Return true when none of the given layer types is requested.
131 return not any(self.params.layers & l for l in layer)
134 def has_feature_layers(self) -> bool:
135 """ Return true if any layer other than ADDRESS or POI is requested.
137 return self.layer_enabled(DataLayer.RAILWAY, DataLayer.MANMADE, DataLayer.NATURAL)
140 def _add_geometry_columns(self, sql: SaLambdaSelect, col: SaColumn) -> SaSelect:
143 if self.params.geometry_simplification > 0.0:
144 col = sa.func.ST_SimplifyPreserveTopology(col, self.params.geometry_simplification)
146 if self.params.geometry_output & GeometryFormat.GEOJSON:
147 out.append(sa.func.ST_AsGeoJSON(col).label('geometry_geojson'))
148 if self.params.geometry_output & GeometryFormat.TEXT:
149 out.append(sa.func.ST_AsText(col).label('geometry_text'))
150 if self.params.geometry_output & GeometryFormat.KML:
151 out.append(sa.func.ST_AsKML(col).label('geometry_kml'))
152 if self.params.geometry_output & GeometryFormat.SVG:
153 out.append(sa.func.ST_AsSVG(col).label('geometry_svg'))
155 return sql.add_columns(*out)
158 def _filter_by_layer(self, table: SaFromClause) -> SaColumn:
159 if self.layer_enabled(DataLayer.MANMADE):
161 if self.layer_disabled(DataLayer.RAILWAY):
162 exclude.append('railway')
163 if self.layer_disabled(DataLayer.NATURAL):
164 exclude.extend(('natural', 'water', 'waterway'))
165 return table.c.class_.not_in(tuple(exclude))
168 if self.layer_enabled(DataLayer.RAILWAY):
169 include.append('railway')
170 if self.layer_enabled(DataLayer.NATURAL):
171 include.extend(('natural', 'water', 'waterway'))
172 return table.c.class_.in_(tuple(include))
175 async def _find_closest_street_or_poi(self, distance: float) -> Optional[SaRow]:
176 """ Look up the closest rank 26+ place in the database, which
177 is closer than the given distance.
179 t = self.conn.t.placex
181 # PostgreSQL must not get the distance as a parameter because
182 # there is a danger it won't be able to proberly estimate index use
183 # when used with prepared statements
184 diststr = sa.text(f"{distance}")
186 sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
187 .where(t.c.geometry.ST_DWithin(WKT_PARAM, diststr))
188 .where(t.c.indexed_status == 0)
189 .where(t.c.linked_place_id == None)
190 .where(sa.or_(sa.not_(t.c.geometry.is_area()),
191 t.c.centroid.ST_Distance(WKT_PARAM) < diststr))
192 .order_by('distance')
195 if self.has_geometries():
196 sql = self._add_geometry_columns(sql, t.c.geometry)
198 restrict: List[Union[SaColumn, Callable[[], SaColumn]]] = []
200 if self.layer_enabled(DataLayer.ADDRESS):
201 max_rank = min(29, self.max_rank)
202 restrict.append(lambda: no_index(t.c.rank_address).between(26, max_rank))
203 if self.max_rank == 30:
204 restrict.append(lambda: _is_address_point(t))
205 if self.layer_enabled(DataLayer.POI) and self.max_rank == 30:
206 restrict.append(lambda: sa.and_(no_index(t.c.rank_search) == 30,
207 t.c.class_.not_in(('place', 'building')),
208 sa.not_(t.c.geometry.is_line_like())))
209 if self.has_feature_layers():
210 restrict.append(sa.and_(no_index(t.c.rank_search).between(26, MAX_RANK_PARAM),
211 no_index(t.c.rank_address) == 0,
212 self._filter_by_layer(t)))
217 sql = sql.where(sa.or_(*restrict))
219 return (await self.conn.execute(sql, self.bind_params)).one_or_none()
222 async def _find_housenumber_for_street(self, parent_place_id: int) -> Optional[SaRow]:
223 t = self.conn.t.placex
225 sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
226 .where(t.c.geometry.ST_DWithin(WKT_PARAM, 0.001))
227 .where(t.c.parent_place_id == parent_place_id)
228 .where(_is_address_point(t))
229 .where(t.c.indexed_status == 0)
230 .where(t.c.linked_place_id == None)
231 .order_by('distance')
234 if self.has_geometries():
235 sql = self._add_geometry_columns(sql, t.c.geometry)
237 return (await self.conn.execute(sql, self.bind_params)).one_or_none()
240 async def _find_interpolation_for_street(self, parent_place_id: Optional[int],
241 distance: float) -> Optional[SaRow]:
242 t = self.conn.t.osmline
244 sql: Any = sa.lambda_stmt(lambda:
246 t.c.linegeo.ST_Distance(WKT_PARAM).label('distance'),
247 _locate_interpolation(t))
248 .where(t.c.linegeo.ST_DWithin(WKT_PARAM, distance))
249 .where(t.c.startnumber != None)
250 .order_by('distance')
253 if parent_place_id is not None:
254 sql += lambda s: s.where(t.c.parent_place_id == parent_place_id)
256 def _wrap_query(base_sql: SaLambdaSelect) -> SaSelect:
257 inner = base_sql.subquery('ipol')
259 return sa.select(inner.c.place_id, inner.c.osm_id,
260 inner.c.parent_place_id, inner.c.address,
261 _interpolated_housenumber(inner),
262 _interpolated_position(inner),
263 inner.c.postcode, inner.c.country_code,
268 if self.has_geometries():
269 sub = sql.subquery('geom')
270 sql = self._add_geometry_columns(sa.select(sub), sub.c.centroid)
272 return (await self.conn.execute(sql, self.bind_params)).one_or_none()
275 async def _find_tiger_number_for_street(self, parent_place_id: int) -> Optional[SaRow]:
276 t = self.conn.t.tiger
278 def _base_query() -> SaSelect:
280 t.c.linegeo.ST_Distance(WKT_PARAM).label('distance'),
281 _locate_interpolation(t))\
282 .where(t.c.linegeo.ST_DWithin(WKT_PARAM, 0.001))\
283 .where(t.c.parent_place_id == parent_place_id)\
284 .order_by('distance')\
288 return sa.select(inner.c.place_id,
289 inner.c.parent_place_id,
290 _interpolated_housenumber(inner),
291 _interpolated_position(inner),
295 sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
297 if self.has_geometries():
298 sub = sql.subquery('geom')
299 sql = self._add_geometry_columns(sa.select(sub), sub.c.centroid)
301 return (await self.conn.execute(sql, self.bind_params)).one_or_none()
304 async def lookup_street_poi(self) -> Tuple[Optional[SaRow], RowFunc]:
305 """ Find a street or POI/address for the given WKT point.
307 log().section('Reverse lookup on street/address level')
309 parent_place_id = None
311 row = await self._find_closest_street_or_poi(distance)
312 row_func: RowFunc = nres.create_from_placex_row
313 log().var_dump('Result (street/building)', row)
315 # If the closest result was a street, but an address was requested,
316 # check for a housenumber nearby which is part of the street.
318 if self.max_rank > 27 \
319 and self.layer_enabled(DataLayer.ADDRESS) \
320 and row.rank_address <= 27:
322 parent_place_id = row.place_id
323 log().comment('Find housenumber for street')
324 addr_row = await self._find_housenumber_for_street(parent_place_id)
325 log().var_dump('Result (street housenumber)', addr_row)
327 if addr_row is not None:
329 row_func = nres.create_from_placex_row
330 distance = addr_row.distance
331 elif row.country_code == 'us' and parent_place_id is not None:
332 log().comment('Find TIGER housenumber for street')
333 addr_row = await self._find_tiger_number_for_street(parent_place_id)
334 log().var_dump('Result (street Tiger housenumber)', addr_row)
336 if addr_row is not None:
337 row_func = cast(RowFunc,
338 functools.partial(nres.create_from_tiger_row,
339 osm_type=row.osm_type,
343 distance = row.distance
345 # Check for an interpolation that is either closer than our result
346 # or belongs to a close street found.
347 if self.max_rank > 27 and self.layer_enabled(DataLayer.ADDRESS):
348 log().comment('Find interpolation for street')
349 addr_row = await self._find_interpolation_for_street(parent_place_id,
351 log().var_dump('Result (street interpolation)', addr_row)
352 if addr_row is not None:
354 row_func = nres.create_from_osmline_row
359 async def _lookup_area_address(self) -> Optional[SaRow]:
360 """ Lookup large addressable areas for the given WKT point.
362 log().comment('Reverse lookup by larger address area features')
363 t = self.conn.t.placex
365 def _base_query() -> SaSelect:
366 # The inner SQL brings results in the right order, so that
367 # later only a minimum of results needs to be checked with ST_Contains.
368 inner = sa.select(t, sa.literal(0.0).label('distance'))\
369 .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
370 .where(t.c.geometry.intersects(WKT_PARAM))\
371 .where(snfn.select_index_placex_geometry_reverse_lookuppolygon('placex'))\
372 .order_by(sa.desc(t.c.rank_search))\
376 return _select_from_placex(inner, False)\
377 .where(inner.c.geometry.ST_Contains(WKT_PARAM))\
378 .order_by(sa.desc(inner.c.rank_search))\
381 sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
382 if self.has_geometries():
383 sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
385 address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
386 log().var_dump('Result (area)', address_row)
388 if address_row is not None and address_row.rank_search < self.max_rank:
389 log().comment('Search for better matching place nodes inside the area')
391 address_rank = address_row.rank_search
392 address_id = address_row.place_id
394 def _place_inside_area_query() -> SaSelect:
397 t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
398 .where(t.c.rank_search > address_rank)\
399 .where(t.c.rank_search <= MAX_RANK_PARAM)\
400 .where(t.c.indexed_status == 0)\
401 .where(snfn.select_index_placex_geometry_reverse_lookupplacenode('placex'))\
403 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
404 .intersects(WKT_PARAM))\
405 .order_by(sa.desc(t.c.rank_search))\
409 touter = t.alias('outer')
410 return _select_from_placex(inner, False)\
411 .join(touter, touter.c.geometry.ST_Contains(inner.c.geometry))\
412 .where(touter.c.place_id == address_id)\
413 .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
414 .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
417 sql = sa.lambda_stmt(_place_inside_area_query)
418 if self.has_geometries():
419 sql = self._add_geometry_columns(sql, sa.literal_column('places.geometry'))
421 place_address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
422 log().var_dump('Result (place node)', place_address_row)
424 if place_address_row is not None:
425 return place_address_row
430 async def _lookup_area_others(self) -> Optional[SaRow]:
431 t = self.conn.t.placex
433 inner = sa.select(t, t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
434 .where(t.c.rank_address == 0)\
435 .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
436 .where(t.c.name != None)\
437 .where(t.c.indexed_status == 0)\
438 .where(t.c.linked_place_id == None)\
439 .where(self._filter_by_layer(t))\
441 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
442 .intersects(WKT_PARAM))\
443 .order_by(sa.desc(t.c.rank_search))\
447 sql = _select_from_placex(inner, False)\
448 .where(sa.or_(sa.not_(inner.c.geometry.is_area()),
449 inner.c.geometry.ST_Contains(WKT_PARAM)))\
450 .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
453 if self.has_geometries():
454 sql = self._add_geometry_columns(sql, inner.c.geometry)
456 row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
457 log().var_dump('Result (non-address feature)', row)
462 async def lookup_area(self) -> Optional[SaRow]:
463 """ Lookup large areas for the current search.
465 log().section('Reverse lookup by larger area features')
467 if self.layer_enabled(DataLayer.ADDRESS):
468 address_row = await self._lookup_area_address()
472 if self.has_feature_layers():
473 other_row = await self._lookup_area_others()
477 return _get_closest(address_row, other_row)
480 async def lookup_country(self) -> Optional[SaRow]:
481 """ Lookup the country for the current search.
483 log().section('Reverse lookup by country code')
484 t = self.conn.t.country_grid
485 sql: SaLambdaSelect = sa.select(t.c.country_code).distinct()\
486 .where(t.c.geometry.ST_Contains(WKT_PARAM))
488 ccodes = tuple((r[0] for r in await self.conn.execute(sql, self.bind_params)))
489 log().var_dump('Country codes', ccodes)
494 t = self.conn.t.placex
495 if self.max_rank > 4:
496 log().comment('Search for place nodes in country')
498 def _base_query() -> SaSelect:
501 t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
502 .where(t.c.rank_search > 4)\
503 .where(t.c.rank_search <= MAX_RANK_PARAM)\
504 .where(t.c.indexed_status == 0)\
505 .where(t.c.country_code.in_(ccodes))\
506 .where(snfn.select_index_placex_geometry_reverse_lookupplacenode('placex'))\
508 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
509 .intersects(WKT_PARAM))\
510 .order_by(sa.desc(t.c.rank_search))\
514 return _select_from_placex(inner, False)\
515 .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
516 .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
519 sql = sa.lambda_stmt(_base_query)
520 if self.has_geometries():
521 sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
523 address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
524 log().var_dump('Result (addressable place node)', address_row)
528 if address_row is None:
529 # Still nothing, then return a country with the appropriate country code.
530 sql = sa.lambda_stmt(lambda: _select_from_placex(t)\
531 .where(t.c.country_code.in_(ccodes))\
532 .where(t.c.rank_address == 4)\
533 .where(t.c.rank_search == 4)\
534 .where(t.c.linked_place_id == None)\
535 .order_by('distance')\
538 if self.has_geometries():
539 sql = self._add_geometry_columns(sql, t.c.geometry)
541 address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
546 async def lookup(self, coord: AnyPoint) -> Optional[nres.ReverseResult]:
547 """ Look up a single coordinate. Returns the place information,
548 if a place was found near the coordinates or None otherwise.
550 log().function('reverse_lookup', coord=coord, params=self.params)
553 self.bind_params['wkt'] = f'POINT({coord[0]} {coord[1]})'
555 row: Optional[SaRow] = None
556 row_func: RowFunc = nres.create_from_placex_row
558 if self.max_rank >= 26:
559 row, tmp_row_func = await self.lookup_street_poi()
561 row_func = tmp_row_func
562 if row is None and self.max_rank > 4:
563 row = await self.lookup_area()
564 if row is None and self.layer_enabled(DataLayer.ADDRESS):
565 row = await self.lookup_country()
567 result = row_func(row, nres.ReverseResult)
568 if result is not None:
569 assert row is not None
570 result.distance = row.distance
571 if hasattr(row, 'bbox'):
572 result.bbox = Bbox.from_wkb(row.bbox)
573 await nres.add_result_details(self.conn, [result], self.params)