Coverage for pygeodesy/datums.py: 94%

239 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-05 16:22 -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.04' 

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 if name and not name.startswith(_MINUS_): 

166 self._register(Transforms, name) 

167 

168 def __eq__(self, other): 

169 '''Compare this and an other transform. 

170 

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

172 

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

174 ''' 

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

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

177 

178 def __hash__(self): 

179 return self._hash # memoized 

180 

181 def __iter__(self): 

182 '''Yield the attribute values. 

183 ''' 

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

185 yield x 

186 

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

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

189 

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

191 ''' 

192 m = _MODS.ellipsoidalBase 

193 _xinstanceof(m.CartesianEllipsoidalBase, point=point) 

194 return point.toTransform(self) 

195 

196 def __neg__(self): 

197 return self.inverse() 

198 

199 @Property_RO 

200 def _hash(self): 

201 return hash(x for x in self) 

202 

203 def items(self): 

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

205 ''' 

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

207 yield n, getattr(self, n) 

208 

209 def inverse(self, name=NN): 

210 '''Return the inverse of this transform. 

211 

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

213 

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

215 

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

217 ''' 

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

219 n = name or _minus(self.name) 

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

221 

222 @Property_RO 

223 def isunity(self): 

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

225 ''' 

226 return not any(s for s in self) 

227 

228 def _rps2(self, s_): 

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

230 ''' 

231 return (_RPS * s_), s_ 

232 

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

234 '''Return this transform as a string. 

235 

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

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

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

239 this transform's name. 

240 

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

242 ''' 

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

244 's1', 'rx', 'ry', 'rz', 

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

246 

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

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

249 

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

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

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

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

254 

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

256 ''' 

257 if self.isunity: 

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

259 

260 xyz1 = x, y, z, _1_0 

261 s1 = self.s1 

262 if inverse: 

263 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

271 name=self.name) 

272 

273 

274class Transforms(_NamedEnum): 

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

276 to accommodate the L{_LazyNamedEnumItem} properties. 

277 ''' 

278 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

280 ''' 

281 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

282 

283Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

286Transforms._assert( 

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

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

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

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

291 s=_F( 1.2727)), 

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

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

294 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

312 Identity = _lazy(_Identity_), 

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

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

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

316 s=_F( -8.15)), 

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

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

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

320 s=_F( -8.15)), 

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

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

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

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

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

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

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

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

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

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

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

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

333 s=_F(-0.0015)), 

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

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

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

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

338 s=_F( 20.4894)), 

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

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

341 WGS84 = _lazy(_WGS84_), # unity 

342) 

343 

344 

345class Datum(_NamedEnumItem): 

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

347 ''' 

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

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

350 

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

352 '''New L{Datum}. 

353 

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

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

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

357 

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

359 

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

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

362 not a L{Transform}. 

363 ''' 

364 self._ellipsoid = ellipsoid or Datum._ellipsoid 

365 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

366 

367 self._transform = transform or Datum._transform 

368 _xinstanceof(Transform, transform=self.transform) 

369 

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

371 

372 def __eq__(self, other): 

373 '''Compare this and an other datum. 

374 

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

376 

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

378 ''' 

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

380 self.ellipsoid == other.ellipsoid and 

381 self.transform == other.transform) 

382 

383 def __hash__(self): 

384 return self._hash # memoized 

385 

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

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

388 

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

390 ''' 

391 m = _MODS.ellipsoidalBase 

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

393 return point.toDatum(self) 

394 

395 def ecef(self, Ecef=None): 

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

397 

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

399 

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

401 

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

403 

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

405 ''' 

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

407 

408 @Property_RO 

409 def ellipsoid(self): 

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

411 ''' 

412 return self._ellipsoid 

413 

414 @Property_RO 

415 def exactTM(self): 

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

417 ''' 

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

419 

420 @Property_RO 

421 def _hash(self): 

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

423 

424 @property_RO 

425 def isEllipsoidal(self): 

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

427 ''' 

428 return self.ellipsoid.isEllipsoidal 

429 

430 @property_RO 

431 def isOblate(self): 

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

433 ''' 

434 return self.ellipsoid.isOblate 

435 

436 @property_RO 

437 def isProlate(self): 

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

439 ''' 

440 return self.ellipsoid.isProlate 

441 

442 @property_RO 

443 def isSpherical(self): 

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

445 ''' 

