gistools/proj/project/
lcca.rs

1use crate::proj::{
2    CoordinateStep, Proj, ProjectCoordinates, TransformCoordinates, enfn, inv_mlfn, mlfn,
3};
4use alloc::{rc::Rc, vec::Vec};
5use core::cell::RefCell;
6use libm::{atan2, cos, fabs, sin, sqrt, tan};
7
8// Lambert Conformal Conic Alternative
9// -----------------------------------
10//
11// This is Gerald Evenden's 2003 implementation of an alternative
12// "almost" LCC, which has been in use historically, but which
13// should NOT be used for new projects - i.e: use this implementation
14// if you need interoperability with old data represented in this
15// projection, but not in any other case.
16//
17// The code was originally discussed on the PROJ.4 mailing list in
18// a thread archived over at
19//
20// http://lists.maptools.org/pipermail/proj/2003-March/000644.html
21//
22// It was discussed again in the thread starting at
23//
24// http://lists.maptools.org/pipermail/proj/2017-October/007828.html
25// and continuing at
26// http://lists.maptools.org/pipermail/proj/2017-November/007831.html
27//
28// which prompted Clifford J. Mugnier to add these clarifying notes:
29//
30// The French Army Truncated Cubic Lambert (partially conformal) Conic
31// projection is the Legal system for the projection in France between
32// the late 1800s and 1948 when the French Legislature changed the law
33// to recognize the fully conformal version.
34//
35// It was (might still be in one or two North African prior French
36// Colonies) used in North Africa in Algeria, Tunisia, & Morocco, as
37// well as in Syria during the Levant.
38//
39// Last time I have seen it used was about 30+ years ago in
40// Algeria when it was used to define Lease Block boundaries for
41// Petroleum Exploration & Production.
42//
43// (signed)
44//
45// Clifford J. Mugnier, c.p., c.m.s.
46// Chief of Geodesy
47// LSU Center for GeoInformatics
48// Dept. of Civil Engineering
49// LOUISIANA STATE UNIVERSITY
50
51const MAX_ITER: usize = 10;
52const DEL_TOL: f64 = 1e-12;
53
54/// Lambert Conformal Conic Alternative variables
55#[derive(Debug, Default, Clone, PartialEq)]
56pub struct LccaData {
57    en: Vec<f64>,
58    r0: f64,
59    l: f64,
60    m0: f64,
61    c: f64,
62}
63
64/// Lambert Conformal Conic Alternative Projection
65#[derive(Debug, Clone, PartialEq)]
66pub struct LambertConformalConicAlternativeProjection {
67    proj: Rc<RefCell<Proj>>,
68    store: RefCell<LccaData>,
69}
70impl ProjectCoordinates for LambertConformalConicAlternativeProjection {
71    fn code(&self) -> i64 {
72        -1
73    }
74    fn name(&self) -> &'static str {
75        "Lambert Conformal Conic Alternative"
76    }
77    fn names() -> &'static [&'static str] {
78        &["Lambert Conformal Conic Alternative", "lcca"]
79    }
80}
81impl CoordinateStep for LambertConformalConicAlternativeProjection {
82    fn new(proj: Rc<RefCell<Proj>>) -> Self {
83        let mut store = LccaData::default();
84        {
85            let proj = &mut proj.borrow_mut();
86            store.en = enfn(proj.n);
87            if proj.phi0 == 0. {
88                panic!("Invalid value for lat_0: it should be different from 0.");
89            }
90            store.l = sin(proj.phi0);
91            store.m0 = mlfn(proj.phi0, store.l, cos(proj.phi0), &store.en);
92            let s2p0 = store.l * store.l;
93            let mut r0 = 1. / (1. - proj.es * s2p0);
94            let n0 = sqrt(r0);
95            r0 *= proj.one_es * n0;
96            let tan0 = tan(proj.phi0);
97            store.r0 = n0 / tan0;
98            store.c = 1. / (6. * r0 * n0);
99        }
100        LambertConformalConicAlternativeProjection { proj, store: store.into() }
101    }
102    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
103        lcca_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
104    }
105    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
106        lcca_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
107    }
108}
109
110fn f_s(s: f64, c: f64) -> f64 {
111    // func to compute dr
112    s * (1. + s * s * c)
113}
114
115fn f_sp(s: f64, c: f64) -> f64 {
116    // deriv of fs
117    1. + 3. * s * s * c
118}
119
120/// Lambert Conformal Conic Alternative Ellipsoidal forward project
121pub fn lcca_e_forward<P: TransformCoordinates>(lcca: &mut LccaData, proj: &Proj, p: &mut P) {
122    let s = mlfn(p.phi(), sin(p.phi()), cos(p.phi()), &lcca.en) - lcca.m0;
123    let dr = f_s(s, lcca.c);
124    let r = lcca.r0 - dr;
125    let lam_mul_l = p.lam() * lcca.l;
126    p.set_x(proj.k0 * (r * sin(lam_mul_l)));
127    p.set_y(proj.k0 * (lcca.r0 - r * cos(lam_mul_l)));
128}
129
130/// Lambert Conformal Conic Alternative Ellipsoidal inverse project
131pub fn lcca_e_inverse<P: TransformCoordinates>(lcca: &mut LccaData, proj: &Proj, p: &mut P) {
132    let x = p.x() / proj.k0;
133    let y = p.y() / proj.k0;
134    let theta = atan2(x, lcca.r0 - y);
135    let dr = y - x * tan(0.5 * theta);
136    p.set_lam(theta / lcca.l);
137    let mut s = dr;
138    let mut i = MAX_ITER;
139    while i > 0 {
140        let dif = f_s(s, lcca.c) - dr;
141        s -= dif / f_sp(s, lcca.c);
142        if fabs(dif) < DEL_TOL {
143            break;
144        }
145        i -= 1;
146    }
147    if i == 0 {
148        panic!("Coordinate outside projection domain");
149    }
150    p.set_phi(inv_mlfn(s + lcca.m0, &lcca.en));
151}
gistools/proj/project/lcca.rs

gistools/proj/project/
lcca.rs
