Coverage for pygeodesy/datums.py: 97%
182 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -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, neg_, _xinstanceof
68from pygeodesy.constants import _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
70# from pygeodesy.errors import _IsnotError # from .fmath
71from pygeodesy.fmath import fdot, fmean, Fmt, _IsnotError
72from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, _cartesian_, \
73 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _ellipsoid_, \
74 _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \
75 _Krassowsky1940_, _NAD27_, _NAD83_, _name_, _s_, _Sphere_, \
76 _spherical_, _sx_, _sy_, _sz_, _transform_, _tx_, _ty_, _tz_, \
77 _UNDER_, _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.03.29'
91_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)
92_BD72_ = 'BD72'
93_DHDN_ = 'DHDN'
94_earth_ = 'earth'
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) and
177 self.tx == other.tx and
178 self.ty == other.ty and
179 self.tz == other.tz and
180 self.rx == other.rx and
181 self.ry == other.ry and
182 self.rz == other.rz and
183 self.s == other.s)
185 def __matmul__(self, other): # PYCHOK Python 3.5+
186 '''Helmert-transform cartesian B{C{other}}.
188 @raise TypeError: Invalid B{C{other}}.
189 '''
190 try: # only CartesianBase
191 return other.toTransform(self)
192 except AttributeError:
193 pass
194 raise _IsnotError(_cartesian_, other=other)
196 def inverse(self, name=NN):
197 '''Return the inverse of this transform.
199 @kwarg name: Optional, unique name (C{str}).
201 @return: Inverse (Transform).
203 @raise NameError: Transform with that B{C{name}} already exists.
204 '''
205 return Transform(name=name or (self.name + _Inverse_),
206 tx=-self.tx, ty=-self.ty, tz=-self.tz,
207 sx=-self.sx, sy=-self.sy, sz=-self.sz, s=-self.s)
209 def toStr(self, prec=5, name=NN, **unused): # PYCHOK expected
210 '''Return this transform as a string.
212 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
213 @kwarg name: Override name (C{str}) or C{None} to exclude
214 this transform's name.
216 @return: Transform attributes (C{str}).
217 '''
218 return self._instr(name, prec, _tx_, _ty_, _tz_,
219 'rx', 'ry', 'rz', _s_, 's1',
220 _sx_, _sy_, _sz_)
222 def transform(self, x, y, z, inverse=False):
223 '''Transform a (geocentric) Cartesian point, forward or inverse.
225 @arg x: X coordinate (C{meter}).
226 @arg y: Y coordinate (C{meter}).
227 @arg z: Z coordinate (C{meter}).
228 @kwarg inverse: Optional direction, forward or inverse (C{bool}).
230 @return: A L{Vector3Tuple}C{(x, y, z)}, transformed.
231 '''
232 xyz1 = x, y, z, _1_0
233 s1 = self.s1
234 if inverse:
235 xyz1 = neg_(*xyz1)
236 s1 -= _2_0 # == -(1 - s * 1e-6)) == -(1 - (s1 - 1))
237 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,
238 # x * .rz + y * .s1 - z * .rx + .ty,
239 # -x * .ry + y * .rx + z * .s1 + .tz)
240 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),
241 fdot(xyz1, self.rz, s1, -self.rx, self.ty),
242 fdot(xyz1, -self.ry, self.rx, s1, self.tz),
243 name=self.name)
246class Transforms(_NamedEnum):
247 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class
248 to accommodate the L{_LazyNamedEnumItem} properties.
249 '''
250 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):
251 '''(INTERNAL) Instantiate the C{Transform}.
252 '''
253 return Transform(**name_tx_ty_tz_sx_sy_sz_s)
255Transforms = Transforms(Transform) # PYCHOK singleton
256'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''
257# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84
258Transforms._assert(
259 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),
260 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84
261 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>
262 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),
263 s=_F( 1.2727)),
264 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),
265 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),
266 s=_F( -8.3)),
267 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),
268 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),
269 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),
270 s=_F( -9.82)), # Germany
271 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),
272 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)
273 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>
274 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
275 sz=_F( 0.156), s=_F(-1.2)),
276 Identity = _lazy(_Identity_),
277 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
278 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),
279 s=_F( -8.15)),
280 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),
281 # XXX rotation signs may be opposite, to be checked
282 sx=_F( -1.042), sy=_F( -0.214), sz=_F( -0.631),
283 s=_F( -1.1)),
284 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
285 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
286 s=_F( -2.423)), # spelling
287 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),
288 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),
289 s=_F( -2.423)), # spelling
290 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),
291 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),
292 s=_F( -2.423)), # Austria
293 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),
294 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),
295 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),
296 s=_F(-0.0015)),
297 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify
298 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),
299 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),
300 s=_F( 20.4894)),
301 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),
302 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),
303 WGS84 = _lazy(_WGS84_), # unity
304)
307class Datum(_NamedEnumItem):
308 '''Ellipsoid and transform parameters for an earth model.
309 '''
310 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})
311 _transform = Transforms.WGS84 # default transform (L{Transform})
313 def __init__(self, ellipsoid, transform=None, name=NN):
314 '''New L{Datum}.
316 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
317 @kwarg transform: Optional transform (L{Transform}).
318 @kwarg name: Optional, unique name (C{str}).
320 @raise NameError: Datum with that B{C{name}} already exists.
322 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}
323 nor L{Ellipsoid2} or B{C{transform}} is
324 not a L{Transform}.
325 '''
326 self._ellipsoid = ellipsoid or Datum._ellipsoid
327 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)
329 self._transform = transform or Datum._transform
330 _xinstanceof(Transform, transform=self.transform)
332 self._register(Datums, name or self.transform.name or self.ellipsoid.name)
334 def __eq__(self, other):
335 '''Compare this and an other datum.
337 @arg other: The other datum (L{Datum}).
339 @return: C{True} if equal, C{False} otherwise.
340 '''
341 return self is other or (isinstance(other, Datum) and
342 self.ellipsoid == other.ellipsoid and
343 self.transform == other.transform)
345 def __matmul__(self, other): # PYCHOK Python 3.5+
346 '''Convert cartesian or ellipsoidal B{C{other}} to this datum.
348 @raise TypeError: Invalid B{C{other}}.
349 '''
350 try: # only CartesianBase and EllipsoidalLatLonBase
351 return other.toDatum(self)
352 except AttributeError:
353 pass
354 raise _IsnotError(_cartesian_, _ellipsoidal_, other=other)
356 def ecef(self, Ecef=None):
357 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
359 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
361 @return: An ECEF converter for this C{datum}.
363 @raise TypeError: Invalid B{C{Ecef}}.
365 @see: Module L{pygeodesy.ecef}.
366 '''
367 return _MODS.ecef._4Ecef(self, Ecef)
369 @Property_RO
370 def ellipsoid(self):
371 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
372 '''
373 return self._ellipsoid
375 @Property_RO
376 def exactTM(self):
377 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).
378 '''
379 return _MODS.etm.ExactTransverseMercator(datum=self)
381 @Property_RO
382 def isEllipsoidal(self):
383 '''Check whether this datum is ellipsoidal (C{bool}).
384 '''
385 return self._ellipsoid.isEllipsoidal
387 @Property_RO
388 def isOblate(self):
389 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).
390 '''
391 return self._ellipsoid.isOblate
393 @Property_RO
394 def isProlate(self):
395 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).
396 '''
397 return self._ellipsoid.isProlate
399 @Property_RO
400 def isSpherical(self):
401 '''Check whether this datum is (near-)spherical (C{bool}).
402 '''
403 return self._ellipsoid.isSpherical
405 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected
406 '''Return this datum as a string.
408 @kwarg sep: Separator to join (C{str}).
409 @kwarg name: Override name (C{str}) or C{None} to exclude
410 this datum's name.
412 @return: Datum attributes (C{str}).
413 '''
414 t = [] if name is None else \
415 [Fmt.EQUAL(name=repr(name or self.named))]
416 for a in (_ellipsoid_, _transform_):
417 v = getattr(self, a)
418 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))
419 return sep.join(t)
421 @Property_RO
422 def transform(self):
423 '''Get this datum's transform (L{Transform}).
424 '''
425 return self._transform
428def _En2(earth, name):
429 '''(INTERNAL) Helper for C{_ellipsoid} and C{_ellipsoidal_datum}.
430 '''
431 if isinstance(earth, (Ellipsoid, Ellipsoid2)):
432 E = earth
433 n = _UNDER(name or E.name)
434 elif isinstance(earth, Datum):
435 E = earth.ellipsoid
436 n = _UNDER(name or earth.name)
437 elif isinstance(earth, a_f2Tuple):
438 n = _UNDER(name or earth.name)
439 E = Ellipsoid(earth.a, earth.b, name=n)
440 elif islistuple(earth, minum=2):
441 a, f = earth[:2]
442 n = _UNDER(name or getattr(earth, _name_, NN))
443 E = Ellipsoid(a, f=f, name=n)
444 else:
445 E, n = None, NN
446 return E, n
449def _a_ellipsoid(a_ellipsoid, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...
450 '''(INTERNAL) Get an ellipsoid from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}},
451 an L{Ellipsoid} or L{Ellipsoid2} from L{Datum} or C{a_f2Tuple}.
452 '''
453 if f is None:
454 E, _ = _En2(a_ellipsoid, name)
455 if raiser and not E:
456 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_ellipsoid})
457 else:
458 E = Ellipsoid2(a_ellipsoid, f, name=name)
459 return E
462def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):
463 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},
464 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.
466 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.
467 '''
468 if isinstance(earth, Datum):
469 d = earth
470 else:
471 E, n = _En2(earth, name)
472 if not E:
473 n = raiser or _earth_
474 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})
475 d = Datum(E, transform=Transforms.Identity, name=n)
476 if raiser and not d.isEllipsoidal:
477 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})
478 return d
481def _mean_radius(radius, *lats):
482 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},
483 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.
484 '''
485 if isscalar(radius):
486 r = Radius_(radius, low=0, Error=TypeError)
487 else:
488 E = _ellipsoidal_datum(radius).ellipsoid
489 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean
490 return r
493def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):
494 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},
495 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.
497 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.
498 '''
499 if isscalar(earth):
500 E = Datums.Sphere.ellipsoid
501 if earth == E.a == E.b and not name:
502 d = Datums.Sphere
503 else:
504 r = Radius_(earth, Error=Error) # invalid datum
505 n = _UNDER(name)
506 E = Ellipsoid(r, r, name=n)
507 d = Datum(E, transform=Transforms.Identity, name=n)
508 else:
509 d = _ellipsoidal_datum(earth, Error=Error, name=name)
510 if raiser and not d.isSpherical:
511 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})
512 return d
515class Datums(_NamedEnum):
516 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class
517 to accommodate the L{_LazyNamedEnumItem} properties.
518 '''
519 def _Lazy(self, ellipsoid_name, transform_name, name=NN):
520 '''(INTERNAL) Instantiate the L{Datum}.
521 '''
522 return Datum(Ellipsoids.get(ellipsoid_name),
523 Transforms.get(transform_name), name=name)
525Datums = Datums(Datum) # PYCHOK singleton
526'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''
527# Datums with associated ellipsoid and Helmert transform parameters
528# to convert from WGS84 into the given datum. More are available at
529# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and
530# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.
531Datums._assert(
532 # Belgian Datum 1972, based on Hayford ellipsoid.
533 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>
534 # <https://SpatialReference.org/ref/sr-org/7718/html/>
535 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),
537 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>
538 # <https://WikiPedia.org/wiki/Helmert_transformation>
539 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),
541 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>
542 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),
544 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>
545# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?
546# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?
547# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),
548 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?
550 # <https://WikiPedia.org/wiki/GRS_80>
551 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),
553 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2
554# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),
555 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),
557 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
558 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?
559 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?
561 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>
562 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),
564 # <https://WikiPedia.org/wiki/Helmert_transformation>
565 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),
567 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>
568 # (If you *really* must convert WGS84<->NAD83, you need more than this!)
569 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),
571 # Nouvelle Triangulation Francaise (Paris) XXX verify
572 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),
574 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>
575 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),
577 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>
578 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),
580 # XXX psuedo-ellipsoids for spherical LatLon
581 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),
583 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>
584 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),
586 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>
587 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),
589 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),
590)
592_WGS84 = Datums.WGS84 # PYCHOK exported internally
594if __name__ == '__main__':
596 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
598 # __doc__ of this file, force all into registery
599 for r in (Datums, Transforms):
600 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)
601 print(_NLATvar_.join(i.strip(_COMMA_) for i in t))
603# **) MIT License
604#
605# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
606#
607# Permission is hereby granted, free of charge, to any person obtaining a
608# copy of this software and associated documentation files (the "Software"),
609# to deal in the Software without restriction, including without limitation
610# the rights to use, copy, modify, merge, publish, distribute, sublicense,
611# and/or sell copies of the Software, and to permit persons to whom the
612# Software is furnished to do so, subject to the following conditions:
613#
614# The above copyright notice and this permission notice shall be included
615# in all copies or substantial portions of the Software.
616#
617# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
618# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
619# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
620# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
621# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
622# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
623# OTHER DEALINGS IN THE SOFTWARE.