gistools/proj/project/
poly.rs

1use crate::proj::{
2    _msfn, CoordinateStep, POLYCONIC, Proj, ProjMethod, ProjectCoordinates, TransformCoordinates,
3    enfn, mlfn,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::cell::RefCell;
7use libm::{asin, cos, fabs, sin, sqrt, tan};
8
9/// Polyconic Variables
10#[derive(Debug, Default, Clone, PartialEq)]
11pub struct PolyData {
12    ml0: f64,
13    en: Vec<f64>,
14}
15
16const TOL: f64 = 1e-10;
17const CONV: f64 = 1e-10;
18const N_ITER: usize = 10;
19const I_ITER: usize = 20;
20const ITOL: f64 = 1e-12;
21
22/// # Polyconic (American)
23///
24/// **Classification**: Pseudoconical
25///
26/// **Available forms**: Forward and inverse, spherical and ellipsoidal
27///
28/// **Defined area**: Global
29///
30/// **Alias**: poly
31///
32/// **Domain**: 2D
33///
34/// **Input type**: Geodetic coordinates
35///
36/// **Output type**: Projected coordinates
37///
38/// ## Projection String
39/// ```ini
40/// +proj=poly
41/// ```
42///
43/// ## Required Parameters
44/// - None
45///
46/// ## Optional Parameters
47/// - `+lon_0=<value>`: Central meridian.
48/// - `+ellps=<value>`: Ellipsoid used.
49/// - `+R=<value>`: Radius of the projection sphere.
50/// - `+x_0=<value>`: False easting.
51/// - `+y_0=<value>`: False northing.
52///
53/// ![Polyconic (American)](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/poly.png?raw=true)
54#[derive(Debug, Clone, PartialEq)]
55pub struct PolyconicProjection {
56    proj: Rc<RefCell<Proj>>,
57    store: RefCell<PolyData>,
58    method: ProjMethod,
59}
60impl ProjectCoordinates for PolyconicProjection {
61    fn code(&self) -> i64 {
62        POLYCONIC
63    }
64    fn name(&self) -> &'static str {
65        "Polyconic"
66    }
67    fn names() -> &'static [&'static str] {
68        &["Polyconic", "Polyconic (American)", "poly"]
69    }
70}
71impl CoordinateStep for PolyconicProjection {
72    fn new(proj: Rc<RefCell<Proj>>) -> Self {
73        let mut store = PolyData::default();
74        let method: ProjMethod;
75        {
76            let proj = &mut proj.borrow_mut();
77
78            method = if proj.es != 0.0 {
79                store.en = enfn(proj.n);
80                store.ml0 = mlfn(proj.phi0, sin(proj.phi0), cos(proj.phi0), &store.en);
81                ProjMethod::Ellipsoidal
82            } else {
83                store.ml0 = -proj.phi0;
84                ProjMethod::Spheroidal
85            }
86        }
87
88        PolyconicProjection { proj, store: store.into(), method }
89    }
90    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
91        if self.method == ProjMethod::Spheroidal {
92            poly_s_forward(&self.store.borrow(), &self.proj.borrow(), p);
93        } else {
94            poly_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
95        }
96    }
97    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
98        if self.method == ProjMethod::Spheroidal {
99            poly_s_inverse(&self.proj.borrow(), p);
100        } else {
101            poly_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
102        }
103    }
104}
105
106/// Polyconic Ellipsoidal forward project
107pub fn poly_e_forward<P: TransformCoordinates>(poly: &PolyData, proj: &Proj, p: &mut P) {
108    if fabs(p.phi()) <= TOL {
109        p.set_x(p.lam());
110        p.set_y(-poly.ml0);
111    } else {
112        let sp = sin(p.phi());
113        let cp = cos(p.phi());
114        let ms = if fabs(cp) > TOL { _msfn(sp, cp, proj.es) / sp } else { 0. };
115        p.set_lam(p.lam() * sp);
116        p.set_x(ms * sin(p.lam()));
117        p.set_y((mlfn(p.phi(), sp, cp, &poly.en) - poly.ml0) + ms * (1. - cos(p.lam())));
118    }
119}
120
121/// Polyconic Spheroidal forward project
122pub fn poly_s_forward<P: TransformCoordinates>(poly: &PolyData, proj: &Proj, p: &mut P) {
123    if fabs(p.phi()) <= TOL {
124        p.set_x(p.lam());
125        p.set_y(poly.ml0);
126    } else {
127        let cot = 1. / tan(p.phi());
128        let e = p.lam() * sin(p.phi());
129        p.set_x(sin(e) * cot);
130        p.set_y(p.phi() - proj.phi0 + cot * (1. - cos(e)));
131    }
132}
133
134/// Polyconic Ellipsoidal inverse project
135pub fn poly_e_inverse<P: TransformCoordinates>(poly: &PolyData, proj: &Proj, p: &mut P) {
136    let x = p.x();
137    let mut y = p.y();
138    let mut phi;
139    let lam;
140    y += poly.ml0;
141    if fabs(y) <= TOL {
142        lam = x;
143        phi = 0.;
144    } else {
145        let r = y * y + x * x;
146        phi = y;
147        let mut i = I_ITER;
148        while i > 0 {
149            let sp = sin(phi);
150            let cp = cos(phi);
151            let s2ph = sp * cp;
152            if fabs(cp) < ITOL {
153                panic!("Coordinate outside projection domain");
154            }
155            let mut mlp = sqrt(1. - proj.es * sp * sp);
156            let c = sp * mlp / cp;
157            let ml = mlfn(phi, sp, cp, &poly.en);
158            let mlb = ml * ml + r;
159            mlp = proj.one_es / (mlp * mlp * mlp);
160            let d_phi = (ml + ml + c * mlb - 2. * y * (c * ml + 1.))
161                / (proj.es * s2ph * (mlb - 2. * y * ml) / c
162                    + 2. * (y - ml) * (c * mlp - 1. / s2ph)
163                    - mlp
164                    - mlp);
165            phi += d_phi;
166            if fabs(d_phi) <= ITOL {
167                break;
168            }
169            i -= 1;
170        }
171        if i == 0 {
172            panic!("Coordinate outside projection domain");
173        }
174        let c = sin(phi);
175        lam = asin(x * tan(phi) * sqrt(1. - proj.es * c * c)) / sin(phi);
176    }
177
178    p.set_lam(lam);
179    p.set_phi(phi);
180}
181
182/// Polyconic Spheroidal inverse project
183pub fn poly_s_inverse<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
184    let x = p.x();
185    let y = proj.phi0 + p.y();
186    if fabs(y) <= TOL {
187        p.set_lam(x);
188        p.set_phi(0.);
189    } else {
190        p.set_phi(y);
191        let b = x * x + y * y;
192        let mut i = N_ITER;
193        loop {
194            let tp = tan(p.phi());
195            let dphi = (y * (p.phi() * tp + 1.) - p.phi() - 0.5 * (p.phi() * p.phi() + b) * tp)
196                / ((p.phi() - y) / tp - 1.);
197            p.set_phi(p.phi() - dphi);
198            if fabs(dphi) <= CONV {
199                break;
200            }
201            i -= 1;
202            if i == 0 {
203                panic!("Coordinate outside projection domain");
204            }
205        }
206        p.set_lam(asin(x * tan(p.phi())) / sin(p.phi()));
207    }
208}
gistools/proj/project/poly.rs

gistools/proj/project/
poly.rs
