gistools/proj/project/
stere.rs

1use crate::proj::{
2    CoordinateStep, EPS10, LATITUDE_STD_PARALLEL, POLAR_STEREOGRAPHIC_VARIANT_A,
3    POLAR_STEREOGRAPHIC_VARIANT_B, POLAR_STEREOGRAPHIC_VARIANT_C, Proj, ProjMethod, ProjMode,
4    ProjectCoordinates, SOUTH, TransformCoordinates, tsfn,
5};
6use alloc::rc::Rc;
7use core::{
8    cell::RefCell,
9    f64::consts::{FRAC_PI_2, FRAC_PI_4},
10};
11use libm::{asin, atan, atan2, cos, fabs, hypot, pow, sin, sqrt, tan};
12
13/// Stereographic Variables
14#[derive(Debug, Default, Clone, PartialEq)]
15pub struct StereData {
16    phits: f64,
17    sin_x1: f64,
18    cos_x1: f64,
19    akm1: f64,
20    mode: ProjMode,
21}
22
23const TOL: f64 = 1e-8;
24const NITER: usize = 8;
25const CONV: f64 = 1e-10;
26
27fn ssfn_(phit: f64, mut sinphi: f64, eccen: f64) -> f64 {
28    sinphi *= eccen;
29    tan(0.5 * (FRAC_PI_2 + phit)) * pow((1. - sinphi) / (1. + sinphi), 0.5 * eccen)
30}
31
32/// general stereographic initialization
33fn stere_setup(proj: &mut Proj, store: &mut StereData) -> ProjMethod {
34    let t = fabs(proj.phi0);
35    if fabs(t - FRAC_PI_2) < EPS10 {
36        store.mode = if proj.phi0 < 0. { ProjMode::SPole } else { ProjMode::NPole };
37    } else {
38        store.mode = if t > EPS10 { ProjMode::Obliq } else { ProjMode::Equit };
39    }
40    store.phits = fabs(store.phits);
41
42    if proj.es != 0.0 {
43        match store.mode {
44            ProjMode::NPole | ProjMode::SPole => {
45                if fabs(store.phits - FRAC_PI_2) < EPS10 {
46                    store.akm1 = 2. * proj.k0
47                        / sqrt(pow(1. + proj.e, 1. + proj.e) * pow(1. - proj.e, 1. - proj.e));
48                } else {
49                    let mut t = sin(store.phits);
50                    store.akm1 = cos(store.phits) / tsfn(store.phits, t, proj.e);
51                    t *= proj.e;
52                    store.akm1 /= sqrt(1. - t * t);
53                }
54            }
55            ProjMode::Equit | ProjMode::Obliq => {
56                let mut t = sin(proj.phi0);
57                let x = 2. * atan(ssfn_(proj.phi0, t, proj.e)) - FRAC_PI_2;
58                t *= proj.e;
59                store.akm1 = 2. * proj.k0 * cos(proj.phi0) / sqrt(1. - t * t);
60                store.sin_x1 = sin(x);
61                store.cos_x1 = cos(x);
62            }
63        }
64        ProjMethod::Ellipsoidal
65    } else {
66        match store.mode {
67            ProjMode::Obliq => {
68                store.sin_x1 = sin(proj.phi0);
69                store.cos_x1 = cos(proj.phi0);
70                store.akm1 = 2. * proj.k0;
71            }
72            ProjMode::Equit => {
73                store.akm1 = 2. * proj.k0;
74            }
75            ProjMode::SPole | ProjMode::NPole => {
76                store.akm1 = if fabs(store.phits - FRAC_PI_2) >= EPS10 {
77                    cos(store.phits) / tan(FRAC_PI_4 - 0.5 * store.phits)
78                } else {
79                    2. * proj.k0
80                };
81            }
82        }
83
84        ProjMethod::Spheroidal
85    }
86}
87
88/// # Stereographic
89///
90/// **Classification**: Azimuthal
91///
92/// **Available forms**: Forward and inverse, spherical and ellipsoidal
93///
94/// **Defined area**: Global
95///
96/// **Alias**: stere
97///
98/// **Domain**: 2D
99///
100/// **Input type**: Geodetic coordinates
101///
102/// **Output type**: Projected coordinates
103///
104/// ## Projection String
105/// ```ini
106/// +proj=stere +lat_0=90 +latTs=75
107/// ```
108///
109/// Note:
110/// This projection method gives different results than the :ref:`sterea`
111/// method in the non-polar cases (i.e. the oblique and equatorial case). The later
112/// projection method is the one referenced by EPSG as "Oblique Stereographic".
113///
114/// ## Required Parameters
115/// - None
116///
117/// ## Optional Parameters
118/// - `+lat_0=<value>`: Latitude of origin.
119/// - `+latTs=<value>`: Latitude where scale is not distorted.
120/// - `+k_0=<value>`: Scale factor.
121/// - `+lon_0=<value>`: Central meridian.
122/// - `+ellps=<value>`: Ellipsoid used.
123/// - `+R=<value>`: Radius of the projection sphere.
124/// - `+x_0=<value>`: False easting.
125/// - `+y_0=<value>`: False northing.
126///
127/// ![Stereographic](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/stere.png?raw=true)
128#[derive(Debug, Clone, PartialEq)]
129pub struct StereographicProjection {
130    proj: Rc<RefCell<Proj>>,
131    store: RefCell<StereData>,
132    method: ProjMethod,
133}
134impl ProjectCoordinates for StereographicProjection {
135    fn code(&self) -> i64 {
136        -1
137    }
138    fn name(&self) -> &'static str {
139        "Stereographic"
140    }
141    fn names() -> &'static [&'static str] {
142        &[
143            "Stereographic",
144            "Polar_Stereographic",
145            "StereographicSouthPole",
146            "Stereographic_South_Pole",
147            "Stereographic South Pole",
148            "Polar Stereographic (variant B)",
149            "stere",
150        ]
151    }
152}
153impl CoordinateStep for StereographicProjection {
154    fn new(proj: Rc<RefCell<Proj>>) -> Self {
155        let mut store = StereData::default();
156        let method: ProjMethod;
157        {
158            let proj = &mut proj.borrow_mut();
159            store.phits = if let Some(lat_ts) = proj.params.get(&LATITUDE_STD_PARALLEL) {
160                lat_ts.f64()
161            } else {
162                FRAC_PI_2
163            };
164
165            method = stere_setup(proj, &mut store);
166        }
167        StereographicProjection { proj, store: store.into(), method }
168    }
169    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
170        if self.method == ProjMethod::Spheroidal {
171            stere_s_forward(&mut self.store.borrow_mut(), p);
172        } else {
173            stere_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
174        }
175    }
176    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
177        if self.method == ProjMethod::Spheroidal {
178            stere_s_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
179        } else {
180            stere_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
181        }
182    }
183}
184
185/// Polar Stereographic Variant A Projection
186pub type PolarStereographicVariantAProjection =
187    UniversalPolarStereographicProjection<POLAR_STEREOGRAPHIC_VARIANT_A>;
188/// Polar Stereographic Variant B Projection
189pub type PolarStereographicVariantBProjection =
190    UniversalPolarStereographicProjection<POLAR_STEREOGRAPHIC_VARIANT_B>;
191/// Polar Stereographic Variant C Projection
192pub type PolarStereographicVariantCProjection =
193    UniversalPolarStereographicProjection<POLAR_STEREOGRAPHIC_VARIANT_C>;
194
195/// Stereographic Projection
196#[derive(Debug, Clone, PartialEq)]
197pub struct UniversalPolarStereographicProjection<const C: i64> {
198    proj: Rc<RefCell<Proj>>,
199    store: RefCell<StereData>,
200    method: ProjMethod,
201}
202impl<const C: i64> ProjectCoordinates for UniversalPolarStereographicProjection<C> {
203    fn code(&self) -> i64 {
204        C
205    }
206    fn name(&self) -> &'static str {
207        "Universal Polar Stereographic"
208    }
209    fn names() -> &'static [&'static str] {
210        &[
211            "Polar Stereographic",
212            "Universal Polar Stereographic",
213            "Polar Stereographic (variant A)",
214            "Polar Stereographic (variant B)",
215            "Polar Stereographic (variant C)",
216        ]
217    }
218}
219impl<const C: i64> CoordinateStep for UniversalPolarStereographicProjection<C> {
220    fn new(proj: Rc<RefCell<Proj>>) -> Self {
221        let mut store = StereData::default();
222        let method: ProjMethod;
223        {
224            let proj = &mut proj.borrow_mut();
225            proj.phi0 = if proj.params.contains_key(&SOUTH) { -FRAC_PI_2 } else { FRAC_PI_2 };
226            if proj.es == 0.0 {
227                panic!("Invalid value for es: only ellipsoidal formulation supported");
228            }
229            proj.k0 = 0.994;
230            proj.x0 = 2000000.;
231            proj.y0 = 2000000.;
232            store.phits = FRAC_PI_2;
233            proj.lam0 = 0.;
234
235            method = stere_setup(proj, &mut store);
236        }
237        UniversalPolarStereographicProjection { proj, store: store.into(), method }
238    }
239    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
240        if self.method == ProjMethod::Spheroidal {
241            stere_s_forward(&mut self.store.borrow_mut(), p);
242        } else {
243            stere_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
244        }
245    }
246    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
247        if self.method == ProjMethod::Spheroidal {
248            stere_s_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
249        } else {
250            stere_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
251        }
252    }
253}
254
255/// Stereographic Ellipsoidal forward project
256pub fn stere_e_forward<P: TransformCoordinates>(stere: &mut StereData, proj: &Proj, p: &mut P) {
257    let mut x;
258    let y;
259    let mut coslam = cos(p.lam());
260    let sinlam = sin(p.lam());
261    let mut sinphi = sin(p.phi());
262    let mut sin_x = 0.;
263    let mut cos_x = 0.;
264    if stere.mode == ProjMode::Obliq || stere.mode == ProjMode::Equit {
265        let x = 2. * atan(ssfn_(p.phi(), sinphi, proj.e)) - FRAC_PI_2;
266        sin_x = sin(x);
267        cos_x = cos(x);
268    }
269
270    match stere.mode {
271        ProjMode::Obliq => {
272            let denom = stere.cos_x1 * (1. + stere.sin_x1 * sin_x + stere.cos_x1 * cos_x * coslam);
273            if denom == 0. {
274                panic!("Coordinate outside projection domain");
275            }
276            let a = stere.akm1 / denom;
277            y = a * (stere.cos_x1 * sin_x - stere.sin_x1 * cos_x * coslam);
278            x = a * cos_x;
279        }
280        ProjMode::Equit => {
281            // avoid zero division
282            let mut a = 0.;
283            if 1. + cos_x * coslam == 0.0 {
284                y = f64::MAX;
285            } else {
286                a = stere.akm1 / (1. + cos_x * coslam);
287                y = a * sin_x;
288            }
289            x = a * cos_x;
290        }
291        ProjMode::SPole => {
292            p.set_phi(-p.phi());
293            coslam = -coslam;
294            sinphi = -sinphi;
295            if fabs(p.phi() - FRAC_PI_2) < 1e-15 {
296                x = 0.;
297            } else {
298                x = stere.akm1 * tsfn(p.phi(), sinphi, proj.e);
299            }
300            y = -x * coslam;
301        }
302        ProjMode::NPole => {
303            if fabs(p.phi() - FRAC_PI_2) < 1e-15 {
304                x = 0.;
305            } else {
306                x = stere.akm1 * tsfn(p.phi(), sinphi, proj.e);
307            }
308            y = -x * coslam;
309        }
310    }
311
312    x *= sinlam;
313    p.set_x(x);
314    p.set_y(y);
315}
316
317/// Stereographic Spheroidal forward project
318pub fn stere_s_forward<P: TransformCoordinates>(stere: &mut StereData, p: &mut P) {
319    let sinphi = sin(p.phi());
320    let cosphi = cos(p.phi());
321    let mut coslam = cos(p.lam());
322    let sinlam = sin(p.lam());
323    let x;
324    let mut y;
325
326    match stere.mode {
327        ProjMode::Equit => {
328            y = 1. + cosphi * coslam;
329            if y <= EPS10 {
330                panic!("Coordinate outside projection domain");
331            }
332            y = stere.akm1 / y;
333            x = y * cosphi * sinlam;
334            y *= if stere.mode == ProjMode::Equit {
335                sinphi
336            } else {
337                stere.cos_x1 * sinphi - stere.sin_x1 * cosphi * coslam
338            };
339        }
340        ProjMode::Obliq => {
341            y = 1. + stere.sin_x1 * sinphi + stere.cos_x1 * cosphi * coslam;
342            if y <= EPS10 {
343                panic!("Coordinate outside projection domain");
344            }
345            y = stere.akm1 / y;
346            x = y * cosphi * sinlam;
347            y *= if stere.mode == ProjMode::Equit {
348                sinphi
349            } else {
350                stere.cos_x1 * sinphi - stere.sin_x1 * cosphi * coslam
351            };
352        }
353        ProjMode::NPole => {
354            coslam = -coslam;
355            p.set_phi(-p.phi());
356            if fabs(p.phi() - FRAC_PI_2) < TOL {
357                panic!("Coordinate outside projection domain");
358            }
359            y = stere.akm1 * tan(FRAC_PI_4 + 0.5 * p.phi());
360            x = sinlam * y;
361            y *= coslam;
362        }
363        ProjMode::SPole => {
364            if fabs(p.phi() - FRAC_PI_2) < TOL {
365                panic!("Coordinate outside projection domain");
366            }
367            y = stere.akm1 * tan(FRAC_PI_4 + 0.5 * p.phi());
368            x = sinlam * y;
369            y *= coslam;
370        }
371    }
372
373    p.set_x(x);
374    p.set_y(y);
375}
376
377/// Stereographic Ellipsoidal inverse project
378pub fn stere_e_inverse<P: TransformCoordinates>(stere: &mut StereData, proj: &Proj, p: &mut P) {
379    let mut x = p.x();
380    let mut y = p.y();
381    let rho = hypot(p.x(), p.y());
382    let mut phi_l;
383    let mut tp;
384    let halfpi;
385    let halfe;
386
387    match stere.mode {
388        ProjMode::Obliq | ProjMode::Equit => {
389            tp = 2. * atan2(rho * stere.cos_x1, stere.akm1);
390            let cosphi = cos(tp);
391            let sinphi = sin(tp);
392            if rho == 0.0 {
393                phi_l = asin(cosphi * stere.sin_x1);
394            } else {
395                phi_l = asin(cosphi * stere.sin_x1 + (y * sinphi * stere.cos_x1 / rho));
396            }
397
398            tp = tan(0.5 * (FRAC_PI_2 + phi_l));
399            x *= sinphi;
400            y = rho * stere.cos_x1 * cosphi - y * stere.sin_x1 * sinphi;
401            halfpi = FRAC_PI_2;
402            halfe = 0.5 * proj.e;
403        }
404        ProjMode::NPole => {
405            y = -y;
406            tp = -rho / stere.akm1;
407            phi_l = FRAC_PI_2 - 2. * atan(tp);
408            halfpi = -FRAC_PI_2;
409            halfe = -0.5 * proj.e;
410        }
411        ProjMode::SPole => {
412            tp = -rho / stere.akm1;
413            phi_l = FRAC_PI_2 - 2. * atan(tp);
414            halfpi = -FRAC_PI_2;
415            halfe = -0.5 * proj.e;
416        }
417    }
418
419    let mut i = NITER;
420    while i > 0 {
421        let sinphi = proj.e * sin(phi_l);
422        p.set_phi(2. * atan(tp * pow((1. + sinphi) / (1. - sinphi), halfe)) - halfpi);
423        if fabs(phi_l - p.phi()) < CONV {
424            if stere.mode == ProjMode::SPole {
425                p.set_phi(-p.phi());
426            }
427            p.set_lam(if x == 0. && y == 0. { 0. } else { atan2(x, y) });
428            return;
429        }
430        phi_l = p.phi();
431        i -= 1;
432    }
433
434    panic!("Coordinate outside projection domain");
435}
436
437/// Stereographic Spheroidal inverse project
438pub fn stere_s_inverse<P: TransformCoordinates>(stere: &mut StereData, proj: &Proj, p: &mut P) {
439    let rh = hypot(p.x(), p.y());
440    let mut c = 2. * atan(rh / stere.akm1);
441    let sinc = sin(c);
442    let cosc = cos(c);
443    let mut lam = 0.;
444    let phi;
445
446    match stere.mode {
447        ProjMode::Equit => {
448            if fabs(rh) <= EPS10 {
449                phi = 0.;
450            } else {
451                phi = asin(p.y() * sinc / rh);
452            }
453            if cosc != 0. || p.x() != 0. {
454                lam = atan2(p.x() * sinc, cosc * rh);
455            }
456        }
457        ProjMode::Obliq => {
458            if fabs(rh) <= EPS10 {
459                phi = proj.phi0;
460            } else {
461                phi = asin(cosc * stere.sin_x1 + p.y() * sinc * stere.cos_x1 / rh);
462            }
463            c = cosc - stere.sin_x1 * sin(phi);
464            if c != 0. || p.x() != 0. {
465                lam = atan2(p.x() * sinc * stere.cos_x1, c * rh);
466            }
467        }
468        ProjMode::NPole => {
469            p.set_y(-p.y());
470            if fabs(rh) <= EPS10 {
471                phi = proj.phi0;
472            } else {
473                phi = asin(if stere.mode == ProjMode::SPole { -cosc } else { cosc });
474            }
475            lam = if p.x() == 0. && p.y() == 0. { 0. } else { atan2(p.x(), p.y()) };
476        }
477        ProjMode::SPole => {
478            if fabs(rh) <= EPS10 {
479                phi = proj.phi0;
480            } else {
481                phi = asin(if stere.mode == ProjMode::SPole { -cosc } else { cosc });
482            }
483            lam = if p.x() == 0. && p.y() == 0. { 0. } else { atan2(p.x(), p.y()) };
484        }
485    }
486
487    p.set_lam(lam);
488    p.set_phi(phi);
489}
gistools/proj/project/stere.rs

gistools/proj/project/
stere.rs
