dwd_cloud_cover_rs/src/projection.rs

/// Orthographic map projection.
///
/// Projects a geographic point (lat, lon) onto a 2-D plane as seen from directly
/// above the centre point (center_lat, center_lon).
///
/// Returns (x, y) normalised so that an angular distance of 1° from the centre
/// equals 1.0 in output units.  Returns `None` when the point is on the far side
/// of the globe (not visible from the projection centre).

pub struct OrthoProjection {
    /// Centre of projection in radians
    phi0: f64,
    lam0: f64,
    /// Half-extents of the visible window in "projected units" (degrees of arc)
    pub half_width: f64,
    pub half_height: f64,
}

impl OrthoProjection {
    /// `center_lat`, `center_lon` in degrees.
    /// `half_width`, `half_height` are the angular half-extents of the view window
    /// in degrees (e.g. 0.8 and 0.4 to match the original Python script).
    pub fn new(center_lat: f64, center_lon: f64, half_width: f64, half_height: f64) -> Self {
        Self {
            phi0: center_lat.to_radians(),
            lam0: center_lon.to_radians(),
            half_width,
            half_height,
        }
    }

    /// Project `(lat, lon)` in degrees.
    ///
    /// Returns `Some((x, y))` in degrees of arc from the centre (same units as
    /// `half_width`/`half_height`), or `None` if the point is not visible.
    pub fn project(&self, lat: f64, lon: f64) -> Option<(f64, f64)> {
        let phi = lat.to_radians();
        let lam = lon.to_radians();
        let dlam = lam - self.lam0;

        let cos_c = self.phi0.sin() * phi.sin() + self.phi0.cos() * phi.cos() * dlam.cos();
        if cos_c <= 0.0 {
            return None; // back of globe
        }

        // x and y are in radians; convert to degrees so they match the half_width unit
        let x = phi.cos() * dlam.sin();
        let y = self.phi0.cos() * phi.sin() - self.phi0.sin() * phi.cos() * dlam.cos();

        Some((x.to_degrees(), y.to_degrees()))
    }

    /// Inverse of `project`: convert projected (x, y) in degrees-of-arc back to
    /// geographic (lat, lon) in degrees.
    pub fn unproject(&self, x: f64, y: f64) -> (f64, f64) {
        let x_rad = x.to_radians();
        let y_rad = y.to_radians();
        let rho = (x_rad * x_rad + y_rad * y_rad).sqrt();

        if rho < 1e-14 {
            // At the projection centre
            return (self.phi0.to_degrees(), self.lam0.to_degrees());
        }

        let c = rho.asin();
        let sin_c = c.sin();
        let cos_c = c.cos();

        let lat = (cos_c * self.phi0.sin() + y_rad * sin_c * self.phi0.cos() / rho).asin();
        let lon = self.lam0
            + (x_rad * sin_c).atan2(rho * self.phi0.cos() * cos_c - y_rad * self.phi0.sin() * sin_c);

        (lat.to_degrees(), lon.to_degrees())
    }

    /// Inverse of `to_pixel`: convert pixel (col, row) back to projected (x, y)
    /// in degrees-of-arc.
    pub fn from_pixel(&self, col: f64, row: f64, width: u32, height: u32) -> (f64, f64) {
        let x = col / width as f64 * (2.0 * self.half_width) - self.half_width;
        let y = self.half_height - row / height as f64 * (2.0 * self.half_height);
        (x, y)
    }

    /// Convenience: convert pixel (col, row) directly to geographic (lat, lon).
    pub fn pixel_to_geo(&self, col: f64, row: f64, width: u32, height: u32) -> (f64, f64) {
        let (x, y) = self.from_pixel(col, row, width, height);
        self.unproject(x, y)
    }

    /// Convert projected (x, y) in degrees-of-arc to pixel (col, row).
    ///
    /// Returns `None` if the point falls outside the image boundary.
    pub fn to_pixel(&self, x: f64, y: f64, width: u32, height: u32) -> Option<(i32, i32)> {
        // Map [-half_width, +half_width] → [0, width]
        // Map [+half_height, -half_height] → [0, height]  (y flipped: north = up = row 0)
        let col = (x + self.half_width) / (2.0 * self.half_width) * width as f64;
        let row = (-y + self.half_height) / (2.0 * self.half_height) * height as f64;

        let col_i = col as i32;
        let row_i = row as i32;

        if col_i < 0 || col_i >= width as i32 || row_i < 0 || row_i >= height as i32 {
            None
        } else {
            Some((col_i, row_i))
        }
    }
}