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 near-spherical

10 earth models.

11

12@see: 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<https://Web.Archive.org/web/20141009142830/

16 http://earth-info.nga.mil/GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf>}.

17'''

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

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

20

21 from pygeodesy.basics import isscalar, _xinstanceof

22 from pygeodesy.constants import PI_2, R_MA, _2_0

23 from pygeodesy.datums import Datum, _spherical_datum

24 from pygeodesy.dms import clipDegrees, parseDMS2

25 from pygeodesy.errors import _parseX, _ValueError, _xattr, _xkwds

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

27 _northing_, _radius_, _SPACE_, _splituple, _x_, _y_

28# from pygeodesy.lazily import _ALL_LAZY from .named

29 from pygeodesy.named import _NamedBase, _NamedTuple, _ALL_LAZY

30 from pygeodesy.namedTuples import LatLon2Tuple, LatLonDatum3Tuple, PhiLam2Tuple

31 from pygeodesy.props import deprecated_method, Property_RO

32 from pygeodesy.streprs import Fmt, strs, _xzipairs

33 from pygeodesy.units import Easting, Lat, Northing, Radius

34 from pygeodesy.utily import degrees90, degrees180

35

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

37

38__all__ = _ALL_LAZY.webmercator

39__version__ = '23.06.12'

40

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

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

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

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

45

46

47 class EasNorRadius3Tuple(_NamedTuple):

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

49 '''

50 _Names_ = (_easting_, _northing_, _radius_)

51 _Units_ = ( Easting, Northing, Radius)

52

53

54 class WebMercatorError(_ValueError):

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

56 '''

57 pass

58

59

60 class Wm(_NamedBase):

61 '''Web Mercator (WM) coordinate.

62 '''

63 _datum = None # set further below

64 _earths = () # dito

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

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

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

68

69 def __init__(self, x, y, earth=R_MA, name=NN, **radius):

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

71

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

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

74 @kwarg earth: Earth radius (C{meter}), datum or

75 ellipsoid (L{Datum}, L{a_f2Tuple},

76 L{Ellipsoid} or L{Ellipsoid2}).

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

78 @kwarg radius: DEPRECATED, use keyword argument B{C{earth}}.

79

80 @note: WM is strictly defined for spherical and WGS84

81 ellipsoidal earth models only.

82

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

84 '''

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

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

87

88 R = radius.get(_radius_, earth)

89 if R not in Wm._earths:

90 self._datum = _datum(R, _radius_ if radius else _earth_)

91 self._radius = self.datum.ellipsoid.a

92

93 if name:

94 self.name = name

95

96 @Property_RO

97 def datum(self):

98 '''Get the datum (C{Datum}).

99 '''

100 return self._datum

101

102 @Property_RO

103 def ellipsoid(self):

104 '''Get the ellipsoid (C{Ellipsoid}).

105 '''

106 return self.datum.ellipsoid

107

108 @Property_RO

109 def latlon(self):

110 '''Get the lat- and longitude (L{LatLon2Tuple}).

111 '''

112 return self.latlon2()

113

114 def latlon2(self, datum=None):

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

116

117 @kwarg datum: Optional datum (L{Datum}, L{Ellipsoid},

118 L{Ellipsoid2} or L{a_f2Tuple}) or earth

119 radius (C{meter}), overriding this WM's

120 C{radius} and C{datum}.

121

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

123

124 @note: WM is strictly defined for spherical and WGS84

125 ellipsoidal earth models only.

126

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

128

129 @see: Method C{toLatLon} for other return types.

130 '''

131 d = self.datum if datum in (None, self.datum, self.radius) else _datum(datum)

132 E = d.ellipsoid

133 R = self.radius

134 x = self.x / R

135 y = atan(exp(self.y / R)) * _2_0 - PI_2

136 if E.es or E.a != R: # strictly, WGS84 only

137 # <https://Web.Archive.org/web/20141009142830/http://earth-info.nga.mil/

138 # GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf>

139 y = y / R # /= chokes PyChecker

140 y -= E.es_atanh(tanh(y))

141 y *= E.a

142 x *= E.a / R

143

144 return LatLon2Tuple(degrees90(y), degrees180(x), name=self.name)

145

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

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

148

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

150 function L{parseWM}.

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

152 overriding this name.

153

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

155

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

157 '''

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

159 name=name or self.name)

160

161 @deprecated_method

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

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

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

165

166 @Property_RO

167 def philam(self):

168 '''Get the lat- and longitude ((L{PhiLam2Tuple}).

169 '''

170 return PhiLam2Tuple(*map(radians, self.latlon), name=self.name)

