gistools/proj/project/
eqearth.rs

1use crate::proj::{
2    CoordinateStep, EQUAL_EARTH, Proj, ProjectCoordinates, TransformCoordinates,
3    authalic_lat_compute_coeffs, authalic_lat_inverse, authalic_lat_q,
4};
5use alloc::{rc::Rc, vec::Vec};
6use core::cell::RefCell;
7use libm::{asin, cos, fabs, sin, sqrt};
8
9// Equal Earth is a projection inspired by the Robinson projection, but unlike
10// the Robinson projection retains the relative size of areas. The projection
11// was designed in 2018 by Bojan Savric, Tom Patterson and Bernhard Jenny.
12//
13// Publication:
14// Bojan Savric, Tom Patterson & Bernhard Jenny (2018). The Equal Earth map
15// projection, International Journal of Geographical Information Science,
16// DOI: 10.1080/13658816.2018.1504949
17//
18// Port to PROJ by Juernjakob Dugge, 16 August 2018
19// Added ellipsoidal equations by Bojan Savric, 22 August 2018
20
21// A1..A4, polynomial coefficients
22const A1: f64 = 1.340264;
23const A2: f64 = -0.081106;
24const A3: f64 = 0.000893;
25const A4: f64 = 0.003796;
26const M: f64 = 0.8660254037844386; // sqrt(3.0) / 2.0;
27
28// 90° latitude on a sphere with radius 1
29const MAX_Y: f64 = 1.3173627591574;
30const EPS: f64 = 1e-11;
31const MAX_ITER: usize = 12;
32
33/// Equal Earth variables
34#[derive(Debug, Default, Clone, PartialEq)]
35pub struct EqEarth {
36    qp: f64,
37    rqda: f64,
38    apa: Vec<f64>,
39}
40
41/// # Equal Earth
42///
43/// **Classification**: Pseudo cylindrical
44///
45/// **Available forms**: Forward and inverse, spherical and ellipsoidal projection
46///
47/// **Defined area**: Global
48///
49/// **Alias**: eqearth
50///
51/// **Domain**: 2D
52///
53/// **Input type**: Geodetic coordinates
54///
55/// **Output type**: Projected coordinates
56///
57/// ## Projection String
58/// ```ini
59/// +proj=eqearth
60/// ```
61///
62/// ## Usage
63///
64/// The Equal Earth projection is designed for world maps and retains the relative size of areas. It was inspired by the Robinson projection.
65///
66/// Example:
67/// ```bash
68/// $ echo 122 47 | proj +proj=eqearth +R=1
69/// 1.55    0.89
70/// ```
71///
72/// ## Parameters
73///
74/// **Note**: All parameters for this projection are optional.
75///
76/// ### Optional
77/// - `+lon_0` (Central meridian)
78/// - `+ellps` (Ellipsoid name)
79/// - `+R` (Radius of the sphere)
80/// - `+x_0` (False easting)
81/// - `+y_0` (False northing)
82///
83/// ## Further Reading
84/// - [The Equal Earth map projection](https://www.researchgate.net/profile/Bojan_Savric2/publication/326879978_The_Equal_Earth_map_projection/links/5b69d0ae299bf14c6d951b77/The-Equal-Earth-map-projection.pdf) by Bojan Savric, Tom Patterson & Bernhard Jenny (2018)
85///
86/// ![Equal Earth](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/eqearth.png?raw=true)
87#[derive(Debug, Clone, PartialEq)]
88pub struct EqualEarthProjection {
89    proj: Rc<RefCell<Proj>>,
90    store: RefCell<EqEarth>,
91}
92impl ProjectCoordinates for EqualEarthProjection {
93    fn code(&self) -> i64 {
94        EQUAL_EARTH
95    }
96    fn name(&self) -> &'static str {
97        "Equal Earth"
98    }
99    fn names() -> &'static [&'static str] {
100        &["Equal Earth", "EqualEarth", "eqearth"]
101    }
102}
103impl CoordinateStep for EqualEarthProjection {
104    fn new(proj: Rc<RefCell<Proj>>) -> Self {
105        let mut store = EqEarth { rqda: 1.0, ..Default::default() };
106        {
107            let proj = &mut proj.borrow_mut();
108            // Ellipsoidal case
109            if proj.es != 0.0 {
110                store.apa = authalic_lat_compute_coeffs(proj.n); // For auth_lat().
111                store.qp = authalic_lat_q(1.0, proj); // For auth_lat().
112                store.rqda = sqrt(0.5 * store.qp); // Authalic radius divided by major axis
113            }
114        }
115        EqualEarthProjection { proj, store: store.into() }
116    }
117    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
118        eqearth_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
119    }
120    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
121        eqearth_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
122    }
123}
124
125/// Equal Earth Ellipsoidal/spheroidal forward project
126pub fn eqearth_e_forward<P: TransformCoordinates>(eq_earth: &mut EqEarth, proj: &Proj, p: &mut P) {
127    // Spheroidal case, using sine latitude
128    let mut sbeta = sin(p.phi());
129
130    // In the ellipsoidal case, we convert sbeta to sine of authalic latitude
131    if proj.es != 0.0 {
132        sbeta = authalic_lat_q(sbeta, proj) / eq_earth.qp;
133
134        // Rounding error.
135        if fabs(sbeta) > 1. {
136            sbeta = if sbeta > 0. { 1. } else { -1. };
137        }
138    }
139
140    // Equal Earth projection
141    let psi = asin(M * sbeta);
142    let psi2 = psi * psi;
143    let psi6 = psi2 * psi2 * psi2;
144
145    let mut x =
146        p.lam() * cos(psi) / (M * (A1 + 3. * A2 * psi2 + psi6 * (7. * A3 + 9. * A4 * psi2)));
147    let mut y = psi * (A1 + A2 * psi2 + psi6 * (A3 + A4 * psi2));
148
149    // Adjusting x and y for authalic radius
150    x *= eq_earth.rqda;
151    y *= eq_earth.rqda;
152
153    p.set_x(x);
154    p.set_y(y);
155}
156
157/// Equal Earth Ellipsoidal/spheroidal inverse project
158pub fn eqearth_e_inverse<P: TransformCoordinates>(eq_earth: &mut EqEarth, proj: &Proj, p: &mut P) {
159    // Adjusting x and y for authalic radius
160    let x = p.x() / eq_earth.rqda;
161    let mut y = p.y() / eq_earth.rqda;
162
163    // Make sure y is inside valid range
164    y = y.clamp(-MAX_Y, MAX_Y);
165
166    let mut yc = y;
167
168    // Newton-Raphson
169    let mut i = MAX_ITER;
170    while i > 0 {
171        let y2 = yc * yc;
172        let y6 = y2 * y2 * y2;
173
174        let f = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - y;
175        let fder = A1 + 3. * A2 * y2 + y6 * (7. * A3 + 9. * A4 * y2);
176
177        let tol = f / fder;
178        yc -= tol;
179
180        if fabs(tol) < EPS {
181            break;
182        }
183        i -= 1;
184    }
185
186    if i == 0 {
187        panic!("Coordinate outside projection domain");
188    }
189
190    // Longitude
191    let y2 = yc * yc;
192    let y6 = y2 * y2 * y2;
193
194    p.set_lam(M * x * (A1 + 3. * A2 * y2 + y6 * (7. * A3 + 9. * A4 * y2)) / cos(yc));
195
196    // Latitude (for spheroidal case, this is latitude
197    p.set_phi(asin(sin(yc) / M));
198
199    // Ellipsoidal case, converting auth. latitude
200    if proj.es != 0.0 {
201        p.set_phi(authalic_lat_inverse(p.phi(), &eq_earth.apa, proj, eq_earth.qp));
202    }
203}
gistools/proj/project/eqearth.rs

gistools/proj/project/
eqearth.rs
