Coverage for pygeodesy/datums.py: 94%

241 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-25 12:57 -0500

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 

71from pygeodesy.constants import R_M, _float as _F, _0_0, _0_26, _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 

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

80 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \ 

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

82 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \ 

83 _NAD83_, _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 

94__all__ = _ALL_LAZY.datums 

95__version__ = '24.02.24' 

96 

97_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

98_BD72_ = 'BD72' 

99_DHDN_ = 'DHDN' 

100_DHDNE_ = 'DHDNE' 

101_DHDNW_ = 'DHDNW' 

102_ED50_ = 'ED50' 

103_GDA2020_ = 'GDA2020' # in .trf 

104_Identity_ = 'Identity' 

105_Irl1965_ = 'Irl1965' 

106_Irl1975_ = 'Irl1975' 

107_MGI_ = 'MGI' 

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

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

110_NTF_ = 'NTF' 

111_OSGB36_ = 'OSGB36' 

112_Potsdam_ = 'Potsdam' 

113_RPS = radians(_1_0 / _3600_0) # radians per degree-second 

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

115_TokyoJapan_ = 'TokyoJapan' 

116 

117 

118class Transform(_NamedEnumItem): 

119 '''Helmert I{datum} transformation. 

120 

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

122 ''' 

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

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

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

126 

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

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

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

130 

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

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

133 

134 sx = _0_0 # x rotation (C{degree-seconds}) 

135 sy = _0_0 # y rotation (C{degree-seconds}) 

136 sz = _0_0 # z rotation (C{degree-seconds}) 

137 

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

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

140 '''New L{Transform}. 

141 

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

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

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

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

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

147 @kwarg sx: X rotation (C{degree-seconds}). 

148 @kwarg sy: Y rotation (C{degree-seconds}). 

149 @kwarg sz: Z rotation (C{degree-seconds}). 

150 

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

152 ''' 

153 if tx: 

154 self.tx = tx 

155 if ty: 

156 self.ty = ty 

157 if tz: 

158 self.tz = tz 

159 if sx: # secs to rads 

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

161 if sy: 

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

163 if sz: 

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

165 if s: 

166 self.s = s 

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

168 

169 self._register(Transforms, name) 

170 

171 def __eq__(self, other): 

172 '''Compare this and an other transform. 

173 

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

175 

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

177 ''' 

178 return self is other or (isinstance(other, Transform) and all( 

179 a == b for a, b in zip(self, other))) 

180 

181 def __hash__(self): 

182 return hash(x for x in self) 

183 

184 def __iter__(self): 

185 '''Yield the initial attribute values. 

186 ''' 

187 for _, x in self.items(): 

188 yield x 

189 

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

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

192 

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

194 

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

196 

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

198 ''' 

199 _ = _xellipsoidall(point) 

200 return point.toTransform(self) 

201 

202 def __neg__(self): 

203 return self.inverse() 

204 

205# def __sub__(self, other): 

206# _xinstanceof(Transform, other=other) 

207# 

208# def _sub(a, b): 

209# for n in _Names11: 

210# yield getattr(a, n) - getattr(b, n) 

211# 

212# return type(self)(_sub(self, other), name=_MINUS_) # .fsums._sub_op 

213 

214# @Property_RO 

215# def _hash(self): 

216# return hash(x for x in self) 

217 

218 def items(self): 

219 '''Yield each attribute as 2-tuple C{(name, value)}. 

220 ''' 

221 for n in _Names7: 

222 yield n, getattr(self, n) 

223 

224 def inverse(self, name=NN): 

225 '''Return the inverse of this transform. 

226 

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

228 

229 @return: Inverse (L{Transform}). 

230 

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

232 ''' 

233 d = dict((n, neg(v)) for n, v in self.items()) 

234 n = name or _minus(self.name) 

235 return type(self)(name=n, **d) 

236 

237 @Property_RO 

238 def isunity(self): 

239 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84? 

240 ''' 

241 return not any(s for s in self) 

242 

243 def _rps2(self, s_): 

244 '''(INTERNAL) Rotation in C{radians} and C{degree-seconds}. 

245 ''' 

246 return (_RPS * s_), s_ 

247 

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

249 '''Return this transform as a string. 

250 

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

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

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

254 this transform's name. 

255 

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

257 ''' 

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

259 

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

261 '''Transform a (geocentric) position, forward or inverse. 

262 

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

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

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

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

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

268 cartesian class and additional C{B{Vector}} 

269 keyword arguments to return the transformed 

270 point. 

271 

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

273 unless some B{C{Vector_and_kwds}} is specified. 

274 ''' 

275 if self.isunity: 