171

172 @Property_RO

173 def radius(self):

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

175 '''

176 return self._radius

177

178 @deprecated_method

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

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

181

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

183 '''

184 return self.latlon2(datum=datum)

185

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

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

188

189 @kwarg LatLon: Ellipsoidal or sphperical C{LatLon} class to

190 return the geodetic point (C{LatLon}) or C{None}.

191 @kwarg datum: Optional, datum (C{Datum}) overriding this WM's.

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

193 keyword arguments, ignored if

194 C{B{LatLon} is None}.

195

196 @return: This WM coordinate as B{C{LatLon}} or if

197 C{B{LatLon} is None} a L{LatLonDatum3Tuple}.

198

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

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

201 invalid.

202 '''

203 d = datum or self.datum

204 _xinstanceof(Datum, datum=d)

205 r = self.latlon2(datum=d)

206 r = LatLonDatum3Tuple(r.lat, r.lon, d, name=r.name) if LatLon is None else \

207 LatLon(r.lat, r.lon, **_xkwds(LatLon_kwds, datum=d, name=r.name))

208 return r

209

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

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

212

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

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

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

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

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

218

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

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

221 C{scalar}.

222

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

224 '''

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

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

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

228

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

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

231

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

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

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

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

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

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

238

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

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

241

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

243 '''

244 fs = self.x, self.y

245 if isscalar(radius):

246 fs += (radius,)

247 elif radius: # is True:

248 fs += (self.radius,)

249 elif radius not in (None, False):

250 raise WebMercatorError(radius=radius)

251 t = strs(fs, prec=prec)

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

253

254 @Property_RO

255 def x(self):

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

257 '''

258 return self._x

259

260 @Property_RO

261 def y(self):

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

263 '''

264 return self._y

265

266 Wm._datum = _spherical_datum(Wm._radius, name=Wm.__name__, raiser=_radius_) # PYCHOK defaults

267 Wm._earths = (Wm._radius, Wm._datum, Wm._datum.ellipsoid)

268

269

270 def _datum(earth, name=_datum_):

271 '''(INTERNAL) Make a datum from an C{earth} radius, datum or ellipsoid.

272 '''

273 if earth in Wm._earths:

274 return Wm._datum

275 try:

276 return _spherical_datum(earth, name=name)

277 except Exception as x:

278 raise WebMercatorError(name, earth, cause=x)

279

280

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

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

283 consisting of easting, northing and an optional radius.

284

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

286 @kwarg radius: Optional earth radius (C{meter}), needed in

287 case B{C{strWM}} doesn't include C{r}.

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

289 or C{None}.

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

291

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

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

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

295

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

297 '''

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

299 w = _splituple(strWM)

300

301 if len(w) == 2:

302 w += (radius,)

303 elif len(w) != 3:

304 raise ValueError

305 x, y, R = map(float, w)

306

307 return EasNorRadius3Tuple(x, y, R, name=name) if Wm is None else \

308 Wm(x, y, earth=R, name=name)

309

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

311 strWM=strWM, Error=WebMercatorError)

312

313

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

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

316

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

318 spherical) geodetic C{LatLon} point.

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

320 @kwarg earth: Earth radius (C{meter}), datum or ellipsoid

321 (L{Datum}, L{a_f2Tuple}, L{Ellipsoid} or

322 L{Ellipsoid2}), overridden by B{C{latlon}}'s

323 datum if present.

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

325 or C{None}.

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

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

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

329

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

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

332

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

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

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

336 B{C{earth}} is invalid.

337 '''

338 if _radius_ in Wm_kwds: # remove DEPRECATED, radius

339 d = _datum(Wm_kwds.pop(_radius_), _radius_)

340 else:

341 d = _datum(earth, _earth_)

342 try:

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

344 y = clipDegrees(y, _LatLimit)

345 d = _xattr(latlon, datum=d)

346 n = name or _xattr(latlon, name=NN)

347 except AttributeError:

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

349 n = name

350

351 E = d.ellipsoid

352 R = E.a

353 s = sin(radians(y))

354 y = atanh(s) # == log(tand((90 + lat) / 2)) == log(tanPI_2_2(radians(lat)))

355 if E.es:

356 y -= E.es_atanh(s) # strictly, WGS84 only

357 y *= R

358 x = R * radians(x)

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

360 Wm(x, y, **_xkwds(Wm_kwds, earth=d, name=n))

361 return r

362

363# **) MIT License

364#

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

366#

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

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

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

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

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

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

373#

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

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

376#

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

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

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

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

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

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

383# OTHER DEALINGS IN THE SOFTWARE.