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, neg_, _xinstanceof

68 from pygeodesy.constants import _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 # from .fmath

71 from pygeodesy.fmath import fdot, fmean, Fmt, _IsnotError

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

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

74 _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \

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

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

77 _UNDER_, _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.03.29'

90

91_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

92_BD72_ = 'BD72'

93_DHDN_ = 'DHDN'

94_earth_ = 'earth'

95_ED50_ = 'ED50'

96_GDA2020_ = 'GDA2020'

97_Identity_ = 'Identity'

98_Inverse_ = 'Inverse'

99_Irl1965_ = 'Irl1965'

100_Irl1975_ = 'Irl1975'

101_MGI_ = 'MGI'

102_NTF_ = 'NTF'

103_OSGB36_ = 'OSGB36'

104_Potsdam_ = 'Potsdam'

105_TokyoJapan_ = 'TokyoJapan'

106

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

108

109

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) and

177 self.tx == other.tx and

178 self.ty == other.ty and

179 self.tz == other.tz and

180 self.rx == other.rx and

181 self.ry == other.ry and

182 self.rz == other.rz and

183 self.s == other.s)

184

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

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

187

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

189 '''

190 try: # only CartesianBase

191 return other.toTransform(self)

192 except AttributeError:

193 pass

194 raise _IsnotError(_cartesian_, other=other)

195

196 def inverse(self, name=NN):

197 '''Return the inverse of this transform.

198

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

200

201 @return: Inverse (Transform).

202

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

204 '''

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

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

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

208

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

210 '''Return this transform as a string.

211

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

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

214 this transform's name.

215

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

217 '''

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

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

220 _sx_, _sy_, _sz_)

221

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

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

224

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

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

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

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

229

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

231 '''

232 xyz1 = x, y, z, _1_0

233 s1 = self.s1

234 if inverse:

235 xyz1 = neg_(*xyz1)

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

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

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

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

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

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

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

243 name=self.name)

244

245

246 class Transforms(_NamedEnum):

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

248 to accommodate the L{_LazyNamedEnumItem} properties.

249 '''

250 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

252 '''

253 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

254

255 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

258 Transforms._assert(

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

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

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

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

263 s=_F( 1.2727)),

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

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

266 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

276 Identity = _lazy(_Identity_),

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

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

279 s=_F( -8.15)),

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

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

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

283 s=_F( -1.1)),

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

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

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

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

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

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

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

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

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

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

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

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

296 s=_F(-0.0015)),

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

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

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

300 s=_F( 20.4894)),

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

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

303 WGS84 = _lazy(_WGS84_), # unity

304)

305

306

307 class Datum(_NamedEnumItem):

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

309 '''

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

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

312

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

314 '''New L{Datum}.

315

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

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

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

319

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

321

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

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

324 not a L{Transform}.

325 '''

326 self._ellipsoid = ellipsoid or Datum._ellipsoid

327 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

328

329 self._transform = transform or Datum._transform

330 _xinstanceof(Transform, transform=self.transform)

331

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

333

334 def __eq__(self, other):

335 '''Compare this and an other datum.

336

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

338

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

340 '''

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

342 self.ellipsoid == other.ellipsoid and

343 self.transform == other.transform)

344

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

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

347

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

349 '''

350 try: # only CartesianBase and EllipsoidalLatLonBase

351 return other.toDatum(self)

352 except AttributeError:

353 pass

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

355

356 def ecef(self, Ecef=None):

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

358

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

360

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

362

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

364

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

366 '''

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

368

369 @Property_RO

370 def ellipsoid(self):

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

372 '''

373 return self._ellipsoid

374

375 @Property_RO

376 def exactTM(self):

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

378 '''

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

380

381 @Property_RO

382 def isEllipsoidal(self):

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

384 '''

385 return self._ellipsoid.isEllipsoidal

386

387 @Property_RO

388 def isOblate(self):

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

390 '''

391 return self._ellipsoid.isOblate

392

393 @Property_RO

394 def isProlate(self):

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

396 '''

397 return self._ellipsoid.isProlate

398

399 @Property_RO

400 def isSpherical(self):

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

402 '''

403 return self._ellipsoid.isSpherical

404

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

406 '''Return this datum as a string.

407

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

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

410 this datum's name.

411

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

413 '''

414 t = [] if name is None else \

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

416 for a in (_ellipsoid_, _transform_):

417 v = getattr(self, a)

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

419 return sep.join(t)

420

421 @Property_RO

422 def transform(self):

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

424 '''

425 return self._transform

426

427

428 def _En2(earth, name):

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

430 '''

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

432 E = earth

433 n = _UNDER(name or E.name)

434 elif isinstance(earth, Datum):

435 E = earth.ellipsoid

436 n = _UNDER(name or earth.name)

437 elif isinstance(earth, a_f2Tuple):

438 n = _UNDER(name or earth.name)

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

440 elif islistuple(earth, minum=2):

441 a, f = earth[:2]

442 n = _UNDER(name or getattr(earth, _name_, NN))

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

444 else:

445 E, n = None, NN

446 return E, n

447

448

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

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

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

452 '''

453 if f is None:

454 E, _ = _En2(a_ellipsoid, name)

455 if raiser and not E:

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

457 else:

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

459 return E

460

461

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

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

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

465

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

467 '''

468 if isinstance(earth, Datum):

469 d = earth

470 else:

471 E, n = _En2(earth, name)

472 if not E:

473 n = raiser or _earth_

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

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

476 if raiser and not d.isEllipsoidal:

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

478 return d

479

480

481 def _mean_radius(radius, *lats):

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

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

484 '''

485 if isscalar(radius):

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

487 else:

488 E = _ellipsoidal_datum(radius).ellipsoid

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

490 return r

491

492

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

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

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

496

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

498 '''

499 if isscalar(earth):

500 E = Datums.Sphere.ellipsoid

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

502 d = Datums.Sphere

503 else:

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

505 n = _UNDER(name)

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

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

508 else:

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

510 if raiser and not d.isSpherical:

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

512 return d

513

514

515 class Datums(_NamedEnum):

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

517 to accommodate the L{_LazyNamedEnumItem} properties.

518 '''

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

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

521 '''

522 return Datum(Ellipsoids.get(ellipsoid_name),

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

524

525 Datums = Datums(Datum) # PYCHOK singleton

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

527# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

531 Datums._assert(

532 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

536

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

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

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

540

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

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

543

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

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

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

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

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

549

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

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

552

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

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

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

556

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

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

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

560

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

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

563

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

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

566

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

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

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

570

571 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

573

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

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

576

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

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

579

580 # XXX psuedo-ellipsoids for spherical LatLon

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

582

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

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

585

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

587 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

588

589 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

590)

591

592_WGS84 = Datums.WGS84 # PYCHOK exported internally

593

594 if __name__ == '__main__':

595

596 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

597

598 # __doc__ of this file, force all into registery

599 for r in (Datums, Transforms):

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

601 print(_NLATvar_.join(i.strip(_COMMA_) for i in t))

602

603# **) MIT License

604#

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

606#

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

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

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

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

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

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

613#

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

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

616#

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

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

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

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

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

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

623# OTHER DEALINGS IN THE SOFTWARE.