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

15 or below the earth’s surface. Latitude is measured in degrees from the equator, lomgitude from

16 the International Reference Meridian and height in meters above an ellipsoid based on the given

17 datum. The datum in turn is based on a reference ellipsoid and tied to geodetic survey

18 reference points.

19

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

21 previously various other reference ellipsoids and datum references were used.

22

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

24 datum, q.v. U{"A Guide to Coordinate Systems in Great Britain", Section 6

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

26

27@var Datums.BD72: Datum(name='BD72', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.BD72)

28@var Datums.DHDN: Datum(name='DHDN', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.DHDN)

29@var Datums.ED50: Datum(name='ED50', ellipsoid=Ellipsoids.Intl1924, transform=Transforms.ED50)

30@var Datums.GDA2020: Datum(name='GDA2020', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)

31@var Datums.GRS80: Datum(name='GRS80', ellipsoid=Ellipsoids.GRS80, transform=Transforms.WGS84)

32@var Datums.Irl1975: Datum(name='Irl1975', ellipsoid=Ellipsoids.AiryModified, transform=Transforms.Irl1975)

33@var Datums.Krassovski1940: Datum(name='Krassovski1940', ellipsoid=Ellipsoids.Krassovski1940, transform=Transforms.Krassovski1940)

34@var Datums.Krassowsky1940: Datum(name='Krassowsky1940', ellipsoid=Ellipsoids.Krassowsky1940, transform=Transforms.Krassowsky1940)

35@var Datums.MGI: Datum(name='MGI', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.MGI)

36@var Datums.NAD27: Datum(name='NAD27', ellipsoid=Ellipsoids.Clarke1866, transform=Transforms.NAD27)

37@var Datums.NAD83: Datum(name='NAD83', ellipsoid=Ellipsoids.GRS80, transform=Transforms.NAD83)

38@var Datums.NTF: Datum(name='NTF', ellipsoid=Ellipsoids.Clarke1880IGN, transform=Transforms.NTF)

39@var Datums.OSGB36: Datum(name='OSGB36', ellipsoid=Ellipsoids.Airy1830, transform=Transforms.OSGB36)

40@var Datums.Potsdam: Datum(name='Potsdam', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841)

41@var Datums.Sphere: Datum(name='Sphere', ellipsoid=Ellipsoids.Sphere, transform=Transforms.WGS84)

42@var Datums.TokyoJapan: Datum(name='TokyoJapan', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.TokyoJapan)

43@var Datums.WGS72: Datum(name='WGS72', ellipsoid=Ellipsoids.WGS72, transform=Transforms.WGS72)

44@var Datums.WGS84: Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84)

45

46@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)

47@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)

48@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)

49@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)

50@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)

51@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)

52@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)

53@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)

54@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)

55@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)

56@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)

57@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)

58@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)

59@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)

60@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)

61@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)

62@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)

63@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)

64'''

65# make sure int/int division yields float quotient, see .basics

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

67

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

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

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

71 _EWGS84, Vector3Tuple

72 from pygeodesy.errors import _IsnotError, _TypeError, _xattr

73 from pygeodesy.fmath import fdot, fmean, Fmt

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

75 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

76 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \

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

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

79 _WGS72_, _WGS84_, _under

80 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

81 from pygeodesy.named import _NamedEnum, _NamedEnumItem, Property_RO, property_RO, \

82 _lazyNamedEnumItem as _lazy

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

84# from pygeodesy.props import Property_RO, property_RO # from .named

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

86 from pygeodesy.units import radians, Radius_

87

88# from math import radians # from .units

89

90__all__ = _ALL_LAZY.datums

91__version__ = '23.11.05'

92

93_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

94_BD72_ = 'BD72'

95_DHDN_ = 'DHDN'

96_DHDNE_ = 'DHDNE'

97_DHDNW_ = 'DHDNW'

98_ED50_ = 'ED50'

99_GDA2020_ = 'GDA2020'

100_Identity_ = 'Identity'

101_Inverse_ = 'Inverse'

102_Irl1965_ = 'Irl1965'

103_Irl1975_ = 'Irl1975'

104_MGI_ = 'MGI'

105_NTF_ = 'NTF'

106_OSGB36_ = 'OSGB36'

107_Potsdam_ = 'Potsdam'

108_RPS = radians(_1_0 / _3600_0) # radians per degree_second

109_TokyoJapan_ = 'TokyoJapan'

110

111

112 def _rps2(s_):

113 '''(INTERNAL) Rotation in C{radians} and C{degree seconds}.

