Coverage for pygeodesy/sphericalBase.py: 94%
251 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-02 13:46 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-02 13:46 -0500
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private spherical base classes C{CartesianSphericalBase} and
5C{LatLonSphericalBase} for L{sphericalNvector} and L{sphericalTrigonometry}.
7A pure Python implementation of geodetic (lat-/longitude) functions,
8transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016}
9and published under the same MIT Licence**, see
10U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.
11'''
12# make sure int/int division yields float quotient, see .basics
13from __future__ import division as _; del _ # PYCHOK semicolon
15from pygeodesy.basics import _copysign, isbool, isinstanceof, isscalar, map1
16from pygeodesy.cartesianBase import CartesianBase, Bearing2Tuple
17from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \
18 _0_0, _0_5, _1_0, _180_0, _360_0, \
19 _over, isnear0, isnon0
20from pygeodesy.datums import Datums, _earth_ellipsoid, _spherical_datum
21from pygeodesy.errors import IntersectionError, _ValueError, \
22 _xattr, _xError
23from pygeodesy.fmath import favg, fdot, hypot, sqrt_a
24from pygeodesy.interns import NN, _COMMA_, _concentric_, _datum_, \
25 _distant_, _exceed_PI_radians_, _name_, \
26 _near_, _radius_, _too_
27from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed
28from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
29# from pygeodesy.namedTuples import Bearing2Tuple # from .cartesianBase
30from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs
31from pygeodesy.props import deprecated_method, property_doc_, \
32 property_RO, _update_all
33# from pygeodesy.streprs import Fmt, _xattrs # from .nvectorBase
34from pygeodesy.units import Bearing, Bearing_, Radians_, Radius, \
35 Radius_, Scalar_, _100km
36from pygeodesy.utily import acos1, asin1, atan2b, atan2d, degrees90, \
37 degrees180, sincos2, sincos2d, _unrollon, \
38 tanPI_2_2, wrapPI
40from math import cos, fabs, log, sin, sqrt
42__all__ = _ALL_LAZY.sphericalBase
43__version__ = '23.12.01'
46class CartesianSphericalBase(CartesianBase):
47 '''(INTERNAL) Base class for spherical C{Cartesian}s.
48 '''
49 _datum = Datums.Sphere # L{Datum}
51 def intersections2(self, rad1, other, rad2, radius=R_M):
52 '''Compute the intersection points of two circles each defined
53 by a center point and a radius.
55 @arg rad1: Radius of the this circle (C{meter} or C{radians},
56 see B{C{radius}}).
57 @arg other: Center of the other circle (C{Cartesian}).
58 @arg rad2: Radius of the other circle (C{meter} or C{radians},
59 see B{C{radius}}).
60 @kwarg radius: Mean earth radius (C{meter} or C{None} if both
61 B{C{rad1}} and B{C{rad2}} are given in C{radians}).
63 @return: 2-Tuple of the intersection points, each C{Cartesian}.
64 For abutting circles, the intersection points are the
65 same C{Cartesian} instance, aka the I{radical center}.
67 @raise IntersectionError: Concentric, antipodal, invalid or
68 non-intersecting circles.
70 @raise TypeError: If B{C{other}} is not C{Cartesian}.
72 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}.
74 @see: U{Calculating intersection of two Circles
75 <https://GIS.StackExchange.com/questions/48937/
76 calculating-intersection-of-two-circles>} and method
77 or function C{trilaterate3d2}.
78 '''
79 x1, x2 = self, self.others(other)
80 r1, r2, x = _rads3(rad1, rad2, radius)
81 if x:
82 x1, x2 = x2, x1
83 try:
84 n, q = x1.cross(x2), x1.dot(x2)
85 n2, q1 = n.length2, (_1_0 - q**2)
86 if n2 < EPS or isnear0(q1):
87 raise ValueError(_near_(_concentric_))
88 c1, c2 = cos(r1), cos(r2)
89 x0 = x1.times((c1 - q * c2) / q1).plus(
90 x2.times((c2 - q * c1) / q1))
91 n1 = _1_0 - x0.length2
92 if n1 < EPS:
93 raise ValueError(_too_(_distant_))
94 except ValueError as x:
95 raise IntersectionError(center=self, rad1=rad1,
96 other=other, rad2=rad2, cause=x)
97 n = n.times(sqrt(n1 / n2))
98 if n.length > EPS:
99 x1 = x0.plus(n)
100 x2 = x0.minus(n)
101 else: # abutting circles
102 x1 = x2 = x0
104 return (_xattrs(x1, self, _datum_, _name_),
105 _xattrs(x2, self, _datum_, _name_))
107 @property_RO
108 def sphericalCartesian(self):
109 '''Get this C{Cartesian}'s spherical class.
110 '''
111 return type(self)
114class LatLonSphericalBase(LatLonBase):
115 '''(INTERNAL) Base class for spherical C{LatLon}s.
116 '''
117 _datum = Datums.Sphere # spherical L{Datum}
118 _napieradius = _100km
120 def __init__(self, latlonh, lon=None, height=0, datum=None, wrap=False, name=NN):
121 '''Create a spherical C{LatLon} point frome the given lat-, longitude and
122 height on the given datum.
124 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or
125 a previous C{LatLon} instance provided C{B{lon}=None}.
126 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or
127 C(None), indicating B{C{latlonh}} is a C{LatLon}.
128 @kwarg height: Optional height above (or below) the earth surface (C{meter},
129 same units as the datum's ellipsoid axes or radius).
130 @kwarg datum: Optional, spherical datum to use (L{Datum}, L{Ellipsoid},
131 L{Ellipsoid2}, L{a_f2Tuple}) or earth radius in C{meter},
132 conventionally).
133 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}
134 (C{bool}).
135 @kwarg name: Optional name (string).
137 @raise TypeError: If B{C{latlonh}} is not a C{LatLon} or B{C{datum}} not
138 spherical.
139 '''
140 LatLonBase.__init__(self, latlonh, lon=lon, height=height, wrap=wrap, name=name)
141 if datum not in (None, self.datum):
142 self.datum = datum
144 def bearingTo2(self, other, wrap=False, raiser=False):
145 '''Return the initial and final bearing (forward and reverse
146 azimuth) from this to an other point.
148 @arg other: The other point (C{LatLon}).
149 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
150 B{C{other}} point (C{bool}).
152 @return: A L{Bearing2Tuple}C{(initial, final)}.
154 @raise TypeError: The B{C{other}} point is not spherical.
156 @see: Methods C{initialBearingTo} and C{finalBearingTo}.
157 '''
158 # .initialBearingTo is inside .-Nvector and .-Trigonometry
159 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo
160 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser)
161 return Bearing2Tuple(i, f, name=self.name)
163 @property_doc_(''' this point's datum (L{Datum}).''')
164 def datum(self):
165 '''Get this point's datum (L{Datum}).
166 '''
167 return self._datum
169 @datum.setter # PYCHOK setter!
170 def datum(self, datum):
171 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid},
172 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius).
174 @raise TypeError: If B{C{datum}} invalid or not not spherical.
175 '''
176 d = _spherical_datum(datum, name=self.name, raiser=_datum_)
177 if self._datum != d:
178 _update_all(self)
179 self._datum = d
181 def finalBearingTo(self, other, wrap=False, raiser=False):
182 '''Return the final bearing (reverse azimuth) from this to
183 an other point.
185 @arg other: The other point (spherical C{LatLon}).
186 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
187 the B{C{other}} point (C{bool}).
189 @return: Final bearing (compass C{degrees360}).
191 @raise TypeError: The B{C{other}} point is not spherical.
193 @example:
195 >>> p = LatLon(52.205, 0.119)
196 >>> q = LatLon(48.857, 2.351)
197 >>> b = p.finalBearingTo(q) # 157.9
198 '''
199 p = self.others(other)
200 if wrap:
201 p = _unrollon(self, p, wrap=wrap)
202 # final bearing is the reverse of the other, initial one
203 b = p.initialBearingTo(self, wrap=False, raiser=raiser) + _180_0
204 return b if b < 360 else (b - _360_0)
206 def intersecant2(self, circle, point, other, radius=R_M, exact=False, # PYCHOK signature
207 height=None, wrap=False):
208 '''Compute the intersections of a circle and a (great circle) line
209 given as two points or as a point and bearing.
211 @arg circle: Radius of the circle centered at this location (C{meter},
212 same units as B{C{radius}}) or a point on the circle
213 (this C{LatLon}).
214 @arg point: A point on the (great circle) line (this C{LatLon}).
215 @arg other: An other point I{on} (this {LatLon}) or the bearing at
216 B{C{point}} I{of} the (great circle) line (compass
217 C{degrees}).
218 @kwarg radius: Mean earth radius (C{meter}, conventionally).
219 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
220 destination and distance, if C{False} use the basic
221 rhumb methods (C{bool}) or if C{None} use the I{great
222 circle} methods.
223 @kwarg height: Optional height for the intersection points (C{meter},
224 conventionally) or C{None} for interpolated heights.
225 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the points
226 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}).
228 @return: 2-Tuple of the intersection points (representing a chord), each
229 an instance of the B{C{point}} class. Both points are the same
230 instance if the (great circle) line is tangent to the circle.
232 @raise IntersectionError: The circle and line do not intersect.
234 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
235 or B{C{other}} invalid.
237 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
238 B{C{exact}}, B{C{height}} or B{C{napieradius}}.
239 '''
240 p = self.others(point=point)
241 try:
242 return _intersecant2(self, circle, p, other, radius=radius, exact=exact,
243 height=height, wrap=wrap)
244 except (TypeError, ValueError) as x:
245 raise _xError(x, center=self, circle=circle, point=point, other=other,
246 radius=radius, exact=exact, height=height, wrap=wrap)
248 def maxLat(self, bearing):
249 '''Return the maximum latitude reached when travelling on a great circle
250 on given bearing from this point based on Clairaut's formula.
252 The maximum latitude is independent of longitude and the same for all
253 points on a given latitude.
255 Negate the result for the minimum latitude (on the Southern hemisphere).
257 @arg bearing: Initial bearing (compass C{degrees360}).
259 @return: Maximum latitude (C{degrees90}).
261 @raise ValueError: Invalid B{C{bearing}}.
262 '''
263 r = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi)))
264 return degrees90(r)
266 def minLat(self, bearing):
267 '''Return the minimum latitude reached when travelling on a great circle
268 on given bearing from this point.
270 @arg bearing: Initial bearing (compass C{degrees360}).
272 @return: Minimum latitude (C{degrees90}).
274 @see: Method L{maxLat} for more details.
276 @raise ValueError: Invalid B{C{bearing}}.
277 '''
278 return -self.maxLat(bearing)
280 def _mpr(self, radius=R_M, exact=None): # meter per radian
281 if exact and not isscalar(radius): # see .rhumb.ekx.Rhumb._mpr
282 radius = _earth_ellipsoid(radius)._Lpr
283 return radius
285 @property_doc_(''' the I{Napier} radius to apply spherical trigonometry.''')
286 def napieradius(self):
287 '''Get the I{Napier} radius (C{meter}, conventionally).
288 '''
289 return self._napieradius
291 @napieradius.setter # PYCHOK setter!
292 def napieradius(self, radius):
293 '''Set this I{Napier} radius (C{meter}, conventionally) or C{0}.
295 In methods L{intersecant2} and L{rhumbIntersecant2}, I{Napier}'s
296 spherical trigonometry is applied if the circle radius exceeds
297 the I{Napier} radius, otherwise planar trigonometry is used.
299 @raise UnitError: Invalid B{C{radius}}.
300 '''
301 self._napieradius = Radius(napieradius=radius or 0)
303# def nearestTo(self, point, other, **radius_exact_height_wrap): # PYCHOK signature
304# p = self.others(point=point)
305# try:
306# p, q = _intersecant2(self, p, p, other, **radius_exact_height_wrap)
307# except (TypeError, ValueError) as x:
308# raise _xError(x, this=self, point=point, other=other, **radius_exact_height_wrap)
309# return p.midpointTo(q)
311 def parse(self, strllh, height=0, sep=_COMMA_, name=NN):
312 '''Parse a string representing a similar, spherical C{LatLon}
313 point, consisting of C{"lat, lon[, height]"}.
315 @arg strllh: Lat, lon and optional height (C{str}),
316 see function L{pygeodesy.parse3llh}.
317 @kwarg height: Optional, default height (C{meter}).
318 @kwarg sep: Optional separator (C{str}).
319 @kwarg name: Optional instance name (C{str}),
320 overriding this name.
322 @return: The similar point (spherical C{LatLon}).
324 @raise ParseError: Invalid B{C{strllh}}.
325 '''
326 t = _MODS.dms.parse3llh(strllh, height=height, sep=sep)
327 r = self.classof(*t)
328 if name:
329 r.rename(name)
330 return r
332 @property_RO
333 def _radius(self):
334 '''(INTERNAL) Get this sphere's radius.
335 '''
336 return self.datum.ellipsoid.equatoradius
338 def _rhumbs3(self, other, wrap, r=False): # != .latlonBase._rhumbx3
339 '''(INTERNAL) Rhumb_ helper function.
341 @arg other: The other point (spherical C{LatLon}).
342 '''
343 p = self.others(other, up=2)
344 if wrap:
345 p = _unrollon(self, p, wrap=wrap)
346 a2, b2 = p.philam
347 a1, b1 = self.philam
348 # if |db| > 180 take shorter rhumb
349 # line across the anti-meridian
350 db = wrapPI(b2 - b1)
351 dp = _logPI_2_2(a2, a1)
352 da = a2 - a1
353 if r:
354 # on Mercator projection, longitude distances shrink
355 # by latitude; the 'stretch factor' q becomes ill-
356 # conditioned along E-W line (0/0); use an empirical
357 # tolerance to avoid it
358 q = (da / dp) if fabs(dp) > EPS else cos(a1)
359 da = hypot(da, q * db) # angular distance radians
360 return da, db, dp
362 def rhumbAzimuthTo(self, other, radius=R_M, exact=False, wrap=False, b360=False):
363 '''Return the azimuth (bearing) of a rhumb line (loxodrome) between
364 this and an other (spherical) point.
366 @arg other: The other point (spherical C{LatLon}).
367 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
368 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
369 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),
370 default C{False} for backward compatibility.
371 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
372 B{C{other}} point (C{bool}).
373 @kwarg b360: If C{True}, return the azimuth in the bearing range.
375 @return: Rhumb azimuth (compass C{degrees180} or C{degrees360}).
377 @raise TypeError: The B{C{other}} point is incompatible or
378 B{C{radius}} is invalid.
380 @example:
382 >>> p = LatLon(51.127, 1.338)
383 >>> q = LatLon(50.964, 1.853)
384 >>> b = p.rhumbBearingTo(q) # 116.7
385 '''
386 if exact: # use series, always
387 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False, # Krüger
388 radius=radius, wrap=wrap, b360=b360)
389 else:
390 _, db, dp = self._rhumbs3(other, wrap)
391 z = (atan2b if b360 else atan2d)(db, dp) # see .rhumbBase.RhumbBase.Inverse
392 return z
394 @deprecated_method
395 def rhumbBearingTo(self, other): # unwrapped
396 '''DEPRECATED, use method C{.rhumbAzimuthTo}.'''
397 return self.rhumbAzimuthTo(other, b360=True) # [0..360)
399 def rhumbDestination(self, distance, azimuth, radius=R_M, height=None, exact=False):
400 '''Return the destination point having travelled the given distance from
401 this point along a rhumb line (loxodrome) of the given azimuth.
403 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}),
404 may be negative if C{B{exact}=True}.
405 @arg azimuth: Azimuth (bearing) of the rhumb line (compass C{degrees}).
406 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
407 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
408 C{B{exact}=True}.
409 @kwarg height: Optional height, overriding the default height (C{meter}.
410 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),
411 default C{False} for backward compatibility.
413 @return: The destination point (spherical C{LatLon}).
415 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}}, B{C{radius}}
416 or B{C{height}}.
418 @example:
420 >>> p = LatLon(51.127, 1.338)
421 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E
422 '''
423 if exact: # use series, always
424 r = LatLonBase.rhumbDestination(self, distance, azimuth, exact=False, # Krüger
425 radius=radius, height=height)
426 else: # radius=None from .rhumbMidpointTo
427 if radius in (None, self._radius):
428 d, r = self.datum, radius
429 else:
430 d = _spherical_datum(radius, raiser=_radius_) # spherical only
431 r = d.ellipsoid.equatoradius
432 r = _m2radians(distance, r, low=-EPS) # distance=0 from .rhumbMidpointTo
434 a1, b1 = self.philam
435 sb, cb = sincos2(Bearing_(azimuth)) # radians
437 da = r * cb
438 a2 = a1 + da
439 # normalize latitude if past pole
440 if fabs(a2) > PI_2:
441 a2 = _copysign(PI, a2) - a2
443 dp = _logPI_2_2(a2, a1)
444 # q becomes ill-conditioned on E-W course 0/0
445 q = cos(a1) if isnear0(dp) else (da / dp)
446 b2 = b1 if isnear0(q) else (b1 + r * sb / q)
448 h = self._heigHt(height)
449 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h)
450 return r
452 def rhumbDistanceTo(self, other, radius=R_M, exact=False, wrap=False):
453 '''Return the distance from this to an other point along
454 a rhumb line (loxodrome).
456 @arg other: The other point (spherical C{LatLon}).
457 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
458 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if
459 C{B{exact}=True}.
460 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),
461 default C{False} for backward compatibility.
462 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
463 B{C{other}} point (C{bool}).
465 @return: Distance (C{meter}, the same units as B{C{radius}}
466 or C{radians} if B{C{radius}} is C{None}).
468 @raise TypeError: The B{C{other}} point is incompatible.
470 @raise ValueError: Invalid B{C{radius}}.
472 @example:
474 >>> p = LatLon(51.127, 1.338)
475 >>> q = LatLon(50.964, 1.853)
476 >>> d = p.rhumbDistanceTo(q) # 403100
477 '''
478 if exact: # use series, always
479 r = LatLonBase.rhumbDistanceTo(self, other, exact=False, # Krüger
480 radius=radius, wrap=wrap)
481 if radius is None: # angular distance in radians
482 r = r / self._radius # /= chokes PyChecker
483 else:
484 # see <https://www.EdWilliams.org/avform.htm#Rhumb>
485 r, _, _ = self._rhumbs3(other, wrap, r=True)
486 if radius is not None:
487 r *= Radius(radius)
488 return r
490 def rhumbIntersecant2(self, circle, point, other, radius=R_M, exact=True, # PYCHOK signature
491 height=None, wrap=False):
492 '''Compute the intersections of a circle and a rhumb line given as two
493 points and as a point and azimuth.
495 @arg circle: Radius of the circle centered at this location (C{meter},
496 same units as B{C{radius}}) or a point on the circle
497 (this C{LatLon}).
498 @arg point: The rhumb line's start point (this C{LatLon}).
499 @arg other: An other point (this I{on} C{LatLon}) or the azimuth I{of}
500 (compass C{degrees}) the rhumb line.
501 @kwarg radius: Mean earth radius (C{meter}, conventionally).
502 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
503 destination and distance, if C{False} use the basic
504 rhumb methods (C{bool}) or if C{None} use the I{great
505 circle} methods.
506 @kwarg height: Optional height for the intersection points (C{meter},
507 conventionally) or C{None}.
508 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the points
509 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}).
511 @return: 2-Tuple of the intersection points (representing a chord),
512 each an instance of this class. For a tangent line, both
513 points are the same instance, wrapped or I{normalized}.
515 @raise IntersectionError: The circle and line do not intersect.
517 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}
518 or B{C{other}} invalid.
520 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
521 B{C{exact}} or B{C{height}}.
522 '''
523 m = LatLonBase.rhumbIntersecant2 if exact else \
524 LatLonSphericalBase.intersecant2
525 return m(self, circle, point, other, radius=radius, exact=exact,
526 height=height, wrap=wrap)
528 def rhumbMidpointTo(self, other, height=None, radius=R_M, exact=False,
529 fraction=_0_5, wrap=False):
530 '''Return the (loxodromic) midpoint on the rhumb line between
531 this and an other point.
533 @arg other: The other point (spherical LatLon).
534 @kwarg height: Optional height, overriding the mean height (C{meter}).
535 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},
536 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).
537 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),
538 default C{False} for backward compatibility.
539 @kwarg fraction: Midpoint location from this point (C{scalar}), may
540 be negative if C{B{exact}=True}.
541 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{other}}
542 point (C{bool}).
544 @return: The (mid)point at the given B{C{fraction}} along the rhumb
545 line (spherical C{LatLon}).
547 @raise TypeError: The B{C{other}} point is incompatible.
549 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}
551 @example:
553 >>> p = LatLon(51.127, 1.338)
554 >>> q = LatLon(50.964, 1.853)
555 >>> m = p.rhumb_midpointTo(q)
556 >>> m.toStr() # '51.0455°N, 001.5957°E'
557 '''
558 if exact: # use series, always
559 r = LatLonBase.rhumbMidpointTo(self, other, exact=False, # Krüger
560 radius=radius, height=height,
561 fraction=fraction, wrap=wrap)
562 elif fraction is not _0_5:
563 f = Scalar_(fraction=fraction) # low=_0_0
564 r, db, dp = self._rhumbs3(other, wrap, r=True) # radians
565 z = atan2b(db, dp)
566 h = self._havg(other, f=f, h=height)
567 r = self.rhumbDestination(r * f, z, radius=None, height=h)
569 else: # for backward compatibility, unwrapped
570 # see <https://MathForum.org/library/drmath/view/51822.html>
571 a1, b1 = self.philam
572 a2, b2 = self.others(other).philam
574 if fabs(b2 - b1) > PI:
575 b1 += PI2 # crossing anti-meridian
577 a3 = favg(a1, a2)
578 b3 = favg(b1, b2)
580 f1 = tanPI_2_2(a1)
581 if isnon0(f1):
582 f2 = tanPI_2_2(a2)
583 f = f2 / f1
584 if isnon0(f):
585 f = log(f)
586 if isnon0(f):
587 f3 = tanPI_2_2(a3)
588 b3 = fdot(map1(log, f1, f2, f3),
589 -b2, b1, b2 - b1) / f
591 d = self.datum if radius in (None, self._radius) else \
592 _spherical_datum(radius, name=self.name, raiser=_radius_)
593 h = self._havg(other, h=height)
594 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h)
595 return r
597 @property_RO
598 def sphericalLatLon(self):
599 '''Get this C{LatLon}'s spherical class.
600 '''
601 return type(self)
603 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature
604 '''Convert this point to C{Nvector} components, I{including
605 height}.
607 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}
608 keyword arguments, ignored if
609 C{B{Nvector} is None}.
611 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}
612 if B{C{Nvector}} is C{None}.
614 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.
615 '''
616 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)
619def _intersecant2(c, r, p, b, radius=R_M, exact=False, height=None, wrap=False):
620 # (INTERNAL) Intersect a circle and line, see L{intersecant2}
621 # above, separated to allow callers to embellish any exceptions
623 if wrap:
624 p = _unrollon(c, p, wrap=wrap)
625 nonexact = exact is None
627 if not isinstanceof(r, c.__class__, p.__class__):
628 r = Radius_(circle=r)
629 elif nonexact:
630 r = c.distanceTo(r, radius=radius, wrap=wrap)
631 elif isbool(exact):
632 r = c.rhumbDistanceTo(r, radius=radius, exact=exact, wrap=wrap)
633 else:
634 raise _ValueError(exact=exact)
636 if not isinstanceof(b, c.__class__, p.__class__):
637 b = Bearing(b)
638 elif nonexact:
639 b = p.initialBearingTo(b, wrap=wrap)
640 else:
641 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap,
642 b360=True)
644 d = p.distanceTo(c, radius=radius) if nonexact else \
645 p.rhumbDistanceTo(c, radius=radius, exact=exact)
646 if d > EPS0:
647 n = _xattr(c, napieradius=0)
648 a = p.initialBearingTo(c) if nonexact else \
649 p.rhumbAzimuthTo(c, radius=radius, exact=exact, b360=True)
650 s, c = sincos2d(b - a) # Napier's sin(A), cos(A)
651 if r > n:
652 # Napier's right spherical triangle rules (R2) and (R1)
653 # <https://WikiPedia.org/wiki/Spherical_trigonometry>
654 m = p._mpr(radius=radius, exact=exact) # meter per radian
655 if fabs(c) > EPS0:
656 d = d / m # /= chokes PyChecker
657 a = asin1(sin(d) * fabs(s)) # Napier's a
658 c = _copysign(cos(a), c)
659 d = acos1(cos(d) / c) * m
660 a *= m # meter
661 else: # point and chord center coincident
662 a, d = d, 0
663 c = cos(a / m)
664 h = (acos1(cos(r / m) / c) * m) if a < r else 0
665 else: # distance from the chord center to ...
666 a = fabs(d * s) # ... the cicle center ...
667 d *= c # ... and to the point
668 h = sqrt_a(r, a) if a < r else 0 # half chord length
669 if a > r:
670 raise IntersectionError(_too_(Fmt.distant(a)))
671 else:
672 d, h = 0, r # point and circle center coincident
674 _intersecant1, kwds = (p.destination, {}) if nonexact else \
675 (p.rhumbDestination, dict(exact=exact))
676 kwds.update(radius=radius, height=height)
677 t = (_intersecant1(d + h, b, **kwds),)
678 if h:
679 t += (_intersecant1(d - h, b, **kwds),)
680 else: # same instance twice
681 t *= 2
682 return t
685def _logPI_2_2(a2, a1):
686 '''(INTERNAL) C{log} of C{tanPI_2_2}'s quotient.
687 '''
688 return log(_over(tanPI_2_2(a2), tanPI_2_2(a1)))
691def _m2radians(distance, radius, low=EPS): # PYCHOK in .spherical*
692 '''(INTERNAL) Distance in C{meter} to angular distance in C{radians}.
694 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}.
695 '''
696 r = float(distance)
697 if radius:
698 r = r / Radius_(radius=radius) # /= chokes PyChecker
699 if low is not None:
700 # small near0 values from .rhumbDestination not exact OK
701 r = _0_0 if low < 0 and r < 0 else Radians_(r, low=low)
702 # _0_0 if low < 0 and low < r < 0 else Radians_(r, low=low)
703 return r
706def _radians2m(rad, radius):
707 '''(INTERNAL) Angular distance in C{radians} to distance in C{meter}.
708 '''
709 if radius is not None: # not in (None, _0_0)
710 rad *= R_M if radius is R_M else Radius(radius)
711 return rad
714def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry
715 '''(INTERNAL) Convert radii to radians.
716 '''
717 r1 = Radius_(rad1=rad1)
718 r2 = Radius_(rad2=rad2)
719 if radius is not None: # convert radii to radians
720 r1 = _m2radians(r1, radius)
721 r2 = _m2radians(r2, radius)
723 x = r1 < r2
724 if x:
725 r1, r2 = r2, r1
726 if r1 > PI:
727 raise IntersectionError(rad1=rad1, rad2=rad2,
728 txt=_exceed_PI_radians_)
729 return r1, r2, x
732__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase)
734# **) MIT License
735#
736# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
737#
738# Permission is hereby granted, free of charge, to any person obtaining a
739# copy of this software and associated documentation files (the "Software"),
740# to deal in the Software without restriction, including without limitation
741# the rights to use, copy, modify, merge, publish, distribute, sublicense,
742# and/or sell copies of the Software, and to permit persons to whom the
743# Software is furnished to do so, subject to the following conditions:
744#
745# The above copyright notice and this permission notice shall be included
746# in all copies or substantial portions of the Software.
747#
748# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
749# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
750# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
751# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
752# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
753# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
754# OTHER DEALINGS IN THE SOFTWARE.