Coverage for pygeodesy/azimuthal.py: 98%
312 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-02 08:40 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-02 08:40 -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.03.21'
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 @Property_RO
137 def flattening(self):
138 '''Get the geodesic's flattening (C{float}).
139 '''
140 return self.datum.ellipsoid.f
142 def forward(self, lat, lon, name=NN): # PYCHOK no cover
143 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
144 '''
145 notOverloaded(self, lat, lon, name=name)
147 def _forward(self, lat, lon, name, _k_t_2):
148 '''(INTERNAL) Azimuthal (spherical) forward C{lat, lon} to C{x, y}.
149 '''
150 lat, lon = Lat_(lat), Lon_(lon)
151 sa, ca, sb, cb = sincos2d_(lat, lon - self.lon0)
152 s0, c0 = self._sc0
154 cb *= ca
155 k, t = _k_t_2(s0 * sa + c0 * cb)
156 if t:
157 r = k * self.radius
158 y = r * (c0 * sa - s0 * cb)
159 e, n, z, _ = _enzh4(r * sb * ca, y, None)
160 else: # 0 or 180
161 e = n = z = _0_0
163 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
164 name=name or self.name)
165 return t
167 def _forwards(self, *lls):
168 '''(INTERNAL) One or more C{.forward} calls, see .ellipsoidalBaseDI.
169 '''
170 _fwd = self.forward
171 for ll in lls:
172 yield _fwd(ll.lat, ll.lon)
174 @Property_RO
175 def lat0(self):
176 '''Get the center latitude (C{degrees90}).
177 '''
178 return self._latlon0.lat
180 @property
181 def latlon0(self):
182 '''Get the center lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}) in (C{degrees90}, C{degrees180}).
183 '''
184 return self._latlon0
186 @latlon0.setter # PYCHOK setter!
187 def latlon0(self, latlon0):
188 '''Set the center lat- and longitude (C{LatLon}, L{LatLon2Tuple} or L{LatLon4Tuple}).
190 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or ellipsoidal mismatch
191 of B{C{latlon0}} and this projection.
192 '''
193 B = _LLEB if self.datum.isEllipsoidal else _LLB
194 _xinstanceof(B, LatLon2Tuple, LatLon4Tuple, latlon0=latlon0)
195 if hasattr(latlon0, _datum_):
196 _xdatum(self.datum, latlon0.datum, Error=AzimuthalError)
197 self.reset(latlon0.lat, latlon0.lon)
199 @Property_RO
200 def lon0(self):
201 '''Get the center longitude (C{degrees180}).
202 '''
203 return self._latlon0.lon
205 @deprecated_Property_RO
206 def majoradius(self): # PYCHOK no cover
207 '''DEPRECATED, use property C{equatoradius}.'''
208 return self.equatoradius
210 @Property_RO
211 def radius(self):
212 '''Get this projection's mean radius of curvature (C{meter}).
213 '''
214 return self.datum.ellipsoid.rocMean(self.lat0)
216 def reset(self, lat0, lon0):
217 '''Set or reset the center point of this azimuthal projection.
219 @arg lat0: Center point latitude (C{degrees90}).
220 @arg lon0: Center point longitude (C{degrees180}).
222 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
223 '''
224 _update_all(self) # zap caches
225 self._reset(lat0, lon0)
227 def _reset(self, lat0, lon0):
228 '''(INTERNAL) Update the center point.
229 '''
230 self._latlon0 = LatLon2Tuple(Lat_(lat0=lat0, Error=AzimuthalError),
231 Lon_(lon0=lon0, Error=AzimuthalError))
232 self._sc0 = sincos2d(self.lat0)
234 def reverse(self, x, y, name=NN, **LatLon_and_kwds): # PYCHOK no cover
235 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
236 '''
237 notOverloaded(self, x, y, name=name, **LatLon_and_kwds)
239 def _reverse(self, x, y, name, _c, lea, LatLon=None, **LatLon_kwds):
240 '''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
241 '''
242 e, n, z, r = _enzh4(x, y)
244 c = _c(r / self.radius)
245 if c is None:
246 lat, lon = self.latlon0
247 k, z = _1_0, _0_0
248 else:
249 s0, c0 = self._sc0
250 sc, cc = sincos2(c)
251 k = c / sc
252 s = s0 * cc
253 if r > EPS0:
254 s += c0 * sc * (n / r)
255 lat = degrees(asin1(s))
256 if lea or fabs(c0) > EPS:
257 d = atan2d(e * sc, r * c0 * cc - n * s0 * sc)
258 else:
259 d = atan2d(e, (n if s0 < 0 else -n))
260 lon = _norm180(self.lon0 + d)
262 if LatLon is None:
263 t = Azimuthal7Tuple(e, n, lat, lon, z, k, self.datum,
264 name=name or self.name)
265 else:
266 t = self._toLatLon(lat, lon, LatLon, LatLon_kwds, name)
267 return t
269 def _reverse2(self, x_t, *y):
270 '''(INTERNAL) See iterating functions .ellipsoidalBaseDI._intersect3,
271 .ellipsoidalBaseDI._intersects2 and .ellipsoidalBaseDI._nearestOne.
272 '''
273 t = self.reverse(x_t, *y) if y else self.reverse(x_t.x, x_t.y) # LatLon=None
274 d = euclid(t.lat - self.lat0, t.lon - self.lon0) # degrees
275 return t, d
277 def _toLatLon(self, lat, lon, LatLon, LatLon_kwds, name):
278 '''(INTERNAL) Check B{C{LatLon}} and return an instance.
279 '''
280 kwds = _xkwds(LatLon_kwds, datum=self.datum)
281 r = self._xnamed(LatLon(lat, lon, **kwds), name=name) # handle .classof
282 B = _LLEB if self.datum.isEllipsoidal else _LLB
283 _xinstanceof(B, LatLon=r)
284 return r
286 def toRepr(self, prec=6, **unused): # PYCHOK expected
287 '''Return a string representation of this projection.
289 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
291 @return: This projection as C{"<classname>(lat0, lon0, ...)"}
292 (C{str}).
293 '''
294 return _fstrLL0(self, prec, True)
296 def toStr(self, prec=6, sep=_SPACE_, **unused): # PYCHOK expected
297 '''Return a string representation of this projection.
299 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
300 @kwarg sep: Separator to join (C{str}).
302 @return: This projection as C{"lat0 lon0"} (C{str}).
303 '''
304 t = _fstrLL0(self, prec, False)
305 return t if sep is None else sep.join(t)
308class AzimuthalError(_ValueError):
309 '''An azimuthal L{Equidistant}, L{EquidistantKarney}, L{Gnomonic},
310 L{LambertEqualArea}, L{Orthographic}, L{Stereographic} or
311 L{Azimuthal7Tuple} issue.
312 '''
313 pass
316class Azimuthal7Tuple(_NamedTuple):
317 '''7-Tuple C{(x, y, lat, lon, azimuth, scale, datum)}, in C{meter}, C{meter},
318 C{degrees90}, C{degrees180}, compass C{degrees}, C{scalar} and C{Datum}
319 where C{(x, y)} is the easting and northing of a projected point, C{(lat,
320 lon)} the geodetic location, C{azimuth} the azimuth, clockwise from true
321 North and C{scale} is the projection scale, either C{1 / reciprocal} or
322 C{1} or C{-1} in the L{Equidistant} case.
323 '''
324 _Names_ = (_x_, _y_, _lat_, _lon_, _azimuth_, _scale_, _datum_)
325 _Units_ = ( Easting, Northing, Lat_, Lon_, Bearing, Scalar, _Pass)
327 def antipodal(self, azimuth=None):
328 '''Return this tuple with the antipodal C{lat} and C{lon}.
330 @kwarg azimuth: Optional azimuth, overriding the current azimuth
331 (C{compass degrees360}).
332 '''
333 a = antipode(self.lat, self.lon) # PYCHOK named
334 z = self.azimuth if azimuth is None else Bearing(azimuth=azimuth) # PYCHOK named
335 return _NamedTuple.dup(self, lat=a.lat, lon=a.lon, azimuth=z)
338class Equidistant(_AzimuthalBase):
339 '''Azimuthal equidistant projection for the sphere***, see U{Snyder, pp 195-197
340 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
341 <https://MathWorld.Wolfram.com/AzimuthalEquidistantProjection.html>}.
343 @note: Results from this L{Equidistant} and an L{EquidistantExact},
344 L{EquidistantGeodSolve} or L{EquidistantKarney} projection
345 C{may differ} by 10% or more. For an example, see method
346 C{testDiscrepancies} in module C{testAzimuthal.py}.
347 '''
348 if _FOR_DOCS:
349 __init__ = _AzimuthalBase.__init__
351 def forward(self, lat, lon, name=NN):
352 '''Convert a geodetic location to azimuthal equidistant east- and northing.
354 @arg lat: Latitude of the location (C{degrees90}).
355 @arg lon: Longitude of the location (C{degrees180}).
356 @kwarg name: Optional name for the location (C{str}).
358 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
359 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
360 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
361 The C{scale} of the projection is C{1} in I{radial} direction and
362 is C{1 / reciprocal} in the direction perpendicular to this.
364 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
366 @note: The C{scale} will be C{-1} if B{C{(lat, lon)}} is antipodal to
367 the projection center C{(lat0, lon0)}.
368 '''
369 def _k_t(c):
370 k = _N_1_0 if c < 0 else _1_0
371 t = fabs(c) < EPS1
372 if t:
373 a = acos(c)
374 s = sin(a)
375 if s:
376 k = a / s
377 return k, t
379 return self._forward(lat, lon, name, _k_t)
381 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
382 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
384 @arg x: Easting of the location (C{meter}).
385 @arg y: Northing of the location (C{meter}).
386 @kwarg name: Optional name for the location (C{str}).
387 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
388 additional B{C{LatLon}} keyword arguments,
389 ignored if C{B{LatLon} is None} or not given.
391 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
392 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
394 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
395 in the range C{[-180..180] degrees}. The C{scale} of the
396 projection is C{1} in I{radial} direction, C{azimuth} clockwise
397 from true North and is C{1 / reciprocal} in the direction
398 perpendicular to this.
399 '''
400 def _c(c):
401 return c if c > EPS else None
403 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
406def equidistant(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
407 '''Return an L{EquidistantExact}, L{EquidistantGeodSolve} or (if I{Karney}'s
408 U{geographiclib<https://PyPI.org/project/geographiclib>} package is
409 installed) an L{EquidistantKarney}, otherwise an L{Equidistant} instance.
411 @arg lat0: Latitude of center point (C{degrees90}).
412 @arg lon0: Longitude of center point (C{degrees180}).
413 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
414 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
415 radius (C{meter}).
416 @kwarg exact: Return an L{EquidistantExact} instance.
417 @kwarg geodsolve: Return an L{EquidistantGeodSolve} instance.
418 @kwarg name: Optional name for the projection (C{str}).
420 @return: An L{EquidistantExact}, L{EquidistantGeodSolve},
421 L{EquidistantKarney} or L{Equidistant} instance.
423 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or (spherical) B{C{datum}}.
425 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
426 or I{Karney}'s wrapped C{Geodesic}.
428 @raise TypeError: Invalid B{C{datum}}.
429 '''
431 E = EquidistantExact if exact else (EquidistantGeodSolve if geodsolve else Equidistant)
432 if E is Equidistant:
433 try:
434 return EquidistantKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
435 except ImportError:
436 pass
437 return E(lat0, lon0, datum=datum, name=name) # PYCHOK types
440class _AzimuthalGeodesic(_AzimuthalBase):
441 '''(INTERNAL) Base class for azimuthal projections using U{Karney Geodesic
442 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>} or
443 exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
444 '''
445 _mask = 0
447 @Property_RO
448 def geodesic(self): # PYCHOK no cover
449 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
450 '''
451 notOverloaded(self)
453 def _7Tuple(self, e, n, r, M=None, name=NN):
454 '''(INTERNAL) Return an C{Azimuthal7Tuple}.
455 '''
456 s = M if M is not None else ( # reciprocal, azimuthal scale
457 (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0)
458 z = _umod_360(r.azi2) # -180 <= r.azi2 < 180 ... 0 <= z < 360
459 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum,
460 name=name or self.name)
463class _EquidistantBase(_AzimuthalGeodesic):
464 '''(INTERNAL) Base for classes L{EquidistantExact}, L{EquidistantGeodSolve}
465 and L{EquidistantKarney}.
466 '''
467 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
468 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or
469 L{EquidistantKarney} projection.
470 '''
471 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
473 g = self.geodesic
474 # g.STANDARD = g.AZIMUTH | g.DISTANCE | g.LATITUDE | g.LONGITUDE
475 self._mask = g.REDUCEDLENGTH | g.STANDARD # | g.LONG_UNROLL
477 def forward(self, lat, lon, name=NN):
478 '''Convert an (ellipsoidal) geodetic location to azimuthal equidistant east- and northing.
480 @arg lat: Latitude of the location (C{degrees90}).
481 @arg lon: Longitude of the location (C{degrees180}).
482 @kwarg name: Optional name for the location (C{str}).
484 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
485 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
486 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
487 The C{scale} of the projection is C{1} in I{radial} direction and
488 is C{1 / reciprocal} in the direction perpendicular to this.
490 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
492 @note: A call to C{.forward} followed by a call to C{.reverse} will return
493 the original C{lat, lon} to within roundoff.
494 '''
495 r = self.geodesic.Inverse(self.lat0, self.lon0,
496 Lat_(lat), Lon_(lon), self._mask)
497 x, y = sincos2d(r.azi1)
498 return self._7Tuple(x * r.s12, y * r.s12, r, name=name)
500 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
501 '''Convert an azimuthal equidistant location to (ellipsoidal) geodetic lat- and longitude.
503 @arg x: Easting of the location (C{meter}).
504 @arg y: Northing of the location (C{meter}).
505 @kwarg name: Optional name for the location (C{str}).
506 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
507 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
508 arguments, ignored if C{B{LatLon} is None}.
510 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
511 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
513 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
514 in the range C{[-180..180] degrees}. The scale of the projection
515 is C{1} in I{radial} direction, C{azimuth} clockwise from true
516 North and is C{1 / reciprocal} in the direction perpendicular
517 to this.
518 '''
519 e, n, z, s = _enzh4(x, y)
521 r = self.geodesic.Direct(self.lat0, self.lon0, z, s, self._mask)
522 return self._7Tuple(e, n, r, name=name) if LatLon is None else \
523 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
526class EquidistantExact(_EquidistantBase):
527 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
528 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
529 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
531 An azimuthal equidistant projection is centered at an arbitrary position on the ellipsoid.
532 For a point in projected space C{(x, y)}, the geodesic distance from the center position
533 is C{hypot(x, y)} and the C{azimuth} of the geodesic from the center point is C{atan2(x, y)},
534 clockwise from true North.
536 The C{.forward} and C{.reverse} methods also return the C{azimuth} of the geodesic at C{(x,
537 y)} and the C{scale} in the azimuthal direction which, together with the basic properties
538 of the projection, serve to specify completely the local affine transformation between
539 geographic and projected coordinates.
540 '''
541 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
542 '''New azimuthal L{EquidistantExact} projection.
544 @arg lat0: Latitude of center point (C{degrees90}).
545 @arg lon0: Longitude of center point (C{degrees180}).
546 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
547 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
548 radius (C{meter}).
549 @kwarg name: Optional name for the projection (C{str}).
551 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
552 '''
553 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
555 if _FOR_DOCS:
556 forward = _EquidistantBase.forward
557 reverse = _EquidistantBase.reverse
559 @Property_RO
560 def geodesic(self):
561 '''Get this projection's exact geodesic (L{GeodesicExact}).
562 '''
563 return self.datum.ellipsoid.geodesicx
566class EquidistantGeodSolve(_EquidistantBase):
567 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
568 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
569 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and intended
570 I{for testing purposes only}.
572 @see: L{EquidistantExact} and module L{geodsolve}.
573 '''
574 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
575 '''New azimuthal L{EquidistantGeodSolve} projection.
577 @arg lat0: Latitude of center point (C{degrees90}).
578 @arg lon0: Longitude of center point (C{degrees180}).
579 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
580 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
581 radius (C{meter}).
582 @kwarg name: Optional name for the projection (C{str}).
584 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
585 '''
586 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
588 if _FOR_DOCS:
589 forward = _EquidistantBase.forward
590 reverse = _EquidistantBase.reverse
592 @Property_RO
593 def geodesic(self):
594 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
595 '''
596 return self.datum.ellipsoid.geodsolve
599class EquidistantKarney(_EquidistantBase):
600 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
601 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
602 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
604 @see: L{EquidistantExact}.
605 '''
606 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
607 '''New azimuthal L{EquidistantKarney} projection.
609 @arg lat0: Latitude of center point (C{degrees90}).
610 @arg lon0: Longitude of center point (C{degrees180}).
611 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
612 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
613 radius (C{meter}).
614 @kwarg name: Optional name for the projection (C{str}).
616 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
617 not installed or not found.
619 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
620 '''
621 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
623 if _FOR_DOCS:
624 forward = _EquidistantBase.forward
625 reverse = _EquidistantBase.reverse
627 @Property_RO
628 def geodesic(self):
629 '''Get this projection's I{wrapped} U{Karney Geodesic
630 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>},
631 provided package U{geographiclib
632 <https://PyPI.org/project/geographiclib>} is installed.
633 '''
634 return self.datum.ellipsoid.geodesic
637_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,
638 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI
641class Gnomonic(_AzimuthalBase):
642 '''Azimuthal gnomonic projection for the sphere***, see U{Snyder, pp 164-168
643 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
644 <https://MathWorld.Wolfram.com/GnomonicProjection.html>}.
645 '''
646 if _FOR_DOCS:
647 __init__ = _AzimuthalBase.__init__
649 def forward(self, lat, lon, name=NN):
650 '''Convert a geodetic location to azimuthal equidistant east- and northing.
652 @arg lat: Latitude of the location (C{degrees90}).
653 @arg lon: Longitude of the location (C{degrees180}).
654 @kwarg name: Optional name for the location (C{str}).
656 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
657 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
658 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
659 The C{scale} of the projection is C{1} in I{radial} direction and
660 is C{1 / reciprocal} in the direction perpendicular to this.
662 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
663 '''
664 def _k_t(c):
665 t = c > EPS
666 k = (_1_0 / c) if t else _1_0
667 return k, t
669 return self._forward(lat, lon, name, _k_t)
671 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
672 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
674 @arg x: Easting of the location (C{meter}).
675 @arg y: Northing of the location (C{meter}).
676 @kwarg name: Optional name for the location (C{str}).
677 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
678 additional B{C{LatLon}} keyword arguments,
679 ignored if C{B{LatLon} is None} or not given.
681 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
682 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
684 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
685 in the range C{[-180..180] degrees}. The C{scale} of the
686 projection is C{1} in I{radial} direction, C{azimuth} clockwise
687 from true North and C{1 / reciprocal} in the direction
688 perpendicular to this.
689 '''
690 def _c(c):
691 return atan(c) if c > EPS else None
693 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
696def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
697 '''Return a L{GnomonicExact} or (if I{Karney}'s U{geographiclib
698 <https://PyPI.org/project/geographiclib>} package is installed)
699 a L{GnomonicKarney}, otherwise a L{Gnomonic} instance.
701 @arg lat0: Latitude of center point (C{degrees90}).
702 @arg lon0: Longitude of center point (C{degrees180}).
703 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
704 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
705 radius (C{meter}).
706 @kwarg exact: Return a L{GnomonicExact} instance.
707 @kwarg geodsolve: Return a L{GnomonicGeodSolve} instance.
708 @kwarg name: Optional name for the projection (C{str}).
710 @return: A L{GnomonicExact}, L{GnomonicGeodSolve},
711 L{GnomonicKarney} or L{Gnomonic} instance.
713 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or
714 (spherical) B{C{datum}}.
716 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
717 or I{Karney}'s wrapped C{Geodesic}.
719 @raise TypeError: Invalid B{C{datum}}.
720 '''
721 G = GnomonicExact if exact else (GnomonicGeodSolve if geodsolve else Gnomonic)
722 if G is Gnomonic:
723 try:
724 return GnomonicKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
725 except ImportError:
726 pass
727 return G(lat0, lon0, datum=datum, name=name) # PYCHOK types
730class _GnomonicBase(_AzimuthalGeodesic):
731 '''(INTERNAL) Base for classes L{GnomonicExact}, L{GnomonicGeodSolve}
732 and L{GnomonicKarney}.
733 '''
734 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
735 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection.
736 '''
737 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
739 g = self.geodesic
740 self._mask = g.ALL # | g.LONG_UNROLL
742 def forward(self, lat, lon, name=NN, raiser=True): # PYCHOK signature
743 '''Convert an (ellipsoidal) geodetic location to azimuthal gnomonic east-
744 and northing.
746 @arg lat: Latitude of the location (C{degrees90}).
747 @arg lon: Longitude of the location (C{degrees180}).
748 @kwarg name: Optional name for the location (C{str}).
749 @kwarg raiser: Do or don't throw an error (C{bool}) if
750 the location lies over the horizon.
752 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
753 with easting C{x} and northing C{y} in C{meter} and C{lat} and
754 C{lon} in C{degrees} and C{azimuth} clockwise from true North.
755 The C{scale} of the projection is C{1 / reciprocal**2} in I{radial}
756 direction and C{1 / reciprocal} in the direction perpendicular to
757 this. Both C{x} and C{y} will be C{NAN} if the (geodetic) location
758 lies over the horizon and C{B{raiser}=False}.
760 @raise AzimuthalError: Invalid B{C{lat}}, B{C{lon}} or the location lies
761 over the horizon and C{B{raiser}=True}.
762 '''
763 self._iteration = 0
765 r = self.geodesic.Inverse(self.lat0, self.lon0,
766 Lat_(lat), Lon_(lon), self._mask)
767 M = r.M21
768 if M > EPS0:
769 q = r.m12 / M # .M12
770 e, n = sincos2d(r.azi1)
771 e *= q
772 n *= q
773 elif raiser: # PYCHOK no cover
774 raise AzimuthalError(lat=lat, lon=lon, txt=_over_horizon_)
775 else: # PYCHOK no cover
776 e = n = NAN
778 t = self._7Tuple(e, n, r, M=M, name=name)
779 t._iteraton = 0
780 return t
782 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
783 '''Convert an azimuthal gnomonic location to (ellipsoidal) geodetic lat- and longitude.
785 @arg x: Easting of the location (C{meter}).
786 @arg y: Northing of the location (C{meter}).
787 @kwarg name: Optional name for the location (C{str}).
788 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
789 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
790 arguments, ignored if C{B{LatLon} is None}.
792 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
793 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
795 @raise AzimuthalError: No convergence.
797 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
798 in the range C{[-180..180] degrees}. The C{azimuth} is clockwise
799 from true North. The scale is C{1 / reciprocal**2} in C{radial}
800 direction and C{1 / reciprocal} in the direction perpendicular
801 to this.
802 '''
803 e, n, z, q = _enzh4(x, y)
805 d = a = self.equatoradius
806 s = a * atan(q / a)
807 if q > a: # PYCHOK no cover
808 def _d(r, q):
809 return (r.M12 - q * r.m12) * r.m12 # negated
811 q = _1_0 / q
812 else: # little == True
813 def _d(r, q): # PYCHOK _d
814 return (q * r.M12 - r.m12) * r.M12 # negated
816 a *= _EPS_K
817 m = self._mask
818 g = self.geodesic
820 P = g.Line(self.lat0, self.lon0, z, m | g.LINE_OFF).Position
821 _S2 = Fsum(s).fsum2_
822 for i in range(1, _TRIPS):
823 r = P(s, m)
824 if fabs(d) < a:
825 break
826 s, d = _S2(_d(r, q))
827 else: # PYCHOK no cover
828 self._iteration = _TRIPS
829 raise AzimuthalError(Fmt.no_convergence(d, a),
830 txt=unstr(self.reverse, x, y))
832 t = self._7Tuple(e, n, r, M=r.M12, name=name) if LatLon is None else \
833 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
834 t._iteration = self._iteration = i
835 return t
838class GnomonicExact(_GnomonicBase):
839 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
840 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
841 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
843 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/
844 classGeographicLib_1_1Gnomonic.html>}, especially the B{Warning}.
845 '''
846 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
847 '''New azimuthal L{GnomonicExact} projection.
849 @arg lat0: Latitude of center point (C{degrees90}).
850 @arg lon0: Longitude of center point (C{degrees180}).
851 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
852 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
853 radius (C{meter}).
854 @kwarg name: Optional name for the projection (C{str}).
856 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
857 '''
858 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
860 if _FOR_DOCS:
861 forward = _GnomonicBase.forward
862 reverse = _GnomonicBase.reverse
864 @Property_RO
865 def geodesic(self):
866 '''Get this projection's exact geodesic (L{GeodesicExact}).
867 '''
868 return self.datum.ellipsoid.geodesicx
871class GnomonicGeodSolve(_GnomonicBase):
872 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
873 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
874 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and
875 intended I{for testing purposes only}.
877 @see: L{GnomonicExact} and module L{geodsolve}.
878 '''
879 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
880 '''New azimuthal L{GnomonicGeodSolve} projection.
882 @arg lat0: Latitude of center point (C{degrees90}).
883 @arg lon0: Longitude of center point (C{degrees180}).
884 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
885 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
886 radius (C{meter}).
887 @kwarg name: Optional name for the projection (C{str}).
889 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
890 '''
891 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
893 if _FOR_DOCS:
894 forward = _GnomonicBase.forward
895 reverse = _GnomonicBase.reverse
897 @Property_RO
898 def geodesic(self):
899 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
900 '''
901 return self.datum.ellipsoid.geodsolve
904class GnomonicKarney(_GnomonicBase):
905 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
906 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
907 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
909 @see: L{GnomonicExact}.
910 '''
911 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
912 '''New azimuthal L{GnomonicKarney} projection.
914 @arg lat0: Latitude of center point (C{degrees90}).
915 @arg lon0: Longitude of center point (C{degrees180}).
916 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
917 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
918 radius (C{meter}).
919 @kwarg name: Optional name for the projection (C{str}).
921 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
922 not installed or not found.
924 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
925 '''
926 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
928 if _FOR_DOCS:
929 forward = _GnomonicBase.forward
930 reverse = _GnomonicBase.reverse
932 @Property_RO
933 def geodesic(self):
934 '''Get this projection's I{wrapped} U{Karney Geodesic
935 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}, provided package
936 U{geographiclib<https://PyPI.org/project/geographiclib>} is installed.
937 '''
938 return self.datum.ellipsoid.geodesic
941class LambertEqualArea(_AzimuthalBase):
942 '''Lambert-equal-area projection for the sphere*** (aka U{Lambert zenithal equal-area
943 projection<https://WikiPedia.org/wiki/Lambert_azimuthal_equal-area_projection>}, see
944 U{Snyder, pp 185-187<https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
945 <https://MathWorld.Wolfram.com/LambertAzimuthalEqual-AreaProjection.html>}.
946 '''
947 if _FOR_DOCS:
948 __init__ = _AzimuthalBase.__init__
950 def forward(self, lat, lon, name=NN):
951 '''Convert a geodetic location to azimuthal Lambert-equal-area east- and northing.
953 @arg lat: Latitude of the location (C{degrees90}).
954 @arg lon: Longitude of the location (C{degrees180}).
955 @kwarg name: Optional name for the location (C{str}).
957 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
958 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
959 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
960 The C{scale} of the projection is C{1} in I{radial} direction and
961 is C{1 / reciprocal} in the direction perpendicular to this.
963 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
964 '''
965 def _k_t(c):
966 c += _1_0
967 t = c > EPS0
968 k = sqrt(_2_0 / c) if t else _1_0
969 return k, t
971 return self._forward(lat, lon, name, _k_t)
973 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
974 '''Convert an azimuthal Lambert-equal-area location to geodetic lat- and longitude.
976 @arg x: Easting of the location (C{meter}).
977 @arg y: Northing of the location (C{meter}).
978 @kwarg name: Optional name for the location (C{str}).
979 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
980 additional B{C{LatLon}} keyword arguments,
981 ignored if C{B{LatLon} is None} or not given.
983 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
984 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
986 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
987 in the range C{[-180..180] degrees}. The C{scale} of the
988 projection is C{1} in I{radial} direction, C{azimuth} clockwise
989 from true North and is C{1 / reciprocal} in the direction
990 perpendicular to this.
991 '''
992 def _c(c):
993 c *= _0_5
994 return (asin1(c) * _2_0) if c > EPS else None
996 return self._reverse(x, y, name, _c, True, **LatLon_and_kwds)
999class Orthographic(_AzimuthalBase):
1000 '''Orthographic projection for the sphere***, see U{Snyder, pp 148-153
1001 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1002 <https://MathWorld.Wolfram.com/OrthographicProjection.html>}.
1003 '''
1004 if _FOR_DOCS:
1005 __init__ = _AzimuthalBase.__init__
1007 def forward(self, lat, lon, name=NN):
1008 '''Convert a geodetic location to azimuthal orthographic east- and northing.
1010 @arg lat: Latitude of the location (C{degrees90}).
1011 @arg lon: Longitude of the location (C{degrees180}).
1012 @kwarg name: Optional name for the location (C{str}).
1014 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1015 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1016 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1017 The C{scale} of the projection is C{1} in I{radial} direction and
1018 is C{1 / reciprocal} in the direction perpendicular to this.
1020 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1021 '''
1022 def _k_t(c):
1023 return _1_0, (c >= 0)
1025 return self._forward(lat, lon, name, _k_t)
1027 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1028 '''Convert an azimuthal orthographic location to geodetic lat- and longitude.
1030 @arg x: Easting of the location (C{meter}).
1031 @arg y: Northing of the location (C{meter}).
1032 @kwarg name: Optional name for the location (C{str}).
1033 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1034 additional B{C{LatLon}} keyword arguments,
1035 ignored if C{B{LatLon} is None} or not given.
1037 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1038 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1040 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1041 in the range C{[-180..180] degrees}. The C{scale} of the
1042 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1043 from true North and is C{1 / reciprocal} in the direction
1044 perpendicular to this.
1045 '''
1046 def _c(c):
1047 return asin1(c) if c > EPS else None
1049 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1052class Stereographic(_AzimuthalBase):
1053 '''Stereographic projection for the sphere***, see U{Snyder, pp 157-160
1054 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1055 <https://MathWorld.Wolfram.com/StereographicProjection.html>}.
1056 '''
1057 _k0 = _1_0 # central scale factor (C{scalar})
1058 _k02 = _2_0 # double ._k0
1060 if _FOR_DOCS:
1061 __init__ = _AzimuthalBase.__init__
1063 def forward(self, lat, lon, name=NN):
1064 '''Convert a geodetic location to azimuthal stereographic east- and northing.
1066 @arg lat: Latitude of the location (C{degrees90}).
1067 @arg lon: Longitude of the location (C{degrees180}).
1068 @kwarg name: Optional name for the location (C{str}).
1070 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1071 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1072 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1073 The C{scale} of the projection is C{1} in I{radial} direction and
1074 is C{1 / reciprocal} in the direction perpendicular to this.
1076 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1077 '''
1078 def _k_t(c):
1079 c += _1_0
1080 t = isnon0(c)
1081 k = (self._k02 / c) if t else _1_0
1082 return k, t
1084 return self._forward(lat, lon, name, _k_t)
1086 @property_doc_(''' optional, central scale factor (C{scalar}).''')
1087 def k0(self):
1088 '''Get the central scale factor (C{scalar}).
1089 '''
1090 return self._k0
1092 @k0.setter # PYCHOK setter!
1093 def k0(self, factor):
1094 '''Set the central scale factor (C{scalar}).
1095 '''
1096 n = Stereographic.k0.fget.__name__
1097 self._k0 = Scalar_(factor, name=n, low=EPS, high=2) # XXX high=1, 2, other?
1098 self._k02 = self._k0 * _2_0
1100 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1101 '''Convert an azimuthal stereographic location to geodetic lat- and longitude.
1103 @arg x: Easting of the location (C{meter}).
1104 @arg y: Northing of the location (C{meter}).
1105 @kwarg name: Optional name for the location (C{str}).
1106 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1107 additional B{C{LatLon}} keyword arguments,
1108 ignored if C{B{LatLon} is None} or not given.
1110 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1111 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1113 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1114 in the range C{[-180..180] degrees}. The C{scale} of the
1115 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1116 from true North and is C{1 / reciprocal} in the direction
1117 perpendicular to this.
1118 '''
1119 def _c(c):
1120 return (_2_0 * atan2(c, self._k02)) if c > EPS else None
1122 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1125__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase)
1127# **) MIT License
1128#
1129# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1130#
1131# Permission is hereby granted, free of charge, to any person obtaining a
1132# copy of this software and associated documentation files (the "Software"),
1133# to deal in the Software without restriction, including without limitation
1134# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1135# and/or sell copies of the Software, and to permit persons to whom the
1136# Software is furnished to do so, subject to the following conditions:
1137#
1138# The above copyright notice and this permission notice shall be included
1139# in all copies or substantial portions of the Software.
1140#
1141# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1142# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1143# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1144# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1145# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1146# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1147# OTHER DEALINGS IN THE SOFTWARE.