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.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

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 __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

245 class 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

254 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

257 Transforms._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

306 class 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

427 def _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

448 def _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

461 def _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

480 def _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

494 def _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

516 class 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

526 Datums = 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>.

532 Datums._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

595 if __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.