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.18'

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 self._hash # memoized

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 '''Convert an I{ellipsoidal} B{C{point}} with the Helmert transform.

192

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

194 '''

195 _ = _xellipsoidall(point)

196 return point.toTransform(self)

197

198 def __neg__(self):

199 return self.inverse()

200

201# def __sub__(self, other):

202# _xinstanceof(Transform, other=other)

203#

204# def _sub(a, b):

205# for n in _Names11:

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

207#

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

209

210 @Property_RO

211 def _hash(self):

212 return hash(x for x in self)

213

214 def items(self):

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

216 '''

217 for n in _Names7:

218 yield n, getattr(self, n)

219

220 def inverse(self, name=NN):

221 '''Return the inverse of this transform.

222

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

224

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

226

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

228 '''

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

230 n = name or _minus(self.name)

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

232

233 @Property_RO

234 def isunity(self):

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

236 '''

237 return not any(s for s in self)

238

239 def _rps2(self, s_):

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

241 '''

242 return (_RPS * s_), s_

243

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

245 '''Return this transform as a string.

246

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

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

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

250 this transform's name.

251

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

253 '''

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

255

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

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

258

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

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

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

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

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

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

265 keyword arguments to return the transformed

266 point.

267

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

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

270 '''

271 if self.isunity:

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

273 else:

274 xyz1 = x, y, z, _1_0

275 s1 = self.s1

276 if inverse:

277 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

285 name=self.name)

286 if Vector_and_kwds:

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

288 if V:

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

290 return r

291

292

293 class Transforms(_NamedEnum):

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

295 to accommodate the L{_LazyNamedEnumItem} properties.

296 '''

297 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

299 '''

300 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

301

302 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

305 Transforms._assert(

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

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

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

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

310 s=_F( 1.2727)),

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

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

313 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

331 Identity = _lazy(_Identity_),

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

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

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

335 s=_F( -8.15)),

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

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

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

339 s=_F( -8.15)),

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

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

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

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

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

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

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

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

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

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

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

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

352 s=_F(-0.0015)),

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

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

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

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

357 s=_F( 20.4894)),

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

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

360 WGS84 = _lazy(_WGS84_), # unity

361)

362

363

364 class Datum(_NamedEnumItem):

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

366 '''

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

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

369

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

371 '''New L{Datum}.

372

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

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

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

376

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

378

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

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

381 not a L{Transform}.

382 '''

383 self._ellipsoid = ellipsoid or Datum._ellipsoid

384 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

385

386 self._transform = transform or Datum._transform

387 _xinstanceof(Transform, transform=self.transform)

388

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

390

391 def __eq__(self, other):

392 '''Compare this and an other datum.

393

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

395

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

397 '''

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

399 self.ellipsoid == other.ellipsoid and

400 self.transform == other.transform)

401

402 def __hash__(self):

403 return self._hash # memoized

404

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

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

407

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

409 '''

410 _ = _xellipsoidall(point)

411 return point.toDatum(self)

412

413 def ecef(self, Ecef=None):

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

415

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

417

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

419

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

421

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

423 '''

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

425

426 @Property_RO

427 def ellipsoid(self):

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

429 '''

430 return self._ellipsoid

431

432 @Property_RO

433 def exactTM(self):

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

435 '''

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

437

438 @Property_RO

439 def _hash(self):

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

441

442 @property_RO

443 def isEllipsoidal(self):

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

445 '''

446 return self.ellipsoid.isEllipsoidal

447

448 @property_RO

449 def isOblate(self):

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

451 '''

452 return self.ellipsoid.isOblate

453

454 @property_RO

455 def isProlate(self):

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

457 '''

458 return self.ellipsoid.isProlate

459

460 @property_RO

461 def isSpherical(self):

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

463 '''

464 return self.ellipsoid.isSpherical

465

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

467 '''Return this datum as a string.

468

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

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

471 this datum's name.

472

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

474 '''

475 t = [] if name is None else \

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

477 for a in (_ellipsoid_, _transform_):

478 v = getattr(self, a)

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

480 return sep.join(t)

481

482 @Property_RO

483 def transform(self):

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

485 '''

486 return self._transform

487

488

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

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

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

492

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

494 '''

495 if f is not None:

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

497 if raiser and not E:

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

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

500 return

501 elif isinstance(a_earth, Datum):

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

503 else:

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

505 if raiser and not E:

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

507 if D is None:

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

509 inst._datum = D

510

511

512 def _earth_ellipsoid(earth, *name_raiser):

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

514 '''

515 return Ellipsoids.Sphere if earth is R_M else (

516 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

517 _spherical_datum(earth, *name_raiser).ellipsoid)

518

519

520 def _ED2(radius, name):

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

522 '''

523 D = Datums.Sphere

524 E = D.ellipsoid

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

526 n = _under(name)

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

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

529 return E, D

530

531

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

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

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

535

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

537 '''

538 if isinstance(earth, Datum):

539 D = earth

540 else:

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

542 if not E:

543 n = raiser or _earth_

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

545 if D is None:

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

547 if raiser and not D.isEllipsoidal:

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

549 return D

550

551

552 def _EnD3(earth, name):

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

554 '''

555 D = None

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

557 E = earth

558 n = _under(name or E.name)

559 elif isinstance(earth, Datum):

560 E = earth.ellipsoid

561 n = _under(name or earth.name)

562 D = earth

563 elif _isRadius(earth):

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

565 n = E.name

566 elif isinstance(earth, a_f2Tuple):

567 n = _under(name or earth.name)

568 E = earth.ellipsoid(name=n)

569 elif islistuple(earth, minum=2):

570 E = Ellipsoids.Sphere

571 a, f = earth[:2]

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

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

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

575 else:

576 n = E.name

577 D = Datums.Sphere

578 else:

579 E, n = None, NN

580 return E, n, D

581

582

583 def _mean_radius(radius, *lats):

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

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

586 '''

587 if radius is R_M:

588 r = radius

589 elif _isRadius(radius):

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

591 else:

592 E = _ellipsoidal_datum(radius).ellipsoid

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

594 return r

595

596

597 def _minus(name): # in .trf

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

599 '''

600 m = _MINUS_

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

602

603

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

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

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

607

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

609 '''

610 if _isRadius(earth):

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

612 else:

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

614 if raiser and not D.isSpherical:

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

616 return D

617

618

619 class Datums(_NamedEnum):

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

621 to accommodate the L{_LazyNamedEnumItem} properties.

622 '''

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

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

625 '''

626 return Datum(Ellipsoids.get(ellipsoid_name),

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

628

629 Datums = Datums(Datum) # PYCHOK singleton

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

631# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

635 Datums._assert(

636 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

640

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

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

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

644

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

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

647

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

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

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

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

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

653

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

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

656

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

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

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

660

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

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

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

664

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

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

667

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

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

670

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

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

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

674

675 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

677

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

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

680

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

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

683

684 # XXX psuedo-ellipsoids for spherical LatLon

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

686

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

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

689

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

691 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

692

693 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

694)

695

696_WGS84 = Datums.WGS84

697 assert _WGS84.ellipsoid is _EWGS84

698# assert _WGS84.transform.isunity

699

700 if __name__ == '__main__':

701

702 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

703 from pygeodesy.lazily import printf

704

705 # __doc__ of this file, force all into registery

706 for r in (Datums, Transforms):

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

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

709

710# **) MIT License

711#

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

713#

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

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

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

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

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

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

720#

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

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

723#

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

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

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

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

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

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

730# OTHER DEALINGS IN THE SOFTWARE.