Coverage for pygeodesy/webmercator.py: 98%

118 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-21 13:14 -0400

1 

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

3 

4u'''Web Mercator (WM) projection. 

5 

6Classes L{Wm} and L{WebMercatorError} and functions L{parseWM} and L{toWm}. 

7 

8Pure Python implementation of a U{Web Mercator<https://WikiPedia.org/wiki/Web_Mercator>} 

9(aka I{Pseudo-Mercator}) class and conversion functions for spherical and 

10near-spherical earth models. 

11 

12References U{Google Maps / Bing Maps Spherical Mercator Projection 

13<https://AlastairA.WordPress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection>}, 

14U{Geomatics Guidance Note 7, part 2<https://www.IOGP.org/wp-content/uploads/2019/09/373-07-02.pdf>} 

15and U{Implementation Practice Web Mercator Map Projection 

16<https://Earth-Info.NGA.mil/GandG/wgs84/web_mercator/%28U%29%20NGA_SIG_0011_1.0.0_WEBMERC.pdf>}. 

17''' 

18 

19from pygeodesy.basics import isscalar, issubclassof 

20from pygeodesy.constants import PI_2, R_M, R_MA, R_MB 

21from pygeodesy.datums import _ALL_LAZY, _ellipsoidal_datum 

22from pygeodesy.dms import clipDegrees, parseDMS2 

23from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB 

24from pygeodesy.errors import _IsnotError, _parseX, _TypeError, \ 

25 _ValueError, _xkwds 

26from pygeodesy.interns import NN, _COMMASPACE_, _easting_, \ 

27 _ellipsoidal_, _northing_, _radius_, \ 

28 _SPACE_, _splituple 

29from pygeodesy.interns import _x_, _y_ # PYCHOK used! 

30# from pygeodesy.lazily import _ALL_LAZY from .datums 

31from pygeodesy.named import _NamedBase, _NamedTuple 

32from pygeodesy.namedTuples import LatLon2Tuple, PhiLam2Tuple 

33from pygeodesy.props import deprecated_method, Property_RO 

34from pygeodesy.streprs import Fmt, strs, _xzipairs 

35from pygeodesy.units import Easting, Lam_, Lat, Lon, Northing, Phi_, \ 

36 Radius, Radius_ 

37from pygeodesy.utily import degrees90, degrees180 

38 

39from math import atan, atanh, exp, radians, sin, tanh 

40 

41__all__ = _ALL_LAZY.webmercator 

42__version__ = '22.09.12' 

43 

44# _FalseEasting = 0 # false Easting (C{meter}) 

45# _FalseNorthing = 0 # false Northing (C{meter}) 

46_LatLimit = Lat(limit=85.051129) # latitudinal limit (C{degrees}) 

47# _LonOrigin = 0 # longitude of natural origin (C{degrees}) 

48 

49 

50class EasNorRadius3Tuple(_NamedTuple): 

51 '''3-Tuple C{(easting, northing, radius)}, all in C{meter}. 

52 ''' 

53 _Names_ = (_easting_, _northing_, _radius_) 

54 _Units_ = ( Easting, Northing, Radius) 

55 

56 

57class WebMercatorError(_ValueError): 

58 '''Web Mercator (WM) parser or L{Wm} issue. 

59 ''' 

60 pass 

61 

62 

63class Wm(_NamedBase): 

64 '''Web Mercator (WM) coordinate. 

65 ''' 

66 _radius = R_MA # earth radius (C{meter}) 

67 _x = 0 # Easting (C{meter}) 

68 _y = 0 # Northing (C{meter}) 

69 

70 def __init__(self, x, y, radius=R_MA, name=NN): 

71 '''New L{Wm} Web Mercator (WM) coordinate. 

72 

73 @arg x: Easting from central meridian (C{meter}). 

74 @arg y: Northing from equator (C{meter}). 

75 @kwarg radius: Optional earth radius (C{meter}). 

76 @kwarg name: Optional name (C{str}). 

77 

78 @raise WebMercatorError: Invalid B{C{x}}, B{C{y}} or B{C{radius}}. 

79 

80 @example: 

81 

82 >>> import pygeodesy 

83 >>> w = pygeodesy.Wm(448251, 5411932) 

84 ''' 

