1

2# -*- coding: utf-8 -*-

3

4 u'''Datums and transformations thereof.

5

6 Classes L{Datum} and L{Transform} and registries L{Datums} and L{Transforms}, respectively.

7

8 Pure Python implementation of geodesy tools for ellipsoidal earth models, including datums

9 and ellipsoid parameters for different geographic coordinate systems and methods for

10 converting between them and to cartesian coordinates. Transcoded from JavaScript originals by

11 I{(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>}.

13

14 Historical geodetic datums: a latitude/longitude point defines a geographic location on, above

15 or below the earth’s surface. Latitude is measured in degrees from the equator, lomgitude from

16 the International Reference Meridian and height in meters above an ellipsoid based on the given

17 datum. The datum in turn is based on a reference ellipsoid and tied to geodetic survey

18 reference points.

19

20 Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems), but

21 previously various other reference ellipsoids and datum references were used.

22

23 The UK Ordnance Survey National Grid References are still based on the otherwise historical OSGB36

24 datum, 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>}.

26

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)

45

46@var Transforms.BD72: Transform(name='BD72', tx=106.87, ty=-52.298, tz=103.72, s1=1.0, rx=-1.6317e-06, ry=-2.2154e-06, rz=-8.9311e-06, s=1.2727, sx=-0.33657, sy=-0.45696, sz=-1.8422)

47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, s1=0.99999, rx=-5.0421e-06, ry=-1.6968e-06, rz=1.4932e-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.5682e-07, rz=-7.0686e-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=1.0, rx=2.618e-07, ry=-2.7634e-07, rz=1.356e-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.99999, rx=-9.7932e-07, ry=-2.1817e-07, rz=1.1902e-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=1.0, rx=0.0, ry=0.0, rz=7.5631e-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.6, tz=-564.56, s1=0.99999, rx=5.0518e-06, ry=1.0375e-06, rz=3.0592e-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.6, tz=-564.56, s1=0.99999, rx=5.0518e-06, ry=1.0375e-06, rz=3.0592e-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=1.0, rx=-9.6963e-08, ry=1.2605e-06, rz=6.3026e-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=1.0, rx=-9.6963e-08, ry=1.2605e-06, rz=6.3026e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)

58@var Transforms.MGI: Transform(name='MGI', tx=-577.33, ty=-90.129, tz=-463.92, s1=1.0, rx=2.4905e-05, ry=7.1462e-06, rz=2.5681e-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.2945e-07, ry=1.6484e-09, rz=5.333e-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.45, ty=125.16, tz=-542.06, s1=1.0, rx=-7.2819e-07, ry=-1.1975e-06, rz=-4.0826e-06, s=20.489, 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.6859e-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

68 from __future__ import division as _; del _ # PYCHOK semicolon

69

70 from pygeodesy.basics import islistuple, map2, neg, _xinstanceof

71 from 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

74 from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, _EWGS84, \

75 Vector3Tuple

76 from pygeodesy.errors import _IsnotError, _TypeError, _xattr, _xellipsoidall, _xkwds, \

77 _xkwds_pop2

78 from pygeodesy.fmath import fdot, fmean, Fmt

79 from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _Bessel1841_, \

80 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \

81 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \

82 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \

83 _NAD83_, _s_, _Sphere_, _spherical_, _transform_, \

84 _UNDER_, _WGS72_, _WGS84_, _under

85 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

86 from pygeodesy.named import _NamedEnum, _NamedEnumItem, _lazyNamedEnumItem as _lazy

87# from pygeodesy.namedTuples import Vector3Tuple # from .ellipsoids

88 from pygeodesy.props import Property_RO, property_RO

89# from pygeodesy.streprs import Fmt # from .fmath

90 from pygeodesy.units import _isRadius, Radius_, radians

91

92# from math import radians # from .units

93

94__all__ = _ALL_LAZY.datums

95__version__ = '24.02.24'

96

97_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

98_BD72_ = 'BD72'

99_DHDN_ = 'DHDN'

100_DHDNE_ = 'DHDNE'

101_DHDNW_ = 'DHDNW'

102_ED50_ = 'ED50'

103_GDA2020_ = 'GDA2020' # in .trf

104_Identity_ = 'Identity'

105_Irl1965_ = 'Irl1965'

106_Irl1975_ = 'Irl1975'

107_MGI_ = 'MGI'

108_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf

109_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:]

110_NTF_ = 'NTF'

111_OSGB36_ = 'OSGB36'

112_Potsdam_ = 'Potsdam'

113_RPS = radians(_1_0 / _3600_0) # radians per degree-second

114_S1_S = 1.e-6 # in .trf

115_TokyoJapan_ = 'TokyoJapan'

116

117

118 class Transform(_NamedEnumItem):

119 '''Helmert I{datum} transformation.

