geodesy/inner_op/
laea.rs

1//! Lambert azimuthal equal area: EPSG coordinate operation method 9820, implemented
2//! following [IOGP, 2019](crate::Bibliography::Iogp19), pp. 78-80
3use crate::authoring::*;
4
5use std::f64::consts::FRAC_PI_2;
6const EPS10: f64 = 1e-10;
7
8// ----- F O R W A R D -----------------------------------------------------------------
9
10fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
11    // Oblique aspect: [IOGP, 2019](crate::Bibliography::Iogp19), pp. 78-80
12    let Ok(xi_0) = op.params.real("xi_0") else {
13        return 0;
14    };
15    let Ok(qp) = op.params.real("qp") else {
16        return 0;
17    };
18    let Ok(rq) = op.params.real("rq") else {
19        return 0;
20    };
21    let Ok(d) = op.params.real("d") else { return 0 };
22
23    let oblique = op.params.boolean("oblique");
24    let north_polar = op.params.boolean("north_polar");
25    let south_polar = op.params.boolean("south_polar");
26
27    let lon_0 = op.params.real("lon_0").unwrap_or(0.).to_radians();
28    let x_0 = op.params.real("x_0").unwrap_or(0.);
29    let y_0 = op.params.real("y_0").unwrap_or(0.);
30    let ellps = op.params.ellps(0);
31    let e = ellps.eccentricity();
32    let a = ellps.semimajor_axis();
33
34    let (sin_xi_0, cos_xi_0) = xi_0.sin_cos();
35
36    let mut successes = 0_usize;
37    let n = operands.len();
38
39    // The polar aspects are fairly simple
40    if north_polar || south_polar {
41        let sign = if north_polar { -1.0 } else { 1.0 };
42        for i in 0..n {
43            let (lon, lat) = operands.xy(i);
44            let (sin_lon, cos_lon) = (lon - lon_0).sin_cos();
45
46            let q = ancillary::qs(lat.sin(), e);
47            let rho = a * (qp + sign * q).sqrt();
48
49            let easting = x_0 + rho * sin_lon;
50            let northing = y_0 + sign * rho * cos_lon;
51            operands.set_xy(i, easting, northing);
52            successes += 1;
53        }
54        return successes;
55    }
56
57    // Either equatorial or oblique aspects
58    for i in 0..n {
59        let (lon, lat) = operands.xy(i);
60        let (sin_lon, cos_lon) = (lon - lon_0).sin_cos();
61
62        // Authalic latitude, 𝜉
63        let xi = (ancillary::qs(lat.sin(), e) / qp).asin();
64        let (sin_xi, cos_xi) = xi.sin_cos();
65
66        let b = if oblique {
67            let factor = 1.0 + sin_xi_0 * sin_xi + (cos_xi_0 * cos_xi * cos_lon);
68            rq * (2.0 / factor).sqrt()
69        } else {
70            1.0
71        };
72
73        let easting = x_0 + (b * d) * (cos_xi * sin_lon);
74        let northing = y_0 + (b / d) * (cos_xi_0 * sin_xi - sin_xi_0 * cos_xi * cos_lon);
75        operands.set_xy(i, easting, northing);
76        successes += 1;
77    }
78
79    successes
80}
81
82// ----- I N V E R S E -----------------------------------------------------------------
83
84fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
85    // Oblique aspect: [IOGP, 2019](crate::Bibliography::Iogp19), pp. 78-80
86    let Ok(xi_0) = op.params.real("xi_0") else {
87        return 0;
88    };
89    let Ok(rq) = op.params.real("rq") else {
90        return 0;
91    };
92    let Ok(d) = op.params.real("d") else { return 0 };
93    let Ok(authalic) = op.params.fourier_coefficients("authalic") else {
94        return 0;
95    };
96
97    let north_polar = op.params.boolean("north_polar");
98    let south_polar = op.params.boolean("south_polar");
99
100    let lon_0 = op.params.real("lon_0").unwrap_or(0.).to_radians();
101    let lat_0 = op.params.real("lat_0").unwrap_or(0.).to_radians();
102    let x_0 = op.params.real("x_0").unwrap_or(0.);
103    let y_0 = op.params.real("y_0").unwrap_or(0.);
104
105    let ellps = op.params.ellps(0);
106    let a = ellps.semimajor_axis();
107    let es = ellps.eccentricity_squared();
108    let e = es.sqrt();
109
110    let (sin_xi_0, cos_xi_0) = xi_0.sin_cos();
111
112    let mut successes = 0_usize;
113    let n = operands.len();
114
115    // The polar aspects are not quite as simple as in the forward case
116    if north_polar || south_polar {
117        let sign = if north_polar { -1.0 } else { 1.0 };
118        for i in 0..n {
119            let (x, y) = operands.xy(i);
120            let rho = (x - x_0).hypot(y - y_0);
121
122            // The authalic latitude is a bit convoluted
123            let denom = a * a * (1.0 - ((1.0 - es) / (2.0 * e)) * ((1.0 - e) / (1.0 + e)).ln());
124            let xi = (-sign) * (1.0 - rho * rho / denom);
125
126            let lon = lon_0 + (x - x_0).atan2(sign * (y - y_0));
127            let lat = ellps.latitude_authalic_to_geographic(xi, &authalic);
128            operands.set_xy(i, lon, lat);
129            successes += 1;
130        }
131        return successes;
132    }
133
134    // Either equatorial or oblique aspects
135    for i in 0..n {
136        let (x, y) = operands.xy(i);
137        let rho = ((x - x_0) / d).hypot(d * (y - y_0));
138
139        // A bit of reality hardening ported from the PROJ implementation
140        if rho < EPS10 {
141            operands.set_xy(i, lon_0, lat_0);
142            successes += 1;
143            continue;
144        }
145
146        // Another case of PROJ reality hardening
147        let asin_argument = 0.5 * rho / rq;
148        if asin_argument.abs() > 1.0 {
149            debug!("LAEA: ({x}, {y}) outside domain");
150            operands.set_xy(i, f64::NAN, f64::NAN);
151            continue;
152        }
153
154        let c = 2.0 * asin_argument.asin();
155        let (sin_c, cos_c) = c.sin_cos();
156
157        // The authalic latitude, 𝜉
158        let xi = (cos_c * sin_xi_0 + (d * (y - y_0) * sin_c * cos_xi_0) / rho).asin();
159
160        let lat = ellps.latitude_authalic_to_geographic(xi, &authalic);
161
162        let num = (x - x_0) * sin_c;
163        let denom = d * rho * cos_xi_0 * cos_c - d * d * (y - y_0) * sin_xi_0 * sin_c;
164        let lon = num.atan2(denom) + lon_0;
165        operands.set_xy(i, lon, lat);
166        successes += 1;
167    }
168
169    successes
170}
171
172// ----- C O N S T R U C T O R ---------------------------------------------------------
173
174#[rustfmt::skip]
175pub const GAMUT: [OpParameter; 6] = [
176    OpParameter::Flag { key: "inv" },
177    OpParameter::Text { key: "ellps", default: Some("GRS80") },
178
179    OpParameter::Real { key: "lat_0", default: Some(0_f64) },
180    OpParameter::Real { key: "lon_0", default: Some(0_f64) },
181
182    OpParameter::Real { key: "x_0",   default: Some(0_f64) },
183    OpParameter::Real { key: "y_0",   default: Some(0_f64) },
184];
185
186pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
187    let def = &parameters.instantiated_as;
188    let mut params = ParsedParameters::new(parameters, &GAMUT)?;
189
190    let lat_0 = params.real("lat_0").unwrap_or(0.).to_radians();
191
192    if lat_0.is_nan() {
193        warn!("LAEA: Bad central latitude!");
194        return Err(Error::BadParam("lat_0".to_string(), def.clone()));
195    }
196
197    let t = lat_0.abs();
198    if t > FRAC_PI_2 + EPS10 {
199        warn!("LAEA: Bad central latitude!");
200        return Err(Error::BadParam("lat_0".to_string(), def.clone()));
201    }
202
203    let polar = (t - FRAC_PI_2).abs() < EPS10;
204    let north = polar && (t > 0.0);
205    let equatorial = !polar && t < EPS10;
206    let oblique = !polar && !equatorial;
207    match (polar, equatorial, north) {
208        (true, _, true) => params.boolean.insert("north_polar"),
209        (true, _, false) => params.boolean.insert("south_polar"),
210        (_, true, _) => params.boolean.insert("equatorial"),
211        _ => params.boolean.insert("oblique"),
212    };
213
214    // --- Precompute some latitude invariant factors ---
215
216    let ellps = params.ellps(0);
217    let a = ellps.semimajor_axis();
218    let es = ellps.eccentricity_squared();
219    let e = es.sqrt();
220    let (sin_phi_0, cos_phi_0) = lat_0.sin_cos();
221
222    // qs for the central parallel
223    let q0 = ancillary::qs(sin_phi_0, e);
224    // qs for the North Pole
225    let qp = ancillary::qs(1.0, e);
226    // Authalic latitude of the central parallel - 𝛽₀ in the IOGP text
227    let xi_0 = (q0 / qp).asin();
228    // Rq in the IOGP text
229    let rq = a * (0.5 * qp).sqrt();
230    // D in the IOGP text
231    let d = if oblique {
232        a * (cos_phi_0 / (1.0 - es * sin_phi_0 * sin_phi_0).sqrt()) / (rq * xi_0.cos())
233    } else if equatorial {
234        rq.recip()
235    } else {
236        a
237    };
238
239    params.real.insert("xi_0", xi_0);
240    params.real.insert("q0", q0);
241    params.real.insert("qp", qp);
242    params.real.insert("rq", rq);
243    params.real.insert("d", d);
244
245    let authalic = ellps.coefficients_for_authalic_latitude_computations();
246    params.fourier_coefficients.insert("authalic", authalic);
247
248    let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv)));
249    Ok(Op {
250        descriptor,
251        params,
252        steps: None,
253    })
254}
255
256// ----- T E S T S ---------------------------------------------------------------------
257
258#[cfg(test)]
259mod tests {
260    use super::*;
261    use float_eq::assert_float_eq;
262
263    #[test]
264    fn laea_oblique() -> Result<(), Error> {
265        let mut ctx = Minimal::default();
266
267        // ETRS-LAEA grid definition
268        let op = ctx.op("laea ellps=GRS80 lat_0=52 lon_0=10  x_0=4321000 y_0=3210000")?;
269
270        // The test point from IOGP
271        let p = Coor2D::geo(50.0, 5.0);
272        let geo = [p];
273        let p = Coor2D::raw(3962799.45, 2999718.85);
274        let projected = [p];
275
276        let mut operands = geo;
277
278        // Forward
279        ctx.apply(op, Fwd, &mut operands)?;
280        assert_float_eq!(operands[0].0, projected[0].0, abs_all <= 0.01);
281        assert!((operands[0][0] - 3962799.45).abs() < 0.01);
282        assert!((operands[0][1] - 2999718.85).abs() < 0.01);
283        ctx.apply(op, Inv, &mut operands)?;
284        assert!((operands[0][0].to_degrees() - 5.0).abs() < 1e-12);
285        assert!((operands[0][1].to_degrees() - 50.).abs() < 1e-12);
286
287        let p = Coor4D::raw(1e30, 1e30, 0.0, 0.0);
288        let mut operands = [p];
289        ctx.apply(op, Inv, &mut operands)?;
290        assert!(operands[0][0].is_nan());
291
292        // Missing test points for the polar aspects
293
294        Ok(())
295    }
296
297    // Test the "if rho < EPS10" branch in the inverse case.
298    // From this issue: https://github.com/busstoptaktik/geodesy/issues/89
299    // reported by @maximkaaa
300    #[test]
301    fn origin() {
302        let mut ctx = Minimal::new();
303        let op = ctx
304            .op("laea lon_0=10 lat_0=52 x_0=4321000 y_0=3210000")
305            .unwrap();
306        let mut data = [Coor2D::geo(52.0, 10.0)];
307        let clone = data;
308        ctx.apply(op, Fwd, &mut data).unwrap();
309        ctx.apply(op, Inv, &mut data).unwrap();
310        assert_eq!(data, clone);
311    }
312}
geodesy/inner_op/laea.rs

geodesy/inner_op/
laea.rs
