miniproj_ops/ops/
lambert_conic_conformal.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};
6use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
7
8#[derive(Copy, Clone, Debug)]
9pub struct LambertConic2SPParams {
10    /// longitude of false origin
11    lon_orig: f64,
12    /// latitude of false origin
13    lat_orig: f64,
14    /// latitude of 1st standard parallel
15    lat_p1: f64,
16    /// latitude of 2nd standard parallel
17    lat_p2: f64,
18    /// easting at false origin
19    false_e: f64,
20    /// easting at false northing
21    false_n: f64,
22}
23
24impl LambertConic2SPParams {
25    pub fn new(
26        lon_orig: f64,
27        lat_orig: f64,
28        lat_p1: f64,
29        lat_p2: f64,
30        false_e: f64,
31        false_n: f64,
32    ) -> Self {
33        Self {
34            lat_orig,
35            lon_orig,
36            lat_p1,
37            lat_p2,
38            false_e,
39            false_n,
40        }
41    }
42
43    /// Get longitude of false origin, radians.
44    pub fn lon_orig(&self) -> f64 {
45        self.lon_orig
46    }
47
48    /// Get latitude of false origin, radians.
49    pub fn lat_orig(&self) -> f64 {
50        self.lat_orig
51    }
52
53    /// Get latitude of 1st standard parallel.
54    pub fn lat_p1(&self) -> f64 {
55        self.lat_p1
56    }
57
58    /// Get latitude of 2nd standard parallel.
59    pub fn lat_p2(&self) -> f64 {
60        self.lat_p2
61    }
62
63    /// Get easting at false origin.
64    pub fn false_e(&self) -> f64 {
65        self.false_e
66    }
67
68    /// Get northing at false origin.
69    pub fn false_n(&self) -> f64 {
70        self.false_n
71    }
72}
73
74/// EPSG:9802: Lambert Conic Conformal (2SP) .
75#[allow(non_snake_case)]
76#[derive(Copy, Clone, Debug)]
77pub struct LambertConic2SPProjection {
78    pub ellipsoid_e: f64,
79    pub ellipsoid_a: f64,
80
81    pub lon_orig: f64,
82    pub lat_orig: f64,
83
84    pub false_e: f64,
85    pub false_n: f64,
86
87    pub n: f64,
88    pub r_F: f64,
89    pub F: f64,
90}
91
92impl LambertConic2SPProjection {
93    const MAX_ITERATIONS: usize = 4;
94
95    #[allow(non_snake_case)]
96    pub fn new(ell: &Ellipsoid, params: &LambertConic2SPParams) -> Self {
97        let n;
98        let F;
99        let r_F;
100        if params.lat_p1() == params.lat_p2() {
101            let m_O = params.lat_p1().cos()
102                / (1f64 - ell.e_squared() * params.lat_p1().sin().powi(2)).sqrt();
103
104            let t_O = (FRAC_PI_4 - params.lat_p1() / 2f64).tan()
105                / ((1f64 - ell.e() * params.lat_p1().sin())
106                    / (1f64 + ell.e() * params.lat_p1().sin()))
107                .powf(ell.e() / 2f64);
108            n = params.lat_p1().sin();
109            F = m_O / (n * t_O.powf(n));
110            r_F = ell.a() * F * t_O.powf(n);
111        } else {
112            let m1 = params.lat_p1().cos()
113                / (1f64 - ell.e_squared() * params.lat_p1().sin().powi(2)).sqrt();
114            let m2 = params.lat_p2().cos()
115                / (1f64 - ell.e_squared() * params.lat_p2().sin().powi(2)).sqrt();
116
117            let t1 = (FRAC_PI_4 - params.lat_p1() / 2f64).tan()
118                / ((1f64 - ell.e() * params.lat_p1().sin())
119                    / (1f64 + ell.e() * params.lat_p1().sin()))
120                .powf(ell.e() / 2f64);
121            let t2 = (FRAC_PI_4 - params.lat_p2() / 2f64).tan()
122                / ((1f64 - ell.e() * params.lat_p2().sin())
123                    / (1f64 + ell.e() * params.lat_p2().sin()))
124                .powf(ell.e() / 2f64);
125            let t_F = (FRAC_PI_4 - params.lat_orig() / 2f64).tan()
126                / ((1f64 - ell.e() * params.lat_orig().sin())
127                    / (1f64 + ell.e() * params.lat_orig().sin()))
128                .powf(ell.e() / 2f64);
129            n = (m1.ln() - m2.ln()) / (t1.ln() - t2.ln());
130            F = m1 / (n * t1.powf(n));
131            r_F = ell.a() * F * t_F.powf(n);
132        }
133        Self {
134            ellipsoid_e: ell.e(),
135            ellipsoid_a: ell.a(),
136
137            lon_orig: params.lon_orig(),
138            lat_orig: params.lat_orig(),
139
140            false_e: params.false_e(),
141            false_n: params.false_n(),
142
143            n,
144            r_F,
145            F,
146        }
147    }
148}
149
150impl Projection for LambertConic2SPProjection {
151    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – May 2022
152    /// longitude & latitude in radians
153    #[allow(non_snake_case)]
154    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
155        let t = (FRAC_PI_4 - latitude / 2f64).tan()
156            / ((1f64 - self.ellipsoid_e * latitude.sin())
157                / (1f64 + self.ellipsoid_e * latitude.sin()))
158            .powf(self.ellipsoid_e / 2f64);
159
160        let theta = self.n * (longitude - self.lon_orig);
161
162        let r = self.ellipsoid_a * self.F * t.powf(self.n);
163        (
164            self.false_e + r * theta.sin(),
165            self.false_n + self.r_F - r * theta.cos(),
166        )
167    }
168
169    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – May 2022
170    /// longitude & latitude in radians
171    #[allow(non_snake_case)]
172    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
173        let theta_ = (self.n.signum() * (easting - self.false_e))
174            .atan2(self.n.signum() * (self.r_F - (northing - self.false_n)));
175        let r_ = self.n.signum()
176            * ((easting - self.false_e).powi(2) + (self.r_F - (northing - self.false_n)).powi(2))
177                .sqrt();
178        let t_ = (r_ / (self.ellipsoid_a * self.F)).powf(1f64 / self.n);
179        let mut phi = FRAC_PI_2 - 2.0 * (t_.atan());
180        for _ in 0..Self::MAX_ITERATIONS {
181            phi = FRAC_PI_2
182                - 2.0
183                    * (t_
184                        * ((1f64 - self.ellipsoid_e * phi.sin())
185                            / (1f64 + self.ellipsoid_e * phi.sin()))
186                        .powf(self.ellipsoid_e / 2f64))
187                    .atan()
188        }
189        (theta_ / self.n + self.lon_orig, phi)
190    }
191}
192
193impl PseudoSerialize for LambertConic2SPProjection {
194    fn to_constructed(&self) -> String {
195        format!(
196            r"LambertConic2SPProjection{{
197    ellipsoid_e: {}f64,
198    ellipsoid_a: {}f64,
199    lon_orig: {}f64,
200    lat_orig: {}f64,
201    false_e: {}f64,
202    false_n: {}f64,
203    n: {}f64,
204    r_F: {}f64,
205    F: {}f64,
206}}",
207            self.ellipsoid_e,
208            self.ellipsoid_a,
209            self.lon_orig,
210            self.lat_orig,
211            self.false_e,
212            self.false_n,
213            self.n,
214            self.r_F,
215            self.F
216        )
217    }
218}
219
220impl DbContstruct for LambertConic2SPProjection {
221    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
222        let params = LambertConic2SPParams::new(
223            params
224                .iter()
225                .find_map(|(c, v)| if *c == 8822 { Some(*v) } else { None })
226                .unwrap(),
227            params
228                .iter()
229                .find_map(|(c, v)| if *c == 8821 { Some(*v) } else { None })
230                .unwrap(),
231            params
232                .iter()
233                .find_map(|(c, v)| if *c == 8823 { Some(*v) } else { None })
234                .unwrap(),
235            params
236                .iter()
237                .find_map(|(c, v)| if *c == 8824 { Some(*v) } else { None })
238                .unwrap(),
239            params
240                .iter()
241                .find_map(|(c, v)| if *c == 8826 { Some(*v) } else { None })
242                .unwrap(),
243            params
244                .iter()
245                .find_map(|(c, v)| if *c == 8827 { Some(*v) } else { None })
246                .unwrap(),
247        );
248        Self::new(ellipsoid, &params)
249    }
250}
251
252impl GetterContstruct for LambertConic2SPProjection {
253    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
254    where
255        G: FnMut(u32) -> Option<f64>,
256    {
257        let params = LambertConic2SPParams::new(
258            getter(8822)?,
259            getter(8821)?,
260            getter(8823)?,
261            getter(8824)?,
262            getter(8826)?,
263            getter(8827)?,
264        );
265        Some(Self::new(ellipsoid, &params))
266    }
267}
268
269#[derive(Copy, Clone, Debug)]
270pub struct LambertConic1SPAParams {
271    /// longitude of false origin
272    lon_nat_orig: f64,
273    /// latitude of false origin
274    lat_nat_orig: f64,
275    /// scale factor at natural origin
276    k_nat_orig: f64,
277    /// false easting
278    false_e: f64,
279    /// false northing
280    false_n: f64,
281}
282
283impl LambertConic1SPAParams {
284    pub fn new(
285        lon_nat_orig: f64,
286        lat_nat_orig: f64,
287        k_nat_orig: f64,
288        false_e: f64,
289        false_n: f64,
290    ) -> Self {
291        Self {
292            lon_nat_orig,
293            lat_nat_orig,
294            k_nat_orig,
295            false_e,
296            false_n,
297        }
298    }
299
300    /// Get longitude of natural origin, radians.
301    pub fn lon_nat_orig(&self) -> f64 {
302        self.lon_nat_orig
303    }
304
305    /// Get latitude of natural origin, radians.
306    pub fn lat_nat_orig(&self) -> f64 {
307        self.lat_nat_orig
308    }
309
310    /// Get latitude of 2nd standard parallel.
311    pub fn k_nat_orig(&self) -> f64 {
312        self.k_nat_orig
313    }
314
315    /// Get easting at false origin.
316    pub fn false_e(&self) -> f64 {
317        self.false_e
318    }
319
320    /// Get northing at false origin.
321    pub fn false_n(&self) -> f64 {
322        self.false_n
323    }
324}
325
326/// EPSG:9801: Lambert Conic Conformal (2SP).
327#[allow(non_snake_case)]
328#[derive(Copy, Clone, Debug)]
329pub struct LambertConic1SPAProjection {
330    pub false_e: f64,
331    pub false_n: f64,
332
333    pub r_O: f64,
334    pub lon_O: f64,
335
336    pub n: f64,
337    pub t_r_fac: f64,
338    pub ellipsoid_e: f64,
339}
340
341impl LambertConic1SPAProjection {
342    const MAX_ITERATIONS: usize = 4;
343
344    #[allow(non_snake_case)]
345    pub fn new(ell: &Ellipsoid, params: &LambertConic1SPAParams) -> Self {
346        let m_O = params.lat_nat_orig().cos()
347            / (1f64 - ell.e_squared() * params.lat_nat_orig().sin().powi(2)).sqrt();
348        let t_O = (FRAC_PI_4 - params.lat_nat_orig() / 2f64).tan()
349            / ((1f64 - ell.e() * params.lat_nat_orig().sin())
350                / (1f64 + ell.e() * params.lat_nat_orig().sin()))
351            .powf(ell.e() / 2f64);
352        let n = params.lat_nat_orig.sin();
353        let F = m_O / (n * t_O.powf(n));
354        let r_O = ell.a() * F * t_O.powf(n) * params.k_nat_orig();
355        Self {
356            false_e: params.false_e(),
357            false_n: params.false_n(),
358
359            r_O,
360            lon_O: params.lon_nat_orig(),
361            n,
362            t_r_fac: ell.a() * F * params.k_nat_orig(),
363            ellipsoid_e: ell.e(),
364        }
365    }
366}
367
368impl Projection for LambertConic1SPAProjection {
369    fn projected_to_rad(&self, x: f64, y: f64) -> (f64, f64) {
370        let theta_ = (self.n.signum() * (x - self.false_e))
371            .atan2(self.n.signum() * (self.r_O - (y - self.false_n)));
372        let r_ = self.n.signum()
373            * ((x - self.false_e).powi(2) + (self.r_O - (y - self.false_n)).powi(2)).sqrt();
374        let t_ = (r_ / self.t_r_fac).powf(1f64 / self.n);
375        let mut phi = FRAC_PI_2 - 2f64 * t_.atan();
376        for _ in 0..Self::MAX_ITERATIONS {
377            phi = FRAC_PI_2
378                - 2f64
379                    * (t_
380                        * ((1f64 - self.ellipsoid_e * phi.sin())
381                            / (1f64 + self.ellipsoid_e * phi.sin()))
382                        .powf(self.ellipsoid_e / 2f64))
383                    .atan();
384        }
385        (theta_ / self.n + self.lon_O, phi)
386    }
387
388    fn rad_to_projected(&self, lon: f64, lat: f64) -> (f64, f64) {
389        let t = (FRAC_PI_4 - lat / 2f64).tan()
390            / ((1f64 - self.ellipsoid_e * lat.sin()) / (1f64 + self.ellipsoid_e * lat.sin()))
391                .powf(self.ellipsoid_e / 2f64);
392        let r = self.t_r_fac * t.powf(self.n);
393        let theta = self.n * (lon - self.lon_O);
394        (
395            self.false_e + r * theta.sin(),
396            self.false_n + self.r_O - r * theta.cos(),
397        )
398    }
399}
400
401impl PseudoSerialize for LambertConic1SPAProjection {
402    fn to_constructed(&self) -> String {
403        format!(
404            "LambertConic1SPAProjection {{
405    false_e: {}f64,
406    false_n: {}f64,
407    r_O: {}f64,
408    lon_O: {}f64,
409    n: {}f64,
410    t_r_fac: {}f64,
411    ellipsoid_e: {}f64
412}}
413",
414            self.false_e,
415            self.false_n,
416            self.r_O,
417            self.lon_O,
418            self.n,
419            self.t_r_fac,
420            self.ellipsoid_e
421        )
422    }
423}
424
425impl DbContstruct for LambertConic1SPAProjection {
426    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
427        let params = LambertConic1SPAParams::new(
428            params
429                .iter()
430                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
431                .unwrap(),
432            params
433                .iter()
434                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
435                .unwrap(),
436            params
437                .iter()
438                .find_map(|(c, v)| if *c == 8805 { Some(*v) } else { None })
439                .unwrap(),
440            params
441                .iter()
442                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
443                .unwrap(),
444            params
445                .iter()
446                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
447                .unwrap(),
448        );
449        Self::new(ellipsoid, &params)
450    }
451}
452
453impl GetterContstruct for LambertConic1SPAProjection {
454    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
455    where
456        G: FnMut(u32) -> Option<f64>,
457    {
458        let params = LambertConic1SPAParams::new(
459            getter(8802)?,
460            getter(8801)?,
461            getter(8805)?,
462            getter(8806)?,
463            getter(8807)?,
464        );
465        Some(Self::new(ellipsoid, &params))
466    }
467}
468
469pub fn direct_projection_2sp(params: &[(u32, f64)], ell: Ellipsoid) -> String {
470    LambertConic2SPProjection::from_database_params(params, &ell).to_constructed()
471}
472
473pub fn direct_projection_1sp_a(params: &[(u32, f64)], ell: Ellipsoid) -> String {
474    LambertConic1SPAProjection::from_database_params(params, &ell).to_constructed()
475}
476
477#[cfg(test)]
478mod tests {
479
480    use crate::ellipsoid::Ellipsoid;
481    use crate::lambert_conic_conformal::*;
482    use crate::traits::*;
483
484    #[test]
485    fn lambert_conic_2sp_consistency() {
486        let ell = Ellipsoid::from_a_f_inv(6378160.0, 298.25);
487        let params = LambertConic2SPParams::new(
488            145f64.to_radians(),
489            37f64.to_radians(),
490            36f64.to_radians(),
491            38f64.to_radians(),
492            2_500_000.0,
493            4_500_000.0,
494        );
495
496        let projection = LambertConic2SPProjection::new(&ell, &params);
497        let easting_goal = 2477968.963;
498        let northing_goal = 4416742.535;
499        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
500        let (easting, northing) = projection.deg_to_projected(lon, lat);
501
502        eprintln!("easting: {easting_goal} - {easting}");
503        eprintln!("northing: {northing_goal} - {northing}");
504
505        assert!((easting - easting_goal).abs() < 0.001);
506
507        assert!((northing - northing_goal).abs() < 0.001);
508    }
509
510    #[test]
511    fn lambert_conic_1sp_a_consistency() {
512        let ell = Ellipsoid::from_a_f_inv(6378206.400, 294.97870);
513        let params = LambertConic1SPAParams::new(
514            18f64.to_radians(),
515            -77f64.to_radians(),
516            1.0,
517            2_500_000.0,
518            1_500_000.0,
519        );
520
521        let projection = LambertConic1SPAProjection::new(&ell, &params);
522        let easting_goal = 255966.58;
523        let northing_goal = 142493.51;
524        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
525        let (easting, northing) = projection.deg_to_projected(lon, lat);
526
527        eprintln!("easting: {easting_goal} - {easting}");
528        eprintln!("northing: {northing_goal} - {northing}");
529
530        assert!((easting - easting_goal).abs() < 0.001);
531
532        assert!((northing - northing_goal).abs() < 0.001);
533    }
534}