Coverage for pygeodesy/rhumbBase.py: 96%
317 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-20 13:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-09-20 13:43 -0400
2# -*- coding: utf-8 -*-
4u'''Base classes C{RhumbBase} and C{RhumbLineBase}, pure Python version of I{Karney}'s C++
5classes U{Rhumb<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Rhumb.html>}
6and U{RhumbLine<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1RhumbLine.html>}
7from I{GeographicLib versions 2.0} and I{2.2} and I{Karney}'s C++ example U{Rhumb intersect
8<https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}.
10Class L{RhumbLine} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively
11find the intersection of two rhumb lines, respectively the nearest point on a rumb line along a
12geodesic or perpendicular rhumb line from an other point.
14For more details, see the C++ U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>}
15documentation, especially the U{Class List<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>},
16the background information on U{Rhumb lines<https://GeographicLib.SourceForge.io/C++/doc/rhumb.html>},
17the utily U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online
18rhumb line calculations<https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve>}.
20Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2014-2023) and licensed under the MIT/X11
21License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
22'''
23# make sure int/int division yields float quotient
24from __future__ import division as _; del _ # PYCHOK semicolon
26# from pygeodesy.basics import unsigned0, _xinstanceof # from .karney
27from pygeodesy.constants import EPS, EPS1, INT0, _EPSqrt as _TOL, NAN, \
28 _0_0, _0_01, _1_0, _90_0, _over
29# from pygeodesy.datums import _spherical_datum # from .formy
30# from pygeodesy.ellipsoids import _EWGS84 # from .karney
31from pygeodesy.errors import IntersectionError, itemsorted, RhumbError, \
32 _xdatum, _xkwds, _Xorder
33# from pygeodesy.etm import ExactTransverseMercator # _MODS
34from pygeodesy.fmath import euclid, favg, fabs, Fsum
35from pygeodesy.formy import opposing, _spherical_datum
36# from pygeodesy.fsums import Fsum # from .fmath
37from pygeodesy.interns import NN, _coincident_, _COMMASPACE_, _Dash, \
38 _intersection_, _no_, _parallel_, _under
39from pygeodesy.karney import _atan2d, Caps, _CapsBase, _diff182, _fix90, \
40 GDict, _norm180, _EWGS84, unsigned0, _xinstanceof
41# from pygeodesy.ktm import KTransverseMercator, _AlpCoeffs # _MODS
42from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS
43# from pygeodesy.named import notOverloaded # _MODS
44from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple
45from pygeodesy.props import deprecated_method, Property, Property_RO, property_RO, \
46 _update_all
47from pygeodesy.streprs import Fmt, pairs, unstr
48from pygeodesy.units import Float_, Lat, Lon, Meter, Int # PYCHOK shared
49from pygeodesy.utily import atan2d, _azireversed, _loneg, sincos2d, sincos2d_, \
50 _unrollon, _Wrap, sincos2_ # PYCHOK shared
51from pygeodesy.vector3d import _intersect3d3, Vector3d # in .intersection2 below
53# from math import fabs # from .fmath
55__all__ = ()
56__version__ = '23.09.20'
58_anti_ = _Dash('anti')
59_rls = [] # instances of C{RbumbLine...} to be updated
60_TRIPS = 65 # .intersection2, .nearestOn4, 19+
63class _Lat(Lat):
64 '''(INTERNAL) Latitude B{C{lat}}.
65 '''
66 def __init__(self, *lat, **Error_name):
67 kwds = _xkwds(Error_name, clip=0, Error=RhumbError)
68 Lat.__new__(_Lat, *lat, **kwds)
71class _Lon(Lon):
72 '''(INTERNAL) Longitude B{C{lon}}.
73 '''
74 def __init__(self, *lon, **Error_name):
75 kwds = _xkwds(Error_name, clip=0, Error=RhumbError)
76 Lon.__new__(_Lon, *lon, **kwds)
79def _update_all_rls(r):
80 '''(INTERNAL) Zap cached/memoized C{Property[_RO]}s
81 of any C{RhumbLine} instances tied to the given
82 C{Rhumb} instance B{C{r}}.
83 '''
84 # _xinstanceof(_MODS.rhumbaux.RhumbAux, _MODS.rhumbx.Rhumb, r=r)
85 _update_all(r)
86 for rl in _rls: # PYCHOK use weakref?
87 if rl._rhumb is r:
88 _update_all(rl)
91class RhumbBase(_CapsBase):
92 '''(INTERNAL) Base class for C{rhumbaux.RhumbAux} and C{rhumbx.Rhumb}.
93 '''
94 _E = _EWGS84
95 _exact = True
96 _f_max = _0_01
97 _mTM = 6 # see .TMorder
99 def __init__(self, a_earth, f, exact, name):
100 '''New C{rhumbaux.RhumbAux} or C{rhumbx.Rhum}.
101 '''
102 if f is not None:
103 self.ellipsoid = a_earth, f
104 elif a_earth not in (_EWGS84, None):
105 self.ellipsoid = a_earth
106 if not exact:
107 self.exact = False
108 if name:
109 self.name = name
111 @Property_RO
112 def a(self):
113 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).
114 '''
115 return self.ellipsoid.a
117 equatoradius = a
119 def ArcDirect(self, lat1, lon1, azi12, a12, outmask=Caps.LATITUDE_LONGITUDE):
120 '''Solve the I{direct rhumb} problem, optionally with area.
122 @arg lat1: Latitude of the first point (C{degrees90}).
123 @arg lon1: Longitude of the first point (C{degrees180}).
124 @arg azi12: Azimuth of the rhumb line (compass C{degrees}).
125 @arg a12: Angle along the rhumb line from the given to the
126 destination point (C{degrees}), can be negative.
128 @return: L{GDict} with 2 up to 8 items C{lat2, lon2, a12, S12,
129 lat1, lon1, azi12, s12} with the destination point's
130 latitude C{lat2} and longitude C{lon2} in C{degrees},
131 the rhumb angle C{a12} in C{degrees} and area C{S12}
132 under the rhumb line in C{meter} I{squared}.
134 @raise ImportError: Package C{numpy} not found or not installed,
135 only required for area C{S12} when C{B{exact}
136 is True} and L{RhumbAux}.
138 @note: If B{C{a12}} is large enough that the rhumb line crosses
139 a pole, the longitude of the second point is indeterminate
140 and C{NAN} is returned for C{lon2} and area C{S12}.
142 @note: If the given point is a pole, the cosine of its latitude is
143 taken to be C{sqrt(L{EPS})}. This position is extremely
144 close to the actual pole and allows the calculation to be
145 carried out in finite terms.
146 '''
147 s12 = a12 * self._mpd
148 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask)
150 @Property_RO
151 def b(self):
152 '''Get the C{ellipsoid}'s polar radius, semi-axis (C{meter}).
153 '''
154 return self.ellipsoid.b
156 polaradius = b
158 def _Direct(self, ll1, azi12, s12, **outmask):
159 '''(INTERNAL) Short-cut version, see .latlonBase.
160 '''
161 return self.Direct(ll1.lat, ll1.lon, azi12, s12, **outmask)
163 def Direct(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE):
164 '''Solve the I{direct rhumb} problem, optionally with area.
166 @arg lat1: Latitude of the first point (C{degrees90}).
167 @arg lon1: Longitude of the first point (C{degrees180}).
168 @arg azi12: Azimuth of the rhumb line (compass C{degrees}).
169 @arg s12: Distance along the rhumb line from the given to
170 the destination point (C{meter}), can be negative.
172 @return: L{GDict} with 2 up to 8 items C{lat2, lon2, a12, S12,
173 lat1, lon1, azi12, s12} with the destination point's
174 latitude C{lat2} and longitude C{lon2} in C{degrees},
175 the rhumb angle C{a12} in C{degrees} and area C{S12}
176 under the rhumb line in C{meter} I{squared}.
178 @raise ImportError: Package C{numpy} not found or not installed,
179 only required for area C{S12} when C{B{exact}
180 is True} and L{RhumbAux}.
182 @note: If B{C{s12}} is large enough that the rhumb line crosses
183 a pole, the longitude of the second point is indeterminate
184 and C{NAN} is returned for C{lon2} and area C{S12}.
186 @note: If the given point is a pole, the cosine of its latitude is
187 taken to be C{sqrt(L{EPS})}. This position is extremely
188 close to the actual pole and allows the calculation to be
189 carried out in finite terms.
190 '''
191 a12 = _over(s12, self._mpd)
192 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask)
194 def Direct8(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA):
195 '''Like method L{Rhumb.Direct} but returning a L{Rhumb8Tuple} with area C{S12}.
196 '''
197 return self.Direct(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple()
199 def _DirectLine(self, ll1, azi12, **caps_name):
200 '''(INTERNAL) Short-cut version, see .latlonBase.
201 '''
202 return self.DirectLine(ll1.lat, ll1.lon, azi12, **caps_name)
204 def DirectLine(self, lat1, lon1, azi12, **caps_name):
205 '''Define a C{RhumbLine} in terms of the I{direct} rhumb
206 problem to compute several points on a single rhumb line.
208 @arg lat1: Latitude of the first point (C{degrees90}).
209 @arg lon1: Longitude of the first point (C{degrees180}).
210 @arg azi12: Azimuth of the rhumb line (compass C{degrees}).
211 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
212 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of
213 L{Caps} values specifying the required capabilities.
214 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}
215 should I{not} be reflected in this rhumb line.
217 @return: A C{RhumbLine...} instance and invoke its method
218 C{.Position} to compute each point.
220 @note: Updates to this rhumb are reflected in the returned
221 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.
222 '''
223 return self._RhumbLine(self, lat1, lon1, azi12, **caps_name)
225 Line = DirectLine # synonyms
227 def _DirectRhumb(self, lat1, lon1, azi12, a12, s12, outmask):
228 '''(INTERNAL) See methods C{.ArcDirect} and C{.Direct}.
229 '''
230 rl = self._RhumbLine(self, lat1, lon1, azi12, caps=Caps.LINE_OFF,
231 name=self.name)
232 return rl._Position(a12, s12, outmask | self._debug) # lat2, lon2, S12
234 @Property
235 def ellipsoid(self):
236 '''Get this rhumb's ellipsoid (L{Ellipsoid}).
237 '''
238 return self._E
240 @ellipsoid.setter # PYCHOK setter!
241 def ellipsoid(self, a_earth_f):
242 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or
243 L{a_f2Tuple}) or (equatorial) radius and flattening (2-tuple C{(a, f)}).
245 @raise RhumbError: If C{abs(B{f}} exceeds non-zero C{f_max} and C{exact=False}.
246 '''
247 E = _spherical_datum(a_earth_f, Error=RhumbError).ellipsoid
248 if self._E != E:
249 self._exactest(self.exact, E, self.f_max)
250 _update_all_rls(self)
251 self._E = E
253 @Property
254 def exact(self):
255 '''Get the I{exact} option (C{bool}).
256 '''
257 return self._exact
259 @exact.setter # PYCHOK setter!
260 def exact(self, exact):
261 '''Set the I{exact} option (C{bool}). If C{True}, use I{exact} rhumb
262 expressions, otherwise a series expansion (accurate for oblate or
263 prolate ellipsoids with C{abs(flattening)} below C{f_max}.
265 @raise RhumbError: If C{B{exact}=False} and C{abs(flattening})
266 exceeds non-zero C{f_max}.
268 @see: Option U{B{-s}<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>}
269 and U{ACCURACY<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html#ACCURACY>}.
270 '''
271 x = bool(exact)
272 if self._exact != x:
273 self._exactest(x, self.ellipsoid, self.f_max)
274 _update_all_rls(self)
275 self._exact = x
277 def _exactest(self, exact, ellipsoid, f_max):
278 # Helper for property setters C{ellipsoid}, C{exact} and C{f_max}
279 if fabs(ellipsoid.f) > f_max > 0 and not exact:
280 raise RhumbError(exact=exact, f=ellipsoid.f, f_max=f_max)
282 @Property_RO
283 def f(self):
284 '''Get the C{ellipsoid}'s flattening (C{float}).
285 '''
286 return self.ellipsoid.f
288 flattening = f
290 @property
291 def f_max(self):
292 '''Get the I{max.} flattening (C{float}).
293 '''
294 return self._f_max
296 @f_max.setter # PYCHOK setter!
297 def f_max(self, f_max): # PYCHOK no cover
298 '''Set the I{max.} flattening, not to exceed (C{float}).
300 @raise RhumbError: If C{exact=False} and C{abs(flattening})
301 exceeds non-zero C{f_max}.
302 '''
303 f = Float_(f_max=f_max, low=_0_0, high=EPS1)
304 if self._f_max != f:
305 self._exactest(self.exact, self.ellipsoid, f)
306 self._f_max = f
308 def _Inverse(self, ll1, ll2, wrap, **outmask):
309 '''(INTERNAL) Short-cut version, see .latlonBase.
310 '''
311 if wrap:
312 ll2 = _unrollon(ll1, _Wrap.point(ll2))
313 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask)
315 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH_DISTANCE):
316 '''Solve the I{inverse rhumb} problem.
318 @arg lat1: Latitude of the first point (C{degrees90}).
319 @arg lon1: Longitude of the first point (C{degrees180}).
320 @arg lat2: Latitude of the second point (C{degrees90}).
321 @arg lon2: Longitude of the second point (C{degrees180}).
323 @return: L{GDict} with 4 to 8 items C{lat1, lon1, lat2, lon2,
324 azi12, s12, a12, S12, }, the rhumb line's azimuth
325 C{azi12} in compass C{degrees} between C{-180} and
326 C{+180}, the rhumb distance C{s12} and rhumb angle
327 C{a12} between both points in C{meter} respectively
328 C{degrees} and the area C{S12} under the rhumb line
329 in C{meter} I{squared}.
331 @raise ImportError: Package C{numpy} not found or not installed,
332 only required for L{RhumbAux} area C{S12}
333 when C{B{exact} is True}.
335 @note: The shortest rhumb line is found. If the end points are
336 on opposite meridians, there are two shortest rhumb lines
337 and the East-going one is chosen.
339 @note: If either point is a pole, the cosine of its latitude is
340 taken to be C{sqrt(L{EPS})}. This position is extremely
341 close to the actual pole and allows the calculation to be
342 carried out in finite terms.
343 '''
344 r = GDict(lat1=lat1, lon1=lon1, lat2=lat2, lon2=lon2, name=self.name)
345 Cs = Caps
346 if (outmask & Cs.AZIMUTH_DISTANCE_AREA):
347 lon12, _ = _diff182(lon1, lon2, K_2_0=True)
348 y, x, s1, s2 = self._Inverse4(lon12, r, outmask)
349 if (outmask & Cs.AZIMUTH):
350 r.set_(azi12=_atan2d(y, x))
351 if (outmask & Cs.AREA):
352 S12 = self._S12d(s1, s2, lon12)
353 r.set_(S12=unsigned0(S12)) # like .gx
354 return r
356 def _Inverse4(self, lon12, r, outmask): # PYCHOK no cover
357 '''I{Must be overloaded}, see function C{notOverloaded}.
358 '''
359 _MODS.named.notOverloaded(self, lon12, r, outmask)
361 def Inverse8(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA):
362 '''Like method L{Rhumb.Inverse} but returning a L{Rhumb8Tuple} with area C{S12}.
363 '''
364 return self.Inverse(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple()
366 def _InverseLine(self, ll1, ll2, wrap, **caps_name):
367 '''(INTERNAL) Short-cut version, see .latlonBase.
368 '''
369 if wrap:
370 ll2 = _unrollon(ll1, _Wrap.point(ll2))
371 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name)
373 def InverseLine(self, lat1, lon1, lat2, lon2, **caps_name):
374 '''Define a C{RhumbLine} in terms of the I{inverse} rhumb problem.
376 @arg lat1: Latitude of the first point (C{degrees90}).
377 @arg lon1: Longitude of the first point (C{degrees180}).
378 @arg lat2: Latitude of the second point (C{degrees90}).
379 @arg lon2: Longitude of the second point (C{degrees180}).
380 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
381 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of
382 L{Caps} values specifying the required capabilities.
383 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}
384 should I{not} be reflected in this rhumb line.
386 @return: A C{RhumbLine...} instance and invoke its method
387 C{ArcPosition} or C{Position} to compute points.
389 @note: Updates to this rhumb are reflected in the returned
390 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.
391 '''
392 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH)
393 return self._RhumbLine(self, lat1, lon1, r.azi12, **caps_name)
395 @Property_RO
396 def _mpd(self): # PYCHOK no cover
397 '''I{Must be overloaded}, see function C{notOverloaded}.
398 '''
399 _MODS.named.notOverloaded(self)
401 @property_RO
402 def _RhumbLine(self): # PYCHOK no cover
403 '''I{Must be overloaded}, see function C{notOverloaded}.
404 '''
405 _MODS.named.notOverloaded(self)
407 def _S12d(self, s1, s2, lon): # PYCHOK no cover
408 '''I{Must be overloaded}, see function C{notOverloaded}.
409 '''
410 _MODS.named.notOverloaded(self, s1, s2, lon)
412 @Property
413 def TMorder(self):
414 '''Get the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
415 '''
416 return self._mTM
418 @TMorder.setter # PYCHOK setter!
419 def TMorder(self, order):
420 '''Set the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
422 @note: Setting C{TMorder} turns property C{exact} off, but only
423 for L{Rhumb} instances.
424 '''
425 m = _Xorder(_MODS.ktm._AlpCoeffs, RhumbError, TMorder=order)
426 if self._mTM != m:
427 _update_all_rls(self)
428 self._mTM = m
429 if self.exact and isinstance(self, _MODS.rhumbx.Rhumb):
430 self.exact = False
432 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK no cover
433 '''I{Must be overloaded}, see function C{notOverloaded}.
434 '''
435 _MODS.named.notOverloaded(self, prec=prec, sep=sep)
438class RhumbLineBase(_CapsBase):
439 '''(INTERNAL) Base class for C{rhumbaux.RhumbLineAux} and C{rhumbx.RhumbLine}.
440 '''
441 _azi12 = _0_0
442 _calp = _1_0
443# _caps = 0
444# _debug = 0
445# _lat1 = _0_0
446# _lon1 = _0_0
447# _lon12 = _0_0
448 _Rhumb = RhumbBase # compatible C{Rhumb} class
449 _rhumb = None # C{Rhumb} instance
450 _salp = _0_0
451 _talp = _0_0
453 def __init__(self, rhumb, lat1, lon1, azi12, caps=Caps.STANDARD, name=NN):
454 '''New C{RhumbLine}.
455 '''
456 _xinstanceof(self._Rhumb, rhumb=rhumb)
458 self._lat1 = _Lat(lat1=_fix90(lat1))
459 self._lon1 = _Lon(lon1= lon1)
460 self._lon12 = _norm180(self._lon1)
461 if azi12: # non-zero, non-None
462 self.azi12 = _norm180(azi12)
464 n = name or rhumb.name
465 if n:
466 self.name=n
468 self._caps = caps
469 self._debug |= (caps | rhumb._debug) & Caps._DEBUG_DIRECT_LINE
470 if (caps & Caps.LINE_OFF): # copy to avoid updates
471 self._rhumb = rhumb.copy(deep=False, name=_under(rhumb.name))
472 else:
473 self._rhumb = rhumb
474 _rls.append(self)
476 def __del__(self): # XXX use weakref?
477 if _rls: # may be empty or None
478 try: # PYCHOK no cover
479 _rls.remove(self)
480 except (TypeError, ValueError):
481 pass
482 self._rhumb = None
483 # _update_all(self) # throws TypeError during Python 2 cleanup
485 def ArcPosition(self, a12, outmask=Caps.LATITUDE_LONGITUDE):
486 '''Compute a point at a given angular distance on this rhumb line.
488 @arg a12: The angle along this rhumb line from its origin to the
489 point (C{degrees}), can be negative.
490 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
491 the quantities to be returned.
493 @return: L{GDict} with 4 to 8 items C{azi12, a12, s12, S12, lat2,
494 lon2, lat1, lon1} with latitude C{lat2} and longitude
495 C{lon2} of the point in C{degrees}, the rhumb distance
496 C{s12} in C{meter} from the start point of and the area
497 C{S12} under this rhumb line in C{meter} I{squared}.
499 @raise ImportError: Package C{numpy} not found or not installed,
500 only required for L{RhumbLineAux} area C{S12}
501 when C{B{exact} is True}.
503 @note: If B{C{a12}} is large enough that the rhumb line crosses a
504 pole, the longitude of the second point is indeterminate and
505 C{NAN} is returned for C{lon2} and area C{S12}.
507 If the first point is a pole, the cosine of its latitude is
508 taken to be C{sqrt(L{EPS})}. This position is extremely
509 close to the actual pole and allows the calculation to be
510 carried out in finite terms.
511 '''
512 return self._Position(a12, self.degrees2m(a12), outmask)
514 @Property
515 def azi12(self):
516 '''Get this rhumb line's I{azimuth} (compass C{degrees}).
517 '''
518 return self._azi12
520 @azi12.setter # PYCHOK setter!
521 def azi12(self, azi12):
522 '''Set this rhumb line's I{azimuth} (compass C{degrees}).
523 '''
524 z = _norm180(azi12)
525 if self._azi12 != z:
526 if self._rhumb:
527 _update_all(self)
528 self._azi12 = z
529 self._salp, self._calp = t = sincos2d(z) # no NEG0
530 self._talp = _over(*t)
532 @property_RO
533 def azi12_sincos2(self): # PYCHOK no cover
534 '''Get the sine and cosine of this rhumb line's I{azimuth} (2-tuple C{(sin, cos)}).
535 '''
536 return self._scalp, self._calp
538 def degrees2m(self, angle):
539 '''Convert an angular distance along this rhumb line to C{meter}.
541 @arg angle: Angular distance (C{degrees}).
543 @return: Distance (C{meter}).
544 '''
545 return float(angle) * self.rhumb._mpd
547 @deprecated_method
548 def distance2(self, lat, lon): # PYCHOK no cover
549 '''DEPRECATED, use method L{RhumbLineAux.Inverse} or L{RhumbLine.Inverse}.
551 @return: A L{Distance2Tuple}C{(distance, initial)} with the C{distance}
552 in C{meter} and C{initial} bearing (azimuth) in C{degrees}.
553 '''
554 r = self.Inverse(lat, lon)
555 return Distance2Tuple(r.s12, r.azi12)
557 @property_RO
558 def ellipsoid(self):
559 '''Get this rhumb line's ellipsoid (L{Ellipsoid}).
560 '''
561 return self.rhumb.ellipsoid
563 @property_RO
564 def exact(self):
565 '''Get this rhumb line's I{exact} option (C{bool}).
566 '''
567 return self.rhumb.exact
569 def intersection2(self, other, tol=_TOL, **eps):
570 '''I{Iteratively} find the intersection of this and an other rhumb line.
572 @arg other: The other rhumb line (C{RhumbLine}).
573 @kwarg tol: Tolerance for longitudinal convergence and parallel
574 error (C{degrees}).
575 @kwarg eps: Tolerance for L{pygeodesy.intersection3d3} (C{EPS}).
577 @return: A L{LatLon2Tuple}{(lat, lon)} with the C{lat}- and
578 C{lon}gitude of the intersection point.
580 @raise IntersectionError: No convergence for this B{C{tol}} or
581 no intersection for an other reason.
583 @see: Methods C{distance2} and C{nearestOn4} and function
584 L{pygeodesy.intersection3d3}.
586 @note: Each iteration involves a round trip to this rhumb line's
587 L{ExactTransverseMercator} or L{KTransverseMercator}
588 projection and invoking function L{pygeodesy.intersection3d3}
589 in that domain.
590 '''
591 _xinstanceof(RhumbLineBase, other=other)
592 _xdatum(self.rhumb, other.rhumb, Error=RhumbError)
593 try:
594 if other is self:
595 raise ValueError(_coincident_)
596 # make invariants and globals locals
597 _s_3d, s_az = self._xTM3d, self.azi12
598 _o_3d, o_az = other._xTM3d, other.azi12
599 p = opposing(s_az, o_az, margin=tol)
600 if p is not None: # == t in (True, False)
601 raise ValueError(_anti_(_parallel_) if p else _parallel_)
602 _diff = euclid # approximate length
603 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3
604 _LL2T = LatLon2Tuple
605 _xTMr = self.xTM.reverse # ellipsoidal or spherical
606 # use halfway point as initial estimate
607 p = _LL2T(favg(self.lat1, other.lat1),
608 favg(self.lon1, other.lon1))
609 for i in range(1, _TRIPS):
610 v = _i3d3(_s_3d(p), s_az, # point + bearing
611 _o_3d(p), o_az, useZ=False, **eps)[0]
612 t = _xTMr(v.x, v.y, lon0=p.lon) # PYCHOK Reverse4Tuple
613 d = _diff(t.lon - p.lon, t.lat) # PYCHOK t.lat + p.lat - p.lat
614 p = _LL2T(t.lat + p.lat, t.lon) # PYCHOK t.lon + p.lon = lon0
615 if d < tol: # 19 trips
616 return _LL2T(p.lat, p.lon, iteration=i, # PYCHOK p...
617 name=self.intersection2.__name__)
618 except Exception as x:
619 raise IntersectionError(_no_(_intersection_), cause=x)
620 t = unstr(self.intersection2, other, **eps)
621 raise IntersectionError(Fmt.no_convergence(d, tol), txt=t)
623 def Inverse(self, lat2, lon2, wrap=False):
624 '''Return the rhumb angle, distance, azimuth, I{reverse} azimuth, etc. of
625 a rhumb line between the given point and this rhumb line's start point.
627 @arg lat2: Latitude of the point (C{degrees}).
628 @arg lon2: Longitude of the points (C{degrees}).
629 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
630 and B{C{lon2}} (C{bool}).
632 @return: L{GDict} with 8 items C{a12, s12, azi12, azi21, lat1, lon1,
633 lat2, lon2}, the rhumb angle C{a12} and rhumb distance C{s12}
634 between both points in C{degrees} respectively C{meter}, the
635 rhumb line's azimuth C{azi12} and I{reverse} azimuth C{azi21}
636 both in compass C{degrees} between C{-180} and C{+180}.
637 '''
638 if wrap:
639 _, lat2, lon2 = _Wrap.latlon3(self.lon1, _fix90(lat2), lon2)
640 r = self.rhumb.Inverse(self.lat1, self.lon1, lat2, lon2)
641 r.set_(azi21=_azireversed(r.azi12))
642 return r
644 @Property_RO
645 def isLoxodrome(self):
646 '''Is this rhumb line a loxodrome (C{True} if non-meridional and
647 non-parallel, C{False} if parallel or C{None} if meridional)?
649 @see: I{Osborne's} U{2.5 Rumb lines and loxodromes
650 <https://Zenodo.org/record/35392>}, page 37.
651 '''
652 return bool(self._salp) if self._calp else None
654 @Property_RO
655 def lat1(self):
656 '''Get this rhumb line's latitude (C{degrees90}).
657 '''
658 return self._lat1
660 @Property_RO
661 def lon1(self):
662 '''Get this rhumb line's longitude (C{degrees180}).
663 '''
664 return self._lon1
666 @Property_RO
667 def latlon1(self):
668 '''Get this rhumb line's lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).
669 '''
670 return LatLon2Tuple(self.lat1, self.lon1)
672 def m2degrees(self, distance):
673 '''Convert a distance along this rhumb line to an angular distance.
675 @arg distance: Distance (C{meter}).
677 @return: Angular distance (C{degrees}).
678 '''
679 return _over(float(distance), self.rhumb._mpd)
681 @property_RO
682 def _mu1(self): # PYCHOK no cover
683 '''I{Must be overloaded}, see function C{notOverloaded}.
684 '''
685 _MODS.named.notOverloaded(self)
687 def _mu2lat(self, mu2): # PYCHOK no cover
688 '''I{Must be overloaded}, see function C{notOverloaded}.
689 '''
690 _MODS.named.notOverloaded(self, mu2)
692 def nearestOn4(self, lat, lon, tol=_TOL, exact=None, eps=EPS, est=None):
693 '''I{Iteratively} locate the point on this rhumb line nearest to the
694 given point, in part transcoded from I{Karney}'s C++ U{rhumb-intercept
695 <https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}.
697 @arg lat: Latitude of the point (C{degrees}).
698 @arg lon: Longitude of the point (C{degrees}).
699 @kwarg tol: Longitudinal convergence tolerance (C{degrees}) or the
700 distance tolerance (C(meter)) when C{B{exact} is None},
701 respectively C{not None}.
702 @kwarg exact: If C{None}, use a rhumb line perpendicular to this rhumb
703 line, otherwise use an I{exact} C{Geodesic...} from the
704 given point perpendicular to this rhumb line (C{bool} or
705 C{Geodesic...}), see method L{Ellipsoid.geodesic_}.
706 @kwarg eps: Optional tolerance for L{pygeodesy.intersection3d3}
707 (C{EPS}), used only if C{B{exact} is None}.
708 @kwarg est: Optional, initial estimate for the distance C{s13} of
709 the intersection I{along} this rhumb line (C{meter}),
710 used only if C{B{exact} is not None}.
712 @return: A L{NearestOn4Tuple}C{(lat, lon, distance, normal)} with
713 the C{lat}- and C{lon}gitude of the nearest point on and
714 the C{distance} in C{meter} to this rhumb line and the
715 C{normal} azimuth at the intersection.
717 @raise ImportError: I{Karney}'s U{geographiclib
718 <https://PyPI.org/project/geographiclib>}
719 package not found or not installed.
721 @raise IntersectionError: No convergence for this B{C{eps}} or
722 no intersection for an other reason.
724 @see: Methods C{distance2} and C{intersection2} and function
725 L{pygeodesy.intersection3d3}.
726 '''
727 if exact is None:
728 z = _norm180(self.azi12 + _90_0) # perpendicular azimuth
729 rl = RhumbLineBase(self.rhumb, lat, lon, z, caps=Caps.LINE_OFF)
730 p = self.intersection2(rl, tol=tol, eps=eps)
731 r = rl.Inverse(p.lat, p.lon)
732 r = NearestOn4Tuple(p.lat, p.lon, r.s12, z,
733 iteration=p.iteration)
734 else: # C{rhumb-intercept}
735 azi = self.azi12
736 szi = self._salp
737 E = self.ellipsoid
738 _gI = E.geodesic_(exact=exact).Inverse
739 Cs = Caps
740 gm = Cs._AZIMUTH_LATITUDE_LONGITUDE | Cs._REDUCEDLENGTH_GEODESICSCALE
741 if est is None: # get an estimate from the perpendicular
742 r = _gI(self.lat1, self.lon1, lat, lon, outmask=Cs.AZIMUTH_DISTANCE)
743 d, _ = _diff182(r.azi2, azi, K_2_0=True)
744 _, c = sincos2d(d)
745 s12 = c * r.s12 # signed
746 else:
747 s12 = Meter(est=est)
748 try:
749 tol = Float_(tol=tol, low=EPS, high=None)
750 # def _over(p, q): # see @note at RhumbLine[Aux].Position
751 # return (p / (q or _copysign(tol, q))) if isfinite(q) else NAN
753 _S12 = Fsum(s12).fsum2_
754 for i in range(1, _TRIPS): # suffix 1 == C++ 2, 2 == C++ 3
755 p = self.Position(s12) # outmask=Cs.LATITUDE_LONGITUDE
756 r = _gI(lat, lon, p.lat2, p.lon2, outmask=gm)
757 d, _ = _diff182(azi, r.azi2, K_2_0=True)
758 s, c, s2, c2 = sincos2d_(d, r.lat2)
759 c2 *= E.rocPrimeVertical(r.lat2) # aka rocTransverse
760 s *= _over(s2 * szi, c2) - _over(s * r.M21, r.m12)
761 s12, t = _S12(c / s) # XXX _over?
762 if fabs(t) < tol: # or fabs(c) < EPS
763 break
764 r = NearestOn4Tuple(r.lat2, r.lon2, s12, r.azi2,
765 iteration=i)
766 except Exception as x: # Fsum(NAN) Value-, ZeroDivisionError
767 t = unstr(self.nearestOn4, lat, lon, tol=tol, exact=exact,
768 iteration=i, eps=eps, est=est)
769 raise IntersectionError(t, cause=x)
770 return r
772 def Position(self, s12, outmask=Caps.LATITUDE_LONGITUDE):
773 '''Compute a point at a given distance on this rhumb line.
775 @arg s12: The distance along this rhumb line from its origin to
776 the point (C{meters}), can be negative.
777 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying
778 the quantities to be returned.
780 @return: L{GDict} with 4 to 8 items C{azi12, a12, s12, S12, lat2,
781 lon2, lat1, lon1} with latitude C{lat2} and longitude
782 C{lon2} of the point in C{degrees}, the rhumb angle C{a12}
783 in C{degrees} from the start point of and the area C{S12}
784 under this rhumb line in C{meter} I{squared}.
786 @raise ImportError: Package C{numpy} not found or not installed,
787 only required for L{RhumbLineAux} area C{S12}
788 when C{B{exact} is True}.
790 @note: If B{C{s12}} is large enough that the rhumb line crosses a
791 pole, the longitude of the second point is indeterminate and
792 C{NAN} is returned for C{lon2} and area C{S12}.
794 If the first point is a pole, the cosine of its latitude is
795 taken to be C{sqrt(L{EPS})}. This position is extremely
796 close to the actual pole and allows the calculation to be
797 carried out in finite terms.
798 '''
799 return self._Position(self.m2degrees(s12), s12, outmask)
801 def _Position(self, a12, s12, outmask):
802 '''(INTERNAL) C{Arc-/Position} helper.
803 '''
804 r = GDict(azi12=self.azi12, a12=a12, s12=s12, name=self.name)
805 Cs = Caps
806 if (outmask & Cs.LATITUDE_LONGITUDE_AREA):
807 if a12 or s12:
808 mu12 = self._calp * a12
809 mu2 = self._mu1 + mu12
810 if fabs(mu2) > 90: # past pole
811 mu2 = _norm180(mu2) # reduce to [-180, 180)
812 if fabs(mu2) > 90: # point on anti-meridian
813 mu2 = _norm180(_loneg(mu2))
814 lat2 = self._mu2lat(mu2)
815 lon2 = S12 = NAN
816 else:
817 lat2, lon2, S1, S2 = self._Position4(a12, mu2, s12, mu12)
818 if (outmask & Cs.AREA):
819 S12 = self.rhumb._S12d(S1, S2, lon2)
820 S12 = unsigned0(S12) # like .gx
821# else:
822# S12 = None # unused
823 if (outmask & Cs.LONGITUDE):
824 if (outmask & Cs.LONG_UNROLL):
825 lon2 += self.lon1
826 else:
827 lon2 = _norm180(self._lon12 + lon2)
828 else: # coincident or outmask 0
829 lat2, lon2 = self.latlon1
830 S12 = _0_0
832 if (outmask & Cs.LATITUDE):
833 r.set_(lat2=lat2, lat1=self.lat1)
834 if (outmask & Cs.LONGITUDE):
835 r.set_(lon2=lon2, lon1=self.lon1)
836 if (outmask & Cs.AREA):
837 r.set_(S12=S12)
838 return r
840 def _Position4(self, a12, mu2, s12, mu12): # PYCHOK no cover
841 '''I{Must be overloaded}, see function C{notOverloaded}.
842 '''
843 _MODS.named.notOverloaded(self, a12, s12, mu2, mu12)
845 @Property_RO
846 def rhumb(self):
847 '''Get this rhumb line's rhumb (L{RhumbAux} or L{Rhumb}).
848 '''
849 return self._rhumb
851 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
852 '''Return this C{RhumbLine} as string.
854 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
855 Trailing zero decimals are stripped for B{C{prec}} values
856 of 1 and above, but kept for negative B{C{prec}} values.
857 @kwarg sep: Separator to join (C{str}).
859 @return: C{RhumbLine} (C{str}).
860 '''
861 d = dict(rhumb=self.rhumb, lat1=self.lat1, lon1=self.lon1,
862 azi12=self.azi12, exact=self.exact,
863 TMorder=self.TMorder, xTM=self.xTM)
864 return sep.join(pairs(itemsorted(d, asorted=False), prec=prec))
866 @property_RO
867 def TMorder(self):
868 '''Get this rhumb line's I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
869 '''
870 return self.rhumb.TMorder
872 @Property_RO
873 def xTM(self):
874 '''Get this rhumb line's I{Transverse Mercator} projection (L{ExactTransverseMercator}
875 if I{exact} and I{ellipsoidal}, otherwise L{KTransverseMercator} for C{TMorder}).
876 '''
877 E = self.ellipsoid
878 # ExactTransverseMercator doesn't handle spherical earth models
879 return _MODS.etm.ExactTransverseMercator(E) if self.exact and E.isEllipsoidal else \
880 _MODS.ktm.KTransverseMercator(E, TMorder=self.TMorder)
882 def _xTM3d(self, latlon0, z=INT0, V3d=Vector3d):
883 '''(INTERNAL) C{xTM.forward} this C{latlon1} to C{V3d} with B{C{latlon0}}
884 as current intersection estimate and central meridian.
885 '''
886 t = self.xTM.forward(self.lat1 - latlon0.lat, self.lon1, lon0=latlon0.lon)
887 return V3d(t.easting, t.northing, z)
890__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase)
892if __name__ == '__main__':
894 from pygeodesy import printf, Rhumb as R, RhumbAux as A
896 A = A(_EWGS84).Line(30, 0, 45)
897 R = R(_EWGS84).Line(30, 0, 45)
899 for i in range(1, 10):
900 s = .5e6 + 1e6 / i
901 a = A.Position(s).lon2
902 r = R.Position(s).lon2
903 e = (fabs(a - r) / a) if a else 0
904 printf('Position.lon2 %.14f vs %.14f, diff %g', r, a, e)
906 for exact in (None, False, True):
907 for est in (None, 1e6):
908 a = A.nearestOn4(60, 0, exact=exact, est=est)
909 r = R.nearestOn4(60, 0, exact=exact, est=est)
910 printf('%s, iteration=%s, exact=%s, est=%s\n%s, iteration=%s',
911 a.toRepr(), a.iteration, exact, est,
912 r.toRepr(), r.iteration, nl=1)
914# % python3 -m pygeodesy.rhumbBase
916# Position.lon2 11.61455846901637 vs 11.61455846901637, diff 1.52942e-16
917# Position.lon2 7.58982302826842 vs 7.58982302826842, diff 2.34045e-16
918# Position.lon2 6.28526067416369 vs 6.28526067416369, diff 2.82623e-16
919# Position.lon2 5.63938995325146 vs 5.63938995325146, diff 1.57495e-16
920# Position.lon2 5.25385527435707 vs 5.25385527435707, diff 0
921# Position.lon2 4.99764604290380 vs 4.99764604290380, diff 8.88597e-16
922# Position.lon2 4.81503363740473 vs 4.81503363740473, diff 1.84459e-16
923# Position.lon2 4.67828821748836 vs 4.67828821748835, diff 5.69553e-16
924# Position.lon2 4.57205667906283 vs 4.57205667906283, diff 1.94262e-16
926# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=None
927# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=3278714.911694, normal=135.0), iteration=9
929# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=1000000.0
930# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=3278714.911694, normal=135.0), iteration=9
932# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=False, est=None
933# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5
935# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=False, est=1000000.0
936# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7
938# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=True, est=None
939# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5
941# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=True, est=1000000.0
942# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7
944# **) MIT License
945#
946# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
947#
948# Permission is hereby granted, free of charge, to any person obtaining a
949# copy of this software and associated documentation files (the "Software"),
950# to deal in the Software without restriction, including without limitation
951# the rights to use, copy, modify, merge, publish, distribute, sublicense,
952# and/or sell copies of the Software, and to permit persons to whom the
953# Software is furnished to do so, subject to the following conditions:
954#
955# The above copyright notice and this permission notice shall be included
956# in all copies or substantial portions of the Software.
957#
958# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
959# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
960# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
961# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
962# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
963# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
964# OTHER DEALINGS IN THE SOFTWARE.