gistools/proj/project/
eqdc.rs

1use crate::proj::{
2    _msfn, CoordinateStep, EPS10, EQUIDISTANT_CONIC, LATITUDE_OF_FIRST_STANDARD_PARALLEL,
3    LATITUDE_OF_SECOND_STANDARD_PARALLEL, Proj, ProjValue, ProjectCoordinates,
4    TransformCoordinates, enfn, inv_mlfn, mlfn,
5};
6use alloc::{rc::Rc, vec::Vec};
7use core::{cell::RefCell, f64::consts::FRAC_PI_2};
8use libm::{atan2, cos, fabs, hypot, sin};
9
10/// Equidistant Conic Variables
11#[derive(Debug, Default, Clone, PartialEq)]
12pub struct EqdcData {
13    phi1: f64,
14    phi2: f64,
15    n: f64,
16    rho: f64,
17    rho0: f64,
18    c: f64,
19    en: Vec<f64>,
20    ellips: bool,
21}
22
23/// # Equidistant Conic
24///
25/// **Classification**: Conic
26///
27/// **Available forms**: Forward and inverse, ellipsoidal
28///
29/// **Defined area**: Global
30///
31/// **Alias**: eqdc
32///
33/// **Domain**: 2D
34///
35/// **Input type**: Geodetic coordinates
36///
37/// **Output type**: Projected coordinates
38///
39/// ## Projection String
40/// ```ini
41/// +proj=eqdc +lat_1=55 +lat_2=60
42/// ```
43///
44/// ## Parameters
45///
46/// ### Required
47/// - `+lat_1` (Latitude of the first standard parallel)
48/// - `+lat_2` (Latitude of the second standard parallel)
49///
50/// ### Optional
51/// - `+lon_0` (Central meridian)
52/// - `+ellps` (Ellipsoid name)
53/// - `+R` (Radius of the sphere)
54/// - `+x_0` (False easting)
55/// - `+y_0` (False northing)
56///
57/// ![Equidistant Conic](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/eqdc.png?raw=true)
58#[derive(Debug, Clone, PartialEq)]
59pub struct EquidistantConicProjection {
60    proj: Rc<RefCell<Proj>>,
61    store: RefCell<EqdcData>,
62}
63impl ProjectCoordinates for EquidistantConicProjection {
64    fn code(&self) -> i64 {
65        EQUIDISTANT_CONIC
66    }
67    fn name(&self) -> &'static str {
68        "Equidistant Conic"
69    }
70    fn names() -> &'static [&'static str] {
71        &["Equidistant Conic", "Equidistant_Conic", "eqdc"]
72    }
73}
74impl CoordinateStep for EquidistantConicProjection {
75    fn new(proj: Rc<RefCell<Proj>>) -> Self {
76        let mut store = EqdcData::default();
77
78        {
79            let proj = &mut proj.borrow_mut();
80
81            store.phi1 = proj
82                .params
83                .get(&LATITUDE_OF_FIRST_STANDARD_PARALLEL)
84                .unwrap_or(&ProjValue::default())
85                .f64();
86            store.phi2 = proj
87                .params
88                .get(&LATITUDE_OF_SECOND_STANDARD_PARALLEL)
89                .unwrap_or(&ProjValue::default())
90                .f64();
91
92            if fabs(store.phi1) > FRAC_PI_2 {
93                panic!("Invalid value for lat_1: |lat_1| should be <= 90°");
94            }
95
96            if fabs(store.phi2) > FRAC_PI_2 {
97                panic!("Invalid value for lat_2: |lat_2| should be <= 90°");
98            }
99            if fabs(store.phi1 + store.phi2) < EPS10 {
100                panic!("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
101            }
102
103            store.en = enfn(proj.n);
104
105            let mut sinphi = sin(store.phi1);
106            store.n = sinphi;
107            let mut cosphi = cos(store.phi1);
108            let secant = fabs(store.phi1 - store.phi2) >= EPS10;
109            store.ellips = proj.es > 0.;
110            if store.ellips {
111                let m1 = _msfn(sinphi, cosphi, proj.es);
112                let ml1 = mlfn(store.phi1, sinphi, cosphi, &store.en);
113                if secant {
114                    // secant cone
115                    sinphi = sin(store.phi2);
116                    cosphi = cos(store.phi2);
117                    let ml2 = mlfn(store.phi2, sinphi, cosphi, &store.en);
118                    if ml1 == ml2 {
119                        panic!("Eccentricity too close to 1");
120                    }
121                    store.n = (m1 - _msfn(sinphi, cosphi, proj.es)) / (ml2 - ml1);
122                    if store.n == 0. {
123                        // Not quite, but es is very close to 1...
124                        panic!("Invalid value for eccentricity");
125                    }
126                }
127                store.c = ml1 + m1 / store.n;
128                store.rho0 = store.c - mlfn(proj.phi0, sin(proj.phi0), cos(proj.phi0), &store.en);
129            } else {
130                if secant {
131                    store.n = (cosphi - cos(store.phi2)) / (store.phi2 - store.phi1);
132                }
133                if store.n == 0. {
134                    panic!("Invalid value for lat_1 and lat_2: lat_1 + lat_2 should be > 0");
135                }
136                store.c = store.phi1 + cos(store.phi1) / store.n;
137                store.rho0 = store.c - proj.phi0;
138            }
139        }
140
141        EquidistantConicProjection { proj, store: store.into() }
142    }
143    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
144        eqdc_e_forward(&mut self.store.borrow_mut(), p);
145    }
146    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
147        eqdc_e_inverse(&mut self.store.borrow_mut(), p);
148    }
149}
150
151/// Equidistant Conic Ellipsoidal forward project
152pub fn eqdc_e_forward<P: TransformCoordinates>(eqdc: &mut EqdcData, p: &mut P) {
153    eqdc.rho = eqdc.c
154        - (if eqdc.ellips { mlfn(p.phi(), sin(p.phi()), cos(p.phi()), &eqdc.en) } else { p.phi() });
155    let lam_mul_n = p.lam() * eqdc.n;
156    p.set_x(eqdc.rho * sin(lam_mul_n));
157    p.set_y(eqdc.rho0 - eqdc.rho * cos(lam_mul_n));
158}
159
160/// Equidistant Conic Ellipsoidal inverse project
161pub fn eqdc_e_inverse<P: TransformCoordinates>(eqdc: &mut EqdcData, p: &mut P) {
162    let mut x = p.x();
163    let mut y = p.y();
164    y = eqdc.rho0 - y;
165    eqdc.rho = hypot(x, y);
166
167    if eqdc.rho != 0.0 {
168        if eqdc.n < 0. {
169            eqdc.rho = -eqdc.rho;
170            x = -x;
171            y = -y;
172        }
173        p.set_phi(eqdc.c - eqdc.rho);
174        if eqdc.ellips {
175            p.set_phi(inv_mlfn(p.phi(), &eqdc.en));
176        }
177        p.set_lam(atan2(x, y) / eqdc.n);
178    } else {
179        p.set_lam(0.);
180        p.set_phi(if eqdc.n > 0. { FRAC_PI_2 } else { -FRAC_PI_2 });
181    }
182}
gistools/proj/project/eqdc.rs

gistools/proj/project/
eqdc.rs
