1use crate::coordinate::{CartesianCoordinate, EllipticalCoordinate, Latitude, Longitude};
6use crate::geo::ellipsoid::Ellipsoid;
7use crate::proj::Projection;
8use core::f64::consts::FRAC_PI_4;
9use irox_units::units::angle::Angle;
10use irox_units::units::length::Length;
11use std::f64::consts::FRAC_PI_2;
12
13fn m(shape: &Ellipsoid, phi: &Angle) -> f64 {
14 let phi = phi.as_radians().value();
15 let sin2 = phi.sin().powi(2);
16 let e2 = shape.first_eccentricity_squared();
17 let cos = phi.cos();
18 cos / (1. - e2 * sin2).sqrt()
19}
20fn t(shape: &Ellipsoid, phi: &Angle) -> f64 {
21 let phi = phi.as_radians().value();
22 let u = (FRAC_PI_4 - phi / 2.).tan();
23 let e = shape.first_eccentricity();
24 let sinphi = phi.sin();
25 let low = ((1. - e * sinphi) / (1. + e * sinphi)).powf(e / 2.);
26 u / low
27}
28#[derive(Default)]
29pub struct LambertConformalConicBuilder {
30 center: EllipticalCoordinate,
31 shape: Ellipsoid,
32 first_parallel: Latitude,
33 second_parallel: Latitude,
34 false_northing: Length,
35 false_easting: Length,
36}
37impl LambertConformalConicBuilder {
38 #[must_use]
39 pub fn with_center(mut self, center: EllipticalCoordinate) -> Self {
40 self.center = center;
41 self
42 }
43 #[must_use]
44 pub fn with_shape(mut self, shape: Ellipsoid) -> Self {
45 self.shape = shape;
46 self
47 }
48 #[must_use]
49 pub fn with_first_parallel(mut self, first_parallel: Latitude) -> Self {
50 self.first_parallel = first_parallel;
51 self
52 }
53 #[must_use]
54 pub fn with_second_parallel(mut self, second_parallel: Latitude) -> Self {
55 self.second_parallel = second_parallel;
56 self
57 }
58 #[must_use]
59 pub fn with_false_northing(mut self, false_northing: Length) -> Self {
60 self.false_northing = false_northing;
61 self
62 }
63 #[must_use]
64 pub fn with_false_easting(mut self, false_easting: Length) -> Self {
65 self.false_easting = false_easting;
66 self
67 }
68 #[must_use]
69 pub fn build(self) -> LambertConformalConic {
70 let phi0 = self.center.get_latitude();
71 let phi1 = self.first_parallel;
72 let phi2 = self.second_parallel;
73
74 let a = self.shape.semi_major_axis_a().as_meters().value();
75
76 let t0 = t(&self.shape, phi0);
77 let t1 = t(&self.shape, &phi1);
78 let t2 = t(&self.shape, &phi2);
79
80 let m1 = m(&self.shape, &phi1);
81 let m2 = m(&self.shape, &phi2);
82
83 let n = (m1.ln() - m2.ln()) / (t1.ln() - t2.ln());
84 let f = m1 / (n * t1.powf(n));
85 let p0 = a * f * t0.powf(n);
86 LambertConformalConic {
87 a,
88 n,
89 f,
90 p0,
91 center: self.center,
92 false_easting: self.false_easting,
93 false_northing: self.false_northing,
94 shape: self.shape,
95 }
96 }
97}
98
99pub struct LambertConformalConic {
100 center: EllipticalCoordinate,
101 shape: Ellipsoid,
102 false_northing: Length,
103 false_easting: Length,
104 p0: f64,
105 a: f64,
106 n: f64,
107 f: f64,
108}
109
110impl Projection for LambertConformalConic {
111 fn get_center_coords(&self) -> &EllipticalCoordinate {
112 &self.center
113 }
114
115 fn project_to_cartesian(&self, coord: &EllipticalCoordinate) -> CartesianCoordinate {
116 let n = self.n;
117 let p0 = self.p0;
118 let a = self.a;
119 let f = self.f;
120 let lam0 = self.center.get_longitude();
121 let phi = coord.get_latitude();
122 let lam = coord.get_longitude();
123 let t = t(&self.shape, phi);
124 let theta = n * (lam.as_radians() - lam0.as_radians());
125
126 let p = a * f * t.powf(n);
127
128 let x = Length::new_meters(p * theta.sin());
129 let y = Length::new_meters(p0 - p * theta.cos());
130
131 CartesianCoordinate::new(
132 x + self.false_easting,
133 y + self.false_northing,
134 Length::ZERO,
135 )
136 }
137
138 fn project_to_elliptical(&self, coord: &CartesianCoordinate) -> EllipticalCoordinate {
139 let x = coord.get_x().as_meters().value();
140 let y = coord.get_y().as_meters().value();
141
142 let p = (x.powi(2) + (self.p0 - y).powi(2)).sqrt() * self.n.signum();
143 let t = (p / (self.a * self.f)).powf(1.0 / self.n);
144 let theta = x.atan2(self.p0 - y);
145 let lam = theta / self.n + self.center.get_longitude().as_radians().value();
146 let e = self.shape.first_eccentricity();
147
148 let mut phi = FRAC_PI_2 - 2. * t.atan();
149 loop {
150 let phit = FRAC_PI_2
151 - 2. * (t * ((1. - e * phi.sin()) / (1. + e * phi.sin())).powf(e / 2.)).atan();
152 let eps = (phit - phi).abs();
153 phi = phit;
154 if eps < 1e-10 {
155 break;
156 }
157 }
158
159 let lat = Latitude(Angle::new_radians(phi));
160 let lon = Longitude(Angle::new_radians(lam));
161
162 EllipticalCoordinate::new(lat, lon, self.shape.into())
163 }
164}
165
166#[cfg(test)]
167mod test {
168 use crate::coordinate::{CartesianCoordinate, EllipticalCoordinate, Longitude};
169 use crate::geo::standards::StandardShapes;
170 use crate::lcc::{LambertConformalConicBuilder, Latitude};
171 use crate::proj::Projection;
172 use irox_tools::assert_eq_eps;
173 use irox_units::units::angle::Angle;
174 use irox_units::units::length::Length;
175
176 #[test]
177 pub fn test_lambert_conic() {
178 let lcc = LambertConformalConicBuilder::default()
179 .with_first_parallel(Latitude(Angle::new_dms(33, 00, 0.0)))
180 .with_second_parallel(Latitude(Angle::new_dms(45, 00, 0.0)))
181 .with_center(EllipticalCoordinate::new(
182 Latitude(Angle::new_dms(23, 0, 0.0)),
183 Longitude(Angle::new_dms(-96, 00, 0.0)),
184 StandardShapes::NAD27.into(),
185 ))
186 .with_shape(StandardShapes::NAD27.as_ellipsoid())
187 .build();
188 let xyz = lcc.project_to_cartesian(&EllipticalCoordinate::new(
189 Latitude(Angle::new_dms(35, 00, 0.0)),
190 Longitude(Angle::new_dms(-75, 00, 0.0)),
191 StandardShapes::NAD27.into(),
192 ));
193 assert_eq_eps!(1894410.9, xyz.get_x().as_meters().value(), 0.1);
194 assert_eq_eps!(1564649.5, xyz.get_y().as_meters().value(), 0.1);
195
196 let lla = lcc.project_to_elliptical(&CartesianCoordinate::new(
197 Length::new_meters(1894410.9),
198 Length::new_meters(1564649.5),
199 Length::ZERO,
200 ));
201 assert_eq_eps!(35.0, lla.get_latitude().as_degrees().value(), 1e-6);
202 assert_eq_eps!(-75.0, lla.get_longitude().as_degrees().value(), 1e-6);
203 }
204}