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