Coverage for pygeodesy/datums.py: 94%

253 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-15 16:36 -0400

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

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

68from __future__ import division as _; del _ # PYCHOK semicolon 

69 

70from pygeodesy.basics import islistuple, map2, neg, _xinstanceof, _zip 

71from 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 

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

75 Vector3Tuple 

76from pygeodesy.errors import _IsnotError, _TypeError, _xattr, _xellipsoidall, _xkwds, \ 

77 _xkwds_pop2 

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

79from pygeodesy.internals import _passarg, _under 

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

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

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

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

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

85 _WGS72_, _WGS84_ 

86from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

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

89from pygeodesy.props import Property_RO, property_RO 

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

91from pygeodesy.units import _isRadius, Radius_, radians 

92 

93# from math import radians # from .units 

94# import operator as _operator # from .fmath 

95 

96__all__ = _ALL_LAZY.datums 

97__version__ = '24.05.14' 

98 

99_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

100_BD72_ = 'BD72' 

101_DHDN_ = 'DHDN' 

102_DHDNE_ = 'DHDNE' 

103_DHDNW_ = 'DHDNW' 

104_ED50_ = 'ED50' 

105_GDA2020_ = 'GDA2020' # in .trf 

106_Identity_ = 'Identity' 

107_Irl1965_ = 'Irl1965' 

108_Irl1975_ = 'Irl1975' 

109_MGI_ = 'MGI' 

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

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

112_NTF_ = 'NTF' 

113_OSGB36_ = 'OSGB36' 

114_Potsdam_ = 'Potsdam' 

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

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

117_TokyoJapan_ = 'TokyoJapan' 

118 

119 

120class Transform(_NamedEnumItem): 

121 '''Helmert I{datum} transformation. 

122 

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

124 ''' 

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

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

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

128 

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

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

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

132 

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

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

135 

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

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

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

139 

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

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

142 '''New L{Transform}. 

143 

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

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

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

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

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

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

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

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

152 

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

154 ''' 

155 if tx: 

156 self.tx = tx 

157 if ty: 

158 self.ty = ty 

159 if tz: 

160 self.tz = tz 

161 if s: 

162 self.s = s 

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

164 if sx: # secs to rads 

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

166 if sy: 

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

168 if sz: 

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

170 

171 self._register(Transforms, name) 

172 

173 def __eq__(self, other): 

174 '''Compare this and an other transform. 

175 

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

177 

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

179 ''' 

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

181 and _equall(other, self)) 

182 

183 def __hash__(self): 

184 return hash(tuple(self)) 

185 

186 def __iter__(self): 

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

188 ''' 

189 for n in _Names7: 

190 yield getattr(self, n) 

191 

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

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

194 

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

196 

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

198 

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

200 ''' 

201 _ = _xellipsoidall(point) 

202 return point.toTransform(self) 

203 

204 def __neg__(self): 

205 return self.inverse() 

206 

207 def inverse(self, name=NN): 

208 '''Return the inverse of this transform. 

209 

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

211 

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

213 ''' 

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

215 n = name or _negastr(self.name) 

216 if n: 

217 r.name = n # unregistered 

218 return r 

219 

220 @Property_RO 

221 def isunity(self): 

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

223 WGS84 with translation, scale and rotation all zero? 

224 ''' 

225 return not any(self) 

226 

227 def items(self, inverse=False): 

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

229 

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

231 ''' 

232 _p = neg if inverse else _passarg 

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

234 yield n, _p(x) 

235 

236 def _rps2(self, s_): 

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

238 ''' 

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

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

241 return (_RPS * s_), s_ 

242 

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

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

245 ''' 

246 Transform.isunity._update(self) 

247 self.s1 = s1 

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

249 return s 

250 

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

252 '''Return this transform as a string. 

253 

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

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

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

257 this transform's name. 

258 

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

260 ''' 

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

262 

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

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

265 

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

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

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

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

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

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

272 arguments to return the transformed position. 

273 

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

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

276 ''' 

277 if self.isunity: 

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

279 else: 

280 xyz1 = x, y, z, _1_0 

281 s1 = self.s1 

282 if inverse: 

283 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

291 name=self.name) 

292 if Vector_and_kwds: 

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

294 if V: 

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

296 return r 

297 

298 

299class Transforms(_NamedEnum): 

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

301 to accommodate the L{_LazyNamedEnumItem} properties. 

302 ''' 

303 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz): 

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

305 ''' 

306 return Transform(**name_tx_ty_tz_s_sx_sy_sz) 

307 

308Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

311Transforms._assert( 

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

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

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

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

316 

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

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

319 

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

321 

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

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

324 

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

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

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

328 

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

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

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

332 

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

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

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

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

337 sz=_F( 0.156)), 

338 Identity = _lazy(_Identity_), 

339 

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

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

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

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

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

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

346 

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

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

349 

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

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

352 

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

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

355 

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

357 

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

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

360 

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

362 

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

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

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

366 

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

368 

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

370 

371 WGS84 = _lazy(_WGS84_), # unity 

372) 

373 

374 

375class Datum(_NamedEnumItem): 

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

377 ''' 

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

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

380 

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

382 '''New L{Datum}. 

383 

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

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

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

387 

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

389 

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

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

392 not a L{Transform}. 

393 ''' 

394 self._ellipsoid = ellipsoid or Datum._ellipsoid 

395 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

396 

397 self._transform = transform or Datum._transform 

398 _xinstanceof(Transform, transform=self.transform) 

399 

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

401 

402 def __eq__(self, other): 

403 '''Compare this and an other datum. 

404 

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

406 

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

408 ''' 

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

