Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%
330 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-15 16:36 -0400
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private, ellipsoidal Direct/Inverse geodesy base
5class C{LatLonEllipsoidalBaseDI} and functions.
6'''
7# make sure int/int division yields float quotient, see .basics
8from __future__ import division as _; del _ # PYCHOK semicolon
10from pygeodesy.basics import isLatLon, issubclassof, map2
11from pygeodesy.constants import EPS, MAX, PI, PI2, PI_4, isnear0, isnear1, \
12 _EPSqrt as _TOL, _0_0, _0_5, _1_5, _3_0
13# from pygeodesy.dms import F_DMS # _MODS
14from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, _TOL_M, property_RO
15from pygeodesy.errors import _AssertionError, IntersectionError, _IsnotError, \
16 _or, _ValueError, _xellipsoidal, _xError, _xkwds_not
17from pygeodesy.fmath import favg, fmean_
18from pygeodesy.fsums import Fmt, fsumf_
19from pygeodesy.formy import _isequalTo, opposing, _radical2
20# from pygeodesy.internals import _dunder_nameof # from .lazily
21from pygeodesy.interns import _antipodal_, _concentric_, _ellipsoidal_, \
22 _exceed_PI_radians_, _low_, _near_, \
23 _SPACE_, _too_
24from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, \
25 _dunder_nameof
26from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \
27 Intersection3Tuple, NearestOn2Tuple, \
28 NearestOn8Tuple, _LL4Tuple
29# from pygeodesy.props import property_RO # from .ellipsoidalBase
30# from pygeodesy.streprs import Fmt # from .fsums
31from pygeodesy.units import _fi_j2, _isDegrees, _isHeight, _isRadius, \
32 Radius_, Scalar
33from pygeodesy.utily import m2km, unroll180, _unrollon, _unrollon3, \
34 _Wrap, wrap360
36from math import degrees, radians
38__all__ = _ALL_LAZY.ellipsoidalBaseDI
39__version__ = '24.05.13'
41_polar__ = 'polar?'
42_B2END = _1_5 # _intersect3 bearing to pseudo-end point factor
43_TRIPS = 33 # _intersect3, _intersects2, _nearestOn interations, 6..9 sufficient?
46class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):
47 '''(INTERNAL) Base class for C{ellipsoidal*.LatLon} classes
48 with I{overloaded} C{Direct} and C{Inverse} methods.
49 '''
51 def bearingTo2(self, other, wrap=False):
52 '''Compute the initial and final bearing (forward and reverse
53 azimuth) from this to an other point, using this C{Inverse}
54 method. See methods L{initialBearingTo} and L{finalBearingTo}
55 for more details.
57 @arg other: The other point (C{LatLon}).
58 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
59 B{C{other}} point (C{bool}).
61 @return: A L{Bearing2Tuple}C{(initial, final)}.
63 @raise TypeError: The B{C{other}} point is not C{LatLon}.
65 @raise ValueError: If this and the B{C{other}} point's L{Datum}
66 ellipsoids are not compatible.
67 '''
68 r = self._Inverse(other, wrap)
69 return Bearing2Tuple(r.initial, r.final, name=self.name)
71 def destination(self, distance, bearing, height=None):
72 '''Compute the destination point after having travelled for
73 the given distance from this point along a geodesic given
74 by an initial bearing, using this C{Direct} method. See
75 method L{destination2} for more details.
77 @arg distance: Distance (C{meter}).
78 @arg bearing: Initial bearing in (compass C{degrees360}).
79 @kwarg height: Optional height, overriding the default
80 height (C{meter}, same units as C{distance}).
82 @return: The destination point (C{LatLon}).
83 '''
84 return self._Direct(distance, bearing, self.classof, height).destination
86 def destination2(self, distance, bearing, height=None):
87 '''Compute the destination point and the final bearing (reverse
88 azimuth) after having travelled for the given distance from
89 this point along a geodesic given by an initial bearing,
90 using this C{Direct} method.
92 The distance must be in the same units as this point's datum
93 axes, conventionally C{meter}. The distance is measured on
94 the surface of the ellipsoid, ignoring this point's height.
96 The initial and final bearing (forward and reverse azimuth)
97 are in compass C{degrees360}.
99 The destination point's height and datum are set to this
100 point's height and datum, unless the former is overridden.
102 @arg distance: Distance (C{meter}).
103 @arg bearing: Initial bearing (compass C{degrees360}).
104 @kwarg height: Optional height, overriding the default
105 height (C{meter}, same units as C{distance}).
107 @return: A L{Destination2Tuple}C{(destination, final)}.
108 '''
109 r = self._Direct(distance, bearing, self.classof, height)
110 return self._xnamed(r)
112 def _Direct(self, distance, bearing, LL, height): # overloaded by I{Vincenty}
113 '''(INTERNAL) I{Karney}'s C{Direct} method.
115 @return: A L{Destination2Tuple}C{(destination, final)} or
116 a L{Destination3Tuple}C{(lat, lon, final)} if
117 B{C{LL}} is C{None}.
118 '''
119 g = self.geodesic
120 r = g.Direct3(self.lat, self.lon, bearing, distance)
121 if LL:
122 r = self._Direct2Tuple(LL, height, r)
123 return r
125 def _Direct2Tuple(self, LL, height, r):
126 '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}.
127 '''
128 h = self.height if height is None else height
129 d = LL(*_Wrap.latlon(r.lat, r.lon), height=h, **_xkwds_not(None,
130 datum=self.datum, name=self.name,
131 epoch=self.epoch, reframe=self.reframe))
132 return Destination2Tuple(d, wrap360(r.final))
134 def distanceTo(self, other, wrap=False, **unused): # ignore radius=R_M
135 '''Compute the distance between this and an other point along
136 a geodesic, using this C{Inverse} method. See method
137 L{distanceTo3} for more details.
139 @arg other: The other point (C{LatLon}).
140 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
141 B{C{other}} point (C{bool}).
143 @return: Distance (C{meter}).
145 @raise TypeError: The B{C{other}} point is not C{LatLon}.
147 @raise ValueError: If this and the B{C{other}} point's L{Datum}
148 ellipsoids are not compatible.
149 '''
150 return self._Inverse(other, wrap, azis=False).distance
152 def distanceTo3(self, other, wrap=False):
153 '''Compute the distance, the initial and final bearing along
154 a geodesic between this and an other point, using this
155 C{Inverse} method.
157 The distance is in the same units as this point's datum axes,
158 conventionally meter. The distance is measured on the surface
159 of the ellipsoid, ignoring this point's height.
161 The initial and final bearing (forward and reverse azimuth)
162 are in compass C{degrees360} from North.
164 @arg other: Destination point (C{LatLon}).
165 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
166 B{C{other}} point (C{bool}).
168 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
170 @raise TypeError: The B{C{other}} point is not C{LatLon}.
172 @raise ValueError: If this and the B{C{other}} point's L{Datum}
173 ellipsoids are not compatible.
174 '''
175 return self._xnamed(self._Inverse(other, wrap))
177 def finalBearingOn(self, distance, bearing):
178 '''Compute the final bearing (reverse azimuth) after having
179 travelled for the given distance along a geodesic given by
180 an initial bearing from this point, using this C{Direct}
181 method. See method L{destination2} for more details.
183 @arg distance: Distance (C{meter}).
184 @arg bearing: Initial bearing (compass C{degrees360}).
186 @return: Final bearing (compass C{degrees360}).
187 '''
188 return self._Direct(distance, bearing, None, None).final
190 def finalBearingTo(self, other, wrap=False):
191 '''Compute the final bearing (reverse azimuth) after having
192 travelled along a geodesic from this point to an other
193 point, using this C{Inverse} method. See method
194 L{distanceTo3} for more details.
196 @arg other: The other point (C{LatLon}).
197 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
198 the B{C{other}} point (C{bool}).
200 @return: Final bearing (compass C{degrees360}).
202 @raise TypeError: The B{C{other}} point is not C{LatLon}.
204 @raise ValueError: If this and the B{C{other}} point's L{Datum}
205 ellipsoids are not compatible.
206 '''
207 return self._Inverse(other, wrap).final
209 @property_RO
210 def geodesic(self): # overloaded by I{Karney}'s, N/A for I{Vincenty}
211 '''N/A, invalid (C{None} I{always}).
212 '''
213 return None # PYCHOK no cover
215 def _g_gl_p3(self, point, other, exact, wrap):
216 '''(INTERNAL) Helper for methods C{.intersecant2} and C{.plumbTo}.
217 '''
218 p = _unrollon(self, self.others(point=point), wrap=wrap)
219 g = self.datum.ellipsoid.geodesic_(exact=exact)
220 gl = g._DirectLine( p, other) if _isDegrees(other) else \
221 g._InverseLine(p, self.others(other), wrap)
222 return g, gl, p
224 def initialBearingTo(self, other, wrap=False):
225 '''Compute the initial bearing (forward azimuth) to travel
226 along a geodesic from this point to an other point,
227 using this C{Inverse} method. See method L{distanceTo3}
228 for more details.
230 @arg other: The other point (C{LatLon}).
231 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
232 the B{C{other}} point (C{bool}).
234 @return: Initial bearing (compass C{degrees360}).
236 @raise TypeError: The B{C{other}} point is not C{LatLon}.
238 @raise ValueError: If this and the B{C{other}} point's L{Datum}
239 ellipsoids are not compatible.
240 '''
241 return self._Inverse(other, wrap).initial
243 def intermediateTo(self, other, fraction, height=None, wrap=False):
244 '''Return the point at given fraction along the geodesic between
245 this and an other point, using this C{Direct} and C{Inverse}
246 methods.
248 @arg other: The other point (C{LatLon}).
249 @arg fraction: Fraction between both points (C{scalar}, 0.0
250 at this and 1.0 at the other point.
251 @kwarg height: Optional height, overriding the fractional
252 height (C{meter}).
253 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
254 B{C{other}} point (C{bool}).
256 @return: Intermediate point (C{LatLon}).
258 @raise TypeError: The B{C{other}} point is not C{LatLon}.
260 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.
262 @raise ValueError: If this and the B{C{other}} point's L{Datum}
263 ellipsoids are not compatible.
265 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and
266 C{rhumbMidpointTo}.
267 '''
268 f = Scalar(fraction=fraction)
269 if isnear0(f):
270 r = self
271 elif isnear1(f) and not wrap:
272 r = self.others(other)
273 else: # negative fraction OK
274 t = self.distanceTo3(other, wrap=wrap)
275 h = self._havg(other, f=f, h=height)
276 r = self.destination(t.distance * f, t.initial, height=h)
277 return r
279 def intersecant2(self, circle, point, other, exact=False, height=None, # PYCHOK signature
280 wrap=False, tol=_TOL):
281 '''Compute the intersections of a circle and a geodesic (line) given as
282 two points or as a point and bearing.
284 @arg circle: Radius of the circle centered at this location (C{meter},
285 conventionally) or a point on the circle (this C{LatLon}).
286 @arg point: A location on the geodesic (this C{LatLon}).
287 @arg other: An other point on the geodesic (this C{LatLon}) or the
288 (forward) bearing at the B{C{point}} (compass C{degrees}).
289 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}), see
290 method L{Ellipsoid.geodesic_}.
291 @kwarg height: Optional height for the intersection points (C{meter},
292 conventionally) or C{None} for interpolated heights.
293 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}},
294 B{C{circle}} and/or B{C{other}} (C{bool}).
295 @kwarg tol: Convergence tolerance (C{scalar}).
297 @return: 2-Tuple of the intersection points (representing a geodesic chord),
298 each an instance of this C{LatLon}. Both points are the same
299 instance if the geodesic (line) is tangential to the circle.
301 @raise IntersectionError: The circle and geodesic do not intersect.
303 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
304 or B{C{other}} invalid.
306 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{exact}} or
307 B{C{height}}.
309 @see: Method L{rhumbIntersecant2<LatLonBase.rhumbIntersecant2>}.
310 '''
311 try:
312 g, gl, p = self._g_gl_p3(point, other, exact, wrap)
313 r = Radius_(circle=circle) if _isRadius(circle) else \
314 g._Inverse(self, self.others(circle=circle), wrap).s12
316 P, Q = _MODS.geodesicw._Intersecant2(gl, self.lat, self.lon, r, tol=tol,
317 form=_MODS.dms.F_DMS)
318 return self._intersecend2(p, other, wrap, height, g, P, Q,
319 self.intersecant2)
320 except (TypeError, ValueError) as x:
321 raise _xError(x, center=self, circle=circle, point=point, other=other,
322 exact=exact, wrap=wrap)
324 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty}
325 '''(INTERNAL) I{Karney}'s C{Inverse} method.
327 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
328 '''
329 _ = self.ellipsoids(other)
330 g = self.geodesic
331 _, lon = unroll180(self.lon, other.lon, wrap=wrap)
332 return g.Inverse3(self.lat, self.lon, other.lat, lon)
334 def nearestOn8(self, points, closed=False, height=None, wrap=False,
335 equidistant=None, tol=_TOL_M):
336 '''I{Iteratively} locate the point on a path or polygon closest
337 to this point.
339 @arg points: The path or polygon points (C{LatLon}[]).
340 @kwarg closed: Optionally, close the polygon (C{bool}).
341 @kwarg height: Optional height, overriding the height of this and
342 all other points (C{meter}, conventionally). If
343 B{C{height}} is C{None}, the height of each point
344 is taken into account for distances.
345 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
346 B{C{points}} (C{bool}).
348 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start,
349 end, initial, final)} with C{distance} in C{meter},
350 conventionally and with the C{closest}, the C{start}
351 the C{end} point each an instance of this C{LatLon}.
353 @raise PointsError: Insufficient number of B{C{points}}.
355 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid.
357 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible
358 or no convergence for the given B{C{tol}}.
360 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}.
361 '''
362 _d3 = self.distanceTo3 # Distance3Tuple
363 _n3 = _nearestOn3
364 try:
365 Ps = self.PointsIter(points, loop=1, wrap=wrap)
366 p1 = c = s = e = Ps[0]
367 _ = self.ellipsoids(p1)
368 c3 = _d3(c, wrap=wrap) # XXX wrap=False?
370 except (TypeError, ValueError) as x:
371 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol,
372 closed=closed, height=height, wrap=wrap)
374 # get the azimuthal equidistant projection, once
375 A = _Equidistant00(equidistant, c)
376 b = _Box(c, c3.distance)
377 m = f = i = 0 # p1..p2 == points[i]..[j]
379 kwds = dict(within=True, height=height, tol=tol,
380 LatLon=self.classof, # this LatLon
381 datum=self.datum, epoch=self.epoch, reframe=self.reframe)
382 try:
383 for j, p2 in Ps.enumerate(closed=closed):
384 if wrap and j != 0: # not Ps.looped
385 p2 = _unrollon(p1, p2)
386 # skip edge if no overlap with box around closest
387 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon):
388 p, t, _ = _n3(self, p1, p2, A, **kwds)
389 d3 = _d3(p, wrap=False) # already unrolled
390 if d3.distance < c3.distance:
391 c3, c, s, e, f = d3, p, p1, p2, (i + t)
392 b = _Box(c, c3.distance)
393 m = max(m, c.iteration)
394 p1, i = p2, j
396 except (TypeError, ValueError) as x:
397 raise _xError(x, Fmt.SQUARE(points=i), p1,
398 Fmt.SQUARE(points=j), p2, this=self, tol=tol,
399 closed=closed, height=height, wrap=wrap)
401 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6
403 n = self.nearestOn8.__name__
404 c.rename(n)
405 if s is not c:
406 s = s.copy(name=n)
407 if e is not c:
408 e = e.copy(name=n)
409 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final,
410 iteration=m) # ._iteration for tests
412 def plumbTo(self, point, other, exact=False, height=None, # PYCHOK signature
413 wrap=False, tol=_TOL):
414 '''Compute the I{perpendicular} intersection of a geodesic (line) with
415 a geodesic from this point.
417 @arg point: A location (C{LatLon}) on the geodesic (line).
418 @arg other: An other point (C{LatLon}) on the geodesic (line) or the
419 (forward) bearing at the B{C{point}} (compass C{degrees}).
420 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}),
421 see method L{Ellipsoid.geodesic_}.
422 @kwarg height: Optional height for the intersection point (C{meter},
423 conventionally) or C{None} for an interpolated height.
424 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
425 and/or B{C{other}} (C{bool}).
426 @kwarg tol: Convergence tolerance (C{scalar}).
428 @return: The intersection point, an instance of this C{LatLon}.
430 @raise TypeError: If B{C{point}} or B{C{other}} not this C{LatLon}.
432 @raise UnitError: Invalid B{C{other}}, B{C{exact}} or B{C{height}}.
433 '''
434 try:
435 g, gl, p = self._g_gl_p3(point, other, exact, wrap)
437 P = _MODS.geodesicw._PlumbTo(gl, self.lat, self.lon, tol=tol)
438 h = self._havg(p, h=height)
439 p = self.classof(P.lat2, P.lon2, datum=self.datum, height=h) # name=n
440 p._iteration = P.iteration
441 except (TypeError, ValueError) as x:
442 raise _xError(x, plumb=self, point=point, other=other,
443 exact=exact, wrap=wrap)
444 return p
447class _Box(object):
448 '''Bounding box around a C{LatLon} point.
450 @see: Function C{_box4} in .clipy.py.
451 '''
452 def __init__(self, center, distance):
453 '''New L{_Box} around point.
455 @arg center: The center point (C{LatLon}).
456 @arg distance: Radius, half-size of the box
457 (C{meter}, conventionally)
458 '''
459 m = Radius_(distance=distance)
460 E = center.ellipsoid() # XXX see above
461 d = E.m2degrees(m) + _0_5 # some margin
462 self._N = center.lat + d
463 self._S = center.lat - d
464 self._E = center.lon + d
465 self._W = center.lon - d
467 def overlaps(self, lat1, lon1, lat2, lon2):
468 '''Check whether this box overlaps an other box.
470 @arg lat1: Latitude of a box corner (C{degrees}).
471 @arg lon1: Longitude of a box corner (C{degrees}).
472 @arg lat2: Latitude of the opposing corner (C{degrees}).
473 @arg lon2: Longitude of the opposing corner (C{degrees}).
475 @return: C{True} if there is some overlap, C{False}
476 otherwise (C{bool}).
477 '''
478 if lat1 > lat2:
479 lat1, lat2 = lat2, lat1
480 if lat1 > self._N or lat2 < self._S:
481 return False
482 if lon1 > lon2:
483 lon1, lon2 = lon2, lon1
484 if lon1 > self._E or lon2 < self._W:
485 return False
486 return True
489class _Tol(object):
490 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}.
491 '''
492 _deg = 0 # tol in degrees
493 _lat = 0
494 _m = 0 # tol in meter
495 _min = MAX # degrees
496 _prev = None
497 _r = 0
499 def __init__(self, tol_m, E, lat, *lats):
500 '''New L{_Tol}.
502 @arg tol_m: Tolerance (C{meter}, only).
503 @arg E: Earth ellipsoid (L{Ellipsoid}).
504 @arg lat: Latitude (C{degrees}).
505 @arg lats: Additional latitudes (C{degrees}).
506 '''
507 self._lat = fmean_(lat, *lats) if lats else lat
508 self._r = max(EPS, E.rocMean(self._lat))
509 self._m = max(EPS, tol_m)
510 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!
512 @property_RO
513 def degrees(self):
514 '''Get this tolerance in C{degrees}.
515 '''
516 return self._deg
518 def degrees2m(self, deg):
519 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.
520 '''
521 return self.radius * radians(deg) / PI2 # avoid degrees2m!
523 def degError(self, Error=_ValueError):
524 '''Compose an error with the C{deg}rees minimum.
525 '''
526 return self.mError(self.degrees2m(self._min), Error=Error)
528 def done(self, deg):
529 '''Check C{deg} vs. tolerance and previous value.
530 '''
531 if deg < self._deg or deg == self._prev:
532 return True
533 self._min = min(self._min, deg)
534 self._prev = deg
535 return False
537 @property_RO
538 def lat(self):
539 '''Get the mean latitude in C{degrees}.
540 '''
541 return self._lat
543 def mError(self, m, Error=_ValueError):
544 '''Compose an error with B{C{m}}eter minimum.
545 '''
546 t = _SPACE_(Fmt.tolerance(self.meter), _too_(_low_))
547 if m2km(m) > self.meter:
548 t = _or(t, _antipodal_, _near_(_polar__))
549 return Error(Fmt.no_convergence(m), txt=t)
551 @property_RO
552 def meter(self):
553 '''Get this tolerance in C{meter}.
554 '''
555 return self._m
557 @property_RO
558 def radius(self):
559 '''Get the earth radius in C{meter}.
560 '''
561 return self._r
563 def reset(self):
564 '''Reset tolerances.
565 '''
566 self._min = MAX
567 self._prev = None
570def _Equidistant00(equidistant, p1):
571 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
572 '''
573 if equidistant is None or not callable(equidistant):
574 equidistant = p1.Equidistant
575 else:
576 Eqs = _MODS.azimuthal._Equidistants
577 if not issubclassof(equidistant, *Eqs): # PYCHOK no cover
578 t = map2(_dunder_nameof, Eqs)
579 raise _IsnotError(*t, equidistant=equidistant)
580 return equidistant(0, 0, p1.datum)
583def intersecant2(center, circle, point, other, **exact_height_wrap_tol):
584 '''Compute the intersections of a circle and a geodesic given as two points
585 or as a point and (forward) bearing.
587 @arg center: Center of the circle (C{LatLon}).
588 @arg circle: Radius of the circle (C{meter}, conventionally) or a point on
589 the circle (C{LatLon}, as B{C{center}}).
590 @arg point: A point of the geodesic (C{LatLon}, as B{C{center}}).
591 @arg other: An other point of the geodesic (C{LatLon}, as B{C{center}}) or
592 the (forward) bearing at the B{C{point}} (compass C{degrees}).
593 @kwarg exact_height_wrap_tol: Optional keyword arguments, see below.
595 @raise NotImplementedError: Method C{intersecant2} not available.
597 @raise TypeError: If B{C{center}}, B{C{point}} or B{C{circle}} or B{C{other}}
598 points not ellipsoidal or not compatible with B{C{center}}.
600 @see: Method C{LatLon.intersecant2} of class L{ellipsoidalExact.LatLon},
601 L{ellipsoidalKarney.LatLon} or L{ellipsoidalVincenty.LatLon}.
602 '''
603 if not isLatLon(center, ellipsoidal=True): # isinstance(center, LatLonEllipsoidalBase)
604 raise _IsnotError(_ellipsoidal_, center=center)
605 return center.intersecant2(circle, point, other, **exact_height_wrap_tol)
608def _intersect3(s1, end1, s2, end2, height=None, wrap=False, # MCCABE 16 was=True
609 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
610 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method
611 L{intersection3}, separated to allow callers to embellish any exceptions.
612 '''
613 _LLS = _MODS.sphericalTrigonometry.LatLon
614 _si = _MODS.sphericalTrigonometry._intersect
615 _vi3 = _MODS.vector3d._intersect3d3
617 def _b_d(s, e, w, t, h=_0_0):
618 # compute opposing and distance
619 t = s.classof(t.lat, t.lon, height=h, name=t.name)
620 t = s.distanceTo3(t, wrap=w) # Distance3Tuple
621 b = opposing(e, t.initial) # "before" start
622 return b, t.distance
624 def _b_e(s, e, w, t):
625 # compute an end point along the initial bearing about
626 # 1.5 times the distance to the gu-/estimate, at least
627 # 1/8 and at most 3/8 of the earth perimeter like the
628 # radians in .sphericalTrigonometry._int3d2 and bearing
629 # comparison in .sphericaltrigonometry._intb
630 b, d = _b_d(s, e, w, t, h=t.height)
631 m = s.ellipsoid().R2x * PI_4 # authalic exact
632 d = min(max(d * _B2END, m), m * _3_0)
633 e = s.destination(d, e)
634 return b, (_unrollon(s, e) if w else e)
636 def _e_ll(s, e, w, **end):
637 # return 2-tuple (end, False if bearing else True)
638 ll = not _isDegrees(e)
639 if ll:
640 e = s.others(**end)
641 if w: # unroll180 == .karney._unroll2
642 e = _unrollon(s, e)
643 return e, ll
645 def _o(o, b, n, s, t, e):
646 # determine C{o}utside before, on or after start point
647 if not o: # intersection may be on start
648 if _isequalTo(s, t, eps=e.degrees):
649 return o
650 return -n if b else n
652 E = s1.ellipsoids(s2)
654 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
655 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
657 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),
658 s2.lat, (e2.lat if ll2 else s2.lat))
660 # get the azimuthal equidistant projection
661 A = _Equidistant00(equidistant, s1)
663 # gu-/estimate initial intersection, spherically ...
664 t = _si(_LLS(s1), (_LLS(e1) if ll1 else e1),
665 _LLS(s2), (_LLS(e2) if ll2 else e2),
666 height=height, wrap=False, LatLon=_LLS) # unrolled already
667 h, n = t.height, t.name
669 if not ll1:
670 b1, e1 = _b_e(s1, e1, wrap, t)
671 if not ll2:
672 b2, e2 = _b_e(s2, e2, wrap, t)
674 # ... and iterate as Karney describes, for references
675 # @see: Function L{ellipsoidalKarney.intersection3}.
676 for i in range(1, _TRIPS):
677 A.reset(t.lat, t.lon) # gu-/estimate as origin
678 # convert start and end points to projection
679 # space and compute an intersection there
680 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),
681 eps=e.meter, useZ=False)
682 # convert intersection back to geodetic
683 t, d = A._reverse2(v)
684 if e.done(d): # below tol or unchanged?
685 break
686 else:
687 raise e.degError(Error=IntersectionError)
689 # like .sphericalTrigonometry._intersect, if this intersection
690 # is "before" the first point, use the antipodal intersection
691 if not (ll1 or ll2): # end1 and end2 are an initial bearing
692 b1, _ = _b_d(s1, end1, wrap, t)
693 if b1:
694 t = t.antipodal()
695 b1 = not b1
696 b2, _ = _b_d(s2, end2, wrap, t)
698 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1,
699 iteration=i, name=n)
700 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),
701 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))
704def _intersection3(start1, end1, start2, end2, height=None, wrap=False, # was=True
705 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
706 '''(INTERNAL) Iteratively compute the intersection point of two lines,
707 each defined by two (ellipsoidal) points or an (ellipsoidal) start
708 point and an initial bearing from North.
709 '''
710 s1 = _xellipsoidal(start1=start1)
711 s2 = s1.others(start2=start2)
712 try:
713 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,
714 equidistant=equidistant, tol=tol,
715 LatLon=LatLon, **LatLon_kwds)
716 except (TypeError, ValueError) as x:
717 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
720def _intersections2(center1, radius1, center2, radius2, height=None, wrap=False, # was=True
721 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
722 '''(INTERNAL) Iteratively compute the intersection points of two circles,
723 each defined by an (ellipsoidal) center point and a radius.
724 '''
725 c1 = _xellipsoidal(center1=center1)
726 c2 = c1.others(center2=center2)
727 try:
728 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,
729 equidistant=equidistant, tol=tol,
730 LatLon=LatLon, **LatLon_kwds)
731 except (TypeError, ValueError) as x:
732 raise _xError(x, center1=center1, radius1=radius1,
733 center2=center2, radius2=radius2)
736def _intersects2(c1, radius1, c2, radius2, height=None, wrap=False, # MCCABE 16 was=True
737 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
738 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}
739 above, separated to allow callers to embellish any exceptions.
740 '''
741 _LLS = _MODS.sphericalTrigonometry.LatLon
742 _si2 = _MODS.sphericalTrigonometry._intersects2
743 _vi2 = _MODS.vector3d._intersects2
745 def _ll4(t, h, n, c):
746 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c,
747 iteration=t.iteration, name=n)
749 r1 = Radius_(radius1=radius1)
750 r2 = Radius_(radius2=radius2)
752 E = c1.ellipsoids(c2)
753 # get the azimuthal equidistant projection
754 A = _Equidistant00(equidistant, c1)
756 if r1 < r2:
757 c1, c2 = c2, c1
758 r1, r2 = r2, r1
760 if r1 > (min(E.b, E.a) * PI):
761 raise _ValueError(_exceed_PI_radians_)
763 if wrap: # unroll180 == .karney._unroll2
764 c2 = _unrollon(c1, c2)
766 # distance between centers and radii are
767 # measured along the ellipsoid's surface
768 m = c1.distanceTo(c2, wrap=False) # meter
769 if m < max(r1 - r2, EPS):
770 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
771 if fsumf_(r1, r2, -m) < 0:
772 raise IntersectionError(_too_(Fmt.distant(m)))
774 f = _radical2(m, r1, r2).ratio # "radical fraction"
775 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))
777 # gu-/estimate initial intersections, spherically ...
778 t1, t2 = _si2(_LLS(c1), r1, _LLS(c2), r2, radius=e.radius,
779 height=height, wrap=False, too_d=m) # unrolled already
780 h, n = t1.height, t1.name
782 # ... and iterate as Karney describes, for references
783 # @see: Function L{ellipsoidalKarney.intersections2}.
784 ts, ta = [], None
785 for t in ((t1,) if t1 is t2 else (t1, t2)):
786 for i in range(1, _TRIPS):
787 A.reset(t.lat, t.lon) # gu-/estimate as origin
788 # convert centers to projection space and
789 # compute the intersections there
790 t1, t2 = A._forwards(c1, c2)
791 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,
792 t2, r2, # XXX * t2.scale?,
793 sphere=False, too_d=m)
794 # convert intersections back to geodetic
795 if v1 is v2: # abutting
796 t, d = A._reverse2(v1)
797 else: # consider the closer intersection
798 t1, d1 = A._reverse2(v1)
799 t2, d2 = A._reverse2(v2)
800 t, d = (t1, d1) if d1 < d2 else (t2, d2)
801 if e.done(d): # below tol or unchanged?
802 t._iteration = i # _NamedTuple._iteration
803 ts.append(t)
804 if v1 is v2: # abutting
805 ta = t
806 break
807 else:
808 raise e.degError(Error=IntersectionError)
809 e.reset()
811 if ta: # abutting circles
812 pass # PYCHOK no cover
813 elif len(ts) == 2:
814 return (_ll4(ts[0], h, n, c1),
815 _ll4(ts[1], h, n, c2))
816 elif len(ts) == 1: # PYCHOK no cover
817 ta = ts[0] # assume abutting
818 else: # PYCHOK no cover
819 raise _AssertionError(ts=ts)
820 r = _ll4(ta, h, n, c1)
821 return r, r
824def _nearestOn2(p, point1, point2, within=True, height=None, wrap=False, # was=True
825 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):
826 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,
827 separated to allow callers to embellish any exceptions.
828 '''
829 p1 = p.others(point1=point1)
830 p2 = p.others(point2=point2)
832 _ = p.ellipsoids(p1)
833# E = p.ellipsoids(p2) # done in _nearestOn3
835 # get the azimuthal equidistant projection
836 A = _Equidistant00(equidistant, p)
838 p1, p2, _ = _unrollon3(p, p1, p2, wrap) # XXX don't unroll?
839 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height,
840 tol=tol, **LatLon_and_kwds)
841 return NearestOn2Tuple(r, f)
844def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M,
845 LatLon=None, **LatLon_kwds):
846 # Only in function C{_nearestOn2} and method C{nearestOn8} above
847 _LLS = _MODS.sphericalNvector.LatLon
848 _V3d = _MODS.vector3d.Vector3d
849 _vn2 = _MODS.vector3d._nearestOn2
851 E = p.ellipsoids(p2)
852 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)
854 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
855 # using sphericalNvector.LatLon.nearestOn for within=False support
856 t = _LLS(p).nearestOn(_LLS(p1), _LLS(p2), within=within,
857 height=height)
858 n, h = t.name, t.height
859 if height is None:
860 h1 = p1.height # use heights as pseudo-Z in projection space
861 h2 = p2.height # to be included in the closest function
862 h0 = favg(h1, h2)
863 else: # ignore heights in distances, Z=0
864 h0 = h1 = h2 = _0_0
866 # ... and iterate to find the closest (to the origin with .z
867 # to interpolate height) as Karney describes, for references
868 # @see: Function L{ellipsoidalKarney.nearestOn}.
869 vp, f = _V3d(_0_0, _0_0, h0), None
870 for i in range(1, _TRIPS):
871 A.reset(t.lat, t.lon) # gu-/estimate as origin
872 # convert points to projection space and compute
873 # the nearest one (and its height) there
874 s, t = A._forwards(p1, p2)
875 v, f = _vn2(vp, _V3d(s.x, s.y, h1),
876 _V3d(t.x, t.y, h2), within=within)
877 # convert nearest one back to geodetic
878 t, d = A._reverse2(v)
879 if e.done(d): # below tol or unchanged?
880 break
881 else:
882 raise e.degError()
884 if height is None:
885 h = v.z # nearest
886 elif _isHeight(height):
887 h = height
888 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,
889 iteration=i, name=n)
890 return r, f, e # fraction or None
893__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2)
895# **) MIT License
896#
897# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
898#
899# Permission is hereby granted, free of charge, to any person obtaining a
900# copy of this software and associated documentation files (the "Software"),
901# to deal in the Software without restriction, including without limitation
902# the rights to use, copy, modify, merge, publish, distribute, sublicense,
903# and/or sell copies of the Software, and to permit persons to whom the
904# Software is furnished to do so, subject to the following conditions:
905#
906# The above copyright notice and this permission notice shall be included
907# in all copies or substantial portions of the Software.
908#
909# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
910# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
911# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
912# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
913# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
914# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
915# OTHER DEALINGS IN THE SOFTWARE.