1

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

3

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

5

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

7

8 Pure 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

10 near-spherical earth models.

11

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

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

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

15 and 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

19 from pygeodesy.basics import isscalar, issubclassof

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

21 from pygeodesy.datums import _ellipsoidal_datum, _ALL_LAZY

22 from pygeodesy.dms import clipDegrees, parseDMS2

23 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB

24 from pygeodesy.errors import _parseX, _TypeError, _ValueError, _xkwds

25 from pygeodesy.interns import NN, _COMMASPACE_, _datum_, _easting_, \

26 _northing_, _radius_, _SPACE_, _splituple

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

28# from pygeodesy.lazily import _ALL_LAZY from .datums

29 from pygeodesy.named import _NamedBase, _NamedTuple

30 from pygeodesy.namedTuples import LatLon2Tuple, PhiLam2Tuple

31 from pygeodesy.props import deprecated_method, Property_RO

32 from pygeodesy.streprs import Fmt, strs, _xzipairs

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

34 Radius, Radius_

35 from pygeodesy.utily import degrees90, degrees180

36

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

38

39__all__ = _ALL_LAZY.webmercator

40__version__ = '23.05.26'

41

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

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

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

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

46

47

48 class EasNorRadius3Tuple(_NamedTuple):

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

50 '''

51 _Names_ = (_easting_, _northing_, _radius_)

52 _Units_ = ( Easting, Northing, Radius)

53

54

55 class WebMercatorError(_ValueError):

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

57 '''

58 pass

59

60

61 class Wm(_NamedBase):

62 '''Web Mercator (WM) coordinate.

63 '''

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

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

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

67

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

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

70

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

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

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

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

75

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

77

78 @example:

79

80 >>> import pygeodesy

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

82 '''

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

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

85 if radius != R_MA:

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

87

88 if name:

89 self.name = name

90

91 @Property_RO

92 def latlon(self):

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

94 '''

95 return self.latlon2()

96

97 def latlon2(self, datum=None):

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

99

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

101 L{Ellipsoid}, L{Ellipsoid2} or

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

103

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

105

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

107

108 @see: Method C{toLatLon}.

109 '''

110 r = self.radius

111 x = self.x / r

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

113 if datum is not None:

114 E = _ellipsoidal_datum(datum, name=self.name, raiser=_datum_).ellipsoid

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

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

117 y = y / r

118 if E.e:

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

120 y *= E.a

121 x *= E.a / r

122

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

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

125

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

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

128

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

130 function L{parseWM}.

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

132 overriding this name.

133

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

135

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

137 '''

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

139 name=name or self.name)

140

141 @deprecated_method

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

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

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

145

146 @Property_RO

147 def philam(self):

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

149 '''

150 r = self.latlon

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

152

153 @Property_RO

154 def radius(self):

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

156 '''

157 return self._radius

158

159 @deprecated_method

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

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

162

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

164 '''

165 return self.latlon2(datum=datum)

166

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

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

169

170 @arg LatLon: Ellipsoidal class to return the geodetic

171 point (C{LatLon}).

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

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

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

175 keyword arguments.

176

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

178

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

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

181 invalid or not ellipsoidal.

182

183 @example:

184

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

186 >>> from pygeodesy import sphericalTrigonometry as sT

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

188 '''

189 e = issubclassof(LatLon, _LLEB)

190 if e and datum:

191 kwds = _xkwds(LatLon_kwds, datum=datum)

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

193 kwds = LatLon_kwds

194 datum = None

195 else:

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

197

198 r = self.latlon2(datum=datum)

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

200

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

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

203

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

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

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

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

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

209

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

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

212 C{scalar}.

213

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

215 '''

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

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

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

219

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

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

222

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

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

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

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

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

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

229

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

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

232

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

234

235 @example:

236

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

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

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

240 '''

241 fs = self.x, self.y

242 if isscalar(radius):

243 fs += (radius,)

244 elif radius: # is True:

245 fs += (self.radius,)

246 elif radius not in (None, False):

247 raise WebMercatorError(radius=radius)

248 t = strs(fs, prec=prec)

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

250

251 @Property_RO

252 def x(self):

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

254 '''

255 return self._x

256

257 @Property_RO

258 def y(self):

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

260 '''

261 return self._y

262

263

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

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

266 consisting of easting, northing and an optional radius.

267

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

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

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

271 or C{None}.

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

273

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

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

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

277

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

279

280 @example:

281

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

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

284 '''

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

286 w = _splituple(strWM)

287

288 if len(w) == 2:

289 w += (radius,)

290 elif len(w) != 3:

291 raise ValueError

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

293

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

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

296

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

298 strWM=strWM, Error=WebMercatorError)

299

300

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

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

303

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

305 spherical) geodetic C{LatLon} point.

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

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

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

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

310 or C{None}.

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

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

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

314

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

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

317

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

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

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

321 B{C{radius}} is invalid.

322

323 @example:

324

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

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

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

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

329 '''

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

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

332 r = Radius(radius)

333 try:

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

335 if isinstance(latlon, _LLEB):

336 E = latlon.datum.ellipsoid

337 e, r = E.e, E.a

338 y = clipDegrees(y, _LatLimit)

339 n = name or latlon.name

340 except AttributeError:

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

342

343 s = sin(radians(y))

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

345 if e:

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

347 y *= r

348 x = r * radians(x)

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

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

351 return r

352

353# **) MIT License

354#

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

356#

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

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

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

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

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

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

363#

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

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

366#

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

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

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

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

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

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

373# OTHER DEALINGS IN THE SOFTWARE.