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, _xellipsoidall, _xkwds, _xkwds_pop2

77 from pygeodesy.fmath import fdot, fmean, Fmt, _operator

78 from pygeodesy.internals import _passarg, _under

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

80 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

81 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _MINUS_, \

82 _Krassovski1940_, _Krassowsky1940_, _NAD27_, _NAD83_, _s_, \

83 _PLUS_, _Sphere_, _spherical_, _transform_, _UNDER_, \

84 _WGS72_, _WGS84_

85 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

86 from pygeodesy.named import _lazyNamedEnumItem as _lazy, _name__, _name2__, _NamedEnum, \

87 _NamedEnumItem

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

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):

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__(**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): # 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: Optional, override C{B{name}=NN} (C{str}) or

257 C{None} to exclude this transform's name.

258

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

260 '''

261 name, _ = _name2__(**name) # name=None

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

263

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

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

266

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

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

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

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

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

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

273 arguments to return the transformed position.

274

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

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

277 '''

278 if self.isunity:

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

280 else:

281 xyz1 = x, y, z, _1_0

282 s1 = self.s1

283 if inverse:

284 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

292 name=self.name)

293 if Vector_and_kwds:

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

295 if V:

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

297 return r

298

299

300 class Transforms(_NamedEnum):

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

302 to accommodate the L{_LazyNamedEnumItem} properties.

303 '''

304 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz):

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

306 '''

307 return Transform(**name_tx_ty_tz_s_sx_sy_sz)

308

309 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

312 Transforms._assert(

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

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

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

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

317

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

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

320

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

322

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

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

325

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

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

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

329

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

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

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

333

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

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

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

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

338 sz=_F( 0.156)),

339 Identity = _lazy(_Identity_),

340

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

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

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

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

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

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

347

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

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

350

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

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

353

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

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

356

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

358

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

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

361

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

363

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

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

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

367

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

369

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

371

372 WGS84 = _lazy(_WGS84_), # unity

373)

374

375

376 class Datum(_NamedEnumItem):

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

378 '''

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

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

381

382 def __init__(self, ellipsoid, transform=None, **name):

383 '''New L{Datum}.

384

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

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

387 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

388

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

390

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

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

393 not a L{Transform}.

394 '''

395 self._ellipsoid = ellipsoid or Datum._ellipsoid

396 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

397

398 self._transform = transform or Datum._transform

399 _xinstanceof(Transform, transform=self.transform)

400

401 self._register(Datums, _name__(name) or self.transform.name

402 or self.ellipsoid.name)

403

404 def __eq__(self, other):

405 '''Compare this and an other datum.

406

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

408

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

410 '''

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

412 self.ellipsoid == other.ellipsoid and

413 self.transform == other.transform)

414

415 def __hash__(self):

416 return self._hash # memoized

417

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

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

420

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

422 '''

423 _ = _xellipsoidall(point)

424 return point.toDatum(self)

425

426 def ecef(self, Ecef=None):

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

428

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

430

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

432

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

434

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

436 '''

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

438

439 @Property_RO

440 def ellipsoid(self):

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

442 '''

443 return self._ellipsoid

444

445 @Property_RO

446 def exactTM(self):

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

448 '''

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

450

451 @Property_RO

452 def _hash(self):

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

454

455 @property_RO

456 def isEllipsoidal(self):

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

458 '''

459 return self.ellipsoid.isEllipsoidal

460

461 @property_RO

462 def isOblate(self):

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

464 '''

465 return self.ellipsoid.isOblate

466

467 @property_RO

468 def isProlate(self):

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

470 '''

471 return self.ellipsoid.isProlate

472

473 @property_RO

474 def isSpherical(self):

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

476 '''

477 return self.ellipsoid.isSpherical

478

479 def toStr(self, sep=_COMMASPACE_, **name): # PYCHOK expected

480 '''Return this datum as a string.

481

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

483 @kwarg name: Optional, override C{B{name}=NN} (C{str}) or

484 C{None} to exclude this datum's name.

485

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