120

121 @see: L{TransformXform<trf.TransformXform>}.

122 '''

123 tx = _0_0 # x translation (C{meter})

124 ty = _0_0 # y translation (C{meter})

125 tz = _0_0 # z translation (C{meter})

126

127 rx = _0_0 # x rotation (C{radians})

128 ry = _0_0 # y rotation (C{radians})

129 rz = _0_0 # z rotation (C{radians})

130

131 s = _0_0 # scale ppm (C{float})

132 s1 = _1_0 # scale + 1 (C{float})

133

134 sx = _0_0 # x rotation (C{degree-seconds})

135 sy = _0_0 # y rotation (C{degree-seconds})

136 sz = _0_0 # z rotation (C{degree-seconds})

137

138 def __init__(self, name=NN, tx=0, ty=0, tz=0,

139 sx=0, sy=0, sz=0, s=0):

140 '''New L{Transform}.

141

142 @kwarg name: Optional, unique name (C{str}).

143 @kwarg tx: X translation (C{meter}).

144 @kwarg ty: Y translation (C{meter}).

145 @kwarg tz: Z translation (C{meter}).

146 @kwarg s: Scale (C{float}), ppm.

147 @kwarg sx: X rotation (C{degree-seconds}).

148 @kwarg sy: Y rotation (C{degree-seconds}).

149 @kwarg sz: Z rotation (C{degree-seconds}).

150

151 @raise NameError: Transform with that B{C{name}} already exists.

152 '''

153 if tx:

154 self.tx = tx

155 if ty:

156 self.ty = ty

157 if tz:

158 self.tz = tz

159 if sx: # secs to rads

160 self.rx, self.sx = self._rps2(sx)

161 if sy:

162 self.ry, self.sy = self._rps2(sy)

163 if sz:

164 self.rz, self.sz = self._rps2(sz)

165 if s:

166 self.s = s

167 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1)

168

169 self._register(Transforms, name)

170

171 def __eq__(self, other):

172 '''Compare this and an other transform.

173

174 @arg other: The other transform (L{Transform}).

175

176 @return: C{True} if equal, C{False} otherwise.

177 '''

178 return self is other or (isinstance(other, Transform) and all(

179 a == b for a, b in zip(self, other)))

180

181 def __hash__(self):

182 return hash(x for x in self)

183

184 def __iter__(self):

185 '''Yield the initial attribute values.

186 '''

187 for _, x in self.items():

188 yield x

189

190 def __matmul__(self, point): # PYCHOK Python 3.5+

191 '''Transform an I{ellipsoidal} B{C{point}} with this Helmert.

192

193 @return: A transformed copy of B{C{point}}.

194

195 @raise TypeError: Invalid B{C{point}}.

196

197 @see: Method C{B{point}.toTransform}.

198 '''

199 _ = _xellipsoidall(point)

200 return point.toTransform(self)

201

202 def __neg__(self):

203 return self.inverse()

204

205# def __sub__(self, other):

206# _xinstanceof(Transform, other=other)

207#

208# def _sub(a, b):

209# for n in _Names11:

210# yield getattr(a, n) - getattr(b, n)

211#

212# return type(self)(_sub(self, other), name=_MINUS_) # .fsums._sub_op

213

214# @Property_RO

215# def _hash(self):

216# return hash(x for x in self)

217

218 def items(self):

219 '''Yield each attribute as 2-tuple C{(name, value)}.

