rustial_engine/

geometry_ops.rs

// ---------------------------------------------------------------------------
//! # Geometry operations -- spatial computations on geographic geometries
//!
//! Provides core spatial analysis functions that operate on the engine's
//! [`Geometry`] / [`GeoCoord`] types:
//!
//! - **Bounding box** -- [`geometry_bbox`], [`Geometry::bbox`]
//! - **Centroid** -- [`geometry_centroid`], [`Geometry::centroid`]
//! - **Area** -- [`polygon_area`], [`ring_area`], [`Geometry::area`]
//! - **Length** -- [`linestring_length`], [`Geometry::length`]
//! - **Bearing** -- [`bearing`], [`initial_bearing`]
//! - **Interpolation** -- [`interpolate_great_circle`]
//!
//! All distance/area calculations use the Haversine formula on the
//! WGS-84 sphere (radius 6,378,137 m) for a good balance of speed
//! and accuracy.
// ---------------------------------------------------------------------------

use crate::geometry::{Geometry, Polygon};
use rustial_math::{GeoBounds, GeoCoord};

/// Mean radius of the Earth in meters (WGS-84 semi-major axis).
const EARTH_RADIUS: f64 = 6_378_137.0;

// ---------------------------------------------------------------------------
// Bounding box
// ---------------------------------------------------------------------------

/// Compute the geographic bounding box of a geometry.
///
/// Returns `None` for empty geometries.
pub fn geometry_bbox(geometry: &Geometry) -> Option<GeoBounds> {
    let mut min_lat = f64::INFINITY;
    let mut max_lat = f64::NEG_INFINITY;
    let mut min_lon = f64::INFINITY;
    let mut max_lon = f64::NEG_INFINITY;

    fn visit(
        coord: &GeoCoord,
        min_lat: &mut f64,
        max_lat: &mut f64,
        min_lon: &mut f64,
        max_lon: &mut f64,
    ) {
        if coord.lat < *min_lat { *min_lat = coord.lat; }
        if coord.lat > *max_lat { *max_lat = coord.lat; }
        if coord.lon < *min_lon { *min_lon = coord.lon; }
        if coord.lon > *max_lon { *max_lon = coord.lon; }
    }

    fn visit_coords(
        coords: &[GeoCoord],
        min_lat: &mut f64,
        max_lat: &mut f64,
        min_lon: &mut f64,
        max_lon: &mut f64,
    ) {
        for c in coords {
            visit(c, min_lat, max_lat, min_lon, max_lon);
        }
    }

    fn visit_geometry(
        g: &Geometry,
        min_lat: &mut f64,
        max_lat: &mut f64,
        min_lon: &mut f64,
        max_lon: &mut f64,
    ) {
        match g {
            Geometry::Point(p) => visit(&p.coord, min_lat, max_lat, min_lon, max_lon),
            Geometry::LineString(ls) => visit_coords(&ls.coords, min_lat, max_lat, min_lon, max_lon),
            Geometry::Polygon(p) => {
                visit_coords(&p.exterior, min_lat, max_lat, min_lon, max_lon);
                for hole in &p.interiors {
                    visit_coords(hole, min_lat, max_lat, min_lon, max_lon);
                }
            }
            Geometry::MultiPoint(mp) => {
                for p in &mp.points {
                    visit(&p.coord, min_lat, max_lat, min_lon, max_lon);
                }
            }
            Geometry::MultiLineString(mls) => {
                for ls in &mls.lines {
                    visit_coords(&ls.coords, min_lat, max_lat, min_lon, max_lon);
                }
            }
            Geometry::MultiPolygon(mp) => {
                for p in &mp.polygons {
                    visit_coords(&p.exterior, min_lat, max_lat, min_lon, max_lon);
                }
            }
            Geometry::GeometryCollection(gc) => {
                for g in gc {
                    visit_geometry(g, min_lat, max_lat, min_lon, max_lon);
                }
            }
        }
    }

    visit_geometry(geometry, &mut min_lat, &mut max_lat, &mut min_lon, &mut max_lon);

    if min_lat.is_infinite() {
        None
    } else {
        Some(GeoBounds::from_coords(min_lon, min_lat, max_lon, max_lat))
    }
}

// ---------------------------------------------------------------------------
// Area
// ---------------------------------------------------------------------------