114 '''

115 return (_RPS * s_), s_

116

117

118 class Transform(_NamedEnumItem):

119 '''Helmert transformation.

120

121 @see: L{Helmert7Tuple}.

122 '''

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

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

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

126

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

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

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

130

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

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

133

134 sx = _0_0 # x rotation (degree seconds)

135 sy = _0_0 # y rotation (degree seconds)

136 sz = _0_0 # z rotation (degree seconds)

137

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

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

140 '''New L{Transform}.

141

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

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

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

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

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

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

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

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

150

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

152 '''

153 if tx:

154 self.tx = tx

155 if ty:

156 self.ty = ty

157 if tz:

158 self.tz = tz

159 if sx: # secs to rads

160 self.rx, self.sx = _rps2(sx)

161 if sy:

162 self.ry, self.sy = _rps2(sy)

163 if sz:

164 self.rz, self.sz = _rps2(sz)

165 if s:

166 self.s = s

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

168

169 self._register(Transforms, name)

170

171 def __eq__(self, other):

172 '''Compare this and an other transform.

173

174 @arg other: The other transform (L{Transform}).

175

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

177 '''

178 return self is other or (isinstance(other, Transform)

179 and self.tx == other.tx

180 and self.ty == other.ty

181 and self.tz == other.tz

182 and self.rx == other.rx

183 and self.ry == other.ry

184 and self.rz == other.rz

185 and self.s == other.s)

186

187 def __hash__(self):

188 return self._hash # memoized

189

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

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

192

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

194 '''

195 try: # only CartesianBase

196 return other.toTransform(self)

197 except AttributeError:

198 pass

199 raise _IsnotError(_cartesian_, other=other)

200

201 @Property_RO

202 def _hash(self):

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

204 self.tx, self.ty, self.tz))

205

206 def inverse(self, name=NN):

207 '''Return the inverse of this transform.

208

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

210

211 @return: Inverse (Transform).

212

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

214 '''

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

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

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

218

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

220 '''Return this transform as a string.

221

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

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

224 this transform's name.

225

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

227 '''

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

229 'rx', 'ry', 'rz', _s_, 's1',

230 _sx_, _sy_, _sz_)

231

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

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

234

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

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

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

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

239

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

241 '''

242 xyz1 = x, y, z, _1_0

243 s1 = self.s1

244 if inverse:

245 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

253 name=self.name)

254

255

256 class Transforms(_NamedEnum):

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

258 to accommodate the L{_LazyNamedEnumItem} properties.

259 '''

260 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