85 self._x = Easting( x=x, Error=WebMercatorError) 

86 self._y = Northing(y=y, Error=WebMercatorError) 

87 if radius != R_MA: 

88 self._radius = Radius_(radius, Error=WebMercatorError) 

89 

90 if name: 

91 self.name = name 

92 

93 @Property_RO 

94 def latlon(self): 

95 '''Get the lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}). 

96 ''' 

97 return self.latlon2() 

98 

99 def latlon2(self, datum=None): 

100 '''Convert this WM coordinate to a lat- and longitude. 

101 

102 @kwarg datum: Optional, ellipsoidal datum (L{Datum}, 

103 L{Ellipsoid}, L{Ellipsoid2} or 

104 L{a_f2Tuple}) or C{None}. 

105 

106 @return: A L{LatLon2Tuple}C{(lat, lon)}. 

107 

108 @raise TypeError: Invalid or non-ellipsoidal B{C{datum}}. 

109 

110 @see: Method C{toLatLon}. 

111 ''' 

112 r = self.radius 

113 x = self.x / r 

114 y = 2 * atan(exp(self.y / r)) - PI_2 

115 if datum is not None: 

116 E = _ellipsoidal_datum(datum, name=self.name).ellipsoid 

117 if not E.isEllipsoidal: 

118 raise _IsnotError(_ellipsoidal_, datum=datum) 

119 # <https://Earth-Info.NGA.mil/GandG/wgs84/web_mercator/ 

120 # %28U%29%20NGA_SIG_0011_1.0.0_WEBMERC.pdf> 

121 y = y / r 

122 if E.e: 

123 y -= E.e * atanh(E.e * tanh(y)) # == E.es_atanh(tanh(y)) 

124 y *= E.a 

125 x *= E.a / r 

126 

127 return LatLon2Tuple(Lat(degrees90( y)), 

128 Lon(degrees180(x)), name=self.name) 

129 

130 def parse(self, strWM, name=NN): 

131 '''Parse a string to a similar L{Wm} instance. 

132 

133 @arg strWM: The WM coordinate (C{str}), see 

134 function L{parseWM}. 

135 @kwarg name: Optional instance name (C{str}), 

136 overriding this name. 

137 

138 @return: The similar instance (L{Wm}). 

139 

140 @raise WebMercatorError: Invalid B{C{strWM}}. 

141 ''' 

142 return parseWM(strWM, radius=self.radius, Wm=self.classof, 

143 name=name or self.name) 

144 

145 @deprecated_method 

146 def parseWM(self, strWM, name=NN): # PYCHOK no cover 

147 '''DEPRECATED, use method L{Wm.parse}.''' 

148 return self.parse(strWM, name=name) 

149 

150 @Property_RO 

151 def philam(self): 

152 '''Get the lat- and longitude ((L{PhiLam2Tuple}C{(phi, lam)}). 

153 ''' 

154 r = self.latlon 

155 return PhiLam2Tuple(Phi_(r.lat), Lam_(r.lon), name=r.name) 

156 

157 @Property_RO 

158 def radius(self): 

159 '''Get the earth radius (C{meter}). 

160 ''' 

161 return self._radius 

162 

163 @deprecated_method 

164 def to2ll(self, datum=None): # PYCHOK no cover 

165 '''DEPRECATED, use method C{latlon2}. 

166 

167 @return: A L{LatLon2Tuple}C{(lat, lon)}. 

168 ''' 

169 return self.latlon2(datum=datum) 

170 

171 def toLatLon(self, LatLon, datum=None, **LatLon_kwds): 

