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.869, ty=-52.2978, tz=103.724, s1=1.0, rx=-1.63174e-06, ry=-2.21538e-06, rz=-8.93114e-06, s=1.2727, sx=-0.33657, sy=-0.456955, sz=-1.84218)

47@var Transforms.Bessel1841: Transform(name='Bessel1841', tx=-582, ty=-105, tz=-414, s1=0.999992, rx=-5.04206e-06, ry=-1.69685e-06, rz=1.49323e-05, s=-8.3, sx=-1.04, sy=-0.35, sz=3.08)

48@var Transforms.Clarke1866: Transform(name='Clarke1866', tx=8.0, ty=-160, tz=-176, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

49@var Transforms.DHDN: Transform(name='DHDN', tx=-591.28, ty=-81.35, tz=-396.39, s1=0.99999, rx=7.1607e-06, ry=-3.56823e-07, rz=-7.06858e-06, s=-9.82, sx=1.477, sy=-0.0736, sz=-1.458)

50@var Transforms.DHDNE: Transform(name='DHDNE', tx=-612.4, ty=-77, tz=-440.2, s1=0.999997, rx=2.61799e-07, ry=-2.76344e-07, rz=1.35602e-05, s=-2.55, sx=0.054, sy=-0.057, sz=2.797)

51@var Transforms.DHDNW: Transform(name='DHDNW', tx=-598.1, ty=-73.7, tz=-418.2, s1=0.999993, rx=-9.79324e-07, ry=-2.18166e-07, rz=1.19022e-05, s=-6.7, sx=-0.202, sy=-0.045, sz=2.455)

52@var Transforms.ED50: Transform(name='ED50', tx=89.5, ty=93.8, tz=123.1, s1=0.999999, rx=0.0, ry=0.0, rz=7.56309e-07, s=-1.2, sx=0.0, sy=0.0, sz=0.156)

53@var Transforms.Identity: Transform(name='Identity', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

54@var Transforms.Irl1965: Transform(name='Irl1965', tx=-482.53, ty=130.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-06, s=-8.15, sx=1.042, sy=0.214, sz=0.631)

55@var Transforms.Irl1975: Transform(name='Irl1975', tx=-482.53, ty=130.596, tz=-564.557, s1=0.999992, rx=5.05176e-06, ry=1.0375e-06, rz=3.05917e-06, s=-8.15, sx=1.042, sy=0.214, sz=0.631)

56@var Transforms.Krassovski1940: Transform(name='Krassovski1940', tx=-24, ty=123.0, tz=94.0, s1=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)

57@var Transforms.Krassowsky1940: Transform(name='Krassowsky1940', tx=-24, ty=123.0, tz=94.0, s1=0.999998, rx=-9.69627e-08, ry=1.26052e-06, rz=6.30258e-07, s=-2.423, sx=-0.02, sy=0.26, sz=0.13)

58@var Transforms.MGI: Transform(name='MGI', tx=-577.326, ty=-90.129, tz=-463.92, s1=0.999998, rx=2.49049e-05, ry=7.14615e-06, rz=2.56806e-05, s=-2.423, sx=5.137, sy=1.474, sz=5.297)

59@var Transforms.NAD27: Transform(name='NAD27', tx=8.0, ty=-160, tz=-176, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

60@var Transforms.NAD83: Transform(name='NAD83', tx=1.004, ty=-1.91, tz=-0.515, s1=1.0, rx=1.29445e-07, ry=1.64837e-09, rz=5.33295e-08, s=-0.0015, sx=0.0267, sy=0.00034, sz=0.011)

61@var Transforms.NTF: Transform(name='NTF', tx=-168, ty=-60, tz=320.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

62@var Transforms.OSGB36: Transform(name='OSGB36', tx=-446.448, ty=125.157, tz=-542.06, s1=1.00002, rx=-7.2819e-07, ry=-1.19749e-06, rz=-4.08262e-06, s=20.4894, sx=-0.1502, sy=-0.247, sz=-0.8421)

63@var Transforms.TokyoJapan: Transform(name='TokyoJapan', tx=148.0, ty=-507, tz=-685, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

64@var Transforms.WGS72: Transform(name='WGS72', tx=0.0, ty=0.0, tz=-4.5, s1=1.0, rx=0.0, ry=0.0, rz=2.68587e-06, s=-0.22, sx=0.0, sy=0.0, sz=0.554)

65@var Transforms.WGS84: Transform(name='WGS84', tx=0.0, ty=0.0, tz=0.0, s1=1.0, rx=0.0, ry=0.0, rz=0.0, s=0.0, sx=0.0, sy=0.0, sz=0.0)

66'''

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

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

69

70 from pygeodesy.basics import islistuple, map2, neg, _xinstanceof

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

72# from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase as _CEB, \

73# LatLonEllipsoidalBase as _LLEB # MODS

74 from pygeodesy.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, _EWGS84, \

75 Vector3Tuple

76 from pygeodesy.errors import _IsnotError, _TypeError, _xattr

77 from pygeodesy.fmath import fdot, fmean, Fmt

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

79 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

80 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _MINUS_, \

81 _Krassovski1940_, _Krassowsky1940_, _NAD27_, _NAD83_, _s_, \

82 _Sphere_, _spherical_, _sx_, _sy_, _sz_, _transform_, _tx_, \

83 _ty_, _tz_, _UNDER_, _WGS72_, _WGS84_, _under

84 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

85 from pygeodesy.named import _NamedEnum, _NamedEnumItem, _lazyNamedEnumItem as _lazy

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

87 from pygeodesy.props import Property_RO, property_RO

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

89 from pygeodesy.units import _isRadius, Radius_, radians

90

91# from math import radians # from .units

92

93__all__ = _ALL_LAZY.datums

94__version__ = '24.02.04'

95

96_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

97_BD72_ = 'BD72'

98_DHDN_ = 'DHDN'

99_DHDNE_ = 'DHDNE'

100_DHDNW_ = 'DHDNW'

101_ED50_ = 'ED50'

102_GDA2020_ = 'GDA2020'

103_Identity_ = 'Identity'

104_Irl1965_ = 'Irl1965'

105_Irl1975_ = 'Irl1975'

106_MGI_ = 'MGI'

107_NTF_ = 'NTF'

108_OSGB36_ = 'OSGB36'

109_Potsdam_ = 'Potsdam'

110_RPS = radians(_1_0 / _3600_0) # radians per degree-second

111_TokyoJapan_ = 'TokyoJapan'

112

113

114 class Transform(_NamedEnumItem):

115 '''Helmert I{datum} transformation.

116

117 @see: L{TransformXform<trf.TransformXform>}.

118 '''

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

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

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

122

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

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

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

126

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

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

129

130 sx = _0_0 # x rotation (C{degree-seconds})

131 sy = _0_0 # y rotation (C{degree-seconds})

132 sz = _0_0 # z rotation (C{degree-seconds})

133

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

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

136 '''New L{Transform}.

137

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

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

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

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

142 @kwarg s: Scale (C{float}), ppm.

143 @kwarg sx: X rotation (C{degree-seconds}).

144 @kwarg sy: Y rotation (C{degree-seconds}).

145 @kwarg sz: Z rotation (C{degree-seconds}).

146

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

148 '''

149 if tx:

150 self.tx = tx

151 if ty:

152 self.ty = ty

153 if tz:

154 self.tz = tz

155 if sx: # secs to rads

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

157 if sy:

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

159 if sz:

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

161 if s:

162 self.s = s

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

164

165 if name and not name.startswith(_MINUS_):

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

176 s == t for s, t in zip(self, other)))

177

178 def __hash__(self):

179 return self._hash # memoized

180

181 def __iter__(self):

182 '''Yield the attribute values.

183 '''

184 for _, x in self.items():

185 yield x

186

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

188 '''Helmert-transform an I{ellipsoidal} cartesian B{C{point}}.

189

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

191 '''

192 m = _MODS.ellipsoidalBase

193 _xinstanceof(m.CartesianEllipsoidalBase, point=point)

194 return point.toTransform(self)

195

196 def __neg__(self):

197 return self.inverse()

198

199 @Property_RO

200 def _hash(self):

201 return hash(x for x in self)

202

203 def items(self):

204 '''Yield each attribute as 2-tuple C{(name, value)}.

205 '''

206 for n in (_tx_, _ty_, _tz_, _s_, _sx_, _sy_, _sz_):

207 yield n, getattr(self, n)

208

209 def inverse(self, name=NN):

210 '''Return the inverse of this transform.

211

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

213

214 @return: Inverse (L{Transform}).

215

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

217 '''

218 d = dict((n, -v) for n, v in self.items())

219 n = name or _minus(self.name)

220 return type(self)(name=n, **d)

221

222 @Property_RO

223 def isunity(self):

224 '''Is this a C{unity, identidy} transform (C{bool}), like WGS84?

225 '''

226 return not any(s for s in self)

227

228 def _rps2(self, s_):

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

230 '''

231 return (_RPS * s_), s_

232

233 def toStr(self, prec=5, fmt=Fmt.g, name=NN, **unused): # PYCHOK expected

234 '''Return this transform as a string.

235

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

237 @kwarg fmt: Optional C{float} format (C{letter}).

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

239 this transform's name.

240

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

242 '''

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

244 's1', 'rx', 'ry', 'rz',

245 _s_, _sx_, _sy_, _sz_, fmt=fmt)

246

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

248 '''Transform a (geocentric) position, forward or inverse.

249

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

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

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

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

254

255 @return: The transformed position (L{Vector3Tuple}C{(x, y, z)}).

256 '''

257 if self.isunity:

258 return Vector3Tuple(x, y, z, name=self.name)

259

260 xyz1 = x, y, z, _1_0

261 s1 = self.s1

262 if inverse:

263 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

271 name=self.name)

272

273

274 class Transforms(_NamedEnum):

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

276 to accommodate the L{_LazyNamedEnumItem} properties.

277 '''

278 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

280 '''

281 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

282

283 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

286 Transforms._assert(

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

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

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

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

291 s=_F( 1.2727)),

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

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

294 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

312 Identity = _lazy(_Identity_),

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

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

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

316 s=_F( -8.15)),

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

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

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

320 s=_F( -8.15)),

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

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

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

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

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

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

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

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

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

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

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

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

333 s=_F(-0.0015)),

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

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

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

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