262 '''

263 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

264

265 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

268 Transforms._assert(

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

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

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

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

273 s=_F( 1.2727)),

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

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

276 s=_F( -8.3)),

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

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

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

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

281 DHDNE = _lazy(_DHDNE_, tx=_F(-612.4), ty=_F(-77.0), tz=_F(-440.2),

282 # <https://EPSG.org/transformation_15869/DHDN-to-WGS-84-3.html>

283 sx=_F( 0.054), sy=_F( -0.057), sz=_F( 2.797),

284 s=_F( -2.55)), # East Germany

285 DHDNW = _lazy(_DHDNW_, tx=_F(-598.1), ty=_F(-73.7), tz=_F(-418.2),

286 # <https://EPSG.org/transformation_1777/DHDN-to-WGS-84-2.html>

287 sx=_F( -0.202), sy=_F( -0.045), sz=_F( 2.455),

288 s=_F( -6.7)), # West Germany

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

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

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

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

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

294 Identity = _lazy(_Identity_),

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

296 # <https://EPSG.org/transformation_1641/TM65-to-WGS-84-2.html>

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

298 s=_F( -8.15)),

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

300 # <https://EPSG.org/transformation_1954/TM75-to-WGS-84-2.html>

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

302 s=_F( -8.15)),

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

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

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

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

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

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

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

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

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

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

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

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

315 s=_F(-0.0015)),

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

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

318 # <https://EPSG.org/transformation_1314/OSGB36-to-WGS-84-6.html>

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

320 s=_F( 20.4894)),

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

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

323 WGS84 = _lazy(_WGS84_), # unity

324)

325

326

327 class Datum(_NamedEnumItem):

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

329 '''

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

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

332

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

334 '''New L{Datum}.

335

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

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

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

339

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

341

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

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

344 not a L{Transform}.

345 '''

346 self._ellipsoid = ellipsoid or Datum._ellipsoid

347 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

348

349 self._transform = transform or Datum._transform

350 _xinstanceof(Transform, transform=self.transform)

351

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

353

354 def __eq__(self, other):

355 '''Compare this and an other datum.

356

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

358

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

360 '''

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

362 self.ellipsoid == other.ellipsoid and

363 self.transform == other.transform)

364

365 def __hash__(self):

366 return self._hash # memoized

367

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

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

370

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

372 '''

373 try: # only CartesianBase and LatLonEllipsoidalBase

374 return other.toDatum(self)

375 except AttributeError:

376 pass

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

378

379 def ecef(self, Ecef=None):

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

381

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

383

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

385

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

387

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

389 '''

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

391

392 @Property_RO

393 def ellipsoid(self):

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

395 '''

396 return self._ellipsoid

397

398 @Property_RO

399 def exactTM(self):

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

401 '''

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

403

404 @Property_RO

405 def _hash(self):

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

407

408 @property_RO

409 def isEllipsoidal(self):

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

411 '''

412 return self.ellipsoid.isEllipsoidal

413

414 @property_RO

415 def isOblate(self):

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

417 '''

418 return self.ellipsoid.isOblate

419

420 @property_RO

421 def isProlate(self):

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

423 '''

424 return self.ellipsoid.isProlate

425

426 @property_RO

427 def isSpherical(self):

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

429 '''

430 return self.ellipsoid.isSpherical

431

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

433 '''Return this datum as a string.

434

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

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

437 this datum's name.

438

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

440 '''

441 t = [] if name is None else \

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

443 for a in (_ellipsoid_, _transform_):

444 v = getattr(self, a)

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

446 return sep.join(t)

447

448 @Property_RO

449 def transform(self):

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

451 '''

452 return self._transform

453

454

455 def _earth_datum(inst, a_earth, f=None, name=NN, raiser=_a_ellipsoid_): # in .karney, .trf, ...

456 '''(INTERNAL) Set C{inst._datum} from C{(B{a_..}, B{f})} or C{B{.._ellipsoid}}

457 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, C{a_f2Tuple} or C{scalar} earth radius).

458

459 @note: C{B{raiser}='a_ellipsoid'} for backward naming compatibility.

460 '''

461 if f is not None:

462 E, n, D = _EnD3((a_earth, f), name)

463 if raiser and not E:

464 raise _TypeError(f=f, **{raiser: a_earth})

465 elif a_earth is _EWGS84 or a_earth in (_EWGS84, _WGS84, None):

466 return

467 elif isinstance(a_earth, Datum):

468 E, n, D = None, NN, a_earth

469 else:

470 E, n, D = _EnD3(a_earth, name)

471 if raiser and not E:

472 _xinstanceof(Ellipsoid, Ellipsoid2, a_f2Tuple, Datum, **{raiser: a_earth})

473 if D is None:

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

475 inst._datum = D

476

477

478 def _earth_ellipsoid(earth, *name_raiser):

479 '''(INTERAL) Return the ellipsoid for the given C{earth} model.

