gistools/proj/convert/
geoc.rs

1use crate::proj::{
2    CoordinateStep, Direction, IoUnits, Proj, ProjectCoordinates, TransformCoordinates,
3};
4use alloc::rc::Rc;
5use core::cell::RefCell;
6use libm::{atan, tan};
7
8/// # Conversion from geographic to geocentric latitude and back.
9///
10/// Convert geographical latitude to geocentric (or the other way round if
11/// direction = PJ_INV)
12///
13/// The conversion involves a call to the tangent function, which goes
14/// through the roof at the poles, so very close (the last centimeter) to the
15/// poles no conversion takes place and the input latitude is copied directly to
16/// the output.
17///
18/// Fortunately, the geocentric latitude converges to the geographical at
19/// the poles, so the difference is negligible.
20///
21/// For the spherical case, the geographical latitude equals the geocentric,
22/// and consequently, the input is copied directly to the output.
23#[derive(Debug, Clone, PartialEq)]
24pub struct GeocentricLatitudeConverter {
25    proj: Rc<RefCell<Proj>>,
26}
27impl ProjectCoordinates for GeocentricLatitudeConverter {
28    fn code(&self) -> i64 {
29        -1
30    }
31
32    fn name(&self) -> &'static str {
33        "geocentric"
34    }
35
36    fn names() -> &'static [&'static str] {
37        &["geocentric", "geoc"]
38    }
39}
40impl CoordinateStep for GeocentricLatitudeConverter {
41    fn new(proj: Rc<RefCell<Proj>>) -> Self {
42        {
43            let proj = &mut proj.borrow_mut();
44            proj.left = IoUnits::RADIANS;
45            proj.right = IoUnits::RADIANS;
46            proj.is_ll = true;
47        }
48        GeocentricLatitudeConverter { proj }
49    }
50    /// Geographical to geocentric
51    fn forward<P: TransformCoordinates>(&self, coords: &mut P) {
52        geocentric_latitude(&self.proj.borrow(), Direction::FWD, coords);
53    }
54    /// Geocentric to geographical
55    fn inverse<P: TransformCoordinates>(&self, coords: &mut P) {
56        geocentric_latitude(&self.proj.borrow(), Direction::INV, coords);
57    }
58}
59
60/// Geocentric latitude conversion function
61pub fn geocentric_latitude<P: TransformCoordinates>(
62    proj: &Proj,
63    direction: Direction,
64    coords: &mut P,
65) {
66    let limit = core::f64::consts::FRAC_PI_2 - 1e-9;
67    let phi = coords.phi();
68    if (phi > limit) || (phi < -limit) || (proj.es == 0.) {
69        return;
70    }
71    if direction == Direction::FWD {
72        coords.set_phi(atan(proj.one_es * tan(phi)));
73    } else {
74        coords.set_phi(atan(proj.rone_es * tan(phi)));
75    }
76}