220 '''

221 for n in _Names7:

222 yield n, getattr(self, n)

223

224 def inverse(self, name=NN):

225 '''Return the inverse of this transform.

226

227 @kwarg name: Optional, unique name (C{str}).

228

229 @return: Inverse (L{Transform}).

230

231 @raise NameError: Transform with that B{C{name}} already exists.

232 '''

233 d = dict((n, neg(v)) for n, v in self.items())

234 n = name or _minus(self.name)

235 return type(self)(name=n, **d)

236

237 @Property_RO

238 def isunity(self):

239 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84?

240 '''

241 return not any(s for s in self)

242

243 def _rps2(self, s_):

244 '''(INTERNAL) Rotation in C{radians} and C{degree-seconds}.

245 '''

246 return (_RPS * s_), s_

247

248 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected

249 '''Return this transform as a string.

250

251 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

252 @kwarg fmt: Optional C{float} format (C{letter}).

253 @kwarg name: Override name (C{str}) or C{None} to exclude

254 this transform's name.

255

256 @return: Transform attributes (C{str}).

257 '''

258 return self._instr(name, prec, *_Names11, fmt=fmt)

259

260 def transform(self, x, y, z, inverse=False, **Vector_and_kwds):

261 '''Transform a (geocentric) position, forward or inverse.

262

263 @arg x: X coordinate (C{meter}).

264 @arg y: Y coordinate (C{meter}).

265 @arg z: Z coordinate (C{meter}).

266 @kwarg inverse: If C{True}, apply the inverse transform (C{bool}).

267 @kwarg Vector_and_kwds: An optional, (3-D) C{B{Vector}=None} or

268 cartesian class and additional C{B{Vector}}

269 keyword arguments to return the transformed

270 point.

271

272 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)})

273 unless some B{C{Vector_and_kwds}} is specified.

274 '''

275 if self.isunity:

276 r = Vector3Tuple(x, y, z, name=self.name)

277 else:

278 xyz1 = x, y, z, _1_0

279 s1 = self.s1

280 if inverse:

281 xyz1 = map2(neg, xyz1)

282 s1 -= _2_0 # = s * 1e-6 - 1 = (s1 - 1) - 1

283 # x', y', z' = (x * .s1 - y * .rz + z * .ry + .tx,

284 # x * .rz + y * .s1 - z * .rx + .ty,

285 # -x * .ry + y * .rx + z * .s1 + .tz)

286 r = Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),

287 fdot(xyz1, self.rz, s1, -self.rx, self.ty),

288 fdot(xyz1, -self.ry, self.rx, s1, self.tz),

289 name=self.name)

290 if Vector_and_kwds:

291 V, kwds = _xkwds_pop2(Vector_and_kwds, Vector=None)

292 if V:

293 r = V(r, **_xkwds(kwds, name=self.name))

294 return r

295

296

297 class Transforms(_NamedEnum):

298 '''(INTERNAL) L{Transform} registry, I{must} be a sub-class

299 to accommodate the L{_LazyNamedEnumItem} properties.

300 '''

301 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

302 '''(INTERNAL) Instantiate the C{Transform}.

303 '''

304 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

305

306 Transforms = Transforms(Transform) # PYCHOK singleton

307'''Some pre-defined L{Transform}s, all I{lazily} instantiated.'''

308# <https://WikiPedia.org/wiki/Helmert_transformation> from WGS84 to ...

309 Transforms._assert(

310 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893),

311 # <https://www.NGI.Be/FR/FR4-4.shtm> ETRS89 == WG84

312 # <https://EPSG.org/transformation_15929/BD72-to-WGS-84-3.html>

313 sx=_F(-0.33657), sy=_F( -0.456955), sz=_F( -1.84218),

314 s=_F( 1.2727)),

315 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0),

316 sx=_F( -1.04), sy=_F( -0.35), sz=_F( 3.08),

317 s=_F( -8.3)),

318 Clarke1866 = _lazy(_Clarke1866_, tx=_F(8), ty=_F(-160), tz=_F(-176)),

319 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39),

320 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458),

321 s=_F( -9.82)), # Germany

322 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2),

323 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>

324 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797),

325 s=_F( -2.55)), # East Germany

326 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2),

