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, _zip

71 from pygeodesy.constants import R_M, _float as _F, _0_0, _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, _operator

79 from pygeodesy.internals import _passarg, _under

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

81 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

82 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _MINUS_, \

83 _Krassovski1940_, _Krassowsky1940_, _NAD27_, _NAD83_, _s_, \

84 _PLUS_, _Sphere_, _spherical_, _transform_, _UNDER_, \

85 _WGS72_, _WGS84_

86 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

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

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

89 from pygeodesy.props import Property_RO, property_RO

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

91 from pygeodesy.units import _isRadius, Radius_, radians

92

93# from math import radians # from .units

94# import operator as _operator # from .fmath

95

96__all__ = _ALL_LAZY.datums

97__version__ = '24.05.14'

98

99_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

100_BD72_ = 'BD72'

101_DHDN_ = 'DHDN'

102_DHDNE_ = 'DHDNE'

103_DHDNW_ = 'DHDNW'

104_ED50_ = 'ED50'

105_GDA2020_ = 'GDA2020' # in .trf

106_Identity_ = 'Identity'

107_Irl1965_ = 'Irl1965'

108_Irl1975_ = 'Irl1975'

109_MGI_ = 'MGI'

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

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

112_NTF_ = 'NTF'

113_OSGB36_ = 'OSGB36'

114_Potsdam_ = 'Potsdam'

115_RPS = radians(_1_0 / _3600_0) # radians per arc-second

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

117_TokyoJapan_ = 'TokyoJapan'

118

119

120 class Transform(_NamedEnumItem):

121 '''Helmert I{datum} transformation.

122

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

124 '''

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

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

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

128

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

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

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

132

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

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

135

136 sx = _0_0 # x rotation (C{arc-seconds})

137 sy = _0_0 # y rotation (C{arc-seconds})

138 sz = _0_0 # z rotation (C{arc-seconds})

139

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

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

142 '''New L{Transform}.

143

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

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

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

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

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

149 @kwarg sx: X rotation (C{arc-seconds}).

150 @kwarg sy: Y rotation (C{arc-seconds}).

151 @kwarg sz: Z rotation (C{arc-seconds}).

152

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

154 '''

155 if tx:

156 self.tx = tx

157 if ty:

158 self.ty = ty

159 if tz:

160 self.tz = tz

161 if s:

162 self.s = s

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

164 if sx: # secs to rads

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

166 if sy:

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

168 if sz:

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

170

171 self._register(Transforms, name)

172

173 def __eq__(self, other):

174 '''Compare this and an other transform.

175

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

177

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

179 '''

180 return self is other or (isinstance(other, Transform)

181 and _equall(other, self))

182

183 def __hash__(self):

184 return hash(tuple(self))

185

186 def __iter__(self):

187 '''Yield the initial attribute values, I{in order}.

188 '''

189 for n in _Names7:

190 yield getattr(self, n)

191

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

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

194

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

196

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

198

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

200 '''

201 _ = _xellipsoidall(point)

202 return point.toTransform(self)

203

204 def __neg__(self):

205 return self.inverse()

206

207 def inverse(self, name=NN):

208 '''Return the inverse of this transform.

209

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

211

212 @return: Inverse (L{Transform}), unregistered.

213 '''

214 r = type(self)(**dict(self.items(inverse=True)))

215 n = name or _negastr(self.name)

216 if n:

217 r.name = n # unregistered

218 return r

219

220 @Property_RO

221 def isunity(self):

222 '''Is this a C{unity, identity} transform (C{bool}), like

223 WGS84 with translation, scale and rotation all zero?

224 '''

225 return not any(self)

226

227 def items(self, inverse=False):

228 '''Yield the initial attributes, each as 2-tuple C{(name, value)}.

229

230 @kwarg inverse: If C{True}, negate the values (C{bool}).

231 '''

232 _p = neg if inverse else _passarg

233 for n, x in _zip(_Names7, self):

234 yield n, _p(x)

235

236 def _rps2(self, s_):

237 '''(INTERNAL) Rotation in C{radians} and C{arc-seconds}.

238 '''

239 # _MR == _RPS * 1.e-3 # radians per milli-arc-second, equ (2)

240 # <https://www.NGS.NOAA.gov/CORS/Articles/SolerSnayASCE.pdf>

241 return (_RPS * s_), s_

242

243 def _s_s1(self, s1): # in .trf

244 '''(INTERNAL) Set C{s1} and C{s}.

