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, Vector3Tuple

70 from pygeodesy.errors import _IsnotError, _xattr

71 from pygeodesy.fmath import fdot, fmean, Fmt

72 from 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

78 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

79 from 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

84 from pygeodesy.units import radians, Radius_

85

86# from math import radians # from .units

87

88__all__ = _ALL_LAZY.datums

89__version__ = '23.08.05'

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

109 def _r_s2(s_):

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

111 '''

112 return _F(_r_s1 * s_), s_

113

114

115 class 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 __hash__(self):

185 return self._hash # memoized

186

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

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

189

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

191 '''

192 try: # only CartesianBase

193 return other.toTransform(self)

194 except AttributeError:

195 pass

196 raise _IsnotError(_cartesian_, other=other)

197

198 @Property_RO

199 def _hash(self):

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

201 self.tx, self.ty, self.tz))

202

203 def inverse(self, name=NN):

204 '''Return the inverse of this transform.

205

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

207

208 @return: Inverse (Transform).

209

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

211 '''

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

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

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

215

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

217 '''Return this transform as a string.

218

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

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

221 this transform's name.

222

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

224 '''

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

226 'rx', 'ry', 'rz', _s_, 's1',

227 _sx_, _sy_, _sz_)

228

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

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

231

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

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

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

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

236

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

238 '''

239 xyz1 = x, y, z, _1_0

240 s1 = self.s1

241 if inverse:

242 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

250 name=self.name)

251

252

253 class Transforms(_NamedEnum):

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

255 to accommodate the L{_LazyNamedEnumItem} properties.

256 '''

257 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

259 '''

260 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

261

262 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

265 Transforms._assert(

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

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

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

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

270 s=_F( 1.2727)),

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

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

273 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

283 Identity = _lazy(_Identity_),

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

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

286 s=_F( -8.15)),

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

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

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

290 s=_F( -1.1)),

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

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

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

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

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

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

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

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

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

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

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

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

303 s=_F(-0.0015)),

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

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

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

307 s=_F( 20.4894)),

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

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

310 WGS84 = _lazy(_WGS84_), # unity

311)

312

313

314 class Datum(_NamedEnumItem):

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

316 '''

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

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

319

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

321 '''New L{Datum}.

322

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

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

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

326

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

328

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

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

331 not a L{Transform}.

332 '''

333 self._ellipsoid = ellipsoid or Datum._ellipsoid

334 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

335

336 self._transform = transform or Datum._transform

337 _xinstanceof(Transform, transform=self.transform)

338

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

340

341 def __eq__(self, other):

342 '''Compare this and an other datum.

343

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

345

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

347 '''

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

349 self.ellipsoid == other.ellipsoid and

350 self.transform == other.transform)

351

352 def __hash__(self):

353 return self._hash # memoized

354

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

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

357

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

359 '''

360 try: # only CartesianBase and EllipsoidalLatLonBase

361 return other.toDatum(self)

362 except AttributeError:

363 pass

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

365

366 def ecef(self, Ecef=None):

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

368

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

370

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

372

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

374

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

376 '''

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

378

379 @Property_RO

380 def ellipsoid(self):

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

382 '''

383 return self._ellipsoid

384

385 @Property_RO

386 def exactTM(self):

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

388 '''

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

390

391 @Property_RO

392 def _hash(self):

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

394

395 @Property_RO

396 def isEllipsoidal(self):

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

398 '''

399 return self.ellipsoid.isEllipsoidal

400

401 @Property_RO

402 def isOblate(self):

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

404 '''

405 return self.ellipsoid.isOblate

406

407 @Property_RO

408 def isProlate(self):

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

410 '''

411 return self.ellipsoid.isProlate

412

413 @Property_RO

414 def isSpherical(self):

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

416 '''

417 return self.ellipsoid.isSpherical

418

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

420 '''Return this datum as a string.

421

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

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

424 this datum's name.

425

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

427 '''

428 t = [] if name is None else \

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

430 for a in (_ellipsoid_, _transform_):

431 v = getattr(self, a)

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

433 return sep.join(t)

434

435 @Property_RO

436 def transform(self):

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

438 '''

439 return self._transform

440

441

442 def _En2(earth, name):

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

444 '''

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

446 E = earth

447 n = _under(name or E.name)

448 elif isinstance(earth, Datum):

449 E = earth.ellipsoid

450 n = _under(name or earth.name)

451 elif isinstance(earth, a_f2Tuple):

452 n = _under(name or earth.name)

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

454 elif islistuple(earth, minum=2):

455 a, f = earth[:2]

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

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

458 else:

459 E, n = None, NN

460 return E, n

461

462

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

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

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

466 '''

467 if f is None:

468 E, _ = _En2(a_ellipsoid, name)

469 if raiser and not E:

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

471 else:

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

473 return E

474

475

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

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

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

479

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

481 '''

482 if isinstance(earth, Datum):

483 d = earth

484 else:

485 E, n = _En2(earth, name)

486 if not E:

487 n = raiser or _earth_

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

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

490 if raiser and not d.isEllipsoidal:

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

492 return d

493

494

495 def _mean_radius(radius, *lats):

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

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

498 '''

499 if radius is R_M:

500 r = radius

501 elif isscalar(radius):

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

503 else:

504 E = _ellipsoidal_datum(radius).ellipsoid

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

506 return r

507

508

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

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

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

512

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

514 '''

515 if isscalar(earth):

516 E = Datums.Sphere.ellipsoid

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

518 d = Datums.Sphere

519 else:

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

521 n = _under(name)

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

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

524 else:

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

526 if raiser and not d.isSpherical:

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

528 return d

529

530

531 class Datums(_NamedEnum):

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

533 to accommodate the L{_LazyNamedEnumItem} properties.

534 '''

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

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

537 '''

538 return Datum(Ellipsoids.get(ellipsoid_name),

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

540

541 Datums = Datums(Datum) # PYCHOK singleton

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

543# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

547 Datums._assert(

548 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

552

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

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

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

556

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

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

559

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

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

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

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

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

565

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

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

568

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

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

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

572

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

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

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

576

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

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

579

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

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

582

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

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

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

586

587 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

589

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

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

592

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

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

595

596 # XXX psuedo-ellipsoids for spherical LatLon

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

598

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

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

601

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

603 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

604

605 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

606)

607

608_WGS84 = Datums.WGS84 # PYCHOK exported internally

609

610 if __name__ == '__main__':

611

612 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

613 from pygeodesy.lazily import printf

614

615 # __doc__ of this file, force all into registery

616 for r in (Datums, Transforms):

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

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

619

620# **) MIT License

621#

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

623#

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

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

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

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

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

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

630#

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

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

633#

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

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

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

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

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

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

640# OTHER DEALINGS IN THE SOFTWARE.