Coverage for pygeodesy/datums.py: 97%

185 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-06-07 08:37 -0400

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Datums and transformations thereof. 

5 

6Classes L{Datum} and L{Transform} and registries L{Datums} and L{Transforms}, respectively. 

7 

8Pure Python implementation of geodesy tools for ellipsoidal earth models, including datums 

9and ellipsoid parameters for different geographic coordinate systems and methods for 

10converting between them and to cartesian coordinates. Transcoded from JavaScript originals by 

11I{(C) Chris Veness 2005-2016} and published under the same MIT Licence**, see U{latlon-ellipsoidal.js 

12<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}. 

13 

14Historical geodetic datums: a latitude/longitude point defines a geographic location on or 

15above/below the earth’s surface, measured in degrees from the equator, from the International 

16Reference Meridian, in meters above the ellipsoid and based on a given datum. The datum in turn 

17is based on a reference ellipsoid and tied to geodetic survey reference points. 

18 

19Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems), 

20but previously various reference ellipsoids and datum references were used. 

21 

22The UK Ordnance Survey National Grid References are still based on the otherwise historical 

23OSGB36 datum, q.v. U{"A Guide to Coordinate Systems in Great Britain", Section 6 

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

25 

26@var Datums.BD72: Datum(name='BD72', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.BD72) 

27@var Datums.DHDN: Datum(name='DHDN', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.DHDN) 

28@var Datums.ED50: Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50) 

29@var Datums.GDA2020: Datum(name='GDA2020', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84) 

30@var Datums.GRS80: Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84) 

31@var Datums.Irl1975: Datum(name='Irl1975', ellipsoid=Ellipsoids.AiryModified, transform=Transforms.Irl1975) 

32@var Datums.Krassovski1940: Datum(name='Krassovski1940', ellipsoid=Ellipsoids.Krassovski1940, transform=Transforms.Krassovski1940) 

33@var Datums.Krassowsky1940: Datum(name='Krassowsky1940', ellipsoid=Ellipsoids.Krassowsky1940, transform=Transforms.Krassowsky1940) 

34@var Datums.MGI: Datum(name='MGI', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.MGI) 

35@var Datums.NAD27: Datum(name='NAD27', ellipsoid=Ellipsoids.Clarke1866, transform=Transforms.NAD27) 

36@var Datums.NAD83: Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83) 

37@var Datums.NTF: Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF) 

38@var Datums.OSGB36: Datum(name='OSGB36', ellipsoid=Ellipsoids.Airy1830, transform=Transforms.OSGB36) 

39@var Datums.Potsdam: Datum(name='Potsdam', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841) 

40@var Datums.Sphere: Datum(name='Sphere', ellipsoid=Ellipsoids.Sphere, transform=Transforms.WGS84) 

41@var Datums.TokyoJapan: Datum(name='TokyoJapan', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.TokyoJapan) 

42@var Datums.WGS72: Datum(name='WGS72', ellipsoid=Ellipsoids.WGS72, transform=Transforms.WGS72) 

43@var Datums.WGS84: Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84) 

44 

45@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) 

46@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) 

47@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) 

48@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) 

49@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) 

50@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) 

51@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) 

52@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) 

53@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) 

54@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) 

55@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) 

56@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) 

57@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) 

58@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) 

59@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) 

60@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) 

61@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) 

62@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) 

63''' 

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

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

66 

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

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

69from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, Vector3Tuple 

70from pygeodesy.errors import _IsnotError, _xattr 

71from pygeodesy.fmath import fdot, fmean, Fmt 

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

73 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _ellipsoid_, \ 

74 _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \ 

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

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

77 _WGS72_, _WGS84_, _UNDER 

78from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

79from pygeodesy.named import _NamedEnum, _NamedEnumItem, \ 

80 _lazyNamedEnumItem as _lazy, Property_RO 

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

82# from pygeodesy.props import Property_RO # from .named 

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

84from pygeodesy.units import radians, Radius_ 

85 

86# from math import radians # from .units 

87 

88__all__ = _ALL_LAZY.datums 

89__version__ = '23.05.31' 

90 

91_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

92_BD72_ = 'BD72' 

93_DHDN_ = 'DHDN' 

94_earth_ = 'earth' 

95_ED50_ = 'ED50' 

96_GDA2020_ = 'GDA2020' 

97_Identity_ = 'Identity' 

98_Inverse_ = 'Inverse' 

99_Irl1965_ = 'Irl1965' 

100_Irl1975_ = 'Irl1975' 

101_MGI_ = 'MGI' 

102_NTF_ = 'NTF' 

103_OSGB36_ = 'OSGB36' 

104_Potsdam_ = 'Potsdam' 

105_TokyoJapan_ = 'TokyoJapan' 

106 

107_r_s1 = radians(1 / _3600_0) # 1 degree second to radians 

108 

109 

110def _r_s2(s_): 

111 '''(INTERNAL) rotation in C{radians} and C{degree seconds}. 