487 '''

488 name, _ = _name2__(**name) # name=None

489 t = [] if name is None else \

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

491 for a in (_ellipsoid_, _transform_):

492 v = getattr(self, a)

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

494 return sep.join(t)

495

496 @Property_RO

497 def transform(self):

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

499 '''

500 return self._transform

501

502

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

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

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

506

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

508 '''

509 if f is not None:

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

511 if raiser and not E:

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

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

514 return

515 elif isinstance(a_earth, Datum):

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

517 else:

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

519 if raiser and not E:

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

521 if D is None:

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

523 inst._datum = D

524

525

526 def _earth_ellipsoid(earth, **name_raiser):

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

528 '''

529 return Ellipsoids.Sphere if earth is R_M else (

530 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

531 _spherical_datum(earth, **name_raiser).ellipsoid)

532

533

534 def _ED2(radius, name):

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

536 '''

537 D = Datums.Sphere

538 E = D.ellipsoid

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

540 n = _under(_name__(name, _or_nameof=D))

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

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

543 return E, D

544

545

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

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

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

549

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

551 '''

552 if isinstance(earth, Datum):

553 D = earth

554 else:

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

556 if not E:

557 n = raiser or _earth_

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

559 if D is None:

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

561 if raiser and not D.isEllipsoidal:

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

563 return D

564

565

566 def _EnD3(earth, name):

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

568 '''

569 D, n = None, _under(_name__(name, _or_nameof=earth))

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

571 E = earth

572 elif isinstance(earth, Datum):

573 E = earth.ellipsoid

574 D = earth

575 elif _isRadius(earth):

576 E, D = _ED2(Radius_(earth), n)

577 n = E.name

578 elif isinstance(earth, a_f2Tuple):

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 E = Ellipsoid(a, f=f, name=n)

585 else:

586 n = E.name

587 D = Datums.Sphere

588 else:

589 E, n = None, NN

590 return E, n, D

591

592

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

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

595 '''

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

597

598

599 def _mean_radius(radius, *lats):

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

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

602 '''

603 if radius is R_M:

604 r = radius

605 elif _isRadius(radius):

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

607 else:

608 E = _ellipsoidal_datum(radius).ellipsoid

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

610 return r

611

612

613 def _negastr(name): # in .trf

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

615 '''

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

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

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

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

620 # def _negastr(name):

621 # n = name.translate(_MINUSxPLUS)

622 # ...

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

624

625

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

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

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

629

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

631 '''

632 if isinstance(earth, Datum):

633 D = earth

634 elif _isRadius(earth):

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

636 else:

637 D = _ellipsoidal_datum(earth, Error=Error, **name)

638 if raiser and not D.isSpherical:

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

640 return D

641

642

643 class Datums(_NamedEnum):

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

645 to accommodate the L{_LazyNamedEnumItem} properties.

646 '''

647 def _Lazy(self, ellipsoid_name, transform_name, **name):

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

649 '''

650 return Datum(Ellipsoids.get(ellipsoid_name),

651 Transforms.get(transform_name), **name)

652

653 Datums = Datums(Datum) # PYCHOK singleton

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

655# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

659 Datums._assert(

660 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

664

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

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

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

668

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

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

671

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

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

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

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

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

677

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

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

680

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

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

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

684

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

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

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

688

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

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

691

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

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

694

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

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

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

698

699 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

701

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

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

704

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

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

707

708 # XXX psuedo-ellipsoids for spherical LatLon

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

710

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

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

713

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

715 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

716

717 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

718)

719

720_WGS84 = Datums.WGS84

721 assert _WGS84.ellipsoid is _EWGS84

722# assert _WGS84.transform.isunity

723

724 if __name__ == '__main__':

725

726 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

727 from pygeodesy import printf

728

729 # __doc__ of this file, force all into registery

730 for r in (Datums, Transforms):

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

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

733

734# **) MIT License

735#

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

737#

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

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

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

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

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

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

744#

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

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

747#

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

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

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

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

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

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

754# OTHER DEALINGS IN THE SOFTWARE.