]> git.openstreetmap.org Git - nominatim.git/blob - nominatim/api/reverse.py
78e3a20c4900827781cfa97eac8c37b0833c57d0
[nominatim.git] / nominatim / api / reverse.py
1 # SPDX-License-Identifier: GPL-3.0-or-later
2 #
3 # This file is part of Nominatim. (https://nominatim.org)
4 #
5 # Copyright (C) 2023 by the Nominatim developer community.
6 # For a full list of authors see the git log.
7 """
8 Implementation of reverse geocoding.
9 """
10 from typing import Optional, List, Callable, Type, Tuple, Dict, Any, cast, Union
11 import functools
12
13 import sqlalchemy as sa
14
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
23 # In SQLAlchemy expression which compare with NULL need to be expressed with
24 # the equal sign.
25 # pylint: disable=singleton-comparison
26
27 RowFunc = Callable[[Optional[SaRow], Type[nres.ReverseResult]], Optional[nres.ReverseResult]]
28
29 WKT_PARAM: SaBind = sa.bindparam('wkt', type_=Geometry)
30 MAX_RANK_PARAM: SaBind = sa.bindparam('max_rank')
31
32 def no_index(expr: SaColumn) -> SaColumn:
33     """ Wrap the given expression, so that the query planner will
34         refrain from using the expression for index lookup.
35     """
36     return sa.func.coalesce(sa.null(), expr) # pylint: disable=not-callable
37
38
39 def _select_from_placex(t: SaFromClause, use_wkt: bool = True) -> SaSelect:
40     """ Create a select statement with the columns relevant for reverse
41         results.
42     """
43     if not use_wkt:
44         distance = t.c.distance
45         centroid = t.c.centroid
46     else:
47         distance = t.c.geometry.ST_Distance(WKT_PARAM)
48         centroid = sa.case((t.c.geometry.is_line_like(), t.c.geometry.ST_ClosestPoint(WKT_PARAM)),
49                            else_=t.c.centroid).label('centroid')
50
51
52     return sa.select(t.c.place_id, t.c.osm_type, t.c.osm_id, t.c.name,
53                      t.c.class_, t.c.type,
54                      t.c.address, t.c.extratags,
55                      t.c.housenumber, t.c.postcode, t.c.country_code,
56                      t.c.importance, t.c.wikipedia,
57                      t.c.parent_place_id, t.c.rank_address, t.c.rank_search,
58                      centroid,
59                      t.c.linked_place_id, t.c.admin_level,
60                      distance.label('distance'),
61                      t.c.geometry.ST_Expand(0).label('bbox'))
62
63
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')
69
70
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
74     return sa.case(
75              (table.c.endnumber == table.c.startnumber, table.c.linegeo.ST_Centroid()),
76               else_=table.c.linegeo.ST_LineInterpolatePoint(rounded_pos)).label('centroid')
77
78
79 def _locate_interpolation(table: SaFromClause) -> SaLabel:
80     """ Given a position, locate the closest point on the line.
81     """
82     return sa.case((table.c.linegeo.is_line_like(),
83                     table.c.linegeo.ST_LineLocatePoint(WKT_PARAM)),
84                    else_=0).label('position')
85
86
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                           sa.func.JsonHasKey(table.c.name, 'addr:housename')))
91
92
93 def _get_closest(*rows: Optional[SaRow]) -> Optional[SaRow]:
94     return min(rows, key=lambda row: 1000 if row is None else row.distance)
95
96
97 class ReverseGeocoder:
98     """ Class implementing the logic for looking up a place from a
99         coordinate.
100     """
101
102     def __init__(self, conn: SearchConnection, params: ReverseDetails,
103                  restrict_to_country_areas: bool = False) -> None:
104         self.conn = conn
105         self.params = params
106         self.restrict_to_country_areas = restrict_to_country_areas
107
108         self.bind_params: Dict[str, Any] = {'max_rank': params.max_rank}
109
110
111     @property
112     def max_rank(self) -> int:
113         """ Return the maximum configured rank.
114         """
115         return self.params.max_rank
116
117
118     def has_geometries(self) -> bool:
119         """ Check if any geometries are requested.
120         """
121         return bool(self.params.geometry_output)
122
123
124     def layer_enabled(self, *layer: DataLayer) -> bool:
125         """ Return true when any of the given layer types are requested.
126         """
127         return any(self.params.layers & l for l in layer)
128
129
130     def layer_disabled(self, *layer: DataLayer) -> bool:
131         """ Return true when none of the given layer types is requested.
132         """
133         return not any(self.params.layers & l for l in layer)
134
135
136     def has_feature_layers(self) -> bool:
137         """ Return true if any layer other than ADDRESS or POI is requested.
138         """
139         return self.layer_enabled(DataLayer.RAILWAY, DataLayer.MANMADE, DataLayer.NATURAL)
140
141
142     def _add_geometry_columns(self, sql: SaLambdaSelect, col: SaColumn) -> SaSelect:
143         out = []
144
145         if self.params.geometry_simplification > 0.0:
146             col = sa.func.ST_SimplifyPreserveTopology(col, self.params.geometry_simplification)
147
148         if self.params.geometry_output & GeometryFormat.GEOJSON:
149             out.append(sa.func.ST_AsGeoJSON(col, 7).label('geometry_geojson'))
150         if self.params.geometry_output & GeometryFormat.TEXT:
151             out.append(sa.func.ST_AsText(col).label('geometry_text'))
152         if self.params.geometry_output & GeometryFormat.KML:
153             out.append(sa.func.ST_AsKML(col, 7).label('geometry_kml'))
154         if self.params.geometry_output & GeometryFormat.SVG:
155             out.append(sa.func.ST_AsSVG(col, 0, 7).label('geometry_svg'))
156
157         return sql.add_columns(*out)
158
159
160     def _filter_by_layer(self, table: SaFromClause) -> SaColumn:
161         if self.layer_enabled(DataLayer.MANMADE):
162             exclude = []
163             if self.layer_disabled(DataLayer.RAILWAY):
164                 exclude.append('railway')
165             if self.layer_disabled(DataLayer.NATURAL):
166                 exclude.extend(('natural', 'water', 'waterway'))
167             return table.c.class_.not_in(tuple(exclude))
168
169         include = []
170         if self.layer_enabled(DataLayer.RAILWAY):
171             include.append('railway')
172         if self.layer_enabled(DataLayer.NATURAL):
173             include.extend(('natural', 'water', 'waterway'))
174         return table.c.class_.in_(tuple(include))
175
176
177     async def _find_closest_street_or_poi(self, distance: float) -> Optional[SaRow]:
178         """ Look up the closest rank 26+ place in the database, which
179             is closer than the given distance.
180         """
181         t = self.conn.t.placex
182
183         # PostgreSQL must not get the distance as a parameter because
184         # there is a danger it won't be able to proberly estimate index use
185         # when used with prepared statements
186         diststr = sa.text(f"{distance}")
187
188         sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
189                 .where(t.c.geometry.ST_DWithin(WKT_PARAM, diststr))
190                 .where(t.c.indexed_status == 0)
191                 .where(t.c.linked_place_id == None)
192                 .where(sa.or_(sa.not_(t.c.geometry.is_area()),
193                               t.c.centroid.ST_Distance(WKT_PARAM) < diststr))
194                 .order_by('distance')
195                 .limit(1))
196
197         if self.has_geometries():
198             sql = self._add_geometry_columns(sql, t.c.geometry)
199
200         restrict: List[Union[SaColumn, Callable[[], SaColumn]]] = []
201
202         if self.layer_enabled(DataLayer.ADDRESS):
203             max_rank = min(29, self.max_rank)
204             restrict.append(lambda: no_index(t.c.rank_address).between(26, max_rank))
205             if self.max_rank == 30:
206                 restrict.append(lambda: _is_address_point(t))
207         if self.layer_enabled(DataLayer.POI) and self.max_rank == 30:
208             restrict.append(lambda: sa.and_(no_index(t.c.rank_search) == 30,
209                                             t.c.class_.not_in(('place', 'building')),
210                                             sa.not_(t.c.geometry.is_line_like())))
211         if self.has_feature_layers():
212             restrict.append(sa.and_(no_index(t.c.rank_search).between(26, MAX_RANK_PARAM),
213                                     no_index(t.c.rank_address) == 0,
214                                     self._filter_by_layer(t)))
215
216         if not restrict:
217             return None
218
219         sql = sql.where(sa.or_(*restrict))
220
221         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
222
223
224     async def _find_housenumber_for_street(self, parent_place_id: int) -> Optional[SaRow]:
225         t = self.conn.t.placex
226
227         sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
228                 .where(t.c.geometry.ST_DWithin(WKT_PARAM, 0.001))
229                 .where(t.c.parent_place_id == parent_place_id)
230                 .where(_is_address_point(t))
231                 .where(t.c.indexed_status == 0)
232                 .where(t.c.linked_place_id == None)
233                 .order_by('distance')
234                 .limit(1))
235
236         if self.has_geometries():
237             sql = self._add_geometry_columns(sql, t.c.geometry)
238
239         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
240
241
242     async def _find_interpolation_for_street(self, parent_place_id: Optional[int],
243                                              distance: float) -> Optional[SaRow]:
244         t = self.conn.t.osmline
245
246         sql: Any = sa.lambda_stmt(lambda:
247                    sa.select(t,
248                              t.c.linegeo.ST_Distance(WKT_PARAM).label('distance'),
249                              _locate_interpolation(t))
250                      .where(t.c.linegeo.ST_DWithin(WKT_PARAM, distance))
251                      .where(t.c.startnumber != None)
252                      .order_by('distance')
253                      .limit(1))
254
255         if parent_place_id is not None:
256             sql += lambda s: s.where(t.c.parent_place_id == parent_place_id)
257
258         def _wrap_query(base_sql: SaLambdaSelect) -> SaSelect:
259             inner = base_sql.subquery('ipol')
260
261             return sa.select(inner.c.place_id, inner.c.osm_id,
262                              inner.c.parent_place_id, inner.c.address,
263                              _interpolated_housenumber(inner),
264                              _interpolated_position(inner),
265                              inner.c.postcode, inner.c.country_code,
266                              inner.c.distance)
267
268         sql += _wrap_query
269
270         if self.has_geometries():
271             sub = sql.subquery('geom')
272             sql = self._add_geometry_columns(sa.select(sub), sub.c.centroid)
273
274         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
275
276
277     async def _find_tiger_number_for_street(self, parent_place_id: int) -> Optional[SaRow]:
278         t = self.conn.t.tiger
279
280         def _base_query() -> SaSelect:
281             inner = sa.select(t,
282                               t.c.linegeo.ST_Distance(WKT_PARAM).label('distance'),
283                               _locate_interpolation(t))\
284                       .where(t.c.linegeo.ST_DWithin(WKT_PARAM, 0.001))\
285                       .where(t.c.parent_place_id == parent_place_id)\
286                       .order_by('distance')\
287                       .limit(1)\
288                       .subquery('tiger')
289
290             return sa.select(inner.c.place_id,
291                              inner.c.parent_place_id,
292                              _interpolated_housenumber(inner),
293                              _interpolated_position(inner),
294                              inner.c.postcode,
295                              inner.c.distance)
296
297         sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
298
299         if self.has_geometries():
300             sub = sql.subquery('geom')
301             sql = self._add_geometry_columns(sa.select(sub), sub.c.centroid)
302
303         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
304
305
306     async def lookup_street_poi(self) -> Tuple[Optional[SaRow], RowFunc]:
307         """ Find a street or POI/address for the given WKT point.
308         """
309         log().section('Reverse lookup on street/address level')
310         distance = 0.006
311         parent_place_id = None
312
313         row = await self._find_closest_street_or_poi(distance)
314         row_func: RowFunc = nres.create_from_placex_row
315         log().var_dump('Result (street/building)', row)
316
317         # If the closest result was a street, but an address was requested,
318         # check for a housenumber nearby which is part of the street.
319         if row is not None:
320             if self.max_rank > 27 \
321                and self.layer_enabled(DataLayer.ADDRESS) \
322                and row.rank_address <= 27:
323                 distance = 0.001
324                 parent_place_id = row.place_id
325                 log().comment('Find housenumber for street')
326                 addr_row = await self._find_housenumber_for_street(parent_place_id)
327                 log().var_dump('Result (street housenumber)', addr_row)
328
329                 if addr_row is not None:
330                     row = addr_row
331                     row_func = nres.create_from_placex_row
332                     distance = addr_row.distance
333                 elif row.country_code == 'us' and parent_place_id is not None:
334                     log().comment('Find TIGER housenumber for street')
335                     addr_row = await self._find_tiger_number_for_street(parent_place_id)
336                     log().var_dump('Result (street Tiger housenumber)', addr_row)
337
338                     if addr_row is not None:
339                         row_func = cast(RowFunc,
340                                         functools.partial(nres.create_from_tiger_row,
341                                                           osm_type=row.osm_type,
342                                                           osm_id=row.osm_id))
343                         row = addr_row
344             else:
345                 distance = row.distance
346
347         # Check for an interpolation that is either closer than our result
348         # or belongs to a close street found.
349         if self.max_rank > 27 and self.layer_enabled(DataLayer.ADDRESS):
350             log().comment('Find interpolation for street')
351             addr_row = await self._find_interpolation_for_street(parent_place_id,
352                                                                  distance)
353             log().var_dump('Result (street interpolation)', addr_row)
354             if addr_row is not None:
355                 row = addr_row
356                 row_func = nres.create_from_osmline_row
357
358         return row, row_func
359
360
361     async def _lookup_area_address(self) -> Optional[SaRow]:
362         """ Lookup large addressable areas for the given WKT point.
363         """
364         log().comment('Reverse lookup by larger address area features')
365         t = self.conn.t.placex
366
367         def _base_query() -> SaSelect:
368             # The inner SQL brings results in the right order, so that
369             # later only a minimum of results needs to be checked with ST_Contains.
370             inner = sa.select(t, sa.literal(0.0).label('distance'))\
371                       .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
372                       .where(t.c.geometry.intersects(WKT_PARAM))\
373                       .where(sa.func.PlacexGeometryReverseLookuppolygon())\
374                       .order_by(sa.desc(t.c.rank_search))\
375                       .limit(50)\
376                       .subquery('area')
377
378             return _select_from_placex(inner, False)\
379                       .where(inner.c.geometry.ST_Contains(WKT_PARAM))\
380                       .order_by(sa.desc(inner.c.rank_search))\
381                       .limit(1)
382
383         sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
384         if self.has_geometries():
385             sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
386
387         address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
388         log().var_dump('Result (area)', address_row)
389
390         if address_row is not None and address_row.rank_search < self.max_rank:
391             log().comment('Search for better matching place nodes inside the area')
392
393             address_rank = address_row.rank_search
394             address_id = address_row.place_id
395
396             def _place_inside_area_query() -> SaSelect:
397                 inner = \
398                     sa.select(t,
399                               t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
400                       .where(t.c.rank_search > address_rank)\
401                       .where(t.c.rank_search <= MAX_RANK_PARAM)\
402                       .where(t.c.indexed_status == 0)\
403                       .where(sa.func.IntersectsReverseDistance(t, WKT_PARAM))\
404                       .order_by(sa.desc(t.c.rank_search))\
405                       .limit(50)\
406                       .subquery('places')
407
408                 touter = t.alias('outer')
409                 return _select_from_placex(inner, False)\
410                     .join(touter, touter.c.geometry.ST_Contains(inner.c.geometry))\
411                     .where(touter.c.place_id == address_id)\
412                     .where(sa.func.IsBelowReverseDistance(inner.c.distance, inner.c.rank_search))\
413                     .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
414                     .limit(1)
415
416             sql = sa.lambda_stmt(_place_inside_area_query)
417             if self.has_geometries():
418                 sql = self._add_geometry_columns(sql, sa.literal_column('places.geometry'))
419
420             place_address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
421             log().var_dump('Result (place node)', place_address_row)
422
423             if place_address_row is not None:
424                 return place_address_row
425
426         return address_row
427
428
429     async def _lookup_area_others(self) -> Optional[SaRow]:
430         t = self.conn.t.placex
431
432         inner = sa.select(t, t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
433                   .where(t.c.rank_address == 0)\
434                   .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
435                   .where(t.c.name != None)\
436                   .where(t.c.indexed_status == 0)\
437                   .where(t.c.linked_place_id == None)\
438                   .where(self._filter_by_layer(t))\
439                   .where(t.c.geometry.intersects(sa.func.ST_Expand(WKT_PARAM, 0.007)))\
440                   .order_by(sa.desc(t.c.rank_search))\
441                   .order_by('distance')\
442                   .limit(50)\
443                   .subquery()
444
445         sql = _select_from_placex(inner, False)\
446                   .where(sa.or_(sa.not_(inner.c.geometry.is_area()),
447                                 inner.c.geometry.ST_Contains(WKT_PARAM)))\
448                   .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
449                   .limit(1)
450
451         if self.has_geometries():
452             sql = self._add_geometry_columns(sql, inner.c.geometry)
453
454         row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
455         log().var_dump('Result (non-address feature)', row)
456
457         return row
458
459
460     async def lookup_area(self) -> Optional[SaRow]:
461         """ Lookup large areas for the current search.
462         """
463         log().section('Reverse lookup by larger area features')
464
465         if self.layer_enabled(DataLayer.ADDRESS):
466             address_row = await self._lookup_area_address()
467         else:
468             address_row = None
469
470         if self.has_feature_layers():
471             other_row = await self._lookup_area_others()
472         else:
473             other_row = None
474
475         return _get_closest(address_row, other_row)
476
477
478     async def lookup_country_codes(self) -> List[str]:
479         """ Lookup the country for the current search.
480         """
481         log().section('Reverse lookup by country code')
482         t = self.conn.t.country_grid
483         sql = sa.select(t.c.country_code).distinct()\
484                 .where(t.c.geometry.ST_Contains(WKT_PARAM))
485
486         ccodes = [cast(str, r[0]) for r in await self.conn.execute(sql, self.bind_params)]
487         log().var_dump('Country codes', ccodes)
488         return ccodes
489
490
491     async def lookup_country(self, ccodes: List[str]) -> Optional[SaRow]:
492         """ Lookup the country for the current search.
493         """
494         if not ccodes:
495             ccodes = await self.lookup_country_codes()
496
497         if not ccodes:
498             return None
499
500         t = self.conn.t.placex
501         if self.max_rank > 4:
502             log().comment('Search for place nodes in country')
503
504             def _base_query() -> SaSelect:
505                 inner = \
506                     sa.select(t,
507                               t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
508                       .where(t.c.rank_search > 4)\
509                       .where(t.c.rank_search <= MAX_RANK_PARAM)\
510                       .where(t.c.indexed_status == 0)\
511                       .where(t.c.country_code.in_(ccodes))\
512                       .where(sa.func.IntersectsReverseDistance(t, WKT_PARAM))\
513                       .order_by(sa.desc(t.c.rank_search))\
514                       .limit(50)\
515                       .subquery('area')
516
517                 return _select_from_placex(inner, False)\
518                     .where(sa.func.IsBelowReverseDistance(inner.c.distance, inner.c.rank_search))\
519                     .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
520                     .limit(1)
521
522             sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
523             if self.has_geometries():
524                 sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
525
526             address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
527             log().var_dump('Result (addressable place node)', address_row)
528         else:
529             address_row = None
530
531         if address_row is None:
532             # Still nothing, then return a country with the appropriate country code.
533             sql = sa.lambda_stmt(lambda: _select_from_placex(t)\
534                       .where(t.c.country_code.in_(ccodes))\
535                       .where(t.c.rank_address == 4)\
536                       .where(t.c.rank_search == 4)\
537                       .where(t.c.linked_place_id == None)\
538                       .order_by('distance')\
539                       .limit(1))
540
541             if self.has_geometries():
542                 sql = self._add_geometry_columns(sql, t.c.geometry)
543
544             address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
545
546         return address_row
547
548
549     async def lookup(self, coord: AnyPoint) -> Optional[nres.ReverseResult]:
550         """ Look up a single coordinate. Returns the place information,
551             if a place was found near the coordinates or None otherwise.
552         """
553         log().function('reverse_lookup', coord=coord, params=self.params)
554
555
556         self.bind_params['wkt'] = f'POINT({coord[0]} {coord[1]})'
557
558         row: Optional[SaRow] = None
559         row_func: RowFunc = nres.create_from_placex_row
560
561         if self.max_rank >= 26:
562             row, tmp_row_func = await self.lookup_street_poi()
563             if row is not None:
564                 row_func = tmp_row_func
565
566         if row is None:
567             if self.restrict_to_country_areas:
568                 ccodes = await self.lookup_country_codes()
569                 if not ccodes:
570                     return None
571             else:
572                 ccodes = []
573
574             if self.max_rank > 4:
575                 row = await self.lookup_area()
576             if row is None and self.layer_enabled(DataLayer.ADDRESS):
577                 row = await self.lookup_country(ccodes)
578
579         result = row_func(row, nres.ReverseResult)
580         if result is not None:
581             assert row is not None
582             result.distance = row.distance
583             if hasattr(row, 'bbox'):
584                 result.bbox = Bbox.from_wkb(row.bbox)
585             await nres.add_result_details(self.conn, [result], self.params)
586
587         return result