gistools/proj/project/
sterea.rs

1use crate::proj::{
2    CoordinateStep, OBLIQUE_STEREOGRAPHIC, Proj, ProjectCoordinates, TransformCoordinates,
3};
4use alloc::rc::Rc;
5use core::{
6    cell::RefCell,
7    f64::consts::{FRAC_PI_2, FRAC_PI_4},
8};
9use libm::{asin, atan, atan2, cos, fabs, hypot, pow, sin, sqrt, tan};
10
11/// Gaussian Variables
12#[derive(Debug, Default, Clone, PartialEq)]
13struct Gauss {
14    c: f64,
15    k: f64,
16    e: f64,
17    ratexp: f64,
18}
19
20/// Oblique Stereographic Alternative Variables
21#[derive(Debug, Default, Clone, PartialEq)]
22pub struct StereaData {
23    phic0: f64,
24    cosc0: f64,
25    sinc0: f64,
26    r2: f64,
27    en: Gauss,
28}
29
30const MAX_ITER: usize = 20;
31const DEL_TOL: f64 = 1e-14;
32
33/// # Oblique Stereographic Alternative
34///
35/// **Classification**: Azimuthal
36///
37/// **Available forms**: Forward and inverse, spherical and ellipsoidal
38///
39/// **Defined area**: Global
40///
41/// **Alias**: sterea
42///
43/// **Domain**: 2D
44///
45/// **Input type**: Geodetic coordinates
46///
47/// **Output type**: Projected coordinates
48///
49/// ## Projection String
50/// ```ini
51/// +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel
52/// ```
53///
54/// ## Note
55/// This projection method, referenced by EPSG as "Oblique Stereographic", is
56/// for example used for the Netherlands "Amersfoort / RD New" projected CRS.
57/// It gives different results than the :ref:`stere` method in the non-polar cases
58/// (i.e. the oblique and equatorial case).
59///
60/// ## Required Parameters
61/// - None
62///
63/// ## Optional Parameters
64/// - `+lat_0=<value>`: Latitude of origin.
65/// - `+lon_0=<value>`: Central meridian.
66/// - `+k=<value>`: Scale factor.
67/// - `+x_0=<value>`: False easting.
68/// - `+y_0=<value>`: False northing.
69/// - `+ellps=<value>`: Ellipsoid used.
70/// - `+R=<value>`: Radius of the projection sphere.
71///
72/// ![Oblique Stereographic Alternative](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/sterea.png?raw=true)
73#[derive(Debug, Clone, PartialEq)]
74pub struct ObliqueStereographicAlternativeProjection {
75    proj: Rc<RefCell<Proj>>,
76    store: RefCell<StereaData>,
77}
78impl ProjectCoordinates for ObliqueStereographicAlternativeProjection {
79    fn code(&self) -> i64 {
80        OBLIQUE_STEREOGRAPHIC
81    }
82    fn name(&self) -> &'static str {
83        "Oblique Stereographic Alternative"
84    }
85    fn names() -> &'static [&'static str] {
86        &["Oblique Stereographic Alternative", "Stereographic_North_Pole", "sterea"]
87    }
88}
89impl CoordinateStep for ObliqueStereographicAlternativeProjection {
90    fn new(proj: Rc<RefCell<Proj>>) -> Self {
91        let mut store = StereaData::default();
92        {
93            let proj = &mut proj.borrow_mut();
94            let mut r: f64 = 0.;
95
96            store.en = gauss_ini(proj.e, proj.phi0, &mut store.phic0, &mut r);
97            store.sinc0 = sin(store.phic0);
98            store.cosc0 = cos(store.phic0);
99            store.r2 = 2. * r;
100        }
101        ObliqueStereographicAlternativeProjection { proj, store: store.into() }
102    }
103    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
104        sterea_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
105    }
106    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
107        sterea_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
108    }
109}
110
111/// Oblique Stereographic Alternative Ellipsoidal forward project
112pub fn sterea_e_forward<P: TransformCoordinates>(sterea: &mut StereaData, proj: &Proj, p: &mut P) {
113    gauss(p, &sterea.en);
114    let sinc = sin(p.phi());
115    let cosc = cos(p.phi());
116    let cosl = cos(p.lam());
117    let denom = 1. + sterea.sinc0 * sinc + sterea.cosc0 * cosc * cosl;
118    if denom == 0.0 {
119        panic!("Coordinate outside projection domain");
120    }
121    let k = proj.k0 * sterea.r2 / denom;
122    p.set_x(k * cosc * sin(p.lam()));
123    p.set_y(k * (sterea.cosc0 * sinc - sterea.sinc0 * cosc * cosl));
124}
125
126/// Oblique Stereographic Alternative Ellipsoidal inverse project
127pub fn sterea_e_inverse<P: TransformCoordinates>(sterea: &mut StereaData, proj: &Proj, p: &mut P) {
128    let x = p.x() / proj.k0;
129    let y = p.y() / proj.k0;
130    let rho = hypot(x, y);
131    if rho != 0.0 {
132        let c = 2. * atan2(rho, sterea.r2);
133        let sinc = sin(c);
134        let cosc = cos(c);
135        p.set_phi(asin(cosc * sterea.sinc0 + y * sinc * sterea.cosc0 / rho));
136        p.set_lam(atan2(x * sinc, rho * sterea.cosc0 * cosc - y * sterea.sinc0 * sinc));
137    } else {
138        p.set_phi(sterea.phic0);
139        p.set_lam(0.);
140    }
141    inv_gauss(p, &sterea.en);
142}
143
144fn srat(esinp: f64, ratexp: f64) -> f64 {
145    pow((1. - esinp) / (1. + esinp), ratexp)
146}
147
148fn gauss_ini(e: f64, phi0: f64, chi: &mut f64, rc: &mut f64) -> Gauss {
149    let mut en = Gauss::default();
150
151    let es = e * e;
152    en.e = e;
153    let sphi = sin(phi0);
154    let mut cphi = cos(phi0);
155    cphi *= cphi;
156    *rc = sqrt(1. - es) / (1. - es * sphi * sphi);
157    en.c = sqrt(1. + es * cphi * cphi / (1. - es));
158    if en.c == 0.0 {
159        panic!("Failed to initialize Gauss projection");
160    }
161    *chi = asin(sphi / en.c);
162    en.ratexp = 0.5 * en.c * e;
163    let srat_val = srat(en.e * sphi, en.ratexp);
164    if srat_val == 0.0 {
165        panic!("Failed to initialize Gauss projection");
166    }
167    if 0.5 * phi0 + FRAC_PI_4 < 1e-10 {
168        en.k = 1.0 / srat_val;
169    } else {
170        en.k = tan(0.5 * *chi + FRAC_PI_4) / (pow(tan(0.5 * phi0 + FRAC_PI_4), en.c) * srat_val);
171    }
172
173    en
174}
175
176fn gauss<P: TransformCoordinates>(elp: &mut P, en: &Gauss) {
177    elp.set_phi(
178        2. * atan(
179            en.k * pow(tan(0.5 * elp.phi() + FRAC_PI_4), en.c)
180                * srat(en.e * sin(elp.phi()), en.ratexp),
181        ) - FRAC_PI_2,
182    );
183    elp.set_lam(en.c * (elp.lam()));
184}
185
186fn inv_gauss<P: TransformCoordinates>(p: &mut P, en: &Gauss) {
187    p.set_lam(p.lam() / en.c);
188    let num = pow(tan(0.5 * p.phi() + FRAC_PI_4) / en.k, 1. / en.c);
189    let mut i = MAX_ITER;
190    while i > 0 {
191        p.set_phi(2. * atan(num * srat(en.e * sin(p.phi()), -0.5 * en.e)) - FRAC_PI_2);
192        if fabs(p.phi() - p.phi()) < DEL_TOL {
193            break;
194        }
195        p.set_phi(p.phi());
196        i -= 1;
197    }
198    // convergence failed
199    if i == 0 {
200        panic!("Coordinate outside projection domain");
201    }
202}
gistools/proj/project/sterea.rs

gistools/proj/project/
sterea.rs
