geodesy/inner_op/
btmerc.rs

1//! Transverse Mercator, following [Bowring (1989)](crate::bibliography::Bibliography::Bow89)
2use crate::authoring::*;
3
4// ----- F O R W A R D -----------------------------------------------------------------
5
6// Forward transverse mercator, following Bowring (1989)
7fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
8    let ellps = op.params.ellps(0);
9    let eps = ellps.second_eccentricity_squared();
10    let lat_0 = op.params.lat(0).to_radians();
11    let lon_0 = op.params.lon(0).to_radians();
12    let x_0 = op.params.x(0);
13    let y_0 = op.params.y(0);
14    let k_0 = op.params.k(0);
15
16    let mut successes = 0_usize;
17    let n = operands.len();
18    for i in 0..n {
19        let mut coord = operands.get_coord(i);
20
21        let lat = coord[1] + lat_0;
22        let (s, c) = lat.sin_cos();
23        let cc = c * c;
24        let ss = s * s;
25
26        let dlon = coord[0] - lon_0;
27        let oo = dlon * dlon;
28
29        #[allow(non_snake_case)]
30        let N = ellps.prime_vertical_radius_of_curvature(lat);
31        let z = eps * dlon.powi(3) * c.powi(5) / 6.;
32        let sd2 = (dlon / 2.).sin();
33
34        let theta_2 = (2. * s * c * sd2 * sd2).atan2(ss + cc * dlon.cos());
35
36        // Easting
37        let sd = dlon.sin();
38        coord[0] = x_0 + k_0 * N * ((c * sd).atanh() + z * (1. + oo * (36. * cc - 29.) / 10.));
39
40        // Northing
41        let m = ellps.meridian_latitude_to_distance(lat);
42        let znos4 = z * N * dlon * s / 4.;
43        let ecc = 4. * eps * cc;
44        coord[1] = y_0 + k_0 * (m + N * theta_2 + znos4 * (9. + ecc + oo * (20. * cc - 11.)));
45        operands.set_coord(i, &coord);
46        successes += 1;
47    }
48
49    successes
50}
51
52// ----- I N V E R S E -----------------------------------------------------------------
53
54// Inverse transverse mercator, following Bowring (1989)
55fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
56    let ellps = op.params.ellps(0);
57    let eps = ellps.second_eccentricity_squared();
58    let lat_0 = op.params.lat(0).to_radians();
59    let lon_0 = op.params.lon(0).to_radians();
60    let x_0 = op.params.x(0);
61    let y_0 = op.params.y(0);
62    let k_0 = op.params.k(0);
63
64    let mut successes = 0_usize;
65    let n = operands.len();
66    for i in 0..n {
67        let mut coord = operands.get_coord(i);
68        // Footpoint latitude, i.e. the latitude of a point on the central meridian
69        // having the same northing as the point of interest
70        let lat = ellps.meridian_distance_to_latitude((coord[1] - y_0) / k_0);
71        let (s, c) = lat.sin_cos();
72        let t = s / c;
73        let cc = c * c;
74        #[allow(non_snake_case)]
75        let N = ellps.prime_vertical_radius_of_curvature(lat);
76        let x = (coord[0] - x_0) / (k_0 * N);
77        let xx = x * x;
78        let theta_4 = x.sinh().atan2(c);
79        let theta_5 = (t * theta_4.cos()).atan();
80
81        // Latitude
82        let xet = xx * xx * eps * t / 24.;
83        coord[1] = lat_0 + (1. + cc * eps) * (theta_5 - xet * (9. - 10. * cc)) - eps * cc * lat;
84
85        // Longitude
86        let approx = lon_0 + theta_4;
87        let coef = eps / 60. * xx * x * c;
88        coord[0] = approx - coef * (10. - 4. * xx / cc + xx * cc);
89        operands.set_coord(i, &coord);
90
91        successes += 1;
92    }
93
94    successes
95}
96
97// ----- C O N S T R U C T O R ---------------------------------------------------------
98
99#[rustfmt::skip]
100pub const GAMUT: [OpParameter; 7] = [
101    OpParameter::Flag { key: "inv" },
102    OpParameter::Text { key: "ellps", default: Some("GRS80") },
103
104    OpParameter::Real { key: "lat_0", default: Some(0_f64) },
105    OpParameter::Real { key: "lon_0", default: Some(0_f64) },
106    OpParameter::Real { key: "x_0",   default: Some(0_f64) },
107    OpParameter::Real { key: "y_0",   default: Some(0_f64) },
108
109    OpParameter::Real { key: "k_0",   default: Some(1_f64) },
110];
111
112pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
113    Op::basic(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT)
114}
115
116#[rustfmt::skip]
117pub const UTM_GAMUT: [OpParameter; 4] = [
118    OpParameter::Flag { key: "inv" },
119    OpParameter::Flag { key: "south" },
120    OpParameter::Text { key: "ellps", default: Some("GRS80") },
121    OpParameter::Natural { key: "zone", default: None },
122];
123
124pub fn utm(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
125    let def = &parameters.instantiated_as;
126    let mut params = ParsedParameters::new(parameters, &UTM_GAMUT)?;
127
128    // The UTM zone should be an integer between 1 and 60
129    let zone = params.natural("zone")?;
130    if !(1..61).contains(&zone) {
131        return Err(Error::General(
132            "UTM: 'zone' must be an integer in the interval 1..60",
133        ));
134    }
135
136    // The scaling factor is 0.9996 by definition of UTM
137    params.real.insert("k_0", 0.9996);
138
139    // The center meridian is determined by the zone
140    params.real.insert("lon_0", -183. + 6. * zone as f64);
141
142    // The base parallel is by definition the equator
143    params.real.insert("lat_0", 0.);
144
145    // The false easting is 500000 m by definition of UTM
146    params.real.insert("x_0", 500_000.);
147
148    // The false northing is 0 m by definition of UTM
149    params.real.insert("y_0", 0.);
150    // or 10_000_000 m if using the southern aspect
151    if params.boolean("south") {
152        params.real.insert("y_0", 10_000_000.0);
153    }
154
155    let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv)));
156
157    Ok(Op {
158        descriptor,
159        params,
160        steps: None,
161    })
162}
163
164// ----- T E S T S ---------------------------------------------------------------------
165
166#[cfg(test)]
167mod tests {
168    use super::*;
169
170    #[test]
171    fn btmerc() -> Result<(), Error> {
172        let mut ctx = Minimal::default();
173        let definition = "btmerc k_0=0.9996 lon_0=9 x_0=500000";
174        let op = ctx.op(definition)?;
175
176        // Validation values from PROJ:
177        // echo 12 55 0 0 | cct -d18 +proj=utm +zone=32 | clip
178        #[rustfmt::skip]
179        let geo = [
180            Coor4D::geo( 55.,  12., 0., 0.),
181            Coor4D::geo(-55.,  12., 0., 0.),
182            Coor4D::geo( 55., -6., 0., 0.),
183            Coor4D::geo(-55., -6., 0., 0.)
184        ];
185
186        #[rustfmt::skip]
187        let projected = [
188            Coor4D::raw( 691_875.632_139_661, 6_098_907.825_005_012, 0., 0.),
189            Coor4D::raw( 691_875.632_139_661,-6_098_907.825_005_012, 0., 0.),
190            Coor4D::raw(-455_673.814_189_040, 6_198_246.671_090_279, 0., 0.),
191            Coor4D::raw(-455_673.814_189_040,-6_198_246.671_090_279, 0., 0.)
192        ];
193
194        let mut operands = geo;
195        ctx.apply(op, Fwd, &mut operands)?;
196        for i in 0..operands.len() {
197            assert!(operands[i].hypot2(&projected[i]) < 5e-3);
198        }
199
200        ctx.apply(op, Inv, &mut operands)?;
201        for i in 0..operands.len() {
202            assert!(operands[i].hypot2(&geo[i]) < 10e-8);
203        }
204        Ok(())
205    }
206
207    #[test]
208    fn butm() -> Result<(), Error> {
209        let mut ctx = Minimal::default();
210        let definition = "butm zone=32";
211        let op = ctx.op(definition)?;
212
213        // Validation values from PROJ:
214        // echo 12 55 0 0 | cct -d18 +proj=utm +zone=32 | clip
215        #[rustfmt::skip]
216        let geo = [
217            Coor4D::geo( 55.,  12., 0., 0.),
218            Coor4D::geo(-55.,  12., 0., 0.),
219            Coor4D::geo( 55., -6., 0., 0.),
220            Coor4D::geo(-55., -6., 0., 0.)
221        ];
222
223        #[rustfmt::skip]
224        let projected = [
225            Coor4D::raw( 691_875.632_139_661, 6_098_907.825_005_012, 0., 0.),
226            Coor4D::raw( 691_875.632_139_661,-6_098_907.825_005_012, 0., 0.),
227            Coor4D::raw(-455_673.814_189_040, 6_198_246.671_090_279, 0., 0.),
228            Coor4D::raw(-455_673.814_189_040,-6_198_246.671_090_279, 0., 0.)
229        ];
230
231        let mut operands = geo;
232        ctx.apply(op, Fwd, &mut operands)?;
233        for i in 0..operands.len() {
234            assert!(operands[i].hypot2(&projected[i]) < 5e-3);
235        }
236
237        ctx.apply(op, Inv, &mut operands)?;
238        for i in 0..operands.len() {
239            assert!(operands[i].hypot2(&geo[i]) < 10e-8);
240        }
241        Ok(())
242    }
243
244    #[test]
245    fn butm_south() -> Result<(), Error> {
246        let mut ctx = Minimal::default();
247        let op = ctx.op("butm zone=32 south")?;
248
249        #[rustfmt::skip]
250        let geo = [
251            Coor4D::geo( 55.,  12., 0., 0.),
252            Coor4D::geo(-55.,  12., 0., 0.),
253            Coor4D::geo( 55., -6., 0., 0.),
254            Coor4D::geo(-55., -6., 0., 0.)
255        ];
256
257        #[rustfmt::skip]
258        let projected = [
259            Coor4D::raw( 691_875.632_139_661, 1e7+6_098_907.825_005_012, 0., 0.),
260            Coor4D::raw( 691_875.632_139_661, 1e7-6_098_907.825_005_012, 0., 0.),
261            Coor4D::raw(-455_673.814_189_040, 1e7+6_198_246.671_090_279, 0., 0.),
262            Coor4D::raw(-455_673.814_189_040, 1e7-6_198_246.671_090_279, 0., 0.)
263        ];
264
265        let mut operands = geo;
266        ctx.apply(op, Fwd, &mut operands)?;
267        for i in 0..operands.len() {
268            assert!(operands[i].hypot2(&projected[i]) < 5e-3);
269        }
270
271        ctx.apply(op, Inv, &mut operands)?;
272        for i in 0..operands.len() {
273            assert!(operands[i].hypot2(&geo[i]) < 10e-8);
274        }
275
276        Ok(())
277    }
278}
geodesy/inner_op/btmerc.rs

geodesy/inner_op/
btmerc.rs
