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