Coverage for pygeodesy / datums.py: 93%

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Datums and transformations thereof. 

5 

6Classes L{Datum} and L{Transform} and registries L{Datums} and L{Transforms}, respectively. 

7 

8Pure Python implementation of geodesy tools for ellipsoidal earth models, including datums 

9and ellipsoid parameters for different geographic coordinate systems and methods for 

10converting between them and to cartesian coordinates. Transcoded from JavaScript originals by 

11I{(C) Chris Veness 2005-2024} 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 

14Historical geodetic datums: a latitude/longitude point defines a geographic location on, above 

15or below the earth’s surface. Latitude is measured in degrees from the equator, lomgitude from 

16the International Reference Meridian and height in meters above an ellipsoid based on the given 

17datum. The datum in turn is based on a reference ellipsoid and tied to geodetic survey 

18reference points. 

19 

20Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems), but 

21previously various other reference ellipsoids and datum references were used. 

22 

23The UK Ordnance Survey National Grid References are still based on the otherwise historical OSGB36 

24datum, 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.Bessel1841: Datum(name='Bessel1841', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841) 

29@var Datums.DHDN: Datum(name='DHDN', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.DHDN) 

30@var Datums.ED50: Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50) 

31@var Datums.GDA2020: Datum(name='GDA2020', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84) 

32@var Datums.GRS80: Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84) 

33@var Datums.Irl1975: Datum(name='Irl1975', ellipsoid=Ellipsoids.AiryModified, transform=Transforms.Irl1975) 

34@var Datums.Krassovski1940: Datum(name='Krassovski1940', ellipsoid=Ellipsoids.Krassovski1940, transform=Transforms.Krassovski1940) 

35@var Datums.Krassowsky1940: Datum(name='Krassowsky1940', ellipsoid=Ellipsoids.Krassowsky1940, transform=Transforms.Krassowsky1940) 

36@var Datums.MGI: Datum(name='MGI', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.MGI) 

37@var Datums.NAD27: Datum(name='NAD27', ellipsoid=Ellipsoids.Clarke1866, transform=Transforms.NAD27) 

38@var Datums.NAD83: Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83) 

39@var Datums.NTF: Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF) 

40@var Datums.OSGB36: Datum(name='OSGB36', ellipsoid=Ellipsoids.Airy1830, transform=Transforms.OSGB36) 

41@var Datums.Potsdam: Datum(name='Potsdam', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841) 

42@var Datums.Sphere: Datum(name='Sphere', ellipsoid=Ellipsoids.Sphere, transform=Transforms.WGS84) 

43@var Datums.TokyoJapan: Datum(name='TokyoJapan', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.TokyoJapan) 

44@var Datums.WGS72: Datum(name='WGS72', ellipsoid=Ellipsoids.WGS72, transform=Transforms.WGS72) 

45@var Datums.WGS84: Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84) 

46 

47@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) 

48@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) 

49@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) 

50@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) 

51@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) 

52@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) 

53@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) 

54@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) 

55@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) 

56@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) 

57@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) 

58@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) 

59@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) 

60@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) 

61@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) 

62@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) 

63@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) 

64@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) 

65@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) 

66@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) 