327 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>

328 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455),

329 s=_F( -6.7)), # West Germany

330 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1),

331 # <https://GeoNet.ESRI.com/thread/36583> sz=_F(-0.156)

332 # <https://GitHub.com/ChrisVeness/geodesy/blob/master/latlon-ellipsoidal.js>

333 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>

334 sz=_F( 0.156), s=_F(-1.2)),

335 Identity = _lazy(_Identity_),

336 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),

337 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>

338 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),

339 s=_F( -8.15)),

340 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557),

341 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>

342 sx=_F( 1.042), sy=_F( 0.214), sz=_F( 0.631),

343 s=_F( -8.15)),

344 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),

345 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),

346 s=_F( -2.423)), # spelling

347 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0),

348 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13),

349 s=_F( -2.423)), # spelling

350 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920),

351 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297),

352 s=_F( -2.423)), # Austria

353 NAD27 = _lazy(_NAD27_, tx=_8_0, ty=_F(-160), tz=_F(-176)),

354 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515),

355 sx=_F( 0.0267), sy=_F( 0.00034), sz=_F( 0.011),

356 s=_F(-0.0015)),

357 NTF = _lazy(_NTF_, tx=_F(-168), ty=_F(-60), tz=_F(320)), # XXX verify

358 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060),

359 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>

360 sx=_F( -0.1502), sy=_F( -0.2470), sz=_F( -0.8421),

361 s=_F( 20.4894)),

362 TokyoJapan = _lazy(_TokyoJapan_, tx=_F(148), ty=_F(-507), tz=_F(-685)),

363 WGS72 = _lazy(_WGS72_, tz=_F(-4.5), sz=_F(0.554), s=_F(-0.22)),

364 WGS84 = _lazy(_WGS84_), # unity

365)

366

367

368 class Datum(_NamedEnumItem):

369 '''Ellipsoid and transform parameters for an earth model.

370 '''

371 _ellipsoid = Ellipsoids.WGS84 # default ellipsoid (L{Ellipsoid}, L{Ellipsoid2})

372 _transform = Transforms.WGS84 # default transform (L{Transform})

373

374 def __init__(self, ellipsoid, transform=None, name=NN):

375 '''New L{Datum}.

376

377 @arg ellipsoid: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).

378 @kwarg transform: Optional transform (L{Transform}).

379 @kwarg name: Optional, unique name (C{str}).

380

381 @raise NameError: Datum with that B{C{name}} already exists.

382

383 @raise TypeError: If B{C{ellipsoid}} is not an L{Ellipsoid}

384 nor L{Ellipsoid2} or B{C{transform}} is

385 not a L{Transform}.

386 '''

387 self._ellipsoid = ellipsoid or Datum._ellipsoid

388 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

389

390 self._transform = transform or Datum._transform

391 _xinstanceof(Transform, transform=self.transform)

392

393 self._register(Datums, name or self.transform.name or self.ellipsoid.name)

394

395 def __eq__(self, other):

396 '''Compare this and an other datum.

397

398 @arg other: The other datum (L{Datum}).

399

400 @return: C{True} if equal, C{False} otherwise.

401 '''

402 return self is other or (isinstance(other, Datum) and

403 self.ellipsoid == other.ellipsoid and

404 self.transform == other.transform)

405

406 def __hash__(self):

407 return self._hash # memoized

408

409 def __matmul__(self, point): # PYCHOK Python 3.5+

410 '''Convert an I{ellipsoidal} B{C{point}} to this datum.

411

412 @raise TypeError: Invalid B{C{point}}.

413 '''

414 _ = _xellipsoidall(point)

415 return point.toDatum(self)

416

417 def ecef(self, Ecef=None):

418 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.

419

420 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.

421

422 @return: An ECEF converter for this C{datum}.

423

424 @raise TypeError: Invalid B{C{Ecef}}.

425

426 @see: Module L{pygeodesy.ecef}.

427 '''

428 return _MODS.ecef._4Ecef(self, Ecef)

429

430 @Property_RO

431 def ellipsoid(self):

432 '''Get this datum's ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).

