Coverage for pygeodesy/webmercator.py: 98%
118 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
2# -*- coding: utf-8 -*-
4u'''Web Mercator (WM) projection.
6Classes L{Wm} and L{WebMercatorError} and functions L{parseWM} and L{toWm}.
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.
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'''
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
39from math import atan, atanh, exp, radians, sin, tanh
41__all__ = _ALL_LAZY.webmercator
42__version__ = '22.09.12'
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})
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)
57class WebMercatorError(_ValueError):
58 '''Web Mercator (WM) parser or L{Wm} issue.
59 '''
60 pass
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})
70 def __init__(self, x, y, radius=R_MA, name=NN):
71 '''New L{Wm} Web Mercator (WM) coordinate.
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}).
78 @raise WebMercatorError: Invalid B{C{x}}, B{C{y}} or B{C{radius}}.
80 @example:
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)
90 if name:
91 self.name = name
93 @Property_RO
94 def latlon(self):
95 '''Get the lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).
96 '''
97 return self.latlon2()
99 def latlon2(self, datum=None):
100 '''Convert this WM coordinate to a lat- and longitude.
102 @kwarg datum: Optional, ellipsoidal datum (L{Datum},
103 L{Ellipsoid}, L{Ellipsoid2} or
104 L{a_f2Tuple}) or C{None}.
106 @return: A L{LatLon2Tuple}C{(lat, lon)}.
108 @raise TypeError: Invalid or non-ellipsoidal B{C{datum}}.
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
127 return LatLon2Tuple(Lat(degrees90( y)),
128 Lon(degrees180(x)), name=self.name)
130 def parse(self, strWM, name=NN):
131 '''Parse a string to a similar L{Wm} instance.
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.
138 @return: The similar instance (L{Wm}).
140 @raise WebMercatorError: Invalid B{C{strWM}}.
141 '''
142 return parseWM(strWM, radius=self.radius, Wm=self.classof,
143 name=name or self.name)
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)
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)
157 @Property_RO
158 def radius(self):
159 '''Get the earth radius (C{meter}).
160 '''
161 return self._radius
163 @deprecated_method
164 def to2ll(self, datum=None): # PYCHOK no cover
165 '''DEPRECATED, use method C{latlon2}.
167 @return: A L{LatLon2Tuple}C{(lat, lon)}.
168 '''
169 return self.latlon2(datum=datum)
171 def toLatLon(self, LatLon, datum=None, **LatLon_kwds):
172 '''Convert this WM coordinate to a geodetic point.
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.
181 @return: Point of this WM coordinate (B{C{LatLon}}).
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.
187 @example:
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)
202 r = self.latlon2(datum=datum)
203 return LatLon(r.lat, r.lon, **_xkwds(kwds, name=r.name))
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.
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.
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}.
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)
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.
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.
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}.
237 @raise WebMercatorError: Invalid B{C{radius}}.
239 @example:
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)
255 @Property_RO
256 def x(self):
257 '''Get the easting (C{meter}).
258 '''
259 return self._x
261 @Property_RO
262 def y(self):
263 '''Get the northing (C{meter}).
264 '''
265 return self._y
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.
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}).
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}.
282 @raise WebMercatorError: Invalid B{C{strWM}}.
284 @example:
286 >>> u = parseWM('448251 5411932')
287 >>> u.toRepr() # [E:448251, N:5411932]
288 '''
289 def _WM(strWM, radius, Wm, name):
290 w = _splituple(strWM)
292 if len(w) == 2:
293 w += (radius,)
294 elif len(w) != 3:
295 raise ValueError
296 x, y, r = map(float, w)
298 return EasNorRadius3Tuple(x, y, r, name=name) if Wm is None else \
299 Wm(x, y, radius=r, name=name)
301 return _parseX(_WM, strWM, radius, Wm, name,
302 strWM=strWM, Error=WebMercatorError)
305def toWm(latlon, lon=None, radius=R_MA, Wm=Wm, name=NN, **Wm_kwds):
306 '''Convert a lat-/longitude point to a WM coordinate.
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}.
319 @return: The WM coordinate (B{C{Wm}}) or if B{C{Wm}} is C{None}
320 an L{EasNorRadius3Tuple}C{(easting, northing, radius)}.
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.
327 @example:
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)
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
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.