480 '''

481 return Ellipsoids.Sphere if earth is R_M else (

482 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

483 _spherical_datum(earth, *name_raiser).ellipsoid)

484

485

486 def _ED2(radius, name):

487 '''(INTERNAL) Helper for C{_EnD3} and C{_spherical_datum}.

488 '''

489 D = Datums.Sphere

490 E = D.ellipsoid

491 if name or radius != E.a: # != E.b

492 n = _under(name)

493 E = Ellipsoid(radius, radius, name=n)

494 D = Datum(E, transform=Transforms.Identity, name=n)

495 return E, D

496

497

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

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

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

501

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

503 '''

504 if isinstance(earth, Datum):

505 D = earth

506 else:

507 E, n, D = _EnD3(earth, name)

508 if not E:

509 n = raiser or _earth_

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

511 if D is None:

512 D = Datum(E, transform=Transforms.Identity, name=_under(n))

513 if raiser and not D.isEllipsoidal:

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

515 return D

516

517

518 def _EnD3(earth, name):

519 '''(INTERNAL) Helper for C{_earth_datum} and C{_ellipsoidal_datum}.

520 '''

521 D = None

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

523 E = earth

524 n = _under(name or E.name)

525 elif isinstance(earth, Datum):

526 E = earth.ellipsoid

527 n = _under(name or earth.name)

528 D = earth

529 elif isscalar(earth):

530 E, D = _ED2(Radius_(earth), name)

531 n = E.name

532 elif isinstance(earth, a_f2Tuple):

533 n = _under(name or earth.name)

534 E = earth.ellipsoid(name=n)

535 elif islistuple(earth, minum=2):

536 E = Ellipsoids.Sphere

537 a, f = earth[:2]

538 if f or a != E.a: # != E.b

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

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

541 else:

542 n = E.name

543 D = Datums.Sphere

544 else:

545 E, n = None, NN

546 return E, n, D

547

548

549 def _mean_radius(radius, *lats):

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

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

552 '''

553 if radius is R_M:

554 r = radius

555 elif isscalar(radius):

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

557 else:

558 E = _ellipsoidal_datum(radius).ellipsoid

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

560 return r

561

562

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

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

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

566

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

568 '''

569 if isscalar(earth):

570 _, D = _ED2(Radius_(earth, Error=Error), name)

571 else:

572 D = _ellipsoidal_datum(earth, Error=Error, name=name)

573 if raiser and not D.isSpherical:

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

575 return D

576

577

578 class Datums(_NamedEnum):

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

580 to accommodate the L{_LazyNamedEnumItem} properties.

581 '''

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

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

584 '''

585 return Datum(Ellipsoids.get(ellipsoid_name),

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

587

588 Datums = Datums(Datum) # PYCHOK singleton

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

590# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

594 Datums._assert(

595 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

599

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

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

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

603

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

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

606

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

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

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

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

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

612

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

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

615

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

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

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

619

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

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

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

623

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

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

626

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

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

629

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

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

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

633

634 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

636

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

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

639

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

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

642

643 # XXX psuedo-ellipsoids for spherical LatLon

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

645

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

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

648

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

650 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

651

652 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

653)

654

655_WGS84 = Datums.WGS84

656 assert _WGS84.ellipsoid is _EWGS84

657

658 if __name__ == '__main__':

659

660 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

661 from pygeodesy.lazily import printf

662

663 # __doc__ of this file, force all into registery

664 for r in (Datums, Transforms):

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

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

667

668# **) MIT License

669#

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

671#

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

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

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

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

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

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

678#

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

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

681#

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

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

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

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

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

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

688# OTHER DEALINGS IN THE SOFTWARE.