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_, _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.05.04' 

102 

103_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

104_BD72_ = 'BD72' 

105_DHDN_ = 'DHDN' 

106_DHDNE_ = 'DHDNE' 

107_DHDNW_ = 'DHDNW' 

108_ED50_ = 'ED50' 

109_GDA2020_ = 'GDA2020' # in .trf 

110_Identity_ = 'Identity' 

111_Irl1965_ = 'Irl1965' 

112_Irl1975_ = 'Irl1975' 

113_MGI_ = 'MGI' 

114_NTF_ = 'NTF' 

115_OSGB36_ = 'OSGB36' 

116_Potsdam_ = 'Potsdam' 

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

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

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

120_TokyoJapan_ = 'TokyoJapan' 

121_uRad = _S1_S # PYCHOK micro-radians to radian 

122 

123 

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

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

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

127 return (_RPS * s_), s_ 

128 

129 

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

131 return r_, (_SPR * r_) 

132 

133 

134class Transform(_NamedEnumItem): 

135 '''Helmert I{datum} transformation. 

136 

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

138 ''' 

139 _Txyzs7 = _Txyzsxyz7 

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

141 

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

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

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

145 

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

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

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

149 

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

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

152 

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

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

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

156 

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

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

159 '''New L{Transform}. 

160 

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

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

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

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

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

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

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

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

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

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

171 

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

173 ''' 

174 if tx: 

175 self.tx = tx 

176 if ty: 

177 self.ty = ty 

178 if tz: 

179 self.tz = tz 

180 if s: 

181 self.s = s 

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

183 if sx: # secs to rads 

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

185 if sy: 

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

187 if sz: 

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

189 

190 self._register(Transforms, name) 

191 

192 def __eq__(self, other): 

193 '''Compare this and an other transform. 

194 

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

196 

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

198 ''' 

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

200 and _equall(other, self)) 

201 

202 def __hash__(self): 

203 return hash(tuple(self)) 

204 

205 def __iter__(self): 

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

207 ''' 

208 for n in self._Txyzs7: 

209 yield getattr(self, n) 

210 

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

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

213 

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

215 

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

217 

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

219 ''' 

220 _ = _xellipsoidall(point) 

221 return point.toTransform(self) 

222 

223 def __neg__(self): 

224 return self.inverse() 

225 

226 def inverse(self, **name): 

227 '''Return the inverse of this transform. 

228 

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

230 

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

232 ''' 

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

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

235 if n: 

236 T.name = n # unregistered 

237 return T 

238 

239 @Property_RO 

240 def isunity(self): 

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

242 WGS84 with translation, scale and rotation all zero? 

243 ''' 

244 return not any(self) 

245 

246 def items(self, inverse=False): 

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

248 

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

250 ''' 

251 _p = neg if inverse else _passarg 

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

253 yield n, _p(x) 

254 

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

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

257 ''' 

258 Transform.isunity._update(self) 

259 self.s1 = s1 

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

261 return s 

262 

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

264 '''Return this transform as a string. 

265 

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

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

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

269 to exclude this transform's name and separater 

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

271 

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

273 ''' 

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

275 

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

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

278 

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

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

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

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

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

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

285 arguments to return the transformed position. 

286 

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

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

289 ''' 

290 if self.isunity: 

291 pass # == inverse 

292 else: 

293 xyz1 = x, y, z, _1_0 

294 s1 = self.s1 

295 if inverse: 

296 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

304 

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

306 

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

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

309 ''' 

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

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

312 if Vector: 

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

314 return r 

315 

316 

317class Similarity(Transform): # in .PyRDNAP 

318 '''Similarity transformation. 

319 ''' 

320 _Txyzs7 = \ 

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

322 

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

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

325 '''New L{Similarity}. 

326 

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

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

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

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

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

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

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

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

335 

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

337 ''' 

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

339 

340 if rx: 

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

342 if ry: 

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

344 if rz: 

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

346 

347 @Property_RO 

348 def _sForward(self): 

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

350 ''' 

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

352 self.ry * _uRad, 

353 self.rz * _uRad) 

354 czsy = cz * sy 

355 szsy = sz * sy 

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

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

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

359 

360 @Property_RO 

361 def _sInverse(self): 

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

363 ''' 

364 return self.inverse()._sForward 

365 

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

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

368 ''' 

369 Similarity._sForward._update(self) 

370 Similarity._sInverse._update(self) 

371 return Transform._s_s1(self, s1) 

372 

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

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

375 

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

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

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

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

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

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

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

383 arguments to return the transformed position. 

384 

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

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

387 ''' 

388 if self.isunity: 

389 pass # == inverse 

390 else: 

391 if xyz0: 

392 x0, y0, z0 = xyz0 

393 x -= x0 

394 y -= y0 

395 z -= z0 

396 else: 

397 x0 = y0 = z0 = _0_0 

398 s1 = self.s1 

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

400 

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

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

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

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

405 

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

407 

408 

409class Transforms(_NamedEnum): 

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

411 to accommodate the L{_LazyNamedEnumItem} properties. 

412 ''' 

413 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz): 

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

415 ''' 

416 return Transform(**name_tx_ty_tz_s_sx_sy_sz) 

417 

418Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

421Transforms._assert( 

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

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

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

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

426 

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

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

429 

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

431 

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

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

434 

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

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

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

438 

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

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

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

442 

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

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

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

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

447 sz=_F( 0.156)), 

448 

449 Identity = _lazy(_Identity_), 

450 

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

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

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

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

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

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

457 

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

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

460 

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

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

463 

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

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

466 

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

468 

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

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

471 

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

473 

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

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

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

477 

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

479 

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

481 

482 WGS84 = _lazy(_WGS84_), # unity 

483) 

484 

485 

486class Datum(_NamedEnumItem): 

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

488 ''' 

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

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

