Coverage for pygeodesy/datums.py: 94%

252 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-06 16:50 -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, _passarg, _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.interns import NN, _a_, _Airy1830_, _AiryModified_, _BAR_, _Bessel1841_, \ 

80 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \ 

81 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \ 

82 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \ 

83 _NAD83_, _PLUS_, _s_, _Sphere_, _spherical_, _transform_, \ 

84 _UNDER_, _WGS72_, _WGS84_, _under 

85from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

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

88from pygeodesy.props import Property_RO, property_RO 

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

90from pygeodesy.units import _isRadius, Radius_, radians 

91 

92# from math import radians # from .units 

93# import operator as _operator # from .fmath 

94 

95__all__ = _ALL_LAZY.datums 

96__version__ = '24.03.07' 

97 

98_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

99_BD72_ = 'BD72' 

100_DHDN_ = 'DHDN' 

101_DHDNE_ = 'DHDNE' 

102_DHDNW_ = 'DHDNW' 

103_ED50_ = 'ED50' 

104_GDA2020_ = 'GDA2020' # in .trf 

105_Identity_ = 'Identity' 

106_Irl1965_ = 'Irl1965' 

107_Irl1975_ = 'Irl1975' 

108_MGI_ = 'MGI' 

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

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

111_NTF_ = 'NTF' 

112_OSGB36_ = 'OSGB36' 

113_Potsdam_ = 'Potsdam' 

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

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

116_TokyoJapan_ = 'TokyoJapan' 

117 

118 

119class Transform(_NamedEnumItem): 

120 '''Helmert I{datum} transformation. 

121 

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

123 ''' 

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

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

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

127 

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

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

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

131 

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

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

134 

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

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

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

138 

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

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

141 '''New L{Transform}. 

142 

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

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

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

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

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

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

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

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

151 

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

153 ''' 

154 if tx: 

155 self.tx = tx 

156 if ty: 

157 self.ty = ty 

158 if tz: 

159 self.tz = tz 

160 if s: 

161 self.s = s 

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

163 if sx: # secs to rads 

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

165 if sy: 

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

167 if sz: 

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

169 

170 self._register(Transforms, name) 

171 

172 def __eq__(self, other): 

173 '''Compare this and an other transform. 

174 

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

176 

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

178 ''' 

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

180 and _equall(other, self)) 

181 

182 def __hash__(self): 

183 return hash(tuple(self)) 

184 

185 def __iter__(self): 

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

187 ''' 

188 for n in _Names7: 

189 yield getattr(self, n) 

190 

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

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

193 

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

195 

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

197 

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

199 ''' 

200 _ = _xellipsoidall(point) 

201 return point.toTransform(self) 

202 

203 def __neg__(self): 

204 return self.inverse() 

205 

206 def inverse(self, name=NN): 

207 '''Return the inverse of this transform. 

208 

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

210 

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

212 ''' 

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

214 n = name or _negastr(self.name) 

215 if n: 

216 r.name = n # unregistered 

217 return r 

218 

219 @Property_RO 

220 def isunity(self): 

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

222 WGS84 with translation, scale and rotation all zero? 

223 ''' 

224 return not any(self) 

225 

226 def items(self, inverse=False): 

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

228 

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

230 ''' 

231 _p = neg if inverse else _passarg 

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

233 yield n, _p(x) 

234 

235 def _rps2(self, s_): 

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

237 ''' 

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

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

240 return (_RPS * s_), s_ 

241 

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

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

244 ''' 

245 Transform.isunity._update(self) 

246 self.s1 = s1 

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

248 return s 

249 

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

251 '''Return this transform as a string. 

252 

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

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

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

256 this transform's name. 

257 

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

259 ''' 

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

261 

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

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

264 

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

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

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

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

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

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

271 arguments to return the transformed position. 

272 

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

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

275 ''' 

276 if self.isunity: 

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

278 else: 

279 xyz1 = x, y, z, _1_0 

280 s1 = self.s1 

281 if inverse: 

282 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

290 name=self.name) 

291 if Vector_and_kwds: 

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

293 if V: 

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

295 return r 

296 

297 

298class Transforms(_NamedEnum): 

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

300 to accommodate the L{_LazyNamedEnumItem} properties. 

301 ''' 

302 def _Lazy(self, **name_tx_ty_tz_s_sx_sy_sz): 

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

304 ''' 

305 return Transform(**name_tx_ty_tz_s_sx_sy_sz) 

306 

307Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

310Transforms._assert( 

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

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

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

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

315 

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

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

318 

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

320 

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

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

323 

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

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

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

327 

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

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

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

331 

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

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

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

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

336 sz=_F( 0.156)), 

337 Identity = _lazy(_Identity_), 

338 

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

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

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

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

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

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

345 

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

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

348 

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

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

351 

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

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

354 

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

356 

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

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

359 

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

361 

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

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

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

365 

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

367 

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

369 

370 WGS84 = _lazy(_WGS84_), # unity 

371) 

372 

373 

374class Datum(_NamedEnumItem): 

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

376 ''' 

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

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

379 

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

381 '''New L{Datum}. 

382 

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

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

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

386 

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

388 

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

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

391 not a L{Transform}. 

392 ''' 

393 self._ellipsoid = ellipsoid or Datum._ellipsoid 

394 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

395 

396 self._transform = transform or Datum._transform 

397 _xinstanceof(Transform, transform=self.transform) 

398 

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

400 

401 def __eq__(self, other): 

402 '''Compare this and an other datum. 

403 

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

405 

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

407 ''' 

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

