Coverage for pygeodesy/datums.py: 94%
241 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-02-12 13:15 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2024-02-12 13:15 -0500
2# -*- coding: utf-8 -*-
4u'''Datums and transformations thereof.
6Classes L{Datum} and L{Transform} and registries L{Datums} and L{Transforms}, respectively.
8Pure Python implementation of geodesy tools for ellipsoidal earth models, including datums
9and ellipsoid parameters for different geographic coordinate systems and methods for
10converting between them and to cartesian coordinates. Transcoded from JavaScript originals by
11I{(C) Chris Veness 2005-2016} and published under the same MIT Licence**, see U{latlon-ellipsoidal.js
12<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
14Historical geodetic datums: a latitude/longitude point defines a geographic location on, above
15or below the earth’s surface. Latitude is measured in degrees from the equator, lomgitude from
16the International Reference Meridian and height in meters above an ellipsoid based on the given
17datum. The datum in turn is based on a reference ellipsoid and tied to geodetic survey
18reference points.
20Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems), but
21previously various other reference ellipsoids and datum references were used.
23The UK Ordnance Survey National Grid References are still based on the otherwise historical OSGB36
24datum, q.v. U{"A Guide to Coordinate Systems in Great Britain", Section 6
25<https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>}.
27@var Datums.BD72: Datum(name='BD72', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.BD72)
28@var Datums.DHDN: Datum(name='DHDN', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.DHDN)
29@var Datums.ED50: Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50)
30@var Datums.GDA2020: Datum(name='GDA2020', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)
31@var Datums.GRS80: Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)
32@var Datums.Irl1975: Datum(name='Irl1975', ellipsoid=Ellipsoids.AiryModified, transform=Transforms.Irl1975)
33@var Datums.Krassovski1940: Datum(name='Krassovski1940', ellipsoid=Ellipsoids.Krassovski1940, transform=Transforms.Krassovski1940)
34@var Datums.Krassowsky1940: Datum(name='Krassowsky1940', ellipsoid=Ellipsoids.Krassowsky1940, transform=Transforms.Krassowsky1940)
35@var Datums.MGI: Datum(name='MGI', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.MGI)
36@var Datums.NAD27: Datum(name='NAD27', ellipsoid=Ellipsoids.Clarke1866, transform=Transforms.NAD27)
37@var Datums.NAD83: Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83)
38@var Datums.NTF: Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF)
39@var Datums.OSGB36: Datum(name='OSGB36', ellipsoid=Ellipsoids.Airy1830, transform=Transforms.OSGB36)
40@var Datums.Potsdam: Datum(name='Potsdam', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841)
41@var Datums.Sphere: Datum(name='Sphere', ellipsoid=Ellipsoids.Sphere, transform=Transforms.WGS84)
42@var Datums.TokyoJapan: Datum(name='TokyoJapan', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.TokyoJapan)
43@var Datums.WGS72: Datum(name='WGS72', ellipsoid=Ellipsoids.WGS72, transform=Transforms.WGS72)
44@var Datums.WGS84: Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84)
46@var Transforms.BD72: Transform(name='BD72', tx=106.869, ty=-52.2978, tz=103.724, s1=1.0, rx=-1.63174e-06, ry=-2.21538e-06, rz=-8.93114e-06, s=1.2727, sx=-0.33657, sy=-0.456955, sz=-1.84218)
47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, s1=0.999992, rx=-5.04206e-06, ry=-1.69685e-06, rz=1.49323e-05, s=-8.3, sx=-1.04, sy=-0.35, sz=3.08)
48@var Transforms.Clarke1866: Transform(name='Clarke1866', tx=8.0, ty=-160, tz=-176, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
49@var Transforms.DHDN: Transform(name='DHDN', tx=-591.28, ty=-81.35, tz=-396.39, s1=0.99999, rx=7.1607e-06, ry=-3.56823e-07, rz=-7.06858e-06, s=-9.82, sx=1.477, sy=-0.0736, sz=-1.458)
50@var Transforms.DHDNE: Transform(name='DHDNE', tx=-612.4, ty=-77, tz=-440.2, s1=0.999997, rx=2.61799e-07, ry=-2.76344e-07, rz=1.35602e-05, s=-2.55, sx=0.054, sy=-0.057, sz=2.797)
51@var Transforms.DHDNW: Transform(name='DHDNW', tx=-598.1, ty=-73.7, tz=-418.2, s1=0.999993, rx=-9.79324e-07, ry=-2.18166e-07, rz=1.19022e-05, s=-6.7, sx=-0.202, sy=-0.045, sz=2.455)
52@var Transforms.ED50: Transform(name='ED50', tx=89.5, ty=93.8, tz=123.1, s1=0.999999, rx=0.0, ry=0.0, rz=7.56309e-07, s=-1.2, sx=0.0, sy=0.0, sz=0.156)
53@var Transforms.Identity: Transform(name='Identity', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
54@var Transforms.Irl1965: Transform(name='Irl1965', tx=-482.53, ty=130.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-06, s=-8.15, sx=1.042, sy=0.214, sz=0.631)
55@var Transforms.Irl1975: Transform(name='Irl1975', tx=-482.53, ty=130.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-06, s=-8.15, sx=1.042, sy=0.214, sz=0.631)
56@var Transforms.Krassovski1940: Transform(name='Krassovski1940', tx=-24, ty=123.0, tz=94.0, s1=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)
57@var Transforms.Krassowsky1940: Transform(name='Krassowsky1940', tx=-24, ty=123.0, tz=94.0, s1=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)
58@var Transforms.MGI: Transform(name='MGI', tx=-577.326, ty=-90.129, tz=-463.92, s1=0.999998, rx=2.49049e-05, ry=7.14615e-06, rz=2.56806e-05, s=-2.423, sx=5.137, sy=1.474, sz=5.297)
59@var Transforms.NAD27: Transform(name='NAD27', tx=8.0, ty=-160, tz=-176, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
60@var Transforms.NAD83: Transform(name='NAD83', tx=1.004, ty=-1.91, tz=-0.515, s1=1.0, rx=1.29445e-07, ry=1.64837e-09, rz=5.33295e-08, s=-0.0015, sx=0.0267, sy=0.00034, sz=0.011)
61@var Transforms.NTF: Transform(name='NTF', tx=-168, ty=-60, tz=320.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
62@var Transforms.OSGB36: Transform(name='OSGB36', tx=-446.448, ty=125.157, tz=-542.06, s1=1.00002, rx=-7.2819e-07, ry=-1.19749e-06, rz=-4.08262e-06, s=20.4894, sx=-0.1502, sy=-0.247, sz=-0.8421)
63@var Transforms.TokyoJapan: Transform(name='TokyoJapan', tx=148.0, ty=-507, tz=-685, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
64@var Transforms.WGS72: Transform(name='WGS72', tx=0.0, ty=0.0, tz=-4.5, s1=1.0, rx=0.0, ry=0.0, rz=2.68587e-06, s=-0.22, sx=0.0, sy=0.0, sz=0.554)
65@var Transforms.WGS84: Transform(name='WGS84', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)
66'''
67# make sure int/int division yields float quotient, see .basics
68from __future__ import division as _; del _ # PYCHOK semicolon
70from pygeodesy.basics import islistuple, map2, neg, _xinstanceof
71from pygeodesy.constants import R_M, _float as _F, _0_0, _0_26, _1_0, _2_0, _8_0, _3600_0
72# from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase as _CEB, \
73# LatLonEllipsoidalBase as _LLEB # MODS
74from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, _EWGS84, \
75 Vector3Tuple
76from pygeodesy.errors import _IsnotError, _TypeError, _xattr
77from pygeodesy.fmath import fdot, fmean, Fmt
78from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, \
79 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \
80 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \
81 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \
82 _NAD83_, _s_, _Sphere_, _spherical_, _transform_, \
83 _UNDER_, _WGS72_, _WGS84_, _under
84from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
85from pygeodesy.named import _NamedEnum, _NamedEnumItem, _lazyNamedEnumItem as _lazy
86# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
87from pygeodesy.props import Property_RO, property_RO
88# from pygeodesy.streprs import Fmt # from .fmath
89from pygeodesy.units import _isRadius, Radius_, radians
91# from math import radians # from .units
93__all__ = _ALL_LAZY.datums
94__version__ = '24.02.12'
96_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
97_BD72_ = 'BD72'
98_DHDN_ = 'DHDN'
99_DHDNE_ = 'DHDNE'
100_DHDNW_ = 'DHDNW'
101_ED50_ = 'ED50'
102_GDA2020_ = 'GDA2020'
103_Identity_ = 'Identity'
104_Irl1965_ = 'Irl1965'
105_Irl1975_ = 'Irl1975'
106_MGI_ = 'MGI'
107_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf
108_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:]
109_NTF_ = 'NTF'
110_OSGB36_ = 'OSGB36'
111_Potsdam_ = 'Potsdam'
112_RPS = radians(_1_0 / _3600_0) # radians per degree-second
113_S1_S = 1.e-6 # in .trf
114_TokyoJapan_ = 'TokyoJapan'
117class Transform(_NamedEnumItem):
118 '''Helmert I{datum} transformation.
120 @see: L{TransformXform<trf.TransformXform>}.
121 '''
122 tx = _0_0 # x translation (C{meter})
123 ty = _0_0 # y translation (C{meter})
124 tz = _0_0 # z translation (C{meter})
126 rx = _0_0 # x rotation (C{radians})
127 ry = _0_0 # y rotation (C{radians})
128 rz = _0_0 # z rotation (C{radians})
130 s = _0_0 # scale ppm (C{float})
131 s1 = _1_0 # scale + 1 (C{float})
133 sx = _0_0 # x rotation (C{degree-seconds})
134 sy = _0_0 # y rotation (C{degree-seconds})
135 sz = _0_0 # z rotation (C{degree-seconds})
137 def __init__(self, name=NN, tx=0, ty=0, tz=0,
138 sx=0, sy=0, sz=0, s=0):
139 '''New L{Transform}.
141 @kwarg name: Optional, unique name (C{str}).
142 @kwarg tx: X translation (C{meter}).
143 @kwarg ty: Y translation (C{meter}).
144 @kwarg tz: Z translation (C{meter}).
145 @kwarg s: Scale (C{float}), ppm.
146 @kwarg sx: X rotation (C{degree-seconds}).
147 @kwarg sy: Y rotation (C{degree-seconds}).
148 @kwarg sz: Z rotation (C{degree-seconds}).
150 @raise NameError: Transform with that B{C{name}} already exists.
151 '''
152 if tx:
153 self.tx = tx
154 if ty:
155 self.ty = ty
156 if tz:
157 self.tz = tz
158 if sx: # secs to rads
159 self.rx, self.sx = self._rps2(sx)
160 if sy:
161 self.ry, self.sy = self._rps2(sy)
162 if sz:
163 self.rz, self.sz = self._rps2(sz)
164 if s:
165 self.s = s
166 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1)
168 self._register(Transforms, name)
170 def __eq__(self, other):
171 '''Compare this and an other transform.
173 @arg other: The other transform (L{Transform}).
175 @return: C{True} if equal, C{False} otherwise.
176 '''
177 return self is other or (isinstance(other, Transform) and all(
178 a == b for a, b in zip(self, other)))
180 def __hash__(self):
181 return self._hash # memoized
183 def __iter__(self):
184 '''Yield the initial attribute values.
185 '''
186 for _, x in self.items():
187 yield x
189 def __matmul__(self, point): # PYCHOK Python 3.5+
190 '''Helmert-transform an I{ellipsoidal} cartesian B{C{point}}.
192 @raise TypeError: Invalid B{C{point}}.
193 '''
194 m = _MODS.ellipsoidalBase
195 _xinstanceof(m.CartesianEllipsoidalBase, point=point)
196 return point.toTransform(self)
198 def __neg__(self):
199 return self.inverse()
201# def __sub__(self, other):
202# _xinstanceof(Transform, other=other)
203#
204# def _sub(a, b):
205# for n in _Names11:
206# yield getattr(a, n) - getattr(b, n)
207#
208# return type(self)(_sub(self, other), name=_MINUS_) # .fsums._sub_op
210 @Property_RO
211 def _hash(self):
212 return hash(x for x in self)
214 def items(self):
215 '''Yield each attribute as 2-tuple C{(name, value)}.
216 '''
217 for n in _Names7:
218 yield n, getattr(self, n)
220 def inverse(self, name=NN):
221 '''Return the inverse of this transform.
223 @kwarg name: Optional, unique name (C{str}).
225 @return: Inverse (L{Transform}).
227 @raise NameError: Transform with that B{C{name}} already exists.
228 '''
229 d = dict((n, neg(v)) for n, v in self.items())
230 n = name or _minus(self.name)
231 return type(self)(name=n, **d)
233 @Property_RO
234 def isunity(self):
235 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84?
236 '''
237 return not any(s for s in self)
239 def _rps2(self, s_):
240 '''(INTERNAL) Rotation in C{radians} and C{degree-seconds}.
241 '''
242 return (_RPS * s_), s_
244 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected
245 '''Return this transform as a string.
247 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
248 @kwarg fmt: Optional C{float} format (C{letter}).
249 @kwarg name: Override name (C{str}) or C{None} to exclude
250 this transform's name.
252 @return: Transform attributes (C{str}).
253 '''
254 return self._instr(name, prec, *_Names11, fmt=fmt)
256 def transform(self, x, y, z, inverse=False):
257 '''Transform a (geocentric) position, forward or inverse.
259 @arg x: X coordinate (C{meter}).
260 @arg y: Y coordinate (C{meter}).
261 @arg z: Z coordinate (C{meter}).
262 @kwarg inverse: Optional direction, forward or inverse (C{bool}).
264 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)}).
265 '''
266 if self.isunity:
267 return Vector3Tuple(x, y, z, name=self.name)
269 xyz1 = x, y, z, _1_0
270 s1 = self.s1
271 if inverse:
272 xyz1 = map2(neg, xyz1)
273 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1
274 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
275 # x * .rz + y * .s1 - z * .rx + .ty,
276 # -x * .ry + y * .rx + z * .s1 + .tz)
277 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
278 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
279 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
280 name=self.name)
283class Transforms(_NamedEnum):
284 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
285 to accommodate the L{_LazyNamedEnumItem} properties.
286 '''
287 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
288 '''(INTERNAL) Instantiate the C{Transform}.
289 '''
290 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
292Transforms = Transforms(Transform) # PYCHOK singleton
293'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
294# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
295Transforms._assert(
296 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
297 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
298 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
299 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
300 s=_F( 1.2727)),
301 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
302 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
303 s=_F( -8.3)),
304 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
305 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
306 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
307 s=_F( -9.82)), # Germany
308 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2),
309 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
310 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797),
311 s=_F( -2.55)), # East Germany
312 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2),
313 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
314 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455),
315 s=_F( -6.7)), # West Germany
316 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
317 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
318 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
319 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
320 sz=_F( 0.156), s=_F(-1.2)),
321 Identity = _lazy(_Identity_),
322 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
323 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
324 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
325 s=_F( -8.15)),
326 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
327 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
328 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
329 s=_F( -8.15)),
330 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
331 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
332 s=_F( -2.423)), # spelling
333 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
334 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
335 s=_F( -2.423)), # spelling
336 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
337 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
338 s=_F( -2.423)), # Austria
339 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
340 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
341 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
342 s=_F(-0.0015)),
343 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
344 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
345 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
346 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
347 s=_F( 20.4894)),
348 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
349 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
350 WGS84 = _lazy(_WGS84_), # unity
351)
354class Datum(_NamedEnumItem):
355 '''Ellipsoid and transform parameters for an earth model.
356 '''
357 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
358 _transform = Transforms.WGS84 # default transform (L{Transform})
360 def __init__(self, ellipsoid, transform=None, name=NN):
361 '''New L{Datum}.
363 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
364 @kwarg transform: Optional transform (L{Transform}).
365 @kwarg name: Optional, unique name (C{str}).
367 @raise NameError: Datum with that B{C{name}} already exists.
369 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
370 nor L{Ellipsoid2} or B{C{transform}} is
371 not a L{Transform}.
372 '''
373 self._ellipsoid = ellipsoid or Datum._ellipsoid
374 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
376 self._transform = transform or Datum._transform
377 _xinstanceof(Transform, transform=self.transform)
379 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
381 def __eq__(self, other):
382 '''Compare this and an other datum.
384 @arg other: The other datum (L{Datum}).
386 @return: C{True} if equal, C{False} otherwise.
387 '''
388 return self is other or (isinstance(other, Datum) and
389 self.ellipsoid == other.ellipsoid and
390 self.transform == other.transform)
392 def __hash__(self):
393 return self._hash # memoized
395 def __matmul__(self, point): # PYCHOK Python 3.5+
396 '''Convert an I{ellipsoidal} B{C{point}} to this datum.
398 @raise TypeError: Invalid B{C{point}}.
399 '''
400 m = _MODS.ellipsoidalBase
401 _xinstanceof(m.CartesianEllipsoidalBase, m.LatLonEllipsoidalBase, point=point)
402 return point.toDatum(self)
404 def ecef(self, Ecef=None):
405 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
407 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
409 @return: An ECEF converter for this C{datum}.
411 @raise TypeError: Invalid B{C{Ecef}}.
413 @see: Module L{pygeodesy.ecef}.
414 '''
415 return _MODS.ecef._4Ecef(self, Ecef)
417 @Property_RO
418 def ellipsoid(self):
419 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
420 '''
421 return self._ellipsoid
423 @Property_RO
424 def exactTM(self):
425 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
426 '''
427 return _MODS.etm.ExactTransverseMercator(datum=self)
429 @Property_RO
430 def _hash(self):
431 return hash(self.ellipsoid) + hash(self.transform)
433 @property_RO
434 def isEllipsoidal(self):
435 '''Check whether this datum is ellipsoidal (C{bool}).
436 '''
437 return self.ellipsoid.isEllipsoidal
439 @property_RO
440 def isOblate(self):
441 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
442 '''
443 return self.ellipsoid.isOblate
445 @property_RO
446 def isProlate(self):
447 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
448 '''
449 return self.ellipsoid.isProlate
451 @property_RO
452 def isSpherical(self):
453 '''Check whether this datum is (near-)spherical (C{bool}).
454 '''
455 return self.ellipsoid.isSpherical
457 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
458 '''Return this datum as a string.
460 @kwarg sep: Separator to join (C{str}).
461 @kwarg name: Override name (C{str}) or C{None} to exclude
462 this datum's name.
464 @return: Datum attributes (C{str}).
465 '''
466 t = [] if name is None else \
467 [Fmt.EQUAL(name=repr(name or self.named))]
468 for a in (_ellipsoid_, _transform_):
469 v = getattr(self, a)
470 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
471 return sep.join(t)
473 @Property_RO
474 def transform(self):
475 '''Get this datum's transform (L{Transform}).
476 '''
477 return self._transform
480def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
481 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
482 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
484 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
485 '''
486 if f is not None:
487 E, n, D = _EnD3((a_earth, f), name)
488 if raiser and not E:
489 raise _TypeError(f=f, **{raiser: a_earth})
490 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None):
491 return
492 elif isinstance(a_earth, Datum):
493 E, n, D = None, NN, a_earth
494 else:
495 E, n, D = _EnD3(a_earth, name)
496 if raiser and not E:
497 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
498 if D is None:
499 D = Datum(E, transform=Transforms.Identity, name=_under(n))
500 inst._datum = D
503def _earth_ellipsoid(earth, *name_raiser):
504 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
505 '''
506 return Ellipsoids.Sphere if earth is R_M else (
507 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
508 _spherical_datum(earth, *name_raiser).ellipsoid)
511def _ED2(radius, name):
512 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
513 '''
514 D = Datums.Sphere
515 E = D.ellipsoid
516 if name or radius != E.a: # != E.b
517 n = _under(name)
518 E = Ellipsoid(radius, radius, name=n)
519 D = Datum(E, transform=Transforms.Identity, name=n)
520 return E, D
523def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
524 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
525 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
527 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
528 '''
529 if isinstance(earth, Datum):
530 D = earth
531 else:
532 E, n, D = _EnD3(earth, name)
533 if not E:
534 n = raiser or _earth_
535 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
536 if D is None:
537 D = Datum(E, transform=Transforms.Identity, name=_under(n))
538 if raiser and not D.isEllipsoidal:
539 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
540 return D
543def _EnD3(earth, name):
544 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
545 '''
546 D = None
547 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
548 E = earth
549 n = _under(name or E.name)
550 elif isinstance(earth, Datum):
551 E = earth.ellipsoid
552 n = _under(name or earth.name)
553 D = earth
554 elif _isRadius(earth):
555 E, D = _ED2(Radius_(earth), name)
556 n = E.name
557 elif isinstance(earth, a_f2Tuple):
558 n = _under(name or earth.name)
559 E = earth.ellipsoid(name=n)
560 elif islistuple(earth, minum=2):
561 E = Ellipsoids.Sphere
562 a, f = earth[:2]
563 if f or a != E.a: # != E.b
564 n = _under(name or _xattr(earth, name=NN))
565 E = Ellipsoid(a, f=f, name=n)
566 else:
567 n = E.name
568 D = Datums.Sphere
569 else:
570 E, n = None, NN
571 return E, n, D
574def _mean_radius(radius, *lats):
575 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
576 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
577 '''
578 if radius is R_M:
579 r = radius
580 elif _isRadius(radius):
581 r = Radius_(radius, low=0, Error=TypeError)
582 else:
583 E = _ellipsoidal_datum(radius).ellipsoid
584 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
585 return r
588def _minus(name): # in .trf
589 '''(INTERNAL) Name of C{inverse} Xform.
590 '''
591 m = _MINUS_
592 return name[len(m):] if name.startswith(m) else NN(m, name)
595def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
596 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
597 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
599 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
600 '''
601 if _isRadius(earth):
602 _, D = _ED2(Radius_(earth, Error=Error), name)
603 else:
604 D = _ellipsoidal_datum(earth, Error=Error, name=name)
605 if raiser and not D.isSpherical:
606 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
607 return D
610class Datums(_NamedEnum):
611 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
612 to accommodate the L{_LazyNamedEnumItem} properties.
613 '''
614 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
615 '''(INTERNAL) Instantiate the L{Datum}.
616 '''
617 return Datum(Ellipsoids.get(ellipsoid_name),
618 Transforms.get(transform_name), name=name)
620Datums = Datums(Datum) # PYCHOK singleton
621'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
622# Datums with associated ellipsoid and Helmert transform parameters
623# to convert from WGS84 into the given datum. More are available at
624# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
625# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
626Datums._assert(
627 # Belgian Datum 1972, based on Hayford ellipsoid.
628 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
629 # <https://SpatialReference.org/ref/sr-org/7718/html/>
630 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
632 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
633 # <https://WikiPedia.org/wiki/Helmert_transformation>
634 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
636 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
637 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
639 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
640# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
641# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
642# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
643 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
645 # <https://WikiPedia.org/wiki/GRS_80>
646 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
648 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
649# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
650 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
652 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
653 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
654 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
656 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
657 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
659 # <https://WikiPedia.org/wiki/Helmert_transformation>
660 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
662 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
663 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
664 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
666 # Nouvelle Triangulation Francaise (Paris) XXX verify
667 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
669 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
670 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
672 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
673 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
675 # XXX psuedo-ellipsoids for spherical LatLon
676 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
678 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
679 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
681 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
682 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
684 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
685)
687_WGS84 = Datums.WGS84
688assert _WGS84.ellipsoid is _EWGS84
689# assert _WGS84.transform.isunity
691if __name__ == '__main__':
693 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
694 from pygeodesy.lazily import printf
696 # __doc__ of this file, force all into registery
697 for r in (Datums, Transforms):
698 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
699 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
701# **) MIT License
702#
703# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
704#
705# Permission is hereby granted, free of charge, to any person obtaining a
706# copy of this software and associated documentation files (the "Software"),
707# to deal in the Software without restriction, including without limitation
708# the rights to use, copy, modify, merge, publish, distribute, sublicense,
709# and/or sell copies of the Software, and to permit persons to whom the
710# Software is furnished to do so, subject to the following conditions:
711#
712# The above copyright notice and this permission notice shall be included
713# in all copies or substantial portions of the Software.
714#
715# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
716# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
717# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
718# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
719# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
720# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
721# OTHER DEALINGS IN THE SOFTWARE.