import math
from typing import List, Sequence

from .agri_field_area_types import FieldArea
from .geo_point_in_polygon_types import GeoPoint
from .math_round_float import round_float
from .math_sin_cos import sin_cos

RADIUS = 6371008.8
RADIANS = 3.141592653589793 / 180
ACRE = 4046.8564224


def _number(value: object) -> bool:
    return isinstance(value, (int, float)) and not isinstance(value, bool) and math.isfinite(value)


def field_area(boundary: Sequence[GeoPoint]) -> FieldArea:
    """Area enclosed by a boundary, by the shoelace formula in the Lambert
    cylindrical equal-area projection (x = R lng, y = R sin lat), which is the
    spherical polygon formula of Chamberlain and Duquette (2007).

    Longitudes are unwrapped from the first point and sines are taken relative
    to the first point's, which keeps a one-hectare field accurate when the
    coordinates are large; sin comes from math.sin-cos so every language adds
    the same doubles in the same order.
    """
    for i, p in enumerate(boundary):
        if not _number(p.lat) or p.lat < -90 or p.lat > 90:
            raise ValueError("point %d: latitude must be a finite number from -90 to 90 degrees" % (i + 1,))
        if not _number(p.lng) or p.lng < -180 or p.lng > 180:
            raise ValueError("point %d: longitude must be a finite number from -180 to 180 degrees" % (i + 1,))
    ring = list(boundary)
    if len(ring) > 1 and ring[0].lat == ring[-1].lat and ring[0].lng == ring[-1].lng:
        ring = ring[:-1]
    if len({(float(p.lat), float(p.lng)) for p in ring}) < 3:
        raise ValueError("a boundary needs at least 3 distinct points")

    n = len(ring)
    x: List[float] = []
    y: List[float] = []
    sin0 = sin_cos(float(ring[0].lat) * RADIANS).sin
    previous = 0.0
    for p in ring:
        d = float(p.lng) - float(ring[0].lng)
        # Take the short way round, so a field across the antimeridian works.
        while d - previous > 180:
            d -= 360
        while d - previous < -180:
            d += 360
        previous = d
        x.append(d * RADIANS)
        y.append(sin_cos(float(p.lat) * RADIANS).sin - sin0)
    total = 0.0
    for i in range(n):
        total += (x[(i + 1) % n] - x[(i + n - 1) % n]) * y[i]
    area = abs(total) * RADIUS * RADIUS / 2
    return FieldArea(
        square_metres=round_float(area, 0),
        hectares=round_float(area / 10000, 4),
        acres=round_float(area / ACRE, 4),
    )