gistools/proj/project/
sinu.rs

1use crate::proj::{
2    CoordinateStep, EPS10, M_VAL, N_VAL, Proj, ProjMethod, ProjectCoordinates,
3    TransformCoordinates, aasin, enfn, inv_mlfn, mlfn,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::{cell::RefCell, f64::consts::FRAC_PI_2};
7use libm::{cos, fabs, sin, sqrt};
8
9const MAX_ITER: usize = 8;
10const LOOP_TOL: f64 = 1e-7;
11
12/// Sinusoidal Variables
13#[derive(Debug, Default, Clone, PartialEq)]
14pub struct SinuData {
15    m: f64,
16    n: f64,
17    c_x: f64,
18    c_y: f64,
19    en: Vec<f64>,
20}
21
22/// for spheres, only
23fn sinu_setup(proj: &mut Proj, sinu: &mut SinuData) {
24    proj.es = 0.;
25    sinu.c_y = sqrt((sinu.m + 1.) / sinu.n);
26    sinu.c_x = sinu.c_y / (sinu.m + 1.);
27}
28
29/// # Sinusoidal (Sanson-Flamsteed)
30///
31/// **Classification**: Pseudocylindrical
32///
33/// **Available forms**: Forward and inverse, spherical and ellipsoidal
34///
35/// **Defined area**: Global
36///
37/// **Alias**: sinu
38///
39/// **Domain**: 2D
40///
41/// **Input type**: Geodetic coordinates
42///
43/// **Output type**: Projected coordinates
44///
45/// ## Projection String
46/// ```ini
47/// +proj=sinu
48/// ```
49///
50/// ## Parameters
51///
52/// All parameters are optional.
53///
54/// - `+lon_0=<value>`: Central meridian.
55/// - `+R=<value>`: Radius of the sphere or semi-major axis of the ellipsoid.
56/// - `+x_0=<value>`: False easting.
57/// - `+y_0=<value>`: False northing.
58///
59/// ## Mathematical Definition
60///
61/// MacBryde and Thomas developed generalized formulas for several of the
62/// pseudocylindricals with sinusoidal meridians. The formulas describing the Sinusoidal
63/// projection are:
64///
65/// Forward projection:
66/// $$x = C\lambda(m+cos\theta) / ( m + 1)$$
67/// $$y = C\theta$$
68///
69/// Inverse projection:
70/// $$\lambda = x \cdot \frac{m + 1}{C \cdot (m + \cos(y / C))}$$
71/// $$\theta = y / C$$
72///
73/// Where:
74/// $$C = \sqrt { (m + 1 ) / n }$$
75///
76/// ## Further Reading
77/// - [Wikipedia](https://en.wikipedia.org/wiki/Sinusoidal_projection)
78///
79/// ![Sinusoidal (Sanson-Flamsteed)](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/sinu.png?raw=true)
80#[derive(Debug, Clone, PartialEq)]
81pub struct SinusoidalProjection {
82    proj: Rc<RefCell<Proj>>,
83    store: RefCell<SinuData>,
84    method: ProjMethod,
85}
86impl ProjectCoordinates for SinusoidalProjection {
87    fn code(&self) -> i64 {
88        -1
89    }
90    fn name(&self) -> &'static str {
91        "Sinusoidal (Sanson-Flamsteed)"
92    }
93    fn names() -> &'static [&'static str] {
94        &["Sinusoidal", "Sinusoidal (Sanson-Flamsteed)", "sinu"]
95    }
96}
97impl CoordinateStep for SinusoidalProjection {
98    fn new(proj: Rc<RefCell<Proj>>) -> Self {
99        let mut store = SinuData::default();
100        let method: ProjMethod;
101        {
102            let proj = &mut proj.borrow_mut();
103
104            store.en = enfn(proj.n);
105
106            method = if proj.es != 0.0 {
107                ProjMethod::Ellipsoidal
108            } else {
109                store.n = 1.;
110                store.m = 0.;
111                sinu_setup(proj, &mut store);
112                ProjMethod::Spheroidal
113            };
114        }
115        SinusoidalProjection { proj, store: store.into(), method }
116    }
117    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
118        if self.method == ProjMethod::Ellipsoidal {
119            sinu_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
120        } else {
121            sinu_s_forward(&self.store.borrow(), p);
122        }
123    }
124    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
125        if self.method == ProjMethod::Ellipsoidal {
126            sinu_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
127        } else {
128            sinu_s_inverse(&self.store.borrow(), p);
129        }
130    }
131}
132
133/// Eckert VI Projection
134#[derive(Debug, Clone, PartialEq)]
135pub struct EckertVIProjection {
136    proj: Rc<RefCell<Proj>>,
137    store: RefCell<SinuData>,
138}
139impl ProjectCoordinates for EckertVIProjection {
140    fn code(&self) -> i64 {
141        -1
142    }
143    fn name(&self) -> &'static str {
144        "Eckert VI"
145    }
146    fn names() -> &'static [&'static str] {
147        &["Eckert VI", "eck6"]
148    }
149}
150impl CoordinateStep for EckertVIProjection {
151    fn new(proj: Rc<RefCell<Proj>>) -> Self {
152        let mut store = SinuData::default();
153        {
154            let proj = &mut proj.borrow_mut();
155
156            store.m = 1.;
157            store.n = 2.570_796_326_794_896_6;
158            sinu_setup(proj, &mut store);
159        }
160        EckertVIProjection { proj, store: store.into() }
161    }
162    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
163        sinu_s_forward(&self.store.borrow(), p);
164    }
165    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
166        sinu_s_inverse(&self.store.borrow(), p);
167    }
168}
169
170/// McBryde-Thomas Flat-Polar Sinusoidal Projection
171#[derive(Debug, Clone, PartialEq)]
172pub struct McBrydeThomasFlatPolarSinusoidalProjection {
173    proj: Rc<RefCell<Proj>>,
174    store: RefCell<SinuData>,
175}
176impl ProjectCoordinates for McBrydeThomasFlatPolarSinusoidalProjection {
177    fn code(&self) -> i64 {
178        -1
179    }
180    fn name(&self) -> &'static str {
181        "McBryde-Thomas Flat-Polar Sinusoidal"
182    }
183    fn names() -> &'static [&'static str] {
184        &["McBryde-Thomas Flat-Polar Sinusoidal", "mbtfps"]
185    }
186}
187impl CoordinateStep for McBrydeThomasFlatPolarSinusoidalProjection {
188    fn new(proj: Rc<RefCell<Proj>>) -> Self {
189        let mut store = SinuData::default();
190        {
191            let proj = &mut proj.borrow_mut();
192
193            store.m = 0.5;
194            store.n = 1.785_398_163_397_448_3;
195            sinu_setup(proj, &mut store);
196        }
197        McBrydeThomasFlatPolarSinusoidalProjection { proj, store: store.into() }
198    }
199    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
200        sinu_s_forward(&self.store.borrow(), p);
201    }
202    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
203        sinu_s_inverse(&self.store.borrow(), p);
204    }
205}
206
207/// General Sinusoidal Series Projection
208#[derive(Debug, Clone, PartialEq)]
209pub struct GeneralSinusoidalSeriesProjection {
210    proj: Rc<RefCell<Proj>>,
211    store: RefCell<SinuData>,
212}
213impl ProjectCoordinates for GeneralSinusoidalSeriesProjection {
214    fn code(&self) -> i64 {
215        -1
216    }
217    fn name(&self) -> &'static str {
218        "General Sinusoidal Series"
219    }
220    fn names() -> &'static [&'static str] {
221        &["General Sinusoidal Series", "General Sinusoidal", "gn_sinu"]
222    }
223}
224impl CoordinateStep for GeneralSinusoidalSeriesProjection {
225    fn new(proj: Rc<RefCell<Proj>>) -> Self {
226        let mut store = SinuData::default();
227        {
228            let proj = &mut proj.borrow_mut();
229
230            if let Some(n) = proj.params.get(&N_VAL) {
231                store.n = n.f64();
232            } else {
233                panic!("Missing parameter n.");
234            }
235            if let Some(m) = proj.params.get(&M_VAL) {
236                store.m = m.f64();
237            } else {
238                panic!("Missing parameter m.");
239            }
240            if store.n <= 0. {
241                panic!("Invalid value for n: it should be > 0.");
242            }
243            if store.m < 0. {
244                panic!("Invalid value for m: it should be >= 0.");
245            }
246
247            sinu_setup(proj, &mut store);
248        }
249        GeneralSinusoidalSeriesProjection { proj, store: store.into() }
250    }
251    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
252        sinu_s_forward(&self.store.borrow(), p);
253    }
254    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
255        sinu_s_inverse(&self.store.borrow(), p);
256    }
257}
258
259/// Sinusoidal Ellipsoidal forward project
260pub fn sinu_e_forward<P: TransformCoordinates>(sinu: &SinuData, proj: &Proj, p: &mut P) {
261    let s = sin(p.phi());
262    let c = cos(p.phi());
263    p.set_y(mlfn(p.phi(), s, c, &sinu.en));
264    p.set_x(p.lam() * c / sqrt(1. - proj.es * s * s));
265}
266
267/// Sinusoidal Ellipsoidal inverse project
268pub fn sinu_e_inverse<P: TransformCoordinates>(sinu: &SinuData, proj: &Proj, p: &mut P) {
269    let x = p.x();
270    let y = p.y();
271
272    p.set_phi(inv_mlfn(y, &sinu.en));
273    let mut s = fabs(p.phi());
274    if s < FRAC_PI_2 {
275        s = sin(p.phi());
276        p.set_lam(x * sqrt(1. - proj.es * s * s) / cos(p.phi()));
277    } else if (s - EPS10) < FRAC_PI_2 {
278        p.set_lam(0.);
279    } else {
280        panic!("Coordinate outside projection domain");
281    }
282}
283
284/// Sinusoidal Spheroidal forward project
285pub fn sinu_s_forward<P: TransformCoordinates>(sinu: &SinuData, p: &mut P) {
286    if sinu.m == 0.0 {
287        p.set_phi(if sinu.n != 1. { aasin(sinu.n * sin(p.phi())) } else { p.phi() });
288    } else {
289        let k = sinu.n * sin(p.phi());
290        let mut i = MAX_ITER;
291        while i > 0 {
292            i -= 1;
293            let v = (sinu.m * p.phi() + sin(p.phi()) - k) / (sinu.m + cos(p.phi()));
294            p.set_phi(p.phi() - v);
295            if fabs(v) < LOOP_TOL {
296                break;
297            }
298        }
299        if i != 0 {
300            panic!("Coordinate outside projection domain");
301        }
302    }
303    p.set_x(sinu.c_x * p.lam() * (sinu.m + cos(p.phi())));
304    p.set_y(sinu.c_y * p.phi());
305}
306
307/// Sinusoidal Spheroidal inverse project
308pub fn sinu_s_inverse<P: TransformCoordinates>(sinu: &SinuData, p: &mut P) {
309    let mut y = p.y();
310
311    y /= sinu.c_y;
312    p.set_phi(if sinu.m != 0.0 {
313        aasin((sinu.m * y + sin(y)) / sinu.n)
314    } else if sinu.n != 1. {
315        aasin(sin(y) / sinu.n)
316    } else {
317        y
318    });
319    p.set_lam(p.x() / (sinu.c_x * (sinu.m + cos(y))));
320}
gistools/proj/project/sinu.rs

gistools/proj/project/
sinu.rs
