Coverage for pygeodesy/datums.py: 97%

184 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-07-12 13:40 -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_, _earth_, \ 

74 _ellipsoid_, _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.06.12' 

90 

91_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_) 

92_BD72_ = 'BD72' 

93_DHDN_ = 'DHDN' 

94_ED50_ = 'ED50' 

95_GDA2020_ = 'GDA2020' 

96_Identity_ = 'Identity' 

97_Inverse_ = 'Inverse' 

98_Irl1965_ = 'Irl1965' 

99_Irl1975_ = 'Irl1975' 

100_MGI_ = 'MGI' 

101_NTF_ = 'NTF' 

102_OSGB36_ = 'OSGB36' 

103_Potsdam_ = 'Potsdam' 

104_TokyoJapan_ = 'TokyoJapan' 

105 

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

107 

108 

109def _r_s2(s_): 

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

111 ''' 

112 return _F(_r_s1 * s_), s_ 

113 

114 

115class Transform(_NamedEnumItem): 

116 '''Helmert transformation. 

117 

118 @see: L{Helmert7Tuple}. 

119 ''' 

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

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

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

123 

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

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

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

127 

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

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

130 

131 sx = _0_0 # x rotation (degree seconds) 

132 sy = _0_0 # y rotation (degree seconds) 

133 sz = _0_0 # z rotation (degree seconds) 

134 

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

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

137 '''New L{Transform}. 

138 

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

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

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

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

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

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

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

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

147 

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

149 ''' 

150 if tx: 

151 self.tx = tx 

152 if ty: 

153 self.ty = ty 

154 if tz: 

155 self.tz = tz 

156 if sx: # secs to rads 

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

158 if sy: 

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

160 if sz: 

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

162 if s: 

163 self.s = s 

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

165 

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) 

176 and self.tx == other.tx 

177 and self.ty == other.ty 

178 and self.tz == other.tz 

179 and self.rx == other.rx 

180 and self.ry == other.ry 

181 and self.rz == other.rz 

182 and self.s == other.s) 

183 

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

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

186 

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

188 ''' 

189 try: # only CartesianBase 

190 return other.toTransform(self) 

191 except AttributeError: 

192 pass 

193 raise _IsnotError(_cartesian_, other=other) 

194 

195 def inverse(self, name=NN): 

196 '''Return the inverse of this transform. 

197 

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

199 

200 @return: Inverse (Transform). 

201 

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

203 ''' 

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

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

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

207 

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

209 '''Return this transform as a string. 

210 

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

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

213 this transform's name. 

214 

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

216 ''' 

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

218 'rx', 'ry', 'rz', _s_, 's1', 

219 _sx_, _sy_, _sz_) 

220 

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

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

223 

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

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

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

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

228 

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

230 ''' 

231 xyz1 = x, y, z, _1_0 

232 s1 = self.s1 

233 if inverse: 

234 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

242 name=self.name) 

243 

244 

245class Transforms(_NamedEnum): 

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

247 to accommodate the L{_LazyNamedEnumItem} properties. 

248 ''' 

249 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

251 ''' 

252 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

253 

254Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

257Transforms._assert( 

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

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

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

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

262 s=_F( 1.2727)), 

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

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

265 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

275 Identity = _lazy(_Identity_), 

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

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

278 s=_F( -8.15)), 

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

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

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

282 s=_F( -1.1)), 

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

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

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

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

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

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

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

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

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

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

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

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

295 s=_F(-0.0015)), 

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

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

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

299 s=_F( 20.4894)), 

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

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

302 WGS84 = _lazy(_WGS84_), # unity 

303) 

304 

305 

306class Datum(_NamedEnumItem): 

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

308 ''' 

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

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

311 

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

313 '''New L{Datum}. 

314 

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

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

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

318 

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

320 

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

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

323 not a L{Transform}. 

324 ''' 

325 self._ellipsoid = ellipsoid or Datum._ellipsoid 

326 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

327 

328 self._transform = transform or Datum._transform 

329 _xinstanceof(Transform, transform=self.transform) 

330 

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

332 

333 def __eq__(self, other): 

334 '''Compare this and an other datum. 

335 

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

337 

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

339 ''' 

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

341 self.ellipsoid == other.ellipsoid and 

342 self.transform == other.transform) 

343 

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

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

346 

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

