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