miniproj_ops/ops/
stereographic.rs

1//This file is licensed under EUPL v1.2 as part of the Digital Earth Viewer
2
3use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
4
5use crate::{ellipsoid::Ellipsoid, traits::GetterContstruct, DbContstruct, PseudoSerialize};
6
7#[derive(Copy, Clone, Debug)]
8pub struct PolarStereographicAParams {
9    /// longitude of natural origin
10    lon_orig: f64,
11    /// latitude of natural origin
12    lat_orig: f64,
13    /// scale factor at natural origin
14    k_orig: f64,
15    /// false easting
16    false_e: f64,
17    /// false northing
18    false_n: f64,
19}
20
21impl PolarStereographicAParams {
22    pub const fn new(
23        lon_orig: f64,
24        lat_orig: f64,
25        k_orig: f64,
26        false_e: f64,
27        false_n: f64,
28    ) -> Self {
29        Self {
30            lat_orig,
31            lon_orig,
32            k_orig,
33            false_e,
34            false_n,
35        }
36    }
37
38    /// longitude of natural origin, radians
39    pub fn lon_orig(&self) -> f64 {
40        self.lon_orig
41    }
42
43    /// latitude of natural origin, radians
44    pub fn lat_orig(&self) -> f64 {
45        self.lat_orig
46    }
47
48    /// scale factor at natural origin
49    pub fn k_orig(&self) -> f64 {
50        self.k_orig
51    }
52
53    /// false easting
54    pub fn false_e(&self) -> f64 {
55        self.false_e
56    }
57
58    /// false northing
59    pub fn false_n(&self) -> f64 {
60        self.false_n
61    }
62}
63
64/// Polar Stereographic coordinate operation.
65#[derive(Copy, Clone, Debug)]
66pub struct PolarStereographicAProjection {
67    pub t_rho_factor: f64,
68
69    pub phi_2_chi_sin_summand_factor: f64,
70    pub phi_4_chi_sin_summand_factor: f64,
71    pub phi_6_chi_sin_summand_factor: f64,
72    pub phi_8_chi_sin_summand_factor: f64,
73
74    //params: &'b PolarStereographicAParams,
75    pub lat_orig: f64,
76    pub lon_orig: f64,
77    pub false_e: f64,
78    pub false_n: f64,
79    //ell: &'a Ellipsoid,
80    pub ell_e: f64,
81}
82
83impl PolarStereographicAProjection {
84    pub fn new(ell: &Ellipsoid, params: &PolarStereographicAParams) -> Self {
85        let t_rho_factor =
86            ((1.0 + ell.e()).powf(1.0 + ell.e()) * (1.0 - ell.e()).powf(1.0 - ell.e())).sqrt()
87                / (2.0 * ell.a() * params.k_orig());
88        let phi_2_chi_sin_summand_factor = ell.e_squared() / 2.0
89            + 5.0 * ell.e_squared().powi(2) / 24.0
90            + ell.e_squared().powi(3) / 12.0
91            + 13.0 * ell.e_squared().powi(4) / 360.0;
92        let phi_4_chi_sin_summand_factor = 7.0 * ell.e_squared().powi(2) / 48.0
93            + 29.0 * ell.e_squared().powi(3) / 240.0
94            + ell.e_squared().powi(4) / 11520.0;
95        let phi_6_chi_sin_summand_factor =
96            7.0 * ell.e_squared().powi(3) / 120.0 + 81.0 * ell.e_squared().powi(4) / 1120.0;
97        let phi_8_chi_sin_summand_factor = 4279.0 * ell.e_squared().powi(4) / 161280.0;
98        Self {
99            t_rho_factor,
100            phi_2_chi_sin_summand_factor,
101            phi_4_chi_sin_summand_factor,
102            phi_6_chi_sin_summand_factor,
103            phi_8_chi_sin_summand_factor,
104
105            lat_orig: params.lat_orig(),
106            lon_orig: params.lon_orig(),
107            false_e: params.false_e(),
108            false_n: params.false_n(),
109
110            ell_e: ell.e(),
111        }
112    }
113}
114
115impl crate::traits::Projection for PolarStereographicAProjection {
116    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
117        if self.lat_orig < 0.0 {
118            // North Pole Case
119            let t = f64::tan(std::f64::consts::FRAC_PI_4 - latitude / 2.0)
120                * ((1.0 + self.ell_e * latitude.sin()) / (1.0 - self.ell_e * latitude.sin()))
121                    .powf(self.ell_e / 2.0);
122            let rho = t / self.t_rho_factor;
123            (
124                self.false_e + rho * f64::sin(longitude - self.lon_orig),
125                self.false_n - rho * f64::cos(longitude - self.lon_orig),
126            )
127        } else {
128            // South Pole Case
129            let t = f64::tan(std::f64::consts::FRAC_PI_4 + latitude / 2.0)
130                / ((1.0 + self.ell_e * latitude.sin()) / (1.0 - self.ell_e * latitude.sin()))
131                    .powf(self.ell_e / 2.0);
132            let rho = t / self.t_rho_factor;
133            (
134                self.false_e + rho * f64::sin(longitude - self.lon_orig),
135                self.false_n + rho * f64::cos(longitude - self.lon_orig),
136            )
137        }
138    }
139
140    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
141        let rho_ = ((easting - self.false_e).powi(2) + (northing - self.false_n).powi(2)).sqrt();
142        let t_ = rho_ * self.t_rho_factor;
143        let chi = if self.lat_orig < 0.0 {
144            // North Pole Case
145            FRAC_PI_2 - 2.0 * t_.atan()
146        } else {
147            // South Pole Case
148            2.0 * t_.atan() - FRAC_PI_2
149        };
150        let phi = chi
151            + self.phi_2_chi_sin_summand_factor * (2.0 * chi).sin()
152            + self.phi_4_chi_sin_summand_factor * (4.0 * chi).sin()
153            + self.phi_6_chi_sin_summand_factor * (6.0 * chi).sin()
154            + self.phi_8_chi_sin_summand_factor * (8.0 * chi).sin();
155        let lambda = /*if easting == self.false_e { //this appears wrong to me so it's commented out. @ me if you think it's right tho.
156            self.lat_orig
157        } else*/ if self.lat_orig < 0.0 { // North Pole Case
158            self.lon_orig + (easting - self.false_e).atan2(self.false_n - northing)
159        } else { // South Pole Case
160            self.lon_orig + (easting - self.false_e).atan2(northing - self.false_n)
161        };
162        (lambda, phi)
163    }
164}
165impl DbContstruct for PolarStereographicAProjection {
166    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
167        let params = PolarStereographicAParams::new(
168            params
169                .iter()
170                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
171                .unwrap(),
172            params
173                .iter()
174                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
175                .unwrap(),
176            params
177                .iter()
178                .find_map(|(c, v)| if *c == 8805 { Some(*v) } else { None })
179                .unwrap(),
180            params
181                .iter()
182                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
183                .unwrap(),
184            params
185                .iter()
186                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
187                .unwrap(),
188        );
189        Self::new(ellipsoid, &params)
190    }
191}
192
193impl GetterContstruct for PolarStereographicAProjection {
194    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
195    where
196        G: FnMut(u32) -> Option<f64>,
197    {
198        let params = PolarStereographicAParams::new(
199            getter(8802)?,
200            getter(8801)?,
201            getter(8805)?,
202            getter(8806)?,
203            getter(8807)?,
204        );
205        Some(Self::new(ellipsoid, &params))
206    }
207}
208
209impl PseudoSerialize for PolarStereographicAProjection {
210    fn to_constructed(&self) -> String {
211        format!(
212            r"PolarStereographicAProjection{{
213    t_rho_factor: {}f64,
214    phi_2_chi_sin_summand_factor: {}f64,
215    phi_4_chi_sin_summand_factor: {}f64,
216    phi_6_chi_sin_summand_factor: {}f64,
217    phi_8_chi_sin_summand_factor: {}f64,
218    
219    lat_orig: {}f64,
220    lon_orig: {}f64,
221    false_e: {}f64,
222    false_n: {}f64,
223
224    ell_e: {}f64
225}}",
226            self.t_rho_factor,
227            self.phi_2_chi_sin_summand_factor,
228            self.phi_4_chi_sin_summand_factor,
229            self.phi_6_chi_sin_summand_factor,
230            self.phi_8_chi_sin_summand_factor,
231            self.lat_orig,
232            self.lon_orig,
233            self.false_e,
234            self.false_n,
235            self.ell_e
236        )
237    }
238}
239pub fn direct_projection_a(params: &[(u32, f64)], ell: Ellipsoid) -> String {
240    PolarStereographicAProjection::from_database_params(params, &ell).to_constructed()
241}
242
243#[derive(Copy, Clone, Debug)]
244pub struct ObliqueStereographicParams {
245    // Longitude of natural origin
246    lon_orig: f64,
247    // Latitude of natural origin
248    lat_orig: f64,
249    // Scale factor at natural origin
250    k_orig: f64,
251    // False easting
252    false_e: f64,
253    // False northing
254    false_n: f64,
255}
256
257impl ObliqueStereographicParams {
258    pub fn new(lon_orig: f64, lat_orig: f64, k_orig: f64, false_e: f64, false_n: f64) -> Self {
259        assert!(lat_orig > 0f64);
260        Self {
261            lat_orig,
262            lon_orig,
263            k_orig,
264            false_e,
265            false_n,
266        }
267    }
268
269    /// longitude of natural origin, radians
270    pub fn lon_orig(&self) -> f64 {
271        self.lon_orig
272    }
273
274    /// latitude of natural origin, radians
275    pub fn lat_orig(&self) -> f64 {
276        self.lat_orig
277    }
278
279    /// scale factor at natural origin
280    pub fn k_orig(&self) -> f64 {
281        self.k_orig
282    }
283
284    /// false easting
285    pub fn false_e(&self) -> f64 {
286        self.false_e
287    }
288
289    /// false northing
290    pub fn false_n(&self) -> f64 {
291        self.false_n
292    }
293}
294
295#[allow(non_snake_case)]
296#[derive(Copy, Clone, Debug)]
297pub struct ObliqueStereographicProjection {
298    pub false_e: f64,
299    pub false_n: f64,
300    pub chi_O: f64,
301    pub R_k_O_2: f64,
302    pub c: f64,
303    pub ellipsoid_e: f64,
304    pub ellipsoid_e_sq: f64,
305    pub n: f64,
306    pub lon_orig: f64,
307    pub g: f64,
308    pub h: f64,
309}
310
311impl ObliqueStereographicProjection {
312    const MAX_ITERATIONS: usize = 4;
313
314    #[allow(non_snake_case)]
315    pub fn new(ell: &Ellipsoid, params: &ObliqueStereographicParams) -> Self {
316        let rho_O = ell.rho(params.lat_orig());
317        let ny_O = ell.ny(params.lat_orig());
318        let R = (rho_O * ny_O).sqrt();
319        let n = (1f64
320            + ((ell.e_squared() * params.lat_orig().cos().powi(4)) / (1f64 - ell.e_squared())))
321        .sqrt();
322        let S_1 = (1f64 + params.lat_orig().sin()) / (1f64 - params.lat_orig.sin());
323        let S_2 =
324            (1f64 - ell.e() * params.lat_orig().sin()) / (1f64 + ell.e() * params.lat_orig().sin());
325        let w_1 = (S_1 * S_2.powf(ell.e())).powf(n);
326        let chi_OO_sin = (w_1 - 1f64) / (w_1 + 1f64);
327        let c = (n + params.lat_orig().sin()) * (1f64 - chi_OO_sin)
328            / ((n - params.lat_orig().sin()) * (1f64 + chi_OO_sin));
329        let w_2 = c * w_1;
330        let chi_O = ((w_2 - 1f64) / (w_2 + 1f64)).asin();
331        let g = 2f64 * R * params.k_orig() * (FRAC_PI_4 - chi_O / 2f64).tan();
332        let h = 4f64 * R * params.k_orig() * chi_O.tan() + g;
333        Self {
334            false_e: params.false_e(),
335            false_n: params.false_n(),
336            chi_O,
337            R_k_O_2: R * params.k_orig() * 2f64,
338            c,
339            ellipsoid_e: ell.e(),
340            ellipsoid_e_sq: ell.e_squared(),
341            n,
342            lon_orig: params.lon_orig(),
343            g,
344            h,
345        }
346    }
347}
348
349impl crate::traits::Projection for ObliqueStereographicProjection {
350    #[allow(non_snake_case)]
351    fn projected_to_rad(&self, x: f64, y: f64) -> (f64, f64) {
352        let i = (x - self.false_e).atan2(self.h + (y - self.false_n));
353        let j = (x - self.false_e).atan2(self.g - (y - self.false_n)) - i;
354        let chi = self.chi_O
355            + 2f64
356                * (((y - self.false_n) - (x - self.false_e) * (j / 2f64).tan()) / self.R_k_O_2)
357                    .atan();
358        let psi = 0.5 * ((1f64 + chi.sin()) / (self.c * (1f64 - chi.sin()))).ln() / self.n;
359        let mut phi = 2f64 * psi.exp().atan() - FRAC_PI_2;
360        for _ in 0..Self::MAX_ITERATIONS {
361            let psi_ = ((phi / 2f64 + FRAC_PI_4).tan()
362                * ((1f64 - self.ellipsoid_e * phi.sin()) / (1f64 + self.ellipsoid_e * phi.sin()))
363                    .powf(self.ellipsoid_e / 2f64))
364            .ln();
365            phi = phi
366                - (psi_ - psi) * phi.cos() * (1f64 - self.ellipsoid_e_sq * phi.sin().powi(2))
367                    / (1f64 - self.ellipsoid_e_sq);
368        }
369        let DeltaLambda = j + 2f64 * i;
370        (DeltaLambda / self.n + self.lon_orig, phi)
371    }
372
373    #[allow(non_snake_case)]
374    fn rad_to_projected(&self, lon: f64, lat: f64) -> (f64, f64) {
375        let S_a = (1f64 + lat.sin()) / (1f64 - lat.sin());
376        let S_b = (1f64 - self.ellipsoid_e * lat.sin()) / (1f64 + self.ellipsoid_e * lat.sin());
377        let DeltaLambda = self.n * (lon - self.lon_orig);
378        let w = self.c * (S_a * S_b.powf(self.ellipsoid_e)).powf(self.n);
379        let chi = ((w - 1f64) / (w + 1f64)).asin();
380        let B =
381            1f64 + chi.sin() * self.chi_O.sin() + chi.cos() * self.chi_O.cos() * DeltaLambda.cos();
382        (
383            self.false_e + self.R_k_O_2 * chi.cos() * DeltaLambda.sin() / B,
384            self.false_n
385                + self.R_k_O_2
386                    * (chi.sin() * self.chi_O.cos()
387                        - chi.cos() * self.chi_O.sin() * DeltaLambda.cos())
388                    / B,
389        )
390    }
391}
392
393impl DbContstruct for ObliqueStereographicProjection {
394    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
395        let params = ObliqueStereographicParams::new(
396            params
397                .iter()
398                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
399                .unwrap(),
400            params
401                .iter()
402                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
403                .unwrap(),
404            params
405                .iter()
406                .find_map(|(c, v)| if *c == 8805 { Some(*v) } else { None })
407                .unwrap(),
408            params
409                .iter()
410                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
411                .unwrap(),
412            params
413                .iter()
414                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
415                .unwrap(),
416        );
417        Self::new(ellipsoid, &params)
418    }
419}
420
421impl GetterContstruct for ObliqueStereographicProjection {
422    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
423    where
424        G: FnMut(u32) -> Option<f64>,
425    {
426        let params = ObliqueStereographicParams::new(
427            getter(8802)?,
428            getter(8801)?,
429            getter(8805)?,
430            getter(8806)?,
431            getter(8807)?,
432        );
433        Some(Self::new(ellipsoid, &params))
434    }
435}
436
437impl PseudoSerialize for ObliqueStereographicProjection {
438    fn to_constructed(&self) -> String {
439        format!(
440            r"ObliqueStereographicProjection{{
441    false_e: {}f64,
442    false_n: {}f64,
443    chi_O: {}f64,
444    R_k_O_2: {}f64,
445    c: {}f64,
446    ellipsoid_e: {}f64,
447    ellipsoid_e_sq: {}f64,
448    n: {}f64,
449    lon_orig: {}f64,
450    g: {}f64,
451    h: {}f64
452}}",
453            self.false_e,
454            self.false_n,
455            self.chi_O,
456            self.R_k_O_2,
457            self.c,
458            self.ellipsoid_e,
459            self.ellipsoid_e_sq,
460            self.n,
461            self.lon_orig,
462            self.g,
463            self.h
464        )
465    }
466}
467
468pub fn direct_projection_oblique(params: &[(u32, f64)], ell: Ellipsoid) -> String {
469    ObliqueStereographicProjection::from_database_params(params, &ell).to_constructed()
470}
471#[cfg(test)]
472mod tests {
473
474    use crate::ellipsoid::Ellipsoid;
475    use crate::stereographic::*;
476    use crate::traits::*;
477
478    #[test]
479    fn polar_stereographic_a_consistency() {
480        let ell = Ellipsoid::from_a_f_inv(6378137.0, 298.257223563);
481        let params = PolarStereographicAParams::new(
482            0.0f64.to_radians(),
483            -90.0f64.to_radians(),
484            0.994,
485            2_000_000.0,
486            2_000_000.0,
487        );
488
489        let projection = PolarStereographicAProjection::new(&ell, &params);
490        let easting_goal = 3329416.75;
491        let northing_goal = 632668.43;
492        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
493        let (easting, northing) = projection.deg_to_projected(lon, lat);
494
495        eprintln!("easting: {easting_goal} - {easting}");
496        eprintln!("northing: {northing_goal} - {northing}");
497
498        assert!((easting - easting_goal).abs() < 0.01);
499
500        assert!((northing - northing_goal).abs() < 0.01);
501    }
502
503    #[test]
504    fn oblique_stereographic_consistency() {
505        let ell = Ellipsoid::from_a_f_inv(6377397.155, 299.15281);
506        let params = ObliqueStereographicParams::new(
507            0.094032038,
508            0.910296727,
509            0.9999079,
510            155000.0,
511            463000.0,
512        );
513
514        let projection = ObliqueStereographicProjection::new(&ell, &params);
515        let easting_goal = 196105.283;
516        let northing_goal = 557057.739;
517        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
518        eprintln!("lon: {lon}, lat: {lat}");
519        let (easting, northing) = projection.deg_to_projected(lon, lat);
520
521        eprintln!("easting: {easting_goal} - {easting}");
522        eprintln!("northing: {northing_goal} - {northing}");
523
524        assert!((easting - easting_goal).abs() < 0.01);
525
526        assert!((northing - northing_goal).abs() < 0.01);
527    }
528}
miniproj_ops/ops/stereographic.rs

miniproj_ops/ops/
stereographic.rs