348 ''' 

349 try: # only CartesianBase and EllipsoidalLatLonBase 

350 return other.toDatum(self) 

351 except AttributeError: 

352 pass 

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

354 

355 def ecef(self, Ecef=None): 

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

357 

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

359 

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

361 

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

363 

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

365 ''' 

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

367 

368 @Property_RO 

369 def ellipsoid(self): 

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

371 ''' 

372 return self._ellipsoid 

373 

374 @Property_RO 

375 def exactTM(self): 

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

377 ''' 

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

379 

380 @Property_RO 

381 def isEllipsoidal(self): 

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

383 ''' 

384 return self._ellipsoid.isEllipsoidal 

385 

386 @Property_RO 

387 def isOblate(self): 

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

389 ''' 

390 return self._ellipsoid.isOblate 

391 

392 @Property_RO 

393 def isProlate(self): 

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

395 ''' 

396 return self._ellipsoid.isProlate 

397 

398 @Property_RO 

399 def isSpherical(self): 

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

401 ''' 

402 return self._ellipsoid.isSpherical 

403 

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

405 '''Return this datum as a string. 

406 

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

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

409 this datum's name. 

410 

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

412 ''' 

413 t = [] if name is None else \ 

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

415 for a in (_ellipsoid_, _transform_): 

416 v = getattr(self, a) 

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

418 return sep.join(t) 

419 

420 @Property_RO 

421 def transform(self): 

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

423 ''' 

424 return self._transform 

425 

426 

427def _En2(earth, name): 

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

429 ''' 

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

431 E = earth 

432 n = _UNDER(name or E.name) 

433 elif isinstance(earth, Datum): 

434 E = earth.ellipsoid 

435 n = _UNDER(name or earth.name) 

436 elif isinstance(earth, a_f2Tuple): 

437 n = _UNDER(name or earth.name) 

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

439 elif islistuple(earth, minum=2): 

440 a, f = earth[:2] 

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

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

443 else: 

444 E, n = None, NN 

445 return E, n 

446 

447 

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

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

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

451 ''' 

452 if f is None: 

453 E, _ = _En2(a_ellipsoid, name) 

454 if raiser and not E: 

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

456 else: 

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

458 return E 

459 

460 

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

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

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

464 

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

466 ''' 

467 if isinstance(earth, Datum): 

468 d = earth 

469 else: 

470 E, n = _En2(earth, name) 

471 if not E: 

472 n = raiser or _earth_ 

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

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

475 if raiser and not d.isEllipsoidal: 

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

477 return d 

478 

479 

480def _mean_radius(radius, *lats): 

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

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

483 ''' 

484 if radius is R_M: 

485 r = radius 

486 elif isscalar(radius): 

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

488 else: 

489 E = _ellipsoidal_datum(radius).ellipsoid 

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

491 return r 

492 

493 

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

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

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

497 

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

499 ''' 

500 if isscalar(earth): 

501 E = Datums.Sphere.ellipsoid 

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

503 d = Datums.Sphere 

504 else: 

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

506 n = _UNDER(name) 

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

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

509 else: 

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

511 if raiser and not d.isSpherical: 

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

513 return d 

514 

515 

516class Datums(_NamedEnum): 

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

518 to accommodate the L{_LazyNamedEnumItem} properties. 

519 ''' 

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

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

522 ''' 

523 return Datum(Ellipsoids.get(ellipsoid_name), 

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

525 

526Datums = Datums(Datum) # PYCHOK singleton 

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

528# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

532Datums._assert( 

533 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

537 

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

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

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

541 

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

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

544 

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

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

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

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

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

550 

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

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

553 

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

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

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

557 

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

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

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

561 

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

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

564 

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

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

567 

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

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

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

571 

572 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

574 

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

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

577 

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

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

580 

581 # XXX psuedo-ellipsoids for spherical LatLon 

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

583 

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

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

586 

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

588 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

589 

590 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

591) 

592 

593_WGS84 = Datums.WGS84 # PYCHOK exported internally 

594 

595if __name__ == '__main__': 

596 

597 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

598 from pygeodesy.lazily import printf 

599 

600 # __doc__ of this file, force all into registery 

601 for r in (Datums, Transforms): 

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

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

604 

605# **) MIT License 

606# 

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

608# 

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

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

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

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

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

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

615# 

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

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

618# 

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

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

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

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

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

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

625# OTHER DEALINGS IN THE SOFTWARE.