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

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

77 from pygeodesy.fmath import fdot, fmean, Fmt

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

79 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \

80 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \

81 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \

82 _NAD83_, _s_, _Sphere_, _spherical_, _transform_, \

83 _UNDER_, _WGS72_, _WGS84_, _under

84 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

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

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

87 from pygeodesy.props import Property_RO, property_RO

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

89 from pygeodesy.units import _isRadius, Radius_, radians

90

91# from math import radians # from .units

92

93__all__ = _ALL_LAZY.datums

94__version__ = '24.02.12'

95

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_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf

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

109_NTF_ = 'NTF'

110_OSGB36_ = 'OSGB36'

111_Potsdam_ = 'Potsdam'

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

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

114_TokyoJapan_ = 'TokyoJapan'

115

116

117 class Transform(_NamedEnumItem):

118 '''Helmert I{datum} transformation.

119

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

121 '''

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

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

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

125

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

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

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

129

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

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

132

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

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

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

136

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

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

139 '''New L{Transform}.

140

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

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

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

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

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

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

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

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

149

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

151 '''

152 if tx:

153 self.tx = tx

154 if ty:

155 self.ty = ty

156 if tz:

157 self.tz = tz

158 if sx: # secs to rads

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

160 if sy:

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

162 if sz:

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

164 if s:

165 self.s = s

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

167

168 self._register(Transforms, name)

169

170 def __eq__(self, other):

171 '''Compare this and an other transform.

172

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

174

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

176 '''

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

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

179

180 def __hash__(self):

181 return self._hash # memoized

182

183 def __iter__(self):

184 '''Yield the initial attribute values.

185 '''

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

187 yield x

188

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

190 '''Helmert-transform an I{ellipsoidal} cartesian B{C{point}}.

191

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

193 '''

194 m = _MODS.ellipsoidalBase

195 _xinstanceof(m.CartesianEllipsoidalBase, point=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):

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: Optional direction, forward or inverse (C{bool}).

263

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

265 '''

266 if self.isunity:

267 return Vector3Tuple(x, y, z, name=self.name)

268

269 xyz1 = x, y, z, _1_0

270 s1 = self.s1

271 if inverse:

272 xyz1 = map2(neg, xyz1)

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

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

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

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

277 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx),

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

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

280 name=self.name)

281

282

283 class Transforms(_NamedEnum):

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

285 to accommodate the L{_LazyNamedEnumItem} properties.

286 '''

287 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

289 '''

290 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

291

292 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

295 Transforms._assert(

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

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

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

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

300 s=_F( 1.2727)),

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

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

303 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

321 Identity = _lazy(_Identity_),

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

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

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

325 s=_F( -8.15)),

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

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

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

329 s=_F( -8.15)),

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

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

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

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

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

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

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

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

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

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

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

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

342 s=_F(-0.0015)),

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

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

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

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

347 s=_F( 20.4894)),

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

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

350 WGS84 = _lazy(_WGS84_), # unity

351)

352

353

354 class Datum(_NamedEnumItem):

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

356 '''

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

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

359

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

361 '''New L{Datum}.

362

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

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

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

366

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

368

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

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

371 not a L{Transform}.

372 '''

373 self._ellipsoid = ellipsoid or Datum._ellipsoid

374 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

375

376 self._transform = transform or Datum._transform

377 _xinstanceof(Transform, transform=self.transform)

378

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

380

381 def __eq__(self, other):

382 '''Compare this and an other datum.

383

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

385

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

387 '''

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

389 self.ellipsoid == other.ellipsoid and

390 self.transform == other.transform)

391

392 def __hash__(self):

393 return self._hash # memoized

394

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

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

397

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

399 '''

400 m = _MODS.ellipsoidalBase

401 _xinstanceof(m.CartesianEllipsoidalBase, m.LatLonEllipsoidalBase, point=point)

402 return point.toDatum(self)

403

404 def ecef(self, Ecef=None):

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

406

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

408

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

410

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

412

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

414 '''

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

416

417 @Property_RO

418 def ellipsoid(self):

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

420 '''

421 return self._ellipsoid

422

423 @Property_RO

424 def exactTM(self):

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

426 '''

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

428

429 @Property_RO

430 def _hash(self):

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

432

433 @property_RO