172 '''Convert this WM coordinate to a geodetic point. 

173 

174 @arg LatLon: Ellipsoidal class to return the geodetic 

175 point (C{LatLon}). 

176 @kwarg datum: Optional datum for ellipsoidal or C{None} 

177 for spherical B{C{LatLon}} (C{Datum}). 

178 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} 

179 keyword arguments. 

180 

181 @return: Point of this WM coordinate (B{C{LatLon}}). 

182 

183 @raise TypeError: If B{C{LatLon}} and B{C{datum}} are 

184 incompatible or if B{C{datum}} is 

185 invalid or not ellipsoidal. 

186 

187 @example: 

188 

189 >>> w = Wm(448251.795, 5411932.678) 

190 >>> from pygeodesy import sphericalTrigonometry as sT 

191 >>> ll = w.toLatLon(sT.LatLon) # 43°39′11.58″N, 004°01′36.17″E 

192 ''' 

193 e = issubclassof(LatLon, _LLEB) 

194 if e and datum: 

195 kwds = _xkwds(LatLon_kwds, datum=datum) 

196 elif not (e or datum): # and LatLon 

197 kwds = LatLon_kwds 

198 datum = None 

199 else: 

200 raise _TypeError(LatLon=LatLon, datum=datum) 

201 

202 r = self.latlon2(datum=datum) 

203 return LatLon(r.lat, r.lon, **_xkwds(kwds, name=r.name)) 

204 

205 def toRepr(self, prec=3, fmt=Fmt.SQUARE, sep=_COMMASPACE_, radius=False, **unused): # PYCHOK expected 

206 '''Return a string representation of this WM coordinate. 

207 

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

209 @kwarg fmt: Enclosing backets format (C{str}). 

210 @kwarg sep: Optional separator between name:value pairs (C{str}). 

211 @kwarg radius: If C{True} include the radius (C{bool}) or 

212 C{scalar} to override this WM's radius. 

213 

214 @return: This WM as "[x:meter, y:meter]" (C{str}) plus 

215 ", radius:meter]" if B{C{radius}} is C{True} or 

216 C{scalar}. 

217 

218 @raise WebMercatorError: Invalid B{C{radius}}. 

219 ''' 

220 t = self.toStr(prec=prec, sep=None, radius=radius) 

221 n = (_x_, _y_, _radius_)[:len(t)] 

222 return _xzipairs(n, t, sep=sep, fmt=fmt) 

223 

224 def toStr(self, prec=3, fmt=Fmt.F, sep=_SPACE_, radius=False, **unused): # PYCHOK expected 

225 '''Return a string representation of this WM coordinate. 

226 

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

228 @kwarg fmt: The C{float} format (C{str}). 

229 @kwarg sep: Optional separator to join (C{str}) or C{None} 

230 to return an unjoined C{tuple} of C{str}s. 

231 @kwarg radius: If C{True} include the radius (C{bool}) or 

232 C{scalar} to override this WM's radius. 

233 

234 @return: This WM as "meter meter" (C{str}) or as "meter meter 

235 radius" if B{C{radius}} is C{True} or C{scalar}. 

236 

237 @raise WebMercatorError: Invalid B{C{radius}}. 

238 

239 @example: 

240 

241 >>> w = Wm(448251, 5411932.0001) 

242 >>> w.toStr(4) # 448251.0 5411932.0001 

243 >>> w.toStr(sep=', ') # 448251, 5411932 

244 ''' 

245 fs = self.x, self.y 

246 if isscalar(radius): 

247 fs += (radius,) 

248 elif radius: # is True: 

249 fs += (self.radius,) 

250 elif radius not in (None, False): 

251 raise WebMercatorError(radius=radius) 

252 t = strs(fs, prec=prec) 

253 return t if sep is None else sep.join(t) 

254 

255 @Property_RO 

256 def x(self): 

257 '''Get the easting (C{meter}). 

258 ''' 

259 return self._x 

260 

261 @Property_RO 

262 def y(self): 

263 '''Get the northing (C{meter}). 

264 ''' 

265 return self._y 

266 

267 

268def parseWM(strWM, radius=R_MA, Wm=Wm, name=NN): 

