Coverage for pygeodesy/datums.py: 94%

238 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-09 17:20 -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.869, ty=-52.2978, tz=103.724, s1=1.0, rx=-1.63174e-06, ry=-2.21538e-06, rz=-8.93114e-06, s=1.2727, sx=-0.33657, sy=-0.456955, sz=-1.84218) 

47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, s1=0.999992, rx=-5.04206e-06, ry=-1.69685e-06, rz=1.49323e-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.56823e-07, rz=-7.06858e-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=0.999997, rx=2.61799e-07, ry=-2.76344e-07, rz=1.35602e-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.999993, rx=-9.79324e-07, ry=-2.18166e-07, rz=1.19022e-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=0.999999, rx=0.0, ry=0.0, rz=7.56309e-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.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-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.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-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=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-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=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13) 

58@var Transforms.MGI: Transform(name='MGI', tx=-577.326, ty=-90.129, tz=-463.92, s1=0.999998, rx=2.49049e-05, ry=7.14615e-06, rz=2.56806e-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.29445e-07, ry=1.64837e-09, rz=5.33295e-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.448, ty=125.157, tz=-542.06, s1=1.00002, rx=-7.2819e-07, ry=-1.19749e-06, rz=-4.08262e-06, s=20.4894, 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.68587e-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 

77from pygeodesy.fmath import fdot, fmean, Fmt 

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

79 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \ 

80 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _MINUS_, \ 

81 _Krassovski1940_, _Krassowsky1940_, _NAD27_, _NAD83_, _s_, \ 

82 _Sphere_, _spherical_, _sx_, _sy_, _sz_, _transform_, _tx_, \ 

83 _ty_, _tz_, _UNDER_, _WGS72_, _WGS84_, _under 

84from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

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

87from pygeodesy.props import Property_RO, property_RO 

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

89from pygeodesy.units import _isRadius, Radius_, radians 

90 

91# from math import radians # from .units 

92 

93__all__ = _ALL_LAZY.datums 

94__version__ = '24.02.06' 

95 

96_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

97_BD72_ = 'BD72' 

98_DHDN_ = 'DHDN' 

99_DHDNE_ = 'DHDNE' 

100_DHDNW_ = 'DHDNW' 

101_ED50_ = 'ED50' 

102_GDA2020_ = 'GDA2020' 

103_Identity_ = 'Identity' 

104_Irl1965_ = 'Irl1965' 

105_Irl1975_ = 'Irl1975' 

106_MGI_ = 'MGI' 

107_NTF_ = 'NTF' 

108_OSGB36_ = 'OSGB36' 

109_Potsdam_ = 'Potsdam' 

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

111_TokyoJapan_ = 'TokyoJapan' 

112 

113 

114class Transform(_NamedEnumItem): 

115 '''Helmert I{datum} transformation. 

116 

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

118 ''' 

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

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

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

122 

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

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

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

126 

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

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

129 

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

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

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

133 

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

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

136 '''New L{Transform}. 

137 

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

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

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

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

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

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

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

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

146 

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

148 ''' 

149 if tx: 

150 self.tx = tx 

151 if ty: 

152 self.ty = ty 

153 if tz: 

154 self.tz = tz 

155 if sx: # secs to rads 

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

157 if sy: 

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

159 if sz: 

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

161 if s: 

162 self.s = s 

163 self.s1 = _F(s * 1e-6 + _1_0) # normalize ppm to (s + 1) 

164 

165 self._register(Transforms, name) 

166 

167 def __eq__(self, other): 

168 '''Compare this and an other transform. 

169 

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

171 

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

173 ''' 

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

175 s == t for s, t in zip(self, other))) 

176 

177 def __hash__(self): 

178 return self._hash # memoized 

179 

180 def __iter__(self): 

181 '''Yield the attribute values. 

182 ''' 

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

184 yield x 

185 

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

187 '''Helmert-transform an I{ellipsoidal} cartesian B{C{point}}. 

188 

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

190 ''' 

191 m = _MODS.ellipsoidalBase 

192 _xinstanceof(m.CartesianEllipsoidalBase, point=point) 

193 return point.toTransform(self) 

194 

195 def __neg__(self): 

196 return self.inverse() 

197 

198 @Property_RO 

199 def _hash(self): 

