gistools/proj/project/
somerc.rs

1use crate::proj::{CoordinateStep, Proj, ProjectCoordinates, SOMERC, TransformCoordinates, aasin};
2use alloc::rc::Rc;
3use core::{
4    cell::RefCell,
5    f64::consts::{FRAC_PI_2, FRAC_PI_4},
6};
7use libm::{atan, cos, exp, fabs, log, sin, sqrt, tan};
8
9/// Swiss Oblique Mercator variables
10#[derive(Debug, Default, Clone, PartialEq)]
11pub struct SomercData {
12    k: f64,
13    c: f64,
14    hlf_e: f64,
15    k_r: f64,
16    cosp0: f64,
17    sinp0: f64,
18}
19
20const EPS: f64 = 1.0e-10;
21const NITER: usize = 6;
22
23/// # Swiss Oblique Mercator
24///
25/// **Classification**: Oblique Mercator
26///
27/// **Available forms**: Forward and inverse, ellipsoidal only
28///
29/// **Defined area**: Global
30///
31/// **Alias**: somerc
32///
33/// **Domain**: 2D
34///
35/// **Input type**: Geodetic coordinates
36///
37/// **Output type**: Projected coordinates
38///
39/// ## Projection String
40/// ```ini
41/// +proj=somerc
42/// ```
43///
44/// ## Required Parameters
45/// - None
46///
47/// ## Optional Parameters
48/// - `+lon_0=<value>`: Central meridian.
49/// - `+ellps=<value>`: Ellipsoid used.
50/// - `+R=<value>`: Radius of the projection sphere.
51/// - `+k_0=<value>`: Scale factor.
52/// - `+x_0=<value>`: False easting.
53/// - `+y_0=<value>`: False northing.
54///
55/// ## References:
56/// Formules et constantes pour le Calcul pour la
57/// projection cylindrique conforme à axe oblique et pour la transformation entre
58/// des systèmes de référence.
59///
60/// <http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf>
61///
62/// ![Swiss Oblique Mercator](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/somerc.png?raw=true)
63#[derive(Debug, Clone, PartialEq)]
64pub struct SwissOblMercatorProjection {
65    proj: Rc<RefCell<Proj>>,
66    store: RefCell<SomercData>,
67}
68impl ProjectCoordinates for SwissOblMercatorProjection {
69    fn code(&self) -> i64 {
70        SOMERC
71    }
72    fn name(&self) -> &'static str {
73        "Swiss. Obl. Mercator"
74    }
75    fn names() -> &'static [&'static str] {
76        &["Swiss. Obl. Mercator", "somerc"]
77    }
78}
79impl CoordinateStep for SwissOblMercatorProjection {
80    fn new(proj: Rc<RefCell<Proj>>) -> Self {
81        let mut store = SomercData::default();
82        {
83            let proj = &mut proj.borrow_mut();
84            store.hlf_e = 0.5 * proj.e;
85            let mut cp = cos(proj.phi0);
86            cp *= cp;
87            store.c = sqrt(1. + proj.es * cp * cp * proj.rone_es);
88            let mut sp = sin(proj.phi0);
89            store.sinp0 = sp / store.c;
90            let phip0 = aasin(store.sinp0);
91            store.cosp0 = cos(phip0);
92            sp *= proj.e;
93            store.k = log(tan(FRAC_PI_4 + 0.5 * phip0))
94                - store.c
95                    * (log(tan(FRAC_PI_4 + 0.5 * proj.phi0))
96                        - store.hlf_e * log((1. + sp) / (1. - sp)));
97            store.k_r = proj.k0 * sqrt(proj.one_es) / (1. - sp * sp);
98        }
99        SwissOblMercatorProjection { proj, store: store.into() }
100    }
101    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
102        somerc_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
103    }
104    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
105        somerc_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
106    }
107}
108
109/// Equal Earth Ellipsoidal forward project
110pub fn somerc_e_forward<P: TransformCoordinates>(somerc: &SomercData, proj: &Proj, p: &mut P) {
111    let sp = proj.e * sin(p.phi());
112    let phip = 2.
113        * atan(exp(somerc.c
114            * (log(tan(FRAC_PI_4 + 0.5 * p.phi())) - somerc.hlf_e * log((1. + sp) / (1. - sp)))
115            + somerc.k))
116        - FRAC_PI_2;
117    let lamp = somerc.c * p.lam();
118    let cp = cos(phip);
119    let phipp = aasin(somerc.cosp0 * sin(phip) - somerc.sinp0 * cp * cos(lamp));
120    let lampp = aasin(cp * sin(lamp) / cos(phipp));
121    p.set_x(somerc.k_r * lampp);
122    p.set_y(somerc.k_r * log(tan(FRAC_PI_4 + 0.5 * phipp)));
123}
124
125/// Equal Earth Ellipsoidal forward project
126pub fn somerc_e_inverse<P: TransformCoordinates>(somerc: &SomercData, proj: &Proj, p: &mut P) {
127    let phipp = 2. * (atan(exp(p.y() / somerc.k_r)) - FRAC_PI_4);
128    let lampp = p.x() / somerc.k_r;
129    let cp = cos(phipp);
130    let mut phip = aasin(somerc.cosp0 * sin(phipp) + somerc.sinp0 * cp * cos(lampp));
131    let lamp = aasin(cp * sin(lampp) / cos(phip));
132    let con = (somerc.k - log(tan(FRAC_PI_4 + 0.5 * phip))) / somerc.c;
133    let mut i = NITER;
134    while i > 0 {
135        let esp = proj.e * sin(phip);
136        let delp = (con + log(tan(FRAC_PI_4 + 0.5 * phip))
137            - somerc.hlf_e * log((1. + esp) / (1. - esp)))
138            * (1. - esp * esp)
139            * cos(phip)
140            * proj.rone_es;
141        phip -= delp;
142        if fabs(delp) < EPS {
143            break;
144        }
145        i -= 1;
146    }
147    if i != 0 {
148        p.set_phi(phip);
149        p.set_lam(lamp / somerc.c);
150    } else {
151        panic!("Coordinate outside projection domain");
152    }
153}
gistools/proj/project/somerc.rs

gistools/proj/project/
somerc.rs