433 '''

434 return self._ellipsoid

435

436 @Property_RO

437 def exactTM(self):

438 '''Get the C{ExactTM} projection (L{ExactTransverseMercator}).

439 '''

440 return _MODS.etm.ExactTransverseMercator(datum=self)

441

442 @Property_RO

443 def _hash(self):

444 return hash(self.ellipsoid) + hash(self.transform)

445

446 @property_RO

447 def isEllipsoidal(self):

448 '''Check whether this datum is ellipsoidal (C{bool}).

449 '''

450 return self.ellipsoid.isEllipsoidal

451

452 @property_RO

453 def isOblate(self):

454 '''Check whether this datum's ellipsoidal is I{oblate} (C{bool}).

455 '''

456 return self.ellipsoid.isOblate

457

458 @property_RO

459 def isProlate(self):

460 '''Check whether this datum's ellipsoidal is I{prolate} (C{bool}).

461 '''

462 return self.ellipsoid.isProlate

463

464 @property_RO

465 def isSpherical(self):

466 '''Check whether this datum is (near-)spherical (C{bool}).

467 '''

468 return self.ellipsoid.isSpherical

469

470 def toStr(self, sep=_COMMASPACE_, name=NN, **unused): # PYCHOK expected

471 '''Return this datum as a string.

472

473 @kwarg sep: Separator to join (C{str}).

474 @kwarg name: Override name (C{str}) or C{None} to exclude

475 this datum's name.

476

477 @return: Datum attributes (C{str}).

478 '''

479 t = [] if name is None else \

480 [Fmt.EQUAL(name=repr(name or self.named))]

481 for a in (_ellipsoid_, _transform_):

482 v = getattr(self, a)

483 t.append(NN(Fmt.EQUAL(a, v.classname), _s_, _DOT_, v.name))

484 return sep.join(t)

485

486 @Property_RO

487 def transform(self):

488 '''Get this datum's transform (L{Transform}).

489 '''

490 return self._transform

491

492

493 def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...

494 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}

495 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).

496

497 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.

498 '''

499 if f is not None:

500 E, n, D = _EnD3((a_earth, f), name)

501 if raiser and not E:

502 raise _TypeError(f=f, **{raiser: a_earth})

503 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None):

504 return

505 elif isinstance(a_earth, Datum):

506 E, n, D = None, NN, a_earth

507 else:

508 E, n, D = _EnD3(a_earth, name)

509 if raiser and not E:

510 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})

511 if D is None:

512 D = Datum(E, transform=Transforms.Identity, name=_under(n))

513 inst._datum = D

514

515

516 def _earth_ellipsoid(earth, *name_raiser):

517 '''(INTERAL) Return the ellipsoid for the given C{earth} model.

518 '''

519 return Ellipsoids.Sphere if earth is R_M else (

520 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

521 _spherical_datum(earth, *name_raiser).ellipsoid)

522

523

524 def _ED2(radius, name):

525 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.

526 '''

527 D = Datums.Sphere

528 E = D.ellipsoid

529 if name or radius != E.a: # != E.b

530 n = _under(name)

531 E = Ellipsoid(radius, radius, name=n)

532 D = Datum(E, transform=Transforms.Identity, name=n)

533 return E, D

534

535

536 def _ellipsoidal_datum(earth, Error=TypeError, name=NN, raiser=NN):

537 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid} or L{Ellipsoid2},

538 C{a_f2Tuple}, 2-tuple or 2-list B{C{earth}} model.

539

540 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not ellipsoidal.

541 '''

542 if isinstance(earth, Datum):

543 D = earth

544 else:

545 E, n, D = _EnD3(earth, name)

546 if not E:

547 n = raiser or _earth_

548 _xinstanceof(Datum, Ellipsoid, Ellipsoid2, a_f2Tuple, **{n: earth})

549 if D is None:

550 D = Datum(E, transform=Transforms.Identity, name=_under(n))

551 if raiser and not D.isEllipsoidal:

552 raise _IsnotError(_ellipsoidal_, Error=Error, **{raiser: earth})

553 return D

554

555

556 def _EnD3(earth, name):

557 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.

