gistools/proj/convert/
cart.rs

1use crate::proj::{
2    CoordinateStep, IoUnits, Proj, ProjectCoordinates, ProjectionTransform, Step,
3    TransformCoordinates,
4};
5use alloc::rc::Rc;
6use core::{cell::RefCell, f64::consts::FRAC_PI_2};
7use libm::{atan, atan2, cos, fabs, sin, sqrt};
8
9/// # CARTESIAN / GEODETIC CONVERSIONS
10///
11/// This material follows:
12///
13/// Bernhard Hofmann-Wellenhof & Helmut Moritz:
14/// Physical Geodesy, 2nd edition.
15/// Springer, 2005.
16///
17/// chapter 5.6: Coordinate transformations
18/// (HM, below),
19///
20/// and
21///
22/// Wikipedia: Geographic Coordinate Conversion,
23/// <https://en.wikipedia.org/wiki/Geographic_coordinate_conversion>
24///
25/// (WP, below).
26///
27/// The cartesian-to-geodetic conversion is based on Bowring's
28/// celebrated method:
29///
30/// B. R. Bowring:
31/// Transformation from spatial to geographical coordinates
32/// Survey Review 23(181), pp. 323-327, 1976
33///
34/// (BB, below),
35///
36/// but could probably use some TLC from a newer and faster
37/// algorithm:
38///
39/// Toshio Fukushima:
40/// Transformation from Cartesian to Geodetic Coordinates
41/// Accelerated by Halley’s Method
42/// Journal of Geodesy, February 2006
43///
44/// (TF, below).
45///
46/// Close to the poles, we avoid singularities by switching to an
47/// approximation requiring knowledge of the geocentric radius
48/// at the given latitude. For this, we use an adaptation of the
49/// formula given in:
50///
51/// Wikipedia: Earth Radius
52/// <https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude>
53/// (Derivation and commentary at <https://gis.stackexchange.com/q/20200>)
54///
55/// (WP2, below)
56///
57/// These routines are probably not as robust at those in
58/// geocent.c, at least they haven't been through as heavy
59/// use as their geocent sisters. Some care has been taken
60/// to avoid singularities, but extreme cases (e.g. setting
61/// es, the squared eccentricity, to 1), will cause havoc.
62#[derive(Debug, Clone, PartialEq)]
63pub struct CartesianConverter {
64    proj: Rc<RefCell<Proj>>,
65}
66impl ProjectCoordinates for CartesianConverter {
67    fn code(&self) -> i64 {
68        -1
69    }
70
71    fn name(&self) -> &'static str {
72        "cartesian"
73    }
74
75    fn names() -> &'static [&'static str] {
76        &["cartesian", "cart"]
77    }
78}
79impl CoordinateStep for CartesianConverter {
80    fn new(proj: Rc<RefCell<Proj>>) -> Self {
81        {
82            let proj = &mut proj.borrow_mut();
83            proj.left = IoUnits::RADIANS;
84            proj.right = IoUnits::CARTESIAN;
85            proj.is_ll = true;
86        }
87        CartesianConverter { proj }
88    }
89    /// Geographical to geocentric
90    fn forward<P: TransformCoordinates>(&self, coords: &mut P) {
91        cartesian(&self.proj.borrow(), coords);
92    }
93    /// Geocentric to geographical
94    fn inverse<P: TransformCoordinates>(&self, coords: &mut P) {
95        geodetic(&self.proj.borrow(), coords);
96    }
97}
98impl From<CartesianConverter> for ProjectionTransform {
99    fn from(c: CartesianConverter) -> ProjectionTransform {
100        ProjectionTransform {
101            proj: c.proj.clone(),
102            method: Step::Cartesian(c.into()),
103            ..Default::default()
104        }
105    }
106}
107
108/// Return the normal radius of curvature of an ellipsoid
109/// with semimajor axis a and squared eccentricity es.
110pub fn normal_radius_of_curvature(a: f64, es: f64, sinphi: f64) -> f64 {
111    if es == 0. {
112        return a;
113    }
114    // This is from WP.  HM formula 2-149 gives an a,b version
115    a / sqrt(1. - es * sinphi * sinphi)
116}
117
118/// Return the geocentric radius at latitude phi, of an ellipsoid
119/// with semimajor axis a and semiminor axis b.
120///
121/// This is from WP2, but uses hypot() for potentially better
122/// numerical robustness
123pub fn geocentric_radius(a: f64, b_div_a: f64, cosphi: f64, sinphi: f64) -> f64 {
124    // Non-optimized version:
125    // let b = a * b_div_a;
126    // return hypot(a * a * cosphi, b * b * sinphi) /
127    //        hypot(a * cosphi, b * sinphi);
128    let cosphi_squared = cosphi * cosphi;
129    let sinphi_squared = sinphi * sinphi;
130    let b_div_a_squared = b_div_a * b_div_a;
131    let b_div_a_squared_mul_sinphi_squared = b_div_a_squared * sinphi_squared;
132    a * sqrt(
133        (cosphi_squared + b_div_a_squared * b_div_a_squared_mul_sinphi_squared)
134            / (cosphi_squared + b_div_a_squared_mul_sinphi_squared),
135    )
136}
137
138/// Cartesian to geodetic
139pub fn cartesian<P: TransformCoordinates>(proj: &Proj, coords: &mut P) {
140    let cosphi: f64 = cos(coords.phi());
141    let sinphi = sin(coords.phi());
142    let n = normal_radius_of_curvature(proj.a, proj.es, sinphi);
143
144    // HM formula 5-27 (z formula follows WP)
145    let z = coords.z();
146    let lam = coords.lam();
147    coords.set_x((n + z) * cosphi * cos(lam));
148    coords.set_y((n + z) * cosphi * sin(lam));
149    coords.set_z((n * (1. - proj.es) + z) * sinphi);
150}
151
152/// Geodetic to cartesian
153pub fn geodetic<P: TransformCoordinates>(proj: &Proj, coords: &mut P) {
154    let phi;
155    // Normalize (x,y,z) to the unit sphere/ellipsoid.
156    let x_div_a = coords.x() * proj.ra;
157    let y_div_a = coords.y() * proj.ra;
158    let z_div_a = coords.z() * proj.ra;
159
160    // Perpendicular distance from point to Z-axis (HM eq. 5-28)
161    let p_div_a = sqrt(x_div_a * x_div_a + y_div_a * y_div_a);
162
163    let b_div_a = 1. - proj.f; // = proj.b / proj.a
164    let p_div_a_b_div_a = p_div_a * b_div_a;
165    let norm = sqrt(z_div_a * z_div_a + p_div_a_b_div_a * p_div_a_b_div_a);
166    let c;
167    let s;
168    if norm != 0. {
169        let inv_norm = 1.0 / norm;
170        c = p_div_a_b_div_a * inv_norm;
171        s = z_div_a * inv_norm;
172    } else {
173        c = 1.;
174        s = 0.;
175    }
176
177    let y_phi = z_div_a + proj.e2s * b_div_a * s * s * s;
178    let x_phi = p_div_a - proj.es * c * c * c;
179    let norm_phi = sqrt(y_phi * y_phi + x_phi * x_phi);
180    let mut cosphi;
181    let mut sinphi;
182    if norm_phi != 0. {
183        let inv_norm_phi = 1.0 / norm_phi;
184        cosphi = x_phi * inv_norm_phi;
185        sinphi = y_phi * inv_norm_phi;
186    } else {
187        cosphi = 1.;
188        sinphi = 0.;
189    }
190    if x_phi <= 0. {
191        // this happen on non-sphere ellipsoid when x,y,z is very close to 0
192        // there is no single solution to the cart->geodetic conversion in
193        // that case, clamp to -90/90 deg and avoid a discontinuous boundary
194        // near the poles
195        phi = if coords.z() >= 0. { FRAC_PI_2 } else { -FRAC_PI_2 };
196        cosphi = 0.;
197        sinphi = if coords.z() >= 0. { 1. } else { -1. };
198    } else {
199        phi = atan(y_phi / x_phi);
200    }
201    let lam = atan2(y_div_a, x_div_a);
202    let z = if cosphi < 1e-6 {
203        // poleward of 89.99994 deg, we avoid division by zero
204        // by computing the height as the cartesian z value
205        // minus the geocentric radius of the Earth at the given
206        // latitude
207        let r = geocentric_radius(proj.a, b_div_a, cosphi, sinphi);
208        fabs(coords.z()) - r
209    } else {
210        let n = normal_radius_of_curvature(proj.a, proj.es, sinphi);
211        proj.a * p_div_a / cosphi - n
212    };
213
214    coords.set_phi(phi);
215    coords.set_lam(lam);
216    coords.set_z(z);
217}
gistools/proj/convert/cart.rs

gistools/proj/convert/
cart.rs
