rustial_engine/visualization/

geo_grid.rs

//! Georeferenced regular grid descriptor.

use rustial_math::GeoCoord;

/// A georeferenced regular grid anchored to a geographic origin.
///
/// Cell `(0, 0)` is the **north-west** corner. Row index increases
/// southward, column index increases eastward. All dimensions are in
/// meters.
#[derive(Debug, Clone)]
pub struct GeoGrid {
    /// Geographic origin (north-west corner of cell `(0, 0)`).
    pub origin: GeoCoord,
    /// Number of rows (south-ward).
    pub rows: usize,
    /// Number of columns (east-ward).
    pub cols: usize,
    /// Cell width in meters (east-west extent per cell).
    pub cell_width: f64,
    /// Cell height in meters (north-south extent per cell).
    pub cell_height: f64,
    /// Grid rotation in radians (clockwise from north). Default `0.0`.
    pub rotation: f64,
    /// Altitude mode for the grid surface.
    pub altitude_mode: crate::models::AltitudeMode,
}

/// Approximate meters per degree of latitude (WGS-84 mean).
const METERS_PER_DEG_LAT: f64 = 111_320.0;

impl GeoGrid {
    /// Create a new grid with the given dimensions and cell size.
    ///
    /// `origin` is the north-west corner. `cell_width` and `cell_height`
    /// are in meters.
    pub fn new(
        origin: GeoCoord,
        rows: usize,
        cols: usize,
        cell_width: f64,
        cell_height: f64,
    ) -> Self {
        Self {
            origin,
            rows,
            cols,
            cell_width,
            cell_height,
            rotation: 0.0,
            altitude_mode: crate::models::AltitudeMode::ClampToGround,
        }
    }

    /// Total number of cells.
    #[inline]
    pub fn cell_count(&self) -> usize {
        self.rows * self.cols
    }

    /// Geographic coordinate of the centre of cell `(row, col)`.
    ///
    /// Returns `None` if `row >= rows` or `col >= cols`.
    pub fn cell_center(&self, row: usize, col: usize) -> Option<GeoCoord> {
        if row >= self.rows || col >= self.cols {
            return None;
        }
        // Offset from origin in meters (unrotated).
        let dx = (col as f64 + 0.5) * self.cell_width;
        let dy = (row as f64 + 0.5) * self.cell_height;

        // Apply rotation.
        let (sin_r, cos_r) = self.rotation.sin_cos();
        let rx = dx * cos_r - dy * sin_r;
        let ry = dx * sin_r + dy * cos_r;

        Some(offset_geo(&self.origin, rx, ry))
    }

    /// Geographic bounding box of the grid `(north-west, south-east)`.
    pub fn geo_bounds(&self) -> (GeoCoord, GeoCoord) {
        let total_dx = self.cols as f64 * self.cell_width;
        let total_dy = self.rows as f64 * self.cell_height;

        if self.rotation.abs() < 1e-12 {
            let se = offset_geo(&self.origin, total_dx, total_dy);
            return (self.origin, se);
        }

        // Sample corners and compute bounding box.
        let corners = [
            (0.0, 0.0),
            (total_dx, 0.0),
            (0.0, total_dy),
            (total_dx, total_dy),
        ];
        let (sin_r, cos_r) = self.rotation.sin_cos();
        let mut min_lat = f64::MAX;
        let mut max_lat = f64::MIN;
        let mut min_lon = f64::MAX;
        let mut max_lon = f64::MIN;
        for &(dx, dy) in &corners {
            let rx = dx * cos_r - dy * sin_r;
            let ry = dx * sin_r + dy * cos_r;
            let c = offset_geo(&self.origin, rx, ry);
            min_lat = min_lat.min(c.lat);
            max_lat = max_lat.max(c.lat);
            min_lon = min_lon.min(c.lon);
            max_lon = max_lon.max(c.lon);
        }
        (
            GeoCoord::from_lat_lon(max_lat, min_lon),
            GeoCoord::from_lat_lon(min_lat, max_lon),
        )
    }

