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.DHDNE: Transform(name='DHDNE', tx=-612.4, ty=-77, tz=-440.2, rx=0, ry=-0, rz=0.00001, s=-2.55, s1=1, 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, rx=-0, ry=-0, rz=0.00001, s=-6.7, s1=0.99999, 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, rx=0, ry=0, rz=0, s=-1.2, s1=1, sx=0, sy=0, sz=0.156)

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

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

55@var Transforms.Irl1975: Transform(name='Irl1975', 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)

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

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

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

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

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

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

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

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

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

65@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) # .isunity

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.ellipsoids import a_f2Tuple, Ellipsoid, Ellipsoid2, Ellipsoids, \

73 _EWGS84, Vector3Tuple

74 from pygeodesy.errors import _IsnotError, _TypeError, _xattr

75 from pygeodesy.fmath import fdot, fmean, Fmt

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

77 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, _earth_, \

78 _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, _Krassovski1940_, \

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

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

81 _WGS72_, _WGS84_, _under

82 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

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

84 _lazyNamedEnumItem as _lazy

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

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

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

88 from pygeodesy.units import _isRadius, radians, Radius_

89

90# from math import radians # from .units

91

92__all__ = _ALL_LAZY.datums

93__version__ = '24.01.25'

94

95_a_ellipsoid_ = _UNDER_(_a_, _ellipsoid_)

96_BD72_ = 'BD72'

97_DHDN_ = 'DHDN'

98_DHDNE_ = 'DHDNE'

99_DHDNW_ = 'DHDNW'

100_ED50_ = 'ED50'

101_GDA2020_ = 'GDA2020'

102_Identity_ = 'Identity'

103_Inverse_ = 'Inverse'

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 def _rps2(s_):

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

116 '''

117 return (_RPS * s_), s_

118

119

120 class Transform(_NamedEnumItem):

121 '''Helmert transformation.

122

123 @see: L{Helmert7Tuple}.

124 '''

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

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

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

128

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

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

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

132

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

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

135

136 sx = _0_0 # x rotation (degree seconds)

137 sy = _0_0 # y rotation (degree seconds)

138 sz = _0_0 # z rotation (degree seconds)

139

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

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

142 '''New L{Transform}.

143

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

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

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

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

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

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

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

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

152

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

154 '''

155 if tx:

156 self.tx = tx

157 if ty:

158 self.ty = ty

159 if tz:

160 self.tz = tz

161 if sx: # secs to rads

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

163 if sy:

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

165 if sz:

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

167 if s:

168 self.s = s

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

170

171 self._register(Transforms, name)

172

173 def __eq__(self, other):

174 '''Compare this and an other transform.

175

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

177

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

179 '''

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

181 and self.tx == other.tx

182 and self.ty == other.ty

183 and self.tz == other.tz

184 and self.rx == other.rx

185 and self.ry == other.ry

186 and self.rz == other.rz

187 and self.s == other.s)

188

189 def __hash__(self):

190 return self._hash # memoized

191

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

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

194

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

196 '''

197 try: # only CartesianBase

198 return other.toTransform(self)

199 except AttributeError:

200 pass

201 raise _IsnotError(_cartesian_, other=other)

202

203 @Property_RO

204 def _hash(self):

205 return hash((self.tx, self.ty, self.tz,

206 self.rx, self.ry, self.rz, self.s))

207

208 def inverse(self, name=NN):

209 '''Return the inverse of this transform.

210

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

212

213 @return: Inverse (Transform).

214

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

216 '''

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

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

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

220

221 @Property_RO

222 def isunity(self):

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

224 '''

225 return not any((self.tx, self.ty, self.tz,

226 self.rx, self.ry, self.rz, self.s))

227

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

229 '''Return this transform as a string.

230

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

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

233 this transform's name.

234

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

236 '''

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

238 'rx', 'ry', 'rz', _s_, 's1',

239 _sx_, _sy_, _sz_)

240

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

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

243

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

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

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

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

248

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

250 '''

251 if self.isunity:

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

253

254 xyz1 = x, y, z, _1_0

255 s1 = self.s1

256 if inverse:

257 xyz1 = map2(neg, xyz1)

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

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

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

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

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

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

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

265 name=self.name)

266

267

268 class Transforms(_NamedEnum):

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

270 to accommodate the L{_LazyNamedEnumItem} properties.

271 '''

272 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s):

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

274 '''

275 return Transform(**name_tx_ty_tz_sx_sy_sz_s)

276

277 Transforms = Transforms(Transform) # PYCHOK singleton

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

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

280 Transforms._assert(

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

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

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

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

285 s=_F( 1.2727)),

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

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

288 s=_F( -8.3)),

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

306 Identity = _lazy(_Identity_),

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

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

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

310 s=_F( -8.15)),

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

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

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

314 s=_F( -8.15)),

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

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

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

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

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

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

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

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

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

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

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

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

327 s=_F(-0.0015)),

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

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

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

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

332 s=_F( 20.4894)),

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

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

335 WGS84 = _lazy(_WGS84_), # unity

336)

337

338

339 class Datum(_NamedEnumItem):

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

341 '''

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

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

344

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

346 '''New L{Datum}.

347

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

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

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

351

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

353

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

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

356 not a L{Transform}.