200 return hash(x for x in self) 

201 

202 def items(self): 

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

204 ''' 

205 for n in (_tx_, _ty_, _tz_, _s_, _sx_, _sy_, _sz_): 

206 yield n, getattr(self, n) 

207 

208 def inverse(self, name=NN): 

209 '''Return the inverse of this transform. 

210 

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

212 

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

214 

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

216 ''' 

217 d = dict((n, -v) for n, v in self.items()) 

218 n = name or _minus(self.name) 

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

220 

221 @Property_RO 

222 def isunity(self): 

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

224 ''' 

225 return not any(s for s in self) 

226 

227 def _rps2(self, s_): 

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

229 ''' 

230 return (_RPS * s_), s_ 

231 

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

233 '''Return this transform as a string. 

234 

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

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

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

238 this transform's name. 

239 

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

241 ''' 

242 return self._instr(name, prec, _tx_, _ty_, _tz_, 

243 's1', 'rx', 'ry', 'rz', 

244 _s_, _sx_, _sy_, _sz_, fmt=fmt) 

245 

246 def transform(self, x, y, z, inverse=False): 

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

248 

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

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

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

252 @kwarg inverse: Optional direction, forward or inverse (C{bool}). 

253 

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

255 ''' 

256 if self.isunity: 

257 return Vector3Tuple(x, y, z, name=self.name) 

258 

259 xyz1 = x, y, z, _1_0 

260 s1 = self.s1 

261 if inverse: 

262 xyz1 = map2(neg, xyz1) 

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

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

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

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

267 return Vector3Tuple(fdot(xyz1, s1, -self.rz, self.ry, self.tx), 

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

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

270 name=self.name) 

271 

272 

273class Transforms(_NamedEnum): 

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

275 to accommodate the L{_LazyNamedEnumItem} properties. 

276 ''' 

277 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

279 ''' 

280 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

281 

282Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

285Transforms._assert( 

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

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

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

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

290 s=_F( 1.2727)), 

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

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

293 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

311 Identity = _lazy(_Identity_), 

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

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

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

315 s=_F( -8.15)), 

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

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

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

319 s=_F( -8.15)), 

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

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

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

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

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

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

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

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

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

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

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

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

332 s=_F(-0.0015)), 

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

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

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

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

337 s=_F( 20.4894)), 

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

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

340 WGS84 = _lazy(_WGS84_), # unity 

341) 

342 

343 

344class Datum(_NamedEnumItem): 

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

346 ''' 

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

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

349 

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

351 '''New L{Datum}. 

352 

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

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

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

356 

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

358 

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

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

361 not a L{Transform}. 

362 ''' 

363 self._ellipsoid = ellipsoid or Datum._ellipsoid 

364 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

365 

366 self._transform = transform or Datum._transform 

367 _xinstanceof(Transform, transform=self.transform) 

368 

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

370 

371 def __eq__(self, other): 

372 '''Compare this and an other datum. 

373 

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

375 

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

377 ''' 

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

379 self.ellipsoid == other.ellipsoid and 

380 self.transform == other.transform) 

381 

382 def __hash__(self): 

383 return self._hash # memoized 

384 

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

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

387 

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

389 ''' 

390 m = _MODS.ellipsoidalBase 

391 _xinstanceof(m.CartesianEllipsoidalBase, m.LatLonEllipsoidalBase, point=point) 

392 return point.toDatum(self) 

393 

394 def ecef(self, Ecef=None): 

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

396 

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

398 

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

400 

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

402 

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

404 ''' 

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

406 

407 @Property_RO 

408 def ellipsoid(self): 

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

410 ''' 

411 return self._ellipsoid 

412 

413 @Property_RO 

414 def exactTM(self): 

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

416 ''' 

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

418 

419 @Property_RO 

420 def _hash(self): 

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

422 

423 @property_RO 

424 def isEllipsoidal(self): 

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

426 ''' 

427 return self.ellipsoid.isEllipsoidal 

428 

429 @property_RO 

430 def isOblate(self): 

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

432 ''' 

433 return self.ellipsoid.isOblate 

434 

435 @property_RO 

436 def isProlate(self): 

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

438 ''' 

439 return self.ellipsoid.isProlate 

440 

441 @property_RO 

442 def isSpherical(self): 

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