67''' 

68# make sure int/int division yields float quotient, see .basics 

69from __future__ import division as _; del _ # noqa: E702 ; 

70 

71from pygeodesy.basics import _isin, islistuple, map2, neg, _xinstanceof, _zip 

72from pygeodesy.constants import R_M, _float as _F, _0_0, _1_0, _2_0, _8_0, _3600_0 

73# from pygeodesy.ecef import _4Ecef # _MODS 

74# from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase as _CEB, \ 

75# LatLonEllipsoidalBase as _LLEB # _MODS 

76from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, _EWGS84, \ 

77 Vector3Tuple 

78from pygeodesy.errors import _IsnotError, _TypeError, _xellipsoidall, _xkwds_pop2 

79# from pygeodesy.etm import ExactTransverseMercator # _MODS 

80from pygeodesy.fmath import _fdotf, fmean, Fmt, _operator 

81from pygeodesy.internals import _passarg, _under 

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

83 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DMAIN_,_DOT_, \ 

84 _earth_, _ED50_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \ 

85 _Krassovski1940_, _Krassowsky1940_, _MINUS_, _NAD27_, _NAD83_, \ 

86 _PLUS_, _s_, _Sphere_, _spherical_, _transform_, _Txyzsxyz7, \ 

87 _UNDER_, _WGS72_, _WGS84_ 

88from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

90 _NamedEnumItem 

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

92from pygeodesy.props import Property_RO, property_RO 

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

94from pygeodesy.units import _isRadius, Radius_, radians 

95# from pygeodesy.utily import sincos2_ # _MODS 

96 

97# from math import radians # from .units 

98# import operator as _operator # from .fmath 

99 

100__all__ = _ALL_LAZY.datums 

101__version__ = '26.07.21' 

102 

103_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

104_BD72_ = 'BD72' 

105_DHDN_ = 'DHDN' 

106_DHDNE_ = 'DHDNE' 

107_DHDNW_ = 'DHDNW' 

108_GDA2020_ = 'GDA2020' # in .trf 

109_Identity_ = 'Identity' 

110_Irl1965_ = 'Irl1965' 

111_Irl1975_ = 'Irl1975' 

112_MGI_ = 'MGI' 

113_NTF_ = 'NTF' 

114_OSGB36_ = 'OSGB36' 

115_Potsdam_ = 'Potsdam' 

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

117_SPR = _1_0 / _RPS # arc-seconds per radian 

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

119_TokyoJapan_ = 'TokyoJapan' 

120_uRad = _S1_S # PYCHOK micro-radians to radian 

121 

122 

123def _rps2(s_): # to C{radians} and C{arc-seconds}. 

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

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

126 return (_RPS * s_), s_ 

127 

128 

129def _spr2(r_): # to C{micro-radians} and C{micro-arc-seconds} 

130 return r_, (_SPR * r_) 

131 

132 

133class Transform(_NamedEnumItem): 

134 '''Helmert I{datum} transformation. 

135 

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

137 ''' 

138 _Txyzs7 = _Txyzsxyz7 

139 _Txyzs11 = _Txyzsxyz7[:3] + ('s1', 'rx', 'ry', 'rz') + _Txyzsxyz7[3:] 

140 

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

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

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

144 

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

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

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

148 

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

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

151 

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

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

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

155 

156 def __init__(self, name=NN, tx=0, ty=0, tz=0, # _Txyzsxyz7 order 

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

158 '''New L{Transform}. 

159 

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

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

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

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

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

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

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

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

168 @kwarg rx_ry_rz: Optional X, Y and Z rotation (C{micro-radians}), 

169 overriding C{sx}, C{sy} and C{sz}. 

170 

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

172 ''' 

173 if tx: 

174 self.tx = tx 

175 if ty: 

176 self.ty = ty 

177 if tz: 

178 self.tz = tz 

179 if s: 

180 self.s = s 

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

182 if sx: # secs to rads 

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

184 if sy: 

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

186 if sz: 

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

188 

189 self._register(Transforms, name) 

190 

191 def __eq__(self, other): 

192 '''Compare this and an other transform. 

193 

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

195 

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

197 ''' 

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

199 and _equall(other, self)) 

200 

201 def __hash__(self): 

202 return hash(tuple(self)) 

203 

204 def __iter__(self): 

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

206 ''' 

207 for n in self._Txyzs7: 

208 yield getattr(self, n) 

209 

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

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

212 

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

214 

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

216 

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

218 ''' 

219 _ = _xellipsoidall(point) 

220 return point.toTransform(self) 

221 

222 def __neg__(self): 

223 return self.inverse() 

224 

225 def inverse(self, **name): 

226 '''Return the inverse of this transform. 

227 

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

229 

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

