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