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