409 self.ellipsoid == other.ellipsoid and 

410 self.transform == other.transform) 

411 

412 def __hash__(self): 

413 return self._hash # memoized 

414 

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

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

417 

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

419 ''' 

420 _ = _xellipsoidall(point) 

421 return point.toDatum(self) 

422 

423 def ecef(self, Ecef=None): 

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

425 

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

427 

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

429 

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

431 

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

433 ''' 

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

435 

436 @Property_RO 

437 def ellipsoid(self): 

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

439 ''' 

440 return self._ellipsoid 

441 

442 @Property_RO 

443 def exactTM(self): 

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

445 ''' 

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

447 

448 @Property_RO 

449 def _hash(self): 

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

451 

452 @property_RO 

453 def isEllipsoidal(self): 

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

455 ''' 

456 return self.ellipsoid.isEllipsoidal 

457 

458 @property_RO 

459 def isOblate(self): 

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

461 ''' 

462 return self.ellipsoid.isOblate 

463 

464 @property_RO 

465 def isProlate(self): 

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

467 ''' 

468 return self.ellipsoid.isProlate 

469 

470 @property_RO 

471 def isSpherical(self): 

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

473 ''' 

474 return self.ellipsoid.isSpherical 

475 

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

477 '''Return this datum as a string. 

478 

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

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

481 this datum's name. 

482 

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

484 ''' 

485 t = [] if name is None else \ 

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

487 for a in (_ellipsoid_, _transform_): 

488 v = getattr(self, a) 

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

490 return sep.join(t) 

491 

492 @Property_RO 

493 def transform(self): 

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

495 ''' 

496 return self._transform 

497 

498 

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

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

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

502 

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

504 ''' 

505 if f is not None: 

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

507 if raiser and not E: 

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

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

510 return 

511 elif isinstance(a_earth, Datum): 

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

513 else: 

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

515 if raiser and not E: 

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

517 if D is None: 

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

519 inst._datum = D 

520 

521 

522def _earth_ellipsoid(earth, *name_raiser): 

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

524 ''' 

525 return Ellipsoids.Sphere if earth is R_M else ( 

526 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

527 _spherical_datum(earth, *name_raiser).ellipsoid) 

528 

529 

530def _ED2(radius, name): 

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

532 ''' 

533 D = Datums.Sphere 

534 E = D.ellipsoid 

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

536 n = _under(name) 

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

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

539 return E, D 

540 

541 

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

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

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

545 

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

547 ''' 

548 if isinstance(earth, Datum): 

549 D = earth 

550 else: 

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

552 if not E: 

553 n = raiser or _earth_ 

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

555 if D is None: 

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

557 if raiser and not D.isEllipsoidal: 

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

559 return D 

560 

561 

562def _EnD3(earth, name): 

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

564 ''' 

565 D = None 

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

567 E = earth 

568 n = _under(name or E.name) 

569 elif isinstance(earth, Datum): 

570 E = earth.ellipsoid 

571 n = _under(name or earth.name) 

572 D = earth 

573 elif _isRadius(earth): 

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

575 n = E.name 

576 elif isinstance(earth, a_f2Tuple): 

577 n = _under(name or earth.name) 

578 E = earth.ellipsoid(name=n) 

579 elif islistuple(earth, minum=2): 

580 E = Ellipsoids.Sphere 

581 a, f = earth[:2] 

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

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

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

585 else: 

586 n = E.name 

587 D = Datums.Sphere 

588 else: 

589 E, n = None, NN 

590 return E, n, D 

591 

592 

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

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

595 ''' 

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

597 

598 

599def _mean_radius(radius, *lats): 

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

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

602 ''' 

603 if radius is R_M: 

604 r = radius 

605 elif _isRadius(radius): 

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

607 else: 

608 E = _ellipsoidal_datum(radius).ellipsoid 

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

610 return r 

611 

612 

613def _negastr(name): # in .trf 

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

615 ''' 

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

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

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

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

620 # def _negastr(name): 

621 # n = name.translate(_MINUSxPLUS) 

622 # ... 

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

624 

625 

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

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

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

629 

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

631 ''' 

632 if _isRadius(earth): 

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

634 else: 

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

636 if raiser and not D.isSpherical: 

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

638 return D 

639 

640 

641class Datums(_NamedEnum): 

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

643 to accommodate the L{_LazyNamedEnumItem} properties. 

644 ''' 

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

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

647 ''' 

648 return Datum(Ellipsoids.get(ellipsoid_name), 

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

650 

651Datums = Datums(Datum) # PYCHOK singleton 

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

653# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

657Datums._assert( 

658 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

662 

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

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

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

666 

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

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

669 

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

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

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

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

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

675 

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

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

678 

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

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

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

682 

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

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

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

686 

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

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

689 

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

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

692 

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

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

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

696 

697 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

699 

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

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

702 

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

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

705 

706 # XXX psuedo-ellipsoids for spherical LatLon 

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

708 

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

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

711 

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

713 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

714 

715 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

716) 

717 

718_WGS84 = Datums.WGS84 

719assert _WGS84.ellipsoid is _EWGS84 

720# assert _WGS84.transform.isunity 

721 

722if __name__ == '__main__': 

723 

724 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

725 from pygeodesy.lazily import printf 

726 

727 # __doc__ of this file, force all into registery 

728 for r in (Datums, Transforms): 

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

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

731 

732# **) MIT License 

733# 

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

735# 

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

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

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

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

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

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

742# 

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

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

745# 

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

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

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

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

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

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

752# OTHER DEALINGS IN THE SOFTWARE.