Coverage for pygeodesy/geodesicx/gxline.py : 98%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicLineExact.html>}.
Copyright (C) Charles Karney (2012-2021) <Charles@Karney.com> and licensed under the MIT/X11 License. For more information, see U{GeographicLib<https://GeographicLib.SourceForge.io>}. ''' # make sure int/int division yields float quotient
# A copy of comments from Karney's C{GeodesicLineExact.cpp}: # # This is a reformulation of the geodesic problem. The # notation is as follows: # - at a general point (no suffix or 1 or 2 as suffix) # - phi = latitude # - beta = latitude on auxiliary sphere # - omega = longitude on auxiliary sphere # - lambda = longitude # - alpha = azimuth of great circle # - sigma = arc length along great circle # - s = distance # - tau = scaled distance (= sigma at multiples of PI/2) # - at northwards equator crossing # - beta = phi = 0 # - omega = lambda = 0 # - alpha = alpha0 # - sigma = s = 0 # - a 12 suffix means a difference, e.g., s12 = s2 - s1. # - s and c prefixes mean sin and cos
# from pygeodesy.errors import _AssertionError # from .karney _coSeries, _GeodesicBase, \ _sincos12, _TINY _0_0, _1_0, _90_0, _180_0 _copysign, _fix90, GDict, _hypot, \ _norm2, _norm180, _sincos2, _sincos2d
'''(INTERNAL) Zap cached/memoized C{Property}s of any L{GeodesicLineExact} instances tied to the given L{GeodesicExact} instance B{C{gX}}. ''' raise _AssertionError(gX=gX)
'''(INTERNAL) Base class for L{GeodesicLineExact}. '''
'''(INTERNAL) New C{[_]GeodesicLineExact} instance. ''' else: # guard against salp0 underflow, # also -0 is converted to +0
if _debug: # PYCHOK no cover caps |= Caps._DEBUG_LINE & _debug # allow lat, azimuth and unrolling of lon
# Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), # Alt: calp0 = _hypot(sbet1, calp1 * cbet1). The following # is slightly better (consider the case salp1 = 0). # Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1). # sig = 0 is nearest northward crossing of equator. # With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line). # With bet1 = pi/2, alp1 = -pi, sig1 = pi/2 # With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2 # Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1). # With alp0 in (0, pi/2], quadrants for sig and omg coincide. # No atan2(0,0) ambiguity at poles since cbet1 = +epsilon. # With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. # Without normalization we have schi1 = somg1. # _norm2(somg1, comg1) # no need to normalize! # _norm2(schi1?, cchi1) # no need to normalize!
# no need to pre-compute other attrs based on _Caps.X. All are # Property_RO's, computed once and cached/memoized until reset # when C4Order is changed or Elliptic function reset is invoked.
try: # PYCHOK no cover _glXs.remove(self) except (TypeError, ValueError): pass # _update_all(self) # throws TypeError during Python 2 cleanup
'''Get the I{equatorial arc} (C{degrees}), the arc length between the northward equatorial crossing and the first point. '''
'''Get (spherical) arc length from the first to the reference point (C{degrees}).
@see: Method L{SetArc}. '''
'''(INTERNAL) Cached/memoized. '''
'''Find the position on the line given B{C{a12}}.
@arg a12: Spherical arc length from the first point to the second point (C{degrees}). @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying the quantities to be returned.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and arc length C{a12} always included, except when C{a12=NAN}.
@note: By default, C{B{outmask}=STANDARD}, meaning thc C{lat1}, C{lon1}, C{azi1}, C{lat2}, C{lon2}, C{azi2}, C{s12} and C{a12} entries are returned, except when C{a12=NAN}. '''
'''Get the I{equatorial azimuth}, the azimuth of this geodesic line as it crosses the equator in a northward direction (C{degrees90}). '''
'''Get the sine and cosine of the I{equatorial azimuth} (2-tuple C{(sin, cos)}). '''
'''Get the azimuth at the first point (compass C{degrees}). '''
'''Get the sine and cosine of the first point's azimuth (2-tuple C{(sin, cos)}). ''' return self._salp1, self._calp1
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized. '''
'''Get the capabilities (bit-or'ed C{Caps}). ''' return self._caps
'''Check the available capabilities.
@arg caps: Bit-or'ed combination of L{Caps} values for all capabilities to be checked.
@return: C{True} if I{all} B{C{caps}} are available, C{False} otherwise (C{bool}). ''' return _all_caps(self.caps, caps)
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized C{Elliptic} function. ''' # see .gx.GeodesicExact._ef_reset_k2
'''(INTERNAL) Generate a new position along the geodesic.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and arc length C{a12} always included, except when C{a12=NAN}. ''' if not (arcmode or self._caps_DISTANCE_IN): # PYCHOK no cover return r # Uninitialized or impossible distance requested
if self.debug: # PYCHOK no cover outmask |= Caps._DEBUG_LINE & self._debug # f_0_01 = False
if arcmode: # PYCHOK no cover # s12_a12 is spherical arc length E2 = _0_0 sig12 = radians(s12_a12) ssig12, csig12 = _sincos2(sig12) a = abs(s12_a12) a -= floor(a / _180_0) * _180_0 if a == _0_0: ssig12 = _0_0 elif a == _90_0: csig12 = _0_0 else: # s12_a12 is distance # tau2 = tau1 + tau12 # if _K_2_0 and abs(gX.f) > _0_01: # f_0_01 = True # raise NotYetImplementedError()
# sig2 = sig1 + sig12 # sin(bet2) = cos(alp0) * sin(sig2) # Alt: cbet2 = _hypot(csig2, salp0 * ssig2); cbet2 = csig2 = _TINY # tan(alp0) = cos(sig2) * tan(alp2)
# AB1 = _E0 * (E2 - _E1) # s12 = _b * (_E0 * sig12 + AB1) # = _b * _E0 * (sig12 + (E2 - _E1)) # = _b * _E0 * (E2 - _E1 + sig12) else:
if (outmask & Caps._DEBUG_LINE): # PYCHOK no cover r.set_(sig12=sig12, dn2=dn2, E0_b=self._E0_b, E1=self._E1, E2=E2, e2=gX.e2, f1=gX.f1, eFk2=eF.k2, eFa2=eF.alpha2)
-self._H1, sig12) atan2(tchi2, cchi2), -atan2(tchi1, cchi1), sig12) else: if (outmask & Caps._DEBUG_LINE): # PYCHOK no cover r.set_(chi12=chi12, lam12=lam12, H0_e2_f1=self._H0_e2_f1, H1=self._H1, ssig2=ssig2, csig2=csig2)
if (outmask & Caps._DEBUG_LINE): # PYCHOK no cover r.set_(dn1=dn1, D0_k2=self._D0_k2, D1=self._D1, J12=J12, ssig1=ssig1, csig1=csig1, b=gX.b) # Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure # accurate cancellation in the case of coincident points. -dn1 * (ssig1 * csig2), -J12 * (csig1 * csig2))) M21=csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2)
# tan(alp) = tan(alp0) * sec(sig) # tan(alp2-alp1) = (tan(alp2) - tan(alp1)) / (tan(alp2) * tan(alp1) + 1) # = calp0 * salp0 * (csig1 - csig2) / (salp0^2 + calp0^2 * csig1 * csig2) # If csig12 > 0, write # csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1) # else # csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1 # No need to normalize (csig1 * (_1_0 - csig12) + ssig1 * ssig12)) * salp0 * calp0 else: # alp12 = alp2 - alp1, used in atan2 so no need to normalize # We used to include some patch up code that purported to deal # with nearly meridional geodesics properly. However, this turned # out to be wrong once salp1 = -0 was allowed (via InverseLine). # In fact, the calculation of {s,c}alp12 was already correct # (following the IEEE rules for handling signed zeros). So, # the patch up code was unnecessary (as well as dangerous). -A4_e2a2 * self._B41, gX.c2 * atan2(salp12, calp12))) if (outmask & Caps._DEBUG_LINE): # PYCHOK no cover r.set_(salp12=salp12, calp12=calp12, B42=B42, A4_e2a2=A4_e2a2, salp0=salp0, calp0=calp0, B41=self._B41, c2=gX.c2)
lat1=self.lat1, # == _fix90(lat1) lon1=self.lon1 if (outmask & Caps.LONG_UNROLL) else self._lon1_norm180, azi1=_norm180(self.azi1))
'''(INTERNAL) Generate a new position along the geodesic.
@return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, s12, m12, M12, M21, S12)}. '''
'''(INTERNAL) Aka C++ C{GenSetDistance}. ''' else:
'''Get the I{exact} geodesic (L{GeodesicExact}). ''' raise _AssertionError(gX=gX)
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) Cached/memoized. '''
'''Get the latitude of the first point (C{degrees}). '''
'''Get the longitude of the first point (C{degrees}). '''
'''(INTERNAL) Cached/memoized. '''
'''Find the position on the line given B{C{s12}}.
@arg s12: Distance from the first point to the second (C{meter}). @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying the quantities to be returned.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and arc length C{a12} always included, except when C{a12=NAN}.
@note: By default, C{B{outmask}=STANDARD}, meaning thc C{lat1}, C{lon1}, C{azi1}, C{lat2}, C{lon2}, C{azi2}, C{s12} and C{a12} entries are returned, except when C{a12=NAN}.
@note: This L{GeodesicLineExact} instance must have been constructed with capability C{Caps.DISTANCE_IN} set. '''
'''Get the distance from the first to the reference point (C{meter}).
@see: Method L{SetDistance}. '''
'''Set reference point 3 in terms of distance to the first point.
@arg a13: Spherical arc length from the first to the reference point (C{degrees}).
@return: The distance C{s13} (C{meter}) between the first and the reference point. ''' # In case the GeodesicLineExact doesn't have the DISTANCE capability.
'''Set reference point 3 in terms of distance to the first point.
@arg s13: Distance from the first to the reference point (C{meter}).
@return: The arc length C{a13} (C{degrees}) between the first and the reference point or C{NAN}. ''' # _a13 will be_NAN if the GeodesicLineExact # doesn't have the DISTANCE_IN capability
'''(INTERNAL) Cached/memoized. ''' # tau1 = sig1 + B11 # unnecessary because Einv inverts E # return -self._eF.deltaEinv(stau1, ctau1)
'''Return this C{GeodesicExactLine} as string.
@kwarg prec: The C{float} precision, number of decimal digits (0..9). Trailing zero decimals are stripped for B{C{prec}} values of 1 and above, but kept for negative B{C{prec}} values. @kwarg sep: Separator to join (C{str}).
@return: C{GeodesicExactLine} (C{str}). ''' lat1=self.lat1, lon1=self.lon1, azi1=self.azi1, a13=self.a13, s13=self.s13)
# **) MIT License # # Copyright (C) 2016-2022 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |