Coverage for pygeodesy/datums.py: 94%
253 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -0400
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, _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.internals import _passarg, _under
80from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _BAR_, _Bessel1841_, \
81 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \
82 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _MINUS_, \
83 _Krassovski1940_, _Krassowsky1940_, _NAD27_, _NAD83_, _s_, \
84 _PLUS_, _Sphere_, _spherical_, _transform_, _UNDER_, \
85 _WGS72_, _WGS84_
86from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
87from pygeodesy.named import _NamedEnum, _NamedEnumItem, _lazyNamedEnumItem as _lazy
88# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
89from pygeodesy.props import Property_RO, property_RO
90# from pygeodesy.streprs import Fmt # from .fmath
91from pygeodesy.units import _isRadius, Radius_, radians
93# from math import radians # from .units
94# import operator as _operator # from .fmath
96__all__ = _ALL_LAZY.datums
97__version__ = '24.05.14'
99_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
100_BD72_ = 'BD72'
101_DHDN_ = 'DHDN'
102_DHDNE_ = 'DHDNE'
103_DHDNW_ = 'DHDNW'
104_ED50_ = 'ED50'
105_GDA2020_ = 'GDA2020' # in .trf
106_Identity_ = 'Identity'
107_Irl1965_ = 'Irl1965'
108_Irl1975_ = 'Irl1975'
109_MGI_ = 'MGI'
110_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf
111_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:]
112_NTF_ = 'NTF'
113_OSGB36_ = 'OSGB36'
114_Potsdam_ = 'Potsdam'
115_RPS = radians(_1_0 / _3600_0) # radians per arc-second
116_S1_S = 1.e-6 # in .trf
117_TokyoJapan_ = 'TokyoJapan'
120class Transform(_NamedEnumItem):
121 '''Helmert I{datum} transformation.
123 @see: L{TransformXform<trf.TransformXform>}.
124 '''
125 tx = _0_0 # x translation (C{meter})
126 ty = _0_0 # y translation (C{meter})
127 tz = _0_0 # z translation (C{meter})
129 rx = _0_0 # x rotation (C{radians})
130 ry = _0_0 # y rotation (C{radians})
131 rz = _0_0 # z rotation (C{radians})
133 s = _0_0 # scale ppm (C{float})
134 s1 = _1_0 # scale + 1 (C{float})
136 sx = _0_0 # x rotation (C{arc-seconds})
137 sy = _0_0 # y rotation (C{arc-seconds})
138 sz = _0_0 # z rotation (C{arc-seconds})
140 def __init__(self, name=NN, tx=0, ty=0, tz=0,
141 s=0, sx=0, sy=0, sz=0):
142 '''New L{Transform}.
144 @kwarg name: Optional, unique name (C{str}).
145 @kwarg tx: X translation (C{meter}).
146 @kwarg ty: Y translation (C{meter}).
147 @kwarg tz: Z translation (C{meter}).
148 @kwarg s: Scale (C{float}), ppm.
149 @kwarg sx: X rotation (C{arc-seconds}).
150 @kwarg sy: Y rotation (C{arc-seconds}).
151 @kwarg sz: Z rotation (C{arc-seconds}).
153 @raise NameError: Transform with that B{C{name}} already exists.
154 '''
155 if tx:
156 self.tx = tx
157 if ty:
158 self.ty = ty
159 if tz:
160 self.tz = tz
161 if s:
162 self.s = s
163 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1)
164 if sx: # secs to rads
165 self.rx, self.sx = self._rps2(sx)
166 if sy:
167 self.ry, self.sy = self._rps2(sy)
168 if sz:
169 self.rz, self.sz = self._rps2(sz)
171 self._register(Transforms, name)
173 def __eq__(self, other):
174 '''Compare this and an other transform.
176 @arg other: The other transform (L{Transform}).
178 @return: C{True} if equal, C{False} otherwise.
179 '''
180 return self is other or (isinstance(other, Transform)
181 and _equall(other, self))
183 def __hash__(self):
184 return hash(tuple(self))
186 def __iter__(self):
187 '''Yield the initial attribute values, I{in order}.
188 '''
189 for n in _Names7:
190 yield getattr(self, n)
192 def __matmul__(self, point): # PYCHOK Python 3.5+
193 '''Transform an I{ellipsoidal} B{C{point}} with this Helmert.
195 @return: A transformed copy of B{C{point}}.
197 @raise TypeError: Invalid B{C{point}}.
199 @see: Method C{B{point}.toTransform}.
200 '''
201 _ = _xellipsoidall(point)
202 return point.toTransform(self)
204 def __neg__(self):
205 return self.inverse()
207 def inverse(self, name=NN):
208 '''Return the inverse of this transform.
210 @kwarg name: Optional, unique name (C{str}).
212 @return: Inverse (L{Transform}), unregistered.
213 '''
214 r = type(self)(**dict(self.items(inverse=True)))
215 n = name or _negastr(self.name)
216 if n:
217 r.name = n # unregistered
218 return r
220 @Property_RO
221 def isunity(self):
222 '''Is this a C{unity, identity} transform (C{bool}), like
223 WGS84 with translation, scale and rotation all zero?
224 '''
225 return not any(self)
227 def items(self, inverse=False):
228 '''Yield the initial attributes, each as 2-tuple C{(name, value)}.
230 @kwarg inverse: If C{True}, negate the values (C{bool}).
231 '''
232 _p = neg if inverse else _passarg
233 for n, x in _zip(_Names7, self):
234 yield n, _p(x)
236 def _rps2(self, s_):
237 '''(INTERNAL) Rotation in C{radians} and C{arc-seconds}.
238 '''
239 # _MR == _RPS * 1.e-3 # radians per milli-arc-second, equ (2)
240 # <https://www.NGS.NOAA.gov/CORS/Articles/SolerSnayASCE.pdf>
241 return (_RPS * s_), s_
243 def _s_s1(self, s1): # in .trf
244 '''(INTERNAL) Set C{s1} and C{s}.
245 '''
246 Transform.isunity._update(self)
247 self.s1 = s1
248 self.s = s = (s1 - _1_0) / _S1_S
249 return s
251 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected
252 '''Return this transform as a string.
254 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
255 @kwarg fmt: Optional C{float} format (C{letter}).
256 @kwarg name: Override name (C{str}) or C{None} to exclude
257 this transform's name.
259 @return: Transform attributes (C{str}).
260 '''
261 return self._instr(name, prec, *_Names11, fmt=fmt)
263 def transform(self, x, y, z, inverse=False, **Vector_and_kwds):
264 '''Transform a (cartesian) position, forward or inverse.
266 @arg x: X coordinate (C{meter}).
267 @arg y: Y coordinate (C{meter}).
268 @arg z: Z coordinate (C{meter}).
269 @kwarg inverse: If C{True}, apply the inverse transform (C{bool}).
270 @kwarg Vector_and_kwds: An optional, (3-D) C{B{Vector}=None} or
271 cartesian class and additional C{B{Vector}} keyword
272 arguments to return the transformed position.
274 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)})
275 unless some B{C{Vector_and_kwds}} are specified.
276 '''
277 if self.isunity:
278 r = Vector3Tuple(x, y, z, name=self.name) # == inverse
279 else:
280 xyz1 = x, y, z, _1_0
281 s1 = self.s1
282 if inverse:
283 xyz1 = map2(neg, xyz1)
284 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1
285 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
286 # x * .rz + y * .s1 - z * .rx + .ty,
287 # -x * .ry + y * .rx + z * .s1 + .tz)
288 r = Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
289 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
290 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
291 name=self.name)
292 if Vector_and_kwds:
293 V, kwds = _xkwds_pop2(Vector_and_kwds, Vector=None)
294 if V:
295 r = V(r, **_xkwds(kwds, name=self.name))
296 return r
299class Transforms(_NamedEnum):
300 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
301 to accommodate the L{_LazyNamedEnumItem} properties.
302 '''
303 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):
304 '''(INTERNAL) Instantiate the C{Transform}.
305 '''
306 return Transform(**name_tx_ty_tz_s_sx_sy_sz)
308Transforms = Transforms(Transform) # PYCHOK singleton
309'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
310# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
311Transforms._assert(
312 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893), s=_F(1.2727),
313 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
314 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
315 sx=_F( -0.33657), sy=_F( -0.456955), sz=_F( -1.84218)),
317 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0), s=_F(-8.3),
318 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08)),
320 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
322 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39), s=_F(-9.82),
323 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458)), # Germany
325 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2), s=_F(-2.55),
326 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
327 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797)), # East Germany
329 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2), s=_F(-6.7),
330 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
331 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455)), # West Germany
333 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1), s=_F(-1.2),
334 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
335 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
336 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
337 sz=_F( 0.156)),
338 Identity = _lazy(_Identity_),
340 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
341 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
342 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
343 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
344 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
345 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
347 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
348 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
350 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
351 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
353 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920), s=_F(-2.423),
354 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297)), # Austria
356 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
358 NAD83 = _lazy(_NAD83_, tx=_F(1.004), ty=_F(-1.910), tz=_F(-0.515), s=_F(-0.0015),
359 sx=_F(0.0267), sy=_F( 0.00034), sz=_F( 0.011)),
361 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
363 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060), s=_F(20.4894),
364 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
365 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421)),
367 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
369 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), s=_F(-0.22), sz=_F(0.554)),
371 WGS84 = _lazy(_WGS84_), # unity
372)
375class Datum(_NamedEnumItem):
376 '''Ellipsoid and transform parameters for an earth model.
377 '''
378 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
379 _transform = Transforms.WGS84 # default transform (L{Transform})
381 def __init__(self, ellipsoid, transform=None, name=NN):
382 '''New L{Datum}.
384 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
385 @kwarg transform: Optional transform (L{Transform}).
386 @kwarg name: Optional, unique name (C{str}).
388 @raise NameError: Datum with that B{C{name}} already exists.
390 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
391 nor L{Ellipsoid2} or B{C{transform}} is
392 not a L{Transform}.
393 '''
394 self._ellipsoid = ellipsoid or Datum._ellipsoid
395 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
397 self._transform = transform or Datum._transform
398 _xinstanceof(Transform, transform=self.transform)
400 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
402 def __eq__(self, other):
403 '''Compare this and an other datum.
405 @arg other: The other datum (L{Datum}).
407 @return: C{True} if equal, C{False} otherwise.
408 '''
409 return self is other or (isinstance(other, Datum) and
410 self.ellipsoid == other.ellipsoid and
411 self.transform == other.transform)
413 def __hash__(self):
414 return self._hash # memoized
416 def __matmul__(self, point): # PYCHOK Python 3.5+
417 '''Convert an I{ellipsoidal} B{C{point}} to this datum.
419 @raise TypeError: Invalid B{C{point}}.
420 '''
421 _ = _xellipsoidall(point)
422 return point.toDatum(self)
424 def ecef(self, Ecef=None):
425 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
427 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
429 @return: An ECEF converter for this C{datum}.
431 @raise TypeError: Invalid B{C{Ecef}}.
433 @see: Module L{pygeodesy.ecef}.
434 '''
435 return _MODS.ecef._4Ecef(self, Ecef)
437 @Property_RO
438 def ellipsoid(self):
439 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
440 '''
441 return self._ellipsoid
443 @Property_RO
444 def exactTM(self):
445 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
446 '''
447 return _MODS.etm.ExactTransverseMercator(datum=self)
449 @Property_RO
450 def _hash(self):
451 return hash(self.ellipsoid) + hash(self.transform)
453 @property_RO
454 def isEllipsoidal(self):
455 '''Check whether this datum is ellipsoidal (C{bool}).
456 '''
457 return self.ellipsoid.isEllipsoidal
459 @property_RO
460 def isOblate(self):
461 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
462 '''
463 return self.ellipsoid.isOblate
465 @property_RO
466 def isProlate(self):
467 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
468 '''
469 return self.ellipsoid.isProlate
471 @property_RO
472 def isSpherical(self):
473 '''Check whether this datum is (near-)spherical (C{bool}).
474 '''
475 return self.ellipsoid.isSpherical
477 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
478 '''Return this datum as a string.
480 @kwarg sep: Separator to join (C{str}).
481 @kwarg name: Override name (C{str}) or C{None} to exclude
482 this datum's name.
484 @return: Datum attributes (C{str}).
485 '''
486 t = [] if name is None else \
487 [Fmt.EQUAL(name=repr(name or self.named))]
488 for a in (_ellipsoid_, _transform_):
489 v = getattr(self, a)
490 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
491 return sep.join(t)
493 @Property_RO
494 def transform(self):
495 '''Get this datum's transform (L{Transform}).
496 '''
497 return self._transform
500def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
501 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
502 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
504 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
505 '''
506 if f is not None:
507 E, n, D = _EnD3((a_earth, f), name)
508 if raiser and not E:
509 raise _TypeError(f=f, **{raiser: a_earth})
510 elif a_earth in (_EWGS84, _WGS84, None) and inst._datum is _WGS84:
511 return
512 elif isinstance(a_earth, Datum):
513 E, n, D = None, NN, a_earth
514 else:
515 E, n, D = _EnD3(a_earth, name)
516 if raiser and not E:
517 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
518 if D is None:
519 D = Datum(E, transform=Transforms.Identity, name=_under(n))
520 inst._datum = D
523def _earth_ellipsoid(earth, *name_raiser):
524 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
525 '''
526 return Ellipsoids.Sphere if earth is R_M else (
527 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
528 _spherical_datum(earth, *name_raiser).ellipsoid)
531def _ED2(radius, name):
532 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
533 '''
534 D = Datums.Sphere
535 E = D.ellipsoid
536 if name or radius != E.a: # != E.b
537 n = _under(name)
538 E = Ellipsoid(radius, radius, name=n)
539 D = Datum(E, transform=Transforms.Identity, name=n)
540 return E, D
543def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
544 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
545 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
547 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
548 '''
549 if isinstance(earth, Datum):
550 D = earth
551 else:
552 E, n, D = _EnD3(earth, name)
553 if not E:
554 n = raiser or _earth_
555 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
556 if D is None:
557 D = Datum(E, transform=Transforms.Identity, name=_under(n))
558 if raiser and not D.isEllipsoidal:
559 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
560 return D
563def _EnD3(earth, name):
564 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
565 '''
566 D = None
567 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
568 E = earth
569 n = _under(name or E.name)
570 elif isinstance(earth, Datum):
571 E = earth.ellipsoid
572 n = _under(name or earth.name)
573 D = earth
574 elif _isRadius(earth):
575 E, D = _ED2(Radius_(earth), name)
576 n = E.name
577 elif isinstance(earth, a_f2Tuple):
578 n = _under(name or earth.name)
579 E = earth.ellipsoid(name=n)
580 elif islistuple(earth, minum=2):
581 E = Ellipsoids.Sphere
582 a, f = earth[:2]
583 if f or a != E.a: # != E.b
584 n = _under(name or _xattr(earth, name=NN))
585 E = Ellipsoid(a, f=f, name=n)
586 else:
587 n = E.name
588 D = Datums.Sphere
589 else:
590 E, n = None, NN
591 return E, n, D
594def _equall(t1, t2): # in .trf
595 '''(INTERNAL) Return L{Transform} C{t1 == t2}.
596 '''
597 return all(map(_operator.eq, t1, t2))
600def _mean_radius(radius, *lats):
601 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
602 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
603 '''
604 if radius is R_M:
605 r = radius
606 elif _isRadius(radius):
607 r = Radius_(radius, low=0, Error=TypeError)
608 else:
609 E = _ellipsoidal_datum(radius).ellipsoid
610 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
611 return r
614def _negastr(name): # in .trf
615 '''(INTERNAL) Negate a C{Transform/-Xform} name.
616 '''
617 b, m, p = _BAR_, _MINUS_, _PLUS_
618 n = name.replace(m, b).replace(p, m).replace(b, p)
619 # as good and fast as (in Python 3+ only) ...
620 # _MINUSxPLUS = str.maketrans({_MINUS_: _PLUS_, _PLUS_: _MINUS_})
621 # def _negastr(name):
622 # n = name.translate(_MINUSxPLUS)
623 # ...
624 return n.lstrip(p) if n.startswith(p) else NN(m, n)
627def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
628 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
629 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
631 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
632 '''
633 if _isRadius(earth):
634 _, D = _ED2(Radius_(earth, Error=Error), name)
635 else:
636 D = _ellipsoidal_datum(earth, Error=Error, name=name)
637 if raiser and not D.isSpherical:
638 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
639 return D
642class Datums(_NamedEnum):
643 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
644 to accommodate the L{_LazyNamedEnumItem} properties.
645 '''
646 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
647 '''(INTERNAL) Instantiate the L{Datum}.
648 '''
649 return Datum(Ellipsoids.get(ellipsoid_name),
650 Transforms.get(transform_name), name=name)
652Datums = Datums(Datum) # PYCHOK singleton
653'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
654# Datums with associated ellipsoid and Helmert transform parameters
655# to convert from WGS84 into the given datum. More are available at
656# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
657# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
658Datums._assert(
659 # Belgian Datum 1972, based on Hayford ellipsoid.
660 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
661 # <https://SpatialReference.org/ref/sr-org/7718/html/>
662 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
664 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
665 # <https://WikiPedia.org/wiki/Helmert_transformation>
666 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
668 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
669 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
671 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
672# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
673# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
674# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
675 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
677 # <https://WikiPedia.org/wiki/GRS_80>
678 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
680 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
681# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
682 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
684 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
685 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
686 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
688 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
689 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
691 # <https://WikiPedia.org/wiki/Helmert_transformation>
692 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
694 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
695 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
696 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
698 # Nouvelle Triangulation Francaise (Paris) XXX verify
699 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
701 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
702 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
704 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
705 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
707 # XXX psuedo-ellipsoids for spherical LatLon
708 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
710 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
711 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
713 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
714 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
716 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
717)
719_WGS84 = Datums.WGS84
720assert _WGS84.ellipsoid is _EWGS84
721# assert _WGS84.transform.isunity
723if __name__ == '__main__':
725 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
726 from pygeodesy import printf
728 # __doc__ of this file, force all into registery
729 for r in (Datums, Transforms):
730 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
731 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
733# **) MIT License
734#
735# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
736#
737# Permission is hereby granted, free of charge, to any person obtaining a
738# copy of this software and associated documentation files (the "Software"),
739# to deal in the Software without restriction, including without limitation
740# the rights to use, copy, modify, merge, publish, distribute, sublicense,
741# and/or sell copies of the Software, and to permit persons to whom the
742# Software is furnished to do so, subject to the following conditions:
743#
744# The above copyright notice and this permission notice shall be included
745# in all copies or substantial portions of the Software.
746#
747# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
748# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
749# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
750# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
751# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
752# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
753# OTHER DEALINGS IN THE SOFTWARE.