gistools/proj/project/
lcc.rs

1use crate::proj::{
2    _msfn, CoordinateStep, EPS10, LAMBERT_CONFORMAL_CONIC_1SP, LAMBERT_CONFORMAL_CONIC_2SP,
3    LATITUDE_OF_FIRST_STANDARD_PARALLEL, LATITUDE_OF_NATURAL_ORIGIN, LATITUDE_OF_PROJECTION_CENTRE,
4    LATITUDE_OF_SECOND_STANDARD_PARALLEL, LONGITUDE_OF_NATURAL_ORIGIN, Proj, ProjValue,
5    ProjectCoordinates, TransformCoordinates, phi2, tsfn,
6};
7use alloc::rc::Rc;
8use core::{
9    cell::RefCell,
10    f64::consts::{FRAC_PI_2, FRAC_PI_4},
11};
12use libm::{atan, atan2, cos, fabs, hypot, log, pow, sin, tan};
13
14/// Lambert Conformal Conic variables
15#[derive(Debug, Default, Clone, PartialEq)]
16pub struct LccData {
17    phi1: f64,
18    phi2: f64,
19    n: f64,
20    rho0: f64,
21    c: f64,
22}
23
24/// Lambert Conformal Conic projection
25///
26/// See [`LambertConformalConicProjection`] for full documentation.
27pub type LambertConformalConic1SPProjection =
28    LambertConformalConicProjection<LAMBERT_CONFORMAL_CONIC_1SP>;
29/// Lambert Conic Conformal (2SP)
30///
31/// See [`LambertConformalConicProjection`] for full documentation.
32pub type LambertConformalConic2SPProjection =
33    LambertConformalConicProjection<LAMBERT_CONFORMAL_CONIC_2SP>;
34
35/// # Lambert Conformal Conic
36///
37///  Lambert Conformal Conic projection (LCC) is a conic map projection
38/// used for aeronautical charts, portions of the State Plane Coordinate
39/// System, and many national and regional mapping systems. It is one of
40/// seven projections introduced by Johann Heinrich Lambert in 1772.
41///
42/// It has several different forms: with one and two standard parallels
43/// (referred to as 1SP and 2SP in EPSG guidance notes). Additionally we
44/// provide "2SP Michigan" form which is very similar to normal 2SP, but
45/// with a scaling factor on the ellipsoid (given as `k_0` parameter).
46/// It is implemented as per EPSG Guidance Note 7-2 (version 54, August
47/// 2018, page 25). It is used in a few systems in the EPSG database which
48/// justifies adding this otherwise non-standard projection.
49///
50/// **Classification**: Conformal conic
51///
52/// **Available forms**: Forward and inverse, spherical and ellipsoidal
53/// - One or two standard parallels (1SP and 2SP).
54/// - "LCC 2SP Michigan" form can be used by setting the `+k_0` parameter to specify ellipsoid scale.
55///
56/// **Defined area**: Best for regions predominantly east-west in extent and located in the middle north or south latitudes.
57///
58/// **Alias**: lcc
59///
60/// **Domain**: 2D
61///
62/// **Input type**: Geodetic coordinates
63///
64/// **Output type**: Projected coordinates
65///
66/// ## Projection String
67/// ```ini
68/// +proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45
69/// ```
70///
71/// ## Required Parameters
72/// - `+lat_1`: Latitude of the first standard parallel.
73///
74/// ## Optional Parameters
75/// - `+lon_0`: Longitude of projection center. Defaults to `0`.
76/// - `+lat_0`: Latitude of projection center. Defaults to `0`.
77/// - `+lat_2`: Latitude of the second standard parallel.
78/// - `+ellps`: Ellipsoid. Defaults to `WGS84`.
79/// - `+R`: Radius of the sphere.
80/// - `+x_0`: False easting. Defaults to `0`.
81/// - `+y_0`: False northing. Defaults to `0`.
82/// - `+k_0`: Scale factor at natural origin (for LCC 1SP) or ellipsoid scale factor (for LCC 2SP Michigan). Defaults to `1.0`.
83///
84/// ## Further reading
85/// - [Wikipedia on Lambert Conformal Conic](https://en.wikipedia.org/wiki/Lambert_conformal_conic_projection)
86/// - [Wolfram Mathworld on Lambert Conformal Conic](http://mathworld.wolfram.com/LambertConformalConicProjection.html)
87/// - [John P. Snyder "Map projections: A working manual"](https://pubs.er.usgs.gov/publication/pp1395)
88/// - [ArcGIS documentation on "Lambert Conformal Conic"](http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/lambert-conformal-conic.htm)
89/// - [EPSG Guidance Note 7-2](http://www.epsg.org/Guidancenotes.aspx)
90///
91/// ![Lambert Conformal Conic](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/lcc.png?raw=true)
92#[derive(Debug, Clone, PartialEq)]
93pub struct LambertConformalConicProjection<const C: i64> {
94    proj: Rc<RefCell<Proj>>,
95    store: RefCell<LccData>,
96}
97impl<const C: i64> ProjectCoordinates for LambertConformalConicProjection<C> {
98    fn code(&self) -> i64 {
99        C
100    }
101    fn name(&self) -> &'static str {
102        "Lambert Conformal Conic"
103    }
104    fn names() -> &'static [&'static str] {
105        &[
106            "Lambert Conic Conformal (1SP)",
107            "Lambert_Conformal_Conic_1SP",
108            "Lambert Conic Conformal (2SP)",
109            "Lambert_Conformal_Conic_2SP",
110            "Lambert Conic Conformal (LCC)",
111            "Lambert_Conformal_Conic",
112            "Lambert Conformal Conic",
113            "LambertConformalConic",
114            "lcc",
115        ]
116    }
117}
118impl<const C: i64> CoordinateStep for LambertConformalConicProjection<C> {
119    fn new(proj: Rc<RefCell<Proj>>) -> Self {
120        // values: 8801, 8802, 8805, 8806, 8807
121        let mut store = LccData::default();
122        {
123            let proj = &mut proj.borrow_mut();
124            let lat_1 = proj
125                .params
126                .get(&LATITUDE_OF_NATURAL_ORIGIN)
127                .or_else(|| proj.params.get(&LATITUDE_OF_FIRST_STANDARD_PARALLEL))
128                .unwrap_or(&ProjValue::default())
129                .f64();
130            let lat_2 = proj
131                .params
132                .get(&LONGITUDE_OF_NATURAL_ORIGIN)
133                .or_else(|| proj.params.get(&LATITUDE_OF_SECOND_STANDARD_PARALLEL))
134                .unwrap_or(&ProjValue::default())
135                .f64();
136            store.phi1 = lat_1;
137            if lat_2 != 0. {
138                store.phi2 = lat_2;
139            } else {
140                store.phi2 = store.phi1;
141                if proj.params.contains_key(&LATITUDE_OF_PROJECTION_CENTRE) {
142                    proj.phi0 = store.phi1;
143                }
144            }
145
146            if fabs(store.phi1 + store.phi2) < EPS10 {
147                panic!("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
148            }
149
150            let mut sinphi = sin(store.phi1);
151            store.n = sinphi;
152            let cosphi = cos(store.phi1);
153
154            if fabs(cosphi) < EPS10 || fabs(store.phi1) >= FRAC_PI_2 {
155                panic!("Invalid value for lat_1: |lat_1| should be < 90°");
156            }
157            if fabs(cos(store.phi2)) < EPS10 || fabs(store.phi2) >= FRAC_PI_2 {
158                panic!("Invalid value for lat_2: |lat_2| should be < 90°");
159            }
160
161            let secant = fabs(store.phi1 - store.phi2) >= EPS10;
162            if proj.es != 0. {
163                let m1 = _msfn(sinphi, cosphi, proj.es);
164                let ml1 = tsfn(store.phi1, sinphi, proj.e);
165                if secant {
166                    // secant cone
167                    sinphi = sin(store.phi2);
168                    store.n = log(m1 / _msfn(sinphi, cos(store.phi2), proj.es));
169                    if store.n == 0. {
170                        panic!("Invalid value for eccentricity");
171                    }
172                    let ml2 = tsfn(store.phi2, sinphi, proj.e);
173                    let denom = log(ml1 / ml2);
174                    if denom == 0. {
175                        panic!("Invalid value for eccentricity");
176                    }
177                    store.n /= denom;
178                }
179                store.rho0 = m1 * pow(ml1, -store.n) / store.n;
180                store.c = store.rho0;
181                store.rho0 *= if fabs(fabs(proj.phi0) - FRAC_PI_2) < EPS10 {
182                    0.
183                } else {
184                    pow(tsfn(proj.phi0, sin(proj.phi0), proj.e), store.n)
185                };
186            } else {
187                if secant {
188                    store.n = log(cosphi / cos(store.phi2))
189                        / log(tan(FRAC_PI_4 + 0.5 * store.phi2) / tan(FRAC_PI_4 + 0.5 * store.phi1));
190                }
191                if store.n == 0. {
192                    // Likely reason is that phi1 / phi2 are too close to zero.
193                    // Can be reproduced with +proj=lcc +a=1 +lat_2=.0000001
194                    panic!("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
195                }
196                store.c = cosphi * pow(tan(FRAC_PI_4 + 0.5 * store.phi1), store.n) / store.n;
197                store.rho0 = if fabs(fabs(proj.phi0) - FRAC_PI_2) < EPS10 {
198                    0.
199                } else {
200                    store.c * pow(tan(FRAC_PI_4 + 0.5 * proj.phi0), -store.n)
201                };
202            }
203        }
204        LambertConformalConicProjection { proj, store: store.into() }
205    }
206    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
207        lcc_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
208    }
209    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
210        lcc_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
211    }
212}
213
214/// Lambert Conformal Conic Ellipsoidal forward project
215pub fn lcc_e_forward<P: TransformCoordinates>(lcc: &mut LccData, proj: &Proj, p: &mut P) {
216    let rho = if fabs(fabs(p.phi()) - FRAC_PI_2) < EPS10 {
217        if (p.phi() * lcc.n) <= 0. {
218            panic!("Coordinate outside projection domain");
219        }
220        0.
221    } else {
222        lcc.c
223            * (if proj.es != 0. {
224                pow(tsfn(p.phi(), sin(p.phi()), proj.e), lcc.n)
225            } else {
226                pow(tan(FRAC_PI_4 + 0.5 * p.phi()), -lcc.n)
227            })
228    };
229    p.set_lam(p.lam() * lcc.n);
230    p.set_x(proj.k0 * (rho * sin(p.lam())));
231    p.set_y(proj.k0 * (lcc.rho0 - rho * cos(p.lam())));
232}
233
234/// Lambert Conformal Conic Ellipsoidal inverse project
235pub fn lcc_e_inverse<P: TransformCoordinates>(lcc: &mut LccData, proj: &Proj, p: &mut P) {
236    let mut x = p.x() / proj.k0;
237    let mut y = p.y() / proj.k0;
238    let phi;
239    let lam;
240
241    y = lcc.rho0 - y;
242    let mut rho = hypot(x, y);
243    if rho != 0. {
244        if lcc.n < 0. {
245            rho = -rho;
246            x = -x;
247            y = -y;
248        }
249        if proj.es != 0. {
250            phi = phi2(pow(rho / lcc.c, 1. / lcc.n), proj.e);
251            if phi == f64::MAX {
252                panic!("Coordinate outside projection domain");
253            }
254        } else {
255            phi = 2. * atan(pow(lcc.c / rho, 1. / lcc.n)) - FRAC_PI_2;
256        }
257        lam = atan2(x, y) / lcc.n;
258    } else {
259        lam = 0.;
260        phi = if lcc.n > 0. { FRAC_PI_2 } else { -FRAC_PI_2 };
261    }
262
263    p.set_phi(phi);
264    p.set_lam(lam);
265}
gistools/proj/project/lcc.rs

gistools/proj/project/
lcc.rs
