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