Coverage for pygeodesy/ellipsoidalBaseDI.py: 96%
310 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -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.23'
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 self._lat = fmean_(lat, *lats) if lats else lat
407 self._r = max(EPS, E.rocMean(self._lat))
408 self._m = max(EPS, tol_m)
409 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!
411 @property_RO
412 def degrees(self):
413 '''Get this tolerance in C{degrees}.
414 '''
415 return self._deg
417 def degrees2m(self, deg):
418 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.
419 '''
420 return self.radius * radians(deg) / PI2 # avoid degrees2m!
422 def degError(self, Error=_ValueError):
423 '''Compose an error with the C{deg}rees minimum.
424 '''
425 return self.mError(self.degrees2m(self._min), Error=Error)
427 def done(self, deg):
428 '''Check C{deg} vs. tolerance and previous value.
429 '''
430 if deg < self._deg or deg == self._prev:
431 return True
432 self._min = min(self._min, deg)
433 self._prev = deg
434 return False
436 @property_RO
437 def lat(self):
438 '''Get the mean latitude in C{degrees}.
439 '''
440 return self._lat
442 def mError(self, m, Error=_ValueError):
443 '''Compose an error with B{C{m}}eter minimum.
444 '''
445 t = _SPACE_(Fmt.tolerance(self.meter), _too_low_)
446 if m2km(m) > self.meter:
447 t = _or(t, _antipodal_, _near_(_polar__))
448 return Error(Fmt.no_convergence(m), txt=t)
450 @property_RO
451 def meter(self):
452 '''Get this tolerance in C{meter}.
453 '''
454 return self._m
456 @property_RO
457 def radius(self):
458 '''Get the earth radius in C{meter}.
459 '''
460 return self._r
462 def reset(self):
463 '''Reset tolerances.
464 '''
465 self._min = MAX
466 self._prev = None
469def _Equidistant00(equidistant, p1):
470 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.
471 '''
472 if equidistant is None or not callable(equidistant):
473 equidistant = p1.Equidistant
474 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover
475 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants)
476 raise _IsnotError(*t, equidistant=equidistant)
477 return equidistant(0, 0, p1.datum)
480def _intersect3(s1, end1, s2, end2, height=None, wrap=True, # MCCABE 16
481 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
482 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method
483 L{intersection3}, separated to allow callers to embellish any exceptions.
484 '''
485 _LLS = _MODS.sphericalTrigonometry.LatLon
486 _si = _MODS.sphericalTrigonometry._intersect
487 _vi3 = _MODS.vector3d._intersect3d3
489 def _b_d(s, e, w, t, h=_0_0):
490 # compute opposing and distance
491 t = s.classof(t.lat, t.lon, height=h, name=t.name)
492 t = s.distanceTo3(t, wrap=w) # Distance3Tuple
493 b = opposing(e, t.initial) # "before" start
494 return b, t.distance
496 def _b_e(s, e, w, t):
497 # compute an end point along the initial bearing
498 # about 1.5 times the distance to the gu-/estimate, at
499 # least 1/8 and at most 3/8 of the earth perimeter like
500 # radians in .sphericalTrigonometry._int3d2 and bearing
501 # comparison in .sphericalTrigonometry._intb
502 b, d = _b_d(s, e, w, t, h=t.height)
503 m = s.ellipsoid().R2x * PI_4 # authalic exact
504 d = min(max(d * _B2END, m), m * _3_0)
505 e = s.destination(d, e)
506 return b, (_unrollon(s, e) if w else e)
508 def _e_ll(s, e, w, **end):
509 # return 2-tuple (end, False if bearing else True)
510 ll = not isscalar(e)
511 if ll:
512 e = s.others(**end)
513 if w: # unroll180 == .karney._unroll2
514 e = _unrollon(s, e)
515 return e, ll
517 def _llS(ll): # return an C{_LLS} instance
518 return _LLS(ll.lat, ll.lon, height=ll.height)
520 def _o(o, b, n, s, t, e):
521 # determine C{o}utside before, on or after start point
522 if not o: # intersection may be on start
523 if _MODS.latlonBase._isequalTo(s, t, eps=e.degrees):
524 return o
525 return -n if b else n
527 E = s1.ellipsoids(s2)
529 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)
530 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)
532 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),
533 s2.lat, (e2.lat if ll2 else s2.lat))
535 # get the azimuthal equidistant projection
536 A = _Equidistant00(equidistant, s1)
538 # gu-/estimate initial intersection, spherically ...
539 t = _si(_llS(s1), (_llS(e1) if ll1 else e1),
540 _llS(s2), (_llS(e2) if ll2 else e2),
541 height=height, wrap=False, LatLon=_LLS) # unrolled already
542 h, n = t.height, t.name
544 if not ll1:
545 b1, e1 = _b_e(s1, e1, wrap, t)
546 if not ll2:
547 b2, e2 = _b_e(s2, e2, wrap, t)
549 # ... and iterate as Karney describes, for references
550 # @see: Function L{ellipsoidalKarney.intersection3}.
551 for i in range(1, _TRIPS):
552 A.reset(t.lat, t.lon) # gu-/estimate as origin
553 # convert start and end points to projection
554 # space and compute an intersection there
555 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),
556 eps=e.meter, useZ=False)
557 # convert intersection back to geodetic
558 t, d = A._reverse2(v)
559 if e.done(d): # below tol or unchanged?
560 t._iteration = i
561 break
562 else:
563 raise e.degError(Error=IntersectionError)
565 # like .sphericalTrigonometry._intersect, if this intersection
566 # is "before" the first point, use the antipodal intersection
567 if not (ll1 or ll2): # end1 and end2 are an initial bearing
568 b1, _ = _b_d(s1, end1, wrap, t)
569 if b1:
570 t = t.antipodal()
571 b1 = not b1
572 b2, _ = _b_d(s2, end2, wrap, t)
574 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1,
575 iteration=t._iteration, name=n)
576 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),
577 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))
580def _intersection3(start1, end1, start2, end2, height=None, wrap=True,
581 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
582 '''(INTERNAL) Iteratively compute the intersection point of two lines,
583 each defined by two (ellipsoidal) points or an (ellipsoidal) start
584 point and an initial bearing from North.
585 '''
586 s1 = _xellipsoidal(start1=start1)
587 s2 = s1.others(start2=start2)
588 try:
589 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,
590 equidistant=equidistant, tol=tol,
591 LatLon=LatLon, **LatLon_kwds)
592 except (TypeError, ValueError) as x:
593 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
596def _intersections2(center1, radius1, center2, radius2, height=None, wrap=True,
597 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
598 '''(INTERNAL) Iteratively compute the intersection points of two circles,
599 each defined by an (ellipsoidal) center point and a radius.
600 '''
601 c1 = _xellipsoidal(center1=center1)
602 c2 = c1.others(center2=center2)
603 try:
604 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,
605 equidistant=equidistant, tol=tol,
606 LatLon=LatLon, **LatLon_kwds)
607 except (TypeError, ValueError) as x:
608 raise _xError(x, center1=center1, radius1=radius1,
609 center2=center2, radius2=radius2)
612def _intersects2(c1, radius1, c2, radius2, height=None, wrap=True, # MCCABE 16
613 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):
614 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}
615 above, separated to allow callers to embellish any exceptions.
616 '''
617 _LLS = _MODS.sphericalTrigonometry.LatLon
618 _si2 = _MODS.sphericalTrigonometry._intersects2
619 _vi2 = _MODS.vector3d._intersects2
621 def _latlon4(t, h, n, c):
622 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c,
623 iteration=t.iteration, name=n)
625 def _llS(ll): # return an C{_LLS} instance
626 return _LLS(ll.lat, ll.lon, height=ll.height)
628 r1 = Radius_(radius1=radius1)
629 r2 = Radius_(radius2=radius2)
631 E = c1.ellipsoids(c2)
632 # get the azimuthal equidistant projection
633 A = _Equidistant00(equidistant, c1)
635 if r1 < r2:
636 c1, c2 = c2, c1
637 r1, r2 = r2, r1
639 if r1 > (min(E.b, E.a) * PI):
640 raise _ValueError(_exceed_PI_radians_)
642 if wrap: # unroll180 == .karney._unroll2
643 c2 = _unrollon(c1, c2)
645 # distance between centers and radii are
646 # measured along the ellipsoid's surface
647 m = c1.distanceTo(c2, wrap=False) # meter
648 if m < max(r1 - r2, EPS):
649 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
650 if fsum_(r1, r2, -m) < 0:
651 raise IntersectionError(_too_(Fmt.distant(m)))
653 f = _radical2(m, r1, r2).ratio # "radical fraction"
654 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))
656 # gu-/estimate initial intersections, spherically ...
657 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius,
658 height=height, wrap=False, too_d=m) # unrolled already
659 h, n = t1.height, t1.name
661 # ... and iterate as Karney describes, for references
662 # @see: Function L{ellipsoidalKarney.intersections2}.
663 ts, ta = [], None
664 for t in ((t1,) if t1 is t2 else (t1, t2)):
665 for i in range(1, _TRIPS):
666 A.reset(t.lat, t.lon) # gu-/estimate as origin
667 # convert centers to projection space and
668 # compute the intersections there
669 t1, t2 = A._forwards(c1, c2)
670 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,
671 t2, r2, # XXX * t2.scale?,
672 sphere=False, too_d=m)
673 # convert intersections back to geodetic
674 if v1 is v2: # abutting
675 t, d = A._reverse2(v1)
676 else: # consider the closer intersection
677 t1, d1 = A._reverse2(v1)
678 t2, d2 = A._reverse2(v2)
679 t, d = (t1, d1) if d1 < d2 else (t2, d2)
680 if e.done(d): # below tol or unchanged?
681 t._iteration = i # _NamedTuple._iteration
682 ts.append(t)
683 if v1 is v2: # abutting
684 ta = t
685 break
686 else:
687 raise e.degError(Error=IntersectionError)
688 e.reset()
690 if ta: # abutting circles
691 pass # PYCHOK no cover
692 elif len(ts) == 2:
693 return (_latlon4(ts[0], h, n, c1),
694 _latlon4(ts[1], h, n, c2))
695 elif len(ts) == 1: # PYCHOK no cover
696 ta = ts[0] # assume abutting
697 else: # PYCHOK no cover
698 raise _AssertionError(ts=ts)
699 r = _latlon4(ta, h, n, c1)
700 return r, r
703def _nearestOn2(p, point1, point2, within=True, height=None, wrap=True,
704 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):
705 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,
706 separated to allow callers to embellish any exceptions.
707 '''
708 p1 = p.others(point1=point1)
709 p2 = p.others(point2=point2)
711 _ = p.ellipsoids(p1)
712# E = p.ellipsoids(p2) # done in _nearestOn3
714 # get the azimuthal equidistant projection
715 A = _Equidistant00(equidistant, p)
717 if wrap:
718 p1 = _unrollon(p, p1) # XXX do not unroll?
719 p2 = _unrollon(p, p2) # XXX do not unroll?
720 p2 = _unrollon(p1, p2)
722 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height,
723 tol=tol, **LatLon_and_kwds)
724 return NearestOn2Tuple(r, f)
727def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M,
728 LatLon=None, **LatLon_kwds):
729 # Only in function C{_nearestOn2} and method C{nearestOn8} above
730 _LLS = _MODS.sphericalNvector.LatLon
731 _V3d = _MODS.vector3d.Vector3d
732 _vn2 = _MODS.vector3d._nearestOn2
734 def _llS(ll): # return an C{_LLS} instance
735 return _LLS(ll.lat, ll.lon, height=ll.height)
737 E = p.ellipsoids(p2)
738 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)
740 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!
741 # using sphericalNvector.LatLon.nearestOn for within=False support
742 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within,
743 height=height)
744 n, h = t.name, t.height
745 if height is None:
746 h1 = p1.height # use heights as pseudo-Z in projection space
747 h2 = p2.height # to be included in the closest function
748 h0 = favg(h1, h2)
749 else: # ignore heights in distances, Z=0
750 h0 = h1 = h2 = _0_0
752 # ... and iterate to find the closest (to the origin with .z
753 # to interpolate height) as Karney describes, for references
754 # @see: Function L{ellipsoidalKarney.nearestOn}.
755 vp, f = _V3d(_0_0, _0_0, h0), None
756 for i in range(1, _TRIPS):
757 A.reset(t.lat, t.lon) # gu-/estimate as origin
758 # convert points to projection space and compute
759 # the nearest one (and its height) there
760 s, t = A._forwards(p1, p2)
761 v, f = _vn2(vp, _V3d(s.x, s.y, h1),
762 _V3d(t.x, t.y, h2), within=within)
763 # convert nearest one back to geodetic
764 t, d = A._reverse2(v)
765 if e.done(d): # below tol or unchanged?
766 t._iteration = i # _NamedTuple._iteration
767 break
768 else:
769 raise e.degError()
771 if height is None:
772 h = v.z # nearest
773 elif isscalar(height):
774 h = height
775 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,
776 iteration=t.iteration, name=n)
777 return r, f, e # fraction or None
780__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)
782# **) MIT License
783#
784# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
785#
786# Permission is hereby granted, free of charge, to any person obtaining a
787# copy of this software and associated documentation files (the "Software"),
788# to deal in the Software without restriction, including without limitation
789# the rights to use, copy, modify, merge, publish, distribute, sublicense,
790# and/or sell copies of the Software, and to permit persons to whom the
791# Software is furnished to do so, subject to the following conditions:
792#
793# The above copyright notice and this permission notice shall be included
794# in all copies or substantial portions of the Software.
795#
796# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
797# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
798# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
799# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
800# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
801# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
802# OTHER DEALINGS IN THE SOFTWARE.