Coverage for pygeodesy/rhumbBase.py: 96%
254 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'''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.
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 _xinstanceof # from .karney
27from pygeodesy.constants import EPS, EPS1, INT0, _EPSqrt as _TOL, \
28 _0_0, _0_01, _1_0, _90_0, _180_0
29# from pygeodesy.datums import _spherical_datum # _MODS
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
35# from pygeodesy.fsums import Fsum # _MODS
36from pygeodesy.interns import NN, _coincident_, _COMMASPACE_, _intersection_, \
37 _no_, _parallel_, _under
38from pygeodesy.karney import Caps, _CapsBase, _diff182, _fix90, _norm180, \
39 _EWGS84, _xinstanceof
40from pygeodesy.ktm import KTransverseMercator, _AlpCoeffs # PYCHOK used!
41from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS
42# from pygeodesy.named import notOverloaded # _MODS
43from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple
44from pygeodesy.props import Property, Property_RO, property_RO, _update_all
45from pygeodesy.streprs import Fmt, pairs, unstr
46from pygeodesy.units import Float_, Lat, Lon, Meter
47from pygeodesy.utily import sincos2d, sincos2d_, _unrollon, _Wrap
48from pygeodesy.vector3d import _intersect3d3, Vector3d # in .intersection2 below
50# from math import fabs # from .fmath
52__all__ = ()
53__version__ = '23.08.25'
55_rls = [] # instances of C{RbumbLine...} to be updated
56_TRIPS = 65 # .intersection2, .nearestOn4, 18+
59class _Lat(Lat):
60 '''(INTERNAL) Latitude B{C{lat}}.
61 '''
62 def __init__(self, *lat, **Error_name):
63 kwds = _xkwds(Error_name, clip=0, Error=RhumbError)
64 Lat.__new__(_Lat, *lat, **kwds)
67class _Lon(Lon):
68 '''(INTERNAL) Longitude B{C{lon}}.
69 '''
70 def __init__(self, *lon, **Error_name):
71 kwds = _xkwds(Error_name, clip=0, Error=RhumbError)
72 Lon.__new__(_Lon, *lon, **kwds)
75def _update_all_rls(r):
76 '''(INTERNAL) Zap cached/memoized C{Property[_RO]}s
77 of any C{RhumbLine} instances tied to the given
78 C{Rhumb} instance B{C{r}}.
79 '''
80 # _xinstanceof(_MODS.rhumbaux.RhumbAux, _MODS.rhumbx.Rhumb, r=r)
81 _update_all(r)
82 for rl in _rls: # PYCHOK use weakref?
83 if rl._rhumb is r:
84 _update_all(rl)
87class RhumbBase(_CapsBase):
88 '''(INTERNAL) Base class for C{rhumbaux.RhumbAux} and C{rhumbx.Rhumb}.
89 '''
90 _E = _EWGS84
91 _exact = True
92 _f_max = _0_01
93 _mTM = 6 # see .TMorder
95 def __init__(self, a_earth, f, exact, name):
96 '''New C{rhumbaux.RhumbAux} or C{rhumbx.Rhum}.
97 '''
98 if f is not None:
99 self.ellipsoid = a_earth, f
100 elif a_earth not in (_EWGS84, None):
101 self.ellipsoid = a_earth
102 if not exact:
103 self.exact = False
104 if name:
105 self.name = name
107 @Property_RO
108 def a(self):
109 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).
110 '''
111 return self.ellipsoid.a
113 equatoradius = a
115 @Property_RO
116 def b(self):
117 '''Get the C{ellipsoid}'s polar radius, semi-axis (C{meter}).
118 '''
119 return self.ellipsoid.b
121 polaradius = b
123 def _Direct(self, ll1, azi12, s12, **outmask):
124 '''(INTERNAL) Short-cut version, see .latlonBase.
125 '''
126 return self.Direct(ll1.lat, ll1.lon, azi12, s12, **outmask)
128 def Direct(self, lat1, lon1, azi12, s12, outmask=0): # PYCHOK no cover
129 '''I{Must be overloaded}, see function C{notOverloaded}.
130 '''
131 _MODS.named.notOverloaded(self, lat1, lon1, azi12, s12, outmask=outmask)
133 def Direct8(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA):
134 '''Like method L{Rhumb.Direct} but returning a L{Rhumb8Tuple} with area C{S12}.
135 '''
136 return self.Direct(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple()
138 def _DirectLine(self, ll1, azi12, **caps_name):
139 '''(INTERNAL) Short-cut version, see .latlonBase.
140 '''
141 return self.DirectLine(ll1.lat, ll1.lon, azi12, **caps_name)
143 def DirectLine(self, lat1, lon1, azi12, **caps_name):
144 '''Define a C{RhumbLine} in terms of the I{direct} rhumb
145 problem to compute several points on a single rhumb line.
147 @arg lat1: Latitude of the first point (C{degrees90}).
148 @arg lon1: Longitude of the first point (C{degrees180}).
149 @arg azi12: Azimuth of the rhumb line (compass C{degrees}).
150 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
151 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of
152 L{Caps} values specifying the required capabilities.
153 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}
154 should I{not} be reflected in this rhumb line.
156 @return: A C{RhumbLine...} instance and invoke its method
157 C{.Position} to compute each point.
159 @note: Updates to this rhumb are reflected in the returned
160 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.
161 '''
162 return self._RhumbLine(self, lat1=lat1, lon1=lon1, azi12=azi12,
163 **caps_name)
165 Line = DirectLine # synonyms
167 @Property
168 def ellipsoid(self):
169 '''Get this rhumb's ellipsoid (L{Ellipsoid}).
170 '''
171 return self._E
173 @ellipsoid.setter # PYCHOK setter!
174 def ellipsoid(self, a_earth_f):
175 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or
176 L{a_f2Tuple}) or (equatorial) radius and flattening (2-tuple C{(a, f)}).
178 @raise RhumbError: If C{abs(B{f}} exceeds non-zero C{f_max} and C{exact=False}.
179 '''
180 E = _MODS.datums._spherical_datum(a_earth_f, Error=RhumbError).ellipsoid
181 if self._E != E:
182 self._exactest(self.exact, E, self.f_max)
183 _update_all_rls(self)
184 self._E = E
186 @Property
187 def exact(self):
188 '''Get the I{exact} option (C{bool}).
189 '''
190 return self._exact
192 @exact.setter # PYCHOK setter!
193 def exact(self, exact):
194 '''Set the I{exact} option (C{bool}). If C{True}, use I{exact} rhumb
195 expressions, otherwise a series expansion (accurate for oblate or
196 prolate ellipsoids with C{abs(flattening)} below C{f_max}.
198 @raise RhumbError: If C{B{exact}=False} and C{abs(flattening})
199 exceeds non-zero C{f_max}.
201 @see: Option U{B{-s}<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>}
202 and U{ACCURACY<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html#ACCURACY>}.
203 '''
204 x = bool(exact)
205 if self._exact != x:
206 self._exactest(x, self.ellipsoid, self.f_max)
207 _update_all_rls(self)
208 self._exact = x
210 def _exactest(self, exact, ellipsoid, f_max):
211 # Helper for property setters C{ellipsoid}, C{exact} and C{f_max}
212 if fabs(ellipsoid.f) > f_max > 0 and not exact:
213 raise RhumbError(exact=exact, f=ellipsoid.f, f_max=f_max)
215 @Property_RO
216 def f(self):
217 '''Get the C{ellipsoid}'s flattening (C{float}).
218 '''
219 return self.ellipsoid.f
221 flattening = f
223 @property
224 def f_max(self):
225 '''Get the I{max.} flattening (C{float}).
226 '''
227 return self._f_max
229 @f_max.setter # PYCHOK setter!
230 def f_max(self, f_max): # PYCHOK no cover
231 '''Set the I{max.} flattening, not to exceed (C{float}).
233 @raise RhumbError: If C{exact=False} and C{abs(flattening})
234 exceeds non-zero C{f_max}.
235 '''
236 f = Float_(f_max=f_max, low=_0_0, high=EPS1)
237 if self._f_max != f:
238 self._exactest(self.exact, self.ellipsoid, f)
239 self._f_max = f
241 def _Inverse(self, ll1, ll2, wrap, **outmask):
242 '''(INTERNAL) Short-cut version, see .latlonBase.
243 '''
244 if wrap:
245 ll2 = _unrollon(ll1, _Wrap.point(ll2))
246 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask)
248 def Inverse(self, lat1, lon1, lat2, lon2, outmask=0): # PYCHOK no cover
249 '''I{Must be overloaded}, see function C{notOverloaded}.
250 '''
251 _MODS.named.notOverloaded(self, lat1, lon1, lat2, lon2, outmask=outmask)
253 def Inverse8(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA):
254 '''Like method L{Rhumb.Inverse} but returning a L{Rhumb8Tuple} with area C{S12}.
255 '''
256 return self.Inverse(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple()
258 def _InverseLine(self, ll1, ll2, wrap, **caps_name):
259 '''(INTERNAL) Short-cut version, see .latlonBase.
260 '''
261 if wrap:
262 ll2 = _unrollon(ll1, _Wrap.point(ll2))
263 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name)
265 def InverseLine(self, lat1, lon1, lat2, lon2, **caps_name):
266 '''Define a C{RhumbLine} in terms of the I{inverse} rhumb problem.
268 @arg lat1: Latitude of the first point (C{degrees90}).
269 @arg lon1: Longitude of the first point (C{degrees180}).
270 @arg lat2: Latitude of the second point (C{degrees90}).
271 @arg lon2: Longitude of the second point (C{degrees180}).
272 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
273 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of
274 L{Caps} values specifying the required capabilities.
275 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}
276 should I{not} be reflected in this rhumb line.
278 @return: A C{RhumbLine...} instance and invoke its method
279 C{.Position} to compute each point.
281 @note: Updates to this rhumb are reflected in the returned
282 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.
283 '''
284 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH)
285 return self._RhumbLine(self, lat1=lat1, lon1=lon1, azi12=r.azi12,
286 **caps_name)
288 @Property_RO
289 def _RhumbLine(self): # PYCHOK no cover
290 '''I{Must be overloaded}, see function C{notOverloaded}.
291 '''
292 _MODS.named.notOverloaded(self)
294 def _TMorder(self, order):
295 '''(INTERNAL) Set the I{Transverse Mercator} order (C{int}).
297 @note: Setting C{TMorder} turns property C{exact} off.
298 '''
299 x = self.exact
300 n = _Xorder(_MODS.ktm._AlpCoeffs, RhumbError, TMorder=order)
301 if self._mTM != n:
302 _update_all_rls(self)
303 self._mTM = n
304 x = False
305 return x
307 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK no cover
308 '''I{Must be overloaded}, see function C{notOverloaded}.
309 '''
310 _MODS.named.notOverloaded(self, prec=prec, sep=sep)
313class RhumbLineBase(RhumbBase):
314 '''(INTERNAL) Class C{RhumbLine}
315 '''
316 _azi12 = _0_0
317# _caps = 0
318# _debug = 0
319# _lat1 = _0_0
320# _lon1 = _0_0
321# _lon12 = _0_0
322 _salp = _0_0
323 _calp = _1_0
324 _Rhumb = RhumbBase # compatible C{Rhumb} class
325 _rhumb = None # C{Rhumb} instance
327 def __init__(self, rhumb, lat1, lon1, azi12, caps=Caps.STANDARD, name=NN):
328 '''New C{RhumbLine}.
329 '''
330 _xinstanceof(self._Rhumb, rhumb=rhumb)
332 self._lat1 = _Lat(lat1=_fix90(lat1))
333 self._lon1 = _Lon(lon1= lon1)
334 self._lon12 = _norm180(self._lon1)
335 if azi12: # non-zero, non-None
336 self.azi12 = _norm180(azi12)
338 n = name or rhumb.name
339 if n:
340 self.name=n
342 self._caps = caps
343 self._debug |= (caps | rhumb._debug) & Caps._DEBUG_DIRECT_LINE
344 if (caps & Caps.LINE_OFF): # copy to avoid updates
345 self._rhumb = rhumb.copy(deep=False, name=_under(rhumb.name))
346 else:
347 self._rhumb = rhumb
348 _rls.append(self)
350 def __del__(self): # XXX use weakref?
351 if _rls: # may be empty or None
352 try: # PYCHOK no cover
353 _rls.remove(self)
354 except (TypeError, ValueError):
355 pass
356 self._rhumb = None
357 # _update_all(self) # throws TypeError during Python 2 cleanup
359 @Property
360 def azi12(self):
361 '''Get this rhumb line's I{azimuth} (compass C{degrees}).
362 '''
363 return self._azi12
365 @azi12.setter # PYCHOK setter!
366 def azi12(self, azi12):
367 '''Set this rhumb line's I{azimuth} (compass C{degrees}).
368 '''
369 z = _norm180(azi12)
370 if self._azi12 != z:
371 if self._rhumb:
372 _update_all(self)
373 self._azi12 = z
374 self._salp, self._calp = sincos2d(z) # no NEG0
376 @property_RO
377 def azi12_sincos2(self): # PYCHOK no cover
378 '''Get this rhumb line's I{azimuth} sine and cosine (2-tuple C{(sin, cos)}).
379 '''
380 return self._scalp, self._calp
382 def distance2(self, lat, lon):
383 '''Return the distance from and (initial) bearing at the given
384 point to this rhumb line's start point.
386 @arg lat: Latitude of the point (C{degrees}).
387 @arg lon: Longitude of the points (C{degrees}).
389 @return: A L{Distance2Tuple}C{(distance, initial)} with the C{distance}
390 in C{meter} and C{initial} bearing in C{degrees}.
392 @see: Methods C{intersection2} and C{nearestOn4} of L{RhumbLine} and
393 L{RhumbLineAux}.
394 '''
395 r = self.rhumb.Inverse(self.lat1, self.lon1, lat, lon)
396# outmask=Caps.AZIMUTH_DISTANCE)
397 return Distance2Tuple(r.s12, r.azi12)
399 @property_RO
400 def ellipsoid(self):
401 '''Get this rhumb line's ellipsoid (L{Ellipsoid}).
402 '''
403 return self.rhumb.ellipsoid
405 @property_RO
406 def exact(self):
407 '''Get this rhumb line's I{exact} option (C{bool}).
408 '''
409 return self.rhumb.exact
411 def intersection2(self, other, tol=_TOL, **eps):
412 '''I{Iteratively} find the intersection of this and an other rhumb line.
414 @arg other: The other rhumb line (C{RhumbLine}).
415 @kwarg tol: Tolerance for longitudinal convergence (C{degrees}).
416 @kwarg eps: Tolerance for L{pygeodesy.intersection3d3} (C{EPS}).
418 @return: A L{LatLon2Tuple}{(lat, lon)} with the C{lat}- and
419 C{lon}gitude of the intersection point.
421 @raise IntersectionError: No convergence for this B{C{tol}} or
422 no intersection for an other reason.
424 @see: Methods C{distance2} and C{nearestOn4} and function
425 L{pygeodesy.intersection3d3}.
427 @note: Each iteration involves a round trip to this rhumb line's
428 L{ExactTransverseMercator} or L{KTransverseMercator}
429 projection and invoking function L{pygeodesy.intersection3d3}
430 in that domain.
431 '''
432 _xinstanceof(RhumbLineBase, other=other)
433 _xdatum(self.rhumb, other.rhumb, Error=RhumbError)
434 try:
435 if other is self:
436 raise ValueError(_coincident_)
437 # make globals and invariants locals
438 _diff = euclid # approximate length
439 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3
440 _LL2T = LatLon2Tuple
441 _xTMr = self.xTM.reverse # ellipsoidal or spherical
442 _s_3d, s_az = self._xTM3d, self.azi12
443 _o_3d, o_az = other._xTM3d, other.azi12
444 t = fabs(s_az - o_az)
445 if t < EPS or fabs(t - _180_0) < EPS:
446 raise ValueError(_parallel_)
447 # use halfway point as initial estimate
448 p = _LL2T(favg(self.lat1, other.lat1),
449 favg(self.lon1, other.lon1))
450 for i in range(1, _TRIPS):
451 v = _i3d3(_s_3d(p), s_az, # point + bearing
452 _o_3d(p), o_az, useZ=False, **eps)[0]
453 t = _xTMr(v.x, v.y, lon0=p.lon) # PYCHOK Reverse4Tuple
454 d = _diff(t.lon - p.lon, t.lat) # PYCHOK t.lat + p.lat - p.lat
455 p = _LL2T(t.lat + p.lat, t.lon) # PYCHOK t.lon + p.lon = lon0
456 if d < tol:
457 return _LL2T(p.lat, p.lon, iteration=i, # PYCHOK p...
458 name=self.intersection2.__name__)
459 except Exception as x:
460 raise IntersectionError(_no_(_intersection_), cause=x)
461 t = unstr(self.intersection2, tol=tol, **eps)
462 raise IntersectionError(Fmt.no_convergence(d), txt=t)
464 @Property_RO
465 def isLoxodrome(self):
466 '''Is this rhumb line a loxodrome (C{True} if non-meridional and
467 non-parallel, C{False} if parallel or C{None} if meridional)?
469 @see: I{Osborne's} U{2.5 Rumb lines and loxodromes
470 <https://Zenodo.org/record/35392>}, page 37.
471 '''
472 return bool(self._salp) if self._calp else None
474 @Property_RO
475 def lat1(self):
476 '''Get this rhumb line's latitude (C{degrees90}).
477 '''
478 return self._lat1
480 @Property_RO
481 def lon1(self):
482 '''Get this rhumb line's longitude (C{degrees180}).
483 '''
484 return self._lon1
486 @Property_RO
487 def latlon1(self):
488 '''Get this rhumb line's lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).
489 '''
490 return LatLon2Tuple(self.lat1, self.lon1)
492 def _mu22(self, mu12, mu1):
493 mu2 = mu12 + mu1
494 x90 = fabs(mu2) > 90
495 if x90: # reduce to [-180, 180)
496 mu2 = _norm180(mu2)
497 if fabs(mu2) > 90: # point on anti-meridian
498 mu2 = _norm180(_180_0 - mu2)
499 return mu2, x90
501 def nearestOn4(self, lat, lon, tol=_TOL, exact=None, eps=EPS, est=None):
502 '''I{Iteratively} locate the point on this rhumb line nearest to the
503 given point, in part transcoded from I{Karney}'s C++ U{rhumb-intercept
504 <https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}.
506 @arg lat: Latitude of the point (C{degrees}).
507 @arg lon: Longitude of the point (C{degrees}).
508 @kwarg tol: Longitudinal convergence tolerance (C{degrees}) or the
509 distance tolerance (C(meter)) when C{B{exact} is None},
510 respectively C{B{exact} is not None}.
511 @kwarg exact: If C{None}, use a rhumb line perpendicular to this rhumb
512 line, otherwise use an I{exact} C{Geodesic...} from the
513 given point perpendicular to this rhumb line (C{bool} or
514 C{Geodesic...}), see method L{Ellipsoid.geodesic_}.
515 @kwarg eps: Optional tolerance for L{pygeodesy.intersection3d3}
516 (C{EPS}), used only if C{B{exact} is None}.
517 @kwarg est: Optional, initial estimate for the distance C{s13} of
518 the intersection I{along} this rhumb line (C{meter}),
519 used only if C{B{exact} is not None}.
521 @return: A L{NearestOn4Tuple}C{(lat, lon, distance, normal)} with
522 the C{lat}- and C{lon}gitude of the nearest point on and
523 the C{distance} in C{meter} to this rhumb line and the
524 C{normal} azimuth at the intersection.
526 @raise ImportError: I{Karney}'s U{geographiclib
527 <https://PyPI.org/project/geographiclib>}
528 package not found or not installed.
530 @raise IntersectionError: No convergence for this B{C{eps}} or
531 no intersection for an other reason.
533 @see: Methods C{distance2} and C{intersection2} and function
534 L{pygeodesy.intersection3d3}.
535 '''
536 if exact is None:
537 z = _norm180(self.azi12 + _90_0) # perpendicular azimuth
538 rl = RhumbLineBase(self.rhumb, lat, lon, z, caps=Caps.LINE_OFF)
539 p = self.intersection2(rl, tol=tol, eps=eps)
540 t = rl.distance2(p.lat, p.lon)
541 r = NearestOn4Tuple(p.lat, p.lon, t.distance, z,
542 iteration=p.iteration)
543 else: # C{rhumb-intercept}
544 azi = self.azi12
545 szi = self._salp
546 E = self.ellipsoid
547 gX = E.geodesic_(exact=exact)
548 Cs = Caps
549 gm = Cs.AZIMUTH | Cs.LATITUDE_LONGITUDE | Cs.REDUCEDLENGTH | Cs.GEODESICSCALE
550 if est is None: # get an estimate from the perpendicular
551 r = gX.Inverse(self.lat1, self.lon1, lat, lon, outmask=Cs.AZIMUTH_DISTANCE)
552 d, _ = _diff182(r.azi2, azi, K_2_0=True)
553 _, c = sincos2d(d)
554 s12 = c * r.s12 # signed
555 else:
556 s12 = Meter(est=est)
557 try:
558 tol = Float_(tol=tol, low=EPS, high=None)
559 # def _over(p, q): # see @note at RhumbLine[Aux].Position
560 # return (p / (q or _copysign(tol, q))) if isfinite(q) else NAN
562 _s12 = _MODS.fsums.Fsum(s12).fsum2_
563 for i in range(1, _TRIPS): # suffix 1 == C++ 2, 2 == C++ 3
564 p = self.Position(s12) # outmask=Cs.LATITUDE_LONGITUDE
565 r = gX.Inverse(lat, lon, p.lat2, p.lon2, outmask=gm)
566 d, _ = _diff182(azi, r.azi2, K_2_0=True)
567 s, c, s2, c2 = sincos2d_(d, r.lat2)
568 c2 *= E.rocPrimeVertical(r.lat2) # aka rocTransverse
569 s *= (s2 * szi / c2) - (s * r.M21 / r.m12) # XXX _over?
570 s12, t = _s12(c / s) # XXX _over?
571 if fabs(t) < tol: # or fabs(c) < EPS
572 break
573 r = NearestOn4Tuple(r.lat2, r.lon2, s12, r.azi2,
574 iteration=i)
575 except Exception as x: # Fsum(NAN) Value-, ZeroDivisionError
576 raise IntersectionError(_no_(_intersection_), cause=x)
577 return r
579 def Position(self, s12, outmask=0): # PYCHOK no cover
580 '''I{Must be overloaded}, see function C{notOverloaded}.
581 '''
582 _MODS.named.notOverloaded(self, s12, outmask=outmask)
584 @Property_RO
585 def rhumb(self):
586 '''Get this rhumb line's rhumb (L{RhumbAux} or L{Rhumb}.
587 '''
588 return self._rhumb
590 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
591 '''Return this C{RhumbLine} as string.
593 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
594 Trailing zero decimals are stripped for B{C{prec}} values
595 of 1 and above, but kept for negative B{C{prec}} values.
596 @kwarg sep: Separator to join (C{str}).
598 @return: C{RhumbLine} (C{str}).
599 '''
600 d = dict(rhumb=self.rhumb, lat1=self.lat1, lon1=self.lon1,
601 azi12=self.azi12, exact=self.exact,
602 TMorder=self.TMorder, xTM=self.xTM)
603 return sep.join(pairs(itemsorted(d, asorted=False), prec=prec))
605 @property_RO
606 def TMorder(self):
607 '''Get this rhumb line's I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).
608 '''
609 return self.rhumb.TMorder
611 @Property_RO
612 def xTM(self):
613 '''Get this rhumb line's I{Transverse Mercator} projection (L{ExactTransverseMercator}
614 if I{exact} and I{ellipsoidal}, otherwise L{KTransverseMercator} for C{TMorder}).
615 '''
616 E = self.ellipsoid
617 # ExactTransverseMercator doesn't handle spherical earth models
618 return _MODS.etm.ExactTransverseMercator(E) if self.exact and E.isEllipsoidal else \
619 _MODS.ktm.KTransverseMercator(E, TMorder=self.TMorder)
621 def _xTM3d(self, latlon0, z=INT0, V3d=Vector3d):
622 '''(INTERNAL) C{xTM.forward} this C{latlon1} to C{V3d} with B{C{latlon0}}
623 as current intersection estimate and central meridian.
624 '''
625 t = self.xTM.forward(self.lat1 - latlon0.lat, self.lon1, lon0=latlon0.lon)
626 return V3d(t.easting, t.northing, z)
629__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase)
632if __name__ == '__main__':
634 from pygeodesy import printf, RhumbAux as A, Rhumb as R
636 A = A(_EWGS84).Line(30, 0, 45)
637 R = R(_EWGS84).Line(30, 0, 45)
639 for i in range(1, 10):
640 s = .5e6 + 1e6 / i
641 a = A.Position(s).lon2
642 r = R.Position(s).lon2
643 e = (fabs(a - r) / a) if a else 0
644 printf('Position %s vs %s, diff %g', r, a, e)
646 for exact in (None, False, True, None):
647 for est in (None, 1e6):
648 a = A.nearestOn4(60, 0, exact=exact, est=est)
649 r = R.nearestOn4(60, 0, exact=exact, est=est)
650 printf('%s, iteration=%s, exact=%s, est=%s\n%s, iteration=%s',
651 a.toRepr(), a.iteration, exact, est, r.toRepr(), r.iteration, nl=1)
653# % python3 -m pygeodesy.rhumbBase
655# Position 11.614558469016366 vs 11.61455846901637, diff 3.05885e-16
656# Position 7.589823028268422 vs 7.589823028268423, diff 1.17022e-16
657# Position 6.285260674163688 vs 6.285260674163687, diff 1.41311e-16
658# Position 5.639389953251457 vs 5.6393899532514595, diff 4.72486e-16
659# Position 5.253855274357075 vs 5.253855274357073, diff 3.38105e-16
660# Position 4.9976460429038 vs 4.9976460429038045, diff 8.88597e-16
661# Position 4.815033637404729 vs 4.81503363740473, diff 1.84459e-16
662# Position 4.678288217488354 vs 4.678288217488353, diff 1.89851e-16
663# Position 4.572056679062825 vs 4.572056679062826, diff 1.94262e-16
665# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=None
666# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9
668# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=1000000.0
669# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9
671# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=False, est=None
672# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5
674# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=False, est=1000000.0
675# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7
677# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=True, est=None
678# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5
680# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=True, est=1000000.0
681# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7
683# **) MIT License
684#
685# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
686#
687# Permission is hereby granted, free of charge, to any person obtaining a
688# copy of this software and associated documentation files (the "Software"),
689# to deal in the Software without restriction, including without limitation
690# the rights to use, copy, modify, merge, publish, distribute, sublicense,
691# and/or sell copies of the Software, and to permit persons to whom the
692# Software is furnished to do so, subject to the following conditions:
693#
694# The above copyright notice and this permission notice shall be included
695# in all copies or substantial portions of the Software.
696#
697# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
698# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
699# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
700# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
701# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
702# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
703# OTHER DEALINGS IN THE SOFTWARE.