Coverage for pygeodesy/datums.py: 94%

241 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-12 13:15 -0500

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Datums and transformations thereof. 

5 

6Classes L{Datum} and L{Transform} and registries L{Datums} and L{Transforms}, respectively. 

7 

8Pure Python implementation of geodesy tools for ellipsoidal earth models, including datums 

9and ellipsoid parameters for different geographic coordinate systems and methods for 

10converting between them and to cartesian coordinates. Transcoded from JavaScript originals by 

11I{(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 

14Historical geodetic datums: a latitude/longitude point defines a geographic location on, above 

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

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

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

18reference points. 

19 

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

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

22 

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

24datum, 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 

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

69 

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

71from 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 

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

75 Vector3Tuple 

76from pygeodesy.errors import _IsnotError, _TypeError, _xattr 

77from pygeodesy.fmath import fdot, fmean, Fmt 

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

79 _Clarke1866_, _Clarke1880IGN_, _COMMASPACE_, _DOT_, \ 

80 _earth_, _ellipsoid_, _ellipsoidal_, _GRS80_, _Intl1924_, \ 

81 _MINUS_, _Krassovski1940_, _Krassowsky1940_, _NAD27_, \ 

82 _NAD83_, _s_, _Sphere_, _spherical_, _transform_, \ 

83 _UNDER_, _WGS72_, _WGS84_, _under 

84from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

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

87from pygeodesy.props import Property_RO, property_RO 

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

89from pygeodesy.units import _isRadius, Radius_, radians 

90 

91# from math import radians # from .units 

92 

93__all__ = _ALL_LAZY.datums 

94__version__ = '24.02.12' 

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_Names7 = 'tx', 'ty', 'tz', _s_, 'sx', 'sy', 'sz' # in .trf 

108_Names11 = _Names7[:3] + ('s1', 'rx', 'ry', 'rz') + _Names7[3:] 

109_NTF_ = 'NTF' 

110_OSGB36_ = 'OSGB36' 

111_Potsdam_ = 'Potsdam' 

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

113_S1_S = 1.e-6 # in .trf 

114_TokyoJapan_ = 'TokyoJapan' 

115 

116 

117class Transform(_NamedEnumItem): 

118 '''Helmert I{datum} transformation. 

119 

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

121 ''' 

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

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

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

125 

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

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

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

129 

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

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

132 

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

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

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

136 

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

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

139 '''New L{Transform}. 

140 

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

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

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

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

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

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

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

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

149 

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

151 ''' 

152 if tx: 

153 self.tx = tx 

154 if ty: 

155 self.ty = ty 

156 if tz: 

157 self.tz = tz 

158 if sx: # secs to rads 

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

160 if sy: 

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

162 if sz: 

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

164 if s: 

165 self.s = s 

166 self.s1 = _F(s * _S1_S + _1_0) # normalize ppM to (s + 1) 

167 

168 self._register(Transforms, name) 

169 

170 def __eq__(self, other): 

171 '''Compare this and an other transform. 

172 

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

174 

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

176 ''' 

177 return self is other or (isinstance(other, Transform) and all( 

178 a == b for a, b in zip(self, other))) 

179 

180 def __hash__(self): 

181 return self._hash # memoized 

182 

183 def __iter__(self): 

184 '''Yield the initial attribute values. 

185 ''' 

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

187 yield x 

188 

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

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

191 

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

193 ''' 

194 m = _MODS.ellipsoidalBase 

195 _xinstanceof(m.CartesianEllipsoidalBase, point=point) 

196 return point.toTransform(self) 

197 

198 def __neg__(self): 

199 return self.inverse() 

200 

201# def __sub__(self, other): 

202# _xinstanceof(Transform, other=other) 

203# 

204# def _sub(a, b): 

205# for n in _Names11: 

206# yield getattr(a, n) - getattr(b, n) 

207# 

208# return type(self)(_sub(self, other), name=_MINUS_) # .fsums._sub_op 

209 

210 @Property_RO 

211 def _hash(self): 

212 return hash(x for x in self) 

213 

214 def items(self): 

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

216 ''' 

217 for n in _Names7: 

218 yield n, getattr(self, n) 

219 

220 def inverse(self, name=NN): 

221 '''Return the inverse of this transform. 

222 

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

224 

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

226 

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

228 ''' 

229 d = dict((n, neg(v)) for n, v in self.items()) 

230 n = name or _minus(self.name) 

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

232 

233 @Property_RO 

234 def isunity(self): 

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

236 ''' 

237 return not any(s for s in self) 

238 

239 def _rps2(self, s_): 

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

241 ''' 

242 return (_RPS * s_), s_ 

243 

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

245 '''Return this transform as a string. 

246 

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

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

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

250 this transform's name. 

251 

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

253 ''' 

254 return self._instr(name, prec, *_Names11, fmt=fmt) 

255 

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

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

258 

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

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

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

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

263 

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

265 ''' 

266 if self.isunity: 

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

268 

269 xyz1 = x, y, z, _1_0 

270 s1 = self.s1 

271 if inverse: 

272 xyz1 = map2(neg, xyz1) 

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

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

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

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

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

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

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

280 name=self.name) 

281 

282 

283class Transforms(_NamedEnum): 

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

285 to accommodate the L{_LazyNamedEnumItem} properties. 

286 ''' 

287 def _Lazy(self, **name_tx_ty_tz_sx_sy_sz_s): 

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

289 ''' 

290 return Transform(**name_tx_ty_tz_sx_sy_sz_s) 

291 

292Transforms = Transforms(Transform) # PYCHOK singleton 

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

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

295Transforms._assert( 

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

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

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

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

300 s=_F( 1.2727)), 

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

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

303 s=_F( -8.3)), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

321 Identity = _lazy(_Identity_), 

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

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

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

325 s=_F( -8.15)), 

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

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

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

329 s=_F( -8.15)), 

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

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

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

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

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

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

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

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

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

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

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

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

342 s=_F(-0.0015)), 

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

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

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

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

347 s=_F( 20.4894)), 

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

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

350 WGS84 = _lazy(_WGS84_), # unity 

351) 

352 

353 

354class Datum(_NamedEnumItem): 

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

356 ''' 

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

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

359 

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

361 '''New L{Datum}. 

362 

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

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

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

366 

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

368 

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

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

371 not a L{Transform}. 

372 ''' 

373 self._ellipsoid = ellipsoid or Datum._ellipsoid 

374 _xinstanceof(Ellipsoid, ellipsoid=self.ellipsoid) 

375 

376 self._transform = transform or Datum._transform 

377 _xinstanceof(Transform, transform=self.transform) 

378 

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

380 

381 def __eq__(self, other): 

382 '''Compare this and an other datum. 

383 

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

385 

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

387 ''' 

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

389 self.ellipsoid == other.ellipsoid and 

390 self.transform == other.transform) 

391 

392 def __hash__(self): 

393 return self._hash # memoized 

394 

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

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

397 

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

399 ''' 

400 m = _MODS.ellipsoidalBase 

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

402 return point.toDatum(self) 

403 

404 def ecef(self, Ecef=None): 

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

406 

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

408 

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

410 

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

412 

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

414 ''' 

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

416 

417 @Property_RO 

418 def ellipsoid(self): 

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

420 ''' 

421 return self._ellipsoid 

422 

423 @Property_RO 

424 def exactTM(self): 

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

426 ''' 

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

428 

429 @Property_RO 

430 def _hash(self): 

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

432 

433 @property_RO 

434 def isEllipsoidal(self): 

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

436 ''' 

437 return self.ellipsoid.isEllipsoidal 

438 

439 @property_RO 

440 def isOblate(self): 

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

442 ''' 

443 return self.ellipsoid.isOblate 

444 

445 @property_RO 

446 def isProlate(self): 

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

448 ''' 

449 return self.ellipsoid.isProlate 

450 

451 @property_RO 

452 def isSpherical(self): 

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

454 ''' 

455 return self.ellipsoid.isSpherical 

456 

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

458 '''Return this datum as a string. 

459 

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

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

462 this datum's name. 

463 

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

465 ''' 

466 t = [] if name is None else \ 

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

468 for a in (_ellipsoid_, _transform_): 

469 v = getattr(self, a) 

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

471 return sep.join(t) 

472 

473 @Property_RO 

474 def transform(self): 

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

476 ''' 

477 return self._transform 

478 

479 

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

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

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

483 

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

485 ''' 

486 if f is not None: 

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

488 if raiser and not E: 

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

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

491 return 

492 elif isinstance(a_earth, Datum): 

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

494 else: 

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

496 if raiser and not E: 

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

498 if D is None: 

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

500 inst._datum = D 

501 

502 

503def _earth_ellipsoid(earth, *name_raiser): 

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

505 ''' 

506 return Ellipsoids.Sphere if earth is R_M else ( 

507 _EWGS84 if earth is _EWGS84 or earth is _WGS84 else 

508 _spherical_datum(earth, *name_raiser).ellipsoid) 

509 

510 

511def _ED2(radius, name): 

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

513 ''' 

514 D = Datums.Sphere 

515 E = D.ellipsoid 

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

517 n = _under(name) 

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

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

520 return E, D 

521 

522 

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

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

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

526 

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

528 ''' 

529 if isinstance(earth, Datum): 

530 D = earth 

531 else: 

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

533 if not E: 

534 n = raiser or _earth_ 

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

536 if D is None: 

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

538 if raiser and not D.isEllipsoidal: 

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

540 return D 

541 

542 

543def _EnD3(earth, name): 

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

545 ''' 

546 D = None 

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

548 E = earth 

549 n = _under(name or E.name) 

550 elif isinstance(earth, Datum): 

551 E = earth.ellipsoid 

552 n = _under(name or earth.name) 

553 D = earth 

554 elif _isRadius(earth): 

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

556 n = E.name 

557 elif isinstance(earth, a_f2Tuple): 

558 n = _under(name or earth.name) 

559 E = earth.ellipsoid(name=n) 

560 elif islistuple(earth, minum=2): 

561 E = Ellipsoids.Sphere 

562 a, f = earth[:2] 

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

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

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

566 else: 

567 n = E.name 

568 D = Datums.Sphere 

569 else: 

570 E, n = None, NN 

571 return E, n, D 

572 

573 

574def _mean_radius(radius, *lats): 

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

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

577 ''' 

578 if radius is R_M: 

579 r = radius 

580 elif _isRadius(radius): 

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

582 else: 

583 E = _ellipsoidal_datum(radius).ellipsoid 

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

585 return r 

586 

587 

588def _minus(name): # in .trf 

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

590 ''' 

591 m = _MINUS_ 

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

593 

594 

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

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

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

598 

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

600 ''' 

601 if _isRadius(earth): 

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

603 else: 

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

605 if raiser and not D.isSpherical: 

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

607 return D 

608 

609 

610class Datums(_NamedEnum): 

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

612 to accommodate the L{_LazyNamedEnumItem} properties. 

613 ''' 

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

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

616 ''' 

617 return Datum(Ellipsoids.get(ellipsoid_name), 

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

619 

620Datums = Datums(Datum) # PYCHOK singleton 

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

622# Datums with associated ellipsoid and Helmert transform parameters 

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

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

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

626Datums._assert( 

627 # Belgian Datum 1972, based on Hayford ellipsoid. 

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

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

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

631 

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

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

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

635 

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

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

638 

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

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

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

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

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

644 

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

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

647 

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

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

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

651 

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

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

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

655 

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

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

658 

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

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

661 

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

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

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

665 

666 # Nouvelle Triangulation Francaise (Paris) XXX verify 

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

668 

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

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

671 

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

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

674 

675 # XXX psuedo-ellipsoids for spherical LatLon 

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

677 

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

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

680 

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

682 WGS72 = _lazy(_WGS72_, _WGS72_, _WGS72_), 

683 

684 WGS84 = _lazy(_WGS84_, _WGS84_, _WGS84_), 

685) 

686 

687_WGS84 = Datums.WGS84 

688assert _WGS84.ellipsoid is _EWGS84 

689# assert _WGS84.transform.isunity 

690 

691if __name__ == '__main__': 

692 

693 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

694 from pygeodesy.lazily import printf 

695 

696 # __doc__ of this file, force all into registery 

697 for r in (Datums, Transforms): 

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

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

700 

701# **) MIT License 

702# 

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

704# 

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

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

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

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

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

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

711# 

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

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

714# 

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

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

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

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

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

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

721# OTHER DEALINGS IN THE SOFTWARE.