Coverage for pygeodesy/ellipsoidalBaseDI.py: 96%
312 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -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, fsum_
19from pygeodesy.formy import opposing, _radical2
20from pygeodesy.interns import _antipodal_, _concentric_, _exceed_PI_radians_, \
21 _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, Height, Radius_, Scalar
29from pygeodesy.utily import m2km, unroll180, _unrollon, wrap90, wrap180, wrap360
31from math import degrees, radians
33__all__ = _ALL_LAZY.ellipsoidalBaseDI
34__version__ = '23.04.11'
36_polar__ = 'polar?'
37_too_low_ = _too_('low')
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: Wrap and unroll longitudes (C{bool}).
56 @return: A L{Bearing2Tuple}C{(initial, final)}.
58 @raise TypeError: The B{C{other}} point is not C{LatLon}.
60 @raise ValueError: If this and the B{C{other}} point's L{Datum}
61 ellipsoids are not compatible.
62 '''
63 r = self._Inverse(other, wrap)
64 return Bearing2Tuple(r.initial, r.final, name=self.name)
66 def destination(self, distance, bearing, height=None):
67 '''Compute the destination point after having travelled for
68 the given distance from this point along a geodesic given
69 by an initial bearing, using this C{Direct} method. See
70 method L{destination2} for more details.
72 @arg distance: Distance (C{meter}).
73 @arg bearing: Initial bearing in (compass C{degrees360}).
74 @kwarg height: Optional height, overriding the default
75 height (C{meter}, same units as C{distance}).
77 @return: The destination point (C{LatLon}).
78 '''
79 return self._Direct(distance, bearing, self.classof, height).destination
81 def destination2(self, distance, bearing, height=None):
82 '''Compute the destination point and the final bearing (reverse
83 azimuth) after having travelled for the given distance from
84 this point along a geodesic given by an initial bearing,
85 using this C{Direct} method.
87 The distance must be in the same units as this point's datum
88 axes, conventionally C{meter}. The distance is measured on
89 the surface of the ellipsoid, ignoring this point's height.
91 The initial and final bearing (forward and reverse azimuth)
92 are in compass C{degrees360}.
94 The destination point's height and datum are set to this
95 point's height and datum, unless the former is overridden.
97 @arg distance: Distance (C{meter}).
98 @arg bearing: Initial bearing (compass C{degrees360}).
99 @kwarg height: Optional height, overriding the default
100 height (C{meter}, same units as C{distance}).
102 @return: A L{Destination2Tuple}C{(destination, final)}.
103 '''
104 r = self._Direct(distance, bearing, self.classof, height)
105 return self._xnamed(r)
107 def _Direct(self, distance, bearing, LL, height): # overloaded by I{Vincenty}
108 '''(INTERNAL) I{Karney}'s C{Direct} method.
110 @return: A L{Destination2Tuple}C{(destination, final)} or
111 a L{Destination3Tuple}C{(lat, lon, final)} if
112 B{C{LL}} is C{None}.
113 '''
114 g = self.geodesic
115 r = g.Direct3(self.lat, self.lon, bearing, distance)
116 if LL:
117 r = self._Direct2Tuple(LL, height, r)
118 return r
120 def _Direct2Tuple(self, LL, height, r):
121 '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}.
122 '''
123 h = self.height if height is None else height
124 d = LL(wrap90(r.lat), wrap180(r.lon), height=h, datum=self.datum, name=self.name,
125 **_xkwds_not(None, epoch=self.epoch, reframe=self.reframe))
126 return Destination2Tuple(d, wrap360(r.final))
128 def distanceTo(self, other, wrap=False, **unused): # ignore radius=R_M
129 '''Compute the distance between this and an other point
130 along a geodesic, using this C{Inverse} method. See method
131 L{distanceTo3} for more details.
133 @arg other: The other point (C{LatLon}).
134 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
136 @return: Distance (C{meter}).
138 @raise TypeError: The B{C{other}} point is not C{LatLon}.
140 @raise ValueError: If this and the B{C{other}} point's L{Datum}
141 ellipsoids are not compatible.
142 '''
143 return self._Inverse(other, wrap, azis=False).distance
145 def distanceTo3(self, other, wrap=False):
146 '''Compute the distance, the initial and final bearing along
147 a geodesic between this and an other point, using this
148 C{Inverse} method.
150 The distance is in the same units as this point's datum axes,
151 conventionally meter. The distance is measured on the surface
152 of the ellipsoid, ignoring this point's height.
154 The initial and final bearing (forward and reverse azimuth)
155 are in compass C{degrees360} from North.
157 @arg other: Destination point (C{LatLon}).
158 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
160 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
162 @raise TypeError: The B{C{other}} point is not C{LatLon}.
164 @raise ValueError: If this and the B{C{other}} point's L{Datum}
165 ellipsoids are not compatible.
166 '''
167 return self._xnamed(self._Inverse(other, wrap))
169 def finalBearingOn(self, distance, bearing):
170 '''Compute the final bearing (reverse azimuth) after having
171 travelled for the given distance along a geodesic given by
172 an initial bearing from this point, using this C{Direct}
173 method. See method L{destination2} for more details.
175 @arg distance: Distance (C{meter}).
176 @arg bearing: Initial bearing (compass C{degrees360}).
178 @return: Final bearing (compass C{degrees360}).
179 '''
180 return self._Direct(distance, bearing, None, None).final
182 def finalBearingTo(self, other, wrap=False):
183 '''Compute the final bearing (reverse azimuth) after having
184 travelled along a geodesic from this point to an other
185 point, using this C{Inverse} method. See method
186 L{distanceTo3} for more details.
188 @arg other: The other point (C{LatLon}).
189 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
191 @return: Final bearing (compass C{degrees360}).
193 @raise TypeError: The B{C{other}} point is not C{LatLon}.
195 @raise ValueError: If this and the B{C{other}} point's L{Datum}
196 ellipsoids are not compatible.
197 '''
198 return self._Inverse(other, wrap).final
200 @Property_RO
201 def geodesic(self): # overloaded by I{Karney}'s, N/A for I{Vincenty}
202 '''N/A, invalid (C{None} I{always}).
203 '''
204 return None # PYCHOK no cover
206 def initialBearingTo(self, other, wrap=False):
207 '''Compute the initial bearing (forward azimuth) to travel
208 along a geodesic from this point to an other point,
209 using this C{Inverse} method. See method L{distanceTo3}
210 for more details.
212 @arg other: The other point (C{LatLon}).
213 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
215 @return: Initial bearing (compass C{degrees360}).
217 @raise TypeError: The B{C{other}} point is not C{LatLon}.
219 @raise ValueError: If this and the B{C{other}} point's L{Datum}
220 ellipsoids are not compatible.
221 '''
222 return self._Inverse(other, wrap).initial
224 def intermediateTo(self, other, fraction, height=None, wrap=False):
225 '''Return the point at given fraction along the geodesic between
226 this and an other point, using this C{Direct} and C{Inverse}
227 methods.
229 @arg other: The other point (C{LatLon}).
230 @arg fraction: Fraction between both points (C{scalar},
231 0.0 at this and 1.0 at the other point.
232 @kwarg height: Optional height, overriding the fractional
233 height (C{meter}).
234 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
236 @return: Intermediate point (C{LatLon}).
238 @raise TypeError: The B{C{other}} point is not C{LatLon}.
240 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.
242 @raise ValueError: If this and the B{C{other}} point's L{Datum}
243 ellipsoids are not compatible.
245 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and
246 C{rhumbMidpointTo}.
247 '''
248 f = Scalar(fraction=fraction)
249 if isnear0(f):
250 r = self
251 elif isnear1(f) and not wrap:
252 r = self.others(other)
253 else: # negative fraction OK
254 t = self.distanceTo3(other, wrap=wrap)
255 h = self._havg(other, f=f) if height is None else Height(height)
256 r = self.destination(t.distance * f, t.initial, height=h)
257 return r
259 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty}
260 '''(INTERNAL) I{Karney}'s C{Inverse} method.
262 @return: A L{Distance3Tuple}C{(distance, initial, final)}.
264 @raise TypeError: The B{C{other}} point is not L{LatLon}.
266 @raise ValueError: If this and the B{C{other}} point's
267 L{Datum} ellipsoids are not compatible.
268 '''
269 _ = self.ellipsoids(other)
270 g = self.geodesic
271 _, lon = unroll180(self.lon, other.lon, wrap=wrap)
272 return g.Inverse3(self.lat, self.lon, other.lat, lon)
274 def nearestOn8(self, points, closed=False, height=None, wrap=False,
275 equidistant=None, tol=_TOL_M):
276 '''Iteratively locate the point on a path or polygon closest
277 to this point.
279 @arg points: The path or polygon points (C{LatLon}[]).
280 @kwarg closed: Optionally, close the polygon (C{bool}).
281 @kwarg height: Optional height, overriding the height of this and
282 all other points (C{meter}, conventionally). If
283 B{C{height}} is C{None}, the height of each point
284 is taken into account for distances.
286 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start,
287 end, initial, final)} with C{distance} in C{meter},
288 conventionally and with the C{closest}, the C{start}
289 the C{end} point each an instance of this C{LatLon}.
291 @raise PointsError: Insufficient number of B{C{points}}.
293 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid.
295 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible
296 or no convergence for the given B{C{tol}}.
298 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}.
299 '''
300 D3 = self.distanceTo3 # Distance3Tuple
302 try:
303 Ps = self.PointsIter(points, loop=1)
304 p1 = c = s = e = Ps[0]
305 _ = self.ellipsoids(p1)
306 c3 = D3(c, wrap=wrap) # XXX wrap=False?
308 except (TypeError, ValueError) as x:
309 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol,
310 closed=closed, height=height, wrap=wrap)
312 # get the azimuthal equidistant projection, once
313 A = _Equidistant00(equidistant, c)
314 b = _Box(c, c3.distance)
315 m = f = i = 0 # p1..p2 == points[i]..[j]
317 kwds = dict(within=True, height=height, tol=tol,
318 LatLon=self.classof, # this LatLon
319 datum=self.datum, epoch=self.epoch, reframe=self.reframe)
320 try:
321 for j, p2 in Ps.enumerate(closed=closed):
322 if wrap and j != 0:
323 p2 = _unrollon(p1, p2)
324 # skip edge if no overlap with box around closest
325 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon):
326 p, t, _ = _nearestOn3_(self, p1, p2, A, **kwds)
327 d3 = D3(p, wrap=False) # already unrolled
328 if d3.distance < c3.distance:
329 c3, c, s, e, f = d3, p, p1, p2, (i + t)
330 b = _Box(c, c3.distance)
331 m = max(m, c.iteration)
332 p1, i = p2, j
334 except (TypeError, ValueError) as x:
335 raise _xError(x, Fmt.SQUARE(points=i), p1,
336 Fmt.SQUARE(points=j), p2, this=self, tol=tol,
337 closed=closed, height=height, wrap=wrap)
339 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6
341 n = self.nearestOn8.__name__
342 c.rename(n)
343 if s is not c:
344 s = s.copy(name=n)
345 if e is not c:
346 e = e.copy(name=n)
347 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final,
348 iteration=m) # ._iteration for tests
351class _Box(object):
352 '''Bounding box around a C{LatLon} point.
354 @see: Function C{_box4} in .clipy.py.
355 '''
356 def __init__(self, center, distance):
357 '''New L{_Box} around point.
359 @arg center: The center point (C{LatLon}).
360 @arg distance: Radius, half-size of the box
361 (C{meter}, conventionally)
362 '''
363 E = center.ellipsoid()
364 d = degrees(distance / max(E.a, E.b)) + _0_5 # some margin
365 self._N = center.lat + d
366 self._S = center.lat - d
367 self._E = center.lon + d
368 self._W = center.lon - d
370 def overlaps(self, lat1, lon1, lat2, lon2):
371 '''Check whether this box overlaps a line between 2 points.
373 @arg lat1: Latitude of first point (C{degrees}).
374 @arg lon1: Longitude of first point (C{degrees}).
375 @arg lat2: Latitude of second point (C{degrees}).
376 @arg lon2: Longitude of second point (C{degrees}).
378 @return: C{False} if there is certainly no overlap,
379 C{True} otherwise (C{bool}).
380 '''
381 non_ = ((lat1 > self._N or lat2 < self._S) if lat1 < lat2 else
382 (lat2 > self._N or lat1 < self._S)) or \
383 ((lon1 > self._E or lon2 < self._W) if lon1 < lon2 else
384 (lon2 > self._E or lon1 < self._W))
385 return not non_
388class _Tol(object):
389 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}.
390 '''
391 _deg = 0 # tol in degrees
392 _lat = 0
393 _m = 0 # tol in meter
394 _min = MAX # degrees
395 _prev = None
396 _r = 0
398 def __init__(self, tol_m, E, lat, *lats):
399 '''New L{_Tol}.
401 @arg tol_m: Tolerance (C{meter}, only).
402 @arg E: Earth ellipsoid (L{Ellipsoid}).
403 @arg lat: Latitude (C{degrees}).
404 @arg lats: Additional latitudes (C{degrees}).
405 '''
406 if lats:
407 lat = fmean_(lat, *lats)
408 self._lat = lat
409 self._r = max(EPS, E.rocMean(lat))
410 self._m = max(EPS, tol_m)
411 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!
413 @property_RO
414 def degrees(self):
415 '''Get this tolerance in C{degrees}.
416 '''
417 return self._deg
419 def degrees2m(self, deg):
420 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.
421 '''
422 return self.radius * radians(deg) / PI2 # avoid degrees2m!
424 def degError(self, Error=_ValueError):
425 '''Compose an error with the C{deg}rees minimum.
426 '''
427 return self.mError(self.degrees2m(self._min), Error=Error)
429 def done(self, deg):
430 '''Check C{deg} vs. tolerance and previous value.
431 '''
432 if deg < self._deg or deg == self._prev:
433 return True
434 self._min = min(self._min, deg)
435 self._prev = deg
436 return False
438 @property_RO
439 def lat(self):
440 '''Get the mean latitude in C{degrees}.
441 '''
442 return self._lat
444 def mError(self, m, Error=_ValueError):
445 '''Compose an error with B{C{m}}eter minimum.
446 '''
447 t = _SPACE_(Fmt.tolerance(self.meter), _too_low_)
448 if m2km(m) > self.meter:
449 t = _or(t, _antipodal_, _near_(_polar__))
450 return Error(Fmt.no_convergence(m), txt=t)
452 @property_RO
453 def meter(self):
454 '''Get this tolerance in C{meter}.
455 '''
456 return self._m
458 @property_RO
459 def radius(self):
460 '''Get the earth radius in C{meter}.
461 '''
462 return self._r
464 def reset(self):
465 '''Reset tolerances.
466 '''
467 self._min = MAX
468 self._prev = None
471def _Equidistant00(equidistant, p1):
472 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
473 '''
474 if equidistant is None or not callable(equidistant):
475 equidistant = p1.Equidistant
476 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover
477 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants)
478 raise _IsnotError(*t, equidistant=equidistant)
479 return equidistant(0, 0, p1.datum)
482def _intersect3(s1, end1, s2, end2, height=None, wrap=True, # MCCABE 16
483 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
484 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method
485 L{intersection3}, separated to allow callers to embellish any exceptions.
486 '''
487 _LLS = _MODS.sphericalTrigonometry.LatLon
488 _si = _MODS.sphericalTrigonometry._intersect
489 _vi3 = _MODS.vector3d._intersect3d3
491 def _b_d(s, e, w, t, h=_0_0):
492 # compute opposing and distance
493 t = s.classof(t.lat, t.lon, height=h, name=t.name)
494 t = s.distanceTo3(t, wrap=w) # Distance3Tuple
495 b = opposing(e, t.initial) # "before" start
496 return b, t.distance
498 def _b_e(s, e, w, t):
499 # compute an end point along the initial bearing
500 # about 1.5 times the distance to the gu-/estimate, at
501 # least 1/8 and at most 3/8 of the earth perimeter like
502 # radians in .sphericalTrigonometry._int3d2 and bearing
503 # comparison in .sphericalTrigonometry._intb
504 b, d = _b_d(s, e, w, t, h=t.height)
505 m = s.ellipsoid().R2x * PI_4 # authalic exact
506 d = min(max(d * _B2END, m), m * _3_0)
507 e = s.destination(d, e)
508 return b, (_unrollon(s, e) if w else e)
510 def _e_ll(s, e, w, **end):
511 # return 2-tuple (end, False if bearing else True)
512 ll = not isscalar(e)
513 if ll:
514 e = s.others(**end)
515 if w: # unroll180 == .karney._unroll2
516 e = _unrollon(s, e)
517 return e, ll
519 def _llS(ll): # return an C{_LLS} instance
520 return _LLS(ll.lat, ll.lon, height=ll.height)
522 def _o(o, b, n, s, t, e):
523 # determine C{o}utside before, on or after start point
524 if not o: # intersection may be on start
525 if _MODS.latlonBase._isequalTo(s, t, eps=e.degrees):
526 return o
527 return -n if b else n
529 E = s1.ellipsoids(s2)
531 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
532 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
534 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),
535 s2.lat, (e2.lat if ll2 else s2.lat))
537 # get the azimuthal equidistant projection
538 A = _Equidistant00(equidistant, s1)
540 # gu-/estimate initial intersection, spherically ...
541 t = _si(_llS(s1), (_llS(e1) if ll1 else e1),
542 _llS(s2), (_llS(e2) if ll2 else e2),
543 height=height, wrap=False, LatLon=_LLS) # unrolled already
544 h, n = t.height, t.name
546 if not ll1:
547 b1, e1 = _b_e(s1, e1, wrap, t)
548 if not ll2:
549 b2, e2 = _b_e(s2, e2, wrap, t)
551 # ... and iterate as Karney describes, for references
552 # @see: Function L{ellipsoidalKarney.intersection3}.
553 for i in range(1, _TRIPS):
554 A.reset(t.lat, t.lon) # gu-/estimate as origin
555 # convert start and end points to projection
556 # space and compute an intersection there
557 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),
558 eps=e.meter, useZ=False)
559 # convert intersection back to geodetic
560 t, d = A._reverse2(v)
561 if e.done(d): # below tol or unchanged?
562 t._iteration = i
563 break
564 else:
565 raise e.degError(Error=IntersectionError)
567 # like .sphericalTrigonometry._intersect, if this intersection
568 # is "before" the first point, use the antipodal intersection
569 if not (ll1 or ll2): # end1 and end2 are an initial bearing
570 b1, _ = _b_d(s1, end1, wrap, t)
571 if b1:
572 t = t.antipodal()
573 b1 = not b1
574 b2, _ = _b_d(s2, end2, wrap, t)
576 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1,
577 iteration=t._iteration, name=n)
578 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),
579 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))
582def _intersection3(start1, end1, start2, end2, height=None, wrap=True,
583 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
584 '''(INTERNAL) Iteratively compute the intersection point of two lines,
585 each defined by two (ellipsoidal) points or an (ellipsoidal) start
586 point and an initial bearing from North.
587 '''
588 s1 = _xellipsoidal(start1=start1)
589 s2 = s1.others(start2=start2)
590 try:
591 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,
592 equidistant=equidistant, tol=tol,
593 LatLon=LatLon, **LatLon_kwds)
594 except (TypeError, ValueError) as x:
595 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
598def _intersections2(center1, radius1, center2, radius2, height=None, wrap=True,
599 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
600 '''(INTERNAL) Iteratively compute the intersection points of two circles,
601 each defined by an (ellipsoidal) center point and a radius.
602 '''
603 c1 = _xellipsoidal(center1=center1)
604 c2 = c1.others(center2=center2)
605 try:
606 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,
607 equidistant=equidistant, tol=tol,
608 LatLon=LatLon, **LatLon_kwds)
609 except (TypeError, ValueError) as x:
610 raise _xError(x, center1=center1, radius1=radius1,
611 center2=center2, radius2=radius2)
614def _intersects2(c1, radius1, c2, radius2, height=None, wrap=True, # MCCABE 16
615 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
616 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}
617 above, separated to allow callers to embellish any exceptions.
618 '''
619 _LLS = _MODS.sphericalTrigonometry.LatLon
620 _si2 = _MODS.sphericalTrigonometry._intersects2
621 _vi2 = _MODS.vector3d._intersects2
623 def _latlon4(t, h, n, c):
624 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c,
625 iteration=t.iteration, name=n)
627 def _llS(ll): # return an C{_LLS} instance
628 return _LLS(ll.lat, ll.lon, height=ll.height)
630 r1 = Radius_(radius1=radius1)
631 r2 = Radius_(radius2=radius2)
633 E = c1.ellipsoids(c2)
634 # get the azimuthal equidistant projection
635 A = _Equidistant00(equidistant, c1)
637 if r1 < r2:
638 c1, c2 = c2, c1
639 r1, r2 = r2, r1
641 if r1 > (min(E.b, E.a) * PI):
642 raise _ValueError(_exceed_PI_radians_)
644 if wrap: # unroll180 == .karney._unroll2
645 c2 = _unrollon(c1, c2)
647 # distance between centers and radii are
648 # measured along the ellipsoid's surface
649 m = c1.distanceTo(c2, wrap=False) # meter
650 if m < max(r1 - r2, EPS):
651 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
652 if fsum_(r1, r2, -m) < 0:
653 raise IntersectionError(_too_(Fmt.distant(m)))
655 f = _radical2(m, r1, r2).ratio # "radical fraction"
656 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))
658 # gu-/estimate initial intersections, spherically ...
659 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius,
660 height=height, wrap=False, too_d=m) # unrolled already
661 h, n = t1.height, t1.name
663 # ... and iterate as Karney describes, for references
664 # @see: Function L{ellipsoidalKarney.intersections2}.
665 ts, ta = [], None
666 for t in ((t1,) if t1 is t2 else (t1, t2)):
667 for i in range(1, _TRIPS):
668 A.reset(t.lat, t.lon) # gu-/estimate as origin
669 # convert centers to projection space and
670 # compute the intersections there
671 t1, t2 = A._forwards(c1, c2)
672 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,
673 t2, r2, # XXX * t2.scale?,
674 sphere=False, too_d=m)
675 # convert intersections back to geodetic
676 if v1 is v2: # abutting
677 t, d = A._reverse2(v1)
678 else: # consider the closer intersection
679 t1, d1 = A._reverse2(v1)
680 t2, d2 = A._reverse2(v2)
681 t, d = (t1, d1) if d1 < d2 else (t2, d2)
682 if e.done(d): # below tol or unchanged?
683 t._iteration = i # _NamedTuple._iteration
684 ts.append(t)
685 if v1 is v2: # abutting
686 ta = t
687 break
688 else:
689 raise e.degError(Error=IntersectionError)
690 e.reset()
692 if ta: # abutting circles
693 pass # PYCHOK no cover
694 elif len(ts) == 2:
695 return (_latlon4(ts[0], h, n, c1),
696 _latlon4(ts[1], h, n, c2))
697 elif len(ts) == 1: # PYCHOK no cover
698 ta = ts[0] # assume abutting
699 else: # PYCHOK no cover
700 raise _AssertionError(ts=ts)
701 r = _latlon4(ta, h, n, c1)
702 return r, r
705def _nearestOn2(p, point1, point2, within=True, height=None, wrap=True,
706 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):
707 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,
708 separated to allow callers to embellish any exceptions.
709 '''
710 p1 = p.others(point1=point1)
711 p2 = p.others(point2=point2)
713 _ = p.ellipsoids(p1)
714# E = p.ellipsoids(p2) # done in __nearestOn2__
716 # get the azimuthal equidistant projection
717 A = _Equidistant00(equidistant, p)
719 if wrap:
720 p1 = _unrollon(p, p1) # XXX do not unroll?
721 p2 = _unrollon(p, p2) # XXX do not unroll?
722 p2 = _unrollon(p1, p2)
724 r, f, e = _nearestOn3_(p, p1, p2, A, within=within, height=height,
725 tol=tol, **LatLon_and_kwds)
726 return NearestOn2Tuple(r, f)
729def _nearestOn3_(p, p1, p2, A, within=True, height=None, tol=_TOL_M,
730 LatLon=None, **LatLon_kwds):
731 # Only in function C{_nearestOn2} and method C{nearestOn8} above
732 _LLS = _MODS.sphericalNvector.LatLon
733 _V3d = _MODS.vector3d.Vector3d
734 _vnOn2 = _MODS.vector3d._nearestOn2
736 def _llS(ll): # return an C{_LLS} instance
737 return _LLS(ll.lat, ll.lon, height=ll.height)
739 E = p.ellipsoids(p2)
740 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)
742 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
743 # using sphericalNvector.LatLon.nearestOn for within=False support
744 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within,
745 height=height)
746 n, h = t.name, t.height
747 if height is None:
748 h1 = p1.height # use heights as pseudo-Z in projection space
749 h2 = p2.height # to be included in the closest function
750 h0 = favg(h1, h2)
751 else: # ignore heights in distances, Z=0
752 h0 = h1 = h2 = _0_0
754 # ... and iterate to find the closest (to the origin with .z
755 # to interpolate height) as Karney describes, for references
756 # @see: Function L{ellipsoidalKarney.nearestOn}.
757 vp, f = _V3d(_0_0, _0_0, h0), None
758 for i in range(1, _TRIPS):
759 A.reset(t.lat, t.lon) # gu-/estimate as origin
760 # convert points to projection space and compute
761 # the nearest one (and its height) there
762 s, t = A._forwards(p1, p2)
763 v, f = _vnOn2(vp, _V3d(s.x, s.y, h1),
764 _V3d(t.x, t.y, h2), within=within)
765 # convert nearest one back to geodetic
766 t, d = A._reverse2(v)
767 if e.done(d): # below tol or unchanged?
768 t._iteration = i # _NamedTuple._iteration
769 break
770 else:
771 raise e.degError()
773 if height is None:
774 h = v.z # nearest
775 elif isscalar(height):
776 h = height
777 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,
778 iteration=t.iteration, name=n)
779 return r, f, e # fraction or None
782__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)
784# **) MIT License
785#
786# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
787#
788# Permission is hereby granted, free of charge, to any person obtaining a
789# copy of this software and associated documentation files (the "Software"),
790# to deal in the Software without restriction, including without limitation
791# the rights to use, copy, modify, merge, publish, distribute, sublicense,
792# and/or sell copies of the Software, and to permit persons to whom the
793# Software is furnished to do so, subject to the following conditions:
794#
795# The above copyright notice and this permission notice shall be included
796# in all copies or substantial portions of the Software.
797#
798# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
799# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
800# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
801# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
802# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
803# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
804# OTHER DEALINGS IN THE SOFTWARE.