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),
        )