Coverage for pygeodesy/datums.py: 94%
252 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-06 16:50 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-06 16:50 -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, _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_, _BAR_, _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.03.07'
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 arc-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{arc-seconds})
136 sy = _0_0 # y rotation (C{arc-seconds})
137 sz = _0_0 # z rotation (C{arc-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{arc-seconds}).
149 @kwarg sy: Y rotation (C{arc-seconds}).
150 @kwarg sz: Z rotation (C{arc-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 _negastr(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{arc-seconds}.
237 '''
238 # _MR == _RPS * 1.e-3 # radians per milli-arc-second, equ (2)
239 # <https://www.NGS.NOAA.gov/CORS/Articles/SolerSnayASCE.pdf>
240 return (_RPS * s_), s_
242 def _s_s1(self, s1): # in .trf
243 '''(INTERNAL) Set C{s1} and C{s}.
244 '''
245 Transform.isunity._update(self)
246 self.s1 = s1
247 self.s = s = (s1 - _1_0) / _S1_S
248 return s
250 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected
251 '''Return this transform as a string.
253 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
254 @kwarg fmt: Optional C{float} format (C{letter}).
255 @kwarg name: Override name (C{str}) or C{None} to exclude
256 this transform's name.
258 @return: Transform attributes (C{str}).
259 '''
260 return self._instr(name, prec, *_Names11, fmt=fmt)
262 def transform(self, x, y, z, inverse=False, **Vector_and_kwds):
263 '''Transform a (cartesian) position, forward or inverse.
265 @arg x: X coordinate (C{meter}).
266 @arg y: Y coordinate (C{meter}).
267 @arg z: Z coordinate (C{meter}).
268 @kwarg inverse: If C{True}, apply the inverse transform (C{bool}).
269 @kwarg Vector_and_kwds: An optional, (3-D) C{B{Vector}=None} or
270 cartesian class and additional C{B{Vector}} keyword
271 arguments to return the transformed position.
273 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)})
274 unless some B{C{Vector_and_kwds}} are specified.
275 '''
276 if self.isunity:
277 r = Vector3Tuple(x, y, z, name=self.name) # == inverse
278 else:
279 xyz1 = x, y, z, _1_0
280 s1 = self.s1
281 if inverse:
282 xyz1 = map2(neg, xyz1)
283 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1
284 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
285 # x * .rz + y * .s1 - z * .rx + .ty,
286 # -x * .ry + y * .rx + z * .s1 + .tz)
287 r = Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
288 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
289 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
290 name=self.name)
291 if Vector_and_kwds:
292 V, kwds = _xkwds_pop2(Vector_and_kwds, Vector=None)
293 if V:
294 r = V(r, **_xkwds(kwds, name=self.name))
295 return r
298class Transforms(_NamedEnum):
299 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
300 to accommodate the L{_LazyNamedEnumItem} properties.
301 '''
302 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):
303 '''(INTERNAL) Instantiate the C{Transform}.
304 '''
305 return Transform(**name_tx_ty_tz_s_sx_sy_sz)
307Transforms = Transforms(Transform) # PYCHOK singleton
308'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
309# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
310Transforms._assert(
311 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893), s=_F(1.2727),
312 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
313 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
314 sx=_F( -0.33657), sy=_F( -0.456955), sz=_F( -1.84218)),
316 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0), s=_F(-8.3),
317 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08)),
319 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
321 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39), s=_F(-9.82),
322 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458)), # Germany
324 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2), s=_F(-2.55),
325 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
326 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797)), # East Germany
328 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2), s=_F(-6.7),
329 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
330 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455)), # West Germany
332 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1), s=_F(-1.2),
333 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
334 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
335 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
336 sz=_F( 0.156)),
337 Identity = _lazy(_Identity_),
339 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
340 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
341 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
342 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),
343 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
344 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631)),
346 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
347 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
349 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),
350 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling
352 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920), s=_F(-2.423),
353 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297)), # Austria
355 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
357 NAD83 = _lazy(_NAD83_, tx=_F(1.004), ty=_F(-1.910), tz=_F(-0.515), s=_F(-0.0015),
358 sx=_F(0.0267), sy=_F( 0.00034), sz=_F( 0.011)),
360 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
362 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060), s=_F(20.4894),
363 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
364 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421)),
366 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
368 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), s=_F(-0.22), sz=_F(0.554)),
370 WGS84 = _lazy(_WGS84_), # unity
371)
374class Datum(_NamedEnumItem):
375 '''Ellipsoid and transform parameters for an earth model.
376 '''
377 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
378 _transform = Transforms.WGS84 # default transform (L{Transform})
380 def __init__(self, ellipsoid, transform=None, name=NN):
381 '''New L{Datum}.
383 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
384 @kwarg transform: Optional transform (L{Transform}).
385 @kwarg name: Optional, unique name (C{str}).
387 @raise NameError: Datum with that B{C{name}} already exists.
389 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
390 nor L{Ellipsoid2} or B{C{transform}} is
391 not a L{Transform}.
392 '''
393 self._ellipsoid = ellipsoid or Datum._ellipsoid
394 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
396 self._transform = transform or Datum._transform
397 _xinstanceof(Transform, transform=self.transform)
399 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
401 def __eq__(self, other):
402 '''Compare this and an other datum.
404 @arg other: The other datum (L{Datum}).
406 @return: C{True} if equal, C{False} otherwise.
407 '''
408 return self is other or (isinstance(other, Datum) and
409 self.ellipsoid == other.ellipsoid and
410 self.transform == other.transform)
412 def __hash__(self):
413 return self._hash # memoized
415 def __matmul__(self, point): # PYCHOK Python 3.5+
416 '''Convert an I{ellipsoidal} B{C{point}} to this datum.
418 @raise TypeError: Invalid B{C{point}}.
419 '''
420 _ = _xellipsoidall(point)
421 return point.toDatum(self)
423 def ecef(self, Ecef=None):
424 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
426 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
428 @return: An ECEF converter for this C{datum}.
430 @raise TypeError: Invalid B{C{Ecef}}.
432 @see: Module L{pygeodesy.ecef}.
433 '''
434 return _MODS.ecef._4Ecef(self, Ecef)
436 @Property_RO
437 def ellipsoid(self):
438 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
439 '''
440 return self._ellipsoid
442 @Property_RO
443 def exactTM(self):
444 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
445 '''
446 return _MODS.etm.ExactTransverseMercator(datum=self)
448 @Property_RO
449 def _hash(self):
450 return hash(self.ellipsoid) + hash(self.transform)
452 @property_RO
453 def isEllipsoidal(self):
454 '''Check whether this datum is ellipsoidal (C{bool}).
455 '''
456 return self.ellipsoid.isEllipsoidal
458 @property_RO
459 def isOblate(self):
460 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
461 '''
462 return self.ellipsoid.isOblate
464 @property_RO
465 def isProlate(self):
466 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
467 '''
468 return self.ellipsoid.isProlate
470 @property_RO
471 def isSpherical(self):
472 '''Check whether this datum is (near-)spherical (C{bool}).
473 '''
474 return self.ellipsoid.isSpherical
476 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
477 '''Return this datum as a string.
479 @kwarg sep: Separator to join (C{str}).
480 @kwarg name: Override name (C{str}) or C{None} to exclude
481 this datum's name.
483 @return: Datum attributes (C{str}).
484 '''
485 t = [] if name is None else \
486 [Fmt.EQUAL(name=repr(name or self.named))]
487 for a in (_ellipsoid_, _transform_):
488 v = getattr(self, a)
489 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
490 return sep.join(t)
492 @Property_RO
493 def transform(self):
494 '''Get this datum's transform (L{Transform}).
495 '''
496 return self._transform
499def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
500 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
501 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
503 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
504 '''
505 if f is not None:
506 E, n, D = _EnD3((a_earth, f), name)
507 if raiser and not E:
508 raise _TypeError(f=f, **{raiser: a_earth})
509 elif a_earth in (_EWGS84, _WGS84, None) and inst._datum is _WGS84:
510 return
511 elif isinstance(a_earth, Datum):
512 E, n, D = None, NN, a_earth
513 else:
514 E, n, D = _EnD3(a_earth, name)
515 if raiser and not E:
516 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
517 if D is None:
518 D = Datum(E, transform=Transforms.Identity, name=_under(n))
519 inst._datum = D
522def _earth_ellipsoid(earth, *name_raiser):
523 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
524 '''
525 return Ellipsoids.Sphere if earth is R_M else (
526 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
527 _spherical_datum(earth, *name_raiser).ellipsoid)
530def _ED2(radius, name):
531 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
532 '''
533 D = Datums.Sphere
534 E = D.ellipsoid
535 if name or radius != E.a: # != E.b
536 n = _under(name)
537 E = Ellipsoid(radius, radius, name=n)
538 D = Datum(E, transform=Transforms.Identity, name=n)
539 return E, D
542def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
543 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
544 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
546 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
547 '''
548 if isinstance(earth, Datum):
549 D = earth
550 else:
551 E, n, D = _EnD3(earth, name)
552 if not E:
553 n = raiser or _earth_
554 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
555 if D is None:
556 D = Datum(E, transform=Transforms.Identity, name=_under(n))
557 if raiser and not D.isEllipsoidal:
558 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
559 return D
562def _EnD3(earth, name):
563 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
564 '''
565 D = None
566 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
567 E = earth
568 n = _under(name or E.name)
569 elif isinstance(earth, Datum):
570 E = earth.ellipsoid
571 n = _under(name or earth.name)
572 D = earth
573 elif _isRadius(earth):
574 E, D = _ED2(Radius_(earth), name)
575 n = E.name
576 elif isinstance(earth, a_f2Tuple):
577 n = _under(name or earth.name)
578 E = earth.ellipsoid(name=n)
579 elif islistuple(earth, minum=2):
580 E = Ellipsoids.Sphere
581 a, f = earth[:2]
582 if f or a != E.a: # != E.b
583 n = _under(name or _xattr(earth, name=NN))
584 E = Ellipsoid(a, f=f, name=n)
585 else:
586 n = E.name
587 D = Datums.Sphere
588 else:
589 E, n = None, NN
590 return E, n, D
593def _equall(t1, t2): # in .trf
594 '''(INTERNAL) Return L{Transform} C{t1 == t2}.
595 '''
596 return all(map(_operator.eq, t1, t2))
599def _mean_radius(radius, *lats):
600 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
601 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
602 '''
603 if radius is R_M:
604 r = radius
605 elif _isRadius(radius):
606 r = Radius_(radius, low=0, Error=TypeError)
607 else:
608 E = _ellipsoidal_datum(radius).ellipsoid
609 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
610 return r
613def _negastr(name): # in .trf
614 '''(INTERNAL) Negate a C{Transform/-Xform} name.
615 '''
616 b, m, p = _BAR_, _MINUS_, _PLUS_
617 n = name.replace(m, b).replace(p, m).replace(b, p)
618 # as good and fast as (in Python 3+ only) ...
619 # _MINUSxPLUS = str.maketrans({_MINUS_: _PLUS_, _PLUS_: _MINUS_})
620 # def _negastr(name):
621 # n = name.translate(_MINUSxPLUS)
622 # ...
623 return n.lstrip(p) if n.startswith(p) else NN(m, n)
626def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
627 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
628 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
630 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
631 '''
632 if _isRadius(earth):
633 _, D = _ED2(Radius_(earth, Error=Error), name)
634 else:
635 D = _ellipsoidal_datum(earth, Error=Error, name=name)
636 if raiser and not D.isSpherical:
637 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
638 return D
641class Datums(_NamedEnum):
642 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
643 to accommodate the L{_LazyNamedEnumItem} properties.
644 '''
645 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
646 '''(INTERNAL) Instantiate the L{Datum}.
647 '''
648 return Datum(Ellipsoids.get(ellipsoid_name),
649 Transforms.get(transform_name), name=name)
651Datums = Datums(Datum) # PYCHOK singleton
652'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
653# Datums with associated ellipsoid and Helmert transform parameters
654# to convert from WGS84 into the given datum. More are available at
655# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
656# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
657Datums._assert(
658 # Belgian Datum 1972, based on Hayford ellipsoid.
659 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
660 # <https://SpatialReference.org/ref/sr-org/7718/html/>
661 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
663 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
664 # <https://WikiPedia.org/wiki/Helmert_transformation>
665 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
667 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
668 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
670 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
671# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
672# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
673# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
674 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
676 # <https://WikiPedia.org/wiki/GRS_80>
677 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
679 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
680# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
681 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
683 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
684 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
685 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
687 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
688 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
690 # <https://WikiPedia.org/wiki/Helmert_transformation>
691 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
693 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
694 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
695 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
697 # Nouvelle Triangulation Francaise (Paris) XXX verify
698 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
700 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
701 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
703 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
704 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
706 # XXX psuedo-ellipsoids for spherical LatLon
707 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
709 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
710 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
712 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
713 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
715 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
716)
718_WGS84 = Datums.WGS84
719assert _WGS84.ellipsoid is _EWGS84
720# assert _WGS84.transform.isunity
722if __name__ == '__main__':
724 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
725 from pygeodesy.lazily import printf
727 # __doc__ of this file, force all into registery
728 for r in (Datums, Transforms):
729 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
730 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
732# **) MIT License
733#
734# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
735#
736# Permission is hereby granted, free of charge, to any person obtaining a
737# copy of this software and associated documentation files (the "Software"),
738# to deal in the Software without restriction, including without limitation
739# the rights to use, copy, modify, merge, publish, distribute, sublicense,
740# and/or sell copies of the Software, and to permit persons to whom the
741# Software is furnished to do so, subject to the following conditions:
742#
743# The above copyright notice and this permission notice shall be included
744# in all copies or substantial portions of the Software.
745#
746# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
747# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
748# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
749# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
750# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
751# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
752# OTHER DEALINGS IN THE SOFTWARE.