Coverage for pygeodesy/datums.py: 94%
223 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-01-26 16:28 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2024-01-26 16:28 -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.86863, ty=-52.29778, tz=103.72389, rx=-0, ry=-0, rz=-0.00001, s=1.2727, s1=1, sx=-0.33657, sy=-0.45696, sz=-1.84218)
47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, rx=-0.00001, ry=-0, rz=0.00001, s=-8.3, s1=0.99999, sx=-1.04, sy=-0.35, sz=3.08)
48@var Transforms.Clarke1866: Transform(name='Clarke1866', tx=8, ty=-160, tz=-176, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0)
49@var Transforms.DHDN: Transform(name='DHDN', tx=-591.28, ty=-81.35, tz=-396.39, rx=0.00001, ry=-0, rz=-0.00001, s=-9.82, s1=0.99999, sx=1.477, sy=-0.0736, sz=-1.458)
50@var Transforms.DHDNE: Transform(name='DHDNE', tx=-612.4, ty=-77, tz=-440.2, rx=0, ry=-0, rz=0.00001, s=-2.55, s1=1, 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, rx=-0, ry=-0, rz=0.00001, s=-6.7, s1=0.99999, 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, rx=0, ry=0, rz=0, s=-1.2, s1=1, sx=0, sy=0, sz=0.156)
53@var Transforms.Identity: Transform(name='Identity', tx=0, ty=0, tz=0, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0)
54@var Transforms.Irl1965: Transform(name='Irl1965', tx=-482.53, ty=130.596, tz=-564.557, rx=0.00001, ry=0, rz=0, s=-8.15, s1=0.99999, 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, rx=0.00001, ry=0, rz=0, s=-8.15, s1=0.99999, sx=1.042, sy=0.214, sz=0.631)
56@var Transforms.Krassovski1940: Transform(name='Krassovski1940', tx=-24, ty=123, tz=94, rx=-0, ry=0, rz=0, s=-2.423, s1=1, sx=-0.02, sy=0.26, sz=0.13)
57@var Transforms.Krassowsky1940: Transform(name='Krassowsky1940', tx=-24, ty=123, tz=94, rx=-0, ry=0, rz=0, s=-2.423, s1=1, 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, rx=0.00002, ry=0.00001, rz=0.00003, s=-2.423, s1=1, sx=5.137, sy=1.474, sz=5.297)
59@var Transforms.NAD27: Transform(name='NAD27', tx=8, ty=-160, tz=-176, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0)
60@var Transforms.NAD83: Transform(name='NAD83', tx=1.004, ty=-1.91, tz=-0.515, rx=0, ry=0, rz=0, s=-0.0015, s1=1, sx=0.0267, sy=0.00034, sz=0.011)
61@var Transforms.NTF: Transform(name='NTF', tx=-168, ty=-60, tz=320, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0)
62@var Transforms.OSGB36: Transform(name='OSGB36', tx=-446.448, ty=125.157, tz=-542.06, rx=-0, ry=-0, rz=-0, s=20.4894, s1=1.00002, sx=-0.1502, sy=-0.247, sz=-0.8421)
63@var Transforms.TokyoJapan: Transform(name='TokyoJapan', tx=148, ty=-507, tz=-685, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0)
64@var Transforms.WGS72: Transform(name='WGS72', tx=0, ty=0, tz=-4.5, rx=0, ry=0, rz=0, s=-0.22, s1=1, sx=0, sy=0, sz=0.554)
65@var Transforms.WGS84: Transform(name='WGS84', tx=0, ty=0, tz=0, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) # .isunity
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
72from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, \
73 _EWGS84, Vector3Tuple
74from pygeodesy.errors import _IsnotError, _TypeError, _xattr
75from pygeodesy.fmath import fdot, fmean, Fmt
76from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, _cartesian_, \
77 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \
78 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \
79 _Krassowsky1940_, _NAD27_, _NAD83_, _s_, _Sphere_, _spherical_, \
80 _sx_, _sy_, _sz_, _transform_, _tx_, _ty_, _tz_, _UNDER_, \
81 _WGS72_, _WGS84_, _under
82from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
83from pygeodesy.named import _NamedEnum, _NamedEnumItem, Property_RO, property_RO, \
84 _lazyNamedEnumItem as _lazy
85# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
86# from pygeodesy.props import Property_RO, property_RO # from .named
87# from pygeodesy.streprs import Fmt # from .fmath
88from pygeodesy.units import _isRadius, radians, Radius_
90# from math import radians # from .units
92__all__ = _ALL_LAZY.datums
93__version__ = '24.01.25'
95_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
96_BD72_ = 'BD72'
97_DHDN_ = 'DHDN'
98_DHDNE_ = 'DHDNE'
99_DHDNW_ = 'DHDNW'
100_ED50_ = 'ED50'
101_GDA2020_ = 'GDA2020'
102_Identity_ = 'Identity'
103_Inverse_ = 'Inverse'
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'
114def _rps2(s_):
115 '''(INTERNAL) Rotation in C{radians} and C{degree seconds}.
116 '''
117 return (_RPS * s_), s_
120class Transform(_NamedEnumItem):
121 '''Helmert transformation.
123 @see: L{Helmert7Tuple}.
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 (degree seconds)
137 sy = _0_0 # y rotation (degree seconds)
138 sz = _0_0 # z rotation (degree seconds)
140 def __init__(self, name=NN, tx=0, ty=0, tz=0,
141 sx=0, sy=0, sz=0, s=0):
142 '''New L{Transform}.
144 @kwarg name: Optional, unique name (C{str}).
145 @kwarg tx: Optional X translation (C{meter}).
146 @kwarg ty: Optional Y translation (C{meter}).
147 @kwarg tz: Optional Z translation (C{meter}).
148 @kwarg s: Optional scale (C{float}), ppm.
149 @kwarg sx: Optional X rotation (C{degree seconds}).
150 @kwarg sy: Optional Y rotation (C{degree seconds}).
151 @kwarg sz: Optional Z rotation (C{degree 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 sx: # secs to rads
162 self.rx, self.sx = _rps2(sx)
163 if sy:
164 self.ry, self.sy = _rps2(sy)
165 if sz:
166 self.rz, self.sz = _rps2(sz)
167 if s:
168 self.s = s
169 self.s1 = _F(s * 1e-6 + _1_0) # normalize ppm to (s + 1)
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 self.tx == other.tx
182 and self.ty == other.ty
183 and self.tz == other.tz
184 and self.rx == other.rx
185 and self.ry == other.ry
186 and self.rz == other.rz
187 and self.s == other.s)
189 def __hash__(self):
190 return self._hash # memoized
192 def __matmul__(self, other): # PYCHOK Python 3.5+
193 '''Helmert-transform a cartesian B{C{other}}.
195 @raise TypeError: Invalid B{C{other}}.
196 '''
197 try: # only CartesianBase
198 return other.toTransform(self)
199 except AttributeError:
200 pass
201 raise _IsnotError(_cartesian_, other=other)
203 @Property_RO
204 def _hash(self):
205 return hash((self.tx, self.ty, self.tz,
206 self.rx, self.ry, self.rz, self.s))
208 def inverse(self, name=NN):
209 '''Return the inverse of this transform.
211 @kwarg name: Optional, unique name (C{str}).
213 @return: Inverse (Transform).
215 @raise NameError: Transform with that B{C{name}} already exists.
216 '''
217 return Transform(name=name or (self.name + _Inverse_),
218 tx=-self.tx, ty=-self.ty, tz=-self.tz,
219 sx=-self.sx, sy=-self.sy, sz=-self.sz, s=-self.s)
221 @Property_RO
222 def isunity(self):
223 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84?
224 '''
225 return not any((self.tx, self.ty, self.tz,
226 self.rx, self.ry, self.rz, self.s))
228 def toStr(self, prec=5, name=NN, **unused): # PYCHOK expected
229 '''Return this transform as a string.
231 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
232 @kwarg name: Override name (C{str}) or C{None} to exclude
233 this transform's name.
235 @return: Transform attributes (C{str}).
236 '''
237 return self._instr(name, prec, _tx_, _ty_, _tz_,
238 'rx', 'ry', 'rz', _s_, 's1',
239 _sx_, _sy_, _sz_)
241 def transform(self, x, y, z, inverse=False):
242 '''Transform a (geocentric) Cartesian point, forward or inverse.
244 @arg x: X coordinate (C{meter}).
245 @arg y: Y coordinate (C{meter}).
246 @arg z: Z coordinate (C{meter}).
247 @kwarg inverse: Optional direction, forward or inverse (C{bool}).
249 @return: A L{Vector3Tuple}C{(x, y, z)}, transformed.
250 '''
251 if self.isunity:
252 return Vector3Tuple(x, y, z, name=self.name)
254 xyz1 = x, y, z, _1_0
255 s1 = self.s1
256 if inverse:
257 xyz1 = map2(neg, xyz1)
258 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1
259 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
260 # x * .rz + y * .s1 - z * .rx + .ty,
261 # -x * .ry + y * .rx + z * .s1 + .tz)
262 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
263 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
264 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
265 name=self.name)
268class Transforms(_NamedEnum):
269 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
270 to accommodate the L{_LazyNamedEnumItem} properties.
271 '''
272 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
273 '''(INTERNAL) Instantiate the C{Transform}.
274 '''
275 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
277Transforms = Transforms(Transform) # PYCHOK singleton
278'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
279# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...
280Transforms._assert(
281 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
282 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
283 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
284 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
285 s=_F( 1.2727)),
286 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
287 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
288 s=_F( -8.3)),
289 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
290 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
291 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
292 s=_F( -9.82)), # Germany
293 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2),
294 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>
295 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797),
296 s=_F( -2.55)), # East Germany
297 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2),
298 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>
299 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455),
300 s=_F( -6.7)), # West Germany
301 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
302 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
303 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
304 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
305 sz=_F( 0.156), s=_F(-1.2)),
306 Identity = _lazy(_Identity_),
307 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
308 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>
309 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
310 s=_F( -8.15)),
311 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
312 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>
313 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
314 s=_F( -8.15)),
315 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
316 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
317 s=_F( -2.423)), # spelling
318 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
319 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
320 s=_F( -2.423)), # spelling
321 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
322 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
323 s=_F( -2.423)), # Austria
324 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
325 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
326 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
327 s=_F(-0.0015)),
328 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
329 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
330 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>
331 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
332 s=_F( 20.4894)),
333 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
334 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
335 WGS84 = _lazy(_WGS84_), # unity
336)
339class Datum(_NamedEnumItem):
340 '''Ellipsoid and transform parameters for an earth model.
341 '''
342 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
343 _transform = Transforms.WGS84 # default transform (L{Transform})
345 def __init__(self, ellipsoid, transform=None, name=NN):
346 '''New L{Datum}.
348 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
349 @kwarg transform: Optional transform (L{Transform}).
350 @kwarg name: Optional, unique name (C{str}).
352 @raise NameError: Datum with that B{C{name}} already exists.
354 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
355 nor L{Ellipsoid2} or B{C{transform}} is
356 not a L{Transform}.
357 '''
358 self._ellipsoid = ellipsoid or Datum._ellipsoid
359 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
361 self._transform = transform or Datum._transform
362 _xinstanceof(Transform, transform=self.transform)
364 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
366 def __eq__(self, other):
367 '''Compare this and an other datum.
369 @arg other: The other datum (L{Datum}).
371 @return: C{True} if equal, C{False} otherwise.
372 '''
373 return self is other or (isinstance(other, Datum) and
374 self.ellipsoid == other.ellipsoid and
375 self.transform == other.transform)
377 def __hash__(self):
378 return self._hash # memoized
380 def __matmul__(self, other): # PYCHOK Python 3.5+
381 '''Convert cartesian or ellipsoidal B{C{other}} to this datum.
383 @raise TypeError: Invalid B{C{other}}.
384 '''
385 try: # only CartesianBase and LatLonEllipsoidalBase
386 return other.toDatum(self)
387 except AttributeError:
388 pass
389 raise _IsnotError(_cartesian_, _ellipsoidal_, other=other)
391 def ecef(self, Ecef=None):
392 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
394 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
396 @return: An ECEF converter for this C{datum}.
398 @raise TypeError: Invalid B{C{Ecef}}.
400 @see: Module L{pygeodesy.ecef}.
401 '''
402 return _MODS.ecef._4Ecef(self, Ecef)
404 @Property_RO
405 def ellipsoid(self):
406 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
407 '''
408 return self._ellipsoid
410 @Property_RO
411 def exactTM(self):
412 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
413 '''
414 return _MODS.etm.ExactTransverseMercator(datum=self)
416 @Property_RO
417 def _hash(self):
418 return hash(self.ellipsoid) + hash(self.transform)
420 @property_RO
421 def isEllipsoidal(self):
422 '''Check whether this datum is ellipsoidal (C{bool}).
423 '''
424 return self.ellipsoid.isEllipsoidal
426 @property_RO
427 def isOblate(self):
428 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
429 '''
430 return self.ellipsoid.isOblate
432 @property_RO
433 def isProlate(self):
434 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
435 '''
436 return self.ellipsoid.isProlate
438 @property_RO
439 def isSpherical(self):
440 '''Check whether this datum is (near-)spherical (C{bool}).
441 '''
442 return self.ellipsoid.isSpherical
444 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
445 '''Return this datum as a string.
447 @kwarg sep: Separator to join (C{str}).
448 @kwarg name: Override name (C{str}) or C{None} to exclude
449 this datum's name.
451 @return: Datum attributes (C{str}).
452 '''
453 t = [] if name is None else \
454 [Fmt.EQUAL(name=repr(name or self.named))]
455 for a in (_ellipsoid_, _transform_):
456 v = getattr(self, a)
457 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
458 return sep.join(t)
460 @Property_RO
461 def transform(self):
462 '''Get this datum's transform (L{Transform}).
463 '''
464 return self._transform
467def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
468 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}
469 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).
471 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.
472 '''
473 if f is not None:
474 E, n, D = _EnD3((a_earth, f), name)
475 if raiser and not E:
476 raise _TypeError(f=f, **{raiser: a_earth})
477 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None):
478 return
479 elif isinstance(a_earth, Datum):
480 E, n, D = None, NN, a_earth
481 else:
482 E, n, D = _EnD3(a_earth, name)
483 if raiser and not E:
484 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})
485 if D is None:
486 D = Datum(E, transform=Transforms.Identity, name=_under(n))
487 inst._datum = D
490def _earth_ellipsoid(earth, *name_raiser):
491 '''(INTERAL) Return the ellipsoid for the given C{earth} model.
492 '''
493 return Ellipsoids.Sphere if earth is R_M else (
494 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else
495 _spherical_datum(earth, *name_raiser).ellipsoid)
498def _ED2(radius, name):
499 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.
500 '''
501 D = Datums.Sphere
502 E = D.ellipsoid
503 if name or radius != E.a: # != E.b
504 n = _under(name)
505 E = Ellipsoid(radius, radius, name=n)
506 D = Datum(E, transform=Transforms.Identity, name=n)
507 return E, D
510def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
511 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
512 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
514 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
515 '''
516 if isinstance(earth, Datum):
517 D = earth
518 else:
519 E, n, D = _EnD3(earth, name)
520 if not E:
521 n = raiser or _earth_
522 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
523 if D is None:
524 D = Datum(E, transform=Transforms.Identity, name=_under(n))
525 if raiser and not D.isEllipsoidal:
526 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
527 return D
530def _EnD3(earth, name):
531 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.
532 '''
533 D = None
534 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
535 E = earth
536 n = _under(name or E.name)
537 elif isinstance(earth, Datum):
538 E = earth.ellipsoid
539 n = _under(name or earth.name)
540 D = earth
541 elif _isRadius(earth):
542 E, D = _ED2(Radius_(earth), name)
543 n = E.name
544 elif isinstance(earth, a_f2Tuple):
545 n = _under(name or earth.name)
546 E = earth.ellipsoid(name=n)
547 elif islistuple(earth, minum=2):
548 E = Ellipsoids.Sphere
549 a, f = earth[:2]
550 if f or a != E.a: # != E.b
551 n = _under(name or _xattr(earth, name=NN))
552 E = Ellipsoid(a, f=f, name=n)
553 else:
554 n = E.name
555 D = Datums.Sphere
556 else:
557 E, n = None, NN
558 return E, n, D
561def _mean_radius(radius, *lats):
562 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
563 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
564 '''
565 if radius is R_M:
566 r = radius
567 elif _isRadius(radius):
568 r = Radius_(radius, low=0, Error=TypeError)
569 else:
570 E = _ellipsoidal_datum(radius).ellipsoid
571 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
572 return r
575def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
576 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
577 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
579 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
580 '''
581 if _isRadius(earth):
582 _, D = _ED2(Radius_(earth, Error=Error), name)
583 else:
584 D = _ellipsoidal_datum(earth, Error=Error, name=name)
585 if raiser and not D.isSpherical:
586 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
587 return D
590class Datums(_NamedEnum):
591 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
592 to accommodate the L{_LazyNamedEnumItem} properties.
593 '''
594 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
595 '''(INTERNAL) Instantiate the L{Datum}.
596 '''
597 return Datum(Ellipsoids.get(ellipsoid_name),
598 Transforms.get(transform_name), name=name)
600Datums = Datums(Datum) # PYCHOK singleton
601'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
602# Datums with associated ellipsoid and Helmert transform parameters
603# to convert from WGS84 into the given datum. More are available at
604# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
605# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
606Datums._assert(
607 # Belgian Datum 1972, based on Hayford ellipsoid.
608 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
609 # <https://SpatialReference.org/ref/sr-org/7718/html/>
610 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
612 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
613 # <https://WikiPedia.org/wiki/Helmert_transformation>
614 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
616 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
617 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
619 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
620# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
621# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
622# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
623 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
625 # <https://WikiPedia.org/wiki/GRS_80>
626 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
628 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
629# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
630 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
632 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
633 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
634 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
636 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
637 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
639 # <https://WikiPedia.org/wiki/Helmert_transformation>
640 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
642 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
643 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
644 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
646 # Nouvelle Triangulation Francaise (Paris) XXX verify
647 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
649 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
650 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
652 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
653 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
655 # XXX psuedo-ellipsoids for spherical LatLon
656 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
658 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
659 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
661 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
662 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
664 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
665)
667_WGS84 = Datums.WGS84
668assert _WGS84.ellipsoid is _EWGS84
669# assert _WGS84.transform.isunity
671if __name__ == '__main__':
673 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
674 from pygeodesy.lazily import printf
676 # __doc__ of this file, force all into registery
677 for r in (Datums, Transforms):
678 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
679 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
681# **) MIT License
682#
683# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
684#
685# Permission is hereby granted, free of charge, to any person obtaining a
686# copy of this software and associated documentation files (the "Software"),
687# to deal in the Software without restriction, including without limitation
688# the rights to use, copy, modify, merge, publish, distribute, sublicense,
689# and/or sell copies of the Software, and to permit persons to whom the
690# Software is furnished to do so, subject to the following conditions:
691#
692# The above copyright notice and this permission notice shall be included
693# in all copies or substantial portions of the Software.
694#
695# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
696# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
697# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
698# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
699# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
700# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
701# OTHER DEALINGS IN THE SOFTWARE.