Skip to main content

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