Coverage for pygeodesy/datums.py: 94%
248 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-03-03 11:31 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2024-03-03 11:31 -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.87, ty=-52.298, tz=103.72, s1=1.0, rx=-1.6317e-06, ry=-2.2154e-06, rz=-8.9311e-06, s=1.2727, sx=-0.33657, sy=-0.45696, sz=-1.8422)
47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, s1=0.99999, rx=-5.0421e-06, ry=-1.6968e-06, rz=1.4932e-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.5682e-07, rz=-7.0686e-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=1.0, rx=2.618e-07, ry=-2.7634e-07, rz=1.356e-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.99999, rx=-9.7932e-07, ry=-2.1817e-07, rz=1.1902e-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=1.0, rx=0.0, ry=0.0, rz=7.5631e-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.6, tz=-564.56, s1=0.99999, rx=5.0518e-06, ry=1.0375e-06, rz=3.0592e-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.6, tz=-564.56, s1=0.99999, rx=5.0518e-06, ry=1.0375e-06, rz=3.0592e-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=1.0, rx=-9.6963e-08, ry=1.2605e-06, rz=6.3026e-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=1.0, rx=-9.6963e-08, ry=1.2605e-06, rz=6.3026e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)
58@var Transforms.MGI: Transform(name='MGI', tx=-577.33, ty=-90.129, tz=-463.92, s1=1.0, rx=2.4905e-05, ry=7.1462e-06, rz=2.5681e-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.2945e-07, ry=1.6484e-09, rz=5.333e-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.45, ty=125.16, tz=-542.06, s1=1.0, rx=-7.2819e-07, ry=-1.1975e-06, rz=-4.0826e-06, s=20.489, 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.6859e-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, _passarg, _xinstanceof, _zip
71from pygeodesy.constants import R_M, _float as _F, _0_0, _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, _xellipsoidall, _xkwds, \
77 _xkwds_pop2
78from pygeodesy.fmath import fdot, fmean, Fmt, _operator
79from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, \
80 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \
81 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \
82 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \
83 _NAD83_, _PLUS_, _s_, _Sphere_, _spherical_, _transform_, \
84 _UNDER_, _WGS72_, _WGS84_, _under
85from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
86from pygeodesy.named import _NamedEnum, _NamedEnumItem, _lazyNamedEnumItem as _lazy
87# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
88from pygeodesy.props import Property_RO, property_RO
89# from pygeodesy.streprs import Fmt # from .fmath
90from pygeodesy.units import _isRadius, Radius_, radians
92# from math import radians # from .units
93# import operator as _operator # from .fmath
95__all__ = _ALL_LAZY.datums
96__version__ = '24.02.29'
98_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
99_BD72_ = 'BD72'
100_DHDN_ = 'DHDN'
101_DHDNE_ = 'DHDNE'
102_DHDNW_ = 'DHDNW'
103_ED50_ = 'ED50'
104_GDA2020_ = 'GDA2020' # in .trf
105_Identity_ = 'Identity'
106_Irl1965_ = 'Irl1965'
107_Irl1975_ = 'Irl1975'
108_MGI_ = 'MGI'
109_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf
110_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:]
111_NTF_ = 'NTF'
112_OSGB36_ = 'OSGB36'
113_Potsdam_ = 'Potsdam'
114_RPS = radians(_1_0 / _3600_0) # radians per degree-second
115_S1_S = 1.e-6 # in .trf
116_TokyoJapan_ = 'TokyoJapan'
119class Transform(_NamedEnumItem):
120 '''Helmert I{datum} transformation.
122 @see: L{TransformXform<trf.TransformXform>}.
123 '''
124 tx = _0_0 # x translation (C{meter})
125 ty = _0_0 # y translation (C{meter})
126 tz = _0_0 # z translation (C{meter})
128 rx = _0_0 # x rotation (C{radians})
129 ry = _0_0 # y rotation (C{radians})
130 rz = _0_0 # z rotation (C{radians})
132 s = _0_0 # scale ppm (C{float})
133 s1 = _1_0 # scale + 1 (C{float})
135 sx = _0_0 # x rotation (C{degree-seconds})
136 sy = _0_0 # y rotation (C{degree-seconds})
137 sz = _0_0 # z rotation (C{degree-seconds})
139 def __init__(self, name=NN, tx=0, ty=0, tz=0,
140 s=0, sx=0, sy=0, sz=0):
141 '''New L{Transform}.
143 @kwarg name: Optional, unique name (C{str}).
144 @kwarg tx: X translation (C{meter}).
145 @kwarg ty: Y translation (C{meter}).
146 @kwarg tz: Z translation (C{meter}).
147 @kwarg s: Scale (C{float}), ppm.
148 @kwarg sx: X rotation (C{degree-seconds}).
149 @kwarg sy: Y rotation (C{degree-seconds}).
150 @kwarg sz: Z rotation (C{degree-seconds}).
152 @raise NameError: Transform with that B{C{name}} already exists.
153 '''
154 if tx:
155 self.tx = tx
156 if ty:
157 self.ty = ty
158 if tz:
159 self.tz = tz
160 if s:
161 self.s = s
162 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1)
163 if sx: # secs to rads
164 self.rx, self.sx = self._rps2(sx)
165 if sy:
166 self.ry, self.sy = self._rps2(sy)
167 if sz:
168 self.rz, self.sz = self._rps2(sz)
170 self._register(Transforms, name)
172 def __eq__(self, other):
173 '''Compare this and an other transform.
175 @arg other: The other transform (L{Transform}).
177 @return: C{True} if equal, C{False} otherwise.
178 '''
179 return self is other or (isinstance(other, Transform)
180 and _equall(other, self))
182 def __hash__(self):
183 return hash(tuple(self))
185 def __iter__(self):
186 '''Yield the initial attribute values, I{in order}.
187 '''
188 for n in _Names7:
189 yield getattr(self, n)
191 def __matmul__(self, point): # PYCHOK Python 3.5+
192 '''Transform an I{ellipsoidal} B{C{point}} with this Helmert.
194 @return: A transformed copy of B{C{point}}.
196 @raise TypeError: Invalid B{C{point}}.
198 @see: Method C{B{point}.toTransform}.
199 '''
200 _ = _xellipsoidall(point)
201 return point.toTransform(self)
203 def __neg__(self):
204 return self.inverse()
206 def inverse(self, name=NN):
207 '''Return the inverse of this transform.
209 @kwarg name: Optional, unique name (C{str}).
211 @return: Inverse (L{Transform}), unregistered.
212 '''
213 r = type(self)(**dict(self.items(inverse=True)))
214 n = name or _negstr(self.name)
215 if n:
216 r.name = n # unregistered
217 return r
219 @Property_RO
220 def isunity(self):
221 '''Is this a C{unity, identity} transform (C{bool}), like
222 WGS84 with translation, scale and rotation all zero?
223 '''
224 return not any(self)
226 def items(self, inverse=False):
227 '''Yield the initial attributes, each as 2-tuple C{(name, value)}.
229 @kwarg inverse: If C{True}, negate the values (C{bool}).
230 '''
231 _p = neg if inverse else _passarg
232 for n, x in _zip(_Names7, self):
233 yield n, _p(x)
235 def _rps2(self, s_):
236 '''(INTERNAL) Rotation in C{radians} and C{degree-seconds}.
237 '''
238 return (_RPS * s_), s_
240 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected
241 '''Return this transform as a string.
243 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
244 @kwarg fmt: Optional C{float} format (C{letter}).
245 @kwarg name: Override name (C{str}) or C{None} to exclude
246 this transform's name.
248 @return: Transform attributes (C{str}).
249 '''
250 return self._instr(name, prec, *_Names11, fmt=fmt)
252 def transform(self, x, y, z, inverse=False, **Vector_and_kwds):
253 '''Transform a (cartesian) position, forward or inverse.
255 @arg x: X coordinate (C{meter}).
256 @arg y: Y coordinate (C{meter}).
257 @arg z: Z coordinate (C{meter}).
258 @kwarg inverse: If C{True}, apply the inverse transform (C{bool}).
259 @kwarg Vector_and_kwds: An optional, (3-D) C{B{Vector}=None} or
260 cartesian class and additional C{B{Vector}} keyword
261 arguments to return the transformed position.
263 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)})
264 unless some B{C{Vector_and_kwds}} are specified.
265 '''
266 if self.isunity:
267 r = Vector3Tuple(x, y, z, name=self.name) # == inverse
268 else:
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 r = 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)
281 if Vector_and_kwds:
282 V, kwds = _xkwds_pop2(Vector_and_kwds, Vector=None)
283 if V:
284 r = V(r, **_xkwds(kwds, name=self.name))
285 return r
288class Transforms(_NamedEnum):
289 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
290 to accommodate the L{_LazyNamedEnumItem} properties.
291 '''
292 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):
293 '''(INTERNAL) Instantiate the C{Transform}.
294 '''
295 return Transform(**name_tx_ty_tz_s_sx_sy_sz)
297Transforms = Transforms(Transform) # PYCHOK singleton
298'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
299# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
300Transforms._assert(
301 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893), s=_F(1.2727),
302 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
303 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
304 sx=_F( -0.33657), sy=_F( -0.456955), sz=_F( -1.84218)),
306 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0), s=_F(-8.3),
307 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08)),
309 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
311 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39), s=_F(-9.82),
312 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458)), # Germany
314 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2), s=_F(-2.55),
315 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
316 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797)), # East Germany
318 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2), s=_F(-6.7),
319 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
320 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455)), # West Germany
322 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1), s=_F(-1.2),
323 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
324 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
325 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
326 sz=_F( 0.156)),
327 Identity = _lazy(_Identity_),
329 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
330 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
331 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
332 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
333 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
334 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
336 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
337 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
339 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
340 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
342 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920), s=_F(-2.423),
343 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297)), # Austria
345 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
347 NAD83 = _lazy(_NAD83_, tx=_F(1.004), ty=_F(-1.910), tz=_F(-0.515), s=_F(-0.0015),
348 sx=_F(0.0267), sy=_F( 0.00034), sz=_F( 0.011)),
350 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
352 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060), s=_F(20.4894),
353 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
354 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421)),
356 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
358 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), s=_F(-0.22), sz=_F(0.554)),
360 WGS84 = _lazy(_WGS84_), # unity
361)
364class Datum(_NamedEnumItem):
365 '''Ellipsoid and transform parameters for an earth model.
366 '''
367 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
368 _transform = Transforms.WGS84 # default transform (L{Transform})
370 def __init__(self, ellipsoid, transform=None, name=NN):
371 '''New L{Datum}.
373 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
374 @kwarg transform: Optional transform (L{Transform}).
375 @kwarg name: Optional, unique name (C{str}).
377 @raise NameError: Datum with that B{C{name}} already exists.
379 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
380 nor L{Ellipsoid2} or B{C{transform}} is
381 not a L{Transform}.
382 '''
383 self._ellipsoid = ellipsoid or Datum._ellipsoid
384 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
386 self._transform = transform or Datum._transform
387 _xinstanceof(Transform, transform=self.transform)
389 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
391 def __eq__(self, other):
392 '''Compare this and an other datum.
394 @arg other: The other datum (L{Datum}).
396 @return: C{True} if equal, C{False} otherwise.
397 '''
398 return self is other or (isinstance(other, Datum) and
399 self.ellipsoid == other.ellipsoid and
400 self.transform == other.transform)
402 def __hash__(self):
403 return self._hash # memoized
405 def __matmul__(self, point): # PYCHOK Python 3.5+
406 '''Convert an I{ellipsoidal} B{C{point}} to this datum.
408 @raise TypeError: Invalid B{C{point}}.
409 '''
410 _ = _xellipsoidall(point)
411 return point.toDatum(self)
413 def ecef(self, Ecef=None):
414 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
416 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
418 @return: An ECEF converter for this C{datum}.
420 @raise TypeError: Invalid B{C{Ecef}}.
422 @see: Module L{pygeodesy.ecef}.
423 '''
424 return _MODS.ecef._4Ecef(self, Ecef)
426 @Property_RO
427 def ellipsoid(self):
428 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
429 '''
430 return self._ellipsoid
432 @Property_RO
433 def exactTM(self):
434 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
435 '''
436 return _MODS.etm.ExactTransverseMercator(datum=self)
438 @Property_RO
439 def _hash(self):
440 return hash(self.ellipsoid) + hash(self.transform)
442 @property_RO
443 def isEllipsoidal(self):
444 '''Check whether this datum is ellipsoidal (C{bool}).
445 '''
446 return self.ellipsoid.isEllipsoidal
448 @property_RO
449 def isOblate(self):
450 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
451 '''
452 return self.ellipsoid.isOblate
454 @property_RO
455 def isProlate(self):
456 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
457 '''
458 return self.ellipsoid.isProlate
460 @property_RO
461 def isSpherical(self):
462 '''Check whether this datum is (near-)spherical (C{bool}).
463 '''
464 return self.ellipsoid.isSpherical
466 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
467 '''Return this datum as a string.
469 @kwarg sep: Separator to join (C{str}).
470 @kwarg name: Override name (C{str}) or C{None} to exclude
471 this datum's name.
473 @return: Datum attributes (C{str}).
474 '''
475 t = [] if name is None else \
476 [Fmt.EQUAL(name=repr(name or self.named))]
477 for a in (_ellipsoid_, _transform_):
478 v = getattr(self, a)
479 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
480 return sep.join(t)
482 @Property_RO
483 def transform(self):
484 '''Get this datum's transform (L{Transform}).
485 '''
486 return self._transform
489def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
490 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
491 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
493 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
494 '''
495 if f is not None:
496 E, n, D = _EnD3((a_earth, f), name)
497 if raiser and not E:
498 raise _TypeError(f=f, **{raiser: a_earth})
499 elif a_earth in (_EWGS84, _WGS84, None) and inst._datum is _WGS84:
500 return
501 elif isinstance(a_earth, Datum):
502 E, n, D = None, NN, a_earth
503 else:
504 E, n, D = _EnD3(a_earth, name)
505 if raiser and not E:
506 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
507 if D is None:
508 D = Datum(E, transform=Transforms.Identity, name=_under(n))
509 inst._datum = D
512def _earth_ellipsoid(earth, *name_raiser):
513 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
514 '''
515 return Ellipsoids.Sphere if earth is R_M else (
516 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
517 _spherical_datum(earth, *name_raiser).ellipsoid)
520def _ED2(radius, name):
521 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
522 '''
523 D = Datums.Sphere
524 E = D.ellipsoid
525 if name or radius != E.a: # != E.b
526 n = _under(name)
527 E = Ellipsoid(radius, radius, name=n)
528 D = Datum(E, transform=Transforms.Identity, name=n)
529 return E, D
532def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
533 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
534 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
536 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
537 '''
538 if isinstance(earth, Datum):
539 D = earth
540 else:
541 E, n, D = _EnD3(earth, name)
542 if not E:
543 n = raiser or _earth_
544 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
545 if D is None:
546 D = Datum(E, transform=Transforms.Identity, name=_under(n))
547 if raiser and not D.isEllipsoidal:
548 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
549 return D
552def _EnD3(earth, name):
553 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
554 '''
555 D = None
556 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
557 E = earth
558 n = _under(name or E.name)
559 elif isinstance(earth, Datum):
560 E = earth.ellipsoid
561 n = _under(name or earth.name)
562 D = earth
563 elif _isRadius(earth):
564 E, D = _ED2(Radius_(earth), name)
565 n = E.name
566 elif isinstance(earth, a_f2Tuple):
567 n = _under(name or earth.name)
568 E = earth.ellipsoid(name=n)
569 elif islistuple(earth, minum=2):
570 E = Ellipsoids.Sphere
571 a, f = earth[:2]
572 if f or a != E.a: # != E.b
573 n = _under(name or _xattr(earth, name=NN))
574 E = Ellipsoid(a, f=f, name=n)
575 else:
576 n = E.name
577 D = Datums.Sphere
578 else:
579 E, n = None, NN
580 return E, n, D
583def _equall(t1, t2): # in .trf
584 '''(INTERNAL) Return L{Transform} C{t1 == t2}.
585 '''
586 return all(map(_operator.eq, t1, t2))
589def _mean_radius(radius, *lats):
590 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
591 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
592 '''
593 if radius is R_M:
594 r = radius
595 elif _isRadius(radius):
596 r = Radius_(radius, low=0, Error=TypeError)
597 else:
598 E = _ellipsoidal_datum(radius).ellipsoid
599 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
600 return r
603try:
604 _MINUSxPLUS = str.maketrans({_MINUS_: _PLUS_, _PLUS_: _MINUS_})
606 def _negstr(name): # in .trf
607 '''(INTERNAL) Negate a C{Transform/-Xform} name.
608 '''
609 n = name.translate(_MINUSxPLUS)
610 return n.lstrip(_PLUS_) if n.startswith(_PLUS_) else NN(_MINUS_, n)
612except AttributeError: # no str.maketran in Python 2-
613 from pygeodesy.interns import _BAR_
615 def _negstr(name): # PYCHOK in .trf
616 '''(INTERNAL) Negate a C{Transform/-Xform} name.
617 '''
618 b, m, p = _BAR_, _MINUS_, _PLUS_
619 n = name.replace(m, b).replace(p, m).replace(b, p)
620 return n.lstrip(p) if n.startswith(p) else NN(m, n)
623def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
624 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
625 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
627 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
628 '''
629 if _isRadius(earth):
630 _, D = _ED2(Radius_(earth, Error=Error), name)
631 else:
632 D = _ellipsoidal_datum(earth, Error=Error, name=name)
633 if raiser and not D.isSpherical:
634 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
635 return D
638class Datums(_NamedEnum):
639 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
640 to accommodate the L{_LazyNamedEnumItem} properties.
641 '''
642 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
643 '''(INTERNAL) Instantiate the L{Datum}.
644 '''
645 return Datum(Ellipsoids.get(ellipsoid_name),
646 Transforms.get(transform_name), name=name)
648Datums = Datums(Datum) # PYCHOK singleton
649'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
650# Datums with associated ellipsoid and Helmert transform parameters
651# to convert from WGS84 into the given datum. More are available at
652# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
653# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
654Datums._assert(
655 # Belgian Datum 1972, based on Hayford ellipsoid.
656 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
657 # <https://SpatialReference.org/ref/sr-org/7718/html/>
658 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
660 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
661 # <https://WikiPedia.org/wiki/Helmert_transformation>
662 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
664 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
665 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
667 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
668# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
669# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
670# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
671 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
673 # <https://WikiPedia.org/wiki/GRS_80>
674 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
676 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
677# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
678 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
680 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
681 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
682 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
684 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
685 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
687 # <https://WikiPedia.org/wiki/Helmert_transformation>
688 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
690 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
691 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
692 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
694 # Nouvelle Triangulation Francaise (Paris) XXX verify
695 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
697 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
698 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
700 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
701 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
703 # XXX psuedo-ellipsoids for spherical LatLon
704 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
706 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
707 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
709 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
710 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
712 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
713)
715_WGS84 = Datums.WGS84
716assert _WGS84.ellipsoid is _EWGS84
717# assert _WGS84.transform.isunity
719if __name__ == '__main__':
721 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
722 from pygeodesy.lazily import printf
724 # __doc__ of this file, force all into registery
725 for r in (Datums, Transforms):
726 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
727 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
729# **) MIT License
730#
731# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
732#
733# Permission is hereby granted, free of charge, to any person obtaining a
734# copy of this software and associated documentation files (the "Software"),
735# to deal in the Software without restriction, including without limitation
736# the rights to use, copy, modify, merge, publish, distribute, sublicense,
737# and/or sell copies of the Software, and to permit persons to whom the
738# Software is furnished to do so, subject to the following conditions:
739#
740# The above copyright notice and this permission notice shall be included
741# in all copies or substantial portions of the Software.
742#
743# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
744# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
745# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
746# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
747# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
748# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
749# OTHER DEALINGS IN THE SOFTWARE.