Coverage for pygeodesy/datums.py: 95%
194 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -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 or
15above/below the earth’s surface, measured in degrees from the equator, from the International
16Reference Meridian, in meters above the ellipsoid and based on a given datum. The datum in turn
17is based on a reference ellipsoid and tied to geodetic survey reference points.
19Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems),
20but previously various reference ellipsoids and datum references were used.
22The UK Ordnance Survey National Grid References are still based on the otherwise historical
23OSGB36 datum, q.v. U{"A Guide to Coordinate Systems in Great Britain", Section 6
24<https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>}.
26@var Datums.BD72: Datum(name='BD72', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.BD72)
27@var Datums.DHDN: Datum(name='DHDN', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.DHDN)
28@var Datums.ED50: Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50)
29@var Datums.GDA2020: Datum(name='GDA2020', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)
30@var Datums.GRS80: Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)
31@var Datums.Irl1975: Datum(name='Irl1975', ellipsoid=Ellipsoids.AiryModified, transform=Transforms.Irl1975)
32@var Datums.Krassovski1940: Datum(name='Krassovski1940', ellipsoid=Ellipsoids.Krassovski1940, transform=Transforms.Krassovski1940)
33@var Datums.Krassowsky1940: Datum(name='Krassowsky1940', ellipsoid=Ellipsoids.Krassowsky1940, transform=Transforms.Krassowsky1940)
34@var Datums.MGI: Datum(name='MGI', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.MGI)
35@var Datums.NAD27: Datum(name='NAD27', ellipsoid=Ellipsoids.Clarke1866, transform=Transforms.NAD27)
36@var Datums.NAD83: Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83)
37@var Datums.NTF: Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF)
38@var Datums.OSGB36: Datum(name='OSGB36', ellipsoid=Ellipsoids.Airy1830, transform=Transforms.OSGB36)
39@var Datums.Potsdam: Datum(name='Potsdam', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841)
40@var Datums.Sphere: Datum(name='Sphere', ellipsoid=Ellipsoids.Sphere, transform=Transforms.WGS84)
41@var Datums.TokyoJapan: Datum(name='TokyoJapan', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.TokyoJapan)
42@var Datums.WGS72: Datum(name='WGS72', ellipsoid=Ellipsoids.WGS72, transform=Transforms.WGS72)
43@var Datums.WGS84: Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84)
45@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)
46@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)
47@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)
48@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)
49@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)
50@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)
51@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)
52@var Transforms.Irl1975: Transform(name='Irl1975', tx=-482.53, ty=130.596, tz=-564.557, rx=-0.00001, ry=-0, rz=-0, s=-1.1, s1=1, sx=-1.042, sy=-0.214, sz=-0.631)
53@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)
54@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)
55@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)
56@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)
57@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)
58@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)
59@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)
60@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)
61@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)
62@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)
63'''
64# make sure int/int division yields float quotient, see .basics
65from __future__ import division as _; del _ # PYCHOK semicolon
67from pygeodesy.basics import islistuple, isscalar, map2, neg, _xinstanceof
68from pygeodesy.constants import R_M, _float as _F, _0_0, _0_26, _1_0, _2_0, _8_0, _3600_0
69from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, Vector3Tuple
70from pygeodesy.errors import _IsnotError, _xattr
71from pygeodesy.fmath import fdot, fmean, Fmt
72from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, _cartesian_, \
73 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \
74 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \
75 _Krassowsky1940_, _NAD27_, _NAD83_, _s_, _Sphere_, _spherical_, \
76 _sx_, _sy_, _sz_, _transform_, _tx_, _ty_, _tz_, _UNDER_, \
77 _WGS72_, _WGS84_, _under
78from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
79from pygeodesy.named import _NamedEnum, _NamedEnumItem, \
80 _lazyNamedEnumItem as _lazy, Property_RO
81# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
82# from pygeodesy.props import Property_RO # from .named
83# from pygeodesy.streprs import Fmt # from .fmath
84from pygeodesy.units import radians, Radius_
86# from math import radians # from .units
88__all__ = _ALL_LAZY.datums
89__version__ = '23.08.05'
91_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
92_BD72_ = 'BD72'
93_DHDN_ = 'DHDN'
94_ED50_ = 'ED50'
95_GDA2020_ = 'GDA2020'
96_Identity_ = 'Identity'
97_Inverse_ = 'Inverse'
98_Irl1965_ = 'Irl1965'
99_Irl1975_ = 'Irl1975'
100_MGI_ = 'MGI'
101_NTF_ = 'NTF'
102_OSGB36_ = 'OSGB36'
103_Potsdam_ = 'Potsdam'
104_TokyoJapan_ = 'TokyoJapan'
106_r_s1 = radians(1 / _3600_0) # 1 degree second to radians
109def _r_s2(s_):
110 '''(INTERNAL) rotation in C{radians} and C{degree seconds}.
111 '''
112 return _F(_r_s1 * s_), s_
115class Transform(_NamedEnumItem):
116 '''Helmert transformation.
118 @see: L{Helmert7Tuple}.
119 '''
120 tx = _0_0 # x translation (C{meter})
121 ty = _0_0 # y translation (C{meter})
122 tz = _0_0 # z translation (C{meter})
124 rx = _0_0 # x rotation (C{radians})
125 ry = _0_0 # y rotation (C{radians})
126 rz = _0_0 # z rotation (C{radians})
128 s = _0_0 # scale ppm (C{float})
129 s1 = _1_0 # scale + 1 (C{float})
131 sx = _0_0 # x rotation (degree seconds)
132 sy = _0_0 # y rotation (degree seconds)
133 sz = _0_0 # z rotation (degree seconds)
135 def __init__(self, name=NN, tx=0, ty=0, tz=0,
136 sx=0, sy=0, sz=0, s=0):
137 '''New L{Transform}.
139 @kwarg name: Optional, unique name (C{str}).
140 @kwarg tx: Optional X translation (C{meter}).
141 @kwarg ty: Optional Y translation (C{meter}).
142 @kwarg tz: Optional Z translation (C{meter}).
143 @kwarg s: Optional scale (C{float}), ppm.
144 @kwarg sx: Optional X rotation (C{degree seconds}).
145 @kwarg sy: Optional Y rotation (C{degree seconds}).
146 @kwarg sz: Optional Z rotation (C{degree seconds}).
148 @raise NameError: Transform with that B{C{name}} already exists.
149 '''
150 if tx:
151 self.tx = tx
152 if ty:
153 self.ty = ty
154 if tz:
155 self.tz = tz
156 if sx: # secs to rads
157 self.rx, self.sx = _r_s2(sx)
158 if sy:
159 self.ry, self.sy = _r_s2(sy)
160 if sz:
161 self.rz, self.sz = _r_s2(sz)
162 if s:
163 self.s = s
164 self.s1 = _F(s * 1e-6 + _1_0) # normalize ppm to (s + 1)
166 self._register(Transforms, name)
168 def __eq__(self, other):
169 '''Compare this and an other transform.
171 @arg other: The other transform (L{Transform}).
173 @return: C{True} if equal, C{False} otherwise.
174 '''
175 return self is other or (isinstance(other, Transform)
176 and self.tx == other.tx
177 and self.ty == other.ty
178 and self.tz == other.tz
179 and self.rx == other.rx
180 and self.ry == other.ry
181 and self.rz == other.rz
182 and self.s == other.s)
184 def __hash__(self):
185 return self._hash # memoized
187 def __matmul__(self, other): # PYCHOK Python 3.5+
188 '''Helmert-transform a cartesian B{C{other}}.
190 @raise TypeError: Invalid B{C{other}}.
191 '''
192 try: # only CartesianBase
193 return other.toTransform(self)
194 except AttributeError:
195 pass
196 raise _IsnotError(_cartesian_, other=other)
198 @Property_RO
199 def _hash(self):
200 return hash((self.rx, self.ry, self.rz, self.s,
201 self.tx, self.ty, self.tz))
203 def inverse(self, name=NN):
204 '''Return the inverse of this transform.
206 @kwarg name: Optional, unique name (C{str}).
208 @return: Inverse (Transform).
210 @raise NameError: Transform with that B{C{name}} already exists.
211 '''
212 return Transform(name=name or (self.name + _Inverse_),
213 tx=-self.tx, ty=-self.ty, tz=-self.tz,
214 sx=-self.sx, sy=-self.sy, sz=-self.sz, s=-self.s)
216 def toStr(self, prec=5, name=NN, **unused): # PYCHOK expected
217 '''Return this transform as a string.
219 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
220 @kwarg name: Override name (C{str}) or C{None} to exclude
221 this transform's name.
223 @return: Transform attributes (C{str}).
224 '''
225 return self._instr(name, prec, _tx_, _ty_, _tz_,
226 'rx', 'ry', 'rz', _s_, 's1',
227 _sx_, _sy_, _sz_)
229 def transform(self, x, y, z, inverse=False):
230 '''Transform a (geocentric) Cartesian point, forward or inverse.
232 @arg x: X coordinate (C{meter}).
233 @arg y: Y coordinate (C{meter}).
234 @arg z: Z coordinate (C{meter}).
235 @kwarg inverse: Optional direction, forward or inverse (C{bool}).
237 @return: A L{Vector3Tuple}C{(x, y, z)}, transformed.
238 '''
239 xyz1 = x, y, z, _1_0
240 s1 = self.s1
241 if inverse:
242 xyz1 = map2(neg, xyz1)
243 s1 -= _2_0 # = -(1 - s * 1e-6)) = -(1 - (s1 - 1)) = -(2 - s1)
244 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
245 # x * .rz + y * .s1 - z * .rx + .ty,
246 # -x * .ry + y * .rx + z * .s1 + .tz)
247 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
248 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
249 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
250 name=self.name)
253class Transforms(_NamedEnum):
254 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
255 to accommodate the L{_LazyNamedEnumItem} properties.
256 '''
257 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
258 '''(INTERNAL) Instantiate the C{Transform}.
259 '''
260 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
262Transforms = Transforms(Transform) # PYCHOK singleton
263'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
264# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84
265Transforms._assert(
266 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
267 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
268 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
269 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
270 s=_F( 1.2727)),
271 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
272 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
273 s=_F( -8.3)),
274 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
275 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
276 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
277 s=_F( -9.82)), # Germany
278 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
279 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
280 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
281 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
282 sz=_F( 0.156), s=_F(-1.2)),
283 Identity = _lazy(_Identity_),
284 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
285 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
286 s=_F( -8.15)),
287 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
288 # XXX rotation signs may be opposite, to be checked
289 sx=_F( -1.042), sy=_F( -0.214), sz=_F( -0.631),
290 s=_F( -1.1)),
291 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
292 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
293 s=_F( -2.423)), # spelling
294 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
295 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
296 s=_F( -2.423)), # spelling
297 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
298 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
299 s=_F( -2.423)), # Austria
300 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
301 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
302 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
303 s=_F(-0.0015)),
304 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
305 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
306 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
307 s=_F( 20.4894)),
308 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
309 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
310 WGS84 = _lazy(_WGS84_), # unity
311)
314class Datum(_NamedEnumItem):
315 '''Ellipsoid and transform parameters for an earth model.
316 '''
317 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
318 _transform = Transforms.WGS84 # default transform (L{Transform})
320 def __init__(self, ellipsoid, transform=None, name=NN):
321 '''New L{Datum}.
323 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
324 @kwarg transform: Optional transform (L{Transform}).
325 @kwarg name: Optional, unique name (C{str}).
327 @raise NameError: Datum with that B{C{name}} already exists.
329 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
330 nor L{Ellipsoid2} or B{C{transform}} is
331 not a L{Transform}.
332 '''
333 self._ellipsoid = ellipsoid or Datum._ellipsoid
334 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
336 self._transform = transform or Datum._transform
337 _xinstanceof(Transform, transform=self.transform)
339 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
341 def __eq__(self, other):
342 '''Compare this and an other datum.
344 @arg other: The other datum (L{Datum}).
346 @return: C{True} if equal, C{False} otherwise.
347 '''
348 return self is other or (isinstance(other, Datum) and
349 self.ellipsoid == other.ellipsoid and
350 self.transform == other.transform)
352 def __hash__(self):
353 return self._hash # memoized
355 def __matmul__(self, other): # PYCHOK Python 3.5+
356 '''Convert cartesian or ellipsoidal B{C{other}} to this datum.
358 @raise TypeError: Invalid B{C{other}}.
359 '''
360 try: # only CartesianBase and EllipsoidalLatLonBase
361 return other.toDatum(self)
362 except AttributeError:
363 pass
364 raise _IsnotError(_cartesian_, _ellipsoidal_, other=other)
366 def ecef(self, Ecef=None):
367 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
369 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
371 @return: An ECEF converter for this C{datum}.
373 @raise TypeError: Invalid B{C{Ecef}}.
375 @see: Module L{pygeodesy.ecef}.
376 '''
377 return _MODS.ecef._4Ecef(self, Ecef)
379 @Property_RO
380 def ellipsoid(self):
381 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
382 '''
383 return self._ellipsoid
385 @Property_RO
386 def exactTM(self):
387 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
388 '''
389 return _MODS.etm.ExactTransverseMercator(datum=self)
391 @Property_RO
392 def _hash(self):
393 return hash(self.ellipsoid) + hash(self.transform)
395 @Property_RO
396 def isEllipsoidal(self):
397 '''Check whether this datum is ellipsoidal (C{bool}).
398 '''
399 return self.ellipsoid.isEllipsoidal
401 @Property_RO
402 def isOblate(self):
403 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
404 '''
405 return self.ellipsoid.isOblate
407 @Property_RO
408 def isProlate(self):
409 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
410 '''
411 return self.ellipsoid.isProlate
413 @Property_RO
414 def isSpherical(self):
415 '''Check whether this datum is (near-)spherical (C{bool}).
416 '''
417 return self.ellipsoid.isSpherical
419 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
420 '''Return this datum as a string.
422 @kwarg sep: Separator to join (C{str}).
423 @kwarg name: Override name (C{str}) or C{None} to exclude
424 this datum's name.
426 @return: Datum attributes (C{str}).
427 '''
428 t = [] if name is None else \
429 [Fmt.EQUAL(name=repr(name or self.named))]
430 for a in (_ellipsoid_, _transform_):
431 v = getattr(self, a)
432 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
433 return sep.join(t)
435 @Property_RO
436 def transform(self):
437 '''Get this datum's transform (L{Transform}).
438 '''
439 return self._transform
442def _En2(earth, name):
443 '''(INTERNAL) Helper for C{_ellipsoid} and C{_ellipsoidal_datum}.
444 '''
445 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
446 E = earth
447 n = _under(name or E.name)
448 elif isinstance(earth, Datum):
449 E = earth.ellipsoid
450 n = _under(name or earth.name)
451 elif isinstance(earth, a_f2Tuple):
452 n = _under(name or earth.name)
453 E = Ellipsoid(earth.a, earth.b, name=n)
454 elif islistuple(earth, minum=2):
455 a, f = earth[:2]
456 n = _under(name or _xattr(earth, name=NN))
457 E = Ellipsoid(a, f=f, name=n)
458 else:
459 E, n = None, NN
460 return E, n
463def _a_ellipsoid(a_ellipsoid, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
464 '''(INTERNAL) Get an ellipsoid from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}},
465 an L{Ellipsoid} or L{Ellipsoid2} from L{Datum} or C{a_f2Tuple}.
466 '''
467 if f is None:
468 E, _ = _En2(a_ellipsoid, name)
469 if raiser and not E:
470 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_ellipsoid})
471 else:
472 E = Ellipsoid2(a_ellipsoid, f, name=name)
473 return E
476def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
477 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
478 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
480 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
481 '''
482 if isinstance(earth, Datum):
483 d = earth
484 else:
485 E, n = _En2(earth, name)
486 if not E:
487 n = raiser or _earth_
488 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
489 d = Datum(E, transform=Transforms.Identity, name=n)
490 if raiser and not d.isEllipsoidal:
491 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
492 return d
495def _mean_radius(radius, *lats):
496 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
497 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
498 '''
499 if radius is R_M:
500 r = radius
501 elif isscalar(radius):
502 r = Radius_(radius, low=0, Error=TypeError)
503 else:
504 E = _ellipsoidal_datum(radius).ellipsoid
505 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
506 return r
509def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
510 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
511 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
513 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
514 '''
515 if isscalar(earth):
516 E = Datums.Sphere.ellipsoid
517 if earth == E.a == E.b and not name:
518 d = Datums.Sphere
519 else:
520 r = Radius_(earth, Error=Error) # invalid datum
521 n = _under(name)
522 E = Ellipsoid(r, r, name=n)
523 d = Datum(E, transform=Transforms.Identity, name=n)
524 else:
525 d = _ellipsoidal_datum(earth, Error=Error, name=name)
526 if raiser and not d.isSpherical:
527 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
528 return d
531class Datums(_NamedEnum):
532 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
533 to accommodate the L{_LazyNamedEnumItem} properties.
534 '''
535 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
536 '''(INTERNAL) Instantiate the L{Datum}.
537 '''
538 return Datum(Ellipsoids.get(ellipsoid_name),
539 Transforms.get(transform_name), name=name)
541Datums = Datums(Datum) # PYCHOK singleton
542'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
543# Datums with associated ellipsoid and Helmert transform parameters
544# to convert from WGS84 into the given datum. More are available at
545# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
546# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
547Datums._assert(
548 # Belgian Datum 1972, based on Hayford ellipsoid.
549 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
550 # <https://SpatialReference.org/ref/sr-org/7718/html/>
551 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
553 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
554 # <https://WikiPedia.org/wiki/Helmert_transformation>
555 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
557 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
558 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
560 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
561# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
562# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
563# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
564 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
566 # <https://WikiPedia.org/wiki/GRS_80>
567 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
569 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
570# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
571 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
573 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
574 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
575 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
577 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
578 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
580 # <https://WikiPedia.org/wiki/Helmert_transformation>
581 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
583 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
584 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
585 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
587 # Nouvelle Triangulation Francaise (Paris) XXX verify
588 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
590 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
591 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
593 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
594 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
596 # XXX psuedo-ellipsoids for spherical LatLon
597 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
599 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
600 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
602 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
603 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
605 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
606)
608_WGS84 = Datums.WGS84 # PYCHOK exported internally
610if __name__ == '__main__':
612 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
613 from pygeodesy.lazily import printf
615 # __doc__ of this file, force all into registery
616 for r in (Datums, Transforms):
617 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
618 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
620# **) MIT License
621#
622# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
623#
624# Permission is hereby granted, free of charge, to any person obtaining a
625# copy of this software and associated documentation files (the "Software"),
626# to deal in the Software without restriction, including without limitation
627# the rights to use, copy, modify, merge, publish, distribute, sublicense,
628# and/or sell copies of the Software, and to permit persons to whom the
629# Software is furnished to do so, subject to the following conditions:
630#
631# The above copyright notice and this permission notice shall be included
632# in all copies or substantial portions of the Software.
633#
634# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
635# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
636# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
637# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
638# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
639# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
640# OTHER DEALINGS IN THE SOFTWARE.