gistools/proj/project/
cass.rs

1use crate::proj::{
2    CASSINI, CoordinateStep, HYPERBOLIC, Proj, ProjMethod, ProjectCoordinates,
3    TransformCoordinates, enfn, generic_inverse_2d, inv_mlfn, mlfn,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::cell::RefCell;
7use libm::{asin, atan2, cos, sin, sqrt, tan};
8
9const C1: f64 = 0.166_666_666_666_666_66;
10const C2: f64 = 0.008_333_333_333_333_333;
11const C3: f64 = 0.041_666_666_666_666_664;
12const C4: f64 = 0.333_333_333_333_333_3;
13const C5: f64 = 0.066_666_666_666_666_67;
14
15/// Cassini variables
16#[derive(Debug, Default, Clone, PartialEq)]
17pub struct CassData {
18    en: Vec<f64>,
19    m0: f64,
20    hyperbolic: bool,
21}
22
23/// # Cassini (Cassini-Soldner) Projection
24///
25/// Although the Cassini projection has been largely replaced by the Transverse Mercator, it is
26/// still in limited use outside the United States and was one of the major topographic mapping
27/// projections until the early 20th century.
28///
29/// **Classification**: Transverse and oblique cylindrical
30///
31/// **Available forms**: Forward and Inverse, Spherical and ellipsoidal
32///
33/// **Defined area**: Global, but best used near the central meridian with long, narrow areas
34///
35/// **Alias**: `cass`
36///
37/// **Domain**: 2D
38///
39/// **Input type**: Geodetic coordinates
40///
41/// **Output type**: Projected coordinates
42///
43/// ## Projection String
44/// ```ini
45/// +proj=cass
46/// ```
47///
48/// ## Required Parameters
49/// - `lat_0`
50///
51/// ## Optional Parameters
52/// - `lon_0`
53/// - `x_0`
54/// - `y_0`
55/// - `ellps`
56/// - `R`
57/// - `+hyperbolic`: Use modified form of the standard Cassini-Soldner projection known as the Hyperbolic Cassini-Soldner (used in EPSG:3139).
58///
59/// ![Cassini (Cassini-Soldner)](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/cass.png?raw=true)
60#[derive(Debug, Clone, PartialEq)]
61pub struct CassiniProjection {
62    proj: Rc<RefCell<Proj>>,
63    store: RefCell<CassData>,
64    method: ProjMethod,
65}
66impl ProjectCoordinates for CassiniProjection {
67    fn code(&self) -> i64 {
68        CASSINI
69    }
70    fn name(&self) -> &'static str {
71        "Cassini"
72    }
73    fn names() -> &'static [&'static str] {
74        &["Cassini", "Cassini-Soldner", "Cassini_Soldner", "cass"]
75    }
76}
77impl CoordinateStep for CassiniProjection {
78    fn new(proj: Rc<RefCell<Proj>>) -> Self {
79        let mut store = CassData::default();
80        let method: ProjMethod;
81        {
82            let proj = &mut proj.borrow_mut();
83            // Spheroidal?
84            method = if 0. == proj.es {
85                ProjMethod::Spheroidal
86            } else {
87                // otherwise it's ellipsoidal
88                store.en = enfn(proj.n);
89                store.m0 = mlfn(proj.phi0, sin(proj.phi0), cos(proj.phi0), &store.en);
90                if proj.params.contains_key(&HYPERBOLIC) {
91                    store.hyperbolic = true;
92                }
93                ProjMethod::Ellipsoidal
94            };
95        }
96        CassiniProjection { proj, store: store.into(), method }
97    }
98    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
99        if self.method == ProjMethod::Ellipsoidal {
100            cass_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
101        } else {
102            cass_s_forward(&self.proj.borrow(), p);
103        }
104    }
105    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
106        if self.method == ProjMethod::Ellipsoidal {
107            let es = self.proj.borrow().es;
108            cass_e_inverse(self, &self.store.borrow(), es, p);
109        } else {
110            cass_s_inverse(&self.proj.borrow(), p);
111        }
112    }
113}
114
115/// Cassini Ellipsoidal forward project
116pub fn cass_e_forward<P: TransformCoordinates>(cass: &CassData, proj: &Proj, p: &mut P) {
117    let sinphi = sin(p.phi());
118    let cosphi = cos(p.phi());
119    let m = mlfn(p.phi(), sinphi, cosphi, &cass.en);
120
121    let nu_square = 1. / (1. - proj.es * sinphi * sinphi);
122    let nu = sqrt(nu_square);
123    let tanphi = tan(p.phi());
124    let t = tanphi * tanphi;
125    let a = p.lam() * cosphi;
126    let c = proj.es * (cosphi * cosphi) / (1. - proj.es);
127    let a2 = a * a;
128
129    let x = nu * a * (1. - a2 * t * (C1 + (8. - t + 8. * c) * a2 * C2));
130    let mut y = m - cass.m0 + nu * tanphi * a2 * (0.5 + (5. - t + 6. * c) * a2 * C3);
131    if cass.hyperbolic {
132        let rho = nu_square * (1. - proj.es) * nu;
133        y -= y * y * y / (6. * rho * nu);
134    }
135
136    p.set_x(x);
137    p.set_y(y);
138}
139
140/// Cassini Spheroidal forward project
141pub fn cass_s_forward<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
142    p.set_x(asin(cos(p.phi()) * sin(p.lam())));
143    p.set_y(atan2(tan(p.phi()), cos(p.lam())) - proj.phi0);
144}
145
146/// Cassini Ellipsoidal inverse project
147pub fn cass_e_inverse<P: TransformCoordinates>(
148    cass: &CassiniProjection,
149    cass_data: &CassData,
150    es: f64,
151    p: &mut P,
152) {
153    let phi1 = inv_mlfn(cass_data.m0 + p.y(), &cass_data.en);
154    let tanphi1 = tan(phi1);
155    let t1 = tanphi1 * tanphi1;
156    let sinphi1 = sin(phi1);
157    let nu1_square = 1. / (1. - es * sinphi1 * sinphi1);
158    let nu1 = sqrt(nu1_square);
159    let rho1 = nu1_square * (1. - es) * nu1;
160    let d = p.x() / nu1;
161    let d2 = d * d;
162    let mut lp = P::default();
163    lp.set_phi(phi1 - (nu1 * tanphi1 / rho1) * d2 * (0.5 - (1. + 3. * t1) * d2 * C3));
164    lp.set_lam(d * (1. + t1 * d2 * (-C4 + (1. + 3. * t1) * d2 * C5)) / cos(phi1));
165
166    // EPSG guidance note 7-2 suggests a custom approximation for the
167    // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
168    // generic inversion method
169    // Actually use it in the non-hyperbolic case. It enables to make the
170    // 5108.gie roundtripping tests to success, with at most 2 iterations.
171    let delta_xy_tolerance = 1e-12;
172    generic_inverse_2d(p, cass, &mut lp, delta_xy_tolerance);
173
174    p.set_phi(lp.phi());
175    p.set_lam(lp.lam());
176}
177
178/// Cassini Spheroidal inverse project
179pub fn cass_s_inverse<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
180    let dd = p.y() + proj.phi0;
181    p.set_phi(asin(sin(dd) * cos(p.x())));
182    p.set_lam(atan2(tan(p.x()), cos(dd)));
183}
gistools/proj/project/cass.rs

gistools/proj/project/
cass.rs
