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