446 return self.ellipsoid.isSpherical 

447 

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

449 '''Return this datum as a string. 

450 

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

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

453 this datum's name. 

454 

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

456 ''' 

457 t = [] if name is None else \ 

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

459 for a in (_ellipsoid_, _transform_): 

460 v = getattr(self, a) 

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

462 return sep.join(t) 

463 

464 @Property_RO 

465 def transform(self): 

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

467 ''' 

468 return self._transform 

469 

470 

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

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

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

474 

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

476 ''' 

477 if f is not None: 

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

479 if raiser and not E: 

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

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

482 return 

483 elif isinstance(a_earth, Datum): 

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

485 else: 

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

487 if raiser and not E: 

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

489 if D is None: 

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

491 inst._datum = D 

492 

493 

494def _earth_ellipsoid(earth, *name_raiser): 

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

496 ''' 

497 return Ellipsoids.Sphere if earth is R_M else ( 

498 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

499 _spherical_datum(earth, *name_raiser).ellipsoid) 

500 

501 

502def _ED2(radius, name): 

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

504 ''' 

505 D = Datums.Sphere 

506 E = D.ellipsoid 

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

508 n = _under(name) 

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

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

511 return E, D 

512 

513 

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

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

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

517 

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

519 ''' 

520 if isinstance(earth, Datum): 

521 D = earth 

522 else: 

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

524 if not E: 

525 n = raiser or _earth_ 

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

527 if D is None: 

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

529 if raiser and not D.isEllipsoidal: 

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

531 return D 

532 

533 

534def _EnD3(earth, name): 

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

536 ''' 

537 D = None 

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

539 E = earth 

540 n = _under(name or E.name) 

541 elif isinstance(earth, Datum): 

542 E = earth.ellipsoid 

543 n = _under(name or earth.name) 

544 D = earth 

545 elif _isRadius(earth): 

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

547 n = E.name 

548 elif isinstance(earth, a_f2Tuple): 

549 n = _under(name or earth.name) 

550 E = earth.ellipsoid(name=n) 

551 elif islistuple(earth, minum=2): 

552 E = Ellipsoids.Sphere 

553 a, f = earth[:2] 

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

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

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

557 else: 

558 n = E.name 

559 D = Datums.Sphere 

560 else: 

561 E, n = None, NN 

562 return E, n, D 

563 

564 

565def _mean_radius(radius, *lats): 

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

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

568 ''' 

569 if radius is R_M: 

570 r = radius 

571 elif _isRadius(radius): 

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

573 else: 

574 E = _ellipsoidal_datum(radius).ellipsoid 

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

576 return r 

577 

578 

579def _minus(name): # in .trf 

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

581 ''' 

582 m = _MINUS_ 

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

584 

585 

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

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

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

589 

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

591 ''' 

592 if _isRadius(earth): 

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

594 else: 

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

596 if raiser and not D.isSpherical: 

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

598 return D 

599 

600 

601class Datums(_NamedEnum): 

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

603 to accommodate the L{_LazyNamedEnumItem} properties. 

604 ''' 

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

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

607 ''' 

608 return Datum(Ellipsoids.get(ellipsoid_name), 

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

610 

611Datums = Datums(Datum) # PYCHOK singleton 

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

613# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

617Datums._assert( 

618 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

622 

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

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

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

626 

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

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

629 

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

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

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

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

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

635 

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

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

638 

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

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

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

642 

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

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

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

646 

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

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

649 

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

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

652 

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

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

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

656 

657 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

659 

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

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

662 

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

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

665 

666 # XXX psuedo-ellipsoids for spherical LatLon 

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

668 

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

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

671 

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

673 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

674 

675 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

676) 

677 

678_WGS84 = Datums.WGS84 

679assert _WGS84.ellipsoid is _EWGS84 

680# assert _WGS84.transform.isunity 

681 

682if __name__ == '__main__': 

683 

684 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

685 from pygeodesy.lazily import printf 

686 

687 # __doc__ of this file, force all into registery 

688 for r in (Datums, Transforms): 

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

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

691 

692# **) MIT License 

693# 

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

695# 

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

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

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

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

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

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

702# 

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

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

705# 

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

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

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

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

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

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

712# OTHER DEALINGS IN THE SOFTWARE.