491 

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

493 '''New L{Datum}. 

494 

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

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

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

498 

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

500 

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

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

503 not a L{Transform}. 

504 ''' 

505 self._ellipsoid = ellipsoid or Datum._ellipsoid 

506 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

507 

508 self._transform = transform or Datum._transform 

509 _xinstanceof(Transform, transform=self.transform) 

510 

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

512 or self.ellipsoid.name) 

513 

514 def __eq__(self, other): 

515 '''Compare this and an other datum. 

516 

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

518 

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

520 ''' 

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

522 self.ellipsoid == other.ellipsoid and 

523 self.transform == other.transform) 

524 

525 def __hash__(self): 

526 return self._hash # memoized 

527 

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

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

530 

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

532 ''' 

533 _ = _xellipsoidall(point) 

534 return point.toDatum(self) 

535 

536 def ecef(self, Ecef=None): 

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

538 

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

540 

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

542 

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

544 

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

546 ''' 

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

548 

549 @Property_RO 

550 def ellipsoid(self): 

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

552 ''' 

553 return self._ellipsoid 

554 

555 @Property_RO 

556 def exactTM(self): 

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

558 ''' 

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

560 

561 @Property_RO 

562 def _hash(self): 

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

564 

565 @property_RO 

566 def isEllipsoidal(self): 

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

568 ''' 

569 return self.ellipsoid.isEllipsoidal 

570 

571 @property_RO 

572 def isOblate(self): 

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

574 ''' 

575 return self.ellipsoid.isOblate 

576 

577 @property_RO 

578 def isProlate(self): 

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

580 ''' 

581 return self.ellipsoid.isProlate 

582 

583 @property_RO 

584 def isSpherical(self): 

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

586 ''' 

587 return self.ellipsoid.isSpherical 

588 

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

590 '''Return this datum as a string. 

591 

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

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

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

595 

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

597 ''' 

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

599 t = [] if name is None else \ 

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

601 for a in (_ellipsoid_, _transform_): 

602 v = getattr(self, a) 

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

604 return sep.join(t) 

605 

606 @Property_RO 

607 def transform(self): 

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

609 ''' 

610 return self._transform 

611 

612 

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

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

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

616 

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

618 ''' 

619 if f is not None: 

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

621 if raiser and not E: 

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

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

624 return 

625 elif isinstance(a_earth, Datum): 

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

627 else: 

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

629 if raiser and not E: 

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

631 if D is None: 

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

633 inst._datum = D 

634 

635 

636def _earth_ellipsoid(earth, **name_raiser): 

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

638 ''' 

639 return Ellipsoids.Sphere if earth is R_M else ( 

640 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

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

642 

643 

644def _ED2(radius, name): 

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

646 ''' 

647 D = Datums.Sphere 

648 E = D.ellipsoid 

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

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

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

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

653 return E, D 

654 

655 

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

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

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

659 

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

661 ''' 

662 if isinstance(earth, Datum): 

663 D = earth 

664 else: 

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

666 if not E: 

667 n = raiser or _earth_ 

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

669 if D is None: 

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

671 if raiser and not D.isEllipsoidal: 

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

673 return D 

674 

675 

676def _EnD3(earth, name): 

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

678 ''' 

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

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

681 E = earth 

682 elif isinstance(earth, Datum): 

683 E = earth.ellipsoid 

684 D = earth 

685 elif _isRadius(earth): 

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

687 n = E.name 

688 elif isinstance(earth, a_f2Tuple): 

689 E = earth.ellipsoid(name=n) 

690 elif islistuple(earth, minum=2): 

691 E = Ellipsoids.Sphere 

692 a, f = earth[:2] 

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

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

695 else: 

696 n = E.name 

697 D = Datums.Sphere 

698 else: 

699 E, n = None, NN 

700 return E, n, D 

701 

702 

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

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

705 ''' 

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

707 

708 

709def _mean_radius(radius, *lats): 

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

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

712 ''' 

713 if radius is R_M: 

714 r = radius 

715 elif _isRadius(radius): 

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

717 else: 

718 E = _ellipsoidal_datum(radius).ellipsoid 

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

720 return r 

721 

722 

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

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

725 ''' 

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

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

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

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

730 # def _negastr(name): 

731 # n = name.translate(_MINUSxPLUS) 

732 # ... 

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

734 

735 

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

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

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

739 

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

741 ''' 

742 if isinstance(earth, Datum): 

743 D = earth 

744 elif _isRadius(earth): 

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

746 else: 

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

748 if raiser and not D.isSpherical: 

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

750 return D 

751 

752 

753class Datums(_NamedEnum): 

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

755 to accommodate the L{_LazyNamedEnumItem} properties. 

756 ''' 

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

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

759 ''' 

760 return Datum(Ellipsoids.get(ellipsoid_name), 

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

762 

763Datums = Datums(Datum) # PYCHOK singleton 

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

765# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

769Datums._assert( 

770 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

774 

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

776 Bessel1841 = _lazy(_Bessel1841_, _Bessel1841_, _Bessel1841_), 

777 

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

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

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

781 

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

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

784 

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

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

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

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

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

790 

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