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, _passarg, _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.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_, _PLUS_, _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# import operator as _operator # from .fmath

94

95__all__ = _ALL_LAZY.datums

96__version__ = '24.02.29'

97

98_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

99_BD72_ = 'BD72'

100_DHDN_ = 'DHDN'

101_DHDNE_ = 'DHDNE'

102_DHDNW_ = 'DHDNW'

103_ED50_ = 'ED50'

104_GDA2020_ = 'GDA2020' # in .trf

105_Identity_ = 'Identity'

106_Irl1965_ = 'Irl1965'

107_Irl1975_ = 'Irl1975'

108_MGI_ = 'MGI'

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

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

111_NTF_ = 'NTF'

112_OSGB36_ = 'OSGB36'

113_Potsdam_ = 'Potsdam'

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

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

116_TokyoJapan_ = 'TokyoJapan'

117

118

119 class Transform(_NamedEnumItem):

120 '''Helmert I{datum} transformation.

121

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

123 '''

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

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

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

127

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

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

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

131

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

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

134

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

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

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

138

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

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

141 '''New L{Transform}.

142

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

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

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

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

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

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

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

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

151

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

153 '''

154 if tx:

155 self.tx = tx

156 if ty:

157 self.ty = ty

158 if tz:

159 self.tz = tz

160 if s:

161 self.s = s

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

163 if sx: # secs to rads

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

165 if sy:

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

167 if sz:

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

169

170 self._register(Transforms, name)

171

172 def __eq__(self, other):

173 '''Compare this and an other transform.

174

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

176

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

178 '''

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

180 and _equall(other, self))

181

182 def __hash__(self):

183 return hash(tuple(self))

184

185 def __iter__(self):

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

187 '''

188 for n in _Names7:

189 yield getattr(self, n)

190

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

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

193

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

195

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

197

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

199 '''

200 _ = _xellipsoidall(point)

201 return point.toTransform(self)

202

203 def __neg__(self):

204 return self.inverse()

205

206 def inverse(self, name=NN):

207 '''Return the inverse of this transform.

208

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

210

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

212 '''

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

214 n = name or _negstr(self.name)

215 if n:

216 r.name = n # unregistered

217 return r

218

219 @Property_RO

220 def isunity(self):

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

222 WGS84 with translation, scale and rotation all zero?

223 '''

224 return not any(self)

225

226 def items(self, inverse=False):

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

228

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

230 '''

231 _p = neg if inverse else _passarg

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

233 yield n, _p(x)

234

235 def _rps2(self, s_):

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

237 '''

238 return (_RPS * s_), s_

239

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

241 '''Return this transform as a string.

242

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

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

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

246 this transform's name.

247

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

249 '''

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

251

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

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

254

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

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

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

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

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

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

261 arguments to return the transformed position.

262

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

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

265 '''

266 if self.isunity:

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

268 else:

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 r = 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 if Vector_and_kwds:

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

283 if V:

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

285 return r

286

287

288 class Transforms(_NamedEnum):

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

290 to accommodate the L{_LazyNamedEnumItem} properties.

291 '''

292 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):

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

294 '''

295 return Transform(**name_tx_ty_tz_s_sx_sy_sz)

296

297 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

300 Transforms._assert(

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

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

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

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

305

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

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

308

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

310

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

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

313

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

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

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

317

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

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

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

321

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

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

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

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

326 sz=_F( 0.156)),

327 Identity = _lazy(_Identity_),

328

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

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

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

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

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

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

335

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

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

338

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

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

341

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

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

344

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

346

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

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

349

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

351

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

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

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

355

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

357

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

359

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 in (_EWGS84, _WGS84, None) and inst._datum is _WGS84:

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 _equall(t1, t2): # in .trf

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

585 '''

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

587

588

589 def _mean_radius(radius, *lats):

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

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

592 '''

593 if radius is R_M:

594 r = radius

595 elif _isRadius(radius):

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

597 else:

598 E = _ellipsoidal_datum(radius).ellipsoid

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

600 return r

601

602

603 try:

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

605

606 def _negstr(name): # in .trf

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

608 '''

609 n = name.translate(_MINUSxPLUS)

610 return n.lstrip(_PLUS_) if n.startswith(_PLUS_) else NN(_MINUS_, n)

611

612 except AttributeError: # no str.maketran in Python 2-

613 from pygeodesy.interns import _BAR_

614

615 def _negstr(name): # PYCHOK in .trf

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

617 '''

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

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

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

621

622

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

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

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

626

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

628 '''

629 if _isRadius(earth):

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

631 else:

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

633 if raiser and not D.isSpherical:

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

635 return D

636

637

638 class Datums(_NamedEnum):

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

640 to accommodate the L{_LazyNamedEnumItem} properties.

641 '''

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

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

644 '''

645 return Datum(Ellipsoids.get(ellipsoid_name),

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

647

648 Datums = Datums(Datum) # PYCHOK singleton

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

650# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

654 Datums._assert(

655 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

659

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

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

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

663

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

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

666

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

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

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

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

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

672

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

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

675

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

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

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

679

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

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

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

683

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

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

686

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

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

689

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

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

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

693

694 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

696

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

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

699

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

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

702

703 # XXX psuedo-ellipsoids for spherical LatLon

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

705

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

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

708

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

710 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

711

712 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

713)

714

715_WGS84 = Datums.WGS84

716 assert _WGS84.ellipsoid is _EWGS84

717# assert _WGS84.transform.isunity

718

719 if __name__ == '__main__':

720

721 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

722 from pygeodesy.lazily import printf

723

724 # __doc__ of this file, force all into registery

725 for r in (Datums, Transforms):

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

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

728

729# **) MIT License

730#

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

732#

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

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

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

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

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

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

739#

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

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

742#

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

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

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

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

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

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

749# OTHER DEALINGS IN THE SOFTWARE.