gistools/proj/project/
moll.rs

1use crate::proj::{CoordinateStep, Proj, ProjectCoordinates, TransformCoordinates, aasin};
2use alloc::rc::Rc;
3use core::{
4    cell::RefCell,
5    f64::consts::{FRAC_PI_2, PI, TAU},
6};
7use libm::{cos, fabs, sin, sqrt};
8
9const MAX_ITER: usize = 30;
10const LOOP_TOL: f64 = 1e-7;
11
12/// Mollweide Variables
13#[derive(Debug, Default, Clone, PartialEq)]
14pub struct MollData {
15    c_x: f64,
16    c_y: f64,
17    c_p: f64,
18}
19
20fn setup(proj: &mut Proj, p: f64) -> MollData {
21    let mut store = MollData::default();
22    let p2 = p + p;
23
24    proj.es = 0.;
25    let sp = sin(p);
26    let r = sqrt(TAU * sp / (p2 + sin(p2)));
27
28    store.c_x = 2. * r / PI;
29    store.c_y = r / sp;
30    store.c_p = p2 + sin(p2);
31
32    store
33}
34
35/// # Mollweide
36///
37/// **Classification**: Pseudocylindrical
38///
39/// **Available forms**: Forward and inverse, spherical projection
40///
41/// **Defined area**: Global
42///
43/// **Alias**: moll
44///
45/// **Domain**: 2D
46///
47/// **Input type**: Geodetic coordinates
48///
49/// **Output type**: Projected coordinates
50///
51/// ## Projection String
52/// ```ini
53/// +proj=moll
54/// ```
55///
56/// ## Required Parameters
57/// - None
58///
59/// ## Optional Parameters
60/// - `+lon_0`: Longitude of projection center. Defaults to `0`.
61/// - `+R`: Radius of the sphere.
62/// - `+x_0`: False easting. Defaults to `0`.
63/// - `+y_0`: False northing. Defaults to `0`.
64///
65/// ## Further reading
66/// - [Wikipedia on Mollweide Projection](https://en.wikipedia.org/wiki/Mollweide_projection)
67///
68/// ![Mollweide](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/moll.png?raw=true)
69#[derive(Debug, Clone, PartialEq)]
70pub struct MollweideProjection {
71    proj: Rc<RefCell<Proj>>,
72    store: RefCell<MollData>,
73}
74impl ProjectCoordinates for MollweideProjection {
75    fn code(&self) -> i64 {
76        -1
77    }
78    fn name(&self) -> &'static str {
79        "Mollweide"
80    }
81    fn names() -> &'static [&'static str] {
82        &["Mollweide", "moll"]
83    }
84}
85impl CoordinateStep for MollweideProjection {
86    fn new(proj: Rc<RefCell<Proj>>) -> Self {
87        let store = setup(&mut proj.borrow_mut(), FRAC_PI_2);
88        MollweideProjection { proj, store: store.into() }
89    }
90    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
91        moll_s_forward(&self.store.borrow(), p);
92    }
93    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
94        moll_s_inverse(&self.store.borrow(), p);
95    }
96}
97
98/// Wagner IV Projection
99#[derive(Debug, Clone, PartialEq)]
100pub struct WagnerIVProjection {
101    proj: Rc<RefCell<Proj>>,
102    store: RefCell<MollData>,
103}
104impl ProjectCoordinates for WagnerIVProjection {
105    fn code(&self) -> i64 {
106        -1
107    }
108    fn name(&self) -> &'static str {
109        "Wagner IV"
110    }
111    fn names() -> &'static [&'static str] {
112        &["Wagner IV", "wag4"]
113    }
114}
115impl CoordinateStep for WagnerIVProjection {
116    fn new(proj: Rc<RefCell<Proj>>) -> Self {
117        let store = setup(&mut proj.borrow_mut(), PI / 3.);
118        WagnerIVProjection { proj, store: store.into() }
119    }
120    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
121        moll_s_forward(&self.store.borrow(), p);
122    }
123    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
124        moll_s_inverse(&self.store.borrow(), p);
125    }
126}
127
128/// Wagner V Projection
129#[derive(Debug, Clone, PartialEq)]
130pub struct WagnerVProjection {
131    proj: Rc<RefCell<Proj>>,
132    store: RefCell<MollData>,
133}
134impl ProjectCoordinates for WagnerVProjection {
135    fn code(&self) -> i64 {
136        -1
137    }
138    fn name(&self) -> &'static str {
139        "Wagner V"
140    }
141    fn names() -> &'static [&'static str] {
142        &["Wagner V", "wag5"]
143    }
144}
145impl CoordinateStep for WagnerVProjection {
146    fn new(proj: Rc<RefCell<Proj>>) -> Self {
147        proj.borrow_mut().es = 0.0;
148        let store = MollData { c_x: 0.90977, c_y: 1.65014, c_p: 3.00896 };
149        WagnerVProjection { proj, store: store.into() }
150    }
151    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
152        moll_s_forward(&self.store.borrow(), p);
153    }
154    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
155        moll_s_inverse(&self.store.borrow(), p);
156    }
157}
158
159/// Mollweide Spheroidal inverse project
160pub fn moll_s_forward<P: TransformCoordinates>(moll: &MollData, p: &mut P) {
161    let k = moll.c_p * sin(p.phi());
162    let mut i = MAX_ITER;
163    while i > 0 {
164        let v = (p.phi() + sin(p.phi()) - k) / (1. + cos(p.phi()));
165        p.set_phi(p.phi() - v);
166        if fabs(v) < LOOP_TOL {
167            break;
168        }
169        i -= 1;
170    }
171    if i != 0 {
172        p.set_phi(if p.phi() < 0. { -FRAC_PI_2 } else { FRAC_PI_2 });
173    } else {
174        p.set_phi(p.phi() * 0.5);
175    }
176    p.set_x(moll.c_x * p.lam() * cos(p.phi()));
177    p.set_y(moll.c_y * sin(p.phi()));
178}
179
180/// Mollweide Spheroidal inverse project
181pub fn moll_s_inverse<P: TransformCoordinates>(moll: &MollData, p: &mut P) {
182    let mut phi = aasin(p.y() / moll.c_y);
183    let mut lam = p.x() / (moll.c_x * cos(phi));
184    if fabs(lam) < PI {
185        phi += phi;
186        phi = aasin((phi + sin(phi)) / moll.c_p);
187    } else {
188        phi = f64::MAX;
189        lam = f64::MAX;
190    }
191    p.set_phi(phi);
192    p.set_lam(lam);
193}
gistools/proj/project/moll.rs

gistools/proj/project/
moll.rs
