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