Coverage for pygeodesy/datums.py: 94%
241 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-02-25 12:57 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2024-02-25 12:57 -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, _xinstanceof
71from pygeodesy.constants import R_M, _float as _F, _0_0, _0_26, _1_0, _2_0, _8_0, _3600_0
72# from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase as _CEB, \
73# LatLonEllipsoidalBase as _LLEB # MODS
74from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, _EWGS84, \
75 Vector3Tuple
76from pygeodesy.errors import _IsnotError, _TypeError, _xattr, _xellipsoidall, _xkwds, \
77 _xkwds_pop2
78from pygeodesy.fmath import fdot, fmean, Fmt
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_, _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
94__all__ = _ALL_LAZY.datums
95__version__ = '24.02.24'
97_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
98_BD72_ = 'BD72'
99_DHDN_ = 'DHDN'
100_DHDNE_ = 'DHDNE'
101_DHDNW_ = 'DHDNW'
102_ED50_ = 'ED50'
103_GDA2020_ = 'GDA2020' # in .trf
104_Identity_ = 'Identity'
105_Irl1965_ = 'Irl1965'
106_Irl1975_ = 'Irl1975'
107_MGI_ = 'MGI'
108_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf
109_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:]
110_NTF_ = 'NTF'
111_OSGB36_ = 'OSGB36'
112_Potsdam_ = 'Potsdam'
113_RPS = radians(_1_0 / _3600_0) # radians per degree-second
114_S1_S = 1.e-6 # in .trf
115_TokyoJapan_ = 'TokyoJapan'
118class Transform(_NamedEnumItem):
119 '''Helmert I{datum} transformation.
121 @see: L{TransformXform<trf.TransformXform>}.
122 '''
123 tx = _0_0 # x translation (C{meter})
124 ty = _0_0 # y translation (C{meter})
125 tz = _0_0 # z translation (C{meter})
127 rx = _0_0 # x rotation (C{radians})
128 ry = _0_0 # y rotation (C{radians})
129 rz = _0_0 # z rotation (C{radians})
131 s = _0_0 # scale ppm (C{float})
132 s1 = _1_0 # scale + 1 (C{float})
134 sx = _0_0 # x rotation (C{degree-seconds})
135 sy = _0_0 # y rotation (C{degree-seconds})
136 sz = _0_0 # z rotation (C{degree-seconds})
138 def __init__(self, name=NN, tx=0, ty=0, tz=0,
139 sx=0, sy=0, sz=0, s=0):
140 '''New L{Transform}.
142 @kwarg name: Optional, unique name (C{str}).
143 @kwarg tx: X translation (C{meter}).
144 @kwarg ty: Y translation (C{meter}).
145 @kwarg tz: Z translation (C{meter}).
146 @kwarg s: Scale (C{float}), ppm.
147 @kwarg sx: X rotation (C{degree-seconds}).
148 @kwarg sy: Y rotation (C{degree-seconds}).
149 @kwarg sz: Z rotation (C{degree-seconds}).
151 @raise NameError: Transform with that B{C{name}} already exists.
152 '''
153 if tx:
154 self.tx = tx
155 if ty:
156 self.ty = ty
157 if tz:
158 self.tz = tz
159 if sx: # secs to rads
160 self.rx, self.sx = self._rps2(sx)
161 if sy:
162 self.ry, self.sy = self._rps2(sy)
163 if sz:
164 self.rz, self.sz = self._rps2(sz)
165 if s:
166 self.s = s
167 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1)
169 self._register(Transforms, name)
171 def __eq__(self, other):
172 '''Compare this and an other transform.
174 @arg other: The other transform (L{Transform}).
176 @return: C{True} if equal, C{False} otherwise.
177 '''
178 return self is other or (isinstance(other, Transform) and all(
179 a == b for a, b in zip(self, other)))
181 def __hash__(self):
182 return hash(x for x in self)
184 def __iter__(self):
185 '''Yield the initial attribute values.
186 '''
187 for _, x in self.items():
188 yield x
190 def __matmul__(self, point): # PYCHOK Python 3.5+
191 '''Transform an I{ellipsoidal} B{C{point}} with this Helmert.
193 @return: A transformed copy of B{C{point}}.
195 @raise TypeError: Invalid B{C{point}}.
197 @see: Method C{B{point}.toTransform}.
198 '''
199 _ = _xellipsoidall(point)
200 return point.toTransform(self)
202 def __neg__(self):
203 return self.inverse()
205# def __sub__(self, other):
206# _xinstanceof(Transform, other=other)
207#
208# def _sub(a, b):
209# for n in _Names11:
210# yield getattr(a, n) - getattr(b, n)
211#
212# return type(self)(_sub(self, other), name=_MINUS_) # .fsums._sub_op
214# @Property_RO
215# def _hash(self):
216# return hash(x for x in self)
218 def items(self):
219 '''Yield each attribute as 2-tuple C{(name, value)}.
220 '''
221 for n in _Names7:
222 yield n, getattr(self, n)
224 def inverse(self, name=NN):
225 '''Return the inverse of this transform.
227 @kwarg name: Optional, unique name (C{str}).
229 @return: Inverse (L{Transform}).
231 @raise NameError: Transform with that B{C{name}} already exists.
232 '''
233 d = dict((n, neg(v)) for n, v in self.items())
234 n = name or _minus(self.name)
235 return type(self)(name=n, **d)
237 @Property_RO
238 def isunity(self):
239 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84?
240 '''
241 return not any(s for s in self)
243 def _rps2(self, s_):
244 '''(INTERNAL) Rotation in C{radians} and C{degree-seconds}.
245 '''
246 return (_RPS * s_), s_
248 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected
249 '''Return this transform as a string.
251 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
252 @kwarg fmt: Optional C{float} format (C{letter}).
253 @kwarg name: Override name (C{str}) or C{None} to exclude
254 this transform's name.
256 @return: Transform attributes (C{str}).
257 '''
258 return self._instr(name, prec, *_Names11, fmt=fmt)
260 def transform(self, x, y, z, inverse=False, **Vector_and_kwds):
261 '''Transform a (geocentric) position, forward or inverse.
263 @arg x: X coordinate (C{meter}).
264 @arg y: Y coordinate (C{meter}).
265 @arg z: Z coordinate (C{meter}).
266 @kwarg inverse: If C{True}, apply the inverse transform (C{bool}).
267 @kwarg Vector_and_kwds: An optional, (3-D) C{B{Vector}=None} or
268 cartesian class and additional C{B{Vector}}
269 keyword arguments to return the transformed
270 point.
272 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)})
273 unless some B{C{Vector_and_kwds}} is specified.
274 '''
275 if self.isunity:
276 r = Vector3Tuple(x, y, z, name=self.name)
277 else:
278 xyz1 = x, y, z, _1_0
279 s1 = self.s1
280 if inverse:
281 xyz1 = map2(neg, xyz1)
282 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1
283 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
284 # x * .rz + y * .s1 - z * .rx + .ty,
285 # -x * .ry + y * .rx + z * .s1 + .tz)
286 r = Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
287 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
288 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
289 name=self.name)
290 if Vector_and_kwds:
291 V, kwds = _xkwds_pop2(Vector_and_kwds, Vector=None)
292 if V:
293 r = V(r, **_xkwds(kwds, name=self.name))
294 return r
297class Transforms(_NamedEnum):
298 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
299 to accommodate the L{_LazyNamedEnumItem} properties.
300 '''
301 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
302 '''(INTERNAL) Instantiate the C{Transform}.
303 '''
304 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
306Transforms = Transforms(Transform) # PYCHOK singleton
307'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
308# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
309Transforms._assert(
310 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
311 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
312 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
313 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
314 s=_F( 1.2727)),
315 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
316 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
317 s=_F( -8.3)),
318 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
319 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
320 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
321 s=_F( -9.82)), # Germany
322 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2),
323 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
324 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797),
325 s=_F( -2.55)), # East Germany
326 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2),
327 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
328 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455),
329 s=_F( -6.7)), # West Germany
330 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
331 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
332 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
333 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
334 sz=_F( 0.156), s=_F(-1.2)),
335 Identity = _lazy(_Identity_),
336 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
337 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
338 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
339 s=_F( -8.15)),
340 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
341 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
342 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
343 s=_F( -8.15)),
344 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
345 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
346 s=_F( -2.423)), # spelling
347 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
348 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
349 s=_F( -2.423)), # spelling
350 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
351 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
352 s=_F( -2.423)), # Austria
353 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
354 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
355 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
356 s=_F(-0.0015)),
357 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
358 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
359 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
360 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
361 s=_F( 20.4894)),
362 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
363 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
364 WGS84 = _lazy(_WGS84_), # unity
365)
368class Datum(_NamedEnumItem):
369 '''Ellipsoid and transform parameters for an earth model.
370 '''
371 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
372 _transform = Transforms.WGS84 # default transform (L{Transform})
374 def __init__(self, ellipsoid, transform=None, name=NN):
375 '''New L{Datum}.
377 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
378 @kwarg transform: Optional transform (L{Transform}).
379 @kwarg name: Optional, unique name (C{str}).
381 @raise NameError: Datum with that B{C{name}} already exists.
383 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
384 nor L{Ellipsoid2} or B{C{transform}} is
385 not a L{Transform}.
386 '''
387 self._ellipsoid = ellipsoid or Datum._ellipsoid
388 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
390 self._transform = transform or Datum._transform
391 _xinstanceof(Transform, transform=self.transform)
393 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
395 def __eq__(self, other):
396 '''Compare this and an other datum.
398 @arg other: The other datum (L{Datum}).
400 @return: C{True} if equal, C{False} otherwise.
401 '''
402 return self is other or (isinstance(other, Datum) and
403 self.ellipsoid == other.ellipsoid and
404 self.transform == other.transform)
406 def __hash__(self):
407 return self._hash # memoized
409 def __matmul__(self, point): # PYCHOK Python 3.5+
410 '''Convert an I{ellipsoidal} B{C{point}} to this datum.
412 @raise TypeError: Invalid B{C{point}}.
413 '''
414 _ = _xellipsoidall(point)
415 return point.toDatum(self)
417 def ecef(self, Ecef=None):
418 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
420 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
422 @return: An ECEF converter for this C{datum}.
424 @raise TypeError: Invalid B{C{Ecef}}.
426 @see: Module L{pygeodesy.ecef}.
427 '''
428 return _MODS.ecef._4Ecef(self, Ecef)
430 @Property_RO
431 def ellipsoid(self):
432 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
433 '''
434 return self._ellipsoid
436 @Property_RO
437 def exactTM(self):
438 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
439 '''
440 return _MODS.etm.ExactTransverseMercator(datum=self)
442 @Property_RO
443 def _hash(self):
444 return hash(self.ellipsoid) + hash(self.transform)
446 @property_RO
447 def isEllipsoidal(self):
448 '''Check whether this datum is ellipsoidal (C{bool}).
449 '''
450 return self.ellipsoid.isEllipsoidal
452 @property_RO
453 def isOblate(self):
454 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
455 '''
456 return self.ellipsoid.isOblate
458 @property_RO
459 def isProlate(self):
460 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
461 '''
462 return self.ellipsoid.isProlate
464 @property_RO
465 def isSpherical(self):
466 '''Check whether this datum is (near-)spherical (C{bool}).
467 '''
468 return self.ellipsoid.isSpherical
470 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
471 '''Return this datum as a string.
473 @kwarg sep: Separator to join (C{str}).
474 @kwarg name: Override name (C{str}) or C{None} to exclude
475 this datum's name.
477 @return: Datum attributes (C{str}).
478 '''
479 t = [] if name is None else \
480 [Fmt.EQUAL(name=repr(name or self.named))]
481 for a in (_ellipsoid_, _transform_):
482 v = getattr(self, a)
483 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
484 return sep.join(t)
486 @Property_RO
487 def transform(self):
488 '''Get this datum's transform (L{Transform}).
489 '''
490 return self._transform
493def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
494 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
495 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
497 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
498 '''
499 if f is not None:
500 E, n, D = _EnD3((a_earth, f), name)
501 if raiser and not E:
502 raise _TypeError(f=f, **{raiser: a_earth})
503 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None):
504 return
505 elif isinstance(a_earth, Datum):
506 E, n, D = None, NN, a_earth
507 else:
508 E, n, D = _EnD3(a_earth, name)
509 if raiser and not E:
510 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
511 if D is None:
512 D = Datum(E, transform=Transforms.Identity, name=_under(n))
513 inst._datum = D
516def _earth_ellipsoid(earth, *name_raiser):
517 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
518 '''
519 return Ellipsoids.Sphere if earth is R_M else (
520 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
521 _spherical_datum(earth, *name_raiser).ellipsoid)
524def _ED2(radius, name):
525 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
526 '''
527 D = Datums.Sphere
528 E = D.ellipsoid
529 if name or radius != E.a: # != E.b
530 n = _under(name)
531 E = Ellipsoid(radius, radius, name=n)
532 D = Datum(E, transform=Transforms.Identity, name=n)
533 return E, D
536def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
537 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
538 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
540 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
541 '''
542 if isinstance(earth, Datum):
543 D = earth
544 else:
545 E, n, D = _EnD3(earth, name)
546 if not E:
547 n = raiser or _earth_
548 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
549 if D is None:
550 D = Datum(E, transform=Transforms.Identity, name=_under(n))
551 if raiser and not D.isEllipsoidal:
552 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
553 return D
556def _EnD3(earth, name):
557 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
558 '''
559 D = None
560 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
561 E = earth
562 n = _under(name or E.name)
563 elif isinstance(earth, Datum):
564 E = earth.ellipsoid
565 n = _under(name or earth.name)
566 D = earth
567 elif _isRadius(earth):
568 E, D = _ED2(Radius_(earth), name)
569 n = E.name
570 elif isinstance(earth, a_f2Tuple):
571 n = _under(name or earth.name)
572 E = earth.ellipsoid(name=n)
573 elif islistuple(earth, minum=2):
574 E = Ellipsoids.Sphere
575 a, f = earth[:2]
576 if f or a != E.a: # != E.b
577 n = _under(name or _xattr(earth, name=NN))
578 E = Ellipsoid(a, f=f, name=n)
579 else:
580 n = E.name
581 D = Datums.Sphere
582 else:
583 E, n = None, NN
584 return E, n, D
587def _mean_radius(radius, *lats):
588 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
589 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
590 '''
591 if radius is R_M:
592 r = radius
593 elif _isRadius(radius):
594 r = Radius_(radius, low=0, Error=TypeError)
595 else:
596 E = _ellipsoidal_datum(radius).ellipsoid
597 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
598 return r
601def _minus(name): # in .trf
602 '''(INTERNAL) Name of C{inverse} Xform.
603 '''
604 m = _MINUS_
605 return name[len(m):] if name.startswith(m) else NN(m, name)
608def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
609 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
610 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
612 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
613 '''
614 if _isRadius(earth):
615 _, D = _ED2(Radius_(earth, Error=Error), name)
616 else:
617 D = _ellipsoidal_datum(earth, Error=Error, name=name)
618 if raiser and not D.isSpherical:
619 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
620 return D
623class Datums(_NamedEnum):
624 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
625 to accommodate the L{_LazyNamedEnumItem} properties.
626 '''
627 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
628 '''(INTERNAL) Instantiate the L{Datum}.
629 '''
630 return Datum(Ellipsoids.get(ellipsoid_name),
631 Transforms.get(transform_name), name=name)
633Datums = Datums(Datum) # PYCHOK singleton
634'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
635# Datums with associated ellipsoid and Helmert transform parameters
636# to convert from WGS84 into the given datum. More are available at
637# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
638# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
639Datums._assert(
640 # Belgian Datum 1972, based on Hayford ellipsoid.
641 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
642 # <https://SpatialReference.org/ref/sr-org/7718/html/>
643 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
645 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
646 # <https://WikiPedia.org/wiki/Helmert_transformation>
647 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
649 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
650 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
652 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
653# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
654# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
655# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
656 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
658 # <https://WikiPedia.org/wiki/GRS_80>
659 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
661 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
662# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
663 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
665 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
666 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
667 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
669 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
670 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
672 # <https://WikiPedia.org/wiki/Helmert_transformation>
673 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
675 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
676 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
677 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
679 # Nouvelle Triangulation Francaise (Paris) XXX verify
680 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
682 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
683 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
685 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
686 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
688 # XXX psuedo-ellipsoids for spherical LatLon
689 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
691 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
692 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
694 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
695 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
697 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
698)
700_WGS84 = Datums.WGS84
701assert _WGS84.ellipsoid is _EWGS84
702# assert _WGS84.transform.isunity
704if __name__ == '__main__':
706 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
707 from pygeodesy.lazily import printf
709 # __doc__ of this file, force all into registery
710 for r in (Datums, Transforms):
711 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
712 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
714# **) MIT License
715#
716# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
717#
718# Permission is hereby granted, free of charge, to any person obtaining a
719# copy of this software and associated documentation files (the "Software"),
720# to deal in the Software without restriction, including without limitation
721# the rights to use, copy, modify, merge, publish, distribute, sublicense,
722# and/or sell copies of the Software, and to permit persons to whom the
723# Software is furnished to do so, subject to the following conditions:
724#
725# The above copyright notice and this permission notice shall be included
726# in all copies or substantial portions of the Software.
727#
728# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
729# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
730# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
731# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
732# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
733# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
734# OTHER DEALINGS IN THE SOFTWARE.