    /// Find the grid cell `(row, col)` containing a geographic coordinate.
    ///
    /// Returns `None` if the coordinate is outside the grid.
    pub fn cell_at_geo(&self, coord: &GeoCoord) -> Option<(usize, usize)> {
        // Offset from origin in meters.
        let (dx, dy) = geo_offset(&self.origin, coord);

        // Undo rotation.
        let (sin_r, cos_r) = self.rotation.sin_cos();
        let ux = dx * cos_r + dy * sin_r;
        let uy = -dx * sin_r + dy * cos_r;

        if ux < 0.0 || uy < 0.0 {
            return None;
        }

        let col = (ux / self.cell_width) as usize;
        let row = (uy / self.cell_height) as usize;

        if row < self.rows && col < self.cols {
            Some((row, col))
        } else {
            None
        }
    }
}

/// Offset a `GeoCoord` by `dx` meters east and `dy` meters south.
fn offset_geo(origin: &GeoCoord, dx_meters: f64, dy_meters: f64) -> GeoCoord {
    let lat = origin.lat - dy_meters / METERS_PER_DEG_LAT;
    let cos_lat = origin.lat.to_radians().cos().max(1e-10);
    let lon = origin.lon + dx_meters / (METERS_PER_DEG_LAT * cos_lat);
    GeoCoord::from_lat_lon(lat, lon)
}

/// Compute the offset in meters from `origin` to `coord` (`dx` east, `dy` south).
fn geo_offset(origin: &GeoCoord, coord: &GeoCoord) -> (f64, f64) {
    let cos_lat = origin.lat.to_radians().cos().max(1e-10);
    let dx = (coord.lon - origin.lon) * METERS_PER_DEG_LAT * cos_lat;
    let dy = (origin.lat - coord.lat) * METERS_PER_DEG_LAT;
    (dx, dy)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn cell_center_round_trips_with_cell_at_geo() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(51.1, 17.0), 10, 10, 100.0, 100.0);
        for row in 0..grid.rows {
            for col in 0..grid.cols {
                let center = grid.cell_center(row, col).unwrap();
                let (r, c) = grid.cell_at_geo(&center).unwrap();
                assert_eq!((r, c), (row, col), "round-trip failed for ({row}, {col})");
            }
        }
    }

    #[test]
    fn cell_center_out_of_bounds() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(0.0, 0.0), 5, 5, 50.0, 50.0);
        assert!(grid.cell_center(5, 0).is_none());
        assert!(grid.cell_center(0, 5).is_none());
    }

    #[test]
    fn cell_at_geo_outside_grid() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(51.1, 17.0), 5, 5, 100.0, 100.0);
        // Far south of the grid
        assert!(grid.cell_at_geo(&GeoCoord::from_lat_lon(40.0, 17.0)).is_none());
        // West of the grid
        assert!(grid.cell_at_geo(&GeoCoord::from_lat_lon(51.1, 16.0)).is_none());
    }

    #[test]
    fn geo_bounds_no_rotation() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(51.1, 17.0), 10, 10, 100.0, 100.0);
        let (nw, se) = grid.geo_bounds();
        assert!((nw.lat - 51.1).abs() < 1e-6);
        assert!((nw.lon - 17.0).abs() < 1e-6);
        assert!(se.lat < nw.lat);
        assert!(se.lon > nw.lon);
    }

    #[test]
    fn cell_count() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(0.0, 0.0), 3, 7, 10.0, 10.0);
        assert_eq!(grid.cell_count(), 21);
    }

    #[test]
    fn geo_bounds_at_high_latitude() {
        let grid = GeoGrid::new(GeoCoord::from_lat_lon(70.0, 25.0), 5, 5, 200.0, 200.0);
        let (nw, se) = grid.geo_bounds();
        assert!(se.lat < nw.lat);
        assert!(se.lon > nw.lon);
    }
}