245 '''

246 Transform.isunity._update(self)

247 self.s1 = s1

248 self.s = s = (s1 - _1_0) / _S1_S

249 return s

250

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

252 '''Return this transform as a string.

253

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

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

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

257 this transform's name.

258

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

260 '''

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

262

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

264 '''Transform a (cartesian) position, forward or inverse.

265

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

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

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

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

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

271 cartesian class and additional C{B{Vector}} keyword

272 arguments to return the transformed position.

273

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

275 unless some B{C{Vector_and_kwds}} are specified.

276 '''

277 if self.isunity:

278 r = Vector3Tuple(x, y, z, name=self.name) # == inverse

279 else:

280 xyz1 = x, y, z, _1_0

281 s1 = self.s1

282 if inverse:

283 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

291 name=self.name)

292 if Vector_and_kwds:

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

294 if V:

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

296 return r

297

298

299 class Transforms(_NamedEnum):

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

301 to accommodate the L{_LazyNamedEnumItem} properties.

302 '''

303 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):

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

305 '''

306 return Transform(**name_tx_ty_tz_s_sx_sy_sz)

307

308 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

311 Transforms._assert(

312 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893), s=_F(1.2727),

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

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

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

316

317 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0), s=_F(-8.3),

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

319

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

321

322 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39), s=_F(-9.82),

323 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458)), # Germany

324

325 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2), s=_F(-2.55),

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

327 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797)), # East Germany

328

329 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2), s=_F(-6.7),

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

331 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455)), # West Germany

332

333 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1), s=_F(-1.2),

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

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

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

337 sz=_F( 0.156)),

338 Identity = _lazy(_Identity_),

339

340 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),

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

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

343 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), s=_F(-8.15),

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

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

346

347 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),

348 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling

349

350 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), s=_F(-2.423),

351 sx=_F( -0.02), sy=_F( 0.26), sz=_F( 0.13)), # spelling

352

353 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920), s=_F(-2.423),

354 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297)), # Austria

355

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

357

358 NAD83 = _lazy(_NAD83_, tx=_F(1.004), ty=_F(-1.910), tz=_F(-0.515), s=_F(-0.0015),

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

360

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

362

363 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060), s=_F(20.4894),

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

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

366

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

368

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

370

371 WGS84 = _lazy(_WGS84_), # unity

372)

373

374

375 class Datum(_NamedEnumItem):

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

377 '''

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

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

380

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

382 '''New L{Datum}.

383

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

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

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

387

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

389

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

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

392 not a L{Transform}.

393 '''

394 self._ellipsoid = ellipsoid or Datum._ellipsoid

395 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

396

397 self._transform = transform or Datum._transform

398 _xinstanceof(Transform, transform=self.transform)

399

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

401

402 def __eq__(self, other):

403 '''Compare this and an other datum.

404

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

406

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

408 '''

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

410 self.ellipsoid == other.ellipsoid and

411 self.transform == other.transform)

412

413 def __hash__(self):

414 return self._hash # memoized

415

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

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

418

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

420 '''

421 _ = _xellipsoidall(point)

422 return point.toDatum(self)

423

424 def ecef(self, Ecef=None):

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

426

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

428

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

430

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

432

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

434 '''

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

436

437 @Property_RO

438 def ellipsoid(self):

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

440 '''

441 return self._ellipsoid

442

443 @Property_RO

444 def exactTM(self):

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

446 '''

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

448

449 @Property_RO

450 def _hash(self):

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

452

453 @property_RO

454 def isEllipsoidal(self):

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

456 '''

457 return self.ellipsoid.isEllipsoidal

458

459 @property_RO

460 def isOblate(self):

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

462 '''

463 return self.ellipsoid.isOblate

464

465 @property_RO

466 def isProlate(self):

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

468 '''

469 return self.ellipsoid.isProlate

470

471 @property_RO

472 def isSpherical(self):

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

474 '''

475 return self.ellipsoid.isSpherical

476

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

478 '''Return this datum as a string.

479

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

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

482 this datum's name.

483

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

485 '''

486 t = [] if name is None else \

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

488 for a in (_ellipsoid_, _transform_):

489 v = getattr(self, a)

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

491 return sep.join(t)

492

493 @Property_RO