231 ''' 

232 T = type(self)(**dict(self.items(inverse=True))) 

233 n = _name__(**name) or _negastr(self.name) 

234 if n: 

235 T.name = n # unregistered 

236 return T 

237 

238 @Property_RO 

239 def isunity(self): 

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

241 WGS84 with translation, scale and rotation all zero? 

242 ''' 

243 return not any(self) 

244 

245 def items(self, inverse=False): 

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

247 

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

249 ''' 

250 _p = neg if inverse else _passarg 

251 for n, x in _zip(self._Txyzs7, self): 

252 yield n, _p(x) 

253 

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

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

256 ''' 

257 Transform.isunity._update(self) 

258 self.s1 = s1 

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

260 return s 

261 

262 def toStr(self, prec=5, fmt=Fmt.g, **sep_name): # PYCHOK expected 

263 '''Return this transform as a string. 

264 

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

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

267 @kwarg sep_name: Optional C{B{name}=NN} (C{str}) or C{None} 

268 to exclude this transform's name and separater 

269 C{B{sep}=", "} to join the items (C{str}). 

270 

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

272 ''' 

273 return self._instr(*self._Txyzs11, fmt=fmt, prec=prec, **sep_name) 

274 

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

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

277 

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

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

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

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

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

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

284 arguments to return the transformed position. 

285 

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

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

288 ''' 

289 if self.isunity: 

290 pass # == inverse 

291 else: 

292 xyz1 = x, y, z, _1_0 

293 s1 = self.s1 

294 if inverse: 

295 xyz1 = map2(neg, xyz1) 

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

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

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

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

300 x = _fdotf(xyz1, s1, -self.rz, self.ry, self.tx) 

301 y = _fdotf(xyz1, self.rz, s1, -self.rx, self.ty) 

302 z = _fdotf(xyz1, -self.ry, self.rx, s1, self.tz) 

303 

304 return self._V(x, y, z, **Vector_and_kwds) 

305 

306 def _V(self, x, y, z, Vector=None, **kwds): 

307 '''(INTERNAL) Return C{r} as a C{Vector}. 

308 ''' 

309 n, kwds = _xkwds_pop2(kwds, name=self.name) 

310 r = Vector3Tuple(x, y, z, name=n) 

311 if Vector: 

312 r = Vector(r, name=n, **kwds) 

313 return r 

314 

315 

316class Similarity(Transform): # in .PyRDNAP 

317 '''Similarity transformation. 

318 ''' 

319 _Txyzs7 = \ 

320 _Txyzs11 = _Txyzsxyz7[:4] + ('rx', 'ry', 'rz') 

321 

322 def __init__(self, name=NN, tx=0, ty=0, tz=0, # _Txyzsxyz7 order 

323 s=0, rx=0, ry=0, rz=0): 

324 '''New L{Similarity}. 

325 

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

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

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

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

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

331 @kwarg rx: X rotation (C{micro-radians}). 

332 @kwarg ry: Y rotation (C{micro-radians}). 

333 @kwarg rz: Z rotation (C{micro-radians}). 

334 

335 @raise NameError: Similarity with that B{C{name}} already exists. 

336 ''' 

337 Transform.__init__(self, name, tx, ty, tz, s) # _Txyzsxyz7 order 

338 

339 if rx: 

340 self.rx, self.sx = _spr2(rx) 

341 if ry: 

342 self.ry, self.sy = _spr2(ry) 

343 if rz: 

344 self.rz, self.sz = _spr2(rz) 

345 

346 @Property_RO 

347 def _sForward(self): 

348 '''(INTERNAL) Get the forward 3-D similarity transform, [3x4] matrix. 

349 ''' 

350 sx, cx, sy, cy, sz, cz = _MODS.utily.sincos2_(self.rx * _uRad, 

351 self.ry * _uRad, 

352 self.rz * _uRad) 

353 czsy = cz * sy 

354 szsy = sz * sy 

355 return (cz * cy, sz * cx + czsy * sx, sz * sx - czsy * cx, self.tx, 

356 -sz * cy, cz * cx - szsy * sx, cz * sx + szsy * cx, self.ty, 

357 sy, -cy * sx, cy * cx, self.tz) 

358 

359 @Property_RO 

360 def _sInverse(self): 

361 '''(INTERNAL) Get the inverse 3-D similarity transform, [3x4] matrix. 

