gistools/proj/project/
aea.rs

1use super::{ProjectCoordinates, TransformCoordinates};
2use crate::proj::{
3    _msfn, ALBERS_EQUAL_AREA, CoordinateStep, EPS10, LATITUDE_OF_FIRST_STANDARD_PARALLEL,
4    LATITUDE_OF_SECOND_STANDARD_PARALLEL, Proj, ProjValue, SOUTH, authalic_lat_compute_coeffs,
5    authalic_lat_inverse, authalic_lat_q,
6};
7use alloc::{rc::Rc, vec::Vec};
8use core::{cell::RefCell, f64::consts::FRAC_PI_2};
9use libm::{asin, atan2, cos, fabs, hypot, log, sin, sqrt};
10
11const TOL7: f64 = 1.0e-7;
12
13/// Albers Equal Area variables. Shared by aea and leac
14#[derive(Debug, Default, Clone, PartialEq)]
15pub struct AeaData {
16    ec: f64,
17    n: f64,
18    c: f64,
19    dd: f64,
20    n2: f64,
21    rho0: f64,
22    rho: f64,
23    phi1: f64,
24    phi2: f64,
25    ellips: bool,
26    apa: Vec<f64>,
27    qp: f64,
28}
29impl AeaData {
30    fn new(proj: Rc<RefCell<Proj>>) -> Self {
31        let proj = &proj.borrow();
32        let mut store = AeaData::default();
33
34        if let Some(lat_1) = proj.params.get(&LATITUDE_OF_FIRST_STANDARD_PARALLEL) {
35            store.phi1 = lat_1.f64();
36        }
37        if let Some(lat_2) = proj.params.get(&LATITUDE_OF_SECOND_STANDARD_PARALLEL) {
38            store.phi2 = lat_2.f64();
39        }
40
41        if fabs(store.phi1) > FRAC_PI_2 {
42            panic!("Invalid value for lat_1: |lat_1| should be <= 90°");
43        }
44        if fabs(store.phi2) > FRAC_PI_2 {
45            panic!("Invalid value for lat_2: |lat_2| should be <= 90°");
46        }
47        if fabs(store.phi1 + store.phi2) < EPS10 {
48            panic!("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
49        }
50        let mut sinphi = sin(store.phi1);
51        store.n = sinphi;
52        let mut cosphi = cos(store.phi1);
53        let secant = fabs(store.phi1 - store.phi2) >= EPS10;
54        store.ellips = proj.es > 0.;
55        if store.ellips {
56            store.apa = authalic_lat_compute_coeffs(proj.n);
57            store.qp = authalic_lat_q(1.0, proj);
58            let m1: f64 = _msfn(sinphi, cosphi, proj.es);
59            let ml1: f64 = authalic_lat_q(sinphi, proj);
60            if secant {
61                // secant cone
62                sinphi = sin(store.phi2);
63                cosphi = cos(store.phi2);
64                let m2: f64 = _msfn(sinphi, cosphi, proj.es);
65                let ml2: f64 = authalic_lat_q(sinphi, proj);
66                if ml2 == ml1 {
67                    panic!("Invalid value for lat_1 and lat_2: latitudes are too close");
68                }
69
70                store.n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
71                if store.n == 0. {
72                    // Not quite, but es is very close to 1...
73                    panic!("Invalid value for eccentricity");
74                }
75            }
76            store.ec = 1. - 0.5 * proj.one_es * log((1. - proj.e) / (1. + proj.e)) / proj.e;
77            store.c = m1 * m1 + store.n * ml1;
78            store.dd = 1. / store.n;
79            store.rho0 = store.dd * sqrt(store.c - store.n * authalic_lat_q(sin(proj.phi0), proj));
80        } else {
81            if secant {
82                store.n = 0.5 * (store.n + sin(store.phi2));
83            }
84            store.n2 = store.n + store.n;
85            store.c = cosphi * cosphi + store.n2 * sinphi;
86            store.dd = 1. / store.n;
87            store.rho0 = store.dd * sqrt(store.c - store.n2 * sin(proj.phi0));
88        }
89
90        store
91    }
92}
93
94/// # Albers Conic Equal Area Projection
95///
96/// **Classification**: Conic
97///
98/// **Available forms**: Forward and inverse, spherical and ellipsoidal
99///
100/// **Defined area**: Global
101///
102/// **Alias**: `aea`
103///
104/// **Domain**: 2D
105///
106/// **Input type**: Geodetic coordinates
107///
108/// **Output type**: Projected coordinates
109///
110/// ## Projection String
111/// ```sh
112/// +proj=aea +lat_1=29.5 +lat_2=42.5
113/// ```
114///
115/// ## Required Parameters
116/// - `lat1`
117/// - `lat2`
118///
119/// ## Optional Parameters
120/// - `lon0`: Longitude of the central meridian.
121/// - `ellps`: Name of the reference ellipsoid.
122/// - `R`: Radius of the sphere if `ellps` is not specified.
123/// - `x0`: False easting (coordinate offset in the x-direction).
124/// - `y0`: False northing (coordinate offset in the y-direction).
125///
126/// ## References
127/// - <https://en.wikipedia.org/wiki/Albers_projection>
128///
129/// ![Albers Conic Equal Area Projection](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/aea.png?raw=true)
130#[derive(Debug, Clone, PartialEq)]
131pub struct AlbersConicEqualAreaProjection {
132    proj: Rc<RefCell<Proj>>,
133    store: RefCell<AeaData>,
134}
135impl ProjectCoordinates for AlbersConicEqualAreaProjection {
136    fn code(&self) -> i64 {
137        ALBERS_EQUAL_AREA
138    }
139    fn name(&self) -> &'static str {
140        "Albers Conic Equal Area"
141    }
142    fn names() -> &'static [&'static str] {
143        &["Albers Conic Equal Area", "Albers_Conic_Equal_Area", "Albers", "aea", "9822"]
144    }
145}
146impl CoordinateStep for AlbersConicEqualAreaProjection {
147    fn new(proj: Rc<RefCell<Proj>>) -> Self {
148        let store = AeaData::new(proj.clone());
149        AlbersConicEqualAreaProjection { proj, store: store.into() }
150    }
151    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
152        aea_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
153    }
154    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
155        aea_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
156    }
157}
158
159/// # Lambert Equal Area Conic Projection
160///
161/// **Classification**: Conical
162///
163/// **Available forms**: Forward and inverse, spherical and ellipsoidal
164///
165/// **Defined area**: Global
166///
167/// **Alias**: `leac`
168///
169/// **Domain**: 2D
170///
171/// **Input type**: Geodetic coordinates
172///
173/// **Output type**: Projected coordinates
174///
175/// ## Projection String
176/// ```ini
177/// +proj=leac
178/// ```
179///
180/// ## Parameters
181///
182/// **Note**: All parameters are optional for the Lambert Equal Area Conic projection.
183///
184/// ## Required Parameters
185/// - `lat1`: Latitude of the first standard parallel.
186/// - `+south`: Sets the second standard parallel to 90°S. When the flag is off, the second standard parallel is set to 90°N.
187///
188/// ## Optional Parameters
189/// - `lon0`: Longitude of the central meridian.
190/// - `ellps`: Name of the reference ellipsoid.
191/// - `R`: Radius of the sphere if `ellps` is not specified.
192/// - `x0`: False easting (coordinate offset in the x-direction).
193/// - `y0`: False northing (coordinate offset in the y-direction).
194///
195/// ## References
196/// - <https://en.wikipedia.org/wiki/Lambert_conformal_conic_projection> (Note: While the name is similar, this link describes the conformal variant. A specific link for the equal-area conic might be needed)
197///
198/// ![Lambert Equal Area Conic Projection](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/leac.png?raw=true)
199#[derive(Debug, Clone, PartialEq)]
200pub struct LambertEqualAreaConicProjection {
201    proj: Rc<RefCell<Proj>>,
202    store: RefCell<AeaData>,
203}
204impl ProjectCoordinates for LambertEqualAreaConicProjection {
205    fn code(&self) -> i64 {
206        -1
207    }
208    fn name(&self) -> &'static str {
209        "Lambert Equal Area Conic"
210    }
211    fn names() -> &'static [&'static str] {
212        &["Lambert Equal Area Conic", "leac"]
213    }
214}
215impl CoordinateStep for LambertEqualAreaConicProjection {
216    fn new(proj: Rc<RefCell<Proj>>) -> Self {
217        let mut store = AeaData::new(proj.clone());
218        store.phi2 = proj.borrow().params.get(&LATITUDE_OF_FIRST_STANDARD_PARALLEL).unwrap().f64();
219        store.phi1 = if proj.borrow().params.get(&SOUTH).unwrap_or(&ProjValue::default()).bool() {
220            -FRAC_PI_2
221        } else {
222            FRAC_PI_2
223        };
224        LambertEqualAreaConicProjection { proj, store: store.into() }
225    }
226    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
227        aea_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
228    }
229    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
230        aea_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
231    }
232}
233
234/// Ellipsoid/spheroid, forward
235pub fn aea_e_forward<P: TransformCoordinates>(aea: &mut AeaData, proj: &Proj, p: &mut P) {
236    let phi = p.phi();
237    let mut lam = p.lam();
238    aea.rho =
239        aea.c - if aea.ellips { aea.n * authalic_lat_q(sin(phi), proj) } else { aea.n2 * sin(phi) };
240
241    if aea.rho < 0. {
242        // ERROR: PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN
243        return;
244    }
245    aea.rho = aea.dd * sqrt(aea.rho);
246    lam *= aea.n;
247    p.set_x(aea.rho * sin(lam));
248    p.set_y(aea.rho0 - aea.rho * cos(lam));
249}
250
251/// Ellipsoid/spheroid, inverse
252pub fn aea_e_inverse<P: TransformCoordinates>(aea: &mut AeaData, proj: &Proj, p: &mut P) {
253    let mut x = p.x();
254    let mut y = p.y();
255    y = aea.rho0 - y;
256    aea.rho = hypot(x, y);
257    if aea.rho != 0.0 {
258        if aea.n < 0. {
259            aea.rho = -aea.rho;
260            x = -x;
261            y = -y;
262        }
263        p.set_phi(aea.rho / aea.dd);
264        if aea.ellips {
265            let qs = (aea.c - p.phi() * p.phi()) / aea.n;
266            if fabs(aea.ec - fabs(qs)) > TOL7 {
267                if fabs(qs) > 2. {
268                    panic!("Coordinate outside projection domain");
269                }
270                p.set_phi(authalic_lat_inverse(asin(qs / aea.qp), &aea.apa, proj, aea.qp));
271                if p.phi() == f64::INFINITY {
272                    panic!("Coordinate outside projection domain");
273                }
274            } else {
275                p.set_phi(if qs < 0. { -FRAC_PI_2 } else { FRAC_PI_2 });
276            }
277        } else {
278            let qs_div_2 = (aea.c - p.phi() * p.phi()) / aea.n2;
279            if fabs(qs_div_2) <= 1. {
280                p.set_phi(asin(qs_div_2));
281            } else {
282                p.set_phi(if qs_div_2 < 0. { -FRAC_PI_2 } else { FRAC_PI_2 });
283            }
284        }
285        p.set_lam(atan2(x, y) / aea.n);
286    } else {
287        p.set_lam(0.);
288        p.set_phi(if aea.n > 0. { FRAC_PI_2 } else { -FRAC_PI_2 });
289    }
290}
gistools/proj/project/aea.rs

gistools/proj/project/
aea.rs
