Coverage for pygeodesy/azimuthal.py: 98%
312 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -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.07.10'
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 the
442 I{wrapped} U{geodesic.Geodesic and geodesicline.GeodesicLine
443 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} or the
444 I{exact} geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
445 '''
446 _mask = 0
448 @Property_RO
449 def geodesic(self): # PYCHOK no cover
450 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
451 '''
452 notOverloaded(self)
454 def _7Tuple(self, e, n, r, M=None, name=NN):
455 '''(INTERNAL) Return an C{Azimuthal7Tuple}.
456 '''
457 s = M if M is not None else ( # reciprocal, azimuthal scale
458 (r.m12 / r.s12) if r.a12 > _EPS_K else _1_0)
459 z = _umod_360(r.azi2) # -180 <= r.azi2 < 180 ... 0 <= z < 360
460 return Azimuthal7Tuple(e, n, r.lat2, r.lon2, z, s, self.datum,
461 name=name or self.name)
464class _EquidistantBase(_AzimuthalGeodesic):
465 '''(INTERNAL) Base for classes L{EquidistantExact}, L{EquidistantGeodSolve}
466 and L{EquidistantKarney}.
467 '''
468 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
469 '''New azimuthal L{EquidistantExact}, L{EquidistantGeodSolve} or
470 L{EquidistantKarney} projection.
471 '''
472 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
474 g = self.geodesic
475 # g.STANDARD = g.AZIMUTH | g.DISTANCE | g.LATITUDE | g.LONGITUDE
476 self._mask = g.REDUCEDLENGTH | g.STANDARD # | g.LONG_UNROLL
478 def forward(self, lat, lon, name=NN):
479 '''Convert an (ellipsoidal) geodetic location to azimuthal equidistant east- and northing.
481 @arg lat: Latitude of the location (C{degrees90}).
482 @arg lon: Longitude of the location (C{degrees180}).
483 @kwarg name: Optional name for the location (C{str}).
485 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
486 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
487 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
488 The C{scale} of the projection is C{1} in I{radial} direction and
489 is C{1 / reciprocal} in the direction perpendicular to this.
491 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
493 @note: A call to C{.forward} followed by a call to C{.reverse} will return
494 the original C{lat, lon} to within roundoff.
495 '''
496 r = self.geodesic.Inverse(self.lat0, self.lon0,
497 Lat_(lat), Lon_(lon), outmask=self._mask)
498 x, y = sincos2d(r.azi1)
499 return self._7Tuple(x * r.s12, y * r.s12, r, name=name)
501 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
502 '''Convert an azimuthal equidistant location to (ellipsoidal) geodetic lat- and longitude.
504 @arg x: Easting of the location (C{meter}).
505 @arg y: Northing of the location (C{meter}).
506 @kwarg name: Optional name for the location (C{str}).
507 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
508 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
509 arguments, ignored if C{B{LatLon} is None}.
511 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
512 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
514 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
515 in the range C{[-180..180] degrees}. The scale of the projection
516 is C{1} in I{radial} direction, C{azimuth} clockwise from true
517 North and is C{1 / reciprocal} in the direction perpendicular
518 to this.
519 '''
520 e, n, z, s = _enzh4(x, y)
522 r = self.geodesic.Direct(self.lat0, self.lon0, z, s, outmask=self._mask)
523 return self._7Tuple(e, n, r, name=name) if LatLon is None else \
524 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
527class EquidistantExact(_EquidistantBase):
528 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
529 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
530 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
532 An azimuthal equidistant projection is centered at an arbitrary position on the ellipsoid.
533 For a point in projected space C{(x, y)}, the geodesic distance from the center position
534 is C{hypot(x, y)} and the C{azimuth} of the geodesic from the center point is C{atan2(x, y)},
535 clockwise from true North.
537 The C{.forward} and C{.reverse} methods also return the C{azimuth} of the geodesic at C{(x,
538 y)} and the C{scale} in the azimuthal direction which, together with the basic properties
539 of the projection, serve to specify completely the local affine transformation between
540 geographic and projected coordinates.
541 '''
542 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
543 '''New azimuthal L{EquidistantExact} projection.
545 @arg lat0: Latitude of center point (C{degrees90}).
546 @arg lon0: Longitude of center point (C{degrees180}).
547 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
548 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
549 radius (C{meter}).
550 @kwarg name: Optional name for the projection (C{str}).
552 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
553 '''
554 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
556 if _FOR_DOCS:
557 forward = _EquidistantBase.forward
558 reverse = _EquidistantBase.reverse
560 @Property_RO
561 def geodesic(self):
562 '''Get this projection's exact geodesic (L{GeodesicExact}).
563 '''
564 return self.datum.ellipsoid.geodesicx
567class EquidistantGeodSolve(_EquidistantBase):
568 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
569 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
570 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and intended
571 I{for testing purposes only}.
573 @see: L{EquidistantExact} and module L{geodsolve}.
574 '''
575 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
576 '''New azimuthal L{EquidistantGeodSolve} projection.
578 @arg lat0: Latitude of center point (C{degrees90}).
579 @arg lon0: Longitude of center point (C{degrees180}).
580 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
581 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
582 radius (C{meter}).
583 @kwarg name: Optional name for the projection (C{str}).
585 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
586 '''
587 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
589 if _FOR_DOCS:
590 forward = _EquidistantBase.forward
591 reverse = _EquidistantBase.reverse
593 @Property_RO
594 def geodesic(self):
595 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
596 '''
597 return self.datum.ellipsoid.geodsolve
600class EquidistantKarney(_EquidistantBase):
601 '''Azimuthal equidistant projection, a Python version of I{Karney}'s C++ class U{AzimuthalEquidistant
602 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>},
603 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
605 @see: L{EquidistantExact}.
606 '''
607 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
608 '''New azimuthal L{EquidistantKarney} projection.
610 @arg lat0: Latitude of center point (C{degrees90}).
611 @arg lon0: Longitude of center point (C{degrees180}).
612 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
613 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
614 radius (C{meter}).
615 @kwarg name: Optional name for the projection (C{str}).
617 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
618 not installed or not found.
620 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or B{C{datum}}.
621 '''
622 _EquidistantBase.__init__(self, lat0, lon0, datum=datum, name=name)
624 if _FOR_DOCS:
625 forward = _EquidistantBase.forward
626 reverse = _EquidistantBase.reverse
628 @Property_RO
629 def geodesic(self):
630 '''Get this projection's I{wrapped} U{geodesic.Geodesic
631 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
632 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
633 package is installed.
634 '''
635 return self.datum.ellipsoid.geodesic
638_Equidistants = (Equidistant, EquidistantExact, EquidistantGeodSolve,
639 EquidistantKarney) # PYCHOK in .ellipsoidalBaseDI
642class Gnomonic(_AzimuthalBase):
643 '''Azimuthal gnomonic projection for the sphere***, see U{Snyder, pp 164-168
644 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
645 <https://MathWorld.Wolfram.com/GnomonicProjection.html>}.
646 '''
647 if _FOR_DOCS:
648 __init__ = _AzimuthalBase.__init__
650 def forward(self, lat, lon, name=NN):
651 '''Convert a geodetic location to azimuthal equidistant east- and northing.
653 @arg lat: Latitude of the location (C{degrees90}).
654 @arg lon: Longitude of the location (C{degrees180}).
655 @kwarg name: Optional name for the location (C{str}).
657 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
658 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
659 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
660 The C{scale} of the projection is C{1} in I{radial} direction and
661 is C{1 / reciprocal} in the direction perpendicular to this.
663 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
664 '''
665 def _k_t(c):
666 t = c > EPS
667 k = (_1_0 / c) if t else _1_0
668 return k, t
670 return self._forward(lat, lon, name, _k_t)
672 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
673 '''Convert an azimuthal equidistant location to geodetic lat- and longitude.
675 @arg x: Easting of the location (C{meter}).
676 @arg y: Northing of the location (C{meter}).
677 @kwarg name: Optional name for the location (C{str}).
678 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
679 additional B{C{LatLon}} keyword arguments,
680 ignored if C{B{LatLon} is None} or not given.
682 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
683 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
685 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
686 in the range C{[-180..180] degrees}. The C{scale} of the
687 projection is C{1} in I{radial} direction, C{azimuth} clockwise
688 from true North and C{1 / reciprocal} in the direction
689 perpendicular to this.
690 '''
691 def _c(c):
692 return atan(c) if c > EPS else None
694 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
697def gnomonic(lat0, lon0, datum=_WGS84, exact=False, geodsolve=False, name=NN):
698 '''Return a L{GnomonicExact} or (if I{Karney}'s U{geographiclib
699 <https://PyPI.org/project/geographiclib>} package is installed)
700 a L{GnomonicKarney}, otherwise a L{Gnomonic} instance.
702 @arg lat0: Latitude of center point (C{degrees90}).
703 @arg lon0: Longitude of center point (C{degrees180}).
704 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
705 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
706 radius (C{meter}).
707 @kwarg exact: Return a L{GnomonicExact} instance.
708 @kwarg geodsolve: Return a L{GnomonicGeodSolve} instance.
709 @kwarg name: Optional name for the projection (C{str}).
711 @return: A L{GnomonicExact}, L{GnomonicGeodSolve},
712 L{GnomonicKarney} or L{Gnomonic} instance.
714 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}} or
715 (spherical) B{C{datum}}.
717 @raise GeodesicError: Issue with L{GeodesicExact}, L{GeodesicSolve}
718 or I{Karney}'s wrapped C{Geodesic}.
720 @raise TypeError: Invalid B{C{datum}}.
721 '''
722 G = GnomonicExact if exact else (GnomonicGeodSolve if geodsolve else Gnomonic)
723 if G is Gnomonic:
724 try:
725 return GnomonicKarney(lat0, lon0, datum=datum, name=name) # PYCHOK types
726 except ImportError:
727 pass
728 return G(lat0, lon0, datum=datum, name=name) # PYCHOK types
731class _GnomonicBase(_AzimuthalGeodesic):
732 '''(INTERNAL) Base for classes L{GnomonicExact}, L{GnomonicGeodSolve}
733 and L{GnomonicKarney}.
734 '''
735 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
736 '''New azimuthal L{GnomonicExact} or L{GnomonicKarney} projection.
737 '''
738 _AzimuthalGeodesic.__init__(self, lat0, lon0, datum=datum, name=name)
740 g = self.geodesic
741 self._mask = g.ALL # | g.LONG_UNROLL
743 def forward(self, lat, lon, name=NN, raiser=True): # PYCHOK signature
744 '''Convert an (ellipsoidal) geodetic location to azimuthal gnomonic east-
745 and northing.
747 @arg lat: Latitude of the location (C{degrees90}).
748 @arg lon: Longitude of the location (C{degrees180}).
749 @kwarg name: Optional name for the location (C{str}).
750 @kwarg raiser: Do or don't throw an error (C{bool}) if
751 the location lies over the horizon.
753 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
754 with easting C{x} and northing C{y} in C{meter} and C{lat} and
755 C{lon} in C{degrees} and C{azimuth} clockwise from true North.
756 The C{scale} of the projection is C{1 / reciprocal**2} in I{radial}
757 direction and C{1 / reciprocal} in the direction perpendicular to
758 this. Both C{x} and C{y} will be C{NAN} if the (geodetic) location
759 lies over the horizon and C{B{raiser}=False}.
761 @raise AzimuthalError: Invalid B{C{lat}}, B{C{lon}} or the location lies
762 over the horizon and C{B{raiser}=True}.
763 '''
764 self._iteration = 0
766 r = self.geodesic.Inverse(self.lat0, self.lon0,
767 Lat_(lat), Lon_(lon), outmask=self._mask)
768 M = r.M21
769 if M > EPS0:
770 q = r.m12 / M # .M12
771 e, n = sincos2d(r.azi1)
772 e *= q
773 n *= q
774 elif raiser: # PYCHOK no cover
775 raise AzimuthalError(lat=lat, lon=lon, txt=_over_horizon_)
776 else: # PYCHOK no cover
777 e = n = NAN
779 t = self._7Tuple(e, n, r, M=M, name=name)
780 t._iteraton = 0
781 return t
783 def reverse(self, x, y, name=NN, LatLon=None, **LatLon_kwds): # PYCHOK signature
784 '''Convert an azimuthal gnomonic location to (ellipsoidal) geodetic lat- and longitude.
786 @arg x: Easting of the location (C{meter}).
787 @arg y: Northing of the location (C{meter}).
788 @kwarg name: Optional name for the location (C{str}).
789 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
790 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
791 arguments, ignored if C{B{LatLon} is None}.
793 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
794 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
796 @raise AzimuthalError: No convergence.
798 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
799 in the range C{[-180..180] degrees}. The C{azimuth} is clockwise
800 from true North. The scale is C{1 / reciprocal**2} in C{radial}
801 direction and C{1 / reciprocal} in the direction perpendicular
802 to this.
803 '''
804 e, n, z, q = _enzh4(x, y)
806 d = a = self.equatoradius
807 s = a * atan(q / a)
808 if q > a: # PYCHOK no cover
809 def _d(r, q):
810 return (r.M12 - q * r.m12) * r.m12 # negated
812 q = _1_0 / q
813 else: # little == True
814 def _d(r, q): # PYCHOK _d
815 return (q * r.M12 - r.m12) * r.M12 # negated
817 a *= _EPS_K
818 m = self._mask
819 g = self.geodesic
821 _P = g.Line(self.lat0, self.lon0, z, m | g.LINE_OFF).Position
822 _S2 = Fsum(s).fsum2_
823 for i in range(1, _TRIPS):
824 r = _P(s, outmask=m)
825 if fabs(d) < a:
826 break
827 s, d = _S2(_d(r, q))
828 else: # PYCHOK no cover
829 self._iteration = _TRIPS
830 raise AzimuthalError(Fmt.no_convergence(d, a),
831 txt=unstr(self.reverse, x, y))
833 t = self._7Tuple(e, n, r, M=r.M12, name=name) if LatLon is None else \
834 self._toLatLon(r.lat2, r.lon2, LatLon, LatLon_kwds, name)
835 t._iteration = self._iteration = i
836 return t
839class GnomonicExact(_GnomonicBase):
840 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
841 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
842 based on exact geodesic classes L{GeodesicExact} and L{GeodesicLineExact}.
844 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/
845 classGeographicLib_1_1Gnomonic.html>}, especially the B{Warning}.
846 '''
847 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
848 '''New azimuthal L{GnomonicExact} projection.
850 @arg lat0: Latitude of center point (C{degrees90}).
851 @arg lon0: Longitude of center point (C{degrees180}).
852 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
853 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
854 radius (C{meter}).
855 @kwarg name: Optional name for the projection (C{str}).
857 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
858 '''
859 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
861 if _FOR_DOCS:
862 forward = _GnomonicBase.forward
863 reverse = _GnomonicBase.reverse
865 @Property_RO
866 def geodesic(self):
867 '''Get this projection's exact geodesic (L{GeodesicExact}).
868 '''
869 return self.datum.ellipsoid.geodesicx
872class GnomonicGeodSolve(_GnomonicBase):
873 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
874 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
875 based on (exact) geodesic I{wrappers} L{GeodesicSolve} and L{GeodesicLineSolve} and
876 intended I{for testing purposes only}.
878 @see: L{GnomonicExact} and module L{geodsolve}.
879 '''
880 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
881 '''New azimuthal L{GnomonicGeodSolve} projection.
883 @arg lat0: Latitude of center point (C{degrees90}).
884 @arg lon0: Longitude of center point (C{degrees180}).
885 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
886 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
887 radius (C{meter}).
888 @kwarg name: Optional name for the projection (C{str}).
890 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
891 '''
892 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
894 if _FOR_DOCS:
895 forward = _GnomonicBase.forward
896 reverse = _GnomonicBase.reverse
898 @Property_RO
899 def geodesic(self):
900 '''Get this projection's (exact) geodesic (L{GeodesicSolve}).
901 '''
902 return self.datum.ellipsoid.geodsolve
905class GnomonicKarney(_GnomonicBase):
906 '''Azimuthal gnomonic projection, a Python version of I{Karney}'s C++ class U{Gnomonic
907 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>},
908 requiring package U{geographiclib<https://PyPI.org/project/geographiclib>} to be installed.
910 @see: L{GnomonicExact}.
911 '''
912 def __init__(self, lat0, lon0, datum=_WGS84, name=NN):
913 '''New azimuthal L{GnomonicKarney} projection.
915 @arg lat0: Latitude of center point (C{degrees90}).
916 @arg lon0: Longitude of center point (C{degrees180}).
917 @kwarg datum: Optional datum or ellipsoid (L{Datum}, L{Ellipsoid},
918 L{Ellipsoid2} or L{a_f2Tuple}) or I{scalar} earth
919 radius (C{meter}).
920 @kwarg name: Optional name for the projection (C{str}).
922 @raise ImportError: Package U{geographiclib<https://PyPI.org/project/geographiclib>}
923 not installed or not found.
925 @raise AzimuthalError: Invalid B{C{lat0}} or B{C{lon0}}.
926 '''
927 _GnomonicBase.__init__(self, lat0, lon0, datum=datum, name=name)
929 if _FOR_DOCS:
930 forward = _GnomonicBase.forward
931 reverse = _GnomonicBase.reverse
933 @Property_RO
934 def geodesic(self):
935 '''Get this projection's I{wrapped} U{geodesic.Geodesic
936 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
937 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
938 package is installed.
939 '''
940 return self.datum.ellipsoid.geodesic
943class LambertEqualArea(_AzimuthalBase):
944 '''Lambert-equal-area projection for the sphere*** (aka U{Lambert zenithal equal-area
945 projection<https://WikiPedia.org/wiki/Lambert_azimuthal_equal-area_projection>}, see
946 U{Snyder, pp 185-187<https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
947 <https://MathWorld.Wolfram.com/LambertAzimuthalEqual-AreaProjection.html>}.
948 '''
949 if _FOR_DOCS:
950 __init__ = _AzimuthalBase.__init__
952 def forward(self, lat, lon, name=NN):
953 '''Convert a geodetic location to azimuthal Lambert-equal-area east- and northing.
955 @arg lat: Latitude of the location (C{degrees90}).
956 @arg lon: Longitude of the location (C{degrees180}).
957 @kwarg name: Optional name for the location (C{str}).
959 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
960 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
961 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
962 The C{scale} of the projection is C{1} in I{radial} direction and
963 is C{1 / reciprocal} in the direction perpendicular to this.
965 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
966 '''
967 def _k_t(c):
968 c += _1_0
969 t = c > EPS0
970 k = sqrt(_2_0 / c) if t else _1_0
971 return k, t
973 return self._forward(lat, lon, name, _k_t)
975 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
976 '''Convert an azimuthal Lambert-equal-area location to geodetic lat- and longitude.
978 @arg x: Easting of the location (C{meter}).
979 @arg y: Northing of the location (C{meter}).
980 @kwarg name: Optional name for the location (C{str}).
981 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
982 additional B{C{LatLon}} keyword arguments,
983 ignored if C{B{LatLon} is None} or not given.
985 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
986 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
988 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
989 in the range C{[-180..180] degrees}. The C{scale} of the
990 projection is C{1} in I{radial} direction, C{azimuth} clockwise
991 from true North and is C{1 / reciprocal} in the direction
992 perpendicular to this.
993 '''
994 def _c(c):
995 c *= _0_5
996 return (asin1(c) * _2_0) if c > EPS else None
998 return self._reverse(x, y, name, _c, True, **LatLon_and_kwds)
1001class Orthographic(_AzimuthalBase):
1002 '''Orthographic projection for the sphere***, see U{Snyder, pp 148-153
1003 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1004 <https://MathWorld.Wolfram.com/OrthographicProjection.html>}.
1005 '''
1006 if _FOR_DOCS:
1007 __init__ = _AzimuthalBase.__init__
1009 def forward(self, lat, lon, name=NN):
1010 '''Convert a geodetic location to azimuthal orthographic east- and northing.
1012 @arg lat: Latitude of the location (C{degrees90}).
1013 @arg lon: Longitude of the location (C{degrees180}).
1014 @kwarg name: Optional name for the location (C{str}).
1016 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1017 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1018 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1019 The C{scale} of the projection is C{1} in I{radial} direction and
1020 is C{1 / reciprocal} in the direction perpendicular to this.
1022 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1023 '''
1024 def _k_t(c):
1025 return _1_0, (c >= 0)
1027 return self._forward(lat, lon, name, _k_t)
1029 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1030 '''Convert an azimuthal orthographic location to geodetic lat- and longitude.
1032 @arg x: Easting of the location (C{meter}).
1033 @arg y: Northing of the location (C{meter}).
1034 @kwarg name: Optional name for the location (C{str}).
1035 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1036 additional B{C{LatLon}} keyword arguments,
1037 ignored if C{B{LatLon} is None} or not given.
1039 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1040 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1042 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1043 in the range C{[-180..180] degrees}. The C{scale} of the
1044 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1045 from true North and is C{1 / reciprocal} in the direction
1046 perpendicular to this.
1047 '''
1048 def _c(c):
1049 return asin1(c) if c > EPS else None
1051 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1054class Stereographic(_AzimuthalBase):
1055 '''Stereographic projection for the sphere***, see U{Snyder, pp 157-160
1056 <https://Pubs.USGS.gov/pp/1395/report.pdf>} and U{MathWorld-Wolfram
1057 <https://MathWorld.Wolfram.com/StereographicProjection.html>}.
1058 '''
1059 _k0 = _1_0 # central scale factor (C{scalar})
1060 _k02 = _2_0 # double ._k0
1062 if _FOR_DOCS:
1063 __init__ = _AzimuthalBase.__init__
1065 def forward(self, lat, lon, name=NN):
1066 '''Convert a geodetic location to azimuthal stereographic east- and northing.
1068 @arg lat: Latitude of the location (C{degrees90}).
1069 @arg lon: Longitude of the location (C{degrees180}).
1070 @kwarg name: Optional name for the location (C{str}).
1072 @return: An L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}
1073 with easting C{x} and northing C{y} of point in C{meter} and C{lat}
1074 and C{lon} in C{degrees} and C{azimuth} clockwise from true North.
1075 The C{scale} of the projection is C{1} in I{radial} direction and
1076 is C{1 / reciprocal} in the direction perpendicular to this.
1078 @raise AzimuthalError: Invalid B{C{lat}} or B{C{lon}}.
1079 '''
1080 def _k_t(c):
1081 c += _1_0
1082 t = isnon0(c)
1083 k = (self._k02 / c) if t else _1_0
1084 return k, t
1086 return self._forward(lat, lon, name, _k_t)
1088 @property_doc_(''' optional, central scale factor (C{scalar}).''')
1089 def k0(self):
1090 '''Get the central scale factor (C{scalar}).
1091 '''
1092 return self._k0
1094 @k0.setter # PYCHOK setter!
1095 def k0(self, factor):
1096 '''Set the central scale factor (C{scalar}).
1097 '''
1098 n = Stereographic.k0.fget.__name__
1099 self._k0 = Scalar_(factor, name=n, low=EPS, high=2) # XXX high=1, 2, other?
1100 self._k02 = self._k0 * _2_0
1102 def reverse(self, x, y, name=NN, **LatLon_and_kwds):
1103 '''Convert an azimuthal stereographic location to geodetic lat- and longitude.
1105 @arg x: Easting of the location (C{meter}).
1106 @arg y: Northing of the location (C{meter}).
1107 @kwarg name: Optional name for the location (C{str}).
1108 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use and
1109 additional B{C{LatLon}} keyword arguments,
1110 ignored if C{B{LatLon} is None} or not given.
1112 @return: The geodetic (C{LatLon}) or if B{C{LatLon}} is C{None} an
1113 L{Azimuthal7Tuple}C{(x, y, lat, lon, azimuth, scale, datum)}.
1115 @note: The C{lat} will be in the range C{[-90..90] degrees} and C{lon}
1116 in the range C{[-180..180] degrees}. The C{scale} of the
1117 projection is C{1} in I{radial} direction, C{azimuth} clockwise
1118 from true North and is C{1 / reciprocal} in the direction
1119 perpendicular to this.
1120 '''
1121 def _c(c):
1122 return (_2_0 * atan2(c, self._k02)) if c > EPS else None
1124 return self._reverse(x, y, name, _c, False, **LatLon_and_kwds)
1127__all__ += _ALL_DOCS(_AzimuthalBase, _AzimuthalGeodesic, _EquidistantBase, _GnomonicBase)
1129# **) MIT License
1130#
1131# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1132#
1133# Permission is hereby granted, free of charge, to any person obtaining a
1134# copy of this software and associated documentation files (the "Software"),
1135# to deal in the Software without restriction, including without limitation
1136# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1137# and/or sell copies of the Software, and to permit persons to whom the
1138# Software is furnished to do so, subject to the following conditions:
1139#
1140# The above copyright notice and this permission notice shall be included
1141# in all copies or substantial portions of the Software.
1142#
1143# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1144# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1145# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1146# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1147# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1148# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1149# OTHER DEALINGS IN THE SOFTWARE.