]> git.openstreetmap.org Git - nominatim.git/blob - nominatim/api/reverse.py
38f993be549f12b997dc1d90a39e53d44cd2af9e
[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
11
12 import sqlalchemy as sa
13
14 from nominatim.typing import SaColumn, SaSelect, SaFromClause, SaLabel, SaRow,\
15                              SaBind, SaLambdaSelect
16 from nominatim.api.connection import SearchConnection
17 import nominatim.api.results as nres
18 from nominatim.api.logging import log
19 from nominatim.api.types import AnyPoint, DataLayer, ReverseDetails, GeometryFormat, Bbox
20 from nominatim.db.sqlalchemy_types import Geometry
21 import nominatim.db.sqlalchemy_functions as snfn
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                      distance.label('distance'),
60                      t.c.geometry.ST_Expand(0).label('bbox'))
61
62
63 def _interpolated_housenumber(table: SaFromClause) -> SaLabel:
64     return sa.cast(table.c.startnumber
65                     + sa.func.round(((table.c.endnumber - table.c.startnumber) * table.c.position)
66                                     / table.c.step) * table.c.step,
67                    sa.Integer).label('housenumber')
68
69
70 def _interpolated_position(table: SaFromClause) -> SaLabel:
71     fac = sa.cast(table.c.step, sa.Float) / (table.c.endnumber - table.c.startnumber)
72     rounded_pos = sa.func.round(table.c.position / fac) * fac
73     return sa.case(
74              (table.c.endnumber == table.c.startnumber, table.c.linegeo.ST_Centroid()),
75               else_=table.c.linegeo.ST_LineInterpolatePoint(rounded_pos)).label('centroid')
76
77
78 def _locate_interpolation(table: SaFromClause) -> SaLabel:
79     """ Given a position, locate the closest point on the line.
80     """
81     return sa.case((table.c.linegeo.is_line_like(),
82                     table.c.linegeo.ST_LineLocatePoint(WKT_PARAM)),
83                    else_=0).label('position')
84
85
86 def _is_address_point(table: SaFromClause) -> SaColumn:
87     return sa.and_(table.c.rank_address == 30,
88                    sa.or_(table.c.housenumber != None,
89                           table.c.name.has_key('housename')))
90
91
92 def _get_closest(*rows: Optional[SaRow]) -> Optional[SaRow]:
93     return min(rows, key=lambda row: 1000 if row is None else row.distance)
94
95
96 class ReverseGeocoder:
97     """ Class implementing the logic for looking up a place from a
98         coordinate.
99     """
100
101     def __init__(self, conn: SearchConnection, params: ReverseDetails) -> None:
102         self.conn = conn
103         self.params = params
104
105         self.bind_params: Dict[str, Any] = {'max_rank': params.max_rank}
106
107
108     @property
109     def max_rank(self) -> int:
110         """ Return the maximum configured rank.
111         """
112         return self.params.max_rank
113
114
115     def has_geometries(self) -> bool:
116         """ Check if any geometries are requested.
117         """
118         return bool(self.params.geometry_output)
119
120
121     def layer_enabled(self, *layer: DataLayer) -> bool:
122         """ Return true when any of the given layer types are requested.
123         """
124         return any(self.params.layers & l for l in layer)
125
126
127     def layer_disabled(self, *layer: DataLayer) -> bool:
128         """ Return true when none of the given layer types is requested.
129         """
130         return not any(self.params.layers & l for l in layer)
131
132
133     def has_feature_layers(self) -> bool:
134         """ Return true if any layer other than ADDRESS or POI is requested.
135         """
136         return self.layer_enabled(DataLayer.RAILWAY, DataLayer.MANMADE, DataLayer.NATURAL)
137
138
139     def _add_geometry_columns(self, sql: SaLambdaSelect, col: SaColumn) -> SaSelect:
140         out = []
141
142         if self.params.geometry_simplification > 0.0:
143             col = sa.func.ST_SimplifyPreserveTopology(col, self.params.geometry_simplification)
144
145         if self.params.geometry_output & GeometryFormat.GEOJSON:
146             out.append(sa.func.ST_AsGeoJSON(col).label('geometry_geojson'))
147         if self.params.geometry_output & GeometryFormat.TEXT:
148             out.append(sa.func.ST_AsText(col).label('geometry_text'))
149         if self.params.geometry_output & GeometryFormat.KML:
150             out.append(sa.func.ST_AsKML(col).label('geometry_kml'))
151         if self.params.geometry_output & GeometryFormat.SVG:
152             out.append(sa.func.ST_AsSVG(col).label('geometry_svg'))
153
154         return sql.add_columns(*out)
155
156
157     def _filter_by_layer(self, table: SaFromClause) -> SaColumn:
158         if self.layer_enabled(DataLayer.MANMADE):
159             exclude = []
160             if self.layer_disabled(DataLayer.RAILWAY):
161                 exclude.append('railway')
162             if self.layer_disabled(DataLayer.NATURAL):
163                 exclude.extend(('natural', 'water', 'waterway'))
164             return table.c.class_.not_in(tuple(exclude))
165
166         include = []
167         if self.layer_enabled(DataLayer.RAILWAY):
168             include.append('railway')
169         if self.layer_enabled(DataLayer.NATURAL):
170             include.extend(('natural', 'water', 'waterway'))
171         return table.c.class_.in_(tuple(include))
172
173
174     async def _find_closest_street_or_poi(self, distance: float) -> Optional[SaRow]:
175         """ Look up the closest rank 26+ place in the database, which
176             is closer than the given distance.
177         """
178         t = self.conn.t.placex
179
180         # PostgreSQL must not get the distance as a parameter because
181         # there is a danger it won't be able to proberly estimate index use
182         # when used with prepared statements
183         diststr = sa.text(f"{distance}")
184
185         sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
186                 .where(t.c.geometry.ST_DWithin(WKT_PARAM, diststr))
187                 .where(t.c.indexed_status == 0)
188                 .where(t.c.linked_place_id == None)
189                 .where(sa.or_(sa.not_(t.c.geometry.is_area()),
190                               t.c.centroid.ST_Distance(WKT_PARAM) < diststr))
191                 .order_by('distance')
192                 .limit(1))
193
194         if self.has_geometries():
195             sql = self._add_geometry_columns(sql, t.c.geometry)
196
197         restrict: List[SaColumn] = []
198
199         if self.layer_enabled(DataLayer.ADDRESS):
200             restrict.append(no_index(t.c.rank_address).between(26, min(29, self.max_rank)))
201             if self.max_rank == 30:
202                 restrict.append(_is_address_point(t))
203         if self.layer_enabled(DataLayer.POI) and self.max_rank == 30:
204             restrict.append(sa.and_(no_index(t.c.rank_search) == 30,
205                                     t.c.class_.not_in(('place', 'building')),
206                                     sa.not_(t.c.geometry.is_line_like())))
207         if self.has_feature_layers():
208             restrict.append(sa.and_(no_index(t.c.rank_search).between(26, MAX_RANK_PARAM),
209                                     no_index(t.c.rank_address) == 0,
210                                     self._filter_by_layer(t)))
211
212         if not restrict:
213             return None
214
215         sql = sql.where(sa.or_(*restrict))
216
217         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
218
219
220     async def _find_housenumber_for_street(self, parent_place_id: int) -> Optional[SaRow]:
221         t = self.conn.t.placex
222
223         sql: SaLambdaSelect = sa.lambda_stmt(lambda: _select_from_placex(t)
224                 .where(t.c.geometry.ST_DWithin(WKT_PARAM, 0.001))
225                 .where(t.c.parent_place_id == parent_place_id)
226                 .where(_is_address_point(t))
227                 .where(t.c.indexed_status == 0)
228                 .where(t.c.linked_place_id == None)
229                 .order_by('distance')
230                 .limit(1))
231
232         if self.has_geometries():
233             sql = self._add_geometry_columns(sql, t.c.geometry)
234
235         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
236
237
238     async def _find_interpolation_for_street(self, parent_place_id: Optional[int],
239                                              distance: float) -> Optional[SaRow]:
240         t = self.conn.t.osmline
241
242         sql: Any = sa.lambda_stmt(lambda:
243                    sa.select(t,
244                              t.c.linegeo.ST_Distance(WKT_PARAM).label('distance'),
245                              _locate_interpolation(t))
246                      .where(t.c.linegeo.ST_DWithin(WKT_PARAM, distance))
247                      .where(t.c.startnumber != None)
248                      .order_by('distance')
249                      .limit(1))
250
251         if parent_place_id is not None:
252             sql += lambda s: s.where(t.c.parent_place_id == parent_place_id)
253
254         def _wrap_query(base_sql: SaLambdaSelect) -> SaSelect:
255             inner = base_sql.subquery('ipol')
256
257             return sa.select(inner.c.place_id, inner.c.osm_id,
258                              inner.c.parent_place_id, inner.c.address,
259                              _interpolated_housenumber(inner),
260                              _interpolated_position(inner),
261                              inner.c.postcode, inner.c.country_code,
262                              inner.c.distance)
263
264         sql += _wrap_query
265
266         if self.has_geometries():
267             sub = sql.subquery('geom')
268             sql = self._add_geometry_columns(sa.select(sub), sub.c.centroid)
269
270         return (await self.conn.execute(sql, self.bind_params)).one_or_none()
271
272
273     async def _find_tiger_number_for_street(self, parent_place_id: int,
274                                             parent_type: str,
275                                             parent_id: int) -> Optional[SaRow]:
276         t = self.conn.t.tiger
277
278         def _base_query() -> SaSelect:
279             inner = sa.select(t,
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')\
285                       .limit(1)\
286                       .subquery('tiger')
287
288             return sa.select(inner.c.place_id,
289                              inner.c.parent_place_id,
290                              sa.sql.expression.label('osm_type', parent_type),
291                              sa.sql.expression.label('osm_id', parent_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                                                                         row.osm_type,
337                                                                         row.osm_id)
338                     log().var_dump('Result (street Tiger housenumber)', addr_row)
339
340                     if addr_row is not None:
341                         row = addr_row
342                         row_func = nres.create_from_tiger_row
343             else:
344                 distance = row.distance
345
346         # Check for an interpolation that is either closer than our result
347         # or belongs to a close street found.
348         if self.max_rank > 27 and self.layer_enabled(DataLayer.ADDRESS):
349             log().comment('Find interpolation for street')
350             addr_row = await self._find_interpolation_for_street(parent_place_id,
351                                                                  distance)
352             log().var_dump('Result (street interpolation)', addr_row)
353             if addr_row is not None:
354                 row = addr_row
355                 row_func = nres.create_from_osmline_row
356
357         return row, row_func
358
359
360     async def _lookup_area_address(self) -> Optional[SaRow]:
361         """ Lookup large addressable areas for the given WKT point.
362         """
363         log().comment('Reverse lookup by larger address area features')
364         t = self.conn.t.placex
365
366         def _base_query() -> SaSelect:
367             # The inner SQL brings results in the right order, so that
368             # later only a minimum of results needs to be checked with ST_Contains.
369             inner = sa.select(t, sa.literal(0.0).label('distance'))\
370                       .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
371                       .where(t.c.geometry.intersects(WKT_PARAM))\
372                       .where(snfn.select_index_placex_geometry_reverse_lookuppolygon('placex'))\
373                       .order_by(sa.desc(t.c.rank_search))\
374                       .limit(50)\
375                       .subquery('area')
376
377             return _select_from_placex(inner, False)\
378                       .where(inner.c.geometry.ST_Contains(WKT_PARAM))\
379                       .order_by(sa.desc(inner.c.rank_search))\
380                       .limit(1)
381
382         sql: SaLambdaSelect = sa.lambda_stmt(_base_query)
383         if self.has_geometries():
384             sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
385
386         address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
387         log().var_dump('Result (area)', address_row)
388
389         if address_row is not None and address_row.rank_search < self.max_rank:
390             log().comment('Search for better matching place nodes inside the area')
391
392             address_rank = address_row.rank_search
393             address_id = address_row.place_id
394
395             def _place_inside_area_query() -> SaSelect:
396                 inner = \
397                     sa.select(t,
398                               t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
399                       .where(t.c.rank_search > address_rank)\
400                       .where(t.c.rank_search <= MAX_RANK_PARAM)\
401                       .where(t.c.indexed_status == 0)\
402                       .where(snfn.select_index_placex_geometry_reverse_lookupplacenode('placex'))\
403                       .where(t.c.geometry
404                                 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
405                                 .intersects(WKT_PARAM))\
406                       .order_by(sa.desc(t.c.rank_search))\
407                       .limit(50)\
408                       .subquery('places')
409
410                 touter = t.alias('outer')
411                 return _select_from_placex(inner, False)\
412                     .join(touter, touter.c.geometry.ST_Contains(inner.c.geometry))\
413                     .where(touter.c.place_id == address_id)\
414                     .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
415                     .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
416                     .limit(1)
417
418             sql = sa.lambda_stmt(_place_inside_area_query)
419             if self.has_geometries():
420                 sql = self._add_geometry_columns(sql, sa.literal_column('places.geometry'))
421
422             place_address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
423             log().var_dump('Result (place node)', place_address_row)
424
425             if place_address_row is not None:
426                 return place_address_row
427
428         return address_row
429
430
431     async def _lookup_area_others(self) -> Optional[SaRow]:
432         t = self.conn.t.placex
433
434         inner = sa.select(t, t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
435                   .where(t.c.rank_address == 0)\
436                   .where(t.c.rank_search.between(5, MAX_RANK_PARAM))\
437                   .where(t.c.name != None)\
438                   .where(t.c.indexed_status == 0)\
439                   .where(t.c.linked_place_id == None)\
440                   .where(self._filter_by_layer(t))\
441                   .where(t.c.geometry
442                                 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
443                                 .intersects(WKT_PARAM))\
444                   .order_by(sa.desc(t.c.rank_search))\
445                   .limit(50)\
446                   .subquery()
447
448         sql = _select_from_placex(inner, False)\
449                   .where(sa.or_(sa.not_(inner.c.geometry.is_area()),
450                                 inner.c.geometry.ST_Contains(WKT_PARAM)))\
451                   .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
452                   .limit(1)
453
454         if self.has_geometries():
455             sql = self._add_geometry_columns(sql, inner.c.geometry)
456
457         row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
458         log().var_dump('Result (non-address feature)', row)
459
460         return row
461
462
463     async def lookup_area(self) -> Optional[SaRow]:
464         """ Lookup large areas for the current search.
465         """
466         log().section('Reverse lookup by larger area features')
467
468         if self.layer_enabled(DataLayer.ADDRESS):
469             address_row = await self._lookup_area_address()
470         else:
471             address_row = None
472
473         if self.has_feature_layers():
474             other_row = await self._lookup_area_others()
475         else:
476             other_row = None
477
478         return _get_closest(address_row, other_row)
479
480
481     async def lookup_country(self) -> Optional[SaRow]:
482         """ Lookup the country for the current search.
483         """
484         log().section('Reverse lookup by country code')
485         t = self.conn.t.country_grid
486         sql: SaLambdaSelect = sa.select(t.c.country_code).distinct()\
487                 .where(t.c.geometry.ST_Contains(WKT_PARAM))
488
489         ccodes = tuple((r[0] for r in await self.conn.execute(sql, self.bind_params)))
490         log().var_dump('Country codes', ccodes)
491
492         if not ccodes:
493             return None
494
495         t = self.conn.t.placex
496         if self.max_rank > 4:
497             log().comment('Search for place nodes in country')
498
499             def _base_query() -> SaSelect:
500                 inner = \
501                     sa.select(t,
502                               t.c.geometry.ST_Distance(WKT_PARAM).label('distance'))\
503                       .where(t.c.rank_search > 4)\
504                       .where(t.c.rank_search <= MAX_RANK_PARAM)\
505                       .where(t.c.indexed_status == 0)\
506                       .where(t.c.country_code.in_(ccodes))\
507                       .where(snfn.select_index_placex_geometry_reverse_lookupplacenode('placex'))\
508                       .where(t.c.geometry
509                                 .ST_Buffer(sa.func.reverse_place_diameter(t.c.rank_search))
510                                 .intersects(WKT_PARAM))\
511                       .order_by(sa.desc(t.c.rank_search))\
512                       .limit(50)\
513                       .subquery('area')
514
515                 return _select_from_placex(inner, False)\
516                     .where(inner.c.distance < sa.func.reverse_place_diameter(inner.c.rank_search))\
517                     .order_by(sa.desc(inner.c.rank_search), inner.c.distance)\
518                     .limit(1)
519
520             sql = sa.lambda_stmt(_base_query)
521             if self.has_geometries():
522                 sql = self._add_geometry_columns(sql, sa.literal_column('area.geometry'))
523
524             address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
525             log().var_dump('Result (addressable place node)', address_row)
526         else:
527             address_row = None
528
529         if address_row is None:
530             # Still nothing, then return a country with the appropriate country code.
531             sql = sa.lambda_stmt(lambda: _select_from_placex(t)\
532                       .where(t.c.country_code.in_(ccodes))\
533                       .where(t.c.rank_address == 4)\
534                       .where(t.c.rank_search == 4)\
535                       .where(t.c.linked_place_id == None)\
536                       .order_by('distance')\
537                       .limit(1))
538
539             if self.has_geometries():
540                 sql = self._add_geometry_columns(sql, t.c.geometry)
541
542             address_row = (await self.conn.execute(sql, self.bind_params)).one_or_none()
543
544         return address_row
545
546
547     async def lookup(self, coord: AnyPoint) -> Optional[nres.ReverseResult]:
548         """ Look up a single coordinate. Returns the place information,
549             if a place was found near the coordinates or None otherwise.
550         """
551         log().function('reverse_lookup', coord=coord, params=self.params)
552
553
554         self.bind_params['wkt'] = f'POINT({coord[0]} {coord[1]})'
555
556         row: Optional[SaRow] = None
557         row_func: RowFunc = nres.create_from_placex_row
558
559         if self.max_rank >= 26:
560             row, tmp_row_func = await self.lookup_street_poi()
561             if row is not None:
562                 row_func = tmp_row_func
563         if row is None and self.max_rank > 4:
564             row = await self.lookup_area()
565         if row is None and self.layer_enabled(DataLayer.ADDRESS):
566             row = await self.lookup_country()
567
568         result = row_func(row, nres.ReverseResult)
569         if result is not None:
570             assert row is not None
571             result.distance = row.distance
572             if hasattr(row, 'bbox'):
573                 result.bbox = Bbox.from_wkb(row.bbox)
574             await nres.add_result_details(self.conn, [result], self.params)
575
576         return result