Coverage for pygeodesy/azimuthal.py: 99%
320 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -0400
2# -*- coding: utf-8 -*-
4u'''Equidistant, Equal-Area, and other Azimuthal projections.
6Classes L{Equidistant}, L{EquidistantExact}, L{EquidistantGeodSolve},
7L{EquidistantKarney}, L{Gnomonic}, L{GnomonicExact}, L{GnomonicKarney},
8L{LambertEqualArea}, L{Orthographic} and L{Stereographic}, classes
9L{AzimuthalError}, L{Azimuthal7Tuple} and functions L{equidistant}
10and L{gnomonic}.
12L{EquidistantExact} and L{GnomonicExact} are based on exact geodesic classes
13L{GeodesicExact} and L{GeodesicLineExact}, Python versions of I{Charles Karney}'s
14C++ original U{GeodesicExact<https://GeographicLib.SourceForge.io/C++/doc/
15classGeographicLib_1_1GeodesicExact.html>}, respectively U{GeodesicLineExact
16<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>}.
18Using L{EquidistantGeodSolve} requires I{Karney}'s utility U{GeodSolve
19<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} to be
20executable and set in env variable C{PYGEODESY_GEODSOLVE}, see module
21L{geodsolve} for more details.
23L{EquidistantKarney} and L{GnomonicKarney} require I{Karney}'s Python package
24U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
26Other azimuthal classes implement only (***) U{Snyder's FORMULAS FOR THE SPHERE
27<https://Pubs.USGS.gov/pp/1395/report.pdf>} and use those for any datum,
28spherical and ellipsoidal. The radius used for the latter is the ellipsoid's
29I{mean radius of curvature} at the latitude of the projection center point. For
30further justification, see the first paragraph under U{Snyder's FORMULAS FOR THE
31ELLIPSOID, page 197<https://Pubs.USGS.gov/pp/1395/report.pdf>}.
33Page numbers in C{Snyder} references apply to U{John P. Snyder, "Map Projections
34-- A Working Manual", 1987<https://Pubs.USGS.gov/pp/1395/report.pdf>}.
36See also U{here<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection>},
37especially the U{Comparison of the Azimuthal equidistant projection and some
38azimuthal projections centred on 90° N at the same scale, ordered by projection
39altitude in Earth radii<https://WikiPedia.org/wiki/Azimuthal_equidistant_projection
40#/media/File:Comparison_azimuthal_projections.svg>}.
41'''
42# make sure int/int division yields float quotient, see .basics
43from __future__ import division as _; del _ # PYCHOK semicolon
45# from pygeodesy.basics import _xinstanceof # from .ellipsoidalBase
46from pygeodesy.constants import EPS, EPS0, EPS1, NAN, isnon0, _umod_360, \
47 _EPStol, _0_0, _0_1, _0_5, _1_0, _N_1_0, _2_0
48from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \
49 _xinstanceof
50from pygeodesy.datums import _spherical_datum, _WGS84
51from pygeodesy.errors import _ValueError, _xdatum, _xkwds
52from pygeodesy.fmath import euclid, hypot as _hypot, Fsum
53# from pygeodesy.fsums import Fsum # from .fmath
54# from pygeodesy.formy import antipode # _MODS
55from pygeodesy.interns import _azimuth_, _datum_, _lat_, _lon_, _scale_, \
56 _SPACE_, _x_, _y_
57from pygeodesy.karney import _norm180
58from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB
59from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS # ALL_MODS
60from pygeodesy.named import _name__, _name2__, _NamedBase, _NamedTuple, _Pass
61from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple
62from pygeodesy.props import deprecated_Property_RO, Property_RO, \
63 property_doc_, _update_all
64from pygeodesy.streprs import Fmt, _fstrLL0, unstr
65from pygeodesy.units import Bearing, Easting, Lat_, Lon_, Northing, \
66 Scalar, Scalar_
67from pygeodesy.utily import asin1, atan1, atan2b, atan2d, sincos2, \
68 sincos2d, sincos2d_
70from math import acos, atan2, degrees, fabs, sin, sqrt
72__all__ = _ALL_LAZY.azimuthal
73__version__ = '24.05.24'
75_EPS_K = _EPStol * _0_1 # Karney's eps_ or _EPSmin * _0_1?
76_over_horizon_ = 'over horizon'
77_TRIPS = 21 # numit, 4 sufficient
80def _enzh4(x, y, *h):
81 '''(INTERNAL) Return 4-tuple (easting, northing, azimuth, hypot).
82 '''
83 e = Easting( x=x)
84 n = Northing(y=y)
85 z = atan2b(e, n) # (x, y) for azimuth from true North
86 return e, n, z, (h[0] if h else _hypot(e, n))
89class _AzimuthalBase(_NamedBase):
90 '''(INTERNAL) Base class for azimuthal projections.
92 @see: I{Karney}'s C++ class U{AzimuthalEquidistant<https://GeographicLib.SourceForge.io/
93 C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>} and U{Gnomonic
94 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>} or
95 the C{PyGeodesy} versions thereof L{EquidistantKarney} respectively L{GnomonicKarney}.
96 '''
97 _datum = _WGS84 # L{Datum}
98 _latlon0 = LatLon2Tuple(_0_0, _0_0) # lat0, lon0 (L{LatLon2Tuple})
99 _sc0 = _0_0, _1_0 # 2-Tuple C{sincos2d(lat0)}
101 def __init__(self, lat0, lon0, datum=None, **name):
102 '''New azimuthal projection.
104 @arg lat0: Latitude of the center point (C{degrees90}).
105 @arg lon0: Longitude of the center point (C{degrees180}).
106 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
107 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
108 radius (C{meter}).
109 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
111 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
113 @raise TypeError: Invalid B{C{datum}}.
114 '''
115 if datum not in (None, self._datum):
116 self._datum = _spherical_datum(datum, **name)
117 if name:
118 self.name = name
120 if lat0 or lon0: # often both 0
121 self._reset(lat0, lon0)
123 @Property_RO
124 def datum(self):
125 '''Get the datum (L{Datum}).
126 '''
127 return self._datum
129 @Property_RO
130 def equatoradius(self):
131 '''Get the geodesic's equatorial radius, semi-axis (C{meter}).
132 '''
133 return self.datum.ellipsoid.a
135 a = equatoradius
137 @Property_RO
138 def flattening(self):
139 '''Get the geodesic's flattening (C{scalar}).
140 '''
141 return self.datum.ellipsoid.f
143 f = flattening
145 def forward(self, lat, lon, **name): # PYCHOK no cover
146 '''I{Must be overloaded}.'''
147 self._notOverloaded(lat, lon, **name)
149 def _forward(self, lat, lon, name, _k_t_2):
150 '''(INTERNAL) Azimuthal (spherical) forward C{lat, lon} to C{x, y}.
151 '''
152 lat, lon = Lat_(lat), Lon_(lon)
153 sa, ca, sb, cb = sincos2d_(lat, lon - self.lon0)
154 s0, c0 = self._sc0
156 cb *= ca
157 k, t = _k_t_2(s0 * sa + c0 * cb)
158 if t:
159 r = k * self.radius
160 y = r * (c0 * sa - s0 * cb)
161 e, n, z, _ = _enzh4(r * sb * ca, y, None)
162 else: # 0 or 180
163 e = n = z = _0_0
165 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
166 name=self._name__(name))
167 return t
169 def _forwards(self, *lls):
170 '''(INTERNAL) One or more C{.forward} calls, see .ellipsoidalBaseDI.
171 '''
172 _fwd = self.forward
173 for ll in lls:
174 yield _fwd(ll.lat, ll.lon)
176 @Property_RO
177 def lat0(self):
178 '''Get the center latitude (C{degrees90}).
179 '''
180 return self._latlon0.lat
182 @property
183 def latlon0(self):
184 '''Get the center lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}) in (C{degrees90}, C{degrees180}).
185 '''
186 return self._latlon0
188 @latlon0.setter # PYCHOK setter!
189 def latlon0(self, latlon0):
190 '''Set the center lat- and longitude (C{LatLon}, L{LatLon2Tuple} or L{LatLon4Tuple}).
192 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or ellipsoidal mismatch
193 of B{C{latlon0}} and this projection.
194 '''
195 B = _LLEB if self.datum.isEllipsoidal else _LLB
196 _xinstanceof(B, LatLon2Tuple, LatLon4Tuple, latlon0=latlon0)
197 if hasattr(latlon0, _datum_):
198 _xdatum(self.datum, latlon0.datum, Error=AzimuthalError)
199 self.reset(latlon0.lat, latlon0.lon)
201 @Property_RO
202 def lon0(self):
203 '''Get the center longitude (C{degrees180}).
204 '''
205 return self._latlon0.lon
207 @deprecated_Property_RO
208 def majoradius(self): # PYCHOK no cover
209 '''DEPRECATED, use property C{equatoradius}.'''
210 return self.equatoradius
212 @Property_RO
213 def radius(self):
214 '''Get this projection's mean radius of curvature (C{meter}).
215 '''
216 return self.datum.ellipsoid.rocMean(self.lat0)
218 def reset(self, lat0, lon0):
219 '''Set or reset the center point of this azimuthal projection.
221 @arg lat0: Center point latitude (C{degrees90}).
222 @arg lon0: Center point longitude (C{degrees180}).
224 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
225 '''
226 _update_all(self) # zap caches
227 self._reset(lat0, lon0)
229 def _reset(self, lat0, lon0):
230 '''(INTERNAL) Update the center point.
231 '''
232 self._latlon0 = LatLon2Tuple(Lat_(lat0=lat0, Error=AzimuthalError),
233 Lon_(lon0=lon0, Error=AzimuthalError))
234 self._sc0 = sincos2d(self.lat0)
236 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK no cover
237 '''I{Must be overloaded}.'''
238 self._notOverloaded(x, y, LatLon=LatLon, **name_LatLon_kwds)
240 def _reverse(self, x, y, _c, lea, LatLon, **name_LatLon_kwds):
241 '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
242 '''
243 e, n, z, r = _enzh4(x, y)
245 c = _c(r / self.radius)
246 if c is None:
247 lat, lon = self.latlon0
248 k, z = _1_0, _0_0
249 else:
250 s0, c0 = self._sc0
251 sc, cc = sincos2(c)
252 k = c / sc
253 s = s0 * cc
254 if r > EPS0:
255 s += c0 * sc * (n / r)
256 lat = degrees(asin1(s))
257 if lea or fabs(c0) > EPS:
258 d = atan2d(e * sc, r * c0 * cc - n * s0 * sc)
259 else:
260 d = atan2d(e, (n if s0 < 0 else -n))
261 lon = _norm180(self.lon0 + d)
263 if LatLon is None:
264 t, _ = _name2__(name_LatLon_kwds, _or_nameof=self)
265 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum, name=t)
266 else:
267 t = self._toLatLon(lat, lon, LatLon, name_LatLon_kwds)
268 return t
270 def _reverse2(self, x_t, *y):
271 '''(INTERNAL) See iterating functions .ellipsoidalBaseDI._intersect3,
272 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.
273 '''
274 t = self.reverse(x_t, *y) if y else self.reverse(x_t.x, x_t.y) # LatLon=None
275 d = euclid(t.lat - self.lat0, t.lon - self.lon0) # degrees
276 return t, d
278 def _toLatLon(self, lat, lon, LatLon, name_LatLon_kwds):
279 '''(INTERNAL) Check B{C{LatLon}} and return an instance.
280 '''
281 kwds = _xkwds(name_LatLon_kwds, datum=self.datum, _or_nameof=self)
282 r = LatLon(lat, lon, **kwds) # handle .classof
283 B = _LLEB if self.datum.isEllipsoidal else _LLB
284 _xinstanceof(B, LatLon=r)
285 return r
287 def toRepr(self, prec=6, **unused): # PYCHOK expected
288 '''Return a string representation of this projection.
290 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
292 @return: This projection as C{"<classname>(lat0, lon0, ...)"}
293 (C{str}).
294 '''
295 return _fstrLL0(self, prec, True)
297 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected
298 '''Return a string representation of this projection.
300 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
301 @kwarg sep: Separator to join (C{str}).
303 @return: This projection as C{"lat0 lon0"} (C{str}).
304 '''
305 t = _fstrLL0(self, prec, False)
306 return t if sep is None else sep.join(t)
309class AzimuthalError(_ValueError):
310 '''An azimuthal L{Equidistant}, L{EquidistantKarney}, L{Gnomonic},
311 L{LambertEqualArea}, L{Orthographic}, L{Stereographic} or
312 L{Azimuthal7Tuple} issue.
313 '''
314 pass
317class Azimuthal7Tuple(_NamedTuple):
318 '''7-Tuple C{(x, y, lat, lon, azimuth, scale, datum)}, in C{meter}, C{meter},
319 C{degrees90}, C{degrees180}, compass C{degrees}, C{scalar} and C{Datum}
320 where C{(x, y)} is the easting and northing of a projected point, C{(lat,
321 lon)} the geodetic location, C{azimuth} the azimuth, clockwise from true
322 North and C{scale} is the projection scale, either C{1 / reciprocal} or
323 C{1} or C{-1} in the L{Equidistant} case.
324 '''
325 _Names_ = (_x_, _y_, _lat_, _lon_, _azimuth_, _scale_, _datum_)
326 _Units_ = ( Easting, Northing, Lat_, Lon_, Bearing, Scalar, _Pass)
328 def antipodal(self, azimuth=None):
329 '''Return this tuple with the antipodal C{lat} and C{lon}.
331 @kwarg azimuth: Optional azimuth, overriding the current azimuth
332 (C{compass degrees360}).
333 '''
334 a = _MODS.formy.antipode(self.lat, self.lon) # PYCHOK named
335 z = self.azimuth if azimuth is None else Bearing(azimuth=azimuth) # PYCHOK named
336 return _NamedTuple.dup(self, lat=a.lat, lon=a.lon, azimuth=z)
339class Equidistant(_AzimuthalBase):
340 '''Azimuthal equidistant projection for the sphere***, see U{Snyder, pp 195-197
341 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
342 <https://MathWorld.Wolfram.com/AzimuthalEquidistantProjection.html>}.
344 @note: Results from this L{Equidistant} and an L{EquidistantExact},
345 L{EquidistantGeodSolve} or L{EquidistantKarney} projection
346 C{may differ} by 10% or more. For an example, see method
347 C{testDiscrepancies} in module C{testAzimuthal.py}.
348 '''
349 if _FOR_DOCS:
350 __init__ = _AzimuthalBase.__init__
352 def forward(self, lat, lon, **name):
353 '''Convert a geodetic location to azimuthal equidistant east- and northing.
355 @arg lat: Latitude of the location (C{degrees90}).
356 @arg lon: Longitude of the location (C{degrees180}).
357 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
359 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
360 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
361 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
362 The C{scale} of the projection is C{1} in I{radial} direction and
363 is C{1 / reciprocal} in the direction perpendicular to this.
365 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
367 @note: The C{scale} will be C{-1} if B{C{(lat, lon)}} is antipodal to
368 the projection center C{(lat0, lon0)}.
369 '''
370 def _k_t(c):
371 k = _N_1_0 if c < 0 else _1_0
372 t = fabs(c) < EPS1
373 if t:
374 a = acos(c)
375 s = sin(a)
376 if s:
377 k = a / s
378 return k, t
380 return self._forward(lat, lon, name, _k_t)
382 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):
383 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
385 @arg x: Easting of the location (C{meter}).
386 @arg y: Northing of the location (C{meter}).
387 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
388 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} forthe location and
389 optional, additional B{C{LatLon}} keyword arguments,
390 ignored if C{B{LatLon} is None}.
392 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
393 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
395 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
396 in the range C{[-180..180] degrees}. The C{scale} of the
397 projection is C{1} in I{radial} direction, C{azimuth} clockwise
398 from true North and is C{1 / reciprocal} in the direction
399 perpendicular to this.
400 '''
401 def _c(c):
402 return c if c > EPS else None
404 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)
407def equidistant(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, **name):
408 '''Return an L{EquidistantExact}, L{EquidistantGeodSolve} or (if I{Karney}'s
409 U{geographiclib<https://PyPI.org/project/geographiclib>} package is
410 installed) an L{EquidistantKarney}, otherwise an L{Equidistant} instance.
412 @arg lat0: Latitude of center point (C{degrees90}).
413 @arg lon0: Longitude of center point (C{degrees180}).
414 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
415 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
416 radius (C{meter}).
417 @kwarg exact: Return an L{EquidistantExact} instance.
418 @kwarg geodsolve: Return an L{EquidistantGeodSolve} instance.
419 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
421 @return: An L{EquidistantExact}, L{EquidistantGeodSolve},
422 L{EquidistantKarney} or L{Equidistant} instance.
424 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
426 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
427 or I{Karney}'s wrapped C{Geodesic}.
429 @raise TypeError: Invalid B{C{datum}}.
430 '''
432 E = EquidistantExact if exact else (EquidistantGeodSolve if geodsolve else Equidistant)
433 if E is Equidistant:
434 try:
435 return EquidistantKarney(lat0, lon0, datum=datum, **name) # PYCHOK types
436 except ImportError:
437 pass
438 return E(lat0, lon0, datum=datum, **name) # PYCHOK types
441class _AzimuthalGeodesic(_AzimuthalBase):
442 '''(INTERNAL) Base class for azimuthal projections using the
443 I{wrapped} U{geodesic.Geodesic and geodesicline.GeodesicLine
444 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} or the
445 I{exact} geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
446 '''
447 _mask = 0
449 @Property_RO
450 def geodesic(self): # PYCHOK no cover
451 '''I{Must be overloaded}.'''
452 self._notOverloaded()
454 def _7Tuple(self, e, n, r, name_LatLon_kwds, M=None):
455 '''(INTERNAL) Return an C{Azimuthal7Tuple}.
456 '''
457 s = M
458 if s is None: # reciprocal, azimuthal scale
459 s = (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0
460 z = _umod_360(r.azi2) # -180 <= r.azi2 < 180 ... 0 <= z < 360
461 t, _ = _name2__(name_LatLon_kwds, _or_nameof=self)
462 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum, name=t)
465class _EquidistantBase(_AzimuthalGeodesic):
466 '''(INTERNAL) Base for classes L{EquidistantExact}, L{EquidistantGeodSolve}
467 and L{EquidistantKarney}.
468 '''
469 def __init__(self, lat0, lon0, datum=_WGS84, **name):
470 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or
471 L{EquidistantKarney} projection.
472 '''
473 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, **name)
475 g = self.geodesic
476 # g.STANDARD = g.AZIMUTH | g.DISTANCE | g.LATITUDE | g.LONGITUDE
477 self._mask = g.REDUCEDLENGTH | g.STANDARD # | g.LONG_UNROLL
479 def forward(self, lat, lon, **name):
480 '''Convert an (ellipsoidal) geodetic location to azimuthal equidistant east- and northing.
482 @arg lat: Latitude of the location (C{degrees90}).
483 @arg lon: Longitude of the location (C{degrees180}).
484 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
486 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
487 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
488 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
489 The C{scale} of the projection is C{1} in I{radial} direction and
490 is C{1 / reciprocal} in the direction perpendicular to this.
492 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
494 @note: A call to C{.forward} followed by a call to C{.reverse} will return
495 the original C{lat, lon} to within roundoff.
496 '''
497 r = self.geodesic.Inverse(self.lat0, self.lon0,
498 Lat_(lat), Lon_(lon), outmask=self._mask)
499 x, y = sincos2d(r.azi1)
500 return self._7Tuple(x * r.s12, y * r.s12, r, _name__(name))
502 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK signature
503 '''Convert an azimuthal equidistant location to (ellipsoidal) geodetic lat- and longitude.
505 @arg x: Easting of the location (C{meter}).
506 @arg y: Northing of the location (C{meter}).
507 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
508 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and
509 optional, additional B{C{LatLon}} keyword arguments,
510 ignored if C{B{LatLon} is None}.
512 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
513 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
515 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
516 in the range C{[-180..180] degrees}. The scale of the projection
517 is C{1} in I{radial} direction, C{azimuth} clockwise from true
518 North and is C{1 / reciprocal} in the direction perpendicular
519 to this.
520 '''
521 e, n, z, s = _enzh4(x, y)
523 r = self.geodesic.Direct(self.lat0, self.lon0, z, s, outmask=self._mask)
524 return self._7Tuple(e, n, r, name_LatLon_kwds) if LatLon is None else \
525 self._toLatLon(r.lat2, r.lon2, LatLon, name_LatLon_kwds)
528class EquidistantExact(_EquidistantBase):
529 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
530 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
531 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
533 An azimuthal equidistant projection is centered at an arbitrary position on the ellipsoid.
534 For a point in projected space C{(x, y)}, the geodesic distance from the center position
535 is C{hypot(x, y)} and the C{azimuth} of the geodesic from the center point is C{atan2(x, y)},
536 clockwise from true North.
538 The C{.forward} and C{.reverse} methods also return the C{azimuth} of the geodesic at C{(x,
539 y)} and the C{scale} in the azimuthal direction which, together with the basic properties
540 of the projection, serve to specify completely the local affine transformation between
541 geographic and projected coordinates.
542 '''
543 def __init__(self, lat0, lon0, datum=_WGS84, **name):
544 '''New azimuthal L{EquidistantExact} projection.
546 @arg lat0: Latitude of center point (C{degrees90}).
547 @arg lon0: Longitude of center point (C{degrees180}).
548 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
549 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
550 radius (C{meter}).
551 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
553 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
554 '''
555 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name)
557 if _FOR_DOCS:
558 forward = _EquidistantBase.forward
559 reverse = _EquidistantBase.reverse
561 @Property_RO
562 def geodesic(self):
563 '''Get this projection's exact geodesic (L{GeodesicExact}).
564 '''
565 return self.datum.ellipsoid.geodesicx
568class EquidistantGeodSolve(_EquidistantBase):
569 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
570 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
571 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and intended
572 I{for testing purposes only}.
574 @see: L{EquidistantExact} and module L{geodsolve}.
575 '''
576 def __init__(self, lat0, lon0, datum=_WGS84, **name):
577 '''New azimuthal L{EquidistantGeodSolve} projection.
579 @arg lat0: Latitude of center point (C{degrees90}).
580 @arg lon0: Longitude of center point (C{degrees180}).
581 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
582 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
583 radius (C{meter}).
584 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
586 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
587 '''
588 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name)
590 if _FOR_DOCS:
591 forward = _EquidistantBase.forward
592 reverse = _EquidistantBase.reverse
594 @Property_RO
595 def geodesic(self):
596 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
597 '''
598 return self.datum.ellipsoid.geodsolve
601class EquidistantKarney(_EquidistantBase):
602 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
603 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
604 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
606 @see: L{EquidistantExact}.
607 '''
608 def __init__(self, lat0, lon0, datum=_WGS84, **name):
609 '''New azimuthal L{EquidistantKarney} projection.
611 @arg lat0: Latitude of center point (C{degrees90}).
612 @arg lon0: Longitude of center point (C{degrees180}).
613 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
614 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
615 radius (C{meter}).
616 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
618 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
619 not installed or not found.
621 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
622 '''
623 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, **name)
625 if _FOR_DOCS:
626 forward = _EquidistantBase.forward
627 reverse = _EquidistantBase.reverse
629 @Property_RO
630 def geodesic(self):
631 '''Get this projection's I{wrapped} U{geodesic.Geodesic
632 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
633 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
634 package is installed.
635 '''
636 return self.datum.ellipsoid.geodesic
639_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,
640 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI
643class Gnomonic(_AzimuthalBase):
644 '''Azimuthal gnomonic projection for the sphere***, see U{Snyder, pp 164-168
645 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
646 <https://MathWorld.Wolfram.com/GnomonicProjection.html>}.
647 '''
648 if _FOR_DOCS:
649 __init__ = _AzimuthalBase.__init__
651 def forward(self, lat, lon, **name):
652 '''Convert a geodetic location to azimuthal equidistant east- and northing.
654 @arg lat: Latitude of the location (C{degrees90}).
655 @arg lon: Longitude of the location (C{degrees180}).
656 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
658 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
659 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
660 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
661 The C{scale} of the projection is C{1} in I{radial} direction and
662 is C{1 / reciprocal} in the direction perpendicular to this.
664 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
665 '''
666 def _k_t(c):
667 t = c > EPS
668 k = (_1_0 / c) if t else _1_0
669 return k, t
671 return self._forward(lat, lon, name, _k_t)
673 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):
674 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
676 @arg x: Easting of the location (C{meter}).
677 @arg y: Northing of the location (C{meter}).
678 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
679 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and
680 optional, additional B{C{LatLon}} keyword arguments,
681 ignored if C{B{LatLon} is None}.
683 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
684 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
686 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
687 in the range C{[-180..180] degrees}. The C{scale} of the
688 projection is C{1} in I{radial} direction, C{azimuth} clockwise
689 from true North and C{1 / reciprocal} in the direction
690 perpendicular to this.
691 '''
692 def _c(c):
693 return atan1(c) if c > EPS else None
695 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)
698def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, **name):
699 '''Return a L{GnomonicExact} or (if I{Karney}'s U{geographiclib
700 <https://PyPI.org/project/geographiclib>} package is installed)
701 a L{GnomonicKarney}, otherwise a L{Gnomonic} instance.
703 @arg lat0: Latitude of center point (C{degrees90}).
704 @arg lon0: Longitude of center point (C{degrees180}).
705 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
706 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
707 radius (C{meter}).
708 @kwarg exact: Return a L{GnomonicExact} instance.
709 @kwarg geodsolve: Return a L{GnomonicGeodSolve} instance.
710 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
712 @return: A L{GnomonicExact}, L{GnomonicGeodSolve},
713 L{GnomonicKarney} or L{Gnomonic} instance.
715 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or
716 (spherical) B{C{datum}}.
718 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
719 or I{Karney}'s wrapped C{Geodesic}.
721 @raise TypeError: Invalid B{C{datum}}.
722 '''
723 G = GnomonicExact if exact else (GnomonicGeodSolve if geodsolve else Gnomonic)
724 if G is Gnomonic:
725 try:
726 return GnomonicKarney(lat0, lon0, datum=datum, **name) # PYCHOK types
727 except ImportError:
728 pass
729 return G(lat0, lon0, datum=datum, **name) # PYCHOK types
732class _GnomonicBase(_AzimuthalGeodesic):
733 '''(INTERNAL) Base for classes L{GnomonicExact}, L{GnomonicGeodSolve}
734 and L{GnomonicKarney}.
735 '''
736 def __init__(self, lat0, lon0, datum=_WGS84, **name):
737 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection.
738 '''
739 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, **name)
741 g = self.geodesic
742 self._mask = g.ALL # | g.LONG_UNROLL
744 def forward(self, lat, lon, raiser=True, **name): # PYCHOK signature
745 '''Convert an (ellipsoidal) geodetic location to azimuthal gnomonic east-
746 and northing.
748 @arg lat: Latitude of the location (C{degrees90}).
749 @arg lon: Longitude of the location (C{degrees180}).
750 @kwarg raiser: Do or don't throw an error (C{bool}) if
751 the location lies over the horizon.
752 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
754 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
755 with easting C{x} and northing C{y} in C{meter} and C{lat} and
756 C{lon} in C{degrees} and C{azimuth} clockwise from true North.
757 The C{scale} of the projection is C{1 / reciprocal**2} in I{radial}
758 direction and C{1 / reciprocal} in the direction perpendicular to
759 this. Both C{x} and C{y} will be C{NAN} if the (geodetic) location
760 lies over the horizon and C{B{raiser}=False}.
762 @raise AzimuthalError: Invalid B{C{lat}}, B{C{lon}} or the location lies
763 over the horizon and C{B{raiser}=True}.
764 '''
765 self._iteration = 0
767 r = self.geodesic.Inverse(self.lat0, self.lon0,
768 Lat_(lat), Lon_(lon), outmask=self._mask)
769 M = r.M21
770 if M > EPS0:
771 q = r.m12 / M # .M12
772 e, n = sincos2d(r.azi1)
773 e *= q
774 n *= q
775 elif raiser: # PYCHOK no cover
776 raise AzimuthalError(lat=lat, lon=lon, txt=_over_horizon_)
777 else: # PYCHOK no cover
778 e = n = NAN
780 t = self._7Tuple(e, n, r, _name__(name), M=M)
781 t._iteraton = 0
782 return t
784 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds): # PYCHOK signature
785 '''Convert an azimuthal gnomonic location to (ellipsoidal) geodetic lat- and longitude.
787 @arg x: Easting of the location (C{meter}).
788 @arg y: Northing of the location (C{meter}).
789 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
790 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and
791 optional, additional B{C{LatLon}} keyword arguments,
792 ignored if C{B{LatLon} is None}.
794 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
795 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
797 @raise AzimuthalError: No convergence.
799 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
800 in the range C{[-180..180] degrees}. The C{azimuth} is clockwise
801 from true North. The scale is C{1 / reciprocal**2} in C{radial}
802 direction and C{1 / reciprocal} in the direction perpendicular
803 to this.
804 '''
805 e, n, z, q = _enzh4(x, y)
807 d = a = self.equatoradius
808 s = a * atan1(q, a)
809 if q > a: # PYCHOK no cover
810 def _d(r, q):
811 return (r.M12 - q * r.m12) * r.m12 # negated
813 q = _1_0 / q
814 else: # little == True
815 def _d(r, q): # PYCHOK _d
816 return (q * r.M12 - r.m12) * r.M12 # negated
818 a *= _EPS_K
819 m = self._mask
820 g = self.geodesic
822 _P = g.Line(self.lat0, self.lon0, z, caps=m | g.LINE_OFF).Position
823 _S2 = Fsum(s).fsum2f_
824 _abs = fabs
825 for i in range(1, _TRIPS):
826 r = _P(s, outmask=m)
827 if _abs(d) < a:
828 break
829 s, d = _S2(_d(r, q))
830 else: # PYCHOK no cover
831 self._iteration = _TRIPS
832 raise AzimuthalError(Fmt.no_convergence(d, a),
833 txt=unstr(self.reverse, x, y))
835 t = self._7Tuple(e, n, r, name_LatLon_kwds, M=r.M12) if LatLon is None else \
836 self._toLatLon(r.lat2, r.lon2, LatLon, name_LatLon_kwds)
837 t._iteration = self._iteration = i
838 return t
841class GnomonicExact(_GnomonicBase):
842 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
843 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
844 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
846 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/
847 classGeographicLib_1_1Gnomonic.html>}, especially the B{Warning}.
848 '''
849 def __init__(self, lat0, lon0, datum=_WGS84, **name):
850 '''New azimuthal L{GnomonicExact} projection.
852 @arg lat0: Latitude of center point (C{degrees90}).
853 @arg lon0: Longitude of center point (C{degrees180}).
854 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
855 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
856 radius (C{meter}).
857 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
859 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
860 '''
861 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, **name)
863 if _FOR_DOCS:
864 forward = _GnomonicBase.forward
865 reverse = _GnomonicBase.reverse
867 @Property_RO
868 def geodesic(self):
869 '''Get this projection's exact geodesic (L{GeodesicExact}).
870 '''
871 return self.datum.ellipsoid.geodesicx
874class GnomonicGeodSolve(_GnomonicBase):
875 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
876 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
877 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and
878 intended I{for testing purposes only}.
880 @see: L{GnomonicExact} and module L{geodsolve}.
881 '''
882 def __init__(self, lat0, lon0, datum=_WGS84, **name):
883 '''New azimuthal L{GnomonicGeodSolve} projection.
885 @arg lat0: Latitude of center point (C{degrees90}).
886 @arg lon0: Longitude of center point (C{degrees180}).
887 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
888 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
889 radius (C{meter}).
890 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
892 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
893 '''
894 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, **name)
896 if _FOR_DOCS:
897 forward = _GnomonicBase.forward
898 reverse = _GnomonicBase.reverse
900 @Property_RO
901 def geodesic(self):
902 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
903 '''
904 return self.datum.ellipsoid.geodsolve
907class GnomonicKarney(_GnomonicBase):
908 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
909 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
910 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
912 @see: L{GnomonicExact}.
913 '''
914 def __init__(self, lat0, lon0, datum=_WGS84, **name):
915 '''New azimuthal L{GnomonicKarney} projection.
917 @arg lat0: Latitude of center point (C{degrees90}).
918 @arg lon0: Longitude of center point (C{degrees180}).
919 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
920 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
921 radius (C{meter}).
922 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).
924 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
925 not installed or not found.
927 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
928 '''
929 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, **name)
931 if _FOR_DOCS:
932 forward = _GnomonicBase.forward
933 reverse = _GnomonicBase.reverse
935 @Property_RO
936 def geodesic(self):
937 '''Get this projection's I{wrapped} U{geodesic.Geodesic
938 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
939 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
940 package is installed.
941 '''
942 return self.datum.ellipsoid.geodesic
945class LambertEqualArea(_AzimuthalBase):
946 '''Lambert-equal-area projection for the sphere*** (aka U{Lambert zenithal equal-area
947 projection<https://WikiPedia.org/wiki/Lambert_azimuthal_equal-area_projection>}, see
948 U{Snyder, pp 185-187<https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
949 <https://MathWorld.Wolfram.com/LambertAzimuthalEqual-AreaProjection.html>}.
950 '''
951 if _FOR_DOCS:
952 __init__ = _AzimuthalBase.__init__
954 def forward(self, lat, lon, **name):
955 '''Convert a geodetic location to azimuthal Lambert-equal-area east- and northing.
957 @arg lat: Latitude of the location (C{degrees90}).
958 @arg lon: Longitude of the location (C{degrees180}).
959 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
961 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
962 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
963 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
964 The C{scale} of the projection is C{1} in I{radial} direction and
965 is C{1 / reciprocal} in the direction perpendicular to this.
967 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
968 '''
969 def _k_t(c):
970 c += _1_0
971 t = c > EPS0
972 k = sqrt(_2_0 / c) if t else _1_0
973 return k, t
975 return self._forward(lat, lon, name, _k_t)
977 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):
978 '''Convert an azimuthal Lambert-equal-area location to geodetic lat- and longitude.
980 @arg x: Easting of the location (C{meter}).
981 @arg y: Northing of the location (C{meter}).
982 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
983 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location
984 and optional, additional B{C{LatLon}} keyword
985 arguments, ignored if C{B{LatLon} is None}.
987 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
988 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
990 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
991 in the range C{[-180..180] degrees}. The C{scale} of the
992 projection is C{1} in I{radial} direction, C{azimuth} clockwise
993 from true North and is C{1 / reciprocal} in the direction
994 perpendicular to this.
995 '''
996 def _c(c):
997 c *= _0_5
998 return (asin1(c) * _2_0) if c > EPS else None
1000 return self._reverse(x, y, _c, True, LatLon, **name_LatLon_kwds)
1003class Orthographic(_AzimuthalBase):
1004 '''Orthographic projection for the sphere***, see U{Snyder, pp 148-153
1005 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1006 <https://MathWorld.Wolfram.com/OrthographicProjection.html>}.
1007 '''
1008 if _FOR_DOCS:
1009 __init__ = _AzimuthalBase.__init__
1011 def forward(self, lat, lon, **name):
1012 '''Convert a geodetic location to azimuthal orthographic east- and northing.
1014 @arg lat: Latitude of the location (C{degrees90}).
1015 @arg lon: Longitude of the location (C{degrees180}).
1016 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
1018 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1019 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1020 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1021 The C{scale} of the projection is C{1} in I{radial} direction and
1022 is C{1 / reciprocal} in the direction perpendicular to this.
1024 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1025 '''
1026 def _k_t(c):
1027 return _1_0, (c >= 0)
1029 return self._forward(lat, lon, name, _k_t)
1031 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):
1032 '''Convert an azimuthal orthographic location to geodetic lat- and longitude.
1034 @arg x: Easting of the location (C{meter}).
1035 @arg y: Northing of the location (C{meter}).
1036 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
1037 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and
1038 optional, additional B{C{LatLon}} keyword arguments,
1039 ignored if C{B{LatLon} is None}.
1041 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1042 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1044 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1045 in the range C{[-180..180] degrees}. The C{scale} of the
1046 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1047 from true North and is C{1 / reciprocal} in the direction
1048 perpendicular to this.
1049 '''
1050 def _c(c):
1051 return asin1(c) if c > EPS else None
1053 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)
1056class Stereographic(_AzimuthalBase):
1057 '''Stereographic projection for the sphere***, see U{Snyder, pp 157-160
1058 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1059 <https://MathWorld.Wolfram.com/StereographicProjection.html>}.
1060 '''
1061 _k0 = _1_0 # central scale factor (C{scalar})
1062 _k02 = _2_0 # double ._k0
1064 if _FOR_DOCS:
1065 __init__ = _AzimuthalBase.__init__
1067 def forward(self, lat, lon, **name):
1068 '''Convert a geodetic location to azimuthal stereographic east- and northing.
1070 @arg lat: Latitude of the location (C{degrees90}).
1071 @arg lon: Longitude of the location (C{degrees180}).
1072 @kwarg name: Optional C{B{name}=NN} for the location (C{str}).
1074 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1075 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1076 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1077 The C{scale} of the projection is C{1} in I{radial} direction and
1078 is C{1 / reciprocal} in the direction perpendicular to this.
1080 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1081 '''
1082 def _k_t(c):
1083 c += _1_0
1084 t = isnon0(c)
1085 k = (self._k02 / c) if t else _1_0
1086 return k, t
1088 return self._forward(lat, lon, name, _k_t)
1090 @property_doc_(''' optional, central scale factor (C{scalar}).''')
1091 def k0(self):
1092 '''Get the central scale factor (C{scalar}).
1093 '''
1094 return self._k0
1096 @k0.setter # PYCHOK setter!
1097 def k0(self, factor):
1098 '''Set the central scale factor (C{scalar}).
1099 '''
1100 n = Stereographic.k0.fget.__name__ # 'k0', name__=Stereographic.k0.fget
1101 self._k0 = Scalar_(factor, name=n, low=EPS, high=2) # XXX high=1, 2, other?
1102 self._k02 = self._k0 * _2_0
1104 def reverse(self, x, y, LatLon=None, **name_LatLon_kwds):
1105 '''Convert an azimuthal stereographic location to geodetic lat- and longitude.
1107 @arg x: Easting of the location (C{meter}).
1108 @arg y: Northing of the location (C{meter}).
1109 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
1110 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} for the location and
1111 optional, additional B{C{LatLon}} keyword arguments,
1112 ignored if C{B{LatLon} is None}.
1114 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1115 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1117 @note: The C{lat} will be in range C{[-90..90] degrees}, C{lon} in range
1118 C{[-180..180] degrees} and C{azimuth} clockwise from true North.
1119 The C{scale} of the projection is C{1} in I{radial} direction and
1120 is C{1 / reciprocal} in the direction perpendicular to this.
1121 '''
1122 def _c(c):
1123 return (atan2(c, self._k02) * _2_0) if c > EPS else None
1125 return self._reverse(x, y, _c, False, LatLon, **name_LatLon_kwds)
1128__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase)
1130# **) MIT License
1131#
1132# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1133#
1134# Permission is hereby granted, free of charge, to any person obtaining a
1135# copy of this software and associated documentation files (the "Software"),
1136# to deal in the Software without restriction, including without limitation
1137# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1138# and/or sell copies of the Software, and to permit persons to whom the
1139# Software is furnished to do so, subject to the following conditions:
1140#
1141# The above copyright notice and this permission notice shall be included
1142# in all copies or substantial portions of the Software.
1143#
1144# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1145# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1146# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1147# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1148# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1149# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1150# OTHER DEALINGS IN THE SOFTWARE.