338 s=_F( 20.4894)),

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

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

341 WGS84 = _lazy(_WGS84_), # unity

342)

343

344

345 class Datum(_NamedEnumItem):

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

347 '''

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

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

350

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

352 '''New L{Datum}.

353

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

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

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

357

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

359

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

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

362 not a L{Transform}.

363 '''

364 self._ellipsoid = ellipsoid or Datum._ellipsoid

365 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

366

367 self._transform = transform or Datum._transform

368 _xinstanceof(Transform, transform=self.transform)

369

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

371

372 def __eq__(self, other):

373 '''Compare this and an other datum.

374

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

376

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

378 '''

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

380 self.ellipsoid == other.ellipsoid and

381 self.transform == other.transform)

382

383 def __hash__(self):

384 return self._hash # memoized

385

386 def __matmul__(self, point): # PYCHOK Python 3.5+

387 '''Convert an I{ellipsoidal} B{C{point}} to this datum.

388

389 @raise TypeError: Invalid B{C{point}}.

390 '''

391 m = _MODS.ellipsoidalBase

392 _xinstanceof(m.CartesianEllipsoidalBase, m.LatLonEllipsoidalBase, point=point)

393 return point.toDatum(self)

394

395 def ecef(self, Ecef=None):

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

397

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

399

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

401

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

403

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

405 '''

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

