64 lines
2.4 KiB
Python
64 lines
2.4 KiB
Python
import math
|
|
from dataclasses import dataclass
|
|
|
|
from shapely.geometry import Polygon, box
|
|
|
|
EARTH_RADIUS = 6_371_008.8
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class Bounds:
|
|
west: float
|
|
south: float
|
|
east: float
|
|
north: float
|
|
|
|
def __post_init__(self) -> None:
|
|
if not all(math.isfinite(v) for v in self.values):
|
|
raise ValueError("Map coordinates must be finite.")
|
|
if not (-180 <= self.west < self.east <= 180):
|
|
raise ValueError("Longitude bounds must increase within [-180, 180].")
|
|
if self.east - self.west > 180:
|
|
raise ValueError("Antimeridian-crossing maps are not supported.")
|
|
if not (-85 <= self.south < self.north <= 85):
|
|
raise ValueError("Latitude bounds must increase within [-85, 85].")
|
|
|
|
@property
|
|
def values(self) -> tuple[float, float, float, float]:
|
|
return self.west, self.south, self.east, self.north
|
|
|
|
@classmethod
|
|
def corners(cls, lon1: float, lat1: float, lon2: float, lat2: float) -> "Bounds":
|
|
return cls(min(lon1, lon2), min(lat1, lat2), max(lon1, lon2), max(lat1, lat2))
|
|
|
|
def expand(self, metres: float) -> "Bounds":
|
|
if not math.isfinite(metres) or metres <= 0:
|
|
raise ValueError("Collection margin must be a positive number of metres.")
|
|
dy = math.degrees(metres / EARTH_RADIUS)
|
|
dx = dy / math.cos(math.radians((self.south + self.north) / 2))
|
|
return Bounds(self.west - dx, self.south - dy, self.east + dx, self.north + dy)
|
|
|
|
def contains(self, lon: float, lat: float) -> bool:
|
|
return self.west <= lon <= self.east and self.south <= lat <= self.north
|
|
|
|
|
|
DEFAULT_BOUNDS = Bounds(4.784, 52.232, 4.864, 52.258)
|
|
|
|
|
|
class Projection:
|
|
"""Local equirectangular projection: geographic bounds remain rectangular."""
|
|
|
|
def __init__(self, bounds: Bounds):
|
|
self.lon0 = (bounds.west + bounds.east) / 2
|
|
self.lat0 = (bounds.south + bounds.north) / 2
|
|
self.x_factor = EARTH_RADIUS * math.cos(math.radians(self.lat0))
|
|
west, south = self.point(bounds.west, bounds.south)
|
|
east, north = self.point(bounds.east, bounds.north)
|
|
self.rectangle: Polygon = box(west, south, east, north)
|
|
|
|
def point(self, lon: float, lat: float) -> tuple[float, float]:
|
|
return (
|
|
self.x_factor * math.radians(lon - self.lon0),
|
|
EARTH_RADIUS * math.radians(lat - self.lat0),
|
|
)
|