geodesy/inner_op/
merc.rs

1//! Mercator
2use crate::authoring::*;
3
4// ----- F O R W A R D -----------------------------------------------------------------
5
6fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
7    let ellps = op.params.ellps(0);
8    let a = ellps.semimajor_axis();
9    let k_0 = op.params.k(0);
10    let x_0 = op.params.x(0);
11    let y_0 = op.params.y(0);
12    let lat_0 = op.params.lat(0);
13    let lon_0 = op.params.lon(0);
14
15    let mut successes = 0_usize;
16    for i in 0..operands.len() {
17        let (lon, lat) = operands.xy(i);
18
19        let easting = (lon - lon_0) * k_0 * a - x_0;
20        let isometric = ellps.latitude_geographic_to_isometric(lat + lat_0);
21        let northing = a * k_0 * isometric - y_0;
22
23        operands.set_xy(i, easting, northing);
24        successes += 1;
25    }
26
27    successes
28}
29
30// ----- I N V E R S E -----------------------------------------------------------------
31
32fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
33    let ellps = op.params.ellps(0);
34    let a = ellps.semimajor_axis();
35    let k_0 = op.params.k(0);
36    let x_0 = op.params.x(0);
37    let y_0 = op.params.y(0);
38    let lat_0 = op.params.lat(0);
39    let lon_0 = op.params.lon(0);
40
41    let mut successes = 0_usize;
42    for i in 0..operands.len() {
43        let (mut x, mut y) = operands.xy(i);
44
45        // Easting -> Longitude
46        x += x_0;
47        let lon = x / (a * k_0) - lon_0;
48
49        // Northing -> Latitude
50        y += y_0;
51        let psi = y / (a * k_0);
52        let lat = ellps.latitude_isometric_to_geographic(psi) - lat_0;
53        operands.set_xy(i, lon, lat);
54        successes += 1;
55    }
56
57    successes
58}
59
60// ----- C O N S T R U C T O R ---------------------------------------------------------
61
62#[rustfmt::skip]
63pub const GAMUT: [OpParameter; 8] = [
64    OpParameter::Flag { key: "inv" },
65    OpParameter::Text { key: "ellps",  default: Some("GRS80") },
66
67    OpParameter::Real { key: "lat_0",  default: Some(0_f64) },
68    OpParameter::Real { key: "lon_0",  default: Some(0_f64) },
69    OpParameter::Real { key: "x_0",    default: Some(0_f64) },
70    OpParameter::Real { key: "y_0",    default: Some(0_f64) },
71
72    OpParameter::Real { key: "k_0",    default: Some(1_f64) },
73    OpParameter::Real { key: "lat_ts", default: Some(0_f64) },
74];
75
76pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
77    let def = &parameters.definition;
78    let mut params = ParsedParameters::new(parameters, &GAMUT)?;
79    let ellps = params.ellps(0);
80
81    let lat_ts = params.real("lat_ts")?;
82    if lat_ts.abs() > 90. {
83        return Err(Error::General(
84            "Merc: Invalid value for lat_ts: |lat_ts| should be <= 90°",
85        ));
86    }
87
88    // lat_ts trumps k_0
89    if lat_ts != 0.0 {
90        let sc = lat_ts.to_radians().sin_cos();
91        let k_0 = sc.1 / (1. - ellps.eccentricity_squared() * sc.0 * sc.0).sqrt();
92        params.real.insert("k_0", k_0);
93    }
94
95    let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv)));
96    let steps = Vec::<Op>::new();
97    let id = OpHandle::new();
98
99    Ok(Op {
100        descriptor,
101        params,
102        steps,
103        id,
104    })
105}
106
107// ----- T E S T S ---------------------------------------------------------------------
108
109#[cfg(test)]
110mod tests {
111    use super::*;
112
113    #[test]
114    fn merc() -> Result<(), Error> {
115        let mut ctx = Minimal::default();
116        let definition = "merc";
117        let op = ctx.op(definition)?;
118
119        // Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc
120        // followed by quadrant tests from PROJ builtins.gie
121        let geo = [
122            Coor4D::geo(55., 12., 0., 0.),
123            Coor4D::geo(1., 2., 0., 0.),
124            Coor4D::geo(-1., 2., 0., 0.),
125            Coor4D::geo(1., -2., 0., 0.),
126            Coor4D::geo(-1., -2., 0., 0.),
127        ];
128
129        let projected = [
130            Coor4D::raw(1_335_833.889_519_282_8, 7_326_837.714_873_877, 0., 0.),
131            Coor4D::raw(222_638.981_586_547, 110_579.965_218_249, 0., 0.),
132            Coor4D::raw(222_638.981_586_547, -110_579.965_218_249, 0., 0.),
133            Coor4D::raw(-222_638.981_586_547, 110_579.965_218_249, 0., 0.),
134            Coor4D::raw(-222_638.981_586_547, -110_579.965_218_249, 0., 0.),
135        ];
136
137        // Forward
138        let mut operands = geo;
139        ctx.apply(op, Fwd, &mut operands)?;
140        for i in 0..operands.len() {
141            assert!(operands[i].hypot2(&projected[i]) < 20e-9);
142        }
143
144        // Roundtrip
145        ctx.apply(op, Inv, &mut operands)?;
146        for i in 0..operands.len() {
147            assert!(operands[i].hypot2(&geo[i]) < 20e-9);
148        }
149
150        Ok(())
151    }
152
153    #[test]
154    fn merc_lat_ts() -> Result<(), Error> {
155        let mut ctx = Minimal::default();
156        let definition = "merc lat_ts=56";
157        let op = ctx.op(definition)?;
158
159        let geo = [Coor4D::geo(55., 12., 0., 0.)];
160
161        // Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc +lat_ts=56
162        let projected = [Coor4D::raw(
163            748_713.257_925_886_8,
164            4_106_573.862_841_270_4,
165            0.,
166            0.,
167        )];
168
169        // Forward
170        let mut operands = geo;
171        ctx.apply(op, Fwd, &mut operands)?;
172        for i in 0..operands.len() {
173            assert!(operands[i].hypot2(&projected[i]) < 20e-9);
174        }
175
176        // Roundtrip
177        ctx.apply(op, Inv, &mut operands)?;
178        for i in 0..operands.len() {
179            assert!(operands[i].hypot2(&geo[i]) < 20e-9);
180        }
181
182        Ok(())
183    }
184}
geodesy/inner_op/merc.rs

geodesy/inner_op/
merc.rs