357 '''

358 self._ellipsoid = ellipsoid or Datum._ellipsoid

359 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid)

360

361 self._transform = transform or Datum._transform

362 _xinstanceof(Transform, transform=self.transform)

363

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

365

366 def __eq__(self, other):

367 '''Compare this and an other datum.

368

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

370

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

372 '''

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

374 self.ellipsoid == other.ellipsoid and

375 self.transform == other.transform)

376

377 def __hash__(self):

378 return self._hash # memoized

379

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

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

382

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

384 '''

385 try: # only CartesianBase and LatLonEllipsoidalBase

386 return other.toDatum(self)

387 except AttributeError:

388 pass

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

390

391 def ecef(self, Ecef=None):

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

393

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

395

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

397

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

399

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

401 '''

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

403

404 @Property_RO

405 def ellipsoid(self):

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

407 '''

408 return self._ellipsoid

409

410 @Property_RO

411 def exactTM(self):

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

413 '''

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

415

416 @Property_RO

417 def _hash(self):

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

419

420 @property_RO

421 def isEllipsoidal(self):

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

423 '''

424 return self.ellipsoid.isEllipsoidal

425

426 @property_RO

427 def isOblate(self):

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

429 '''

430 return self.ellipsoid.isOblate

431

432 @property_RO

433 def isProlate(self):

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

435 '''

436 return self.ellipsoid.isProlate

437

438 @property_RO

439 def isSpherical(self):

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

441 '''

442 return self.ellipsoid.isSpherical

443

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

445 '''Return this datum as a string.

446

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

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

449 this datum's name.

450

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

452 '''

453 t = [] if name is None else \

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

455 for a in (_ellipsoid_, _transform_):

456 v = getattr(self, a)

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

458 return sep.join(t)

459

460 @Property_RO

461 def transform(self):

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

463 '''

464 return self._transform

465

466

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

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

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

470

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

472 '''

473 if f is not None:

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

475 if raiser and not E:

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

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

478 return

479 elif isinstance(a_earth, Datum):

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

481 else:

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

483 if raiser and not E:

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

485 if D is None:

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

487 inst._datum = D

488

489

490 def _earth_ellipsoid(earth, *name_raiser):

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

492 '''

493 return Ellipsoids.Sphere if earth is R_M else (

494 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else

495 _spherical_datum(earth, *name_raiser).ellipsoid)

496

497

498 def _ED2(radius, name):

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

500 '''

501 D = Datums.Sphere

502 E = D.ellipsoid

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

504 n = _under(name)

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

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

507 return E, D

508

509

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

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

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

513

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

515 '''

516 if isinstance(earth, Datum):

517 D = earth

518 else:

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

520 if not E:

521 n = raiser or _earth_

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

523 if D is None:

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

525 if raiser and not D.isEllipsoidal:

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

527 return D

528

529

530 def _EnD3(earth, name):

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

532 '''

533 D = None

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

535 E = earth

536 n = _under(name or E.name)

537 elif isinstance(earth, Datum):

538 E = earth.ellipsoid

539 n = _under(name or earth.name)

540 D = earth

541 elif _isRadius(earth):

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

543 n = E.name

544 elif isinstance(earth, a_f2Tuple):

545 n = _under(name or earth.name)

546 E = earth.ellipsoid(name=n)

547 elif islistuple(earth, minum=2):

548 E = Ellipsoids.Sphere

549 a, f = earth[:2]

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

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

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

553 else:

554 n = E.name

555 D = Datums.Sphere

556 else:

557 E, n = None, NN

558 return E, n, D

559

560

561 def _mean_radius(radius, *lats):

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

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

564 '''

565 if radius is R_M:

566 r = radius

567 elif _isRadius(radius):

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

569 else:

570 E = _ellipsoidal_datum(radius).ellipsoid

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

572 return r

573

574

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

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

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

578

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

580 '''

581 if _isRadius(earth):

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

583 else:

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

585 if raiser and not D.isSpherical:

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

587 return D

588

589

590 class Datums(_NamedEnum):

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

592 to accommodate the L{_LazyNamedEnumItem} properties.

593 '''

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

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

596 '''

597 return Datum(Ellipsoids.get(ellipsoid_name),

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

599

600 Datums = Datums(Datum) # PYCHOK singleton

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

602# Datums with associated ellipsoid and Helmert transform parameters

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

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

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

606 Datums._assert(

607 # Belgian Datum 1972, based on Hayford ellipsoid.

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

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

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

611

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

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

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

615

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

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

618

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

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

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

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

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

624

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

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

627

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

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

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

631

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

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

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

635

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

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

638

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

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

641

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

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

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

645

646 # Nouvelle Triangulation Francaise (Paris) XXX verify

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

648

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

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

651

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

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

654

655 # XXX psuedo-ellipsoids for spherical LatLon

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

657

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

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

660

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

662 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_),

663

664 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_),

665)

666

667_WGS84 = Datums.WGS84

668 assert _WGS84.ellipsoid is _EWGS84

669# assert _WGS84.transform.isunity

670

671 if __name__ == '__main__':

672

673 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

674 from pygeodesy.lazily import printf

675

676 # __doc__ of this file, force all into registery

677 for r in (Datums, Transforms):

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

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

680

681# **) MIT License

682#

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

684#

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

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

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

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

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

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

691#

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

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

694#

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

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

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

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

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

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

701# OTHER DEALINGS IN THE SOFTWARE.