Initial commit

src/projection.rs

@@ -0,0 +1,71 @@
+/// 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()))
+    }
+
+    /// 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))
+        }
+    }
+}