miniproj_ops/ops/
transverse_mercator.rs

1//This file is licensed under EUPL v1.2 as part of the Digital Earth Viewer
2
3use crate::{
4    ellipsoid::Ellipsoid, traits::GetterContstruct, DbContstruct, Projection, PseudoSerialize,
5};
6
7#[derive(Copy, Clone, Debug)]
8pub struct TransverseMercatorParams {
9    /// longitude of natural origin
10    lon_orig: f64,
11    /// latitude of natural origin
12    lat_orig: f64,
13    /// scale factor at natural origin
14    k_orig: f64,
15    /// false easting
16    false_e: f64,
17    /// false northing
18    false_n: f64,
19}
20
21impl TransverseMercatorParams {
22    pub const fn new(
23        lon_orig: f64,
24        lat_orig: f64,
25        k_orig: f64,
26        false_e: f64,
27        false_n: f64,
28    ) -> Self {
29        Self {
30            lat_orig,
31            lon_orig,
32            k_orig,
33            false_e,
34            false_n,
35        }
36    }
37
38    /// Get longitude of natural origin, radians.
39    pub fn lon_orig(&self) -> f64 {
40        self.lon_orig
41    }
42
43    /// Get latitude of natural origin, radians.
44    pub fn lat_orig(&self) -> f64 {
45        self.lat_orig
46    }
47
48    /// Get scale factor at natural origin.
49    pub fn k_orig(&self) -> f64 {
50        self.k_orig
51    }
52
53    /// Get false easting.
54    pub fn false_e(&self) -> f64 {
55        self.false_e
56    }
57
58    /// Get false northing.
59    pub fn false_n(&self) -> f64 {
60        self.false_n
61    }
62}
63
64/// Transverse Mercator coordinate operation (EPSG:9807).
65#[allow(non_snake_case)]
66#[derive(Copy, Clone, Debug)]
67pub struct TransverseMercatorProjection {
68    pub ellipsoid_e: f64,
69
70    pub lon_orig: f64,
71    pub false_e: f64,
72    pub false_n: f64,
73    pub k_orig: f64,
74
75    pub B: f64,
76    pub h_1: f64,
77    pub h_2: f64,
78    pub h_3: f64,
79    pub h_4: f64,
80    pub M_orig: f64,
81
82    pub h_1_: f64,
83    pub h_2_: f64,
84    pub h_3_: f64,
85    pub h_4_: f64,
86}
87
88impl TransverseMercatorProjection {
89    const MAX_ITERATIONS: usize = 4;
90
91    #[allow(non_snake_case)]
92    pub fn new(ell: &Ellipsoid, params: &TransverseMercatorParams) -> Self {
93        let n = ell.f() / (2.0 - ell.f());
94        let B = (ell.a() / (1.0 + n)) * (1.0 + n.powi(2) / 4.0 + n.powi(4) / 64.0);
95
96        let h_1 = n / 2.0 - (2.0 / 3.0) * n.powi(2)
97            + (5.0 / 16.0) * n.powi(3)
98            + (41.0 / 180.0) * n.powi(4);
99        let h_2 =
100            (13.0 / 48.0) * n.powi(2) - (3.0 / 5.0) * n.powi(3) + (557.0 / 1440.0) * n.powi(4);
101        let h_3 = (61.0 / 240.0) * n.powi(3) - (103.0 / 140.0) * n.powi(4);
102        let h_4 = (49561.0 / 161280.0) * n.powi(4);
103
104        let M_orig = if params.lat_orig() == 0.0 {
105            0.0
106        } else if params.lat_orig() == std::f64::consts::FRAC_PI_2 {
107            B * std::f64::consts::FRAC_PI_2
108        } else if params.lat_orig() == -std::f64::consts::FRAC_PI_2 {
109            -B * std::f64::consts::FRAC_PI_2
110        } else {
111            let Q_orig = params.lat_orig().tan().asinh()
112                - (ell.e() * f64::atanh(ell.e() * params.lat_orig().sin()));
113
114            let beta_orig = Q_orig.sinh().atan();
115            let xi_orig_0 = beta_orig;
116
117            let xi_orig_1 = h_1 * f64::sin(2.0 * xi_orig_0);
118            let xi_orig_2 = h_2 * f64::sin(4.0 * xi_orig_0);
119            let xi_orig_3 = h_3 * f64::sin(6.0 * xi_orig_0);
120            let xi_orig_4 = h_4 * f64::sin(8.0 * xi_orig_0);
121            let xi_orig = xi_orig_0 + xi_orig_1 + xi_orig_2 + xi_orig_3 + xi_orig_4;
122            B * xi_orig
123        };
124
125        let h_1_ = n / 2.0 - (2.0 / 3.0) * n.powi(2) + (37.0 / 96.0) * n.powi(3)
126            - (1.0 / 360.0) * n.powi(4);
127        let h_2_ =
128            (1.0 / 48.0) * n.powi(2) + (1.0 / 15.0) * n.powi(3) - (437.0 / 1440.0) * n.powi(4);
129        let h_3_ = (17.0 / 480.0) * n.powi(3) - (37.0 / 840.0) * n.powi(4);
130        let h_4_ = (4397.0 / 161280.0) * n.powi(4);
131
132        Self {
133            ellipsoid_e: ell.e(),
134            lon_orig: params.lon_orig(),
135            false_e: params.false_e(),
136            false_n: params.false_n(),
137            k_orig: params.k_orig(),
138
139            B,
140            h_1,
141            h_2,
142            h_3,
143            h_4,
144            M_orig,
145
146            h_1_,
147            h_2_,
148            h_3_,
149            h_4_,
150        }
151    }
152}
153
154impl Projection for TransverseMercatorProjection {
155    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
156    /// longitude & latitude in radians
157    #[allow(non_snake_case)]
158    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
159        let Q = latitude.tan().asinh()
160            - (self.ellipsoid_e * f64::atanh(self.ellipsoid_e * latitude.sin()));
161        let beta = Q.sinh().atan();
162        let eta_0 = f64::atanh(beta.cos() * f64::sin(longitude - self.lon_orig));
163        let xi_0 = f64::asin(beta.sin() * eta_0.cosh());
164
165        let xi_1 = self.h_1 * f64::sin(2.0 * xi_0) * f64::cosh(2.0 * eta_0);
166        let xi_2 = self.h_2 * f64::sin(4.0 * xi_0) * f64::cosh(4.0 * eta_0);
167        let xi_3 = self.h_3 * f64::sin(6.0 * xi_0) * f64::cosh(6.0 * eta_0);
168        let xi_4 = self.h_4 * f64::sin(8.0 * xi_0) * f64::cosh(8.0 * eta_0);
169        let xi = xi_0 + xi_1 + xi_2 + xi_3 + xi_4;
170
171        let eta_1 = self.h_1 * f64::cos(2.0 * xi_0) * f64::sinh(2.0 * eta_0);
172        let eta_2 = self.h_2 * f64::cos(4.0 * xi_0) * f64::sinh(4.0 * eta_0);
173        let eta_3 = self.h_3 * f64::cos(6.0 * xi_0) * f64::sinh(6.0 * eta_0);
174        let eta_4 = self.h_4 * f64::cos(8.0 * xi_0) * f64::sinh(8.0 * eta_0);
175        let eta = eta_0 + eta_1 + eta_2 + eta_3 + eta_4;
176
177        (
178            self.false_e + self.k_orig * self.B * eta,
179            self.false_n + self.k_orig * (self.B * xi - self.M_orig),
180        )
181    }
182
183    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
184    /// longitude & latitude in radians
185    #[allow(non_snake_case)]
186    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
187        let eta_ = (easting - self.false_e) / (self.B * self.k_orig);
188        let xi_ = ((northing - self.false_n) + self.k_orig * self.M_orig) / (self.B * self.k_orig);
189
190        let xi_1_ = self.h_1_ * f64::sin(2.0 * xi_) * f64::cosh(2.0 * eta_);
191        let xi_2_ = self.h_2_ * f64::sin(4.0 * xi_) * f64::cosh(4.0 * eta_);
192        let xi_3_ = self.h_3_ * f64::sin(6.0 * xi_) * f64::cosh(6.0 * eta_);
193        let xi_4_ = self.h_4_ * f64::sin(8.0 * xi_) * f64::cosh(8.0 * eta_);
194        let xi_0_ = xi_ - (xi_1_ + xi_2_ + xi_3_ + xi_4_);
195
196        let eta_1_ = self.h_1_ * f64::cos(2.0 * xi_) * f64::sinh(2.0 * eta_);
197        let eta_2_ = self.h_2_ * f64::cos(4.0 * xi_) * f64::sinh(4.0 * eta_);
198        let eta_3_ = self.h_3_ * f64::cos(6.0 * xi_) * f64::sinh(6.0 * eta_);
199        let eta_4_ = self.h_4_ * f64::cos(8.0 * xi_) * f64::sinh(8.0 * eta_);
200        let eta_0_ = eta_ - (eta_1_ + eta_2_ + eta_3_ + eta_4_);
201
202        let beta_ = f64::asin(xi_0_.sin() / eta_0_.cosh());
203        let Q_ = beta_.tan().asinh();
204        let mut Q__ = Q_ + (self.ellipsoid_e * f64::atanh(self.ellipsoid_e * Q_.tanh()));
205        for _ in 0..Self::MAX_ITERATIONS {
206            Q__ = Q_ + (self.ellipsoid_e * f64::atanh(self.ellipsoid_e * Q__.tanh()));
207        }
208
209        (
210            self.lon_orig + f64::asin(eta_0_.tanh() / beta_.cos()),
211            Q__.sinh().atan(),
212        )
213    }
214}
215
216impl PseudoSerialize for TransverseMercatorProjection {
217    fn to_constructed(&self) -> String {
218        format!(
219            r"TransverseMercatorProjection{{
220    ellipsoid_e: {}f64,
221    lon_orig: {}f64,
222    false_e: {}f64,
223    false_n: {}f64,
224    k_orig: {}f64,
225
226    B: {}f64,
227    h_1: {}f64,
228    h_2: {}f64,
229    h_3: {}f64,
230    h_4: {}f64,
231    M_orig: {}f64,
232
233    h_1_: {}f64,
234    h_2_: {}f64,
235    h_3_: {}f64,
236    h_4_: {}f64,
237}}",
238            self.ellipsoid_e,
239            self.lon_orig,
240            self.false_e,
241            self.false_n,
242            self.k_orig,
243            self.B,
244            self.h_1,
245            self.h_2,
246            self.h_3,
247            self.h_4,
248            self.M_orig,
249            self.h_1_,
250            self.h_2_,
251            self.h_3_,
252            self.h_4_
253        )
254    }
255}
256
257impl DbContstruct for TransverseMercatorProjection {
258    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
259        /*
260        ImplementedProjection::new(
261            9807,
262            // lon   lat     k     e     n
263            &[8802, 8801, 8805, 8806, 8807],
264            "TransverseMercatorParams",
265            "TransverseMercatorProjection"
266        ),
267        */
268        let params = TransverseMercatorParams::new(
269            params
270                .iter()
271                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
272                .unwrap(),
273            params
274                .iter()
275                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
276                .unwrap(),
277            params
278                .iter()
279                .find_map(|(c, v)| if *c == 8805 { Some(*v) } else { None })
280                .unwrap(),
281            params
282                .iter()
283                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
284                .unwrap(),
285            params
286                .iter()
287                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
288                .unwrap(),
289        );
290        Self::new(ellipsoid, &params)
291    }
292}
293
294impl GetterContstruct for TransverseMercatorProjection {
295    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
296    where
297        G: FnMut(u32) -> Option<f64>,
298    {
299        let params = TransverseMercatorParams::new(
300            getter(8802)?,
301            getter(8801)?,
302            getter(8805)?,
303            getter(8806)?,
304            getter(8807)?,
305        );
306        Some(Self::new(ellipsoid, &params))
307    }
308}
309
310pub fn direct_projection(params: &[(u32, f64)], ell: Ellipsoid) -> String {
311    TransverseMercatorProjection::from_database_params(params, &ell).to_constructed()
312}
313#[cfg(test)]
314mod tests {
315
316    use crate::ellipsoid::Ellipsoid;
317    use crate::traits::*;
318    use crate::transverse_mercator::*;
319
320    #[test]
321    fn transverse_mercator_consistency() {
322        let wgs_84_ellipsoid = Ellipsoid::from_a_f_inv(6378137.0, 298.257223563);
323        let utm_32_n = TransverseMercatorParams::new(
324            9.0f64.to_radians(),
325            0.0f64.to_radians(),
326            0.9996,
327            500_000.0,
328            0.0,
329        );
330
331        let projection = TransverseMercatorProjection::new(&wgs_84_ellipsoid, &utm_32_n);
332        let easting_goal = 577274.99;
333        let northing_goal = 69740.50;
334        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
335        let (easting, northing) = projection.deg_to_projected(lon, lat);
336
337        eprintln!("easting: {easting_goal} - {easting}");
338        eprintln!("northing: {northing_goal} - {northing}");
339
340        assert!((easting - easting_goal).abs() < 0.01);
341
342        assert!((northing - northing_goal).abs() < 0.01);
343    }
344}
miniproj_ops/ops/transverse_mercator.rs

miniproj_ops/ops/
transverse_mercator.rs