558 '''

559 D = None

560 if isinstance(earth, (Ellipsoid, Ellipsoid2)):

561 E = earth

562 n = _under(name or E.name)

563 elif isinstance(earth, Datum):

564 E = earth.ellipsoid

565 n = _under(name or earth.name)

566 D = earth

567 elif _isRadius(earth):

568 E, D = _ED2(Radius_(earth), name)

569 n = E.name

570 elif isinstance(earth, a_f2Tuple):

571 n = _under(name or earth.name)

572 E = earth.ellipsoid(name=n)

573 elif islistuple(earth, minum=2):

574 E = Ellipsoids.Sphere

575 a, f = earth[:2]

576 if f or a != E.a: # != E.b

577 n = _under(name or _xattr(earth, name=NN))

578 E = Ellipsoid(a, f=f, name=n)

579 else:

580 n = E.name

581 D = Datums.Sphere

582 else:

583 E, n = None, NN

584 return E, n, D

585

586

587 def _mean_radius(radius, *lats):

588 '''(INTERNAL) Compute the mean radius of a L{Datum} from an L{Ellipsoid},

589 L{Ellipsoid2} or scalar earth C{radius} over several latitudes.

590 '''

591 if radius is R_M:

592 r = radius

593 elif _isRadius(radius):

594 r = Radius_(radius, low=0, Error=TypeError)

595 else:

596 E = _ellipsoidal_datum(radius).ellipsoid

597 r = fmean(map(E.Rgeocentric, lats)) if lats else E.Rmean

598 return r

599

600

601 def _minus(name): # in .trf

602 '''(INTERNAL) Name of C{inverse} Xform.

603 '''

604 m = _MINUS_

605 return name[len(m):] if name.startswith(m) else NN(m, name)

606

607

608 def _spherical_datum(earth, Error=TypeError, name=NN, raiser=NN):

609 '''(INTERNAL) Create a L{Datum} from an L{Ellipsoid}, L{Ellipsoid2},

610 C{a_f2Tuple}, 2-tuple, 2-list B{C{earth}} model or C{scalar} radius.

611

612 @kwarg raiser: If not C{NN}, raise an B{C{Error}} if not spherical.

613 '''

614 if _isRadius(earth):

615 _, D = _ED2(Radius_(earth, Error=Error), name)

616 else:

617 D = _ellipsoidal_datum(earth, Error=Error, name=name)

618 if raiser and not D.isSpherical:

619 raise _IsnotError(_spherical_, Error=Error, **{raiser: earth})

620 return D

621

622

623 class Datums(_NamedEnum):

624 '''(INTERNAL) L{Datum} registry, I{must} be a sub-class

625 to accommodate the L{_LazyNamedEnumItem} properties.

626 '''

627 def _Lazy(self, ellipsoid_name, transform_name, name=NN):

628 '''(INTERNAL) Instantiate the L{Datum}.

