gistools/proj/project/
robin.rs

1use crate::proj::{CoordinateStep, Proj, ProjectCoordinates, TransformCoordinates};
2use alloc::rc::Rc;
3use core::{
4    cell::RefCell,
5    f64::consts::{FRAC_PI_2, PI},
6};
7use libm::{fabs, floor, round};
8
9fn v(c: &RobinCoefs, z: f64) -> f64 {
10    c.c0 + z * (c.c1 + z * (c.c2 + z * c.c3))
11}
12fn dv(c: &RobinCoefs, z: f64) -> f64 {
13    c.c1 + 2.0 * z * c.c2 + z * z * 3.0 * c.c3
14}
15
16// note: following terms based upon 5 deg. intervals in degrees.
17//
18// Some background on these coefficients is available at:
19//
20// http://article.gmane.org/gmane.comp.gis.proj-4.devel/6039
21// http://trac.osgeo.org/proj/ticket/113
22
23/// Robin coefficients
24#[derive(Debug, Default, Clone, PartialEq)]
25pub struct RobinCoefs {
26    /// First coefficient
27    pub c0: f64,
28    /// Second coefficient
29    pub c1: f64,
30    /// Third coefficient
31    pub c2: f64,
32    /// Fourth coefficient
33    pub c3: f64,
34}
35const X_COEFS: [RobinCoefs; 19] = [
36    RobinCoefs { c0: 1.0, c1: 2.2199e-17, c2: -7.15515e-05, c3: 3.1103e-06 },
37    RobinCoefs { c0: 0.9986, c1: -0.000482243, c2: -2.4897e-05, c3: -1.3309e-06 },
38    RobinCoefs { c0: 0.9954, c1: -0.00083103, c2: -4.48605e-05, c3: -9.86701e-07 },
39    RobinCoefs { c0: 0.99, c1: -0.00135364, c2: -5.9661e-05, c3: 3.6777e-06 },
40    RobinCoefs { c0: 0.9822, c1: -0.00167442, c2: -4.49547e-06, c3: -5.72411e-06 },
41    RobinCoefs { c0: 0.973, c1: -0.00214868, c2: -9.03571e-05, c3: 1.8736e-08 },
42    RobinCoefs { c0: 0.96, c1: -0.00305085, c2: -9.00761e-05, c3: 1.64917e-06 },
43    RobinCoefs { c0: 0.9427, c1: -0.00382792, c2: -6.53386e-05, c3: -2.6154e-06 },
44    RobinCoefs { c0: 0.9216, c1: -0.00467746, c2: -0.00010457, c3: 4.81243e-06 },
45    RobinCoefs { c0: 0.8962, c1: -0.00536223, c2: -3.23831e-05, c3: -5.43432e-06 },
46    RobinCoefs { c0: 0.8679, c1: -0.00609363, c2: -0.000113898, c3: 3.32484e-06 },
47    RobinCoefs { c0: 0.835, c1: -0.00698325, c2: -6.40253e-05, c3: 9.34959e-07 },
48    RobinCoefs { c0: 0.7986, c1: -0.00755338, c2: -5.00009e-05, c3: 9.35324e-07 },
49    RobinCoefs { c0: 0.7597, c1: -0.00798324, c2: -3.5971e-05, c3: -2.27626e-06 },
50    RobinCoefs { c0: 0.7186, c1: -0.00851367, c2: -7.01149e-05, c3: -8.6303e-06 },
51    RobinCoefs { c0: 0.6732, c1: -0.00986209, c2: -0.000199569, c3: 1.91974e-05 },
52    RobinCoefs { c0: 0.6213, c1: -0.010418, c2: 8.83923e-05, c3: 6.24051e-06 },
53    RobinCoefs { c0: 0.5722, c1: -0.00906601, c2: 0.000182, c3: 6.24051e-06 },
54    RobinCoefs { c0: 0.5322, c1: -0.00677797, c2: 0.000275608, c3: 6.24051e-06 },
55];
56
57const Y_COEFS: [RobinCoefs; 19] = [
58    RobinCoefs { c0: -5.20417e-18, c1: 0.0124, c2: 1.21431e-18, c3: -8.45284e-11 },
59    RobinCoefs { c0: 0.062, c1: 0.0124, c2: -1.26793e-09, c3: 4.22642e-10 },
60    RobinCoefs { c0: 0.124, c1: 0.0124, c2: 5.07171e-09, c3: -1.60604e-09 },
61    RobinCoefs { c0: 0.186, c1: 0.0123999, c2: -1.90189e-08, c3: 6.00152e-09 },
62    RobinCoefs { c0: 0.248, c1: 0.0124002, c2: 7.10039e-08, c3: -2.24e-08 },
63    RobinCoefs { c0: 0.31, c1: 0.0123992, c2: -2.64997e-07, c3: 8.35986e-08 },
64    RobinCoefs { c0: 0.372, c1: 0.0124029, c2: 9.88983e-07, c3: -3.11994e-07 },
65    RobinCoefs { c0: 0.434, c1: 0.0123893, c2: -3.69093e-06, c3: -4.35621e-07 },
66    RobinCoefs { c0: 0.4958, c1: 0.0123198, c2: -1.02252e-05, c3: -3.45523e-07 },
67    RobinCoefs { c0: 0.5571, c1: 0.0121916, c2: -1.54081e-05, c3: -5.82288e-07 },
68    RobinCoefs { c0: 0.6176, c1: 0.0119938, c2: -2.41424e-05, c3: -5.25327e-07 },
69    RobinCoefs { c0: 0.6769, c1: 0.011713, c2: -3.20223e-05, c3: -5.16405e-07 },
70    RobinCoefs { c0: 0.7346, c1: 0.0113541, c2: -3.97684e-05, c3: -6.09052e-07 },
71    RobinCoefs { c0: 0.7903, c1: 0.0109107, c2: -4.89042e-05, c3: -1.04739e-06 },
72    RobinCoefs { c0: 0.8435, c1: 0.0103431, c2: -6.4615e-05, c3: -1.40374e-09 },
73    RobinCoefs { c0: 0.8936, c1: 0.00969686, c2: -6.4636e-05, c3: -8.547e-06 },
74    RobinCoefs { c0: 0.9394, c1: 0.00840947, c2: -0.000192841, c3: -4.2106e-06 },
75    RobinCoefs { c0: 0.9761, c1: 0.00616527, c2: -0.000256, c3: -4.2106e-06 },
76    RobinCoefs { c0: 1.0, c1: 0.00328947, c2: -0.000319159, c3: -4.2106e-06 },
77];
78
79const FXC: f64 = 0.8487;
80const FYC: f64 = 1.3523;
81const C1: f64 = 11.459_155_902_616_464;
82const RC1: f64 = 0.087_266_462_599_716_47;
83const NODES: usize = 18;
84const ONEEPS: f64 = 1.000001;
85const EPS: f64 = 1e-10;
86// Not sure at all of the appropriate number for MAX_ITER...
87const MAX_ITER: usize = 100;
88
89/// # Robinson
90///
91/// **Classification**: Pseudocylindrical
92///
93/// **Available forms**: Forward and inverse, spherical projection
94///
95/// **Defined area**: Global
96///
97/// **Alias**: robin
98///
99/// **Domain**: 2D
100///
101/// **Input type**: Geodetic coordinates
102///
103/// **Output type**: Projected coordinates
104///
105/// ## Projection String
106/// ```ini
107/// +proj=robin
108/// ```
109///
110/// ## Required Parameters
111/// - None
112///
113/// ## Optional Parameters
114/// - `+lon_0=<value>`: Central meridian.
115/// - `+R=<value>`: Radius of the projection sphere.
116/// - `+x_0=<value>`: False easting.
117/// - `+y_0=<value>`: False northing.
118///
119/// ![Robinson](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/robin.png?raw=true)
120#[derive(Debug, Clone, PartialEq)]
121pub struct RobinsonProjection {
122    proj: Rc<RefCell<Proj>>,
123}
124impl ProjectCoordinates for RobinsonProjection {
125    fn code(&self) -> i64 {
126        -1
127    }
128    fn name(&self) -> &'static str {
129        "Robinson"
130    }
131    fn names() -> &'static [&'static str] {
132        &["Robinson", "robin"]
133    }
134}
135impl CoordinateStep for RobinsonProjection {
136    fn new(proj: Rc<RefCell<Proj>>) -> Self {
137        proj.borrow_mut().es = 0.;
138        RobinsonProjection { proj }
139    }
140    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
141        robin_s_forward(p);
142    }
143    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
144        robin_s_inverse(p);
145    }
146}
147
148/// Equal Earth Spheroidal forward project
149pub fn robin_s_forward<P: TransformCoordinates>(p: &mut P) {
150    let mut dphi = fabs(p.phi());
151    let mut i = if f64::is_nan(p.phi()) { -1. } else { round(floor(dphi * C1 + 1e-15)) };
152    if i < 0. {
153        panic!("Coordinate outside projection domain");
154    }
155    if i >= NODES as f64 {
156        i = NODES as f64;
157    }
158    dphi = (dphi - RC1 * i).to_degrees();
159    let x = v(&X_COEFS[i as usize], dphi) * FXC * p.lam();
160    let mut y = v(&Y_COEFS[i as usize], dphi) * FYC;
161    if p.phi() < 0. {
162        y = -y;
163    }
164    p.set_x(x);
165    p.set_y(y);
166}
167
168/// Equal Earth Spheroidal inverse project
169pub fn robin_s_inverse<P: TransformCoordinates>(p: &mut P) {
170    let mut lam = p.x() / FXC;
171    let mut phi = fabs(p.y() / FYC);
172    if phi >= 1. {
173        // simple pathologic cases
174        if phi > ONEEPS {
175            panic!("Coordinate outside projection domain");
176        } else {
177            phi = if p.y() < 0. { -FRAC_PI_2 } else { FRAC_PI_2 };
178            lam /= X_COEFS[NODES].c0;
179        }
180    } else {
181        // general problem
182        // in Y space, reduce to table interval
183        let i = if f64::is_nan(phi) { -1. } else { round(floor(phi * NODES as f64)) };
184        if i < 0. || i >= NODES as f64 {
185            panic!("Coordinate outside projection domain");
186        }
187        let mut i = i as usize;
188        loop {
189            if Y_COEFS[i].c0 > phi {
190                i -= 1;
191            } else if Y_COEFS[i + 1].c0 <= phi {
192                i += 1;
193            } else {
194                break;
195            }
196        }
197        let t_coef = &Y_COEFS[i];
198        // first guess, linear interp
199        let mut t = 5. * (phi - t_coef.c0) / (Y_COEFS[i + 1].c0 - t_coef.c0);
200        let mut iters = MAX_ITER;
201        while iters > 0 {
202            // Newton-Raphson
203            let t1 = (v(t_coef, t) - phi) / dv(t_coef, t);
204            t -= t1;
205            if fabs(t1) < EPS {
206                break;
207            }
208            iters -= 1;
209        }
210        if iters == 0 {
211            panic!("Coordinate outside projection domain");
212        }
213        phi = (5. * (i as f64) + t).to_radians();
214        if p.y() < 0. {
215            phi = -phi;
216        }
217        lam /= v(&X_COEFS[i], t);
218        if fabs(lam) > PI {
219            panic!("Coordinate outside projection domain");
220        }
221    }
222    p.set_phi(phi);
223    p.set_lam(lam);
224}
gistools/proj/project/robin.rs

gistools/proj/project/
robin.rs