/// Signed area of a ring on a sphere using the Haversine-based excess formula.
///
/// Positive for counter-clockwise rings, negative for clockwise.
/// The result is in **square meters** on the WGS-84 sphere.
pub fn ring_area(coords: &[GeoCoord]) -> f64 {
    let n = coords.len();
    if n < 3 {
        return 0.0;
    }

    // Use the spherical excess formula:
    //   A = R² * |Σ (λ[i+1] - λ[i-1]) * sin(φ[i])|
    // This is the planar-approximation formula that works well for small polygons.
    // For a more rigorous approach we use the shoelace formula in radians
    // on the sphere.
    let mut area = 0.0;
    for i in 0..n {
        let j = if i == 0 { n - 1 } else { i - 1 };
        let k = if i == n - 1 { 0 } else { i + 1 };
        let lon_j = coords[j].lon.to_radians();
        let lon_k = coords[k].lon.to_radians();
        let lat_i = coords[i].lat.to_radians();
        area += (lon_k - lon_j) * lat_i.sin();
    }
    area * EARTH_RADIUS * EARTH_RADIUS * 0.5
}

/// Area of a polygon (exterior minus holes) in **square meters**.
///
/// Always returns a non-negative value.
pub fn polygon_area(polygon: &Polygon) -> f64 {
    let mut area = ring_area(&polygon.exterior).abs();
    for hole in &polygon.interiors {
        area -= ring_area(hole).abs();
    }
    area.abs()
}

// ---------------------------------------------------------------------------
// Length
// ---------------------------------------------------------------------------

/// Haversine distance between two geographic coordinates in meters.
pub(crate) fn haversine(a: &GeoCoord, b: &GeoCoord) -> f64 {
    let dlat = (b.lat - a.lat).to_radians();
    let dlon = (b.lon - a.lon).to_radians();
    let lat1 = a.lat.to_radians();
    let lat2 = b.lat.to_radians();
    let h = (dlat * 0.5).sin().powi(2) + lat1.cos() * lat2.cos() * (dlon * 0.5).sin().powi(2);
    2.0 * EARTH_RADIUS * h.sqrt().asin()
}

/// Total length of a linestring in meters (Haversine).
pub fn linestring_length(coords: &[GeoCoord]) -> f64 {
    if coords.len() < 2 {
        return 0.0;
    }
    coords.windows(2).map(|w| haversine(&w[0], &w[1])).sum()
}

// ---------------------------------------------------------------------------
// Centroid
// ---------------------------------------------------------------------------

/// Centroid of a geometry.
///
/// - **Point** -- the point itself.
/// - **LineString** -- midpoint weighted by segment length.
/// - **Polygon** -- area-weighted centroid of the exterior ring.
/// - **Multi*** -- weighted average of component centroids.
/// - **GeometryCollection** -- unweighted mean of non-empty children.
///
/// Returns `None` for empty geometries.
pub fn geometry_centroid(geometry: &Geometry) -> Option<GeoCoord> {
    match geometry {
        Geometry::Point(p) => Some(p.coord),
        Geometry::LineString(ls) => linestring_centroid(&ls.coords),
        Geometry::Polygon(p) => ring_centroid(&p.exterior),
        Geometry::MultiPoint(mp) => {
            if mp.points.is_empty() {
                return None;
            }
            let (sum_lat, sum_lon) = mp.points.iter().fold((0.0, 0.0), |(lat, lon), p| {
                (lat + p.coord.lat, lon + p.coord.lon)
            });
            let n = mp.points.len() as f64;
            Some(GeoCoord::from_lat_lon(sum_lat / n, sum_lon / n))
        }
        Geometry::MultiLineString(mls) => {
            let mut total_len = 0.0;
            let mut wlat = 0.0;
            let mut wlon = 0.0;
            for ls in &mls.lines {
                let len = linestring_length(&ls.coords);
                if let Some(c) = linestring_centroid(&ls.coords) {
                    wlat += c.lat * len;
                    wlon += c.lon * len;
                    total_len += len;
                }
            }
            if total_len == 0.0 {
                return None;
            }
            Some(GeoCoord::from_lat_lon(wlat / total_len, wlon / total_len))
        }
        Geometry::MultiPolygon(mp) => {
            let mut total_area = 0.0;
            let mut wlat = 0.0;
            let mut wlon = 0.0;
            for poly in &mp.polygons {
                let area = polygon_area(poly);
                if let Some(c) = ring_centroid(&poly.exterior) {
                    wlat += c.lat * area;
                    wlon += c.lon * area;
                    total_area += area;
                }
            }
            if total_area == 0.0 {
                return None;
            }
            Some(GeoCoord::from_lat_lon(wlat / total_area, wlon / total_area))
        }
        Geometry::GeometryCollection(gc) => {
            let centroids: Vec<GeoCoord> = gc.iter().filter_map(geometry_centroid).collect();
            if centroids.is_empty() {
                return None;
            }
            let n = centroids.len() as f64;
            let (slat, slon) = centroids.iter().fold((0.0, 0.0), |(lat, lon), c| {
                (lat + c.lat, lon + c.lon)
            });
            Some(GeoCoord::from_lat_lon(slat / n, slon / n))
        }
    }
}

