"""Bounding box operations for geohashes.
This module provides functions for working with geohash bounding boxes,
including calculating the bounding box for a geohash and operations
related to geospatial regions.
"""
from __future__ import annotations
import math
from typing import Any, Iterable, List, NamedTuple, Optional, Tuple, Type, TypeVar
from pygeohash.geohash import decode_exactly, encode
from pygeohash.logging import get_logger
logger = get_logger(__name__)
_FIELD_ORDER = "BoundingBox fields are (min_lat, min_lon, max_lat, max_lon)"
def _validate_point_coordinates(lat: float, lon: float) -> None:
"""Validate coordinates supplied to point-containment operations."""
for name, value, lower, upper in (
("latitude", lat, -90.0, 90.0),
("longitude", lon, -180.0, 180.0),
):
if isinstance(value, bool):
raise ValueError(f"{name} ({value}) must be a number, not a bool")
try:
is_finite = math.isfinite(value)
except TypeError:
is_finite = False
if not is_finite:
raise ValueError(f"{name} ({value}) must be a finite number")
try:
is_in_range = lower <= value <= upper
except TypeError:
is_in_range = False
if not is_in_range:
raise ValueError(f"{name} ({value}) must be between {lower:g} and {upper:g}")
class _BoundingBoxFields(NamedTuple):
"""Field layout for :class:`BoundingBox`; ``typing.NamedTuple`` forbids overriding ``__new__``."""
min_lat: float
min_lon: float
max_lat: float
max_lon: float
_BoundingBoxT = TypeVar("_BoundingBoxT", bound=_BoundingBoxFields)
[docs]
class BoundingBox(_BoundingBoxFields):
"""Named tuple representing a geospatial bounding box.
The fields interleave latitude and longitude, so the order is
``(min_lat, min_lon, max_lat, max_lon)`` rather than the grouped
``(min_lat, max_lat, min_lon, max_lon)``. Coordinates must be finite numbers within
the geographic bounds for their axis; booleans are rejected even though ``bool`` is a
subclass of ``int``. Construction also rejects an inverted box
(``min_lat > max_lat`` or ``min_lon > max_lon``) with a ``ValueError``, which is
what a grouped argument list produces. A degenerate box whose minimum equals its
maximum on either axis is valid. Boxes spanning the antimeridian, which would need
``min_lon > max_lon``, are not supported.
Attributes:
min_lat (float): The minimum (southern) latitude of the box in decimal degrees.
Must be between -90 and 90 and not exceed ``max_lat``.
min_lon (float): The minimum (western) longitude of the box in decimal degrees.
Must be between -180 and 180 and not exceed ``max_lon``.
max_lat (float): The maximum (northern) latitude of the box in decimal degrees,
between -90 and 90.
max_lon (float): The maximum (eastern) longitude of the box in decimal degrees,
between -180 and 180.
Raises:
ValueError: If a coordinate is a boolean, non-finite, outside its geographic
bounds, or the box has ``min_lat > max_lat`` or ``min_lon > max_lon``.
"""
__slots__ = ()
def __new__(cls, min_lat: float, min_lon: float, max_lat: float, max_lon: float) -> "BoundingBox":
for field, value, lower, upper in (
("min_lat", min_lat, -90.0, 90.0),
("min_lon", min_lon, -180.0, 180.0),
("max_lat", max_lat, -90.0, 90.0),
("max_lon", max_lon, -180.0, 180.0),
):
# bool is a subclass of int, so it has to be rejected explicitly.
if isinstance(value, bool):
raise ValueError(f"{field} ({value}) must be a number, not a bool")
try:
is_finite = math.isfinite(value)
except TypeError:
is_finite = False
if not is_finite:
raise ValueError(f"{field} ({value}) must be a finite number")
if not lower <= value <= upper:
raise ValueError(f"{field} ({value}) must be between {lower:g} and {upper:g}")
if min_lat > max_lat:
raise ValueError(f"min_lat ({min_lat}) must not exceed max_lat ({max_lat}); {_FIELD_ORDER}")
if min_lon > max_lon:
raise ValueError(
f"min_lon ({min_lon}) must not exceed max_lon ({max_lon}); "
"boxes spanning the antimeridian are not supported"
)
return super().__new__(cls, min_lat, min_lon, max_lat, max_lon)
# mypy rejects the narrowed self-type against the synthesized namedtuple signature.
@classmethod
def _make(cls: Type[_BoundingBoxT], iterable: Iterable[Any]) -> _BoundingBoxT: # type: ignore[override]
"""Build a box from an iterable of field values, validating the ordering.
Overridden so that ``_make`` and ``_replace`` route through ``__new__`` instead of
``tuple.__new__``, which would bypass validation.
"""
return cls(*iterable)
@classmethod
def _unvalidated(cls, min_lat: float, min_lon: float, max_lat: float, max_lon: float) -> "BoundingBox":
"""Build a box from values the library itself derived, skipping validation."""
return tuple.__new__(cls, (min_lat, min_lon, max_lat, max_lon))
[docs]
def get_bounding_box(geohash: str) -> BoundingBox:
"""Calculate the bounding box for a geohash.
Args:
geohash (str): The geohash string to calculate the bounding box for.
Returns:
BoundingBox: A named tuple containing the minimum and maximum latitude and longitude
values that define the bounding box of the geohash.
Example:
>>> tuple(round(value, 6) for value in get_bounding_box("u4pruyd"))
(57.64801, 10.406799, 57.649384, 10.408173)
Note:
The precision of the coordinates in the bounding box depends on the length
of the geohash. Longer geohashes result in smaller bounding boxes with more
precise coordinates.
"""
# Get the center point and error margins
lat, lon, lat_err, lon_err = decode_exactly(geohash)
return BoundingBox._unvalidated(lat - lat_err, lon - lon_err, lat + lat_err, lon + lon_err)
[docs]
def is_point_in_box(lat: float, lon: float, bbox: BoundingBox) -> bool:
"""Check if a point is within a bounding box.
Args:
lat (float): The latitude of the point to check.
lon (float): The longitude of the point to check.
bbox (BoundingBox): The bounding box to check against.
Returns:
bool: True if the point is within the bounding box, False otherwise.
Raises:
ValueError: If latitude or longitude is a boolean, non-numeric, non-finite,
or outside its geographic bounds.
Example:
>>> bbox = get_bounding_box("u4pruyd")
>>> is_point_in_box(57.649, 10.407, bbox)
True
>>> is_point_in_box(40.0, 10.0, bbox)
False
"""
_validate_point_coordinates(lat, lon)
logger.debug("Checking if point (lat=%f, lon=%f) is in box: %s", lat, lon, bbox)
result = bbox.min_lat <= lat <= bbox.max_lat and bbox.min_lon <= lon <= bbox.max_lon
logger.debug("Point is %s the box", "inside" if result else "outside")
return result
[docs]
def is_point_in_geohash(lat: float, lon: float, geohash: str) -> bool:
"""Check if a point is within a geohash's bounding box.
Args:
lat (float): The latitude of the point to check.
lon (float): The longitude of the point to check.
geohash (str): The geohash to check against.
Returns:
bool: True if the point is within the geohash's bounding box, False otherwise.
Raises:
ValueError: If latitude or longitude is a boolean, non-numeric, non-finite,
or outside its geographic bounds, or if ``geohash`` is invalid.
Example:
>>> is_point_in_geohash(57.649, 10.407, "u4pruyd")
True
>>> is_point_in_geohash(40.0, 10.0, "u4pruyd")
False
"""
logger.debug("Checking if point (lat=%f, lon=%f) is in geohash: %s", lat, lon, geohash)
bbox: BoundingBox = get_bounding_box(geohash)
return is_point_in_box(lat, lon, bbox)
[docs]
def do_boxes_intersect(bbox1: BoundingBox, bbox2: BoundingBox) -> bool:
"""Check if two bounding boxes intersect.
Args:
bbox1 (BoundingBox): The first bounding box.
bbox2 (BoundingBox): The second bounding box.
Returns:
bool: True if the bounding boxes intersect, False otherwise.
Example:
>>> box1 = BoundingBox(10.0, 20.0, 30.0, 40.0)
>>> box2 = BoundingBox(20.0, 30.0, 40.0, 50.0)
>>> do_boxes_intersect(box1, box2)
True
"""
logger.debug("Checking intersection between boxes: %s and %s", bbox1, bbox2)
result = not (
bbox1.max_lat < bbox2.min_lat
or bbox1.min_lat > bbox2.max_lat
or bbox1.max_lon < bbox2.min_lon
or bbox1.min_lon > bbox2.max_lon
)
logger.debug("Boxes %s intersect", "do" if result else "do not")
return result
DEFAULT_MAX_CELLS: int = 5_000_000
def _validate_max_cells(max_cells: Optional[int]) -> None:
if max_cells is None:
return
if isinstance(max_cells, bool) or not isinstance(max_cells, int):
raise ValueError("max_cells must be an integer or None, not a bool or other type")
if max_cells < 1:
raise ValueError("max_cells must be at least 1")
def _check_cell_count(estimated: int, max_cells: Optional[int], precision: int) -> None:
if max_cells is not None and estimated > max_cells:
raise ValueError(
f"Region would enumerate about {estimated} candidate cells at precision {precision}, "
f"which exceeds max_cells={max_cells}. Lower the precision, use a smaller region, "
"or pass max_cells=None to disable the limit."
)
def _box_lattice(bbox: BoundingBox, precision: int) -> Tuple[int, int, int, int, float, float, float, float]:
"""Return ``(first_row, last_row, first_col, last_col, lat_step, lon_step, lat_err, lon_err)``."""
center_lat: float = (bbox.min_lat + bbox.max_lat) / 2
center_lon: float = (bbox.min_lon + bbox.max_lon) / 2
center_geohash: str = encode(center_lat, center_lon, precision)
logger.debug("Center geohash: %s at (lat=%f, lon=%f)", center_geohash, center_lat, center_lon)
_cell_lat, _cell_lon, lat_err, lon_err = decode_exactly(center_geohash)
lat_step: float = 2.0 * lat_err
lon_step: float = 2.0 * lon_err
lat_rows: int = 1 << ((5 * precision) // 2)
lon_cols: int = 1 << ((5 * precision + 1) // 2)
first_row: int = max(0, math.floor((bbox.min_lat + 90.0) / lat_step) - 1)
last_row: int = min(lat_rows - 1, math.floor((bbox.max_lat + 90.0) / lat_step) + 1)
first_col: int = max(0, math.floor((bbox.min_lon + 180.0) / lon_step) - 1)
last_col: int = min(lon_cols - 1, math.floor((bbox.max_lon + 180.0) / lon_step) + 1)
return first_row, last_row, first_col, last_col, lat_step, lon_step, lat_err, lon_err
def _lattice_size(lattice: Tuple[int, int, int, int, float, float, float, float]) -> int:
first_row, last_row, first_col, last_col = lattice[:4]
return (last_row - first_row + 1) * (last_col - first_col + 1)
[docs]
def geohashes_in_box(
bbox: BoundingBox, precision: int = 6, *, max_cells: Optional[int] = DEFAULT_MAX_CELLS
) -> List[str]:
"""Find geohashes that intersect with a given bounding box.
Args:
bbox (BoundingBox): The bounding box to find geohashes for.
precision (int, optional): The precision of the geohashes to return. Defaults to 6.
max_cells (Optional[int]): Keyword-only upper bound on the number of grid
cells the enumeration may visit, estimated from the covering row and
column counts before any allocation. Defaults to 5,000,000; ``None``
disables the check.
Returns:
List[str]: A sorted list of geohashes that intersect with the bounding box.
Example:
>>> box = BoundingBox(57.64, 10.40, 57.65, 10.41)
>>> geohashes_in_box(box, precision=5)
['u4pru']
Raises:
ValueError: If ``precision`` or ``max_cells`` is invalid, or if the estimated
cell count exceeds ``max_cells``; lower the precision, shrink the box, or
pass ``max_cells=None``.
Note:
The number of geohashes returned depends on the size of the bounding box
and the precision requested. Higher precision values will result in more
geohashes for the same bounding box. Cells are enumerated directly on
the geohash grid instead of sampling points, so the result is returned
in a deterministically sorted order that is identical across processes.
"""
_validate_max_cells(max_cells)
logger.debug("Finding geohashes in box %s with precision %d", bbox, precision)
# Building the lattice encodes the box center first, which validates the
# precision argument exactly as the previous sampling implementation did.
lattice = _box_lattice(bbox, precision)
_check_cell_count(_lattice_size(lattice), max_cells, precision)
first_row, last_row, first_col, last_col, lat_step, lon_step, lat_err, lon_err = lattice
# A cell is classified interior (fully inside the box) or exterior (clear
# of it) only when the decision clears this margin on both axes; anything
# closer is tested exactly with do_boxes_intersect. Decoded cell bounds are
# exact multiples of the cell size (the C decoder halves intervals exactly)
# and -90 + row * lat_step reproduces the same values exactly, so 1e-12
# degrees is orders of magnitude above the arithmetic slack while staying
# far below half a cell at precision 12 (~8e-8 degrees).
margin: float = 1e-12
min_lat: float = bbox.min_lat
max_lat: float = bbox.max_lat
min_lon: float = bbox.min_lon
max_lon: float = bbox.max_lon
# Column data is row-independent: precompute each candidate column's
# center longitude and interior flag once, so the row loop below only
# does latitude arithmetic per cell.
col_centers: List[float] = []
col_interior: List[bool] = []
for col in range(first_col, last_col + 1):
cell_min_lon: float = -180.0 + col * lon_step
cell_max_lon: float = cell_min_lon + lon_step
if cell_max_lon < min_lon - margin or cell_min_lon > max_lon + margin:
continue # column entirely outside the box
col_centers.append(cell_min_lon + lon_err)
col_interior.append(cell_min_lon >= min_lon + margin and cell_max_lon <= max_lon - margin)
result: List[str] = []
append = result.append
for row in range(first_row, last_row + 1):
cell_min_lat: float = -90.0 + row * lat_step
cell_max_lat: float = cell_min_lat + lat_step
if cell_max_lat < min_lat - margin or cell_min_lat > max_lat + margin:
continue # row entirely outside the box
row_center_lat: float = cell_min_lat + lat_err
lat_interior: bool = cell_min_lat >= min_lat + margin and cell_max_lat <= max_lat - margin
for center_lon, lon_interior in zip(col_centers, col_interior, strict=True):
cell_geohash: str = encode(row_center_lat, center_lon, precision)
if lat_interior and lon_interior:
append(cell_geohash) # interior cell: fully inside the box
elif do_boxes_intersect(bbox, get_bounding_box(cell_geohash)):
append(cell_geohash)
logger.debug("Found %d intersecting geohashes", len(result))
return sorted(result)
def _distance_to_box(latitude: float, longitude: float, bbox: BoundingBox) -> float:
"""Return the spherical distance to the closest point in a cell."""
from pygeohash.types import EARTH_RADIUS
if latitude == 90.0:
return EARTH_RADIUS * math.radians(90.0 - bbox.max_lat)
if latitude == -90.0:
return EARTH_RADIUS * math.radians(bbox.min_lat + 90.0)
if bbox.min_lon <= longitude <= bbox.max_lon:
delta_lon = 0.0
else:
delta_lon = min(abs((edge - longitude + 180.0) % 360.0 - 180.0) for edge in (bbox.min_lon, bbox.max_lon))
phi = math.radians(latitude)
delta = math.radians(delta_lon)
closest_lat = math.degrees(math.atan2(math.sin(phi), math.cos(phi) * math.cos(delta)))
candidates = (bbox.min_lat, bbox.max_lat, max(bbox.min_lat, min(bbox.max_lat, closest_lat)))
distances = []
for lat in candidates:
target_phi = math.radians(lat)
a = math.sin((target_phi - phi) / 2.0) ** 2 + math.cos(phi) * math.cos(target_phi) * math.sin(delta / 2.0) ** 2
distances.append(2.0 * EARTH_RADIUS * math.asin(math.sqrt(max(0.0, min(1.0, a)))))
return min(distances)
[docs]
def geohashes_in_radius(
latitude: float,
longitude: float,
radius_m: float,
precision: int = 6,
*,
max_cells: Optional[int] = DEFAULT_MAX_CELLS,
) -> List[str]:
"""Find all geohash cells intersecting a spherical circle.
Args:
latitude (float): Center latitude in degrees, between -90 and 90.
longitude (float): Center longitude in degrees, between -180 and 180.
radius_m (float): Finite, positive radius in meters. Booleans are rejected.
precision (int): Cell precision from 1 to 12. Defaults to 6.
max_cells (Optional[int]): Keyword-only upper bound on the total number of
grid cells the enumeration may visit, summed over antimeridian
intervals and estimated before any allocation. Defaults to 5,000,000;
``None`` disables the check.
Returns:
List[str]: Sorted, unique cells whose nearest point is within the radius,
including cells touching the circle. Antimeridian crossings and polar caps
are supported. Distances use the spherical ``EARTH_RADIUS`` constant.
Raises:
ValueError: If coordinates, radius, precision, or ``max_cells`` are invalid,
or if the estimated cell count exceeds ``max_cells``. A circle touching a
pole spans every longitude, so small radii at high precision can exceed
the limit; lower the precision or pass ``max_cells=None``.
Note:
These are candidate cells for proximity queries; filter individual records
by their distance afterwards. Output size grows with radius and precision;
large circles at high precision can require very large lists.
"""
from pygeohash.types import EARTH_RADIUS
_validate_point_coordinates(latitude, longitude)
if isinstance(precision, bool) or not isinstance(precision, int):
raise ValueError("Precision must be an integer, not a bool")
if not 1 <= precision <= 12:
raise ValueError("Precision must be between 1 and 12")
if isinstance(radius_m, bool):
raise ValueError("radius_m must be a number, not a bool")
try:
finite_radius = math.isfinite(radius_m)
except TypeError:
finite_radius = False
if not finite_radius:
raise ValueError("radius_m must be a finite number")
if radius_m <= 0:
raise ValueError("radius_m must be positive")
_validate_max_cells(max_cells)
angular_radius = min(math.pi, radius_m / EARTH_RADIUS)
lat_span = math.degrees(angular_radius)
min_lat = max(-90.0, latitude - lat_span)
max_lat = min(90.0, latitude + lat_span)
if lat_span >= 90.0 - abs(latitude):
intervals = [(-180.0, 180.0)]
else:
lon_span = math.degrees(math.asin(min(1.0, math.sin(angular_radius) / math.cos(math.radians(latitude)))))
west, east = longitude - lon_span, longitude + lon_span
if west <= -180.0:
intervals = [(-180.0, east), (west + 360.0, 180.0)]
elif east >= 180.0:
intervals = [(west, 180.0), (-180.0, east - 360.0)]
else:
intervals = [(west, east)]
boxes = [BoundingBox._unvalidated(min_lat, west, max_lat, east) for west, east in intervals]
_check_cell_count(sum(_lattice_size(_box_lattice(box, precision)) for box in boxes), max_cells, precision)
candidates = set()
for box in boxes:
candidates.update(geohashes_in_box(box, precision, max_cells=None))
return sorted(
cell for cell in candidates if _distance_to_box(latitude, longitude, get_bounding_box(cell)) <= radius_m
)