gistools/proj/parse/
datum.rs

1use crate::proj::{
2    Proj, ProjectionTransform, SRS_WGS84_ESQUARED, SRS_WGS84_SEMIMAJOR, SRS_WGS84_SEMIMINOR,
3    TransformCoordinates,
4};
5use alloc::{vec, vec::Vec};
6use core::f64::consts::{FRAC_PI_2, PI, TAU};
7use libm::{atan, atan2, cos, sin, sqrt};
8
9/// Check if the source and destination datums are not WGS84
10pub fn check_not_wgs84(src: &ProjectionTransform, dest: &ProjectionTransform) -> bool {
11    let src_proj = src.proj.borrow();
12    let dest_proj = dest.proj.borrow();
13    !src_proj.datum_type.is_wgs84(&src_proj.datum_params)
14        || !dest_proj.datum_type.is_wgs84(&dest_proj.datum_params)
15}
16
17/// Datum Type helps quicker assessment of a datum
18#[derive(Debug, Default, Clone, Copy, PartialEq, PartialOrd)]
19pub enum DatumType {
20    /// 3 Parameter Datum
21    Param3 = 1,
22    /// 7 Parameter Datum
23    Param7 = 2,
24    /// Grid Shift
25    GridShift = 3,
26    /// WGS84 (Base case)
27    WGS84 = 4,
28    /// Unknown
29    #[default]
30    NoDatum = 5,
31}
32impl DatumType {
33    /// Check if the datum type is 3 or 7
34    pub fn is_params(&self) -> bool {
35        matches!(self, DatumType::Param3 | DatumType::Param7)
36    }
37
38    /// Check if the datum type is WGS84
39    pub fn is_wgs84(&self, params: &DatumParams) -> bool {
40        matches!(self, DatumType::WGS84) || params.is_wgs84()
41    }
42}
43
44/// Datum Parameters can be either 3 or 7
45#[derive(Debug, Copy, Clone, PartialEq, PartialOrd)]
46pub enum DatumParams {
47    /// 3 parameter datum (translate-x, translate-y, translate-z)
48    Param3(f64, f64, f64),
49    /// 7 parameter datum (translate-x, translate-y, translate-z, rotate-x, rotate-y, rotate-z, scale)
50    Param7(f64, f64, f64, f64, f64, f64, f64),
51}
52impl Default for DatumParams {
53    fn default() -> Self {
54        DatumParams::Param3(0.0, 0.0, 0.0)
55    }
56}
57impl DatumParams {
58    /// Returns the datum parameters as a vector
59    pub fn to_vec(&self) -> Vec<f64> {
60        match self {
61            DatumParams::Param3(x, y, z) => vec![*x, *y, *z],
62            DatumParams::Param7(x, y, z, rx, ry, rz, s) => vec![*x, *y, *z, *rx, *ry, *rz, *s],
63        }
64    }
65    /// Given a vector of datum parameters, returns a DatumParams
66    pub fn from_vec(v: Vec<f64>) -> DatumParams {
67        if v.len() == 3 {
68            DatumParams::Param3(v[0], v[1], v[2])
69        } else {
70            DatumParams::Param7(v[0], v[1], v[2], v[3], v[4], v[5], v[6])
71        }
72    }
73    /// Check if the datum is WGS84
74    pub fn is_wgs84(&self) -> bool {
75        match self {
76            Self::Param3(p3_1, p3_2, p3_3) => *p3_1 == 0. && *p3_2 == 0. && *p3_3 == 0.,
77            Self::Param7(p7_1, p7_2, p7_3, p7_4, p7_5, p7_6, p7_7) => {
78                *p7_1 == 0.
79                    && *p7_2 == 0.
80                    && *p7_3 == 0.
81                    && *p7_4 == 0.
82                    && *p7_5 == 0.
83                    && *p7_6 == 0.
84                    && *p7_7 == 0.
85            }
86        }
87    }
88    /// Geocentric to WGS84
89    /// pj_geocentic_to_wgs84( p )
90    /// p = point to transform in geocentric coordinates (x,y,z)
91    /// point object, nothing fancy, just allows values to be
92    /// passed back and forth by reference rather than by value.
93    /// Other point classes may be used as long as they have
94    /// x and y properties, which will get modified in the transform method.
95    pub fn to_wgs84<P: TransformCoordinates>(&self, p: &mut P) {
96        match self {
97            Self::Param3(p3_x, p3_y, p3_z) => {
98                p.set_x(p.x() + *p3_x);
99                p.set_y(p.y() + *p3_y);
100                p.set_z(p.z() - *p3_z);
101            }
102            Self::Param7(p7_dx, p7_dy, p7_dz, p7_rx, p7_ry, p7_rz, p7_s) => {
103                let z = p.z();
104                p.set_x(p7_s * (p.x() - p7_rz * p.y() + p7_ry * z) + *p7_dx);
105                p.set_y(p7_s * (p7_rz * p.x() + p.y() - p7_rx * z) + *p7_dy);
106                p.set_z(p7_s * (-p7_ry * p.x() + p7_rx * p.y() + z) + *p7_dz);
107            }
108        }
109    }
110
111    /// Geocentric from WGS84
112    /// pj_geocentic_from_wgs84() coordinate system definition,
113    /// point to transform in geocentric coordinates (x,y,z)
114    pub fn from_wgs84<P: TransformCoordinates>(&self, p: &mut P) {
115        match self {
116            Self::Param3(p3_x, p3_y, p3_z) => {
117                p.set_x(p.x() - *p3_x);
118                p.set_y(p.y() - *p3_y);
119                p.set_z(p.z() - *p3_z);
120            }
121            Self::Param7(p7_dx, p7_dy, p7_dz, p7_rx, p7_ry, p7_rz, p7_s) => {
122                let z = p.z();
123                let x_tmp = (p.x() - *p7_dx) / *p7_s;
124                let y_tmp = (p.y() - *p7_dy) / *p7_s;
125                let z_tmp = (z - *p7_dz) / *p7_s;
126                p.set_x(x_tmp + *p7_rz * y_tmp - *p7_ry * z_tmp);
127                p.set_y(-p7_rz * x_tmp + y_tmp + *p7_rx * z_tmp);
128                p.set_z(*p7_ry * x_tmp - *p7_rx * y_tmp + z_tmp);
129            }
130        }
131    }
132}
133
134/// Transforms a point from one datum to another
135/// point - lon-lat WGS84 point to mutate
136/// source - source projection
137/// dest - destination projection
138pub fn datum_transform<P: TransformCoordinates>(point: &mut P, source: &Proj, dest: &Proj) {
139    // Short cut if the datums are identical.
140    if source.datum_params == dest.datum_params
141        || source.datum_type == DatumType::NoDatum
142        || dest.datum_type == DatumType::NoDatum
143    {
144        return;
145    }
146
147    // If this datum requires grid shifts, then apply it to geodetic coordinates.
148    let mut source_a = source.a;
149    let mut source_es = source.es;
150    if source.datum_type == DatumType::GridShift {
151        // source
152        // TODO:
153        // apply_grid_shift(source, false, point);
154        source_a = SRS_WGS84_SEMIMAJOR;
155        source_es = SRS_WGS84_ESQUARED;
156    }
157
158    let mut dest_a = dest.a;
159    let mut dest_b = dest.b;
160    let mut dest_es = dest.es;
161    if dest.datum_type == DatumType::GridShift {
162        dest_a = SRS_WGS84_SEMIMAJOR;
163        dest_b = SRS_WGS84_SEMIMINOR;
164        dest_es = SRS_WGS84_ESQUARED;
165    }
166
167    // Do we need to go through geocentric coordinates?
168    if source_es == dest_es
169        && source_a == dest_a
170        && !source.datum_type.is_params()
171        && !dest.datum_type.is_params()
172    {
173        return;
174    }
175
176    // Convert to geocentric coordinates.
177    geodetic_to_geocentric(point, source_es, source_a);
178    // Convert between datums
179    if source.datum_type.is_params() {
180        // geocentric_to_wgs84(point, source.datum_type, source.datum_params);
181        source.datum_params.to_wgs84(point);
182    }
183    if dest.datum_type.is_params() {
184        // geocentric_from_wgs84(point, dest.datum_type, dest.datum_params);
185        dest.datum_params.from_wgs84(point);
186    }
187    // Convert back to geodetic coordinates.
188    geocentric_to_geodetic(point, dest_es, dest_a, dest_b);
189
190    if dest.datum_type == DatumType::GridShift {
191        // TODO:
192        // apply_grid_shift(dest, true, point);
193    }
194}
195
196/// The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
197/// (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
198/// according to the current ellipsoid parameters.
199///
200/// Latitude  : Geodetic latitude in radians                     (input)
201/// Longitude : Geodetic longitude in radians                    (input)
202/// Height    : Geodetic height, in meters                       (input)
203/// X         : Calculated Geocentric X coordinate, in meters    (output)
204/// Y         : Calculated Geocentric Y coordinate, in meters    (output)
205/// Z         : Calculated Geocentric Z coordinate, in meters    (output)
206/// p - lon-lat WGS84 point
207/// es - eccentricity
208/// a - semi-major axis
209pub fn geodetic_to_geocentric<P: TransformCoordinates>(p: &mut P, es: f64, a: f64) {
210    let mut longitude = p.x();
211    let mut latitude = p.y();
212    let height = p.z(); //Z value not always supplied
213    // Don't blow up if Latitude is just a little out of the value
214    // range as it may just be a rounding issue.  Also removed longitude
215    // test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.
216    if latitude < -FRAC_PI_2 && latitude > -1.001 * FRAC_PI_2 {
217        latitude = -FRAC_PI_2;
218    } else if latitude > FRAC_PI_2 && latitude < 1.001 * FRAC_PI_2 {
219        latitude = FRAC_PI_2;
220    } else if !(-FRAC_PI_2..=FRAC_PI_2).contains(&latitude) {
221        panic!("geocent:lat out of range: {}", latitude);
222    }
223
224    if longitude > PI {
225        longitude -= TAU;
226    }
227    let sin_lat = sin(latitude); /*  sin(latitude)  */
228    let cos_lat = cos(latitude); /*  cos(latitude)  */
229    let sin2_lat = sin_lat * sin_lat; /*  Square of sin(latitude)  */
230    let r_n = a / sqrt(1.0 - es * sin2_lat); /*  Earth radius at location  */
231
232    p.set_x((r_n + height) * cos_lat * cos(longitude));
233    p.set_y((r_n + height) * cos_lat * sin(longitude));
234    p.set_z((r_n * (1. - es) + height) * sin_lat);
235}
236
237/// converts a geocentric point to a geodetic point
238pub fn geocentric_to_geodetic<P: TransformCoordinates>(point: &mut P, es: f64, a: f64, _b: f64) {
239    // local defintions and variables
240    // end-criterium of loop, accuracy of sin(Latitude)
241    let genau = 1e-12;
242    let genau2 = genau * genau;
243    let maxiter = 30;
244
245    let mut rx;
246    let mut rk;
247    let mut rn; /* Earth radius at location */
248    let mut cphi0; /* cos of start or old geodetic latitude in iterations */
249    let mut sphi0; /* sin of start or old geodetic latitude in iterations */
250    let mut cphi; /* cos of searched geodetic latitude */
251    let mut sphi; /* sin of searched geodetic latitude */
252    let mut sdphi; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
253    let mut iter; /* # of continous iteration, max. 30 is always enough (s.a.) */
254
255    let x = point.x();
256    let y = point.y();
257    let z = point.z(); //Z value not always supplied
258    let p = sqrt(x * x + y * y); /* distance between semi-minor axis and location */
259    let rr = sqrt(x * x + y * y + z * z); /* distance between center and location */
260    let longitude = if p / a < genau { 0.0 } else { atan2(y, x) };
261    let mut height;
262
263    // --------------------------------------------------------------
264    // Following iterative algorithm was developped by
265    // "Institut for Erdmessung", University of Hannover, July 1988.
266    // Internet: www.ife.uni-hannover.de
267    // Iterative computation of CPHI,SPHI and Height.
268    // Iteration of CPHI and SPHI to 10**-12 radian resp.
269    // 2*10**-7 arcsec.
270    // --------------------------------------------------------------
271    let ct = z / rr; /* sin of geocentric latitude */
272    let st = p / rr; /* cos of geocentric latitude */
273    rx = 1.0 / sqrt(1.0 - es * (2.0 - es) * st * st);
274    cphi0 = st * (1.0 - es) * rx;
275    sphi0 = ct * rx;
276    iter = 0;
277
278    // loop to find sin(Latitude) resp. Latitude
279    // until |sin(Latitude(iter)-Latitude(iter-1))| < genau
280    loop {
281        iter += 1;
282        rn = a / sqrt(1.0 - es * sphi0 * sphi0);
283
284        // ellipsoidal (geodetic) height
285        height = p * cphi0 + z * sphi0 - rn * (1.0 - es * sphi0 * sphi0);
286
287        rk = (es * rn) / (rn + height);
288        rx = 1.0 / sqrt(1.0 - rk * (2.0 - rk) * st * st);
289        cphi = st * (1.0 - rk) * rx;
290        sphi = ct * rx;
291        sdphi = sphi * cphi0 - cphi * sphi0;
292        cphi0 = cphi;
293        sphi0 = sphi;
294        if sdphi * sdphi <= genau2 && iter >= maxiter {
295            break;
296        }
297    }
298
299    // ellipsoidal (geodetic) latitude
300    let latitude = atan(sphi / cphi.abs());
301
302    point.set_x(longitude);
303    point.set_y(latitude);
304    point.set_z(height);
305}
306
307/// Description of a WGS84 datum
308#[derive(Debug, Copy, Clone, PartialEq)]
309pub struct ToWGS84Datum {
310    datum_params: DatumParams,
311    ellipse: &'static str,
312    datum_name: &'static str,
313}
314
315/// WGS84 Datum
316pub const TO_WGS84: ToWGS84Datum = ToWGS84Datum {
317    datum_params: DatumParams::Param7(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
318    ellipse: "WGS84",
319    datum_name: "WGS84",
320};
321
322/// Swiss Datum
323pub const CH1903: ToWGS84Datum = ToWGS84Datum {
324    datum_params: DatumParams::Param3(674.374, 15.056, 405.346),
325    ellipse: "bessel",
326    datum_name: "swiss",
327};
328
329/// Greek_Geodetic_Reference_System_1987 Datum
330pub const GGRS87: ToWGS84Datum = ToWGS84Datum {
331    datum_params: DatumParams::Param3(-199.87, 74.79, 246.62),
332    ellipse: "GRS80",
333    datum_name: "Greek_Geodetic_Reference_System_1987",
334};
335
336/// North_American_Datum_1983 Datum
337pub const NAD83: ToWGS84Datum = ToWGS84Datum {
338    datum_params: DatumParams::Param3(0.0, 0.0, 0.0),
339    ellipse: "GRS80",
340    datum_name: "North_American_Datum_1983",
341};
342
343/// Potsdam Rauenberg 1950 DHDN Datum
344pub const POTSDAM: ToWGS84Datum = ToWGS84Datum {
345    datum_params: DatumParams::Param7(598.1, 73.7, 418.2, 0.202, 0.045, -2.455, 6.7),
346    ellipse: "bessel",
347    datum_name: "Potsdam Rauenberg 1950 DHDN",
348};
349
350/// Carthage 1934 Tunisia Datum
351pub const CARTHAGE: ToWGS84Datum = ToWGS84Datum {
352    datum_params: DatumParams::Param3(-263.0, 6.0, 431.0),
353    ellipse: "clark80",
354    datum_name: "Carthage 1934 Tunisia",
355};
356
357/// Hermannskogel Datum
358pub const HERMANNSKOGEL: ToWGS84Datum = ToWGS84Datum {
359    datum_params: DatumParams::Param7(577.326, 90.129, 463.919, 5.137, 1.474, 5.297, 2.4232),
360    ellipse: "bessel",
361    datum_name: "Hermannskogel",
362};
363
364/// Militar-Geographische Institut Datum
365pub const MILITARGEOGRAPHISCHE_INSTITUT: ToWGS84Datum = ToWGS84Datum {
366    datum_params: DatumParams::Param7(577.326, 90.129, 463.919, 5.137, 1.474, 5.297, 2.4232),
367    ellipse: "bessel",
368    datum_name: "Militar-Geographische Institut",
369};
370
371/// Irish National Datum
372pub const OSNI52: ToWGS84Datum = ToWGS84Datum {
373    datum_params: DatumParams::Param7(482.53, -130.596, 564.557, -1.042, -0.214, -0.631, 8.15),
374    ellipse: "airy",
375    datum_name: "Irish National",
376};
377
378/// Ireland 1965 Datum
379pub const IRE65: ToWGS84Datum = ToWGS84Datum {
380    datum_params: DatumParams::Param7(482.53, -130.596, 564.557, -1.042, -0.214, -0.631, 8.15),
381    ellipse: "mod_airy",
382    datum_name: "Ireland 1965",
383};
384
385/// Rassadiran Datum
386pub const RASSADIRAN: ToWGS84Datum = ToWGS84Datum {
387    datum_params: DatumParams::Param3(-133.63, -157.5, -158.62),
388    ellipse: "intl",
389    datum_name: "Rassadiran",
390};
391
392/// New Zealand Geodetic Datum 1949 Datum
393pub const NZGD49: ToWGS84Datum = ToWGS84Datum {
394    datum_params: DatumParams::Param7(59.47, -5.04, 187.44, 0.47, -0.1, 1.024, -4.5993),
395    ellipse: "intl",
396    datum_name: "New Zealand Geodetic Datum 1949",
397};
398
399/// Airy 1830 Datum
400pub const OSGB36: ToWGS84Datum = ToWGS84Datum {
401    datum_params: DatumParams::Param7(446.448, -125.157, 542.06, 0.1502, 0.247, 0.8421, -20.4894),
402    ellipse: "airy",
403    datum_name: "Airy 1830",
404};
405
406/// S-JTSK (Ferro) Datum
407pub const S_JTSK: ToWGS84Datum = ToWGS84Datum {
408    datum_params: DatumParams::Param3(589.0, 76.0, 480.0),
409    ellipse: "bessel",
410    datum_name: "S-JTSK (Ferro)",
411};
412
413/// Beduaram Datum
414pub const BEDUARAM: ToWGS84Datum = ToWGS84Datum {
415    datum_params: DatumParams::Param3(-106.0, -87.0, 188.0),
416    ellipse: "clrk80",
417    datum_name: "Beduaram",
418};
419
420/// Gunung Segara Jakarta Datum
421pub const GUNUNG_SEGARA: ToWGS84Datum = ToWGS84Datum {
422    datum_params: DatumParams::Param3(-403.0, 684.0, 41.0),
423    ellipse: "bessel",
424    datum_name: "Gunung Segara Jakarta",
425};
426
427/// Reseau National Belge 1972
428pub const RNB72: ToWGS84Datum = ToWGS84Datum {
429    datum_params: DatumParams::Param7(
430        106.869, -52.2978, 103.724, -0.33657, 0.456955, -1.84218, 1.0,
431    ),
432    ellipse: "intl",
433    datum_name: "Reseau National Belge 1972",
434};
435
436/// Given a name, return the corresponding ellipsoid
437#[cfg_attr(feature = "nightly", coverage(off))]
438pub fn get_datum(name: &str) -> Option<ToWGS84Datum> {
439    // fix name to remove _ and convert to uppercase
440    let name = name.to_uppercase().replace("_", "");
441    match name.as_str() {
442        "WGS84" => Some(TO_WGS84),
443        "CH1903" => Some(CH1903),
444        "GGRS87" => Some(GGRS87),
445        "NAD83" => Some(NAD83),
446        "RASSADIRAN" => Some(RASSADIRAN),
447        "NZGD49" => Some(NZGD49),
448        "OSGB36" => Some(OSGB36),
449        "SJTSK" => Some(S_JTSK),
450        "BEDUARAM" => Some(BEDUARAM),
451        "POTSDAM" => Some(POTSDAM),
452        "CARTHAGE" => Some(CARTHAGE),
453        "HERMANNSKOGEL" => Some(HERMANNSKOGEL),
454        "MILITARGEOGRAPHISCHEINSTITUT" => Some(MILITARGEOGRAPHISCHE_INSTITUT),
455        "OSNI52" => Some(OSNI52),
456        "IRE65" => Some(IRE65),
457        "RNB72" => Some(RNB72),
458        "GUNUNG_SEGARA" => Some(GUNUNG_SEGARA),
459        _ => None,
460    }
461}
gistools/proj/parse/datum.rs

gistools/proj/parse/
datum.rs