/// Centroid of a linestring weighted by segment lengths.
fn linestring_centroid(coords: &[GeoCoord]) -> Option<GeoCoord> {
    if coords.is_empty() {
        return None;
    }
    if coords.len() == 1 {
        return Some(coords[0]);
    }
    let mut total_len = 0.0;
    let mut wlat = 0.0;
    let mut wlon = 0.0;
    for w in coords.windows(2) {
        let len = haversine(&w[0], &w[1]);
        let mid_lat = (w[0].lat + w[1].lat) * 0.5;
        let mid_lon = (w[0].lon + w[1].lon) * 0.5;
        wlat += mid_lat * len;
        wlon += mid_lon * len;
        total_len += len;
    }
    if total_len == 0.0 {
        return Some(coords[0]);
    }
    Some(GeoCoord::from_lat_lon(wlat / total_len, wlon / total_len))
}

/// Centroid of a closed ring (signed-area weighted).
fn ring_centroid(coords: &[GeoCoord]) -> Option<GeoCoord> {
    let n = coords.len();
    if n < 3 {
        return None;
    }
    let mut cx = 0.0;
    let mut cy = 0.0;
    let mut signed_area = 0.0;
    for i in 0..n {
        let j = (i + 1) % n;
        let xi = coords[i].lon;
        let yi = coords[i].lat;
        let xj = coords[j].lon;
        let yj = coords[j].lat;
        let a = xi * yj - xj * yi;
        signed_area += a;
        cx += (xi + xj) * a;
        cy += (yi + yj) * a;
    }
    if signed_area.abs() < 1e-20 {
        return None;
    }
    signed_area *= 0.5;
    cx /= 6.0 * signed_area;
    cy /= 6.0 * signed_area;
    Some(GeoCoord::from_lat_lon(cy, cx))
}

// ---------------------------------------------------------------------------
// Bearing
// ---------------------------------------------------------------------------

/// Initial bearing (forward azimuth) from `a` to `b` in **degrees** (0–360).
///
/// North = 0°, East = 90°, South = 180°, West = 270°.
pub fn bearing(from: &GeoCoord, to: &GeoCoord) -> f64 {
    let lat1 = from.lat.to_radians();
    let lat2 = to.lat.to_radians();
    let dlon = (to.lon - from.lon).to_radians();
    let y = dlon.sin() * lat2.cos();
    let x = lat1.cos() * lat2.sin() - lat1.sin() * lat2.cos() * dlon.cos();
    (y.atan2(x).to_degrees() + 360.0) % 360.0
}

/// Alias for [`bearing`] -- initial (forward) azimuth in degrees.
pub fn initial_bearing(from: &GeoCoord, to: &GeoCoord) -> f64 {
    bearing(from, to)
}

// ---------------------------------------------------------------------------
// Interpolation
// ---------------------------------------------------------------------------