407

408 @Property_RO

409 def ellipsoid(self):

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

411 '''

412 return self._ellipsoid

413

414 @Property_RO

415 def exactTM(self):

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

417 '''

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

419

420 @Property_RO

421 def _hash(self):

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

423

424 @property_RO

425 def isEllipsoidal(self):

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

427 '''

428 return self.ellipsoid.isEllipsoidal

429

430 @property_RO

431 def isOblate(self):

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

433 '''

434 return self.ellipsoid.isOblate

435

436 @property_RO

437 def isProlate(self):

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

439 '''

440 return self.ellipsoid.isProlate

441

442 @property_RO

443 def isSpherical(self):

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

445 '''

446 return self.ellipsoid.isSpherical

447

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

449 '''Return this datum as a string.

450

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

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

453 this datum's name.

454

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

456 '''

457 t = [] if name is None else \

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

459 for a in (_ellipsoid_, _transform_):

460 v = getattr(self, a)

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

462 return sep.join(t)

463

464 @Property_RO

465 def transform(self):

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

467 '''

468 return self._transform

469

470

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

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

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

474

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

476 '''

477 if f is not None:

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

479 if raiser and not E:

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

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

482 return

483 elif isinstance(a_earth, Datum):

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

485 else:

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

487 if raiser and not E:

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