276 r = Vector3Tuple(x, y, z, name=self.name) 

277 else: 

278 xyz1 = x, y, z, _1_0 

279 s1 = self.s1 

280 if inverse: 

281 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

289 name=self.name) 

290 if Vector_and_kwds: 

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

292 if V: 

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

294 return r 

295 

296 

297class Transforms(_NamedEnum): 

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

299 to accommodate the L{_LazyNamedEnumItem} properties. 

300 ''' 

301 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

303 ''' 

304 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

305 

306Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

309Transforms._assert( 

310 BD72 = _lazy(_BD72_, tx=_F(106.868628), ty=_F(-52.297783), tz=_F(103.723893), 

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

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

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

314 s=_F( 1.2727)), 

315 Bessel1841 = _lazy(_Bessel1841_, tx=_F(-582.0), ty=_F(-105.0), tz=_F(-414.0), 

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

317 s=_F( -8.3)), 

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

319 DHDN = _lazy(_DHDN_, tx=_F(-591.28), ty=_F(-81.35), tz=_F(-396.39), 

320 sx=_F( 1.477), sy=_F( -0.0736), sz=_F( -1.458), 

321 s=_F( -9.82)), # Germany 

322 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2), 

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

324 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797), 

325 s=_F( -2.55)), # East Germany 

326 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2), 

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

328 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455), 

329 s=_F( -6.7)), # West Germany 

330 ED50 = _lazy(_ED50_, tx=_F(89.5), ty=_F(93.8), tz=_F(123.1), 

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

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

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

334 sz=_F( 0.156), s=_F(-1.2)), 

335 Identity = _lazy(_Identity_), 

336 Irl1965 = _lazy(_Irl1965_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), 

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

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

339 s=_F( -8.15)), 

340 Irl1975 = _lazy(_Irl1975_, tx=_F(-482.530), ty=_F(130.596), tz=_F(-564.557), 

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

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

343 s=_F( -8.15)), 

344 Krassovski1940 = _lazy(_Krassovski1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), 

345 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13), 

346 s=_F( -2.423)), # spelling 

347 Krassowsky1940 = _lazy(_Krassowsky1940_, tx=_F(-24.0), ty=_F(123.0), tz=_F(94.0), 

348 sx=_F( -0.02), sy= _0_26, sz=_F( 0.13), 

349 s=_F( -2.423)), # spelling 

350 MGI = _lazy(_MGI_, tx=_F(-577.326), ty=_F(-90.129), tz=_F(-463.920), 

351 sx=_F( 5.137), sy=_F( 1.474), sz=_F( 5.297), 

352 s=_F( -2.423)), # Austria 

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

354 NAD83 = _lazy(_NAD83_, tx=_F( 1.004), ty=_F(-1.910), tz=_F(-0.515), 

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

356 s=_F(-0.0015)), 

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

358 OSGB36 = _lazy(_OSGB36_, tx=_F(-446.448), ty=_F(125.157), tz=_F(-542.060), 

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

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

361 s=_F( 20.4894)), 

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

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

364 WGS84 = _lazy(_WGS84_), # unity 

365) 

366 

367 

368class Datum(_NamedEnumItem): 

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

370 ''' 

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

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

373 

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

375 '''New L{Datum}. 

376 

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

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

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

380 

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

382 

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

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

385 not a L{Transform}. 

386 ''' 

387 self._ellipsoid = ellipsoid or Datum._ellipsoid 

388 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

389 

390 self._transform = transform or Datum._transform 

391 _xinstanceof(Transform, transform=self.transform) 

392 

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

394 

395 def __eq__(self, other): 

396 '''Compare this and an other datum. 

397 

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

399 

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

401 ''' 

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

403 self.ellipsoid == other.ellipsoid and 

404 self.transform == other.transform) 

405 

406 def __hash__(self): 

407 return self._hash # memoized 

408 

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

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

411 

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

413 ''' 

414 _ = _xellipsoidall(point) 

415 return point.toDatum(self) 

416 

417 def ecef(self, Ecef=None): 

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

419 

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

421 

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

423 

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

425 

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

427 ''' 

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

429 

430 @Property_RO 

431 def ellipsoid(self): 

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

433 ''' 

434 return self._ellipsoid 

435 

436 @Property_RO 

437 def exactTM(self): 

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

439 ''' 

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

441 

442 @Property_RO 

443 def _hash(self): 

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

445 

446 @property_RO 

447 def isEllipsoidal(self): 

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

449 ''' 

450 return self.ellipsoid.isEllipsoidal 

451 

452 @property_RO 

453 def isOblate(self): 

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

