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