489 if D is None:

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

491 inst._datum = D

492

493

494 def _earth_ellipsoid(earth, *name_raiser):

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

496 '''

497 return Ellipsoids.Sphere if earth is R_M else (

498 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

499 _spherical_datum(earth, *name_raiser).ellipsoid)

500

501

502 def _ED2(radius, name):

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

504 '''

505 D = Datums.Sphere

506 E = D.ellipsoid

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

508 n = _under(name)

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

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

511 return E, D

512

513

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

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

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

517

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

519 '''

520 if isinstance(earth, Datum):

521 D = earth

522 else:

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

524 if not E:

525 n = raiser or _earth_

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

527 if D is None:

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

529 if raiser and not D.isEllipsoidal:

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

531 return D

532

533

534 def _EnD3(earth, name):

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

536 '''

537 D = None

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

539 E = earth

540 n = _under(name or E.name)

541 elif isinstance(earth, Datum):

542 E = earth.ellipsoid

543 n = _under(name or earth.name)

544 D = earth

545 elif _isRadius(earth):

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

547 n = E.name

548 elif isinstance(earth, a_f2Tuple):

549 n = _under(name or earth.name)

550 E = earth.ellipsoid(name=n)

551 elif islistuple(earth, minum=2):

552 E = Ellipsoids.Sphere

553 a, f = earth[:2]

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

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

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

557 else:

558 n = E.name

559 D = Datums.Sphere

560 else:

561 E, n = None, NN

562 return E, n, D

563

564

565 def _mean_radius(radius, *lats):

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

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

568 '''

569 if radius is R_M:

570 r = radius

571 elif _isRadius(radius):

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

573 else:

574 E = _ellipsoidal_datum(radius).ellipsoid

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

576 return r

577

578

579 def _minus(name): # in .trf

580 '''(INTERNAL) Name of C{inverse} Xform.

581 '''

582 m = _MINUS_

583 return name[len(m):] if name.startswith(m) else NN(m, name)

584

585

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

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

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

589

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

591 '''

592 if _isRadius(earth):

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

594 else:

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

596 if raiser and not D.isSpherical:

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

598 return D

599

600

601 class Datums(_NamedEnum):

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

603 to accommodate the L{_LazyNamedEnumItem} properties.

604 '''

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

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

607 '''

608 return Datum(Ellipsoids.get(ellipsoid_name),

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

610

611 Datums = Datums(Datum) # PYCHOK singleton

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

613# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

617 Datums._assert(

618 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

622

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

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

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

626

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

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

629

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

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

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

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

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

635

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

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

638

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

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

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

642

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

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

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

646

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

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

649

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

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

652

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

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

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

656

657 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

659

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

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

662

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

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

665

666 # XXX psuedo-ellipsoids for spherical LatLon

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

668

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

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

671

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

673 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

674

675 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

676)

677

678_WGS84 = Datums.WGS84

679 assert _WGS84.ellipsoid is _EWGS84

680# assert _WGS84.transform.isunity

681

682 if __name__ == '__main__':

683

684 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

685 from pygeodesy.lazily import printf

686

687 # __doc__ of this file, force all into registery

688 for r in (Datums, Transforms):

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

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

691

692# **) MIT License

693#

694# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

695#

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

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

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

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

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

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

702#

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

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

705#

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

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

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

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

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

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

712# OTHER DEALINGS IN THE SOFTWARE.