/// Interpolate a point along the great-circle path between `from` and `to`.
///
/// `fraction` ranges from 0.0 (`from`) to 1.0 (`to`).  Values outside
/// `[0, 1]` extrapolate beyond the endpoints.
pub fn interpolate_great_circle(from: &GeoCoord, to: &GeoCoord, fraction: f64) -> GeoCoord {
    if fraction == 0.0 {
        return *from;
    }
    if fraction == 1.0 {
        return *to;
    }
    let lat1 = from.lat.to_radians();
    let lon1 = from.lon.to_radians();
    let lat2 = to.lat.to_radians();
    let lon2 = to.lon.to_radians();

    let d = {
        let dlat = lat2 - lat1;
        let dlon = lon2 - lon1;
        let a = (dlat * 0.5).sin().powi(2) + lat1.cos() * lat2.cos() * (dlon * 0.5).sin().powi(2);
        2.0 * a.sqrt().asin()
    };

    if d.abs() < 1e-15 {
        return *from;
    }

    let a = ((1.0 - fraction) * d).sin() / d.sin();
    let b = (fraction * d).sin() / d.sin();

    let x = a * lat1.cos() * lon1.cos() + b * lat2.cos() * lon2.cos();
    let y = a * lat1.cos() * lon1.sin() + b * lat2.cos() * lon2.sin();
    let z = a * lat1.sin() + b * lat2.sin();

    let lat = z.atan2((x * x + y * y).sqrt()).to_degrees();
    let lon = y.atan2(x).to_degrees();
    GeoCoord::from_lat_lon(lat, lon)
}

// ---------------------------------------------------------------------------
// Geometry trait extensions
// ---------------------------------------------------------------------------

impl Geometry {
    /// Compute the bounding box of this geometry.
    pub fn bbox(&self) -> Option<GeoBounds> {
        geometry_bbox(self)
    }

    /// Compute the centroid of this geometry.
    pub fn centroid(&self) -> Option<GeoCoord> {
        geometry_centroid(self)
    }

    /// Compute the area of this geometry in square meters.
    ///
    /// Returns 0.0 for non-areal geometries (Points, LineStrings).
    pub fn area(&self) -> f64 {
        match self {
            Geometry::Polygon(p) => polygon_area(p),
            Geometry::MultiPolygon(mp) => mp.polygons.iter().map(polygon_area).sum(),
            _ => 0.0,
        }
    }

