Coverage for pygeodesy/ellipsoidalBaseDI.py: 93%
314 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -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 isscalar, issubclassof
11from pygeodesy.constants import EPS, MAX, PI, PI2, PI_4, isnear0, isnear1, \
12 _0_0, _0_5, _1_5, _3_0
13from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, Property_RO, \
14 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_, _exceed_PI_radians_, \
21 _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, 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.05.15'
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 initialBearingTo(self, other, wrap=False):
212 '''Compute the initial bearing (forward azimuth) to travel
213 along a geodesic from this point to an other point,
214 using this C{Inverse} method. See method L{distanceTo3}
215 for more details.
217 @arg other: The other point (C{LatLon}).
218 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
219 the B{C{other}} point (C{bool}).
221 @return: Initial bearing (compass C{degrees360}).
223 @raise TypeError: The B{C{other}} point is not C{LatLon}.
225 @raise ValueError: If this and the B{C{other}} point's L{Datum}
226 ellipsoids are not compatible.
227 '''
228 return self._Inverse(other, wrap).initial
230 def intermediateTo(self, other, fraction, height=None, wrap=False):
231 '''Return the point at given fraction along the geodesic between
232 this and an other point, using this C{Direct} and C{Inverse}
233 methods.
235 @arg other: The other point (C{LatLon}).
236 @arg fraction: Fraction between both points (C{scalar}, 0.0
237 at this and 1.0 at the other point.
238 @kwarg height: Optional height, overriding the fractional
239 height (C{meter}).
240 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
241 B{C{other}} point (C{bool}).
243 @return: Intermediate point (C{LatLon}).
245 @raise TypeError: The B{C{other}} point is not C{LatLon}.
247 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.
249 @raise ValueError: If this and the B{C{other}} point's L{Datum}
250 ellipsoids are not compatible.
252 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and
253 C{rhumbMidpointTo}.
254 '''
255 f = Scalar(fraction=fraction)
256 if isnear0(f):
257 r = self
258 elif isnear1(f) and not wrap:
259 r = self.others(other)
260 else: # negative fraction OK
261 t = self.distanceTo3(other, wrap=wrap)
262 h = self._havg(other, f=f, h=height)
263 r = self.destination(t.distance * f, t.initial, height=h)
264 return r
266 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty}
267 '''(INTERNAL) I{Karney}'s C{Inverse} method.
269 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
270 '''
271 _ = self.ellipsoids(other)
272 g = self.geodesic
273 _, lon = unroll180(self.lon, other.lon, wrap=wrap)
274 return g.Inverse3(self.lat, self.lon, other.lat, lon)
276 def nearestOn8(self, points, closed=False, height=None, wrap=False,
277 equidistant=None, tol=_TOL_M):
278 '''I{Iteratively} locate the point on a path or polygon closest
279 to this point.
281 @arg points: The path or polygon points (C{LatLon}[]).
282 @kwarg closed: Optionally, close the polygon (C{bool}).
283 @kwarg height: Optional height, overriding the height of this and
284 all other points (C{meter}, conventionally). If
285 B{C{height}} is C{None}, the height of each point
286 is taken into account for distances.
287 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
288 B{C{points}} (C{bool}).
290 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start,
291 end, initial, final)} with C{distance} in C{meter},
292 conventionally and with the C{closest}, the C{start}
293 the C{end} point each an instance of this C{LatLon}.
295 @raise PointsError: Insufficient number of B{C{points}}.
297 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid.
299 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible
300 or no convergence for the given B{C{tol}}.
302 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}.
303 '''
304 _d3 = self.distanceTo3 # Distance3Tuple
305 _n3 = _nearestOn3
306 try:
307 Ps = self.PointsIter(points, loop=1, wrap=wrap)
308 p1 = c = s = e = Ps[0]
309 _ = self.ellipsoids(p1)
310 c3 = _d3(c, wrap=wrap) # XXX wrap=False?
312 except (TypeError, ValueError) as x:
313 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol,
314 closed=closed, height=height, wrap=wrap)
316 # get the azimuthal equidistant projection, once
317 A = _Equidistant00(equidistant, c)
318 b = _Box(c, c3.distance)
319 m = f = i = 0 # p1..p2 == points[i]..[j]
321 kwds = dict(within=True, height=height, tol=tol,
322 LatLon=self.classof, # this LatLon
323 datum=self.datum, epoch=self.epoch, reframe=self.reframe)
324 try:
325 for j, p2 in Ps.enumerate(closed=closed):
326 if wrap and j != 0: # not Ps.looped
327 p2 = _unrollon(p1, p2)
328 # skip edge if no overlap with box around closest
329 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon):
330 p, t, _ = _n3(self, p1, p2, A, **kwds)
331 d3 = _d3(p, wrap=False) # already unrolled
332 if d3.distance < c3.distance:
333 c3, c, s, e, f = d3, p, p1, p2, (i + t)
334 b = _Box(c, c3.distance)
335 m = max(m, c.iteration)
336 p1, i = p2, j
338 except (TypeError, ValueError) as x:
339 raise _xError(x, Fmt.SQUARE(points=i), p1,
340 Fmt.SQUARE(points=j), p2, this=self, tol=tol,
341 closed=closed, height=height, wrap=wrap)
343 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6
345 n = self.nearestOn8.__name__
346 c.rename(n)
347 if s is not c:
348 s = s.copy(name=n)
349 if e is not c:
350 e = e.copy(name=n)
351 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final,
352 iteration=m) # ._iteration for tests
355class _Box(object):
356 '''Bounding box around a C{LatLon} point.
358 @see: Function C{_box4} in .clipy.py.
359 '''
360 def __init__(self, center, distance):
361 '''New L{_Box} around point.
363 @arg center: The center point (C{LatLon}).
364 @arg distance: Radius, half-size of the box
365 (C{meter}, conventionally)
366 '''
367 E = center.ellipsoid()
368 d = degrees(distance / max(E.a, E.b)) + _0_5 # some margin
369 self._N = center.lat + d
370 self._S = center.lat - d
371 self._E = center.lon + d
372 self._W = center.lon - d
374 def overlaps(self, lat1, lon1, lat2, lon2):
375 '''Check whether this box overlaps a line between 2 points.
377 @arg lat1: Latitude of first point (C{degrees}).
378 @arg lon1: Longitude of first point (C{degrees}).
379 @arg lat2: Latitude of second point (C{degrees}).
380 @arg lon2: Longitude of second point (C{degrees}).
382 @return: C{False} if there is certainly no overlap,
383 C{True} otherwise (C{bool}).
384 '''
385 if lat1 > lat2:
386 lat1, lat2 = lat2, lat1
387 if lat1 > self._N or lat2 < self._S:
388 return False
389 if lon1 > lon2:
390 lon1, lon2 = lon2, lon1
391 if lon1 > self._E or lon2 < self._W:
392 return False
393 return True
396class _Tol(object):
397 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}.
398 '''
399 _deg = 0 # tol in degrees
400 _lat = 0
401 _m = 0 # tol in meter
402 _min = MAX # degrees
403 _prev = None
404 _r = 0
406 def __init__(self, tol_m, E, lat, *lats):
407 '''New L{_Tol}.
409 @arg tol_m: Tolerance (C{meter}, only).
410 @arg E: Earth ellipsoid (L{Ellipsoid}).
411 @arg lat: Latitude (C{degrees}).
412 @arg lats: Additional latitudes (C{degrees}).
413 '''
414 self._lat = fmean_(lat, *lats) if lats else lat
415 self._r = max(EPS, E.rocMean(self._lat))
416 self._m = max(EPS, tol_m)
417 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!
419 @property_RO
420 def degrees(self):
421 '''Get this tolerance in C{degrees}.
422 '''
423 return self._deg
425 def degrees2m(self, deg):
426 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.
427 '''
428 return self.radius * radians(deg) / PI2 # avoid degrees2m!
430 def degError(self, Error=_ValueError):
431 '''Compose an error with the C{deg}rees minimum.
432 '''
433 return self.mError(self.degrees2m(self._min), Error=Error)
435 def done(self, deg):
436 '''Check C{deg} vs. tolerance and previous value.
437 '''
438 if deg < self._deg or deg == self._prev:
439 return True
440 self._min = min(self._min, deg)
441 self._prev = deg
442 return False
444 @property_RO
445 def lat(self):
446 '''Get the mean latitude in C{degrees}.
447 '''
448 return self._lat
450 def mError(self, m, Error=_ValueError):
451 '''Compose an error with B{C{m}}eter minimum.
452 '''
453 t = _SPACE_(Fmt.tolerance(self.meter), _too_(_low_))
454 if m2km(m) > self.meter:
455 t = _or(t, _antipodal_, _near_(_polar__))
456 return Error(Fmt.no_convergence(m), txt=t)
458 @property_RO
459 def meter(self):
460 '''Get this tolerance in C{meter}.
461 '''
462 return self._m
464 @property_RO
465 def radius(self):
466 '''Get the earth radius in C{meter}.
467 '''
468 return self._r
470 def reset(self):
471 '''Reset tolerances.
472 '''
473 self._min = MAX
474 self._prev = None
477def _Equidistant00(equidistant, p1):
478 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
479 '''
480 if equidistant is None or not callable(equidistant):
481 equidistant = p1.Equidistant
482 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover
483 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants)
484 raise _IsnotError(*t, equidistant=equidistant)
485 return equidistant(0, 0, p1.datum)
488def _intersect3(s1, end1, s2, end2, height=None, wrap=False, # MCCABE 16 was=True
489 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
490 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method
491 L{intersection3}, separated to allow callers to embellish any exceptions.
492 '''
493 _LLS = _MODS.sphericalTrigonometry.LatLon
494 _si = _MODS.sphericalTrigonometry._intersect
495 _vi3 = _MODS.vector3d._intersect3d3
497 def _b_d(s, e, w, t, h=_0_0):
498 # compute opposing and distance
499 t = s.classof(t.lat, t.lon, height=h, name=t.name)
500 t = s.distanceTo3(t, wrap=w) # Distance3Tuple
501 b = opposing(e, t.initial) # "before" start
502 return b, t.distance
504 def _b_e(s, e, w, t):
505 # compute an end point along the initial bearing
506 # about 1.5 times the distance to the gu-/estimate, at
507 # least 1/8 and at most 3/8 of the earth perimeter like
508 # radians in .sphericalTrigonometry._int3d2 and bearing
509 # comparison in .sphericalTrigonometry._intb
510 b, d = _b_d(s, e, w, t, h=t.height)
511 m = s.ellipsoid().R2x * PI_4 # authalic exact
512 d = min(max(d * _B2END, m), m * _3_0)
513 e = s.destination(d, e)
514 return b, (_unrollon(s, e) if w else e)
516 def _e_ll(s, e, w, **end):
517 # return 2-tuple (end, False if bearing else True)
518 ll = not isscalar(e)
519 if ll:
520 e = s.others(**end)
521 if w: # unroll180 == .karney._unroll2
522 e = _unrollon(s, e)
523 return e, ll
525 def _llS(ll): # return an C{_LLS} instance
526 return _LLS(ll.lat, ll.lon, height=ll.height)
528 def _o(o, b, n, s, t, e):
529 # determine C{o}utside before, on or after start point
530 if not o: # intersection may be on start
531 if _MODS.latlonBase._isequalTo(s, t, eps=e.degrees):
532 return o
533 return -n if b else n
535 E = s1.ellipsoids(s2)
537 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
538 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
540 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),
541 s2.lat, (e2.lat if ll2 else s2.lat))
543 # get the azimuthal equidistant projection
544 A = _Equidistant00(equidistant, s1)
546 # gu-/estimate initial intersection, spherically ...
547 t = _si(_llS(s1), (_llS(e1) if ll1 else e1),
548 _llS(s2), (_llS(e2) if ll2 else e2),
549 height=height, wrap=False, LatLon=_LLS) # unrolled already
550 h, n = t.height, t.name
552 if not ll1:
553 b1, e1 = _b_e(s1, e1, wrap, t)
554 if not ll2:
555 b2, e2 = _b_e(s2, e2, wrap, t)
557 # ... and iterate as Karney describes, for references
558 # @see: Function L{ellipsoidalKarney.intersection3}.
559 for i in range(1, _TRIPS):
560 A.reset(t.lat, t.lon) # gu-/estimate as origin
561 # convert start and end points to projection
562 # space and compute an intersection there
563 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),
564 eps=e.meter, useZ=False)
565 # convert intersection back to geodetic
566 t, d = A._reverse2(v)
567 if e.done(d): # below tol or unchanged?
568 t._iteration = i
569 break
570 else:
571 raise e.degError(Error=IntersectionError)
573 # like .sphericalTrigonometry._intersect, if this intersection
574 # is "before" the first point, use the antipodal intersection
575 if not (ll1 or ll2): # end1 and end2 are an initial bearing
576 b1, _ = _b_d(s1, end1, wrap, t)
577 if b1:
578 t = t.antipodal()
579 b1 = not b1
580 b2, _ = _b_d(s2, end2, wrap, t)
582 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1,
583 iteration=t._iteration, name=n)
584 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),
585 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))
588def _intersection3(start1, end1, start2, end2, height=None, wrap=False, # was=True
589 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
590 '''(INTERNAL) Iteratively compute the intersection point of two lines,
591 each defined by two (ellipsoidal) points or an (ellipsoidal) start
592 point and an initial bearing from North.
593 '''
594 s1 = _xellipsoidal(start1=start1)
595 s2 = s1.others(start2=start2)
596 try:
597 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,
598 equidistant=equidistant, tol=tol,
599 LatLon=LatLon, **LatLon_kwds)
600 except (TypeError, ValueError) as x:
601 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
604def _intersections2(center1, radius1, center2, radius2, height=None, wrap=False, # was=True
605 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
606 '''(INTERNAL) Iteratively compute the intersection points of two circles,
607 each defined by an (ellipsoidal) center point and a radius.
608 '''
609 c1 = _xellipsoidal(center1=center1)
610 c2 = c1.others(center2=center2)
611 try:
612 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,
613 equidistant=equidistant, tol=tol,
614 LatLon=LatLon, **LatLon_kwds)
615 except (TypeError, ValueError) as x:
616 raise _xError(x, center1=center1, radius1=radius1,
617 center2=center2, radius2=radius2)
620def _intersects2(c1, radius1, c2, radius2, height=None, wrap=False, # MCCABE 16 was=True
621 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
622 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}
623 above, separated to allow callers to embellish any exceptions.
624 '''
625 _LLS = _MODS.sphericalTrigonometry.LatLon
626 _si2 = _MODS.sphericalTrigonometry._intersects2
627 _vi2 = _MODS.vector3d._intersects2
629 def _latlon4(t, h, n, c):
630 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c,
631 iteration=t.iteration, name=n)
633 def _llS(ll): # return an C{_LLS} instance
634 return _LLS(ll.lat, ll.lon, height=ll.height)
636 r1 = Radius_(radius1=radius1)
637 r2 = Radius_(radius2=radius2)
639 E = c1.ellipsoids(c2)
640 # get the azimuthal equidistant projection
641 A = _Equidistant00(equidistant, c1)
643 if r1 < r2:
644 c1, c2 = c2, c1
645 r1, r2 = r2, r1
647 if r1 > (min(E.b, E.a) * PI):
648 raise _ValueError(_exceed_PI_radians_)
650 if wrap: # unroll180 == .karney._unroll2
651 c2 = _unrollon(c1, c2)
653 # distance between centers and radii are
654 # measured along the ellipsoid's surface
655 m = c1.distanceTo(c2, wrap=False) # meter
656 if m < max(r1 - r2, EPS):
657 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
658 if fsumf_(r1, r2, -m) < 0:
659 raise IntersectionError(_too_(Fmt.distant(m)))
661 f = _radical2(m, r1, r2).ratio # "radical fraction"
662 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))
664 # gu-/estimate initial intersections, spherically ...
665 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius,
666 height=height, wrap=False, too_d=m) # unrolled already
667 h, n = t1.height, t1.name
669 # ... and iterate as Karney describes, for references
670 # @see: Function L{ellipsoidalKarney.intersections2}.
671 ts, ta = [], None
672 for t in ((t1,) if t1 is t2 else (t1, t2)):
673 for i in range(1, _TRIPS):
674 A.reset(t.lat, t.lon) # gu-/estimate as origin
675 # convert centers to projection space and
676 # compute the intersections there
677 t1, t2 = A._forwards(c1, c2)
678 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,
679 t2, r2, # XXX * t2.scale?,
680 sphere=False, too_d=m)
681 # convert intersections back to geodetic
682 if v1 is v2: # abutting
683 t, d = A._reverse2(v1)
684 else: # consider the closer intersection
685 t1, d1 = A._reverse2(v1)
686 t2, d2 = A._reverse2(v2)
687 t, d = (t1, d1) if d1 < d2 else (t2, d2)
688 if e.done(d): # below tol or unchanged?
689 t._iteration = i # _NamedTuple._iteration
690 ts.append(t)
691 if v1 is v2: # abutting
692 ta = t
693 break
694 else:
695 raise e.degError(Error=IntersectionError)
696 e.reset()
698 if ta: # abutting circles
699 pass # PYCHOK no cover
700 elif len(ts) == 2:
701 return (_latlon4(ts[0], h, n, c1),
702 _latlon4(ts[1], h, n, c2))
703 elif len(ts) == 1: # PYCHOK no cover
704 ta = ts[0] # assume abutting
705 else: # PYCHOK no cover
706 raise _AssertionError(ts=ts)
707 r = _latlon4(ta, h, n, c1)
708 return r, r
711def _nearestOn2(p, point1, point2, within=True, height=None, wrap=False, # was=True
712 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):
713 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,
714 separated to allow callers to embellish any exceptions.
715 '''
716 p1 = p.others(point1=point1)
717 p2 = p.others(point2=point2)
719 _ = p.ellipsoids(p1)
720# E = p.ellipsoids(p2) # done in _nearestOn3
722 # get the azimuthal equidistant projection
723 A = _Equidistant00(equidistant, p)
725 p1, p2, _ = _unrollon3(p, p1, p2, wrap) # XXX don't unroll?
726 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height,
727 tol=tol, **LatLon_and_kwds)
728 return NearestOn2Tuple(r, f)
731def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M,
732 LatLon=None, **LatLon_kwds):
733 # Only in function C{_nearestOn2} and method C{nearestOn8} above
734 _LLS = _MODS.sphericalNvector.LatLon
735 _V3d = _MODS.vector3d.Vector3d
736 _vn2 = _MODS.vector3d._nearestOn2
738 def _llS(ll): # return an C{_LLS} instance
739 return _LLS(ll.lat, ll.lon, height=ll.height)
741 E = p.ellipsoids(p2)
742 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)
744 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
745 # using sphericalNvector.LatLon.nearestOn for within=False support
746 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within,
747 height=height)
748 n, h = t.name, t.height
749 if height is None:
750 h1 = p1.height # use heights as pseudo-Z in projection space
751 h2 = p2.height # to be included in the closest function
752 h0 = favg(h1, h2)
753 else: # ignore heights in distances, Z=0
754 h0 = h1 = h2 = _0_0
756 # ... and iterate to find the closest (to the origin with .z
757 # to interpolate height) as Karney describes, for references
758 # @see: Function L{ellipsoidalKarney.nearestOn}.
759 vp, f = _V3d(_0_0, _0_0, h0), None
760 for i in range(1, _TRIPS):
761 A.reset(t.lat, t.lon) # gu-/estimate as origin
762 # convert points to projection space and compute
763 # the nearest one (and its height) there
764 s, t = A._forwards(p1, p2)
765 v, f = _vn2(vp, _V3d(s.x, s.y, h1),
766 _V3d(t.x, t.y, h2), within=within)
767 # convert nearest one back to geodetic
768 t, d = A._reverse2(v)
769 if e.done(d): # below tol or unchanged?
770 t._iteration = i # _NamedTuple._iteration
771 break
772 else:
773 raise e.degError()
775 if height is None:
776 h = v.z # nearest
777 elif isscalar(height):
778 h = height
779 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,
780 iteration=t.iteration, name=n)
781 return r, f, e # fraction or None
784__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)
786# **) MIT License
787#
788# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
789#
790# Permission is hereby granted, free of charge, to any person obtaining a
791# copy of this software and associated documentation files (the "Software"),
792# to deal in the Software without restriction, including without limitation
793# the rights to use, copy, modify, merge, publish, distribute, sublicense,
794# and/or sell copies of the Software, and to permit persons to whom the
795# Software is furnished to do so, subject to the following conditions:
796#
797# The above copyright notice and this permission notice shall be included
798# in all copies or substantial portions of the Software.
799#
800# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
801# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
802# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
803# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
804# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
805# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
806# OTHER DEALINGS IN THE SOFTWARE.