gistools/proj/project/
oea.rs

1use crate::proj::{
2    CoordinateStep, M_VAL, N_VAL, Proj, ProjValue, ProjectCoordinates, THETA, TransformCoordinates,
3    aacos, aasin, aatan2,
4};
5use alloc::rc::Rc;
6use core::cell::RefCell;
7use libm::{cos, hypot, sin};
8
9/// Oblated Equal Area Variables
10#[derive(Debug, Default, Clone, PartialEq)]
11pub struct OeaData {
12    theta: f64,
13    m: f64,
14    n: f64,
15    two_r_m: f64,
16    two_r_n: f64,
17    rm: f64,
18    rn: f64,
19    hm: f64,
20    hn: f64,
21    cp0: f64,
22    sp0: f64,
23}
24
25/// Oblated Equal Area Projection
26#[derive(Debug, Clone, PartialEq)]
27pub struct OblatedEqualAreaProjection {
28    proj: Rc<RefCell<Proj>>,
29    store: RefCell<OeaData>,
30}
31impl ProjectCoordinates for OblatedEqualAreaProjection {
32    fn code(&self) -> i64 {
33        -1
34    }
35    fn name(&self) -> &'static str {
36        "Oblated Equal Area"
37    }
38    fn names() -> &'static [&'static str] {
39        &["Oblated Equal Area", "oea"]
40    }
41}
42impl CoordinateStep for OblatedEqualAreaProjection {
43    fn new(proj: Rc<RefCell<Proj>>) -> Self {
44        let mut store = OeaData::default();
45        {
46            let proj = &mut proj.borrow_mut();
47
48            store.n = proj.params.get(&N_VAL).unwrap_or(&ProjValue::default()).f64();
49            if store.n <= 0. {
50                panic!("Invalid value for n: it should be > 0");
51            }
52            store.m = proj.params.get(&M_VAL).unwrap_or(&ProjValue::default()).f64();
53            if store.m <= 0. {
54                panic!("Invalid value for m: it should be > 0");
55            }
56            store.theta =
57                proj.params.get(&THETA).unwrap_or(&ProjValue::default()).f64().to_radians();
58
59            store.sp0 = sin(proj.phi0);
60            store.cp0 = cos(proj.phi0);
61            store.rn = 1. / store.n;
62            store.rm = 1. / store.m;
63            store.two_r_n = 2. * store.rn;
64            store.two_r_m = 2. * store.rm;
65            store.hm = 0.5 * store.m;
66            store.hn = 0.5 * store.n;
67            proj.es = 0.;
68        }
69        OblatedEqualAreaProjection { proj, store: store.into() }
70    }
71    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
72        oea_s_forward(&self.store.borrow(), p);
73    }
74    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
75        oea_s_inverse(&self.store.borrow(), p);
76    }
77}
78
79/// Oblated Equal Area Spheroidal forward project
80pub fn oea_s_forward<P: TransformCoordinates>(oae: &OeaData, p: &mut P) {
81    let cp = cos(p.phi());
82    let sp = sin(p.phi());
83    let cl = cos(p.lam());
84    let az = aatan2(cp * sin(p.lam()), oae.cp0 * sp - oae.sp0 * cp * cl) + oae.theta;
85    let shz = sin(0.5 * aacos(oae.sp0 * sp + oae.cp0 * cp * cl));
86    let m = aasin(shz * sin(az));
87    let n = aasin(shz * cos(az) * cos(m) / cos(m * oae.two_r_m));
88    p.set_y(oae.n * sin(n * oae.two_r_n));
89    p.set_x(oae.m * sin(m * oae.two_r_m) * cos(n) / cos(n * oae.two_r_n));
90}
91
92/// Oblated Equal Area Spheroidal inverse project
93pub fn oea_s_inverse<P: TransformCoordinates>(oae: &OeaData, p: &mut P) {
94    let n = oae.hn * aasin(p.y() * oae.rn);
95    let m = oae.hm * aasin(p.x() * oae.rm * cos(n * oae.two_r_n) / cos(n));
96    let xp = 2. * sin(m);
97    let yp = 2. * sin(n) * cos(m * oae.two_r_m) / cos(m);
98    let az = aatan2(xp, yp) - oae.theta;
99    let c_az = cos(az);
100    let z = 2. * aasin(0.5 * hypot(xp, yp));
101    let sz = sin(z);
102    let cz = cos(z);
103    p.set_phi(aasin(oae.sp0 * cz + oae.cp0 * sz * c_az));
104    p.set_lam(aatan2(sz * sin(az), oae.cp0 * cz - oae.sp0 * sz * c_az));
105}
gistools/proj/project/oea.rs

gistools/proj/project/
oea.rs