455 ''' 

456 return self.ellipsoid.isOblate 

457 

458 @property_RO 

459 def isProlate(self): 

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

461 ''' 

462 return self.ellipsoid.isProlate 

463 

464 @property_RO 

465 def isSpherical(self): 

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

467 ''' 

468 return self.ellipsoid.isSpherical 

469 

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

471 '''Return this datum as a string. 

472 

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

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

475 this datum's name. 

476 

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

478 ''' 

479 t = [] if name is None else \ 

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

481 for a in (_ellipsoid_, _transform_): 

482 v = getattr(self, a) 

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

484 return sep.join(t) 

485 

486 @Property_RO 

487 def transform(self): 

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

489 ''' 

490 return self._transform 

491 

492 

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

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

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

496 

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

498 ''' 

499 if f is not None: 

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

501 if raiser and not E: 

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

503 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None): 

504 return 

505 elif isinstance(a_earth, Datum): 

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

507 else: 

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

509 if raiser and not E: 

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

511 if D is None: 

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

513 inst._datum = D 

514 

515 

516def _earth_ellipsoid(earth, *name_raiser): 

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

518 ''' 

519 return Ellipsoids.Sphere if earth is R_M else ( 

520 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

521 _spherical_datum(earth, *name_raiser).ellipsoid) 

522 

523 

524def _ED2(radius, name): 

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

526 ''' 

527 D = Datums.Sphere 

528 E = D.ellipsoid 

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

530 n = _under(name) 

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

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

533 return E, D 

534 

535 

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

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

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

539 

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

541 ''' 

542 if isinstance(earth, Datum): 

543 D = earth 

544 else: 

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

546 if not E: 

547 n = raiser or _earth_ 

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

549 if D is None: 

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

551 if raiser and not D.isEllipsoidal: 

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

553 return D 

554 

555 

556def _EnD3(earth, name): 

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

558 ''' 

559 D = None 

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

561 E = earth 

562 n = _under(name or E.name) 

563 elif isinstance(earth, Datum): 

564 E = earth.ellipsoid 

565 n = _under(name or earth.name) 

566 D = earth 

567 elif _isRadius(earth): 

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

569 n = E.name 

570 elif isinstance(earth, a_f2Tuple): 

571 n = _under(name or earth.name) 

572 E = earth.ellipsoid(name=n) 

573 elif islistuple(earth, minum=2): 

574 E = Ellipsoids.Sphere 

575 a, f = earth[:2] 

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

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

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

579 else: 

580 n = E.name 

581 D = Datums.Sphere 

582 else: 

583 E, n = None, NN 

584 return E, n, D 

585 

586 

587def _mean_radius(radius, *lats): 

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

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

590 ''' 

591 if radius is R_M: 

592 r = radius 

593 elif _isRadius(radius): 

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

595 else: 

596 E = _ellipsoidal_datum(radius).ellipsoid 

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

598 return r 

599 

600 

601def _minus(name): # in .trf 

602 '''(INTERNAL) Name of C{inverse} Xform. 

603 ''' 

604 m = _MINUS_ 

605 return name[len(m):] if name.startswith(m) else NN(m, name) 

606 

607 

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

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

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

611 

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

613 ''' 

614 if _isRadius(earth): 

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

616 else: 

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

618 if raiser and not D.isSpherical: 

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

620 return D 

621 

622 

623class Datums(_NamedEnum): 

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

625 to accommodate the L{_LazyNamedEnumItem} properties. 

626 ''' 

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

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

629 ''' 

630 return Datum(Ellipsoids.get(ellipsoid_name), 

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

632 

633Datums = Datums(Datum) # PYCHOK singleton 

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

635# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

639Datums._assert( 

640 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

644 

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

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

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

648 

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

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

651 

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

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

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

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

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

657 

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

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

660 

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

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

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

664 

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

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

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

668 

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

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

671 

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

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

674 

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

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

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

678 

679 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

681 

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

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

684 

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

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

687 

688 # XXX psuedo-ellipsoids for spherical LatLon 

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

690 

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

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

693 

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

695 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

696 

697 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

698) 

699 

700_WGS84 = Datums.WGS84 

701assert _WGS84.ellipsoid is _EWGS84 

702# assert _WGS84.transform.isunity 

703 

704if __name__ == '__main__': 

705 

706 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

707 from pygeodesy.lazily import printf 

708 

709 # __doc__ of this file, force all into registery 

710 for r in (Datums, Transforms): 

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

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

713 

714# **) MIT License 

715# 

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

717# 

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

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

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

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

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

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

724# 

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

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

727# 

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

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

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

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

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

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

734# OTHER DEALINGS IN THE SOFTWARE.