112 ''' 

113 return _F(_r_s1 * s_), s_ 

114 

115 

116class Transform(_NamedEnumItem): 

117 '''Helmert transformation. 

118 

119 @see: L{Helmert7Tuple}. 

120 ''' 

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

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

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

124 

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

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

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

128 

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

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

131 

132 sx = _0_0 # x rotation (degree seconds) 

133 sy = _0_0 # y rotation (degree seconds) 

134 sz = _0_0 # z rotation (degree seconds) 

135 

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

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

138 '''New L{Transform}. 

139 

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

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

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

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

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

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

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

147 @kwarg sz: Optional Z rotation (C{degree seconds}). 

148 

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

150 ''' 

151 if tx: 

152 self.tx = tx 

153 if ty: 

154 self.ty = ty 

155 if tz: 

156 self.tz = tz 

157 if sx: # secs to rads 

158 self.rx, self.sx = _r_s2(sx) 

159 if sy: 

160 self.ry, self.sy = _r_s2(sy) 

161 if sz: 

162 self.rz, self.sz = _r_s2(sz) 

163 if s: 

164 self.s = s 

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

166 

167 self._register(Transforms, name) 

168 

169 def __eq__(self, other): 

170 '''Compare this and an other transform. 

171 

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

173 

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

175 ''' 

176 return self is other or (isinstance(other, Transform) and 

177 self.tx == other.tx and 

178 self.ty == other.ty and 

179 self.tz == other.tz and 

180 self.rx == other.rx and 

181 self.ry == other.ry and 

182 self.rz == other.rz and 

183 self.s == other.s) 

184 

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

186 '''Helmert-transform cartesian B{C{other}}. 

187 

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

189 ''' 

190 try: # only CartesianBase 

191 return other.toTransform(self) 

192 except AttributeError: 

193 pass 

194 raise _IsnotError(_cartesian_, other=other) 

195 

196 def inverse(self, name=NN): 

197 '''Return the inverse of this transform. 

198 

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

200 

201 @return: Inverse (Transform). 

202 

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

204 ''' 

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

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

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

208 

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

210 '''Return this transform as a string. 

211 

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

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

214 this transform's name. 

215 

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

217 ''' 

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

219 'rx', 'ry', 'rz', _s_, 's1', 

220 _sx_, _sy_, _sz_) 

221 

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

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

224 

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

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

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

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

229 

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

231 ''' 

232 xyz1 = x, y, z, _1_0 

233 s1 = self.s1 

234 if inverse: 

235 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

243 name=self.name) 

244 

245 

246class Transforms(_NamedEnum): 

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

248 to accommodate the L{_LazyNamedEnumItem} properties. 

249 ''' 

250 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

252 ''' 

253 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

254 

255Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

258Transforms._assert( 

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

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

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

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

263 s=_F( 1.2727)), 

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

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

266 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

276 Identity = _lazy(_Identity_), 

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

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

279 s=_F( -8.15)), 

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

281 # XXX rotation signs may be opposite, to be checked 

282 sx=_F( -1.042), sy=_F( -0.214), sz=_F( -0.631), 

283 s=_F( -1.1)), 

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

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

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

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

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

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

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

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

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

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

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

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

296 s=_F(-0.0015)), 

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

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

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

300 s=_F( 20.4894)), 

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

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

303 WGS84 = _lazy(_WGS84_), # unity 

304) 

305 

306 

307class Datum(_NamedEnumItem): 

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

309 ''' 

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

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

312 

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

314 '''New L{Datum}. 

315 

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

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

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

319 

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

321 

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

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

324 not a L{Transform}. 

325 ''' 

326 self._ellipsoid = ellipsoid or Datum._ellipsoid 

327 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

328 

329 self._transform = transform or Datum._transform 

330 _xinstanceof(Transform, transform=self.transform) 

331 

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

333 

334 def __eq__(self, other): 

335 '''Compare this and an other datum. 

336 

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

338 

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

340 ''' 

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

342 self.ellipsoid == other.ellipsoid and 

343 self.transform == other.transform) 

344 

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

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

347 

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

349 ''' 

350 try: # only CartesianBase and EllipsoidalLatLonBase 

351 return other.toDatum(self) 

352 except AttributeError: 

353 pass 

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

355 

356 def ecef(self, Ecef=None): 

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

358 

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

360 

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

362 

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

364 

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

366 ''' 

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

368 

369 @Property_RO 

370 def ellipsoid(self): 

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

372 ''' 

373 return self._ellipsoid 

374 

375 @Property_RO 

376 def exactTM(self): 

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

