Coverage for pygeodesy/datums.py: 95%
195 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-15 09:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-15 09:43 -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, \
70 _EWGS84, Vector3Tuple
71from pygeodesy.errors import _IsnotError, _xattr
72from pygeodesy.fmath import fdot, fmean, Fmt
73from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, _cartesian_, \
74 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \
75 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \
76 _Krassowsky1940_, _NAD27_, _NAD83_, _s_, _Sphere_, _spherical_, \
77 _sx_, _sy_, _sz_, _transform_, _tx_, _ty_, _tz_, _UNDER_, \
78 _WGS72_, _WGS84_, _under
79from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
80from pygeodesy.named import _NamedEnum, _NamedEnumItem, \
81 _lazyNamedEnumItem as _lazy, Property_RO
82# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids
83# from pygeodesy.props import Property_RO # from .named
84# from pygeodesy.streprs import Fmt # from .fmath
85from pygeodesy.units import radians, Radius_
87# from math import radians # from .units
89__all__ = _ALL_LAZY.datums
90__version__ = '23.08.20'
92_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
93_BD72_ = 'BD72'
94_DHDN_ = 'DHDN'
95_ED50_ = 'ED50'
96_GDA2020_ = 'GDA2020'
97_Identity_ = 'Identity'
98_Inverse_ = 'Inverse'
99_Irl1965_ = 'Irl1965'
100_Irl1975_ = 'Irl1975'
101_MGI_ = 'MGI'
102_NTF_ = 'NTF'
103_OSGB36_ = 'OSGB36'
104_Potsdam_ = 'Potsdam'
105_TokyoJapan_ = 'TokyoJapan'
107_r_s1 = radians(1 / _3600_0) # 1 degree second to radians
110def _r_s2(s_):
111 '''(INTERNAL) rotation in C{radians} and C{degree seconds}.
112 '''
113 return _F(_r_s1 * s_), s_
116class Transform(_NamedEnumItem):
117 '''Helmert transformation.
119 @see: L{Helmert7Tuple}.
120 '''
121 tx = _0_0 # x translation (C{meter})
122 ty = _0_0 # y translation (C{meter})
123 tz = _0_0 # z translation (C{meter})
125 rx = _0_0 # x rotation (C{radians})
126 ry = _0_0 # y rotation (C{radians})
127 rz = _0_0 # z rotation (C{radians})
129 s = _0_0 # scale ppm (C{float})
130 s1 = _1_0 # scale + 1 (C{float})
132 sx = _0_0 # x rotation (degree seconds)
133 sy = _0_0 # y rotation (degree seconds)
134 sz = _0_0 # z rotation (degree seconds)
136 def __init__(self, name=NN, tx=0, ty=0, tz=0,
137 sx=0, sy=0, sz=0, s=0):
138 '''New L{Transform}.
140 @kwarg name: Optional, unique name (C{str}).
141 @kwarg tx: Optional X translation (C{meter}).
142 @kwarg ty: Optional Y translation (C{meter}).
143 @kwarg tz: Optional Z translation (C{meter}).
144 @kwarg s: Optional scale (C{float}), ppm.
145 @kwarg sx: Optional X rotation (C{degree seconds}).
146 @kwarg sy: Optional Y rotation (C{degree seconds}).
147 @kwarg sz: Optional Z rotation (C{degree seconds}).
149 @raise NameError: Transform with that B{C{name}} already exists.
150 '''
151 if tx:
152 self.tx = tx
153 if ty:
154 self.ty = ty
155 if tz:
156 self.tz = tz
157 if sx: # secs to rads
158 self.rx, self.sx = _r_s2(sx)
159 if sy:
160 self.ry, self.sy = _r_s2(sy)
161 if sz:
162 self.rz, self.sz = _r_s2(sz)
163 if s:
164 self.s = s
165 self.s1 = _F(s * 1e-6 + _1_0) # normalize ppm to (s + 1)
167 self._register(Transforms, name)
169 def __eq__(self, other):
170 '''Compare this and an other transform.
172 @arg other: The other transform (L{Transform}).
174 @return: C{True} if equal, C{False} otherwise.
175 '''
176 return self is other or (isinstance(other, Transform)
177 and self.tx == other.tx
178 and self.ty == other.ty
179 and self.tz == other.tz
180 and self.rx == other.rx
181 and self.ry == other.ry
182 and self.rz == other.rz
183 and self.s == other.s)
185 def __hash__(self):
186 return self._hash # memoized
188 def __matmul__(self, other): # PYCHOK Python 3.5+
189 '''Helmert-transform a cartesian B{C{other}}.
191 @raise TypeError: Invalid B{C{other}}.
192 '''
193 try: # only CartesianBase
194 return other.toTransform(self)
195 except AttributeError:
196 pass
197 raise _IsnotError(_cartesian_, other=other)
199 @Property_RO
200 def _hash(self):
201 return hash((self.rx, self.ry, self.rz, self.s,
202 self.tx, self.ty, self.tz))
204 def inverse(self, name=NN):
205 '''Return the inverse of this transform.
207 @kwarg name: Optional, unique name (C{str}).
209 @return: Inverse (Transform).
211 @raise NameError: Transform with that B{C{name}} already exists.
212 '''
213 return Transform(name=name or (self.name + _Inverse_),
214 tx=-self.tx, ty=-self.ty, tz=-self.tz,
215 sx=-self.sx, sy=-self.sy, sz=-self.sz, s=-self.s)
217 def toStr(self, prec=5, name=NN, **unused): # PYCHOK expected
218 '''Return this transform as a string.
220 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
221 @kwarg name: Override name (C{str}) or C{None} to exclude
222 this transform's name.
224 @return: Transform attributes (C{str}).
225 '''
226 return self._instr(name, prec, _tx_, _ty_, _tz_,
227 'rx', 'ry', 'rz', _s_, 's1',
228 _sx_, _sy_, _sz_)
230 def transform(self, x, y, z, inverse=False):
231 '''Transform a (geocentric) Cartesian point, forward or inverse.
233 @arg x: X coordinate (C{meter}).
234 @arg y: Y coordinate (C{meter}).
235 @arg z: Z coordinate (C{meter}).
236 @kwarg inverse: Optional direction, forward or inverse (C{bool}).
238 @return: A L{Vector3Tuple}C{(x, y, z)}, transformed.
239 '''
240 xyz1 = x, y, z, _1_0
241 s1 = self.s1
242 if inverse:
243 xyz1 = map2(neg, xyz1)
244 s1 -= _2_0 # = -(1 - s * 1e-6)) = -(1 - (s1 - 1)) = -(2 - s1)
245 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
246 # x * .rz + y * .s1 - z * .rx + .ty,
247 # -x * .ry + y * .rx + z * .s1 + .tz)
248 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
249 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
250 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
251 name=self.name)
254class Transforms(_NamedEnum):
255 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
256 to accommodate the L{_LazyNamedEnumItem} properties.
257 '''
258 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
259 '''(INTERNAL) Instantiate the C{Transform}.
260 '''
261 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
263Transforms = Transforms(Transform) # PYCHOK singleton
264'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
265# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84
266Transforms._assert(
267 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
268 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
269 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
270 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
271 s=_F( 1.2727)),
272 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
273 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
274 s=_F( -8.3)),
275 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
276 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
277 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
278 s=_F( -9.82)), # Germany
279 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
280 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
281 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
282 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
283 sz=_F( 0.156), s=_F(-1.2)),
284 Identity = _lazy(_Identity_),
285 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
286 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
287 s=_F( -8.15)),
288 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
289 # XXX rotation signs may be opposite, to be checked
290 sx=_F( -1.042), sy=_F( -0.214), sz=_F( -0.631),
291 s=_F( -1.1)),
292 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
293 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
294 s=_F( -2.423)), # spelling
295 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
296 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
297 s=_F( -2.423)), # spelling
298 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
299 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
300 s=_F( -2.423)), # Austria
301 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
302 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
303 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
304 s=_F(-0.0015)),
305 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
306 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
307 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
308 s=_F( 20.4894)),
309 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
310 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
311 WGS84 = _lazy(_WGS84_), # unity
312)
315class Datum(_NamedEnumItem):
316 '''Ellipsoid and transform parameters for an earth model.
317 '''
318 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
319 _transform = Transforms.WGS84 # default transform (L{Transform})
321 def __init__(self, ellipsoid, transform=None, name=NN):
322 '''New L{Datum}.
324 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
325 @kwarg transform: Optional transform (L{Transform}).
326 @kwarg name: Optional, unique name (C{str}).
328 @raise NameError: Datum with that B{C{name}} already exists.
330 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
331 nor L{Ellipsoid2} or B{C{transform}} is
332 not a L{Transform}.
333 '''
334 self._ellipsoid = ellipsoid or Datum._ellipsoid
335 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
337 self._transform = transform or Datum._transform
338 _xinstanceof(Transform, transform=self.transform)
340 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
342 def __eq__(self, other):
343 '''Compare this and an other datum.
345 @arg other: The other datum (L{Datum}).
347 @return: C{True} if equal, C{False} otherwise.
348 '''
349 return self is other or (isinstance(other, Datum) and
350 self.ellipsoid == other.ellipsoid and
351 self.transform == other.transform)
353 def __hash__(self):
354 return self._hash # memoized
356 def __matmul__(self, other): # PYCHOK Python 3.5+
357 '''Convert cartesian or ellipsoidal B{C{other}} to this datum.
359 @raise TypeError: Invalid B{C{other}}.
360 '''
361 try: # only CartesianBase and EllipsoidalLatLonBase
362 return other.toDatum(self)
363 except AttributeError:
364 pass
365 raise _IsnotError(_cartesian_, _ellipsoidal_, other=other)
367 def ecef(self, Ecef=None):
368 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
370 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
372 @return: An ECEF converter for this C{datum}.
374 @raise TypeError: Invalid B{C{Ecef}}.
376 @see: Module L{pygeodesy.ecef}.
377 '''
378 return _MODS.ecef._4Ecef(self, Ecef)
380 @Property_RO
381 def ellipsoid(self):
382 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
383 '''
384 return self._ellipsoid
386 @Property_RO
387 def exactTM(self):
388 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
389 '''
390 return _MODS.etm.ExactTransverseMercator(datum=self)
392 @Property_RO
393 def _hash(self):
394 return hash(self.ellipsoid) + hash(self.transform)
396 @Property_RO
397 def isEllipsoidal(self):
398 '''Check whether this datum is ellipsoidal (C{bool}).
399 '''
400 return self.ellipsoid.isEllipsoidal
402 @Property_RO
403 def isOblate(self):
404 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
405 '''
406 return self.ellipsoid.isOblate
408 @Property_RO
409 def isProlate(self):
410 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
411 '''
412 return self.ellipsoid.isProlate
414 @Property_RO
415 def isSpherical(self):
416 '''Check whether this datum is (near-)spherical (C{bool}).
417 '''
418 return self.ellipsoid.isSpherical
420 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
421 '''Return this datum as a string.
423 @kwarg sep: Separator to join (C{str}).
424 @kwarg name: Override name (C{str}) or C{None} to exclude
425 this datum's name.
427 @return: Datum attributes (C{str}).
428 '''
429 t = [] if name is None else \
430 [Fmt.EQUAL(name=repr(name or self.named))]
431 for a in (_ellipsoid_, _transform_):
432 v = getattr(self, a)
433 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
434 return sep.join(t)
436 @Property_RO
437 def transform(self):
438 '''Get this datum's transform (L{Transform}).
439 '''
440 return self._transform
443def _En2(earth, name):
444 '''(INTERNAL) Helper for C{_ellipsoid} and C{_ellipsoidal_datum}.
445 '''
446 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
447 E = earth
448 n = _under(name or E.name)
449 elif isinstance(earth, Datum):
450 E = earth.ellipsoid
451 n = _under(name or earth.name)
452 elif isinstance(earth, a_f2Tuple):
453 n = _under(name or earth.name)
454 E = Ellipsoid(earth.a, earth.b, name=n)
455 elif islistuple(earth, minum=2):
456 a, f = earth[:2]
457 n = _under(name or _xattr(earth, name=NN))
458 E = Ellipsoid(a, f=f, name=n)
459 else:
460 E, n = None, NN
461 return E, n
464def _a_ellipsoid(a_ellipsoid, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
465 '''(INTERNAL) Get an ellipsoid from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}},
466 an L{Ellipsoid} or L{Ellipsoid2} from L{Datum} or C{a_f2Tuple}.
467 '''
468 if f is None:
469 E, _ = _En2(a_ellipsoid, name)
470 if raiser and not E:
471 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_ellipsoid})
472 else:
473 E = Ellipsoid2(a_ellipsoid, f, name=name)
474 return E
477def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
478 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
479 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
481 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
482 '''
483 if isinstance(earth, Datum):
484 d = earth
485 else:
486 E, n = _En2(earth, name)
487 if not E:
488 n = raiser or _earth_
489 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
490 d = Datum(E, transform=Transforms.Identity, name=n)
491 if raiser and not d.isEllipsoidal:
492 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
493 return d
496def _mean_radius(radius, *lats):
497 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
498 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
499 '''
500 if radius is R_M:
501 r = radius
502 elif isscalar(radius):
503 r = Radius_(radius, low=0, Error=TypeError)
504 else:
505 E = _ellipsoidal_datum(radius).ellipsoid
506 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
507 return r
510def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
511 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
512 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
514 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
515 '''
516 if isscalar(earth):
517 E = Datums.Sphere.ellipsoid
518 if earth == E.a == E.b and not name:
519 d = Datums.Sphere
520 else:
521 r = Radius_(earth, Error=Error) # invalid datum
522 n = _under(name)
523 E = Ellipsoid(r, r, name=n)
524 d = Datum(E, transform=Transforms.Identity, name=n)
525 else:
526 d = _ellipsoidal_datum(earth, Error=Error, name=name)
527 if raiser and not d.isSpherical:
528 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
529 return d
532class Datums(_NamedEnum):
533 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
534 to accommodate the L{_LazyNamedEnumItem} properties.
535 '''
536 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
537 '''(INTERNAL) Instantiate the L{Datum}.
538 '''
539 return Datum(Ellipsoids.get(ellipsoid_name),
540 Transforms.get(transform_name), name=name)
542Datums = Datums(Datum) # PYCHOK singleton
543'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
544# Datums with associated ellipsoid and Helmert transform parameters
545# to convert from WGS84 into the given datum. More are available at
546# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
547# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
548Datums._assert(
549 # Belgian Datum 1972, based on Hayford ellipsoid.
550 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
551 # <https://SpatialReference.org/ref/sr-org/7718/html/>
552 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
554 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
555 # <https://WikiPedia.org/wiki/Helmert_transformation>
556 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
558 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
559 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
561 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
562# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
563# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
564# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
565 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
567 # <https://WikiPedia.org/wiki/GRS_80>
568 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
570 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
571# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
572 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
574 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
575 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
576 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
578 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
579 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
581 # <https://WikiPedia.org/wiki/Helmert_transformation>
582 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
584 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
585 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
586 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
588 # Nouvelle Triangulation Francaise (Paris) XXX verify
589 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
591 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
592 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
594 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
595 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
597 # XXX psuedo-ellipsoids for spherical LatLon
598 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
600 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
601 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
603 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
604 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
606 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
607)
609_WGS84 = Datums.WGS84
610assert _WGS84.ellipsoid is _EWGS84
612if __name__ == '__main__':
614 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
615 from pygeodesy.lazily import printf
617 # __doc__ of this file, force all into registery
618 for r in (Datums, Transforms):
619 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
620 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
622# **) MIT License
623#
624# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
625#
626# Permission is hereby granted, free of charge, to any person obtaining a
627# copy of this software and associated documentation files (the "Software"),
628# to deal in the Software without restriction, including without limitation
629# the rights to use, copy, modify, merge, publish, distribute, sublicense,
630# and/or sell copies of the Software, and to permit persons to whom the
631# Software is furnished to do so, subject to the following conditions:
632#
633# The above copyright notice and this permission notice shall be included
634# in all copies or substantial portions of the Software.
635#
636# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
637# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
638# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
639# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
640# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
641# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
642# OTHER DEALINGS IN THE SOFTWARE.