362 ''' 

363 return self.inverse()._sForward 

364 

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

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

367 ''' 

368 Similarity._sForward._update(self) 

369 Similarity._sInverse._update(self) 

370 return Transform._s_s1(self, s1) 

371 

372 def transform(self, x, y, z, xyz0=(), inverse=False, **Vector_and_kwds): # PYCHOK signature 

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

374 

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

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

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

378 @arg xyz0: Optional pivot point (3-tuple C{meter}). 

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

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

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

382 arguments to return the transformed position. 

383 

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

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

386 ''' 

387 if self.isunity: 

388 pass # == inverse 

389 else: 

390 if xyz0: 

391 x0, y0, z0 = xyz0 

392 x -= x0 

393 y -= y0 

394 z -= z0 

395 else: 

396 x0 = y0 = z0 = _0_0 

397 s1 = self.s1 

398 _1xyz1 = _1_0, (x * s1), (y * s1), (z * s1), _1_0 

399 

400 S = self._sInverse if inverse else self._sForward 

401 x = _fdotf(_1xyz1, x0, *S[0:4]) 

402 y = _fdotf(_1xyz1, y0, *S[4:8]) 

403 z = _fdotf(_1xyz1, z0, *S[8:12]) 

404 

405 return self._V(x, y, z, **Vector_and_kwds) 

406 

407 

408class Transforms(_NamedEnum): 

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

410 to accommodate the L{_LazyNamedEnumItem} properties. 

411 ''' 

412 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz): 

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

414 ''' 

415 return Transform(**name_tx_ty_tz_s_sx_sy_sz) 

416 

417Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

420Transforms._assert( 

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

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

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

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

425 

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

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

428 

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

430 

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

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

433 

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

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

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

437 

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

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

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

441 

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

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

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

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

446 sz=_F( 0.156)), 

447 

448 Identity = _lazy(_Identity_), 

449 

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

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

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

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

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

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

456 

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

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

459 

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

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

462 

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

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

465 

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

467 

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

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

470 

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

472 

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

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

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

476 

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

478 

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

480 

481 WGS84 = _lazy(_WGS84_), # unity 

482) 

483 

484 

485class Datum(_NamedEnumItem): 

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

487 ''' 

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

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

490 

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

492 '''New L{Datum}. 

493 

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

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

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

497 

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

499 

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

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

502 not a L{Transform}. 

503 ''' 

504 self._ellipsoid = ellipsoid or Datum._ellipsoid 

505 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

506 

507 self._transform = transform or Datum._transform 

508 _xinstanceof(Transform, transform=self.transform) 

509 

510 self._register(Datums, _name__(name) or self.transform.name # first 

511 or self.ellipsoid.name) 

512 

513 def __eq__(self, other): 

514 '''Compare this and an other datum. 

515 

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

517 

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

519 ''' 

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

521 self.ellipsoid == other.ellipsoid and 

522 self.transform == other.transform) 

523 

524 def __hash__(self): 

525 return self._hash # memoized 

526 

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

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

529 

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

531 ''' 

532 _ = _xellipsoidall(point) 

533 return point.toDatum(self) 

534 

535 def ecef(self, Ecef=None): 

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

537 

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

539 

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

541 

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

543 

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

545 ''' 

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

547 

548 @Property_RO 

549 def ellipsoid(self): 

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

551 ''' 

552 return self._ellipsoid 

553 

554 @Property_RO 

555 def exactTM(self): 

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

557 ''' 

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

559 

560 @Property_RO 

561 def _hash(self): 

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

563 

564 @property_RO 

565 def isEllipsoidal(self): 

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

567 ''' 

568 return self.ellipsoid.isEllipsoidal 

569 

570 @property_RO 

571 def isOblate(self): 

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

573 ''' 

