gistools/proj/project/
bonne.rs

1use crate::proj::{
2    BONNE, CoordinateStep, EPS10, LATITUDE_OF_FIRST_STANDARD_PARALLEL, Proj, ProjMethod,
3    ProjectCoordinates, TransformCoordinates, enfn, inv_mlfn, mlfn,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::{cell::RefCell, f64::consts::FRAC_PI_2};
7use libm::{atan2, copysign, cos, fabs, hypot, sin, sqrt, tan};
8
9/// Bonne Variables
10#[derive(Debug, Default, Clone, PartialEq)]
11pub struct BonneData {
12    phi1: f64,
13    cphi1: f64,
14    am1: f64,
15    m1: f64,
16    en: Vec<f64>,
17}
18
19/// # Bonne (Werner lat_1=90) Projection
20///
21/// **Classification**: Miscellaneous
22///
23/// **Available forms**: Forward and inverse, spherical and ellipsoidal
24///
25/// **Defined area**: Global
26///
27/// **Alias**: `bonne`
28///
29/// **Domain**: 2D
30///
31/// **Input type**: Geodetic coordinates
32///
33/// **Output type**: Projected coordinates
34///
35/// ## Projection String
36/// ```ini
37/// +proj=bonne +lat_1=10
38/// ```
39///
40/// ## Required Parameters
41/// - `lat1`: Latitude of first standard parallel
42///
43/// ## Optional Parameters
44/// - `lon0`: Longitude of origin
45/// - `ellps`: Ellipsoid name
46/// - `R`: Radius of sphere
47/// - `x0`: False easting
48/// - `y0`: False northing
49///
50/// ![Bonne (Werner lat_1=90) Projection](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/bonne.png?raw=true)
51#[derive(Debug, Clone, PartialEq)]
52pub struct BonneProjection {
53    proj: Rc<RefCell<Proj>>,
54    store: RefCell<BonneData>,
55    method: ProjMethod,
56}
57impl ProjectCoordinates for BonneProjection {
58    fn code(&self) -> i64 {
59        BONNE
60    }
61    fn name(&self) -> &'static str {
62        "Bonne"
63    }
64    fn names() -> &'static [&'static str] {
65        &["Bonne (Werner lat_1=90)", "bonne_werner", "Bonne", "bonne"]
66    }
67}
68impl CoordinateStep for BonneProjection {
69    fn new(proj: Rc<RefCell<Proj>>) -> Self {
70        let mut store = BonneData::default();
71        let method: ProjMethod;
72        {
73            let proj = &mut proj.borrow_mut();
74            store.phi1 = proj.params.get(&LATITUDE_OF_FIRST_STANDARD_PARALLEL).unwrap().f64();
75            if fabs(store.phi1) < EPS10 {
76                panic!("Invalid value for lat_1: |lat_1| should be > 0");
77            }
78
79            method = if proj.es != 0.0 {
80                store.en = enfn(proj.n);
81                store.am1 = sin(store.phi1);
82                let c = cos(store.phi1);
83                store.m1 = mlfn(store.phi1, store.am1, c, &store.en);
84                store.am1 = c / (sqrt(1. - proj.es * store.am1 * store.am1) * store.am1);
85                ProjMethod::Ellipsoidal
86            } else {
87                if fabs(store.phi1) + EPS10 >= FRAC_PI_2 {
88                    store.cphi1 = 0.;
89                } else {
90                    store.cphi1 = 1. / tan(store.phi1);
91                }
92                ProjMethod::Spheroidal
93            };
94        }
95        BonneProjection { proj, store: store.into(), method }
96    }
97    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
98        if self.method == ProjMethod::Ellipsoidal {
99            bonne_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
100        } else {
101            bonne_s_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
102        }
103    }
104    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
105        if self.method == ProjMethod::Ellipsoidal {
106            bonne_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
107        } else {
108            bonne_s_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
109        }
110    }
111}
112
113/// Bonne Ellipsoidal Forward
114pub fn bonne_e_forward<P: TransformCoordinates>(bonne: &mut BonneData, proj: &Proj, p: &mut P) {
115    let mut e = sin(p.phi());
116    let c = cos(p.phi());
117    let rh = bonne.am1 + bonne.m1 - mlfn(p.phi(), e, c, &bonne.en);
118    if fabs(rh) > EPS10 {
119        e = c * p.lam() / (rh * sqrt(1. - proj.es * e * e));
120        p.set_x(rh * sin(e));
121        p.set_y(bonne.am1 - rh * cos(e));
122    } else {
123        p.set_x(0.);
124        p.set_y(0.);
125    }
126}
127
128/// Bonne Spheroidal Forward
129pub fn bonne_s_forward<P: TransformCoordinates>(bonne: &mut BonneData, _proj: &Proj, p: &mut P) {
130    let rh = bonne.cphi1 + bonne.phi1 - p.phi();
131    if fabs(rh) > EPS10 {
132        let e = p.lam() * cos(p.phi()) / rh;
133        p.set_x(rh * sin(e));
134        p.set_y(bonne.cphi1 - rh * cos(e));
135    } else {
136        p.set_x(0.);
137        p.set_y(0.);
138    }
139}
140
141/// Bonne Spheroidal Inverse
142pub fn bonne_s_inverse<P: TransformCoordinates>(bonne: &mut BonneData, _proj: &Proj, p: &mut P) {
143    p.set_y(bonne.cphi1 - p.y());
144    let rh = copysign(hypot(p.x(), p.y()), bonne.phi1);
145    let phi = bonne.cphi1 + bonne.phi1 - rh;
146    let lam: f64;
147    let abs_phi = fabs(phi);
148    if abs_phi > FRAC_PI_2 {
149        panic!("Coordinate outside projection domain");
150    }
151    if FRAC_PI_2 - abs_phi <= EPS10 {
152        lam = 0.;
153    } else {
154        let lm = rh / cos(phi);
155        if bonne.phi1 > 0. {
156            lam = lm * atan2(p.x(), p.y());
157        } else {
158            lam = lm * atan2(-p.x(), -p.y());
159        }
160    }
161
162    p.set_phi(phi);
163    p.set_lam(lam);
164}
165
166/// Bonne Ellipsoidal Inverse
167pub fn bonne_e_inverse<P: TransformCoordinates>(bonne: &mut BonneData, proj: &Proj, p: &mut P) {
168    p.set_y(bonne.am1 - p.y());
169    let rh = copysign(hypot(p.x(), p.y()), bonne.phi1);
170    let phi = inv_mlfn(bonne.am1 + bonne.m1 - rh, &bonne.en);
171    let lam: f64;
172    let abs_phi = fabs(phi);
173    if abs_phi < FRAC_PI_2 {
174        let sinphi = sin(phi);
175        let lm = rh * sqrt(1. - proj.es * sinphi * sinphi) / cos(phi);
176        if bonne.phi1 > 0. {
177            lam = lm * atan2(p.x(), p.y());
178        } else {
179            lam = lm * atan2(-p.x(), -p.y());
180        }
181    } else if abs_phi - FRAC_PI_2 <= EPS10 {
182        lam = 0.;
183    } else {
184        panic!("Coordinates are outside the projection domain");
185    }
186    p.set_phi(phi);
187    p.set_lam(lam);
188}
gistools/proj/project/bonne.rs

gistools/proj/project/
bonne.rs