444 ''' 

445 return self.ellipsoid.isSpherical 

446 

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

448 '''Return this datum as a string. 

449 

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

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

452 this datum's name. 

453 

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

455 ''' 

456 t = [] if name is None else \ 

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

458 for a in (_ellipsoid_, _transform_): 

459 v = getattr(self, a) 

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

461 return sep.join(t) 

462 

463 @Property_RO 

464 def transform(self): 

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

466 ''' 

467 return self._transform 

468 

469 

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

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

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

473 

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

475 ''' 

476 if f is not None: 

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

478 if raiser and not E: 

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

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

481 return 

482 elif isinstance(a_earth, Datum): 

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

484 else: 

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

486 if raiser and not E: 

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

488 if D is None: 

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

490 inst._datum = D 

491 

492 

493def _earth_ellipsoid(earth, *name_raiser): 

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

495 ''' 

496 return Ellipsoids.Sphere if earth is R_M else ( 

497 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

498 _spherical_datum(earth, *name_raiser).ellipsoid) 

499 

500 

501def _ED2(radius, name): 

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

503 ''' 

504 D = Datums.Sphere 

505 E = D.ellipsoid 

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

507 n = _under(name) 

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

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

510 return E, D 

511 

512 

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

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

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

516 

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

518 ''' 

519 if isinstance(earth, Datum): 

520 D = earth 

521 else: 

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

523 if not E: 

524 n = raiser or _earth_ 

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

526 if D is None: 

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

528 if raiser and not D.isEllipsoidal: 

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

530 return D 

531 

532 

533def _EnD3(earth, name): 

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

535 ''' 

536 D = None 

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

538 E = earth 

539 n = _under(name or E.name) 

540 elif isinstance(earth, Datum): 

541 E = earth.ellipsoid 

542 n = _under(name or earth.name) 

543 D = earth 

544 elif _isRadius(earth): 

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

546 n = E.name 

547 elif isinstance(earth, a_f2Tuple): 

548 n = _under(name or earth.name) 

549 E = earth.ellipsoid(name=n) 

550 elif islistuple(earth, minum=2): 

551 E = Ellipsoids.Sphere 

552 a, f = earth[:2] 

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

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

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

556 else: 

557 n = E.name 

558 D = Datums.Sphere 

559 else: 

560 E, n = None, NN 

561 return E, n, D 

562 

563 

564def _mean_radius(radius, *lats): 

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

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

567 ''' 

568 if radius is R_M: 

569 r = radius 

570 elif _isRadius(radius): 

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

572 else: 

573 E = _ellipsoidal_datum(radius).ellipsoid 

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

575 return r 

576 

577 

578def _minus(name): # in .trf 

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

580 ''' 

581 m = _MINUS_ 

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

583 

584 

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

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

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

588 

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

590 ''' 

591 if _isRadius(earth): 

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

593 else: 

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

595 if raiser and not D.isSpherical: 

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

597 return D 

598 

599 

600class Datums(_NamedEnum): 

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

602 to accommodate the L{_LazyNamedEnumItem} properties. 

603 ''' 

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

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

606 ''' 

607 return Datum(Ellipsoids.get(ellipsoid_name), 

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

609 

610Datums = Datums(Datum) # PYCHOK singleton 

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

612# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

616Datums._assert( 

617 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

621 

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

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

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

625 

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

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

628 

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

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

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

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

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

634 

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

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

637 

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

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

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

641 

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

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

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

645 

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

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

648 

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

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

651 

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

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

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

655 

656 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

658 

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

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

661 

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

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

664 

665 # XXX psuedo-ellipsoids for spherical LatLon 

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

667 

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

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

670 

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

672 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

673 

674 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

675) 

676 

677_WGS84 = Datums.WGS84 

678assert _WGS84.ellipsoid is _EWGS84 

679# assert _WGS84.transform.isunity 

680 

681if __name__ == '__main__': 

682 

683 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

684 from pygeodesy.lazily import printf 

685 

686 # __doc__ of this file, force all into registery 

687 for r in (Datums, Transforms): 

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

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

690 

691# **) MIT License 

692# 

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

694# 

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

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

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

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

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

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

701# 

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

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

704# 

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

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

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

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

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

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

711# OTHER DEALINGS IN THE SOFTWARE.