Coverage for pygeodesy/azimuthal.py: 99%
316 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -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, \
47 _EPStol, _umod_360, _0_0, _0_1, \
48 _0_5, _1_0, _N_1_0, _2_0
49from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB, \
50 _xinstanceof
51from pygeodesy.datums import _spherical_datum, _WGS84
52from pygeodesy.errors import _ValueError, _xdatum, _xkwds
53from pygeodesy.fmath import euclid, hypot as _hypot, Fsum
54# from pygeodesy.fsums import Fsum # from .fmath
55# from pygeodesy.formy import antipode # _MODS
56from pygeodesy.interns import NN, _azimuth_, _datum_, _lat_, _lon_, \
57 _scale_, _SPACE_, _x_, _y_
58from pygeodesy.karney import _norm180
59from pygeodesy.latlonBase import _MODS, LatLonBase as _LLB
60from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS # ALL_MODS
61from pygeodesy.named import _NamedBase, _NamedTuple, _Pass
62from pygeodesy.namedTuples import LatLon2Tuple, LatLon4Tuple
63from pygeodesy.props import deprecated_Property_RO, Property_RO, \
64 property_doc_, _update_all
65from pygeodesy.streprs import Fmt, _fstrLL0, unstr
66from pygeodesy.units import Bearing, Easting, Lat_, Lon_, Northing, \
67 Scalar, Scalar_
68from pygeodesy.utily import asin1, atan1, atan2b, atan2d, sincos2, \
69 sincos2d, sincos2d_
71from math import acos, atan2, degrees, fabs, sin, sqrt
73__all__ = _ALL_LAZY.azimuthal
74__version__ = '24.04.07'
76_EPS_K = _EPStol * _0_1 # Karney's eps_ or _EPSmin * _0_1?
77_over_horizon_ = 'over horizon'
78_TRIPS = 21 # numit, 4 sufficient
81def _enzh4(x, y, *h):
82 '''(INTERNAL) Return 4-tuple (easting, northing, azimuth, hypot).
83 '''
84 e = Easting( x=x)
85 n = Northing(y=y)
86 z = atan2b(e, n) # (x, y) for azimuth from true North
87 return e, n, z, (h[0] if h else _hypot(e, n))
90class _AzimuthalBase(_NamedBase):
91 '''(INTERNAL) Base class for azimuthal projections.
93 @see: I{Karney}'s C++ class U{AzimuthalEquidistant<https://GeographicLib.SourceForge.io/
94 C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>} and U{Gnomonic
95 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>} or the
96 C{PyGeodesy} versions thereof L{EquidistantKarney} respectively L{GnomonicKarney}.
97 '''
98 _datum = _WGS84 # L{Datum}
99 _latlon0 = LatLon2Tuple(_0_0, _0_0) # lat0, lon0 (L{LatLon2Tuple})
100 _sc0 = _0_0, _1_0 # 2-Tuple C{sincos2d(lat0)}
102 def __init__(self, lat0, lon0, datum=None, name=NN):
103 '''New azimuthal projection.
105 @arg lat0: Latitude of the center point (C{degrees90}).
106 @arg lon0: Longitude of the center point (C{degrees180}).
107 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
108 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
109 radius (C{meter}).
110 @kwarg name: Optional name for the projection (C{str}).
112 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
114 @raise TypeError: Invalid B{C{datum}}.
115 '''
116 if datum not in (None, self._datum):
117 self._datum = _spherical_datum(datum, name=name)
118 if name:
119 self.name = name
121 if lat0 or lon0: # often both 0
122 self._reset(lat0, lon0)
124 @Property_RO
125 def datum(self):
126 '''Get the datum (L{Datum}).
127 '''
128 return self._datum
130 @Property_RO
131 def equatoradius(self):
132 '''Get the geodesic's equatorial radius, semi-axis (C{meter}).
133 '''
134 return self.datum.ellipsoid.a
136 a = equatoradius
138 @Property_RO
139 def flattening(self):
140 '''Get the geodesic's flattening (C{scalar}).
141 '''
142 return self.datum.ellipsoid.f
144 f = flattening
146 def forward(self, lat, lon, name=NN): # PYCHOK no cover
147 '''I{Must be overloaded}.'''
148 self._notOverloaded(lat, lon, name=name)
150 def _forward(self, lat, lon, name, _k_t_2):
151 '''(INTERNAL) Azimuthal (spherical) forward C{lat, lon} to C{x, y}.
152 '''
153 lat, lon = Lat_(lat), Lon_(lon)
154 sa, ca, sb, cb = sincos2d_(lat, lon - self.lon0)
155 s0, c0 = self._sc0
157 cb *= ca
158 k, t = _k_t_2(s0 * sa + c0 * cb)
159 if t:
160 r = k * self.radius
161 y = r * (c0 * sa - s0 * cb)
162 e, n, z, _ = _enzh4(r * sb * ca, y, None)
163 else: # 0 or 180
164 e = n = z = _0_0
166 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
167 name=name or self.name)
168 return t
170 def _forwards(self, *lls):
171 '''(INTERNAL) One or more C{.forward} calls, see .ellipsoidalBaseDI.
172 '''
173 _fwd = self.forward
174 for ll in lls:
175 yield _fwd(ll.lat, ll.lon)
177 @Property_RO
178 def lat0(self):
179 '''Get the center latitude (C{degrees90}).
180 '''
181 return self._latlon0.lat
183 @property
184 def latlon0(self):
185 '''Get the center lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}) in (C{degrees90}, C{degrees180}).
186 '''
187 return self._latlon0
189 @latlon0.setter # PYCHOK setter!
190 def latlon0(self, latlon0):
191 '''Set the center lat- and longitude (C{LatLon}, L{LatLon2Tuple} or L{LatLon4Tuple}).
193 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or ellipsoidal mismatch
194 of B{C{latlon0}} and this projection.
195 '''
196 B = _LLEB if self.datum.isEllipsoidal else _LLB
197 _xinstanceof(B, LatLon2Tuple, LatLon4Tuple, latlon0=latlon0)
198 if hasattr(latlon0, _datum_):
199 _xdatum(self.datum, latlon0.datum, Error=AzimuthalError)
200 self.reset(latlon0.lat, latlon0.lon)
202 @Property_RO
203 def lon0(self):
204 '''Get the center longitude (C{degrees180}).
205 '''
206 return self._latlon0.lon
208 @deprecated_Property_RO
209 def majoradius(self): # PYCHOK no cover
210 '''DEPRECATED, use property C{equatoradius}.'''
211 return self.equatoradius
213 @Property_RO
214 def radius(self):
215 '''Get this projection's mean radius of curvature (C{meter}).
216 '''
217 return self.datum.ellipsoid.rocMean(self.lat0)
219 def reset(self, lat0, lon0):
220 '''Set or reset the center point of this azimuthal projection.
222 @arg lat0: Center point latitude (C{degrees90}).
223 @arg lon0: Center point longitude (C{degrees180}).
225 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
226 '''
227 _update_all(self) # zap caches
228 self._reset(lat0, lon0)
230 def _reset(self, lat0, lon0):
231 '''(INTERNAL) Update the center point.
232 '''
233 self._latlon0 = LatLon2Tuple(Lat_(lat0=lat0, Error=AzimuthalError),
234 Lon_(lon0=lon0, Error=AzimuthalError))
235 self._sc0 = sincos2d(self.lat0)
237 def reverse(self, x, y, name=NN, **LatLon_and_kwds): # PYCHOK no cover
238 '''I{Must be overloaded}.'''
239 self._notOverloaded(x, y, name=name, **LatLon_and_kwds)
241 def _reverse(self, x, y, name, _c, lea, LatLon=None, **LatLon_kwds):
242 '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
243 '''
244 e, n, z, r = _enzh4(x, y)
246 c = _c(r / self.radius)
247 if c is None:
248 lat, lon = self.latlon0
249 k, z = _1_0, _0_0
250 else:
251 s0, c0 = self._sc0
252 sc, cc = sincos2(c)
253 k = c / sc
254 s = s0 * cc
255 if r > EPS0:
256 s += c0 * sc * (n / r)
257 lat = degrees(asin1(s))
258 if lea or fabs(c0) > EPS:
259 d = atan2d(e * sc, r * c0 * cc - n * s0 * sc)
260 else:
261 d = atan2d(e, (n if s0 < 0 else -n))
262 lon = _norm180(self.lon0 + d)
264 if LatLon is None:
265 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
266 name=name or self.name)
267 else:
268 t = self._toLatLon(lat, lon, LatLon, LatLon_kwds, name)
269 return t
271 def _reverse2(self, x_t, *y):
272 '''(INTERNAL) See iterating functions .ellipsoidalBaseDI._intersect3,
273 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.
274 '''
275 t = self.reverse(x_t, *y) if y else self.reverse(x_t.x, x_t.y) # LatLon=None
276 d = euclid(t.lat - self.lat0, t.lon - self.lon0) # degrees
277 return t, d
279 def _toLatLon(self, lat, lon, LatLon, LatLon_kwds, name):
280 '''(INTERNAL) Check B{C{LatLon}} and return an instance.
281 '''
282 kwds = _xkwds(LatLon_kwds, datum=self.datum)
283 r = self._xnamed(LatLon(lat, lon, **kwds), name=name) # handle .classof
284 B = _LLEB if self.datum.isEllipsoidal else _LLB
285 _xinstanceof(B, LatLon=r)
286 return r
288 def toRepr(self, prec=6, **unused): # PYCHOK expected
289 '''Return a string representation of this projection.
291 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
293 @return: This projection as C{"<classname>(lat0, lon0, ...)"}
294 (C{str}).
295 '''
296 return _fstrLL0(self, prec, True)
298 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected
299 '''Return a string representation of this projection.
301 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
302 @kwarg sep: Separator to join (C{str}).
304 @return: This projection as C{"lat0 lon0"} (C{str}).
305 '''
306 t = _fstrLL0(self, prec, False)
307 return t if sep is None else sep.join(t)
310class AzimuthalError(_ValueError):
311 '''An azimuthal L{Equidistant}, L{EquidistantKarney}, L{Gnomonic},
312 L{LambertEqualArea}, L{Orthographic}, L{Stereographic} or
313 L{Azimuthal7Tuple} issue.
314 '''
315 pass
318class Azimuthal7Tuple(_NamedTuple):
319 '''7-Tuple C{(x, y, lat, lon, azimuth, scale, datum)}, in C{meter}, C{meter},
320 C{degrees90}, C{degrees180}, compass C{degrees}, C{scalar} and C{Datum}
321 where C{(x, y)} is the easting and northing of a projected point, C{(lat,
322 lon)} the geodetic location, C{azimuth} the azimuth, clockwise from true
323 North and C{scale} is the projection scale, either C{1 / reciprocal} or
324 C{1} or C{-1} in the L{Equidistant} case.
325 '''
326 _Names_ = (_x_, _y_, _lat_, _lon_, _azimuth_, _scale_, _datum_)
327 _Units_ = ( Easting, Northing, Lat_, Lon_, Bearing, Scalar, _Pass)
329 def antipodal(self, azimuth=None):
330 '''Return this tuple with the antipodal C{lat} and C{lon}.
332 @kwarg azimuth: Optional azimuth, overriding the current azimuth
333 (C{compass degrees360}).
334 '''
335 a = _MODS.formy.antipode(self.lat, self.lon) # PYCHOK named
336 z = self.azimuth if azimuth is None else Bearing(azimuth=azimuth) # PYCHOK named
337 return _NamedTuple.dup(self, lat=a.lat, lon=a.lon, azimuth=z)
340class Equidistant(_AzimuthalBase):
341 '''Azimuthal equidistant projection for the sphere***, see U{Snyder, pp 195-197
342 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
343 <https://MathWorld.Wolfram.com/AzimuthalEquidistantProjection.html>}.
345 @note: Results from this L{Equidistant} and an L{EquidistantExact},
346 L{EquidistantGeodSolve} or L{EquidistantKarney} projection
347 C{may differ} by 10% or more. For an example, see method
348 C{testDiscrepancies} in module C{testAzimuthal.py}.
349 '''
350 if _FOR_DOCS:
351 __init__ = _AzimuthalBase.__init__
353 def forward(self, lat, lon, name=NN):
354 '''Convert a geodetic location to azimuthal equidistant east- and northing.
356 @arg lat: Latitude of the location (C{degrees90}).
357 @arg lon: Longitude of the location (C{degrees180}).
358 @kwarg name: Optional name for the location (C{str}).
360 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
361 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
362 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
363 The C{scale} of the projection is C{1} in I{radial} direction and
364 is C{1 / reciprocal} in the direction perpendicular to this.
366 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
368 @note: The C{scale} will be C{-1} if B{C{(lat, lon)}} is antipodal to
369 the projection center C{(lat0, lon0)}.
370 '''
371 def _k_t(c):
372 k = _N_1_0 if c < 0 else _1_0
373 t = fabs(c) < EPS1
374 if t:
375 a = acos(c)
376 s = sin(a)
377 if s:
378 k = a / s
379 return k, t
381 return self._forward(lat, lon, name, _k_t)
383 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
384 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
386 @arg x: Easting of the location (C{meter}).
387 @arg y: Northing of the location (C{meter}).
388 @kwarg name: Optional name for the location (C{str}).
389 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
390 additional B{C{LatLon}} keyword arguments,
391 ignored if C{B{LatLon} is None} or not given.
393 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
394 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
396 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
397 in the range C{[-180..180] degrees}. The C{scale} of the
398 projection is C{1} in I{radial} direction, C{azimuth} clockwise
399 from true North and is C{1 / reciprocal} in the direction
400 perpendicular to this.
401 '''
402 def _c(c):
403 return c if c > EPS else None
405 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
408def equidistant(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
409 '''Return an L{EquidistantExact}, L{EquidistantGeodSolve} or (if I{Karney}'s
410 U{geographiclib<https://PyPI.org/project/geographiclib>} package is
411 installed) an L{EquidistantKarney}, otherwise an L{Equidistant} instance.
413 @arg lat0: Latitude of center point (C{degrees90}).
414 @arg lon0: Longitude of center point (C{degrees180}).
415 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
416 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
417 radius (C{meter}).
418 @kwarg exact: Return an L{EquidistantExact} instance.
419 @kwarg geodsolve: Return an L{EquidistantGeodSolve} instance.
420 @kwarg name: Optional name for the projection (C{str}).
422 @return: An L{EquidistantExact}, L{EquidistantGeodSolve},
423 L{EquidistantKarney} or L{Equidistant} instance.
425 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
427 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
428 or I{Karney}'s wrapped C{Geodesic}.
430 @raise TypeError: Invalid B{C{datum}}.
431 '''
433 E = EquidistantExact if exact else (EquidistantGeodSolve if geodsolve else Equidistant)
434 if E is Equidistant:
435 try:
436 return EquidistantKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
437 except ImportError:
438 pass
439 return E(lat0, lon0, datum=datum, name=name) # PYCHOK types
442class _AzimuthalGeodesic(_AzimuthalBase):
443 '''(INTERNAL) Base class for azimuthal projections using the
444 I{wrapped} U{geodesic.Geodesic and geodesicline.GeodesicLine
445 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} or the
446 I{exact} geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
447 '''
448 _mask = 0
450 @Property_RO
451 def geodesic(self): # PYCHOK no cover
452 '''I{Must be overloaded}.'''
453 self._notOverloaded()
455 def _7Tuple(self, e, n, r, M=None, name=NN):
456 '''(INTERNAL) Return an C{Azimuthal7Tuple}.
457 '''
458 s = M if M is not None else ( # reciprocal, azimuthal scale
459 (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 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum,
462 name=name or self.name)
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=NN):
470 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or
471 L{EquidistantKarney} projection.
472 '''
473 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=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=NN):
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 name 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, name=NN, LatLon=None, **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 name: Optional name for the location (C{str}).
508 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
509 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
510 arguments, 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=name) if LatLon is None else \
525 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
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=NN):
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 name 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=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=NN):
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 name 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=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=NN):
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 name 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=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=NN):
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 name 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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).
679 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
680 additional B{C{LatLon}} keyword arguments,
681 ignored if C{B{LatLon} is None} or not given.
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, name, _c, False, **LatLon_and_kwds)
698def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
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 name 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=name) # PYCHOK types
727 except ImportError:
728 pass
729 return G(lat0, lon0, datum=datum, name=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=NN):
737 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection.
738 '''
739 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
741 g = self.geodesic
742 self._mask = g.ALL # | g.LONG_UNROLL
744 def forward(self, lat, lon, name=NN, raiser=True): # 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 name: Optional name for the location (C{str}).
751 @kwarg raiser: Do or don't throw an error (C{bool}) if
752 the location lies over the horizon.
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, M=M, name=name)
781 t._iteraton = 0
782 return t
784 def reverse(self, x, y, name=NN, LatLon=None, **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 name: Optional name for the location (C{str}).
790 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
791 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
792 arguments, 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, M=r.M12, name=name) if LatLon is None else \
836 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
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=NN):
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 name 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=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=NN):
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 name 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=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=NN):
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 name 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=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=NN):
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 name 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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).
983 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
984 additional B{C{LatLon}} keyword arguments,
985 ignored if C{B{LatLon} is None} or not given.
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, name, _c, True, **LatLon_and_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=NN):
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 name 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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).
1037 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1038 additional B{C{LatLon}} keyword arguments,
1039 ignored if C{B{LatLon} is None} or not given.
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, name, _c, False, **LatLon_and_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=NN):
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 name 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__
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, name=NN, **LatLon_and_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 name: Optional name for the location (C{str}).
1110 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1111 additional B{C{LatLon}} keyword arguments,
1112 ignored if C{B{LatLon} is None} or not given.
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 the range C{[-90..90] degrees} and C{lon}
1118 in the range C{[-180..180] degrees}. The C{scale} of the
1119 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1120 from true North and is C{1 / reciprocal} in the direction
1121 perpendicular to this.
1122 '''
1123 def _c(c):
1124 return (atan2(c, self._k02) * _2_0) if c > EPS else None
1126 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1129__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase)
1131# **) MIT License
1132#
1133# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1134#
1135# Permission is hereby granted, free of charge, to any person obtaining a
1136# copy of this software and associated documentation files (the "Software"),
1137# to deal in the Software without restriction, including without limitation
1138# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1139# and/or sell copies of the Software, and to permit persons to whom the
1140# Software is furnished to do so, subject to the following conditions:
1141#
1142# The above copyright notice and this permission notice shall be included
1143# in all copies or substantial portions of the Software.
1144#
1145# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1146# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1147# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1148# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1149# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1150# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1151# OTHER DEALINGS IN THE SOFTWARE.