Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%
337 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-24 11:18 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-24 11:18 -0500
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
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, property_RO, _TOL_M
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_
19# from pygeodesy.formy import _isequalTo, opposing, _radical2 # _MODS
20from pygeodesy.interns import _antipodal_, _concentric_, _ellipsoidal_, \
21 _exceed_PI_radians_, _low_, _near_, _SPACE_, _too_
22from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
23from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \
24 Intersection3Tuple, NearestOn2Tuple, \
25 NearestOn8Tuple, _LL4Tuple
26# from pygeodesy.props import property_RO # from .ellipsoidalBase
27# from pygeodesy.streprs import Fmt # from .fsums
28from pygeodesy.units import _fi_j2, _isDegrees, _isHeight, _isRadius, \
29 Radius_, Scalar
30from pygeodesy.utily import m2km, unroll180, _unrollon, _unrollon3, \
31 _Wrap, wrap360
33from math import degrees, radians
35__all__ = _ALL_LAZY.ellipsoidalBaseDI
36__version__ = '23.12.12'
38_polar__ = 'polar?'
39_B2END = _1_5 # _intersect3 bearing to pseudo-end point factor
40_TRIPS = 33 # _intersect3, _intersects2, _nearestOn interations, 6..9 sufficient?
43class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):
44 '''(INTERNAL) Base class for C{ellipsoidal*.LatLon} classes
45 with I{overloaded} C{Direct} and C{Inverse} methods.
46 '''
48 def bearingTo2(self, other, wrap=False):
49 '''Compute the initial and final bearing (forward and reverse
50 azimuth) from this to an other point, using this C{Inverse}
51 method. See methods L{initialBearingTo} and L{finalBearingTo}
52 for more details.
54 @arg other: The other point (C{LatLon}).
55 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
56 B{C{other}} point (C{bool}).
58 @return: A L{Bearing2Tuple}C{(initial, final)}.
60 @raise TypeError: The B{C{other}} point is not C{LatLon}.
62 @raise ValueError: If this and the B{C{other}} point's L{Datum}
63 ellipsoids are not compatible.
64 '''
65 r = self._Inverse(other, wrap)
66 return Bearing2Tuple(r.initial, r.final, name=self.name)
68 def destination(self, distance, bearing, height=None):
69 '''Compute the destination point after having travelled for
70 the given distance from this point along a geodesic given
71 by an initial bearing, using this C{Direct} method. See
72 method L{destination2} for more details.
74 @arg distance: Distance (C{meter}).
75 @arg bearing: Initial bearing in (compass C{degrees360}).
76 @kwarg height: Optional height, overriding the default
77 height (C{meter}, same units as C{distance}).
79 @return: The destination point (C{LatLon}).
80 '''
81 return self._Direct(distance, bearing, self.classof, height).destination
83 def destination2(self, distance, bearing, height=None):
84 '''Compute the destination point and the final bearing (reverse
85 azimuth) after having travelled for the given distance from
86 this point along a geodesic given by an initial bearing,
87 using this C{Direct} method.
89 The distance must be in the same units as this point's datum
90 axes, conventionally C{meter}. The distance is measured on
91 the surface of the ellipsoid, ignoring this point's height.
93 The initial and final bearing (forward and reverse azimuth)
94 are in compass C{degrees360}.
96 The destination point's height and datum are set to this
97 point's height and datum, unless the former is overridden.
99 @arg distance: Distance (C{meter}).
100 @arg bearing: Initial bearing (compass C{degrees360}).
101 @kwarg height: Optional height, overriding the default
102 height (C{meter}, same units as C{distance}).
104 @return: A L{Destination2Tuple}C{(destination, final)}.
105 '''
106 r = self._Direct(distance, bearing, self.classof, height)
107 return self._xnamed(r)
109 def _Direct(self, distance, bearing, LL, height): # overloaded by I{Vincenty}
110 '''(INTERNAL) I{Karney}'s C{Direct} method.
112 @return: A L{Destination2Tuple}C{(destination, final)} or
113 a L{Destination3Tuple}C{(lat, lon, final)} if
114 B{C{LL}} is C{None}.
115 '''
116 g = self.geodesic
117 r = g.Direct3(self.lat, self.lon, bearing, distance)
118 if LL:
119 r = self._Direct2Tuple(LL, height, r)
120 return r
122 def _Direct2Tuple(self, LL, height, r):
123 '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}.
124 '''
125 h = self.height if height is None else height
126 d = LL(*_Wrap.latlon(r.lat, r.lon), height=h, **_xkwds_not(None,
127 datum=self.datum, name=self.name,
128 epoch=self.epoch, reframe=self.reframe))
129 return Destination2Tuple(d, wrap360(r.final))
131 def distanceTo(self, other, wrap=False, **unused): # ignore radius=R_M
132 '''Compute the distance between this and an other point along
133 a geodesic, using this C{Inverse} method. See method
134 L{distanceTo3} for more details.
136 @arg other: The other point (C{LatLon}).
137 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
138 B{C{other}} point (C{bool}).
140 @return: Distance (C{meter}).
142 @raise TypeError: The B{C{other}} point is not C{LatLon}.
144 @raise ValueError: If this and the B{C{other}} point's L{Datum}
145 ellipsoids are not compatible.
146 '''
147 return self._Inverse(other, wrap, azis=False).distance
149 def distanceTo3(self, other, wrap=False):
150 '''Compute the distance, the initial and final bearing along
151 a geodesic between this and an other point, using this
152 C{Inverse} method.
154 The distance is in the same units as this point's datum axes,
155 conventionally meter. The distance is measured on the surface
156 of the ellipsoid, ignoring this point's height.
158 The initial and final bearing (forward and reverse azimuth)
159 are in compass C{degrees360} from North.
161 @arg other: Destination point (C{LatLon}).
162 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
163 B{C{other}} point (C{bool}).
165 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
167 @raise TypeError: The B{C{other}} point is not C{LatLon}.
169 @raise ValueError: If this and the B{C{other}} point's L{Datum}
170 ellipsoids are not compatible.
171 '''
172 return self._xnamed(self._Inverse(other, wrap))
174 def finalBearingOn(self, distance, bearing):
175 '''Compute the final bearing (reverse azimuth) after having
176 travelled for the given distance along a geodesic given by
177 an initial bearing from this point, using this C{Direct}
178 method. See method L{destination2} for more details.
180 @arg distance: Distance (C{meter}).
181 @arg bearing: Initial bearing (compass C{degrees360}).
183 @return: Final bearing (compass C{degrees360}).
184 '''
185 return self._Direct(distance, bearing, None, None).final
187 def finalBearingTo(self, other, wrap=False):
188 '''Compute the final bearing (reverse azimuth) after having
189 travelled along a geodesic from this point to an other
190 point, using this C{Inverse} method. See method
191 L{distanceTo3} for more details.
193 @arg other: The other point (C{LatLon}).
194 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
195 the B{C{other}} point (C{bool}).
197 @return: Final bearing (compass C{degrees360}).
199 @raise TypeError: The B{C{other}} point is not C{LatLon}.
201 @raise ValueError: If this and the B{C{other}} point's L{Datum}
202 ellipsoids are not compatible.
203 '''
204 return self._Inverse(other, wrap).final
206 @property_RO
207 def geodesic(self): # overloaded by I{Karney}'s, N/A for I{Vincenty}
208 '''N/A, invalid (C{None} I{always}).
209 '''
210 return None # PYCHOK no cover
212 def _g_gl_p3(self, point, other, exact, wrap):
213 '''(INTERNAL) Helper for methods C{.intersecant2} and C{.plumbTo}.
214 '''
215 p = _unrollon(self, self.others(point=point), wrap=wrap)
216 g = self.datum.ellipsoid.geodesic_(exact=exact)
217 gl = g._DirectLine( p, other) if _isDegrees(other) else \
218 g._InverseLine(p, self.others(other), wrap)
219 return g, gl, p
221 def initialBearingTo(self, other, wrap=False):
222 '''Compute the initial bearing (forward azimuth) to travel
223 along a geodesic from this point to an other point,
224 using this C{Inverse} method. See method L{distanceTo3}
225 for more details.
227 @arg other: The other point (C{LatLon}).
228 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
229 the B{C{other}} point (C{bool}).
231 @return: Initial bearing (compass C{degrees360}).
233 @raise TypeError: The B{C{other}} point is not C{LatLon}.
235 @raise ValueError: If this and the B{C{other}} point's L{Datum}
236 ellipsoids are not compatible.
237 '''
238 return self._Inverse(other, wrap).initial
240 def intermediateTo(self, other, fraction, height=None, wrap=False):
241 '''Return the point at given fraction along the geodesic between
242 this and an other point, using this C{Direct} and C{Inverse}
243 methods.
245 @arg other: The other point (C{LatLon}).
246 @arg fraction: Fraction between both points (C{scalar}, 0.0
247 at this and 1.0 at the other point.
248 @kwarg height: Optional height, overriding the fractional
249 height (C{meter}).
250 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
251 B{C{other}} point (C{bool}).
253 @return: Intermediate point (C{LatLon}).
255 @raise TypeError: The B{C{other}} point is not C{LatLon}.
257 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.
259 @raise ValueError: If this and the B{C{other}} point's L{Datum}
260 ellipsoids are not compatible.
262 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and
263 C{rhumbMidpointTo}.
264 '''
265 f = Scalar(fraction=fraction)
266 if isnear0(f):
267 r = self
268 elif isnear1(f) and not wrap:
269 r = self.others(other)
270 else: # negative fraction OK
271 t = self.distanceTo3(other, wrap=wrap)
272 h = self._havg(other, f=f, h=height)
273 r = self.destination(t.distance * f, t.initial, height=h)
274 return r
276 def intersecant2(self, circle, point, other, exact=False, height=None, # PYCHOK signature
277 wrap=False, tol=_TOL):
278 '''Compute the intersections of a circle and a geodesic (line) given as
279 two points or as a point and bearing.
281 @arg circle: Radius of the circle centered at this location (C{meter},
282 conventionally) or a point on the circle (this C{LatLon}).
283 @arg point: A location on the geodesic (this C{LatLon}).
284 @arg other: An other point on the geodesic (this C{LatLon}) or the
285 (forward) bearing at the B{C{point}} (compass C{degrees}).
286 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}), see
287 method L{Ellipsoid.geodesic_}.
288 @kwarg height: Optional height for the intersection points (C{meter},
289 conventionally) or C{None} for interpolated heights.
290 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}},
291 B{C{circle}} and/or B{C{other}} (C{bool}).
292 @kwarg tol: Convergence tolerance (C{scalar}).
294 @return: 2-Tuple of the intersection points (representing a geodesic chord),
295 each an instance of this C{LatLon}. Both points are the same
296 instance if the geodesic (line) is tangential to the circle.
298 @raise IntersectionError: The circle and geodesic do not intersect.
300 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
301 or B{C{other}} invalid.
303 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{exact}} or
304 B{C{height}}.
306 @see: Method L{rhumbIntersecant2<LatLonBase.rhumbIntersecant2>}.
307 '''
308 try:
309 g, gl, p = self._g_gl_p3(point, other, exact, wrap)
310 r = Radius_(circle=circle) if _isRadius(circle) else \
311 g._Inverse(self, self.others(circle=circle), wrap).s12
313 P, Q = _MODS.geodesicw._Intersecant2(gl, self.lat, self.lon, r, tol=tol,
314 form=_MODS.dms.F_DMS)
315 return self._intersecend2(p, other, wrap, height, g, P, Q,
316 self.intersecant2)
317 except (TypeError, ValueError) as x:
318 raise _xError(x, center=self, circle=circle, point=point, other=other,
319 exact=exact, wrap=wrap)
321 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty}
322 '''(INTERNAL) I{Karney}'s C{Inverse} method.
324 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
325 '''
326 _ = self.ellipsoids(other)
327 g = self.geodesic
328 _, lon = unroll180(self.lon, other.lon, wrap=wrap)
329 return g.Inverse3(self.lat, self.lon, other.lat, lon)
331 def nearestOn8(self, points, closed=False, height=None, wrap=False,
332 equidistant=None, tol=_TOL_M):
333 '''I{Iteratively} locate the point on a path or polygon closest
334 to this point.
336 @arg points: The path or polygon points (C{LatLon}[]).
337 @kwarg closed: Optionally, close the polygon (C{bool}).
338 @kwarg height: Optional height, overriding the height of this and
339 all other points (C{meter}, conventionally). If
340 B{C{height}} is C{None}, the height of each point
341 is taken into account for distances.
342 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
343 B{C{points}} (C{bool}).
345 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start,
346 end, initial, final)} with C{distance} in C{meter},
347 conventionally and with the C{closest}, the C{start}
348 the C{end} point each an instance of this C{LatLon}.
350 @raise PointsError: Insufficient number of B{C{points}}.
352 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid.
354 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible
355 or no convergence for the given B{C{tol}}.
357 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}.
358 '''
359 _d3 = self.distanceTo3 # Distance3Tuple
360 _n3 = _nearestOn3
361 try:
362 Ps = self.PointsIter(points, loop=1, wrap=wrap)
363 p1 = c = s = e = Ps[0]
364 _ = self.ellipsoids(p1)
365 c3 = _d3(c, wrap=wrap) # XXX wrap=False?
367 except (TypeError, ValueError) as x:
368 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol,
369 closed=closed, height=height, wrap=wrap)
371 # get the azimuthal equidistant projection, once
372 A = _Equidistant00(equidistant, c)
373 b = _Box(c, c3.distance)
374 m = f = i = 0 # p1..p2 == points[i]..[j]
376 kwds = dict(within=True, height=height, tol=tol,
377 LatLon=self.classof, # this LatLon
378 datum=self.datum, epoch=self.epoch, reframe=self.reframe)
379 try:
380 for j, p2 in Ps.enumerate(closed=closed):
381 if wrap and j != 0: # not Ps.looped
382 p2 = _unrollon(p1, p2)
383 # skip edge if no overlap with box around closest
384 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon):
385 p, t, _ = _n3(self, p1, p2, A, **kwds)
386 d3 = _d3(p, wrap=False) # already unrolled
387 if d3.distance < c3.distance:
388 c3, c, s, e, f = d3, p, p1, p2, (i + t)
389 b = _Box(c, c3.distance)
390 m = max(m, c.iteration)
391 p1, i = p2, j
393 except (TypeError, ValueError) as x:
394 raise _xError(x, Fmt.SQUARE(points=i), p1,
395 Fmt.SQUARE(points=j), p2, this=self, tol=tol,
396 closed=closed, height=height, wrap=wrap)
398 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6
400 n = self.nearestOn8.__name__
401 c.rename(n)
402 if s is not c:
403 s = s.copy(name=n)
404 if e is not c:
405 e = e.copy(name=n)
406 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final,
407 iteration=m) # ._iteration for tests
409 def plumbTo(self, point, other, exact=False, height=None, # PYCHOK signature
410 wrap=False, tol=_TOL):
411 '''Compute the I{perpendicular} intersection of a geodesic (line) with
412 a geodesic from this point.
414 @arg point: A location (C{LatLon}) on the geodesic (line).
415 @arg other: An other point (C{LatLon}) on the geodesic (line) or the
416 (forward) bearing at the B{C{point}} (compass C{degrees}).
417 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}),
418 see method L{Ellipsoid.geodesic_}.
419 @kwarg height: Optional height for the intersection point (C{meter},
420 conventionally) or C{None} for an interpolated height.
421 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
422 and/or B{C{other}} (C{bool}).
423 @kwarg tol: Convergence tolerance (C{scalar}).
425 @return: The intersection point, an instance of this C{LatLon}.
427 @raise TypeError: If B{C{point}} or B{C{other}} not this C{LatLon}.
429 @raise UnitError: Invalid B{C{other}}, B{C{exact}} or B{C{height}}.
430 '''
431 try:
432 g, gl, p = self._g_gl_p3(point, other, exact, wrap)
434 P = _MODS.geodesicw._PlumbTo(gl, self.lat, self.lon, tol=tol)
435 h = self._havg(p, h=height)
436 p = self.classof(P.lat2, P.lon2, datum=self.datum, height=h) # name=n
437 p._iteration = P.iteration
438 except (TypeError, ValueError) as x:
439 raise _xError(x, plumb=self, point=point, other=other,
440 exact=exact, wrap=wrap)
441 return p
444class _Box(object):
445 '''Bounding box around a C{LatLon} point.
447 @see: Function C{_box4} in .clipy.py.
448 '''
449 def __init__(self, center, distance):
450 '''New L{_Box} around point.
452 @arg center: The center point (C{LatLon}).
453 @arg distance: Radius, half-size of the box
454 (C{meter}, conventionally)
455 '''
456 m = Radius_(distance=distance)
457 E = center.ellipsoid() # XXX see above
458 d = E.m2degrees(m) + _0_5 # some margin
459 self._N = center.lat + d
460 self._S = center.lat - d
461 self._E = center.lon + d
462 self._W = center.lon - d
464 def overlaps(self, lat1, lon1, lat2, lon2):
465 '''Check whether this box overlaps an other box.
467 @arg lat1: Latitude of a box corner (C{degrees}).
468 @arg lon1: Longitude of a box corner (C{degrees}).
469 @arg lat2: Latitude of the opposing corner (C{degrees}).
470 @arg lon2: Longitude of the opposing corner (C{degrees}).
472 @return: C{True} if there is some overlap, C{False}
473 otherwise (C{bool}).
474 '''
475 if lat1 > lat2:
476 lat1, lat2 = lat2, lat1
477 if lat1 > self._N or lat2 < self._S:
478 return False
479 if lon1 > lon2:
480 lon1, lon2 = lon2, lon1
481 if lon1 > self._E or lon2 < self._W:
482 return False
483 return True
486class _Tol(object):
487 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}.
488 '''
489 _deg = 0 # tol in degrees
490 _lat = 0
491 _m = 0 # tol in meter
492 _min = MAX # degrees
493 _prev = None
494 _r = 0
496 def __init__(self, tol_m, E, lat, *lats):
497 '''New L{_Tol}.
499 @arg tol_m: Tolerance (C{meter}, only).
500 @arg E: Earth ellipsoid (L{Ellipsoid}).
501 @arg lat: Latitude (C{degrees}).
502 @arg lats: Additional latitudes (C{degrees}).
503 '''
504 self._lat = fmean_(lat, *lats) if lats else lat
505 self._r = max(EPS, E.rocMean(self._lat))
506 self._m = max(EPS, tol_m)
507 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!
509 @property_RO
510 def degrees(self):
511 '''Get this tolerance in C{degrees}.
512 '''
513 return self._deg
515 def degrees2m(self, deg):
516 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.
517 '''
518 return self.radius * radians(deg) / PI2 # avoid degrees2m!
520 def degError(self, Error=_ValueError):
521 '''Compose an error with the C{deg}rees minimum.
522 '''
523 return self.mError(self.degrees2m(self._min), Error=Error)
525 def done(self, deg):
526 '''Check C{deg} vs. tolerance and previous value.
527 '''
528 if deg < self._deg or deg == self._prev:
529 return True
530 self._min = min(self._min, deg)
531 self._prev = deg
532 return False
534 @property_RO
535 def lat(self):
536 '''Get the mean latitude in C{degrees}.
537 '''
538 return self._lat
540 def mError(self, m, Error=_ValueError):
541 '''Compose an error with B{C{m}}eter minimum.
542 '''
543 t = _SPACE_(Fmt.tolerance(self.meter), _too_(_low_))
544 if m2km(m) > self.meter:
545 t = _or(t, _antipodal_, _near_(_polar__))
546 return Error(Fmt.no_convergence(m), txt=t)
548 @property_RO
549 def meter(self):
550 '''Get this tolerance in C{meter}.
551 '''
552 return self._m
554 @property_RO
555 def radius(self):
556 '''Get the earth radius in C{meter}.
557 '''
558 return self._r
560 def reset(self):
561 '''Reset tolerances.
562 '''
563 self._min = MAX
564 self._prev = None
567def _Equidistant00(equidistant, p1):
568 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
569 '''
570 if equidistant is None or not callable(equidistant):
571 equidistant = p1.Equidistant
572 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover
573 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants)
574 raise _IsnotError(*t, equidistant=equidistant)
575 return equidistant(0, 0, p1.datum)
578def intersecant2(center, circle, point, other, **exact_height_wrap_tol):
579 '''Compute the intersections of a circle and a geodesic given as two points
580 or as a point and (forward) bearing.
582 @arg center: Center of the circle (C{LatLon}).
583 @arg circle: Radius of the circle (C{meter}, conventionally) or a point on
584 the circle (C{LatLon}, as B{C{center}}).
585 @arg point: A point of the geodesic (C{LatLon}, as B{C{center}}).
586 @arg other: An other point of the geodesic (C{LatLon}, as B{C{center}}) or
587 the (forward) bearing at the B{C{point}} (compass C{degrees}).
588 @kwarg exact_height_wrap_tol: Optional keyword arguments, see below.
590 @raise NotImplementedError: Method C{intersecant2} not available.
592 @raise TypeError: If B{C{center}}, B{C{point}} or B{C{circle}} or B{C{other}}
593 points not ellipsoidal or not compatible with B{C{center}}.
595 @see: Method C{LatLon.intersecant2} of class L{ellipsoidalExact.LatLon},
596 L{ellipsoidalKarney.LatLon} or L{ellipsoidalVincenty.LatLon}.
597 '''
598 if not isLatLon(center, ellipsoidal=True): # isinstance(center, LatLonEllipsoidalBase)
599 raise _IsnotError(_ellipsoidal_, center=center)
600 return center.intersecant2(circle, point, other, **exact_height_wrap_tol)
603def _intersect3(s1, end1, s2, end2, height=None, wrap=False, # MCCABE 16 was=True
604 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
605 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method
606 L{intersection3}, separated to allow callers to embellish any exceptions.
607 '''
608 _opp = _MODS.formy.opposing
609 _LLS = _MODS.sphericalTrigonometry.LatLon
610 _si = _MODS.sphericalTrigonometry._intersect
611 _vi3 = _MODS.vector3d._intersect3d3
613 def _b_d(s, e, w, t, h=_0_0):
614 # compute opposing and distance
615 t = s.classof(t.lat, t.lon, height=h, name=t.name)
616 t = s.distanceTo3(t, wrap=w) # Distance3Tuple
617 b = _opp(e, t.initial) # "before" start
618 return b, t.distance
620 def _b_e(s, e, w, t):
621 # compute an end point along the initial bearing
622 # about 1.5 times the distance to the gu-/estimate, at
623 # least 1/8 and at most 3/8 of the earth perimeter like
624 # radians in .sphericalTrigonometry._int3d2 and bearing
625 # comparison in .sphericalTrigonometry._intb
626 b, d = _b_d(s, e, w, t, h=t.height)
627 m = s.ellipsoid().R2x * PI_4 # authalic exact
628 d = min(max(d * _B2END, m), m * _3_0)
629 e = s.destination(d, e)
630 return b, (_unrollon(s, e) if w else e)
632 def _e_ll(s, e, w, **end):
633 # return 2-tuple (end, False if bearing else True)
634 ll = not _isDegrees(e)
635 if ll:
636 e = s.others(**end)
637 if w: # unroll180 == .karney._unroll2
638 e = _unrollon(s, e)
639 return e, ll
641 def _llS(ll): # return an C{_LLS} instance
642 return _LLS(ll.lat, ll.lon, height=ll.height)
644 def _o(o, b, n, s, t, e):
645 # determine C{o}utside before, on or after start point
646 if not o: # intersection may be on start
647 if _MODS.formy._isequalTo(s, t, eps=e.degrees):
648 return o
649 return -n if b else n
651 E = s1.ellipsoids(s2)
653 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
654 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
656 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),
657 s2.lat, (e2.lat if ll2 else s2.lat))
659 # get the azimuthal equidistant projection
660 A = _Equidistant00(equidistant, s1)
662 # gu-/estimate initial intersection, spherically ...
663 t = _si(_llS(s1), (_llS(e1) if ll1 else e1),
664 _llS(s2), (_llS(e2) if ll2 else e2),
665 height=height, wrap=False, LatLon=_LLS) # unrolled already
666 h, n = t.height, t.name
668 if not ll1:
669 b1, e1 = _b_e(s1, e1, wrap, t)
670 if not ll2:
671 b2, e2 = _b_e(s2, e2, wrap, t)
673 # ... and iterate as Karney describes, for references
674 # @see: Function L{ellipsoidalKarney.intersection3}.
675 for i in range(1, _TRIPS):
676 A.reset(t.lat, t.lon) # gu-/estimate as origin
677 # convert start and end points to projection
678 # space and compute an intersection there
679 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),
680 eps=e.meter, useZ=False)
681 # convert intersection back to geodetic
682 t, d = A._reverse2(v)
683 if e.done(d): # below tol or unchanged?
684 t._iteration = i
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=t._iteration, 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 _latlon4(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 def _llS(ll): # return an C{_LLS} instance
750 return _LLS(ll.lat, ll.lon, height=ll.height)
752 r1 = Radius_(radius1=radius1)
753 r2 = Radius_(radius2=radius2)
755 E = c1.ellipsoids(c2)
756 # get the azimuthal equidistant projection
757 A = _Equidistant00(equidistant, c1)
759 if r1 < r2:
760 c1, c2 = c2, c1
761 r1, r2 = r2, r1
763 if r1 > (min(E.b, E.a) * PI):
764 raise _ValueError(_exceed_PI_radians_)
766 if wrap: # unroll180 == .karney._unroll2
767 c2 = _unrollon(c1, c2)
769 # distance between centers and radii are
770 # measured along the ellipsoid's surface
771 m = c1.distanceTo(c2, wrap=False) # meter
772 if m < max(r1 - r2, EPS):
773 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
774 if fsumf_(r1, r2, -m) < 0:
775 raise IntersectionError(_too_(Fmt.distant(m)))
777 f = _MODS.formy._radical2(m, r1, r2).ratio # "radical fraction"
778 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))
780 # gu-/estimate initial intersections, spherically ...
781 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius,
782 height=height, wrap=False, too_d=m) # unrolled already
783 h, n = t1.height, t1.name
785 # ... and iterate as Karney describes, for references
786 # @see: Function L{ellipsoidalKarney.intersections2}.
787 ts, ta = [], None
788 for t in ((t1,) if t1 is t2 else (t1, t2)):
789 for i in range(1, _TRIPS):
790 A.reset(t.lat, t.lon) # gu-/estimate as origin
791 # convert centers to projection space and
792 # compute the intersections there
793 t1, t2 = A._forwards(c1, c2)
794 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,
795 t2, r2, # XXX * t2.scale?,
796 sphere=False, too_d=m)
797 # convert intersections back to geodetic
798 if v1 is v2: # abutting
799 t, d = A._reverse2(v1)
800 else: # consider the closer intersection
801 t1, d1 = A._reverse2(v1)
802 t2, d2 = A._reverse2(v2)
803 t, d = (t1, d1) if d1 < d2 else (t2, d2)
804 if e.done(d): # below tol or unchanged?
805 t._iteration = i # _NamedTuple._iteration
806 ts.append(t)
807 if v1 is v2: # abutting
808 ta = t
809 break
810 else:
811 raise e.degError(Error=IntersectionError)
812 e.reset()
814 if ta: # abutting circles
815 pass # PYCHOK no cover
816 elif len(ts) == 2:
817 return (_latlon4(ts[0], h, n, c1),
818 _latlon4(ts[1], h, n, c2))
819 elif len(ts) == 1: # PYCHOK no cover
820 ta = ts[0] # assume abutting
821 else: # PYCHOK no cover
822 raise _AssertionError(ts=ts)
823 r = _latlon4(ta, h, n, c1)
824 return r, r
827def _nearestOn2(p, point1, point2, within=True, height=None, wrap=False, # was=True
828 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):
829 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,
830 separated to allow callers to embellish any exceptions.
831 '''
832 p1 = p.others(point1=point1)
833 p2 = p.others(point2=point2)
835 _ = p.ellipsoids(p1)
836# E = p.ellipsoids(p2) # done in _nearestOn3
838 # get the azimuthal equidistant projection
839 A = _Equidistant00(equidistant, p)
841 p1, p2, _ = _unrollon3(p, p1, p2, wrap) # XXX don't unroll?
842 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height,
843 tol=tol, **LatLon_and_kwds)
844 return NearestOn2Tuple(r, f)
847def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M,
848 LatLon=None, **LatLon_kwds):
849 # Only in function C{_nearestOn2} and method C{nearestOn8} above
850 _LLS = _MODS.sphericalNvector.LatLon
851 _V3d = _MODS.vector3d.Vector3d
852 _vn2 = _MODS.vector3d._nearestOn2
854 def _llS(ll): # return an C{_LLS} instance
855 return _LLS(ll.lat, ll.lon, height=ll.height)
857 E = p.ellipsoids(p2)
858 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)
860 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
861 # using sphericalNvector.LatLon.nearestOn for within=False support
862 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within,
863 height=height)
864 n, h = t.name, t.height
865 if height is None:
866 h1 = p1.height # use heights as pseudo-Z in projection space
867 h2 = p2.height # to be included in the closest function
868 h0 = favg(h1, h2)
869 else: # ignore heights in distances, Z=0
870 h0 = h1 = h2 = _0_0
872 # ... and iterate to find the closest (to the origin with .z
873 # to interpolate height) as Karney describes, for references
874 # @see: Function L{ellipsoidalKarney.nearestOn}.
875 vp, f = _V3d(_0_0, _0_0, h0), None
876 for i in range(1, _TRIPS):
877 A.reset(t.lat, t.lon) # gu-/estimate as origin
878 # convert points to projection space and compute
879 # the nearest one (and its height) there
880 s, t = A._forwards(p1, p2)
881 v, f = _vn2(vp, _V3d(s.x, s.y, h1),
882 _V3d(t.x, t.y, h2), within=within)
883 # convert nearest one back to geodetic
884 t, d = A._reverse2(v)
885 if e.done(d): # below tol or unchanged?
886 t._iteration = i # _NamedTuple._iteration
887 break
888 else:
889 raise e.degError()
891 if height is None:
892 h = v.z # nearest
893 elif _isHeight(height):
894 h = height
895 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,
896 iteration=t.iteration, name=n)
897 return r, f, e # fraction or None
900__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2)
902# **) MIT License
903#
904# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
905#
906# Permission is hereby granted, free of charge, to any person obtaining a
907# copy of this software and associated documentation files (the "Software"),
908# to deal in the Software without restriction, including without limitation
909# the rights to use, copy, modify, merge, publish, distribute, sublicense,
910# and/or sell copies of the Software, and to permit persons to whom the
911# Software is furnished to do so, subject to the following conditions:
912#
913# The above copyright notice and this permission notice shall be included
914# in all copies or substantial portions of the Software.
915#
916# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
917# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
918# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
919# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
920# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
921# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
922# OTHER DEALINGS IN THE SOFTWARE.