1

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

3

4 u'''Datums and transformations thereof.

5

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

7

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

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

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

11 I{(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

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

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

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

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

18

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

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

21

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

23 OSGB36 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

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

66

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

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

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

70 _EWGS84, Vector3Tuple

71 from pygeodesy.errors import _IsnotError, _xattr

72 from pygeodesy.fmath import fdot, fmean, Fmt

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

74 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

75 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \

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

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

78 _WGS72_, _WGS84_, _under

79 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

80 from pygeodesy.named import _NamedEnum, _NamedEnumItem, \

81 _lazyNamedEnumItem as _lazy, Property_RO

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

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

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

85 from pygeodesy.units import radians, Radius_

86

87# from math import radians # from .units

88

89__all__ = _ALL_LAZY.datums

90__version__ = '23.08.20'

91

92_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

93_BD72_ = 'BD72'

94_DHDN_ = 'DHDN'

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

110 def _r_s2(s_):

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

112 '''

113 return _F(_r_s1 * s_), s_

114

115

116 class 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)

177 and self.tx == other.tx

178 and self.ty == other.ty

179 and self.tz == other.tz

180 and self.rx == other.rx

181 and self.ry == other.ry

182 and self.rz == other.rz

183 and self.s == other.s)

184

185 def __hash__(self):

186 return self._hash # memoized

187

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

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

190

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

192 '''

193 try: # only CartesianBase

194 return other.toTransform(self)

195 except AttributeError:

196 pass

197 raise _IsnotError(_cartesian_, other=other)

198

199 @Property_RO

200 def _hash(self):

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

202 self.tx, self.ty, self.tz))

203

204 def inverse(self, name=NN):

205 '''Return the inverse of this transform.

206

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

208

209 @return: Inverse (Transform).

210

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

212 '''

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

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

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

216

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

218 '''Return this transform as a string.

219

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

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

222 this transform's name.

223

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

225 '''

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

227 'rx', 'ry', 'rz', _s_, 's1',

228 _sx_, _sy_, _sz_)

229

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

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

232

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

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

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

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

237

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

239 '''

240 xyz1 = x, y, z, _1_0

241 s1 = self.s1

242 if inverse:

243 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

251 name=self.name)

252

253

254 class Transforms(_NamedEnum):

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

256 to accommodate the L{_LazyNamedEnumItem} properties.

257 '''

258 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

260 '''

261 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

262

263 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

266 Transforms._assert(

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

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

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

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

271 s=_F( 1.2727)),

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

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

274 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

284 Identity = _lazy(_Identity_),

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

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

287 s=_F( -8.15)),

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

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

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

291 s=_F( -1.1)),

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

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

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

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

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

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

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

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

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

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

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

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

304 s=_F(-0.0015)),

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

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

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

308 s=_F( 20.4894)),

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

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

311 WGS84 = _lazy(_WGS84_), # unity

312)

313

314

315 class Datum(_NamedEnumItem):

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

317 '''

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

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

320

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

322 '''New L{Datum}.

323

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

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

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

327

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

329

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

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

332 not a L{Transform}.

333 '''

334 self._ellipsoid = ellipsoid or Datum._ellipsoid

335 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

336

337 self._transform = transform or Datum._transform

338 _xinstanceof(Transform, transform=self.transform)

339

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

341

342 def __eq__(self, other):

343 '''Compare this and an other datum.

344

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

346

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

348 '''

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

350 self.ellipsoid == other.ellipsoid and

351 self.transform == other.transform)

352

353 def __hash__(self):

354 return self._hash # memoized

355

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

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

358

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

360 '''

361 try: # only CartesianBase and EllipsoidalLatLonBase

362 return other.toDatum(self)

363 except AttributeError:

364 pass

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

366

367 def ecef(self, Ecef=None):

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

369

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

371

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

373

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

375

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

377 '''

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

379

380 @Property_RO

381 def ellipsoid(self):

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

383 '''

384 return self._ellipsoid

385

386 @Property_RO

387 def exactTM(self):

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

389 '''

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

391

392 @Property_RO

393 def _hash(self):

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

395

396 @Property_RO

397 def isEllipsoidal(self):

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

399 '''

400 return self.ellipsoid.isEllipsoidal

401

402 @Property_RO

403 def isOblate(self):

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

405 '''

406 return self.ellipsoid.isOblate

407

408 @Property_RO

409 def isProlate(self):

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

411 '''

412 return self.ellipsoid.isProlate

413

414 @Property_RO

415 def isSpherical(self):

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

417 '''

418 return self.ellipsoid.isSpherical

419

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

421 '''Return this datum as a string.

422

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

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

425 this datum's name.

426

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

428 '''

429 t = [] if name is None else \

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

431 for a in (_ellipsoid_, _transform_):

432 v = getattr(self, a)

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

434 return sep.join(t)

435

436 @Property_RO

437 def transform(self):

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

439 '''

440 return self._transform

441

442

443 def _En2(earth, name):

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

445 '''

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

447 E = earth

448 n = _under(name or E.name)

449 elif isinstance(earth, Datum):

450 E = earth.ellipsoid

451 n = _under(name or earth.name)

452 elif isinstance(earth, a_f2Tuple):

453 n = _under(name or earth.name)

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

455 elif islistuple(earth, minum=2):

456 a, f = earth[:2]

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

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

459 else:

460 E, n = None, NN

461 return E, n

462

463

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

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

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

467 '''

468 if f is None:

469 E, _ = _En2(a_ellipsoid, name)

470 if raiser and not E:

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

472 else:

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

474 return E

475

476

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

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

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

480

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

482 '''

483 if isinstance(earth, Datum):

484 d = earth

485 else:

486 E, n = _En2(earth, name)

487 if not E:

488 n = raiser or _earth_

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

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

491 if raiser and not d.isEllipsoidal:

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

493 return d

494

495

496 def _mean_radius(radius, *lats):

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

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

499 '''

500 if radius is R_M:

501 r = radius

502 elif isscalar(radius):

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

504 else:

505 E = _ellipsoidal_datum(radius).ellipsoid

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

507 return r

508

509

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

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

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

513

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

515 '''

516 if isscalar(earth):

517 E = Datums.Sphere.ellipsoid

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

519 d = Datums.Sphere

520 else:

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

522 n = _under(name)

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

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

525 else:

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

527 if raiser and not d.isSpherical:

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

529 return d

530

531

532 class Datums(_NamedEnum):

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

534 to accommodate the L{_LazyNamedEnumItem} properties.

535 '''

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

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

538 '''

539 return Datum(Ellipsoids.get(ellipsoid_name),

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

541

542 Datums = Datums(Datum) # PYCHOK singleton

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

544# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

548 Datums._assert(

549 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

553

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

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

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

557

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

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

560

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

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

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

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

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

566

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

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

569

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

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

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

573

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

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

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

577

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

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

580

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

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

583

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

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

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

587

588 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

590

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

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

593

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

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

596

597 # XXX psuedo-ellipsoids for spherical LatLon

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

599

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

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

602

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

604 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

605

606 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

607)

608

609_WGS84 = Datums.WGS84

610 assert _WGS84.ellipsoid is _EWGS84

611

612 if __name__ == '__main__':

613

614 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

615 from pygeodesy.lazily import printf

616

617 # __doc__ of this file, force all into registery

618 for r in (Datums, Transforms):

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

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

621

622# **) MIT License

623#

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

625#

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

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

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

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

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

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

632#

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

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

635#

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

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

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

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

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

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

642# OTHER DEALINGS IN THE SOFTWARE.