gistools/proj/project/
cea.rs

1use crate::proj::{
2    CoordinateStep, LATITUDE_STD_PARALLEL, Proj, ProjMethod, ProjectCoordinates,
3    TransformCoordinates, authalic_lat_compute_coeffs, authalic_lat_inverse, authalic_lat_q,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::{cell::RefCell, f64::consts::FRAC_PI_2};
7use libm::{asin, cos, fabs, sin, sqrt};
8
9/// Equal Area Cylindrical Variables
10#[derive(Debug, Default, Clone, PartialEq)]
11pub struct CeaData {
12    qp: f64,
13    apa: Vec<f64>,
14}
15
16const EPS: f64 = 1e-10;
17
18/// # Equal Area Cylindrical Projection
19///
20/// **Classification**: Cylindrical
21///
22/// **Available forms**: Forward and Inverse, spherical and ellipsoidal
23///
24/// **Defined area**: Global
25///
26/// **Alias**: `cea`
27///
28/// **Domain**: 2D
29///
30/// **Input type**: Geodetic coordinates
31///
32/// **Output type**: Projected coordinates
33///
34/// ## Projection String
35/// ```ini
36/// +proj=cea
37/// ```
38///
39/// ## Named Specializations
40///
41/// The Equal Area Cylindrical projection is sometimes known under other names when instantiated with particular values of the `lat_ts` parameter:
42///
43/// - **Lambert cylindrical equal-area**: `lat_ts = 0`
44/// - **Behrmann**: `lat_ts = 30`
45/// - **Gall-Peters**: `lat_ts = 45`
46///
47/// ## Optional Parameters
48/// - `lat_ts`: Latitude of true scale
49/// - `lon_0`: Longitude of origin
50/// - `ellps`: Ellipsoid
51/// - `R`: Radius of the sphere
52/// - `k_0`: Scaling factor
53/// - `x_0`: False easting
54/// - `y_0`: False northing
55///
56/// **Note**: `lat_ts` and `k_0` are mutually exclusive. If `lat_ts` is specified, it is equivalent to setting `k_0` to:
57///
58/// $$k_0 = \cos(lat_{ts}) / \sqrt{1 - e^2 \sin^2(lat_{ts})}$$
59///
60/// reference:
61/// "Cartographic Projection Procedures for the UNIX Environment-
62/// A User's Manual" by Gerald I. Evenden,
63/// USGS Open File Report 90-284and Release 4 Interim Reports (2003)
64///
65/// ![Equal Area Cylindrical Projection](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/cae.png?raw=true)
66#[derive(Debug, Clone, PartialEq)]
67pub struct EqualAreaCylindricalProjection {
68    proj: Rc<RefCell<Proj>>,
69    store: RefCell<CeaData>,
70    method: ProjMethod,
71}
72impl ProjectCoordinates for EqualAreaCylindricalProjection {
73    fn code(&self) -> i64 {
74        -1
75    }
76    fn name(&self) -> &'static str {
77        "Equal Area Cylindrical"
78    }
79    fn names() -> &'static [&'static str] {
80        &["Equal Area Cylindrical", "cea"]
81    }
82}
83impl CoordinateStep for EqualAreaCylindricalProjection {
84    fn new(proj: Rc<RefCell<Proj>>) -> Self {
85        let mut store = CeaData::default();
86        let method: ProjMethod;
87        {
88            let proj = &mut proj.borrow_mut();
89            let mut t = 0.0;
90
91            if let Some(lat_ts) = proj.params.get(&LATITUDE_STD_PARALLEL) {
92                t = lat_ts.f64();
93                proj.k0 = cos(t);
94                if proj.k0 < 0.0 {
95                    panic!("Invalid value for lat_ts: |lat_ts| should be <= 90°");
96                }
97            }
98            method = if proj.es != 0.0 {
99                t = sin(t);
100                proj.k0 /= sqrt(1. - proj.es * t * t);
101                proj.e = sqrt(proj.es);
102                store.apa = authalic_lat_compute_coeffs(proj.n);
103
104                store.qp = authalic_lat_q(1.0, proj);
105                ProjMethod::Ellipsoidal
106            } else {
107                ProjMethod::Spheroidal
108            };
109        }
110        EqualAreaCylindricalProjection { proj, store: store.into(), method }
111    }
112    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
113        if self.method == ProjMethod::Spheroidal {
114            cea_s_forward(&self.proj.borrow(), p);
115        } else {
116            cea_e_forward(&self.proj.borrow(), p);
117        }
118    }
119    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
120        if self.method == ProjMethod::Spheroidal {
121            cea_s_inverse(&self.proj.borrow(), p);
122        } else {
123            cea_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
124        }
125    }
126}
127
128/// Equal Area Cylindrical Ellipsoidal forward project
129pub fn cea_e_forward<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
130    p.set_x(proj.k0 * p.lam());
131    p.set_y(0.5 * authalic_lat_q(sin(p.phi()), proj) / proj.k0);
132}
133
134/// Equal Area Spheroidal forward project
135pub fn cea_s_forward<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
136    p.set_x(proj.k0 * p.lam());
137    p.set_y(sin(p.phi()) / proj.k0);
138}
139
140/// Equal Area Ellipsoidal inverse project
141pub fn cea_e_inverse<P: TransformCoordinates>(cea: &CeaData, proj: &Proj, p: &mut P) {
142    p.set_phi(authalic_lat_inverse(asin(2. * p.y() * proj.k0 / cea.qp), &cea.apa, proj, cea.qp));
143    p.set_lam(p.x() / proj.k0);
144}
145
146/// Equal Area Spheroidal inverse project
147pub fn cea_s_inverse<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
148    let y = p.y() * proj.k0;
149    let x = p.x();
150    let t = fabs(y);
151    if t - EPS <= 1. {
152        if t >= 1. {
153            p.set_phi(if y < 0. { -FRAC_PI_2 } else { FRAC_PI_2 });
154        } else {
155            p.set_phi(asin(y));
156        }
157        p.set_lam(x / proj.k0);
158    } else {
159        panic!("Coordinate outside projection domain");
160    }
161}
gistools/proj/project/cea.rs

gistools/proj/project/
cea.rs
