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