410 self.ellipsoid == other.ellipsoid and 

411 self.transform == other.transform) 

412 

413 def __hash__(self): 

414 return self._hash # memoized 

415 

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

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

418 

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

420 ''' 

421 _ = _xellipsoidall(point) 

422 return point.toDatum(self) 

423 

424 def ecef(self, Ecef=None): 

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

426 

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

428 

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

430 

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

432 

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

434 ''' 

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

436 

437 @Property_RO 

438 def ellipsoid(self): 

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

440 ''' 

441 return self._ellipsoid 

442 

443 @Property_RO 

444 def exactTM(self): 

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

446 ''' 

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

448 

449 @Property_RO 

450 def _hash(self): 

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

452 

453 @property_RO 

454 def isEllipsoidal(self): 

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

456 ''' 

457 return self.ellipsoid.isEllipsoidal 

458 

459 @property_RO 

460 def isOblate(self): 

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

462 ''' 

463 return self.ellipsoid.isOblate 

464 

465 @property_RO 

466 def isProlate(self): 

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

468 ''' 

469 return self.ellipsoid.isProlate 

470 

471 @property_RO 

472 def isSpherical(self): 

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

474 ''' 

475 return self.ellipsoid.isSpherical 

476 

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

478 '''Return this datum as a string. 

479 

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

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

482 this datum's name. 

483 

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

485 ''' 

486 t = [] if name is None else \ 

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

488 for a in (_ellipsoid_, _transform_): 

489 v = getattr(self, a) 

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

491 return sep.join(t) 

492 

493 @Property_RO 

494 def transform(self): 

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

496 ''' 

497 return self._transform 

498 

499 

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

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

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

503 

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

505 ''' 

506 if f is not None: 

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

508 if raiser and not E: 

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

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

511 return 

512 elif isinstance(a_earth, Datum): 

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

514 else: 

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

516 if raiser and not E: 

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

518 if D is None: 

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

520 inst._datum = D 

521 

522 

523def _earth_ellipsoid(earth, *name_raiser): 

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

525 ''' 

526 return Ellipsoids.Sphere if earth is R_M else ( 

527 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

528 _spherical_datum(earth, *name_raiser).ellipsoid) 

529 

530 

531def _ED2(radius, name): 

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

533 ''' 

534 D = Datums.Sphere 

535 E = D.ellipsoid 

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

537 n = _under(name) 

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

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

540 return E, D 

541 

542 

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

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

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

546 

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

548 ''' 

549 if isinstance(earth, Datum): 

550 D = earth 

551 else: 

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

553 if not E: 

554 n = raiser or _earth_ 

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

556 if D is None: 

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

558 if raiser and not D.isEllipsoidal: 

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

560 return D 

561 

562 

563def _EnD3(earth, name): 

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

565 ''' 

566 D = None 

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

568 E = earth 

569 n = _under(name or E.name) 

570 elif isinstance(earth, Datum): 

571 E = earth.ellipsoid 

572 n = _under(name or earth.name) 

573 D = earth 

574 elif _isRadius(earth): 

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

576 n = E.name 

577 elif isinstance(earth, a_f2Tuple): 

578 n = _under(name or earth.name) 

579 E = earth.ellipsoid(name=n) 

580 elif islistuple(earth, minum=2): 

581 E = Ellipsoids.Sphere 

582 a, f = earth[:2] 

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

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

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

586 else: 

587 n = E.name 

588 D = Datums.Sphere 

589 else: 

590 E, n = None, NN 

591 return E, n, D 

592 

593 

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

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

596 ''' 

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

598 

599 

600def _mean_radius(radius, *lats): 

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

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

603 ''' 

604 if radius is R_M: 

605 r = radius 

606 elif _isRadius(radius): 

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

608 else: 

609 E = _ellipsoidal_datum(radius).ellipsoid 

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

611 return r 

612 

613 

614def _negastr(name): # in .trf 

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

616 ''' 

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

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

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

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

621 # def _negastr(name): 

622 # n = name.translate(_MINUSxPLUS) 

623 # ... 

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

625 

626 

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

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

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

630 

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

632 ''' 

633 if _isRadius(earth): 

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

635 else: 

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

637 if raiser and not D.isSpherical: 

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

639 return D 

640 

641 

642class Datums(_NamedEnum): 

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

644 to accommodate the L{_LazyNamedEnumItem} properties. 

645 ''' 

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

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

648 ''' 

649 return Datum(Ellipsoids.get(ellipsoid_name), 

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

651 

652Datums = Datums(Datum) # PYCHOK singleton 

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

654# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

658Datums._assert( 

659 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

663 

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

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

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

667 

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

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

670 

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

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

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

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

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

676 

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

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

679 

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

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

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

683 

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

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

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

687 

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

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

690 

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

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

693 

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

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

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

697 

698 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

700 

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

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

703 

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

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

706 

707 # XXX psuedo-ellipsoids for spherical LatLon 

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

709 

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

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

712 

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

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

715 

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

717) 

718 

719_WGS84 = Datums.WGS84 

720assert _WGS84.ellipsoid is _EWGS84 

721# assert _WGS84.transform.isunity 

722 

723if __name__ == '__main__': 

724 

725 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

726 from pygeodesy import printf 

727 

728 # __doc__ of this file, force all into registery 

729 for r in (Datums, Transforms): 

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

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

732 

733# **) MIT License 

734# 

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

736# 

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

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

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

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

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

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

743# 

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

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

746# 

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

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

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

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

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

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

753# OTHER DEALINGS IN THE SOFTWARE.