269 '''Parse a string C{"e n [r]"} representing a WM coordinate, 

270 consisting of easting, northing and an optional radius. 

271 

272 @arg strWM: A WM coordinate (C{str}). 

273 @kwarg radius: Optional earth radius (C{meter}). 

274 @kwarg Wm: Optional class to return the WM coordinate (L{Wm}) 

275 or C{None}. 

276 @kwarg name: Optional name (C{str}). 

277 

278 @return: The WM coordinate (B{C{Wm}}) or an 

279 L{EasNorRadius3Tuple}C{(easting, northing, radius)} 

280 if B{C{Wm}} is C{None}. 

281 

282 @raise WebMercatorError: Invalid B{C{strWM}}. 

283 

284 @example: 

285 

286 >>> u = parseWM('448251 5411932') 

287 >>> u.toRepr() # [E:448251, N:5411932] 

288 ''' 

289 def _WM(strWM, radius, Wm, name): 

290 w = _splituple(strWM) 

291 

292 if len(w) == 2: 

293 w += (radius,) 

294 elif len(w) != 3: 

295 raise ValueError 

296 x, y, r = map(float, w) 

297 

298 return EasNorRadius3Tuple(x, y, r, name=name) if Wm is None else \ 

299 Wm(x, y, radius=r, name=name) 

300 

301 return _parseX(_WM, strWM, radius, Wm, name, 

302 strWM=strWM, Error=WebMercatorError) 

303 

304 

305def toWm(latlon, lon=None, radius=R_MA, Wm=Wm, name=NN, **Wm_kwds): 

306 '''Convert a lat-/longitude point to a WM coordinate. 

307 

308 @arg latlon: Latitude (C{degrees}) or an (ellipsoidal or 

309 spherical) geodetic C{LatLon} point. 

310 @kwarg lon: Optional longitude (C{degrees} or C{None}). 

311 @kwarg radius: Optional earth radius (C{meter}), overridden 

312 by the B{C{latlon}}'s equatorial radius. 

313 @kwarg Wm: Optional class to return the WM coordinate (L{Wm}) 

314 or C{None}. 

315 @kwarg name: Optional name (C{str}). 

316 @kwarg Wm_kwds: Optional, additional B{C{Wm}} keyword arguments, 

317 ignored if C{B{Wm} is None}. 

318 

319 @return: The WM coordinate (B{C{Wm}}) or if B{C{Wm}} is C{None} 

320 an L{EasNorRadius3Tuple}C{(easting, northing, radius)}. 

321 

322 @raise ValueError: If B{C{lon}} value is missing, if B{C{latlon}} is not 

323 scalar, if B{C{latlon}} is beyond the valid WM range 

324 and L{pygeodesy.rangerrors} is set to C{True} or if 

325 B{C{radius}} is invalid. 

326 

327 @example: 

328 

329 >>> p = LatLon(48.8582, 2.2945) # 448251.8 5411932.7 

330 >>> w = toWm(p) # 448252 5411933 

331 >>> p = LatLon(13.4125, 103.8667) # 377302.4 1483034.8 

332 >>> w = toWm(p) # 377302 1483035 

333 ''' 

334 e, r, n = 0, radius, name 

335 if r not in (R_MA, R_MB, R_M): 

336 r = Radius(radius) 

337 try: 

338 y, x = latlon.lat, latlon.lon 

339 if isinstance(latlon, _LLEB): 

340 E = latlon.datum.ellipsoid 

341 e, r = E.e, E.a 

342 y = clipDegrees(y, _LatLimit) 

343 n = name or latlon.name 

344 except AttributeError: 

345 y, x = parseDMS2(latlon, lon, clipLat=_LatLimit) 

346 

347 s = sin(radians(y)) 

348 y = atanh(s) # == log(tan((90 + lat) / 2)) == log(tanPI_2_2(radians(lat))) 

349 if e: 

350 y -= e * atanh(e * s) 

351 y *= r 

352 x = r * radians(x) 

353 r = EasNorRadius3Tuple(x, y, r, name=n) if Wm is None else \ 

354 Wm(x, y, **_xkwds(Wm_kwds, radius=r, name=n)) 

355 return r 

356 

357# **) MIT License 

358# 

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

360# 

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

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

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

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

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

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

367# 

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

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

370# 

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

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

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

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

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

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

377# OTHER DEALINGS IN THE SOFTWARE.