Coverage for pygeodesy/azimuthal.py: 98%
314 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-20 13:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-20 13: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 .datums
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 _spherical_datum, _WGS84, _xinstanceof
51# from pygeodesy.datums import _spherical_datum, _WGS84, _xinstanceof
52from pygeodesy.errors import _ValueError, _xdatum, _xkwds
53from pygeodesy.fmath import euclid, Fsum, hypot
54# from pygeodesy.fsums import Fsum # from .fmath
55# from pygeodesy.formy import antipode # from latlonBase
56from pygeodesy.interns import NN, _azimuth_, _datum_, _lat_, _lon_, \
57 _scale_, _SPACE_, _x_, _y_
58from pygeodesy.karney import _norm180
59from pygeodesy.latlonBase import antipode, LatLonBase as _LLB
60from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS
61from pygeodesy.named import _NamedBase, _NamedTuple, notOverloaded, _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_, \
67 Northing, Scalar, Scalar_
68from pygeodesy.utily import asin1, atan2b, atan2d, sincos2, \
69 sincos2d, sincos2d_
71from math import acos, atan, atan2, degrees, fabs, sin, sqrt
73__all__ = _ALL_LAZY.azimuthal
74__version__ = '23.09.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 html/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 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
148 '''
149 notOverloaded(self, lat, lon, name=name)
151 def _forward(self, lat, lon, name, _k_t_2):
152 '''(INTERNAL) Azimuthal (spherical) forward C{lat, lon} to C{x, y}.
153 '''
154 lat, lon = Lat_(lat), Lon_(lon)
155 sa, ca, sb, cb = sincos2d_(lat, lon - self.lon0)
156 s0, c0 = self._sc0
158 cb *= ca
159 k, t = _k_t_2(s0 * sa + c0 * cb)
160 if t:
161 r = k * self.radius
162 y = r * (c0 * sa - s0 * cb)
163 e, n, z, _ = _enzh4(r * sb * ca, y, None)
164 else: # 0 or 180
165 e = n = z = _0_0
167 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
168 name=name or self.name)
169 return t
171 def _forwards(self, *lls):
172 '''(INTERNAL) One or more C{.forward} calls, see .ellipsoidalBaseDI.
173 '''
174 _fwd = self.forward
175 for ll in lls:
176 yield _fwd(ll.lat, ll.lon)
178 @Property_RO
179 def lat0(self):
180 '''Get the center latitude (C{degrees90}).
181 '''
182 return self._latlon0.lat
184 @property
185 def latlon0(self):
186 '''Get the center lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}) in (C{degrees90}, C{degrees180}).
187 '''
188 return self._latlon0
190 @latlon0.setter # PYCHOK setter!
191 def latlon0(self, latlon0):
192 '''Set the center lat- and longitude (C{LatLon}, L{LatLon2Tuple} or L{LatLon4Tuple}).
194 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or ellipsoidal mismatch
195 of B{C{latlon0}} and this projection.
196 '''
197 B = _LLEB if self.datum.isEllipsoidal else _LLB
198 _xinstanceof(B, LatLon2Tuple, LatLon4Tuple, latlon0=latlon0)
199 if hasattr(latlon0, _datum_):
200 _xdatum(self.datum, latlon0.datum, Error=AzimuthalError)
201 self.reset(latlon0.lat, latlon0.lon)
203 @Property_RO
204 def lon0(self):
205 '''Get the center longitude (C{degrees180}).
206 '''
207 return self._latlon0.lon
209 @deprecated_Property_RO
210 def majoradius(self): # PYCHOK no cover
211 '''DEPRECATED, use property C{equatoradius}.'''
212 return self.equatoradius
214 @Property_RO
215 def radius(self):
216 '''Get this projection's mean radius of curvature (C{meter}).
217 '''
218 return self.datum.ellipsoid.rocMean(self.lat0)
220 def reset(self, lat0, lon0):
221 '''Set or reset the center point of this azimuthal projection.
223 @arg lat0: Center point latitude (C{degrees90}).
224 @arg lon0: Center point longitude (C{degrees180}).
226 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
227 '''
228 _update_all(self) # zap caches
229 self._reset(lat0, lon0)
231 def _reset(self, lat0, lon0):
232 '''(INTERNAL) Update the center point.
233 '''
234 self._latlon0 = LatLon2Tuple(Lat_(lat0=lat0, Error=AzimuthalError),
235 Lon_(lon0=lon0, Error=AzimuthalError))
236 self._sc0 = sincos2d(self.lat0)
238 def reverse(self, x, y, name=NN, **LatLon_and_kwds): # PYCHOK no cover
239 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
240 '''
241 notOverloaded(self, x, y, name=name, **LatLon_and_kwds)
243 def _reverse(self, x, y, name, _c, lea, LatLon=None, **LatLon_kwds):
244 '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
245 '''
246 e, n, z, r = _enzh4(x, y)
248 c = _c(r / self.radius)
249 if c is None:
250 lat, lon = self.latlon0
251 k, z = _1_0, _0_0
252 else:
253 s0, c0 = self._sc0
254 sc, cc = sincos2(c)
255 k = c / sc
256 s = s0 * cc
257 if r > EPS0:
258 s += c0 * sc * (n / r)
259 lat = degrees(asin1(s))
260 if lea or fabs(c0) > EPS:
261 d = atan2d(e * sc, r * c0 * cc - n * s0 * sc)
262 else:
263 d = atan2d(e, (n if s0 < 0 else -n))
264 lon = _norm180(self.lon0 + d)
266 if LatLon is None:
267 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
268 name=name or self.name)
269 else:
270 t = self._toLatLon(lat, lon, LatLon, LatLon_kwds, name)
271 return t
273 def _reverse2(self, x_t, *y):
274 '''(INTERNAL) See iterating functions .ellipsoidalBaseDI._intersect3,
275 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.
276 '''
277 t = self.reverse(x_t, *y) if y else self.reverse(x_t.x, x_t.y) # LatLon=None
278 d = euclid(t.lat - self.lat0, t.lon - self.lon0) # degrees
279 return t, d
281 def _toLatLon(self, lat, lon, LatLon, LatLon_kwds, name):
282 '''(INTERNAL) Check B{C{LatLon}} and return an instance.
283 '''
284 kwds = _xkwds(LatLon_kwds, datum=self.datum)
285 r = self._xnamed(LatLon(lat, lon, **kwds), name=name) # handle .classof
286 B = _LLEB if self.datum.isEllipsoidal else _LLB
287 _xinstanceof(B, LatLon=r)
288 return r
290 def toRepr(self, prec=6, **unused): # PYCHOK expected
291 '''Return a string representation of this projection.
293 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
295 @return: This projection as C{"<classname>(lat0, lon0, ...)"}
296 (C{str}).
297 '''
298 return _fstrLL0(self, prec, True)
300 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected
301 '''Return a string representation of this projection.
303 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
304 @kwarg sep: Separator to join (C{str}).
306 @return: This projection as C{"lat0 lon0"} (C{str}).
307 '''
308 t = _fstrLL0(self, prec, False)
309 return t if sep is None else sep.join(t)
312class AzimuthalError(_ValueError):
313 '''An azimuthal L{Equidistant}, L{EquidistantKarney}, L{Gnomonic},
314 L{LambertEqualArea}, L{Orthographic}, L{Stereographic} or
315 L{Azimuthal7Tuple} issue.
316 '''
317 pass
320class Azimuthal7Tuple(_NamedTuple):
321 '''7-Tuple C{(x, y, lat, lon, azimuth, scale, datum)}, in C{meter}, C{meter},
322 C{degrees90}, C{degrees180}, compass C{degrees}, C{scalar} and C{Datum}
323 where C{(x, y)} is the easting and northing of a projected point, C{(lat,
324 lon)} the geodetic location, C{azimuth} the azimuth, clockwise from true
325 North and C{scale} is the projection scale, either C{1 / reciprocal} or
326 C{1} or C{-1} in the L{Equidistant} case.
327 '''
328 _Names_ = (_x_, _y_, _lat_, _lon_, _azimuth_, _scale_, _datum_)
329 _Units_ = ( Easting, Northing, Lat_, Lon_, Bearing, Scalar, _Pass)
331 def antipodal(self, azimuth=None):
332 '''Return this tuple with the antipodal C{lat} and C{lon}.
334 @kwarg azimuth: Optional azimuth, overriding the current azimuth
335 (C{compass degrees360}).
336 '''
337 a = antipode(self.lat, self.lon) # PYCHOK named
338 z = self.azimuth if azimuth is None else Bearing(azimuth=azimuth) # PYCHOK named
339 return _NamedTuple.dup(self, lat=a.lat, lon=a.lon, azimuth=z)
342class Equidistant(_AzimuthalBase):
343 '''Azimuthal equidistant projection for the sphere***, see U{Snyder, pp 195-197
344 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
345 <https://MathWorld.Wolfram.com/AzimuthalEquidistantProjection.html>}.
347 @note: Results from this L{Equidistant} and an L{EquidistantExact},
348 L{EquidistantGeodSolve} or L{EquidistantKarney} projection
349 C{may differ} by 10% or more. For an example, see method
350 C{testDiscrepancies} in module C{testAzimuthal.py}.
351 '''
352 if _FOR_DOCS:
353 __init__ = _AzimuthalBase.__init__
355 def forward(self, lat, lon, name=NN):
356 '''Convert a geodetic location to azimuthal equidistant east- and northing.
358 @arg lat: Latitude of the location (C{degrees90}).
359 @arg lon: Longitude of the location (C{degrees180}).
360 @kwarg name: Optional name for the location (C{str}).
362 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
363 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
364 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
365 The C{scale} of the projection is C{1} in I{radial} direction and
366 is C{1 / reciprocal} in the direction perpendicular to this.
368 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
370 @note: The C{scale} will be C{-1} if B{C{(lat, lon)}} is antipodal to
371 the projection center C{(lat0, lon0)}.
372 '''
373 def _k_t(c):
374 k = _N_1_0 if c < 0 else _1_0
375 t = fabs(c) < EPS1
376 if t:
377 a = acos(c)
378 s = sin(a)
379 if s:
380 k = a / s
381 return k, t
383 return self._forward(lat, lon, name, _k_t)
385 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
386 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
388 @arg x: Easting of the location (C{meter}).
389 @arg y: Northing of the location (C{meter}).
390 @kwarg name: Optional name for the location (C{str}).
391 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
392 additional B{C{LatLon}} keyword arguments,
393 ignored if C{B{LatLon} is None} or not given.
395 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
396 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
398 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
399 in the range C{[-180..180] degrees}. The C{scale} of the
400 projection is C{1} in I{radial} direction, C{azimuth} clockwise
401 from true North and is C{1 / reciprocal} in the direction
402 perpendicular to this.
403 '''
404 def _c(c):
405 return c if c > EPS else None
407 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
410def equidistant(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
411 '''Return an L{EquidistantExact}, L{EquidistantGeodSolve} or (if I{Karney}'s
412 U{geographiclib<https://PyPI.org/project/geographiclib>} package is
413 installed) an L{EquidistantKarney}, otherwise an L{Equidistant} instance.
415 @arg lat0: Latitude of center point (C{degrees90}).
416 @arg lon0: Longitude of center point (C{degrees180}).
417 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
418 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
419 radius (C{meter}).
420 @kwarg exact: Return an L{EquidistantExact} instance.
421 @kwarg geodsolve: Return an L{EquidistantGeodSolve} instance.
422 @kwarg name: Optional name for the projection (C{str}).
424 @return: An L{EquidistantExact}, L{EquidistantGeodSolve},
425 L{EquidistantKarney} or L{Equidistant} instance.
427 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
429 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
430 or I{Karney}'s wrapped C{Geodesic}.
432 @raise TypeError: Invalid B{C{datum}}.
433 '''
435 E = EquidistantExact if exact else (EquidistantGeodSolve if geodsolve else Equidistant)
436 if E is Equidistant:
437 try:
438 return EquidistantKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
439 except ImportError:
440 pass
441 return E(lat0, lon0, datum=datum, name=name) # PYCHOK types
444class _AzimuthalGeodesic(_AzimuthalBase):
445 '''(INTERNAL) Base class for azimuthal projections using the
446 I{wrapped} U{geodesic.Geodesic and geodesicline.GeodesicLine
447 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} or the
448 I{exact} geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
449 '''
450 _mask = 0
452 @Property_RO
453 def geodesic(self): # PYCHOK no cover
454 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
455 '''
456 notOverloaded(self)
458 def _7Tuple(self, e, n, r, M=None, name=NN):
459 '''(INTERNAL) Return an C{Azimuthal7Tuple}.
460 '''
461 s = M if M is not None else ( # reciprocal, azimuthal scale
462 (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0)
463 z = _umod_360(r.azi2) # -180 <= r.azi2 < 180 ... 0 <= z < 360
464 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum,
465 name=name or self.name)
468class _EquidistantBase(_AzimuthalGeodesic):
469 '''(INTERNAL) Base for classes L{EquidistantExact}, L{EquidistantGeodSolve}
470 and L{EquidistantKarney}.
471 '''
472 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
473 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or
474 L{EquidistantKarney} projection.
475 '''
476 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
478 g = self.geodesic
479 # g.STANDARD = g.AZIMUTH | g.DISTANCE | g.LATITUDE | g.LONGITUDE
480 self._mask = g.REDUCEDLENGTH | g.STANDARD # | g.LONG_UNROLL
482 def forward(self, lat, lon, name=NN):
483 '''Convert an (ellipsoidal) geodetic location to azimuthal equidistant east- and northing.
485 @arg lat: Latitude of the location (C{degrees90}).
486 @arg lon: Longitude of the location (C{degrees180}).
487 @kwarg name: Optional name for the location (C{str}).
489 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
490 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
491 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
492 The C{scale} of the projection is C{1} in I{radial} direction and
493 is C{1 / reciprocal} in the direction perpendicular to this.
495 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
497 @note: A call to C{.forward} followed by a call to C{.reverse} will return
498 the original C{lat, lon} to within roundoff.
499 '''
500 r = self.geodesic.Inverse(self.lat0, self.lon0,
501 Lat_(lat), Lon_(lon), outmask=self._mask)
502 x, y = sincos2d(r.azi1)
503 return self._7Tuple(x * r.s12, y * r.s12, r, name=name)
505 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
506 '''Convert an azimuthal equidistant location to (ellipsoidal) geodetic lat- and longitude.
508 @arg x: Easting of the location (C{meter}).
509 @arg y: Northing of the location (C{meter}).
510 @kwarg name: Optional name for the location (C{str}).
511 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
512 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
513 arguments, ignored if C{B{LatLon} is None}.
515 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
516 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
518 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
519 in the range C{[-180..180] degrees}. The scale of the projection
520 is C{1} in I{radial} direction, C{azimuth} clockwise from true
521 North and is C{1 / reciprocal} in the direction perpendicular
522 to this.
523 '''
524 e, n, z, s = _enzh4(x, y)
526 r = self.geodesic.Direct(self.lat0, self.lon0, z, s, outmask=self._mask)
527 return self._7Tuple(e, n, r, name=name) if LatLon is None else \
528 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
531class EquidistantExact(_EquidistantBase):
532 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
533 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
534 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
536 An azimuthal equidistant projection is centered at an arbitrary position on the ellipsoid.
537 For a point in projected space C{(x, y)}, the geodesic distance from the center position
538 is C{hypot(x, y)} and the C{azimuth} of the geodesic from the center point is C{atan2(x, y)},
539 clockwise from true North.
541 The C{.forward} and C{.reverse} methods also return the C{azimuth} of the geodesic at C{(x,
542 y)} and the C{scale} in the azimuthal direction which, together with the basic properties
543 of the projection, serve to specify completely the local affine transformation between
544 geographic and projected coordinates.
545 '''
546 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
547 '''New azimuthal L{EquidistantExact} projection.
549 @arg lat0: Latitude of center point (C{degrees90}).
550 @arg lon0: Longitude of center point (C{degrees180}).
551 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
552 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
553 radius (C{meter}).
554 @kwarg name: Optional name for the projection (C{str}).
556 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
557 '''
558 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
560 if _FOR_DOCS:
561 forward = _EquidistantBase.forward
562 reverse = _EquidistantBase.reverse
564 @Property_RO
565 def geodesic(self):
566 '''Get this projection's exact geodesic (L{GeodesicExact}).
567 '''
568 return self.datum.ellipsoid.geodesicx
571class EquidistantGeodSolve(_EquidistantBase):
572 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
573 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
574 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and intended
575 I{for testing purposes only}.
577 @see: L{EquidistantExact} and module L{geodsolve}.
578 '''
579 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
580 '''New azimuthal L{EquidistantGeodSolve} projection.
582 @arg lat0: Latitude of center point (C{degrees90}).
583 @arg lon0: Longitude of center point (C{degrees180}).
584 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
585 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
586 radius (C{meter}).
587 @kwarg name: Optional name for the projection (C{str}).
589 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
590 '''
591 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
593 if _FOR_DOCS:
594 forward = _EquidistantBase.forward
595 reverse = _EquidistantBase.reverse
597 @Property_RO
598 def geodesic(self):
599 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
600 '''
601 return self.datum.ellipsoid.geodsolve
604class EquidistantKarney(_EquidistantBase):
605 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
606 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
607 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
609 @see: L{EquidistantExact}.
610 '''
611 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
612 '''New azimuthal L{EquidistantKarney} projection.
614 @arg lat0: Latitude of center point (C{degrees90}).
615 @arg lon0: Longitude of center point (C{degrees180}).
616 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
617 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
618 radius (C{meter}).
619 @kwarg name: Optional name for the projection (C{str}).
621 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
622 not installed or not found.
624 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
625 '''
626 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
628 if _FOR_DOCS:
629 forward = _EquidistantBase.forward
630 reverse = _EquidistantBase.reverse
632 @Property_RO
633 def geodesic(self):
634 '''Get this projection's I{wrapped} U{geodesic.Geodesic
635 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
636 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
637 package is installed.
638 '''
639 return self.datum.ellipsoid.geodesic
642_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,
643 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI
646class Gnomonic(_AzimuthalBase):
647 '''Azimuthal gnomonic projection for the sphere***, see U{Snyder, pp 164-168
648 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
649 <https://MathWorld.Wolfram.com/GnomonicProjection.html>}.
650 '''
651 if _FOR_DOCS:
652 __init__ = _AzimuthalBase.__init__
654 def forward(self, lat, lon, name=NN):
655 '''Convert a geodetic location to azimuthal equidistant east- and northing.
657 @arg lat: Latitude of the location (C{degrees90}).
658 @arg lon: Longitude of the location (C{degrees180}).
659 @kwarg name: Optional name for the location (C{str}).
661 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
662 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
663 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
664 The C{scale} of the projection is C{1} in I{radial} direction and
665 is C{1 / reciprocal} in the direction perpendicular to this.
667 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
668 '''
669 def _k_t(c):
670 t = c > EPS
671 k = (_1_0 / c) if t else _1_0
672 return k, t
674 return self._forward(lat, lon, name, _k_t)
676 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
677 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
679 @arg x: Easting of the location (C{meter}).
680 @arg y: Northing of the location (C{meter}).
681 @kwarg name: Optional name for the location (C{str}).
682 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
683 additional B{C{LatLon}} keyword arguments,
684 ignored if C{B{LatLon} is None} or not given.
686 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
687 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
689 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
690 in the range C{[-180..180] degrees}. The C{scale} of the
691 projection is C{1} in I{radial} direction, C{azimuth} clockwise
692 from true North and C{1 / reciprocal} in the direction
693 perpendicular to this.
694 '''
695 def _c(c):
696 return atan(c) if c > EPS else None
698 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
701def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
702 '''Return a L{GnomonicExact} or (if I{Karney}'s U{geographiclib
703 <https://PyPI.org/project/geographiclib>} package is installed)
704 a L{GnomonicKarney}, otherwise a L{Gnomonic} instance.
706 @arg lat0: Latitude of center point (C{degrees90}).
707 @arg lon0: Longitude of center point (C{degrees180}).
708 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
709 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
710 radius (C{meter}).
711 @kwarg exact: Return a L{GnomonicExact} instance.
712 @kwarg geodsolve: Return a L{GnomonicGeodSolve} instance.
713 @kwarg name: Optional name for the projection (C{str}).
715 @return: A L{GnomonicExact}, L{GnomonicGeodSolve},
716 L{GnomonicKarney} or L{Gnomonic} instance.
718 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or
719 (spherical) B{C{datum}}.
721 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
722 or I{Karney}'s wrapped C{Geodesic}.
724 @raise TypeError: Invalid B{C{datum}}.
725 '''
726 G = GnomonicExact if exact else (GnomonicGeodSolve if geodsolve else Gnomonic)
727 if G is Gnomonic:
728 try:
729 return GnomonicKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
730 except ImportError:
731 pass
732 return G(lat0, lon0, datum=datum, name=name) # PYCHOK types
735class _GnomonicBase(_AzimuthalGeodesic):
736 '''(INTERNAL) Base for classes L{GnomonicExact}, L{GnomonicGeodSolve}
737 and L{GnomonicKarney}.
738 '''
739 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
740 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection.
741 '''
742 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
744 g = self.geodesic
745 self._mask = g.ALL # | g.LONG_UNROLL
747 def forward(self, lat, lon, name=NN, raiser=True): # PYCHOK signature
748 '''Convert an (ellipsoidal) geodetic location to azimuthal gnomonic east-
749 and northing.
751 @arg lat: Latitude of the location (C{degrees90}).
752 @arg lon: Longitude of the location (C{degrees180}).
753 @kwarg name: Optional name for the location (C{str}).
754 @kwarg raiser: Do or don't throw an error (C{bool}) if
755 the location lies over the horizon.
757 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
758 with easting C{x} and northing C{y} in C{meter} and C{lat} and
759 C{lon} in C{degrees} and C{azimuth} clockwise from true North.
760 The C{scale} of the projection is C{1 / reciprocal**2} in I{radial}
761 direction and C{1 / reciprocal} in the direction perpendicular to
762 this. Both C{x} and C{y} will be C{NAN} if the (geodetic) location
763 lies over the horizon and C{B{raiser}=False}.
765 @raise AzimuthalError: Invalid B{C{lat}}, B{C{lon}} or the location lies
766 over the horizon and C{B{raiser}=True}.
767 '''
768 self._iteration = 0
770 r = self.geodesic.Inverse(self.lat0, self.lon0,
771 Lat_(lat), Lon_(lon), outmask=self._mask)
772 M = r.M21
773 if M > EPS0:
774 q = r.m12 / M # .M12
775 e, n = sincos2d(r.azi1)
776 e *= q
777 n *= q
778 elif raiser: # PYCHOK no cover
779 raise AzimuthalError(lat=lat, lon=lon, txt=_over_horizon_)
780 else: # PYCHOK no cover
781 e = n = NAN
783 t = self._7Tuple(e, n, r, M=M, name=name)
784 t._iteraton = 0
785 return t
787 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
788 '''Convert an azimuthal gnomonic location to (ellipsoidal) geodetic lat- and longitude.
790 @arg x: Easting of the location (C{meter}).
791 @arg y: Northing of the location (C{meter}).
792 @kwarg name: Optional name for the location (C{str}).
793 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
794 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
795 arguments, ignored if C{B{LatLon} is None}.
797 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
798 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
800 @raise AzimuthalError: No convergence.
802 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
803 in the range C{[-180..180] degrees}. The C{azimuth} is clockwise
804 from true North. The scale is C{1 / reciprocal**2} in C{radial}
805 direction and C{1 / reciprocal} in the direction perpendicular
806 to this.
807 '''
808 e, n, z, q = _enzh4(x, y)
810 d = a = self.equatoradius
811 s = a * atan(q / a)
812 if q > a: # PYCHOK no cover
813 def _d(r, q):
814 return (r.M12 - q * r.m12) * r.m12 # negated
816 q = _1_0 / q
817 else: # little == True
818 def _d(r, q): # PYCHOK _d
819 return (q * r.M12 - r.m12) * r.M12 # negated
821 a *= _EPS_K
822 m = self._mask
823 g = self.geodesic
825 _P = g.Line(self.lat0, self.lon0, z, m | g.LINE_OFF).Position
826 _S2 = Fsum(s).fsum2_
827 for i in range(1, _TRIPS):
828 r = _P(s, outmask=m)
829 if fabs(d) < a:
830 break
831 s, d = _S2(_d(r, q))
832 else: # PYCHOK no cover
833 self._iteration = _TRIPS
834 raise AzimuthalError(Fmt.no_convergence(d, a),
835 txt=unstr(self.reverse, x, y))
837 t = self._7Tuple(e, n, r, M=r.M12, name=name) if LatLon is None else \
838 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
839 t._iteration = self._iteration = i
840 return t
843class GnomonicExact(_GnomonicBase):
844 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
845 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
846 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
848 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/
849 classGeographicLib_1_1Gnomonic.html>}, especially the B{Warning}.
850 '''
851 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
852 '''New azimuthal L{GnomonicExact} projection.
854 @arg lat0: Latitude of center point (C{degrees90}).
855 @arg lon0: Longitude of center point (C{degrees180}).
856 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
857 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
858 radius (C{meter}).
859 @kwarg name: Optional name for the projection (C{str}).
861 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
862 '''
863 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
865 if _FOR_DOCS:
866 forward = _GnomonicBase.forward
867 reverse = _GnomonicBase.reverse
869 @Property_RO
870 def geodesic(self):
871 '''Get this projection's exact geodesic (L{GeodesicExact}).
872 '''
873 return self.datum.ellipsoid.geodesicx
876class GnomonicGeodSolve(_GnomonicBase):
877 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
878 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
879 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and
880 intended I{for testing purposes only}.
882 @see: L{GnomonicExact} and module L{geodsolve}.
883 '''
884 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
885 '''New azimuthal L{GnomonicGeodSolve} projection.
887 @arg lat0: Latitude of center point (C{degrees90}).
888 @arg lon0: Longitude of center point (C{degrees180}).
889 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
890 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
891 radius (C{meter}).
892 @kwarg name: Optional name for the projection (C{str}).
894 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
895 '''
896 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
898 if _FOR_DOCS:
899 forward = _GnomonicBase.forward
900 reverse = _GnomonicBase.reverse
902 @Property_RO
903 def geodesic(self):
904 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
905 '''
906 return self.datum.ellipsoid.geodsolve
909class GnomonicKarney(_GnomonicBase):
910 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
911 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
912 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
914 @see: L{GnomonicExact}.
915 '''
916 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
917 '''New azimuthal L{GnomonicKarney} projection.
919 @arg lat0: Latitude of center point (C{degrees90}).
920 @arg lon0: Longitude of center point (C{degrees180}).
921 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
922 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
923 radius (C{meter}).
924 @kwarg name: Optional name for the projection (C{str}).
926 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
927 not installed or not found.
929 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
930 '''
931 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
933 if _FOR_DOCS:
934 forward = _GnomonicBase.forward
935 reverse = _GnomonicBase.reverse
937 @Property_RO
938 def geodesic(self):
939 '''Get this projection's I{wrapped} U{geodesic.Geodesic
940 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
941 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
942 package is installed.
943 '''
944 return self.datum.ellipsoid.geodesic
947class LambertEqualArea(_AzimuthalBase):
948 '''Lambert-equal-area projection for the sphere*** (aka U{Lambert zenithal equal-area
949 projection<https://WikiPedia.org/wiki/Lambert_azimuthal_equal-area_projection>}, see
950 U{Snyder, pp 185-187<https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
951 <https://MathWorld.Wolfram.com/LambertAzimuthalEqual-AreaProjection.html>}.
952 '''
953 if _FOR_DOCS:
954 __init__ = _AzimuthalBase.__init__
956 def forward(self, lat, lon, name=NN):
957 '''Convert a geodetic location to azimuthal Lambert-equal-area east- and northing.
959 @arg lat: Latitude of the location (C{degrees90}).
960 @arg lon: Longitude of the location (C{degrees180}).
961 @kwarg name: Optional name for the location (C{str}).
963 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
964 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
965 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
966 The C{scale} of the projection is C{1} in I{radial} direction and
967 is C{1 / reciprocal} in the direction perpendicular to this.
969 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
970 '''
971 def _k_t(c):
972 c += _1_0
973 t = c > EPS0
974 k = sqrt(_2_0 / c) if t else _1_0
975 return k, t
977 return self._forward(lat, lon, name, _k_t)
979 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
980 '''Convert an azimuthal Lambert-equal-area location to geodetic lat- and longitude.
982 @arg x: Easting of the location (C{meter}).
983 @arg y: Northing of the location (C{meter}).
984 @kwarg name: Optional name for the location (C{str}).
985 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
986 additional B{C{LatLon}} keyword arguments,
987 ignored if C{B{LatLon} is None} or not given.
989 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
990 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
992 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
993 in the range C{[-180..180] degrees}. The C{scale} of the
994 projection is C{1} in I{radial} direction, C{azimuth} clockwise
995 from true North and is C{1 / reciprocal} in the direction
996 perpendicular to this.
997 '''
998 def _c(c):
999 c *= _0_5
1000 return (asin1(c) * _2_0) if c > EPS else None
1002 return self._reverse(x, y, name, _c, True, **LatLon_and_kwds)
1005class Orthographic(_AzimuthalBase):
1006 '''Orthographic projection for the sphere***, see U{Snyder, pp 148-153
1007 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1008 <https://MathWorld.Wolfram.com/OrthographicProjection.html>}.
1009 '''
1010 if _FOR_DOCS:
1011 __init__ = _AzimuthalBase.__init__
1013 def forward(self, lat, lon, name=NN):
1014 '''Convert a geodetic location to azimuthal orthographic east- and northing.
1016 @arg lat: Latitude of the location (C{degrees90}).
1017 @arg lon: Longitude of the location (C{degrees180}).
1018 @kwarg name: Optional name for the location (C{str}).
1020 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1021 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1022 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1023 The C{scale} of the projection is C{1} in I{radial} direction and
1024 is C{1 / reciprocal} in the direction perpendicular to this.
1026 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1027 '''
1028 def _k_t(c):
1029 return _1_0, (c >= 0)
1031 return self._forward(lat, lon, name, _k_t)
1033 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1034 '''Convert an azimuthal orthographic location to geodetic lat- and longitude.
1036 @arg x: Easting of the location (C{meter}).
1037 @arg y: Northing of the location (C{meter}).
1038 @kwarg name: Optional name for the location (C{str}).
1039 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1040 additional B{C{LatLon}} keyword arguments,
1041 ignored if C{B{LatLon} is None} or not given.
1043 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1044 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1046 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1047 in the range C{[-180..180] degrees}. The C{scale} of the
1048 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1049 from true North and is C{1 / reciprocal} in the direction
1050 perpendicular to this.
1051 '''
1052 def _c(c):
1053 return asin1(c) if c > EPS else None
1055 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1058class Stereographic(_AzimuthalBase):
1059 '''Stereographic projection for the sphere***, see U{Snyder, pp 157-160
1060 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1061 <https://MathWorld.Wolfram.com/StereographicProjection.html>}.
1062 '''
1063 _k0 = _1_0 # central scale factor (C{scalar})
1064 _k02 = _2_0 # double ._k0
1066 if _FOR_DOCS:
1067 __init__ = _AzimuthalBase.__init__
1069 def forward(self, lat, lon, name=NN):
1070 '''Convert a geodetic location to azimuthal stereographic east- and northing.
1072 @arg lat: Latitude of the location (C{degrees90}).
1073 @arg lon: Longitude of the location (C{degrees180}).
1074 @kwarg name: Optional name for the location (C{str}).
1076 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1077 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1078 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1079 The C{scale} of the projection is C{1} in I{radial} direction and
1080 is C{1 / reciprocal} in the direction perpendicular to this.
1082 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1083 '''
1084 def _k_t(c):
1085 c += _1_0
1086 t = isnon0(c)
1087 k = (self._k02 / c) if t else _1_0
1088 return k, t
1090 return self._forward(lat, lon, name, _k_t)
1092 @property_doc_(''' optional, central scale factor (C{scalar}).''')
1093 def k0(self):
1094 '''Get the central scale factor (C{scalar}).
1095 '''
1096 return self._k0
1098 @k0.setter # PYCHOK setter!
1099 def k0(self, factor):
1100 '''Set the central scale factor (C{scalar}).
1101 '''
1102 n = Stereographic.k0.fget.__name__
1103 self._k0 = Scalar_(factor, name=n, low=EPS, high=2) # XXX high=1, 2, other?
1104 self._k02 = self._k0 * _2_0
1106 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1107 '''Convert an azimuthal stereographic location to geodetic lat- and longitude.
1109 @arg x: Easting of the location (C{meter}).
1110 @arg y: Northing of the location (C{meter}).
1111 @kwarg name: Optional name for the location (C{str}).
1112 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1113 additional B{C{LatLon}} keyword arguments,
1114 ignored if C{B{LatLon} is None} or not given.
1116 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1117 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1119 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1120 in the range C{[-180..180] degrees}. The C{scale} of the
1121 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1122 from true North and is C{1 / reciprocal} in the direction
1123 perpendicular to this.
1124 '''
1125 def _c(c):
1126 return (_2_0 * atan2(c, self._k02)) if c > EPS else None
1128 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1131__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase)
1133# **) MIT License
1134#
1135# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1136#
1137# Permission is hereby granted, free of charge, to any person obtaining a
1138# copy of this software and associated documentation files (the "Software"),
1139# to deal in the Software without restriction, including without limitation
1140# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1141# and/or sell copies of the Software, and to permit persons to whom the
1142# Software is furnished to do so, subject to the following conditions:
1143#
1144# The above copyright notice and this permission notice shall be included
1145# in all copies or substantial portions of the Software.
1146#
1147# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1148# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1149# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1150# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1151# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1152# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1153# OTHER DEALINGS IN THE SOFTWARE.