Coverage for pygeodesy/datums.py: 94%

218 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -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.86863, ty=-52.29778, tz=103.72389, rx=-0, ry=-0, rz=-0.00001, s=1.2727, s1=1, sx=-0.33657, sy=-0.45696, sz=-1.84218) 

47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, rx=-0.00001, ry=-0, rz=0.00001, s=-8.3, s1=0.99999, sx=-1.04, sy=-0.35, sz=3.08) 

48@var Transforms.Clarke1866: Transform(name='Clarke1866', tx=8, ty=-160, tz=-176, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

49@var Transforms.DHDN: Transform(name='DHDN', tx=-591.28, ty=-81.35, tz=-396.39, rx=0.00001, ry=-0, rz=-0.00001, s=-9.82, s1=0.99999, sx=1.477, sy=-0.0736, sz=-1.458) 

50@var Transforms.ED50: Transform(name='ED50', tx=89.5, ty=93.8, tz=123.1, rx=0, ry=0, rz=0, s=-1.2, s1=1, sx=0, sy=0, sz=0.156) 

51@var Transforms.Identity: Transform(name='Identity', tx=0, ty=0, tz=0, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

52@var Transforms.Irl1965: Transform(name='Irl1965', tx=-482.53, ty=130.596, tz=-564.557, rx=0.00001, ry=0, rz=0, s=-8.15, s1=0.99999, sx=1.042, sy=0.214, sz=0.631) 

53@var Transforms.Irl1975: Transform(name='Irl1975', tx=-482.53, ty=130.596, tz=-564.557, rx=-0.00001, ry=-0, rz=-0, s=-1.1, s1=1, sx=-1.042, sy=-0.214, sz=-0.631) 

54@var Transforms.Krassovski1940: Transform(name='Krassovski1940', tx=-24, ty=123, tz=94, rx=-0, ry=0, rz=0, s=-2.423, s1=1, sx=-0.02, sy=0.26, sz=0.13) 

55@var Transforms.Krassowsky1940: Transform(name='Krassowsky1940', tx=-24, ty=123, tz=94, rx=-0, ry=0, rz=0, s=-2.423, s1=1, sx=-0.02, sy=0.26, sz=0.13) 

56@var Transforms.MGI: Transform(name='MGI', tx=-577.326, ty=-90.129, tz=-463.92, rx=0.00002, ry=0.00001, rz=0.00003, s=-2.423, s1=1, sx=5.137, sy=1.474, sz=5.297) 

57@var Transforms.NAD27: Transform(name='NAD27', tx=8, ty=-160, tz=-176, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

58@var Transforms.NAD83: Transform(name='NAD83', tx=1.004, ty=-1.91, tz=-0.515, rx=0, ry=0, rz=0, s=-0.0015, s1=1, sx=0.0267, sy=0.00034, sz=0.011) 

59@var Transforms.NTF: Transform(name='NTF', tx=-168, ty=-60, tz=320, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

60@var Transforms.OSGB36: Transform(name='OSGB36', tx=-446.448, ty=125.157, tz=-542.06, rx=-0, ry=-0, rz=-0, s=20.4894, s1=1.00002, sx=-0.1502, sy=-0.247, sz=-0.8421) 

61@var Transforms.TokyoJapan: Transform(name='TokyoJapan', tx=148, ty=-507, tz=-685, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

62@var Transforms.WGS72: Transform(name='WGS72', tx=0, ty=0, tz=-4.5, rx=0, ry=0, rz=0, s=-0.22, s1=1, sx=0, sy=0, sz=0.554) 

63@var Transforms.WGS84: Transform(name='WGS84', tx=0, ty=0, tz=0, rx=0, ry=0, rz=0, s=0, s1=1, sx=0, sy=0, sz=0) 

64''' 

65# make sure int/int division yields float quotient, see .basics 

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

67 

68from pygeodesy.basics import islistuple, isscalar, map2, neg, _xinstanceof 

69from pygeodesy.constants import R_M, _float as _F, _0_0, _0_26, _1_0, _2_0, _8_0, _3600_0 

70from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, \ 

71 _EWGS84, Vector3Tuple 

72from pygeodesy.errors import _IsnotError, _TypeError, _xattr 

73from pygeodesy.fmath import fdot, fmean, Fmt 

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

75 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \ 

76 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \ 

77 _Krassowsky1940_, _NAD27_, _NAD83_, _s_, _Sphere_, _spherical_, \ 

78 _sx_, _sy_, _sz_, _transform_, _tx_, _ty_, _tz_, _UNDER_, \ 

79 _WGS72_, _WGS84_, _under 

80from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

81from pygeodesy.named import _NamedEnum, _NamedEnumItem, Property_RO, property_RO, \ 

82 _lazyNamedEnumItem as _lazy 

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

84# from pygeodesy.props import Property_RO, property_RO # from .named 

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

86from pygeodesy.units import radians, Radius_ 

87 

88# from math import radians # from .units 

89 

90__all__ = _ALL_LAZY.datums 

91__version__ = '23.11.05' 

92 

93_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

94_BD72_ = 'BD72' 

95_DHDN_ = 'DHDN' 

96_DHDNE_ = 'DHDNE' 

97_DHDNW_ = 'DHDNW' 

98_ED50_ = 'ED50' 

99_GDA2020_ = 'GDA2020' 

100_Identity_ = 'Identity' 

101_Inverse_ = 'Inverse' 

102_Irl1965_ = 'Irl1965' 

103_Irl1975_ = 'Irl1975' 

104_MGI_ = 'MGI' 

105_NTF_ = 'NTF' 

106_OSGB36_ = 'OSGB36' 

107_Potsdam_ = 'Potsdam' 

108_RPS = radians(_1_0 / _3600_0) # radians per degree_second 

109_TokyoJapan_ = 'TokyoJapan' 

110 

111 

112def _rps2(s_): 

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

114 ''' 

115 return (_RPS * s_), s_ 

116 

117 

118class Transform(_NamedEnumItem): 

119 '''Helmert transformation. 

120 

121 @see: L{Helmert7Tuple}. 

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 (degree seconds) 

135 sy = _0_0 # y rotation (degree seconds) 

136 sz = _0_0 # z rotation (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: Optional X translation (C{meter}). 

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

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

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

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

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

149 @kwarg sz: Optional 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 = _rps2(sx) 

161 if sy: 

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

163 if sz: 

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

165 if s: 

166 self.s = s 

167 self.s1 = _F(s * 1e-6 + _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) 

179 and self.tx == other.tx 

180 and self.ty == other.ty 

181 and self.tz == other.tz 

182 and self.rx == other.rx 

183 and self.ry == other.ry 

184 and self.rz == other.rz 

185 and self.s == other.s) 

186 

187 def __hash__(self): 

188 return self._hash # memoized 

189 

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

191 '''Helmert-transform a cartesian B{C{other}}. 

192 

193 @raise TypeError: Invalid B{C{other}}. 

194 ''' 

195 try: # only CartesianBase 

196 return other.toTransform(self) 

197 except AttributeError: 

198 pass 

199 raise _IsnotError(_cartesian_, other=other) 

200 

201 @Property_RO 

202 def _hash(self): 

203 return hash((self.rx, self.ry, self.rz, self.s, 

204 self.tx, self.ty, self.tz)) 

205 

206 def inverse(self, name=NN): 

207 '''Return the inverse of this transform. 

208 

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

210 

211 @return: Inverse (Transform). 

212 

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

214 ''' 

215 return Transform(name=name or (self.name + _Inverse_), 

216 tx=-self.tx, ty=-self.ty, tz=-self.tz, 

217 sx=-self.sx, sy=-self.sy, sz=-self.sz, s=-self.s) 

218 

219 def toStr(self, prec=5, name=NN, **unused): # PYCHOK expected 

220 '''Return this transform as a string. 

221 

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

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

224 this transform's name. 

225 

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

227 ''' 

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

229 'rx', 'ry', 'rz', _s_, 's1', 

230 _sx_, _sy_, _sz_) 

231 

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

233 '''Transform a (geocentric) Cartesian point, forward or inverse. 

234 

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

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

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

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

239 

240 @return: A L{Vector3Tuple}C{(x, y, z)}, transformed. 

241 ''' 

242 xyz1 = x, y, z, _1_0 

243 s1 = self.s1 

244 if inverse: 

245 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

253 name=self.name) 

254 

255 

256class Transforms(_NamedEnum): 

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

258 to accommodate the L{_LazyNamedEnumItem} properties. 

259 ''' 

260 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

262 ''' 

263 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

264 

265Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

268Transforms._assert( 

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

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

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

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

273 s=_F( 1.2727)), 

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

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

276 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

294 Identity = _lazy(_Identity_), 

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

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

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

298 s=_F( -8.15)), 

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

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

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

302 s=_F( -8.15)), 

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

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

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

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

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

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

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

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

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

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

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

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

315 s=_F(-0.0015)), 

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

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

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

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

320 s=_F( 20.4894)), 

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

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

323 WGS84 = _lazy(_WGS84_), # unity 

324) 

325 

326 

327class Datum(_NamedEnumItem): 

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

329 ''' 

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

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

332 

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

334 '''New L{Datum}. 

335 

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

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

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

339 

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

341 

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

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

344 not a L{Transform}. 

345 ''' 

346 self._ellipsoid = ellipsoid or Datum._ellipsoid 

347 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

348 

349 self._transform = transform or Datum._transform 

350 _xinstanceof(Transform, transform=self.transform) 

351 

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

353 

354 def __eq__(self, other): 

355 '''Compare this and an other datum. 

356 

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

358 

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

360 ''' 

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

362 self.ellipsoid == other.ellipsoid and 

363 self.transform == other.transform) 

364 

365 def __hash__(self): 

366 return self._hash # memoized 

367 

368 def __matmul__(self, other): # PYCHOK Python 3.5+ 

369 '''Convert cartesian or ellipsoidal B{C{other}} to this datum. 

370 

371 @raise TypeError: Invalid B{C{other}}. 

372 ''' 

373 try: # only CartesianBase and LatLonEllipsoidalBase 

374 return other.toDatum(self) 

375 except AttributeError: 

376 pass 

377 raise _IsnotError(_cartesian_, _ellipsoidal_, other=other) 

378 

379 def ecef(self, Ecef=None): 

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

381 

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

383 

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

385 

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

387 

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

389 ''' 

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

391 

392 @Property_RO 

393 def ellipsoid(self): 

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

395 ''' 

396 return self._ellipsoid 

397 

398 @Property_RO 

399 def exactTM(self): 

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

401 ''' 

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

403 

404 @Property_RO 

405 def _hash(self): 

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

407 

408 @property_RO 

409 def isEllipsoidal(self): 

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

411 ''' 

412 return self.ellipsoid.isEllipsoidal 

413 

414 @property_RO 

415 def isOblate(self): 

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

417 ''' 

418 return self.ellipsoid.isOblate 

419 

420 @property_RO 

421 def isProlate(self): 

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

423 ''' 

424 return self.ellipsoid.isProlate 

425 

426 @property_RO 

427 def isSpherical(self): 

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

429 ''' 

430 return self.ellipsoid.isSpherical 

431 

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

433 '''Return this datum as a string. 

434 

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

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

437 this datum's name. 

438 

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

440 ''' 

441 t = [] if name is None else \ 

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

443 for a in (_ellipsoid_, _transform_): 

444 v = getattr(self, a) 

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

446 return sep.join(t) 

447 

448 @Property_RO 

449 def transform(self): 

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

451 ''' 

452 return self._transform 

453 

454 

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

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

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

458 

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

460 ''' 

461 if f is not None: 

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

463 if raiser and not E: 

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

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

466 return 

467 elif isinstance(a_earth, Datum): 

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

469 else: 

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

471 if raiser and not E: 

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

473 if D is None: 

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

475 inst._datum = D 

476 

477 

478def _earth_ellipsoid(earth, *name_raiser): 

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

480 ''' 

481 return Ellipsoids.Sphere if earth is R_M else ( 

482 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

483 _spherical_datum(earth, *name_raiser).ellipsoid) 

484 

485 

486def _ED2(radius, name): 

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

488 ''' 

489 D = Datums.Sphere 

490 E = D.ellipsoid 

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

492 n = _under(name) 

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

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

495 return E, D 

496 

497 

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

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

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

501 

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

503 ''' 

504 if isinstance(earth, Datum): 

505 D = earth 

506 else: 

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

508 if not E: 

509 n = raiser or _earth_ 

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

511 if D is None: 

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

513 if raiser and not D.isEllipsoidal: 

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

515 return D 

516 

517 

518def _EnD3(earth, name): 

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

520 ''' 

521 D = None 

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

523 E = earth 

524 n = _under(name or E.name) 

525 elif isinstance(earth, Datum): 

526 E = earth.ellipsoid 

527 n = _under(name or earth.name) 

528 D = earth 

529 elif isscalar(earth): 

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

531 n = E.name 

532 elif isinstance(earth, a_f2Tuple): 

533 n = _under(name or earth.name) 

534 E = earth.ellipsoid(name=n) 

535 elif islistuple(earth, minum=2): 

536 E = Ellipsoids.Sphere 

537 a, f = earth[:2] 

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

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

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

541 else: 

542 n = E.name 

543 D = Datums.Sphere 

544 else: 

545 E, n = None, NN 

546 return E, n, D 

547 

548 

549def _mean_radius(radius, *lats): 

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

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

552 ''' 

553 if radius is R_M: 

554 r = radius 

555 elif isscalar(radius): 

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

557 else: 

558 E = _ellipsoidal_datum(radius).ellipsoid 

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

560 return r 

561 

562 

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

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

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

566 

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

568 ''' 

569 if isscalar(earth): 

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

571 else: 

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

573 if raiser and not D.isSpherical: 

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

575 return D 

576 

577 

578class Datums(_NamedEnum): 

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

580 to accommodate the L{_LazyNamedEnumItem} properties. 

581 ''' 

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

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

584 ''' 

585 return Datum(Ellipsoids.get(ellipsoid_name), 

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

587 

588Datums = Datums(Datum) # PYCHOK singleton 

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

590# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

594Datums._assert( 

595 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

599 

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

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

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

603 

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

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

606 

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

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

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

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

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

612 

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

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

615 

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

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

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

619 

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

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

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

623 

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

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

626 

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

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

629 

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

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

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

633 

634 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

636 

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

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

639 

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

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

642 

643 # XXX psuedo-ellipsoids for spherical LatLon 

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

645 

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

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

648 

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

650 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

651 

652 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

653) 

654 

655_WGS84 = Datums.WGS84 

656assert _WGS84.ellipsoid is _EWGS84 

657 

658if __name__ == '__main__': 

659 

660 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

661 from pygeodesy.lazily import printf 

662 

663 # __doc__ of this file, force all into registery 

664 for r in (Datums, Transforms): 

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

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

667 

668# **) MIT License 

669# 

670# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

671# 

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

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

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

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

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

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

678# 

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

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

681# 

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

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

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

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

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

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

688# OTHER DEALINGS IN THE SOFTWARE.