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