629 '''

630 return Datum(Ellipsoids.get(ellipsoid_name),

631 Transforms.get(transform_name), name=name)

632

633 Datums = Datums(Datum) # PYCHOK singleton

634'''Some pre-defined L{Datum}s, all I{lazily} instantiated.'''

635# Datums with associated ellipsoid and Helmert transform parameters

636# to convert from WGS84 into the given datum. More are available at

637# <https://Earth-Info.NGA.mil/GandG/coordsys/datums/NATO_DT.pdf> and

638# <XXX://www.FieldenMaps.info/cconv/web/cconv_params.js>.

639 Datums._assert(

640 # Belgian Datum 1972, based on Hayford ellipsoid.

641 # <https://NL.WikiPedia.org/wiki/Belgian_Datum_1972>

642 # <https://SpatialReference.org/ref/sr-org/7718/html/>

643 BD72 = _lazy(_BD72_, _Intl1924_, _BD72_),

644

645 # Germany <https://WikiPedia.org/wiki/Bessel-Ellipsoid>

646 # <https://WikiPedia.org/wiki/Helmert_transformation>

647 DHDN = _lazy(_DHDN_, _Bessel1841_, _DHDN_),

648

649 # <https://www.Gov.UK/guidance/oil-and-gas-petroleum-operations-notices#pon-4>

650 ED50 = _lazy(_ED50_, _Intl1924_, _ED50_),

651

652 # Australia <https://ICSM.Gov.AU/datum/gda2020-and-gda94-technical-manuals>

653# ADG66 = _lazy(_ADG66_, _ANS_, _WGS84_), # XXX Transform?

654# ADG84 = _lazy(_ADG84_, _ANS_, _WGS84_), # XXX Transform?

655# GDA94 = _lazy(_GDA94_, _GRS80_, _WGS84_),

656 GDA2020 = _lazy(_GDA2020_, _GRS80_, _WGS84_), # XXX Transform?

657

658 # <https://WikiPedia.org/wiki/GRS_80>

659 GRS80 = _lazy(_GRS80_, _GRS80_, _WGS84_),

660

661 # <https://OSI.IE/wp-content/uploads/2015/05/transformations_booklet.pdf> Table 2

662# Irl1975 = _lazy(_Irl1965_, _AiryModified_, _Irl1965_),

663 Irl1975 = _lazy(_Irl1975_, _AiryModified_, _Irl1975_),

664

665 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>

666 Krassovski1940 = _lazy(_Krassovski1940_, _Krassovski1940_, _Krassovski1940_), # XXX spelling?

667 Krassowsky1940 = _lazy(_Krassowsky1940_, _Krassowsky1940_, _Krassowsky1940_), # XXX spelling?

668

669 # Austria <https://DE.WikiPedia.org/wiki/Datum_Austria>

670 MGI = _lazy(_MGI_, _Bessel1841_, _MGI_),

671

672 # <https://WikiPedia.org/wiki/Helmert_transformation>

673 NAD27 = _lazy(_NAD27_, _Clarke1866_, _NAD27_),

674

675 # NAD83 (2009) == WGS84 - <https://www.UVM.edu/giv/resources/WGS84_NAD83.pdf>

676 # (If you *really* must convert WGS84<->NAD83, you need more than this!)

677 NAD83 = _lazy(_NAD83_, _GRS80_, _NAD83_),

678

679 # Nouvelle Triangulation Francaise (Paris) XXX verify

680 NTF = _lazy(_NTF_, _Clarke1880IGN_, _NTF_),

681

682 # <https://www.OrdnanceSurvey.co.UK/docs/support/guide-coordinate-systems-great-britain.pdf>

683 OSGB36 = _lazy(_OSGB36_, _Airy1830_, _OSGB36_),

684

685 # Germany <https://WikiPedia.org/wiki/Helmert_transformation>

686 Potsdam = _lazy(_Potsdam_, _Bessel1841_, _Bessel1841_),

687

688 # XXX psuedo-ellipsoids for spherical LatLon

689 Sphere = _lazy(_Sphere_, _Sphere_, _WGS84_),

690

691 # <https://www.GeoCachingToolbox.com?page=datumEllipsoidDetails>

692 TokyoJapan = _lazy(_TokyoJapan_, _Bessel1841_, _TokyoJapan_),

693

694 # <https://www.ICAO.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf>

695 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

696

697 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

698)

699

700_WGS84 = Datums.WGS84

701 assert _WGS84.ellipsoid is _EWGS84

702# assert _WGS84.transform.isunity

703

704 if __name__ == '__main__':

705

706 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

707 from pygeodesy.lazily import printf

708

709 # __doc__ of this file, force all into registery

710 for r in (Datums, Transforms):

711 t = [NN] + r.toRepr(all=True, asorted=True).split(_NL_)

712 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

713

714# **) MIT License

715#

716# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

717#

718# Permission is hereby granted, free of charge, to any person obtaining a

719# copy of this software and associated documentation files (the "Software"),

720# to deal in the Software without restriction, including without limitation

721# the rights to use, copy, modify, merge, publish, distribute, sublicense,

722# and/or sell copies of the Software, and to permit persons to whom the

723# Software is furnished to do so, subject to the following conditions:

724#

725# The above copyright notice and this permission notice shall be included

726# in all copies or substantial portions of the Software.

727#

728# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

729# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

730# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

731# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

732# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

733# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

734# OTHER DEALINGS IN THE SOFTWARE.