434 def isEllipsoidal(self):

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

436 '''

437 return self.ellipsoid.isEllipsoidal

438

439 @property_RO

440 def isOblate(self):

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

442 '''

443 return self.ellipsoid.isOblate

444

445 @property_RO

446 def isProlate(self):

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

448 '''

449 return self.ellipsoid.isProlate

450

451 @property_RO

452 def isSpherical(self):

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

454 '''

455 return self.ellipsoid.isSpherical

456

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

458 '''Return this datum as a string.

459

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

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

462 this datum's name.

463

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

465 '''

466 t = [] if name is None else \

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

468 for a in (_ellipsoid_, _transform_):

469 v = getattr(self, a)

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

471 return sep.join(t)

472

473 @Property_RO

474 def transform(self):

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

476 '''

477 return self._transform

478

479

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

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

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

483

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

485 '''

486 if f is not None:

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

488 if raiser and not E:

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

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

491 return

492 elif isinstance(a_earth, Datum):

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

494 else:

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

496 if raiser and not E:

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

498 if D is None:

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

500 inst._datum = D

501

502

503 def _earth_ellipsoid(earth, *name_raiser):

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

505 '''

506 return Ellipsoids.Sphere if earth is R_M else (

507 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

508 _spherical_datum(earth, *name_raiser).ellipsoid)

509

510

511 def _ED2(radius, name):

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

513 '''

514 D = Datums.Sphere

515 E = D.ellipsoid

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

517 n = _under(name)

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

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

520 return E, D

521

522

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

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

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

526

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

528 '''

529 if isinstance(earth, Datum):

530 D = earth

531 else:

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

533 if not E:

534 n = raiser or _earth_

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

536 if D is None:

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

538 if raiser and not D.isEllipsoidal:

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

540 return D

541

542

543 def _EnD3(earth, name):

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

545 '''

546 D = None

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

548 E = earth

549 n = _under(name or E.name)

550 elif isinstance(earth, Datum):

551 E = earth.ellipsoid

552 n = _under(name or earth.name)

553 D = earth

554 elif _isRadius(earth):

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

556 n = E.name

557 elif isinstance(earth, a_f2Tuple):

558 n = _under(name or earth.name)

559 E = earth.ellipsoid(name=n)

560 elif islistuple(earth, minum=2):

561 E = Ellipsoids.Sphere

562 a, f = earth[:2]

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

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

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

566 else:

567 n = E.name

568 D = Datums.Sphere

569 else:

570 E, n = None, NN

571 return E, n, D

572

573

574 def _mean_radius(radius, *lats):

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

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

577 '''

578 if radius is R_M:

579 r = radius

580 elif _isRadius(radius):

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

582 else:

583 E = _ellipsoidal_datum(radius).ellipsoid

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

585 return r

586

587

588 def _minus(name): # in .trf

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

590 '''

591 m = _MINUS_

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

593

594

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

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

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

598

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

600 '''

601 if _isRadius(earth):

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

603 else:

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

605 if raiser and not D.isSpherical:

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

607 return D

608

609

610 class Datums(_NamedEnum):

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

612 to accommodate the L{_LazyNamedEnumItem} properties.

613 '''

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

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

616 '''

617 return Datum(Ellipsoids.get(ellipsoid_name),

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

619

620 Datums = Datums(Datum) # PYCHOK singleton

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

622# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

626 Datums._assert(

627 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

631

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

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

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

635

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

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

638

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

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

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

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

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

644

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

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

647

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

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

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

651

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

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

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

655

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

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

658

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

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

661

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

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

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

665

666 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

668

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

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

671

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

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

674

675 # XXX psuedo-ellipsoids for spherical LatLon

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

677

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

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

680

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

682 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

683

684 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

685)

686

687_WGS84 = Datums.WGS84

688 assert _WGS84.ellipsoid is _EWGS84

689# assert _WGS84.transform.isunity

690

691 if __name__ == '__main__':

692

693 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

694 from pygeodesy.lazily import printf

695

696 # __doc__ of this file, force all into registery

697 for r in (Datums, Transforms):

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

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

700

701# **) MIT License

702#

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

704#

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

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

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

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

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

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

711#

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

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

714#

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

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

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

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

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

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

721# OTHER DEALINGS IN THE SOFTWARE.