574 return self.ellipsoid.isOblate 

575 

576 @property_RO 

577 def isProlate(self): 

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

579 ''' 

580 return self.ellipsoid.isProlate 

581 

582 @property_RO 

583 def isSpherical(self): 

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

585 ''' 

586 return self.ellipsoid.isSpherical 

587 

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

589 '''Return this datum as a string. 

590 

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

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

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

594 

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

596 ''' 

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

598 t = [] if name is None else \ 

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

600 for a in (_ellipsoid_, _transform_): 

601 v = getattr(self, a) 

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

603 return sep.join(t) 

604 

605 @Property_RO 

606 def transform(self): 

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

608 ''' 

609 return self._transform 

610 

611 

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

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

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

615 

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

617 ''' 

618 if f is not None: 

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

620 if raiser and not E: 

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

622 elif _isin(a_earth, None, _EWGS84, _WGS84) and inst._datum is _WGS84: 

623 return 

624 elif isinstance(a_earth, Datum): 

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

626 else: 

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

628 if raiser and not E: 

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

630 if D is None: 

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

632 inst._datum = D 

633 

634 

635def _earth_ellipsoid(earth, **name_raiser): 

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

637 ''' 

638 return Ellipsoids.Sphere if earth is R_M else ( 

639 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

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

641 

642 

643def _ED2(radius, name): 

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

645 ''' 

646 D = Datums.Sphere 

647 E = D.ellipsoid 

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

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

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

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

652 return E, D 

653 

654 

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

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

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

658 

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

660 ''' 

661 if isinstance(earth, Datum): 

662 D = earth 

663 else: 

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

665 if not E: 

666 n = raiser or _earth_ 

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

668 if D is None: 

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

670 if raiser and not D.isEllipsoidal: 

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

672 return D 

673 

674 

675def _EnD3(earth, name): 

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

677 ''' 

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

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

680 E = earth 

681 elif isinstance(earth, Datum): 

682 E = earth.ellipsoid 

683 D = earth 

684 elif _isRadius(earth): 

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

686 n = E.name 

687 elif isinstance(earth, a_f2Tuple): 

688 E = earth.ellipsoid(name=n) 

689 elif islistuple(earth, minum=2): 

690 E = Ellipsoids.Sphere 

691 a, f = earth[:2] 

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

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

694 else: 

695 n = E.name 

696 D = Datums.Sphere 

697 else: 

698 E, n = None, NN 

699 return E, n, D 

700 

701 

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

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

704 ''' 

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

706 

707 

708def _mean_radius(radius, *lats): 

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

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

711 ''' 

712 if radius is R_M: 

713 r = radius 

714 elif _isRadius(radius): 

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

716 else: 

717 E = _ellipsoidal_datum(radius).ellipsoid 

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

719 return r 

720 

721 

722def _negastr(name): # in .trf, test/testTrf 

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

724 ''' 

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

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

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

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

729 # def _negastr(name): 

730 # n = name.translate(_MINUSxPLUS) 

731 # ... 

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

733 

734 

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

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

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

738 

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

740 ''' 

741 if isinstance(earth, Datum): 

742 D = earth 

743 elif _isRadius(earth): 

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

745 else: 

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

747 if raiser and not D.isSpherical: 

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

749 return D 

750 

751 

752class Datums(_NamedEnum): 

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

754 to accommodate the L{_LazyNamedEnumItem} properties. 

755 ''' 

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

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

758 ''' 

759 return Datum(Ellipsoids.get(ellipsoid_name), 

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

761 

762Datums = Datums(Datum) # PYCHOK singleton 

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

764# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

768Datums._assert( 

769 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

773 

774 # Netherlands' RD-NAP RijksDriehoeksmeting-NormaalAmsterdamsPeil, ETRS89 

775 Bessel1841 = _lazy(_Bessel1841_, _Bessel1841_, _Bessel1841_), 

776 

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

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

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

780 

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

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

783 

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

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

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

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

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

789 

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

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