1use 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 lon_orig: f64,
12 lat_orig: f64,
14 lat_p1: f64,
16 lat_p2: f64,
18 false_e: f64,
20 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 pub fn lon_orig(&self) -> f64 {
45 self.lon_orig
46 }
47
48 pub fn lat_orig(&self) -> f64 {
50 self.lat_orig
51 }
52
53 pub fn lat_p1(&self) -> f64 {
55 self.lat_p1
56 }
57
58 pub fn lat_p2(&self) -> f64 {
60 self.lat_p2
61 }
62
63 pub fn false_e(&self) -> f64 {
65 self.false_e
66 }
67
68 pub fn false_n(&self) -> f64 {
70 self.false_n
71 }
72}
73
74#[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 #[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 #[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, ¶ms)
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, ¶ms))
266 }
267}
268
269#[derive(Copy, Clone, Debug)]
270pub struct LambertConic1SPAParams {
271 lon_nat_orig: f64,
273 lat_nat_orig: f64,
275 k_nat_orig: f64,
277 false_e: f64,
279 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 pub fn lon_nat_orig(&self) -> f64 {
302 self.lon_nat_orig
303 }
304
305 pub fn lat_nat_orig(&self) -> f64 {
307 self.lat_nat_orig
308 }
309
310 pub fn k_nat_orig(&self) -> f64 {
312 self.k_nat_orig
313 }
314
315 pub fn false_e(&self) -> f64 {
317 self.false_e
318 }
319
320 pub fn false_n(&self) -> f64 {
322 self.false_n
323 }
324}
325
326#[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, ¶ms)
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, ¶ms))
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, ¶ms);
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, ¶ms);
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}