    /// Compute the length of this geometry in meters.
    ///
    /// Returns 0.0 for non-linear geometries (Points, Polygons).
    pub fn length(&self) -> f64 {
        match self {
            Geometry::LineString(ls) => linestring_length(&ls.coords),
            Geometry::MultiLineString(mls) => {
                mls.lines.iter().map(|ls| linestring_length(&ls.coords)).sum()
            }
            _ => 0.0,
        }
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::geometry::{LineString, MultiPoint, Point};

    fn square_polygon() -> Polygon {
        Polygon {
            exterior: vec![
                GeoCoord::from_lat_lon(0.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 1.0),
                GeoCoord::from_lat_lon(1.0, 1.0),
                GeoCoord::from_lat_lon(1.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 0.0),
            ],
            interiors: vec![],
        }
    }

    // -- bbox --

    #[test]
    fn bbox_point() {
        let g = Geometry::Point(Point {
            coord: GeoCoord::from_lat_lon(10.0, 20.0),
        });
        let bb = g.bbox().unwrap();
        assert_eq!(bb.south(), 10.0);
        assert_eq!(bb.north(), 10.0);
        assert_eq!(bb.west(), 20.0);
        assert_eq!(bb.east(), 20.0);
    }

    #[test]
    fn bbox_polygon() {
        let g = Geometry::Polygon(square_polygon());
        let bb = g.bbox().unwrap();
        assert_eq!(bb.south(), 0.0);
        assert_eq!(bb.north(), 1.0);
        assert_eq!(bb.west(), 0.0);
        assert_eq!(bb.east(), 1.0);
    }

    #[test]
    fn bbox_empty_collection() {
        let g = Geometry::GeometryCollection(vec![]);
        assert!(g.bbox().is_none());
    }

    // -- centroid --

    #[test]
    fn centroid_point() {
        let g = Geometry::Point(Point {
            coord: GeoCoord::from_lat_lon(10.0, 20.0),
        });
        let c = g.centroid().unwrap();
        assert_eq!(c.lat, 10.0);
        assert_eq!(c.lon, 20.0);
    }

    #[test]
    fn centroid_square_polygon() {
        let g = Geometry::Polygon(square_polygon());
        let c = g.centroid().unwrap();
        assert!((c.lat - 0.5).abs() < 1e-10);
        assert!((c.lon - 0.5).abs() < 1e-10);
    }

    #[test]
    fn centroid_linestring() {
        let g = Geometry::LineString(LineString {
            coords: vec![
                GeoCoord::from_lat_lon(0.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 2.0),
            ],
        });
        let c = g.centroid().unwrap();
        assert!((c.lat - 0.0).abs() < 1e-10);
        assert!((c.lon - 1.0).abs() < 1e-10);
    }

    #[test]
    fn centroid_multi_point() {
        let g = Geometry::MultiPoint(MultiPoint {
            points: vec![
                Point { coord: GeoCoord::from_lat_lon(0.0, 0.0) },
                Point { coord: GeoCoord::from_lat_lon(2.0, 2.0) },
            ],
        });
        let c = g.centroid().unwrap();
        assert!((c.lat - 1.0).abs() < 1e-10);
        assert!((c.lon - 1.0).abs() < 1e-10);
    }

    // -- area --

    #[test]
    fn area_square_at_equator() {
        let area = polygon_area(&square_polygon());
        // 1° × 1° at equator ≈ 111 km × 111 km ≈ 1.234e10 m²
        assert!(area > 1e10, "area should be > 1e10, got {area}");
        assert!(area < 2e10, "area should be < 2e10, got {area}");
    }

    #[test]
    fn area_point_is_zero() {
        let g = Geometry::Point(Point {
            coord: GeoCoord::from_lat_lon(0.0, 0.0),
        });
        assert_eq!(g.area(), 0.0);
    }

    #[test]
    fn area_polygon_with_hole() {
        let poly = Polygon {
            exterior: vec![
                GeoCoord::from_lat_lon(0.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 2.0),
                GeoCoord::from_lat_lon(2.0, 2.0),
                GeoCoord::from_lat_lon(2.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 0.0),
            ],
            interiors: vec![vec![
                GeoCoord::from_lat_lon(0.5, 0.5),
                GeoCoord::from_lat_lon(0.5, 1.5),
                GeoCoord::from_lat_lon(1.5, 1.5),
                GeoCoord::from_lat_lon(1.5, 0.5),
                GeoCoord::from_lat_lon(0.5, 0.5),
            ]],
        };
        let full = polygon_area(&Polygon {
            exterior: poly.exterior.clone(),
            interiors: vec![],
        });
        let with_hole = polygon_area(&poly);
        assert!(with_hole < full, "area with hole should be less than full area");
        // Hole is 1° × 1°, exterior is 2° × 2° → ratio ≈ 3/4.
        let ratio = with_hole / full;
        assert!((ratio - 0.75).abs() < 0.05, "ratio should be ~0.75, got {ratio}");
    }

    // -- length --

    #[test]
    fn length_equator_one_degree() {
        let g = Geometry::LineString(LineString {
            coords: vec![
                GeoCoord::from_lat_lon(0.0, 0.0),
                GeoCoord::from_lat_lon(0.0, 1.0),
            ],
        });
        let len = g.length();
        // ~111 km at equator.
        assert!(len > 110_000.0);
        assert!(len < 112_000.0);
    }

    #[test]
    fn length_point_is_zero() {
        let g = Geometry::Point(Point {
            coord: GeoCoord::from_lat_lon(0.0, 0.0),
        });
        assert_eq!(g.length(), 0.0);
    }

    // -- bearing --

    #[test]
    fn bearing_due_north() {
        let from = GeoCoord::from_lat_lon(0.0, 0.0);
        let to = GeoCoord::from_lat_lon(1.0, 0.0);
        let b = bearing(&from, &to);
        assert!((b - 0.0).abs() < 0.01 || (b - 360.0).abs() < 0.01);
    }

    #[test]
    fn bearing_due_east() {
        let from = GeoCoord::from_lat_lon(0.0, 0.0);
        let to = GeoCoord::from_lat_lon(0.0, 1.0);
        let b = bearing(&from, &to);
        assert!((b - 90.0).abs() < 0.01, "expected ~90°, got {b}°");
    }

    // -- interpolation --

    #[test]
    fn interpolate_endpoints() {
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(10.0, 10.0);
        let start = interpolate_great_circle(&a, &b, 0.0);
        let end = interpolate_great_circle(&a, &b, 1.0);
        assert_eq!(start.lat, a.lat);
        assert_eq!(start.lon, a.lon);
        assert_eq!(end.lat, b.lat);
        assert_eq!(end.lon, b.lon);
    }

    #[test]
    fn interpolate_midpoint_on_equator() {
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(0.0, 10.0);
        let mid = interpolate_great_circle(&a, &b, 0.5);
        assert!((mid.lat - 0.0).abs() < 0.01);
        assert!((mid.lon - 5.0).abs() < 0.01);
    }
}