494 def transform(self):

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

496 '''

497 return self._transform

498

499

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

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

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

503

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

505 '''

506 if f is not None:

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

508 if raiser and not E:

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

510 elif a_earth in (_EWGS84, _WGS84, None) and inst._datum is _WGS84:

511 return

512 elif isinstance(a_earth, Datum):

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

514 else:

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

516 if raiser and not E:

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

518 if D is None:

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

520 inst._datum = D

521

522

523 def _earth_ellipsoid(earth, *name_raiser):

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

525 '''

526 return Ellipsoids.Sphere if earth is R_M else (

527 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

528 _spherical_datum(earth, *name_raiser).ellipsoid)

529

530

531 def _ED2(radius, name):

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

533 '''

534 D = Datums.Sphere

535 E = D.ellipsoid

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

537 n = _under(name)

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

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

540 return E, D

541

542

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

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

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

546

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

548 '''

549 if isinstance(earth, Datum):

550 D = earth

551 else:

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

553 if not E:

554 n = raiser or _earth_

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

556 if D is None:

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

558 if raiser and not D.isEllipsoidal:

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

560 return D

561

562

563 def _EnD3(earth, name):

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

565 '''

566 D = None

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

568 E = earth

569 n = _under(name or E.name)

570 elif isinstance(earth, Datum):

571 E = earth.ellipsoid

572 n = _under(name or earth.name)

573 D = earth

574 elif _isRadius(earth):

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

576 n = E.name

577 elif isinstance(earth, a_f2Tuple):

578 n = _under(name or earth.name)

579 E = earth.ellipsoid(name=n)

580 elif islistuple(earth, minum=2):

581 E = Ellipsoids.Sphere

582 a, f = earth[:2]

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

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

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

586 else:

587 n = E.name

588 D = Datums.Sphere

589 else:

590 E, n = None, NN

591 return E, n, D

592

593

594 def _equall(t1, t2): # in .trf

595 '''(INTERNAL) Return L{Transform} C{t1 == t2}.

596 '''

597 return all(map(_operator.eq, t1, t2))

598

599

600 def _mean_radius(radius, *lats):

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

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

603 '''

604 if radius is R_M:

605 r = radius

606 elif _isRadius(radius):

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

608 else:

609 E = _ellipsoidal_datum(radius).ellipsoid

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

611 return r

612

613

614 def _negastr(name): # in .trf

615 '''(INTERNAL) Negate a C{Transform/-Xform} name.

616 '''

617 b, m, p = _BAR_, _MINUS_, _PLUS_

618 n = name.replace(m, b).replace(p, m).replace(b, p)

619 # as good and fast as (in Python 3+ only) ...

620 # _MINUSxPLUS = str.maketrans({_MINUS_: _PLUS_, _PLUS_: _MINUS_})

621 # def _negastr(name):

622 # n = name.translate(_MINUSxPLUS)

623 # ...

624 return n.lstrip(p) if n.startswith(p) else NN(m, n)

625

626

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

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

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

630

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

632 '''

633 if _isRadius(earth):

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

635 else:

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

637 if raiser and not D.isSpherical:

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

639 return D

640

641

642 class Datums(_NamedEnum):

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

644 to accommodate the L{_LazyNamedEnumItem} properties.

645 '''

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

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

648 '''

649 return Datum(Ellipsoids.get(ellipsoid_name),

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

651

652 Datums = Datums(Datum) # PYCHOK singleton

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

654# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

658 Datums._assert(

659 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

663

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

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

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

667

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

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

670

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

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

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

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

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

676

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

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

679

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

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

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

683

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

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

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

687

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

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

690

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

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

693

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

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

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

697

698 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

700

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

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

703

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

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

706

707 # XXX psuedo-ellipsoids for spherical LatLon

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

709

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

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

712

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

714 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

715

716 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

717)

718

719_WGS84 = Datums.WGS84

720 assert _WGS84.ellipsoid is _EWGS84

721# assert _WGS84.transform.isunity

722

723 if __name__ == '__main__':

724

725 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

726 from pygeodesy import printf

727

728 # __doc__ of this file, force all into registery

729 for r in (Datums, Transforms):

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

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

732

733# **) MIT License

734#

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

736#

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

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

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

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

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

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

743#

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

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

746#

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

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

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

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

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

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

753# OTHER DEALINGS IN THE SOFTWARE.