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