378 ''' 

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

380 

381 @Property_RO 

382 def isEllipsoidal(self): 

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

384 ''' 

385 return self._ellipsoid.isEllipsoidal 

386 

387 @Property_RO 

388 def isOblate(self): 

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

390 ''' 

391 return self._ellipsoid.isOblate 

392 

393 @Property_RO 

394 def isProlate(self): 

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

396 ''' 

397 return self._ellipsoid.isProlate 

398 

399 @Property_RO 

400 def isSpherical(self): 

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

402 ''' 

403 return self._ellipsoid.isSpherical 

404 

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

406 '''Return this datum as a string. 

407 

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

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

410 this datum's name. 

411 

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

413 ''' 

414 t = [] if name is None else \ 

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

416 for a in (_ellipsoid_, _transform_): 

417 v = getattr(self, a) 

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

419 return sep.join(t) 

420 

421 @Property_RO 

422 def transform(self): 

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

424 ''' 

425 return self._transform 

426 

427 

428def _En2(earth, name): 

429 '''(INTERNAL) Helper for C{_ellipsoid} and C{_ellipsoidal_datum}. 

430 ''' 

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

432 E = earth 

433 n = _UNDER(name or E.name) 

434 elif isinstance(earth, Datum): 

435 E = earth.ellipsoid 

436 n = _UNDER(name or earth.name) 

437 elif isinstance(earth, a_f2Tuple): 

438 n = _UNDER(name or earth.name) 

439 E = Ellipsoid(earth.a, earth.b, name=n) 

440 elif islistuple(earth, minum=2): 

441 a, f = earth[:2] 

442 n = _UNDER(name or _xattr(earth, name=NN)) 

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

444 else: 

445 E, n = None, NN 

446 return E, n 

447 

448 

449def _a_ellipsoid(a_ellipsoid, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ... 

450 '''(INTERNAL) Get an ellipsoid from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}, 

451 an L{Ellipsoid} or L{Ellipsoid2} from L{Datum} or C{a_f2Tuple}. 

452 ''' 

453 if f is None: 

454 E, _ = _En2(a_ellipsoid, name) 

455 if raiser and not E: 

456 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_ellipsoid}) 

457 else: 

458 E = Ellipsoid2(a_ellipsoid, f, name=name) 

459 return E 

460 

461 

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

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

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

465 

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

467 ''' 

468 if isinstance(earth, Datum): 

469 d = earth 

470 else: 

471 E, n = _En2(earth, name) 

472 if not E: 

473 n = raiser or _earth_ 

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

475 d = Datum(E, transform=Transforms.Identity, name=n) 

476 if raiser and not d.isEllipsoidal: 

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

478 return d 

479 

480 

481def _mean_radius(radius, *lats): 

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

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

484 ''' 

485 if radius is R_M: 

486 r = radius 

487 elif isscalar(radius): 

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

489 else: 

490 E = _ellipsoidal_datum(radius).ellipsoid 

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

492 return r 

493 

494 

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

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

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

498 

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

500 ''' 

501 if isscalar(earth): 

502 E = Datums.Sphere.ellipsoid 

503 if earth == E.a == E.b and not name: 

504 d = Datums.Sphere 

505 else: 

506 r = Radius_(earth, Error=Error) # invalid datum 

507 n = _UNDER(name) 

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

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

510 else: 

511 d = _ellipsoidal_datum(earth, Error=Error, name=name) 

512 if raiser and not d.isSpherical: 

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

514 return d 

515 

516 

517class Datums(_NamedEnum): 

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

519 to accommodate the L{_LazyNamedEnumItem} properties. 

520 ''' 

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

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

523 ''' 

524 return Datum(Ellipsoids.get(ellipsoid_name), 

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

526 

527Datums = Datums(Datum) # PYCHOK singleton 

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

529# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

533Datums._assert( 

534 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

538 

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

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

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

542 

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

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

545 

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

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

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

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

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

551 

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

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

554 

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

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

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

558 

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

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

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

562 

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

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

565 

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

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

568 

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

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

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

572 

573 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

575 

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

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

578 

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

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

581 

582 # XXX psuedo-ellipsoids for spherical LatLon 

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

584 

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

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

587 

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

589 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

590 

591 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

592) 

593 

594_WGS84 = Datums.WGS84 # PYCHOK exported internally 

595 

596if __name__ == '__main__': 

597 

598 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

599 from pygeodesy.lazily import printf 

600 

601 # __doc__ of this file, force all into registery 

602 for r in (Datums, Transforms): 

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

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

605 

606# **) MIT License 

607# 

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

609# 

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

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

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

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

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

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

616# 

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

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

619# 

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

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

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

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

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

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

626# OTHER DEALINGS IN THE SOFTWARE.