proj_core/
datum.rs

1use crate::ellipsoid::{self, Ellipsoid};
2use crate::error::{Error, Result};
3use crate::grid::GridDefinition;
4use smallvec::SmallVec;
5
6/// A geodetic datum, defined by a reference ellipsoid and its relationship to WGS84.
7#[derive(Debug, Clone)]
8pub struct Datum {
9    /// The reference ellipsoid.
10    ellipsoid: Ellipsoid,
11    /// Explicit relationship from this datum to WGS84.
12    to_wgs84: DatumToWgs84,
13}
14
15impl Datum {
16    /// Create a datum from an ellipsoid and an explicit path to WGS84.
17    pub fn new(ellipsoid: Ellipsoid, to_wgs84: DatumToWgs84) -> Result<Self> {
18        to_wgs84.validate()?;
19        Ok(Self {
20            ellipsoid,
21            to_wgs84,
22        })
23    }
24
25    const fn new_unchecked(ellipsoid: Ellipsoid, to_wgs84: DatumToWgs84) -> Self {
26        Self {
27            ellipsoid,
28            to_wgs84,
29        }
30    }
31
32    /// Return the reference ellipsoid.
33    pub const fn ellipsoid(&self) -> Ellipsoid {
34        self.ellipsoid
35    }
36
37    /// Return the explicit relationship from this datum to WGS84.
38    pub const fn to_wgs84(&self) -> &DatumToWgs84 {
39        &self.to_wgs84
40    }
41
42    /// Returns true if this datum is WGS84 or functionally identical (no Helmert shift needed).
43    pub fn is_wgs84_compatible(&self) -> bool {
44        matches!(self.to_wgs84, DatumToWgs84::Identity)
45    }
46
47    /// Returns true if this datum has a known path to WGS84.
48    pub fn has_known_wgs84_transform(&self) -> bool {
49        !matches!(self.to_wgs84, DatumToWgs84::Unknown)
50    }
51
52    /// Returns true if this datum's WGS84 path uses one or more horizontal grids.
53    pub fn uses_grid_shift(&self) -> bool {
54        self.to_wgs84.uses_grid_shift()
55    }
56
57    /// Return the Helmert parameters for this datum's path to WGS84, when available.
58    pub fn helmert_to_wgs84(&self) -> Option<&HelmertParams> {
59        match &self.to_wgs84 {
60            DatumToWgs84::Helmert(params) => Some(params),
61            DatumToWgs84::Identity | DatumToWgs84::GridShift(_) | DatumToWgs84::Unknown => None,
62        }
63    }
64
65    pub fn approximate_helmert_to(&self, target: &Datum) -> Option<HelmertParams> {
66        if self.uses_grid_shift() || target.uses_grid_shift() {
67            return None;
68        }
69        if !self.has_known_wgs84_transform() || !target.has_known_wgs84_transform() {
70            return None;
71        }
72
73        let source_to_wgs84 = self
74            .helmert_to_wgs84()
75            .copied()
76            .unwrap_or_else(HelmertParams::identity);
77        let wgs84_to_target = target
78            .helmert_to_wgs84()
79            .map(HelmertParams::inverse)
80            .unwrap_or_else(HelmertParams::identity);
81
82        source_to_wgs84.compose_approx(&wgs84_to_target).ok()
83    }
84
85    /// Returns true if two datums are the same (same ellipsoid, same Helmert parameters).
86    pub fn same_datum(&self, other: &Datum) -> bool {
87        let same_ellipsoid =
88            (self.ellipsoid.semi_major_axis() - other.ellipsoid.semi_major_axis()).abs() < 1e-6
89                && (self.ellipsoid.flattening() - other.ellipsoid.flattening()).abs() < 1e-12;
90
91        match (&self.to_wgs84, &other.to_wgs84) {
92            (DatumToWgs84::Identity, DatumToWgs84::Identity) => same_ellipsoid,
93            (DatumToWgs84::Helmert(a), DatumToWgs84::Helmert(b)) => {
94                same_ellipsoid && a.approx_eq(b)
95            }
96            (DatumToWgs84::GridShift(a), DatumToWgs84::GridShift(b)) => same_ellipsoid && a == b,
97            (DatumToWgs84::Unknown, DatumToWgs84::Unknown) => false,
98            _ => false,
99        }
100    }
101}
102
103/// Explicit WGS84 relationship for a datum.
104#[derive(Debug, Clone, PartialEq)]
105pub enum DatumToWgs84 {
106    /// The datum can be treated as WGS84-compatible in the current model.
107    Identity,
108    /// The datum requires the provided Helmert transform to reach WGS84.
109    Helmert(HelmertParams),
110    /// The datum requires horizontal grid interpolation to reach WGS84.
111    GridShift(Box<DatumGridShift>),
112    /// The datum's path to WGS84 is not known.
113    Unknown,
114}
115
116impl DatumToWgs84 {
117    pub fn uses_grid_shift(&self) -> bool {
118        matches!(self, DatumToWgs84::GridShift(shift) if shift.uses_grid_shift())
119    }
120
121    pub fn validate(&self) -> Result<()> {
122        match self {
123            Self::Helmert(params) => params.validate(),
124            Self::Identity | Self::GridShift(_) | Self::Unknown => Ok(()),
125        }
126    }
127}
128
129/// Ordered PROJ-style datum grid list.
130#[derive(Debug, Clone, PartialEq)]
131pub struct DatumGridShift {
132    entries: SmallVec<[DatumGridShiftEntry; 4]>,
133}
134
135impl DatumGridShift {
136    pub fn new(entries: SmallVec<[DatumGridShiftEntry; 4]>) -> Self {
137        Self { entries }
138    }
139
140    pub fn from_vec(entries: Vec<DatumGridShiftEntry>) -> Self {
141        Self {
142            entries: SmallVec::from_vec(entries),
143        }
144    }
145
146    pub fn entries(&self) -> &[DatumGridShiftEntry] {
147        &self.entries
148    }
149
150    pub fn uses_grid_shift(&self) -> bool {
151        self.entries
152            .iter()
153            .any(|entry| matches!(entry, DatumGridShiftEntry::Grid { .. }))
154    }
155}
156
157/// One entry from a datum grid list.
158#[derive(Debug, Clone, PartialEq)]
159pub enum DatumGridShiftEntry {
160    /// Try this horizontal grid. Optional grids may be missing from providers.
161    Grid {
162        definition: GridDefinition,
163        optional: bool,
164    },
165    /// PROJ's `null` grid: stop grid lookup and apply no shift.
166    Null,
167}
168
169/// 7-parameter Helmert (Bursa-Wolf) transformation parameters.
170///
171/// Defines the transformation from one datum to WGS84 geocentric coordinates:
172/// ```text
173/// [X']   [dx]         [  1  -rz   ry] [X]
174/// [Y'] = [dy] + (1+ds)[  rz   1  -rx] [Y]
175/// [Z']   [dz]         [ -ry  rx   1 ] [Z]
176/// ```
177#[derive(Debug, Clone, Copy, PartialEq)]
178pub struct HelmertParams {
179    /// X-axis translation in meters.
180    dx: f64,
181    /// Y-axis translation in meters.
182    dy: f64,
183    /// Z-axis translation in meters.
184    dz: f64,
185    /// X-axis rotation in arc-seconds.
186    rx: f64,
187    /// Y-axis rotation in arc-seconds.
188    ry: f64,
189    /// Z-axis rotation in arc-seconds.
190    rz: f64,
191    /// Scale difference in parts-per-million (ppm).
192    ds: f64,
193}
194
195impl HelmertParams {
196    /// Create a 7-parameter transformation.
197    pub fn new(dx: f64, dy: f64, dz: f64, rx: f64, ry: f64, rz: f64, ds: f64) -> Result<Self> {
198        let params = Self::new_unchecked(dx, dy, dz, rx, ry, rz, ds);
199        params.validate()?;
200        Ok(params)
201    }
202
203    /// Create a translation-only (3-parameter) transformation.
204    pub fn translation(dx: f64, dy: f64, dz: f64) -> Result<Self> {
205        Self::new(dx, dy, dz, 0.0, 0.0, 0.0, 0.0)
206    }
207
208    /// Identity transformation.
209    pub const fn identity() -> Self {
210        Self::new_unchecked(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
211    }
212
213    const fn translation_unchecked(dx: f64, dy: f64, dz: f64) -> Self {
214        Self::new_unchecked(dx, dy, dz, 0.0, 0.0, 0.0, 0.0)
215    }
216
217    const fn new_unchecked(dx: f64, dy: f64, dz: f64, rx: f64, ry: f64, rz: f64, ds: f64) -> Self {
218        Self {
219            dx,
220            dy,
221            dz,
222            rx,
223            ry,
224            rz,
225            ds,
226        }
227    }
228
229    pub const fn dx(&self) -> f64 {
230        self.dx
231    }
232
233    pub const fn dy(&self) -> f64 {
234        self.dy
235    }
236
237    pub const fn dz(&self) -> f64 {
238        self.dz
239    }
240
241    pub const fn rx(&self) -> f64 {
242        self.rx
243    }
244
245    pub const fn ry(&self) -> f64 {
246        self.ry
247    }
248
249    pub const fn rz(&self) -> f64 {
250        self.rz
251    }
252
253    pub const fn ds(&self) -> f64 {
254        self.ds
255    }
256
257    pub fn validate(&self) -> Result<()> {
258        if self.dx.is_finite()
259            && self.dy.is_finite()
260            && self.dz.is_finite()
261            && self.rx.is_finite()
262            && self.ry.is_finite()
263            && self.rz.is_finite()
264            && self.ds.is_finite()
265        {
266            return Ok(());
267        }
268
269        Err(Error::InvalidDefinition(
270            "Helmert parameters must be finite".into(),
271        ))
272    }
273
274    /// Return the inverse parameters (WGS84 → this datum).
275    pub fn inverse(&self) -> Self {
276        Self {
277            dx: -self.dx,
278            dy: -self.dy,
279            dz: -self.dz,
280            rx: -self.rx,
281            ry: -self.ry,
282            rz: -self.rz,
283            ds: -self.ds,
284        }
285    }
286
287    pub fn compose_approx(&self, next: &Self) -> Result<Self> {
288        Self::new(
289            self.dx + next.dx,
290            self.dy + next.dy,
291            self.dz + next.dz,
292            self.rx + next.rx,
293            self.ry + next.ry,
294            self.rz + next.rz,
295            self.ds + next.ds,
296        )
297    }
298
299    fn approx_eq(&self, other: &Self) -> bool {
300        (self.dx - other.dx).abs() < 1e-6
301            && (self.dy - other.dy).abs() < 1e-6
302            && (self.dz - other.dz).abs() < 1e-6
303            && (self.rx - other.rx).abs() < 1e-9
304            && (self.ry - other.ry).abs() < 1e-9
305            && (self.rz - other.rz).abs() < 1e-9
306            && (self.ds - other.ds).abs() < 1e-9
307    }
308}
309
310// ---------------------------------------------------------------------------
311// Well-known datums
312// ---------------------------------------------------------------------------
313
314/// WGS 84 datum.
315pub const WGS84: Datum = Datum::new_unchecked(ellipsoid::WGS84, DatumToWgs84::Identity);
316
317/// NAD83 datum (functionally identical to WGS84 for sub-meter work).
318pub const NAD83: Datum = Datum::new_unchecked(ellipsoid::GRS80, DatumToWgs84::Identity);
319
320/// NAD27 datum (Clarke 1866 ellipsoid).
321/// Helmert parameters from EPSG dataset (approximate continental US average).
322pub const NAD27: Datum = Datum::new_unchecked(
323    ellipsoid::CLARKE1866,
324    DatumToWgs84::Helmert(HelmertParams::translation_unchecked(-8.0, 160.0, 176.0)),
325);
326
327/// ETRS89 datum (European Terrestrial Reference System 1989).
328/// Functionally identical to WGS84 for most purposes.
329pub const ETRS89: Datum = Datum::new_unchecked(ellipsoid::GRS80, DatumToWgs84::Identity);
330
331/// OSGB36 datum (Ordnance Survey Great Britain 1936).
332pub const OSGB36: Datum = Datum::new_unchecked(
333    ellipsoid::AIRY1830,
334    DatumToWgs84::Helmert(HelmertParams::new_unchecked(
335        446.448, -125.157, 542.060, 0.1502, 0.2470, 0.8421, -20.4894,
336    )),
337);
338
339/// Pulkovo 1942 datum (used in Russia and former Soviet states).
340pub const PULKOVO1942: Datum = Datum::new_unchecked(
341    ellipsoid::KRASSOWSKY,
342    DatumToWgs84::Helmert(HelmertParams::translation_unchecked(23.92, -141.27, -80.9)),
343);
344
345/// ED50 datum (European Datum 1950).
346pub const ED50: Datum = Datum::new_unchecked(
347    ellipsoid::INTL1924,
348    DatumToWgs84::Helmert(HelmertParams::translation_unchecked(-87.0, -98.0, -121.0)),
349);
350
351/// Tokyo datum (used in Japan).
352pub const TOKYO: Datum = Datum::new_unchecked(
353    ellipsoid::BESSEL1841,
354    DatumToWgs84::Helmert(HelmertParams::translation_unchecked(
355        -146.414, 507.337, 680.507,
356    )),
357);
358
359#[cfg(test)]
360mod tests {
361    use super::*;
362
363    #[test]
364    fn wgs84_is_wgs84_compatible() {
365        assert!(WGS84.is_wgs84_compatible());
366        assert!(NAD83.is_wgs84_compatible());
367        assert!(ETRS89.is_wgs84_compatible());
368    }
369
370    #[test]
371    fn nad27_is_not_wgs84_compatible() {
372        assert!(!NAD27.is_wgs84_compatible());
373        assert!(!OSGB36.is_wgs84_compatible());
374    }
375
376    #[test]
377    fn same_datum_identity() {
378        assert!(WGS84.same_datum(&WGS84));
379        assert!(NAD27.same_datum(&NAD27));
380    }
381
382    #[test]
383    fn different_datums() {
384        assert!(!WGS84.same_datum(&NAD27));
385        assert!(!NAD27.same_datum(&OSGB36));
386    }
387
388    #[test]
389    fn unknown_datums_are_not_collapsed_by_ellipsoid() {
390        let a = Datum::new(ellipsoid::WGS84, DatumToWgs84::Unknown).unwrap();
391        let b = Datum::new(ellipsoid::WGS84, DatumToWgs84::Unknown).unwrap();
392
393        assert!(!a.same_datum(&b));
394    }
395
396    #[test]
397    fn helmert_inverse_negates() {
398        let h = HelmertParams::new(1.0, 2.0, 3.0, 0.1, 0.2, 0.3, 0.5).unwrap();
399        let inv = h.inverse();
400        assert_eq!(inv.dx(), -1.0);
401        assert_eq!(inv.rx(), -0.1);
402        assert_eq!(inv.ds(), -0.5);
403    }
404
405    #[test]
406    fn helmert_params_reject_non_finite_values() {
407        let err = HelmertParams::new(f64::NAN, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0).unwrap_err();
408        assert!(matches!(err, Error::InvalidDefinition(_)), "got {err}");
409
410        let err = HelmertParams::new(0.0, 0.0, 0.0, 0.0, 0.0, f64::INFINITY, 0.0).unwrap_err();
411        assert!(err.to_string().contains("Helmert parameters"), "{err}");
412    }
413}
proj_core/datum.rs

proj_core/
datum.rs
