Coverage for pygeodesy/formy.py: 99%
421 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-23 12:10 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-23 12:10 -0400
2# -*- coding: utf-8 -*-
4u'''Formulary of basic geodesy functions and approximations.
5'''
6# make sure int/int division yields float quotient, see .basics
7from __future__ import division as _; del _ # PYCHOK semicolon
9# from pygeodesy.basics import isscalar # from .fsums
10from pygeodesy.constants import EPS, EPS0, EPS1, PI, PI2, PI3, PI_2, R_M, \
11 _umod_PI2, float0_, isnon0, remainder, \
12 _0_0, _0_125, _0_25, _0_5, _1_0, _2_0, \
13 _4_0, _32_0, _90_0, _180_0, _360_0
14from pygeodesy.datums import Datum, Ellipsoid, _ellipsoidal_datum, \
15 _mean_radius, _spherical_datum, _WGS84, _EWGS84
16# from pygeodesy.ellipsoids import Ellipsoid, _EWGS84 # from .datums
17from pygeodesy.errors import IntersectionError, LimitError, limiterrors, \
18 _TypeError, _ValueError, _xattr, _xError, \
19 _xkwds, _xkwds_pop
20from pygeodesy.fmath import euclid, hypot, hypot2, sqrt0
21from pygeodesy.fsums import fsumf_, isscalar
22from pygeodesy.interns import NN, _delta_, _distant_, _SPACE_, _too_
23from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
24from pygeodesy.named import _NamedTuple, _xnamed, Fmt, unstr
25from pygeodesy.namedTuples import Bearing2Tuple, Distance4Tuple, \
26 Intersection3Tuple, LatLon2Tuple, \
27 PhiLam2Tuple, Vector3Tuple
28# from pygeodesy.streprs import Fmt, unstr # from .named
29from pygeodesy.units import Bearing, Degrees_, Distance, Distance_, Height, \
30 Lam_, Lat, Lon, Meter_, Phi_, Radians, Radians_, \
31 Radius, Radius_, Scalar, _100km
32from pygeodesy.utily import acos1, atan2b, atan2d, degrees2m, m2degrees, \
33 tan_2, sincos2, sincos2_, sincos2d_, _Wrap
35from contextlib import contextmanager
36from math import asin, atan, atan2, cos, degrees, fabs, radians, sin, sqrt # pow
38__all__ = _ALL_LAZY.formy
39__version__ = '23.08.20'
41_D2_R2 = (PI / _180_0)**2 # degrees- to radians-squared
42_ratio_ = 'ratio'
43_xline_ = 'xline'
46def _anti2(a, b, n_2, n, n2):
47 '''(INTERNAL) Helper for C{antipode} and C{antipode_}.
48 '''
49 r = remainder(a, n) if fabs(a) > n_2 else a
50 if r == a:
51 r = -r
52 b += n
53 if fabs(b) > n:
54 b = remainder(b, n2)
55 return float0_(r, b)
58def antipode(lat, lon, name=NN):
59 '''Return the antipode, the point diametrically opposite
60 to a given point in C{degrees}.
62 @arg lat: Latitude (C{degrees}).
63 @arg lon: Longitude (C{degrees}).
64 @kwarg name: Optional name (C{str}).
66 @return: A L{LatLon2Tuple}C{(lat, lon)}.
68 @see: Functions L{antipode_} and L{normal} and U{Geosphere
69 <https://CRAN.R-Project.org/web/packages/geosphere/geosphere.pdf>}.
70 '''
71 return LatLon2Tuple(*_anti2(lat, lon, _90_0, _180_0, _360_0), name=name)
74def antipode_(phi, lam, name=NN):
75 '''Return the antipode, the point diametrically opposite
76 to a given point in C{radians}.
78 @arg phi: Latitude (C{radians}).
79 @arg lam: Longitude (C{radians}).
80 @kwarg name: Optional name (C{str}).
82 @return: A L{PhiLam2Tuple}C{(phi, lam)}.
84 @see: Functions L{antipode} and L{normal_} and U{Geosphere
85 <https://CRAN.R-Project.org/web/packages/geosphere/geosphere.pdf>}.
86 '''
87 return PhiLam2Tuple(*_anti2(phi, lam, PI_2, PI, PI2), name=name)
90def bearing(lat1, lon1, lat2, lon2, **final_wrap):
91 '''Compute the initial or final bearing (forward or reverse
92 azimuth) between a (spherical) start and end point.
94 @arg lat1: Start latitude (C{degrees}).
95 @arg lon1: Start longitude (C{degrees}).
96 @arg lat2: End latitude (C{degrees}).
97 @arg lon2: End longitude (C{degrees}).
98 @kwarg final_wrap: Optional keyword arguments for function
99 L{pygeodesy.bearing_}.
101 @return: Initial or final bearing (compass C{degrees360}) or
102 zero if start and end point coincide.
103 '''
104 r = bearing_(Phi_(lat1=lat1), Lam_(lon1=lon1),
105 Phi_(lat2=lat2), Lam_(lon2=lon2), **final_wrap)
106 return degrees(r)
109def bearing_(phi1, lam1, phi2, lam2, final=False, wrap=False):
110 '''Compute the initial or final bearing (forward or reverse azimuth)
111 between a (spherical) start and end point.
113 @arg phi1: Start latitude (C{radians}).
114 @arg lam1: Start longitude (C{radians}).
115 @arg phi2: End latitude (C{radians}).
116 @arg lam2: End longitude (C{radians}).
117 @kwarg final: Return final bearing if C{True}, initial otherwise (C{bool}).
118 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{phi2}} and
119 B{C{lam2}} (C{bool}).
121 @return: Initial or final bearing (compass C{radiansPI2}) or zero if start
122 and end point coincide.
124 @see: U{Bearing<https://www.Movable-Type.co.UK/scripts/latlong.html>}, U{Course
125 between two points<https://www.EdWilliams.org/avform147.htm#Crs>} and
126 U{Bearing Between Two Points<https://web.Archive.org/web/20020630205931/
127 https://MathForum.org/library/drmath/view/55417.html>}.
128 '''
129 db, phi2, lam2 = _Wrap.philam3(lam1, phi2, lam2, wrap)
130 if final: # swap plus PI
131 phi1, lam1, phi2, lam2, db = phi2, lam2, phi1, lam1, -db
132 r = PI3
133 else:
134 r = PI2
135 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(phi1, phi2, db)
137 x = ca1 * sa2 - sa1 * ca2 * cdb
138 y = sdb * ca2
139 return _umod_PI2(atan2(y, x) + r) # .utily.wrapPI2
142def _bearingTo2(p1, p2, wrap=False): # for points.ispolar, sphericalTrigonometry.areaOf
143 '''(INTERNAL) Compute initial and final bearing.
144 '''
145 try: # for LatLon_ and ellipsoidal LatLon
146 return p1.bearingTo2(p2, wrap=wrap)
147 except AttributeError:
148 pass
149 # XXX spherical version, OK for ellipsoidal ispolar?
150 a1, b1 = p1.philam
151 a2, b2 = p2.philam
152 return Bearing2Tuple(degrees(bearing_(a1, b1, a2, b2, final=False, wrap=wrap)),
153 degrees(bearing_(a1, b1, a2, b2, final=True, wrap=wrap)),
154 name=_bearingTo2.__name__)
157def compassAngle(lat1, lon1, lat2, lon2, adjust=True, wrap=False):
158 '''Return the angle from North for the direction vector M{(lon2 - lon1,
159 lat2 - lat1)} between two points.
161 Suitable only for short, not near-polar vectors up to a few hundred
162 Km or Miles. Use function L{pygeodesy.bearing} for longer vectors.
164 @arg lat1: From latitude (C{degrees}).
165 @arg lon1: From longitude (C{degrees}).
166 @arg lat2: To latitude (C{degrees}).
167 @arg lon2: To longitude (C{degrees}).
168 @kwarg adjust: Adjust the longitudinal delta by the cosine of the
169 mean latitude (C{bool}).
170 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
171 and B{C{lon2}} (C{bool}).
173 @return: Compass angle from North (C{degrees360}).
175 @note: Courtesy of Martin Schultz.
177 @see: U{Local, flat earth approximation
178 <https://www.EdWilliams.org/avform.htm#flat>}.
179 '''
180 d_lon, lat2, lon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)
181 if adjust: # scale delta lon
182 d_lon *= _scale_deg(lat1, lat2)
183 return atan2b(d_lon, lat2 - lat1)
186def cosineAndoyerLambert(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):
187 '''Compute the distance between two (ellipsoidal) points using the
188 U{Andoyer-Lambert correction<https://NavLib.net/wp-content/uploads/2013/10/
189 admiralty-manual-of-navigation-vol-1-1964-english501c.pdf>} of the U{Law of
190 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.
192 @arg lat1: Start latitude (C{degrees}).
193 @arg lon1: Start longitude (C{degrees}).
194 @arg lat2: End latitude (C{degrees}).
195 @arg lon2: End longitude (C{degrees}).
196 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
197 L{Ellipsoid2} or L{a_f2Tuple}) to use.
198 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
199 B{C{lat2}} and B{C{lon2}} (C{bool}).
201 @return: Distance (C{meter}, same units as the B{C{datum}}'s
202 ellipsoid axes or C{radians} if B{C{datum}} is C{None}).
204 @raise TypeError: Invalid B{C{datum}}.
206 @see: Functions L{cosineAndoyerLambert_}, L{cosineForsytheAndoyerLambert},
207 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},
208 L{flatPolar}, L{haversine}, L{thomas} and L{vincentys} and method
209 L{Ellipsoid.distance2}.
210 '''
211 return _dE(cosineAndoyerLambert_, datum, wrap, lat1, lon1, lat2, lon2)
214def cosineAndoyerLambert_(phi2, phi1, lam21, datum=_WGS84):
215 '''Compute the I{angular} distance between two (ellipsoidal) points using the
216 U{Andoyer-Lambert correction<https://NavLib.net/wp-content/uploads/2013/10/
217 admiralty-manual-of-navigation-vol-1-1964-english501c.pdf>} of the U{Law of
218 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.
220 @arg phi2: End latitude (C{radians}).
221 @arg phi1: Start latitude (C{radians}).
222 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
223 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
224 L{Ellipsoid2} or L{a_f2Tuple}) to use.
226 @return: Angular distance (C{radians}).
228 @raise TypeError: Invalid B{C{datum}}.
230 @see: Functions L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert_},
231 L{cosineLaw_}, L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
232 L{flatPolar_}, L{haversine_}, L{thomas_} and L{vincentys_} and U{Geodesy-PHP
233 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/Distance/
234 AndoyerLambert.php>}.
235 '''
236 s2, c2, s1, c1, r, c21 = _sincosa6(phi2, phi1, lam21)
237 if isnon0(c1) and isnon0(c2):
238 E = _ellipsoidal(datum, cosineAndoyerLambert_)
239 if E.f: # ellipsoidal
240 r2 = atan2(E.b_a * s2, c2)
241 r1 = atan2(E.b_a * s1, c1)
242 s2, c2, s1, c1 = sincos2_(r2, r1)
243 r = acos1(s1 * s2 + c1 * c2 * c21)
244 if r:
245 sr, _, sr_2, cr_2 = sincos2_(r, r * _0_5)
246 if isnon0(sr_2) and isnon0(cr_2):
247 s = (sr + r) * ((s1 - s2) / sr_2)**2
248 c = (sr - r) * ((s1 + s2) / cr_2)**2
249 r += (c - s) * E.f * _0_125
250 return r
253def cosineForsytheAndoyerLambert(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):
254 '''Compute the distance between two (ellipsoidal) points using the
255 U{Forsythe-Andoyer-Lambert correction<https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} of
256 the U{Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
257 formula.
259 @arg lat1: Start latitude (C{degrees}).
260 @arg lon1: Start longitude (C{degrees}).
261 @arg lat2: End latitude (C{degrees}).
262 @arg lon2: End longitude (C{degrees}).
263 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
264 L{Ellipsoid2} or L{a_f2Tuple}) to use.
265 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
266 B{C{lat2}} and B{C{lon2}} (C{bool}).
268 @return: Distance (C{meter}, same units as the B{C{datum}}'s
269 ellipsoid axes or C{radians} if B{C{datum}} is C{None}).
271 @raise TypeError: Invalid B{C{datum}}.
273 @see: Functions L{cosineForsytheAndoyerLambert_}, L{cosineAndoyerLambert},
274 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},
275 L{flatPolar}, L{haversine}, L{thomas} and L{vincentys} and method
276 L{Ellipsoid.distance2}.
277 '''
278 return _dE(cosineForsytheAndoyerLambert_, datum, wrap, lat1, lon1, lat2, lon2)
281def cosineForsytheAndoyerLambert_(phi2, phi1, lam21, datum=_WGS84):
282 '''Compute the I{angular} distance between two (ellipsoidal) points using the
283 U{Forsythe-Andoyer-Lambert correction<https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} of
284 the U{Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
285 formula.
287 @arg phi2: End latitude (C{radians}).
288 @arg phi1: Start latitude (C{radians}).
289 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
290 @kwarg datum: Datum (L{Datum}) or ellipsoid to use (L{Ellipsoid},
291 L{Ellipsoid2} or L{a_f2Tuple}).
293 @return: Angular distance (C{radians}).
295 @raise TypeError: Invalid B{C{datum}}.
297 @see: Functions L{cosineForsytheAndoyerLambert}, L{cosineAndoyerLambert_},
298 L{cosineLaw_}, L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
299 L{flatPolar_}, L{haversine_}, L{thomas_} and L{vincentys_} and U{Geodesy-PHP
300 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/
301 Distance/ForsytheCorrection.php>}.
302 '''
303 s2, c2, s1, c1, r, _ = _sincosa6(phi2, phi1, lam21)
304 if r and isnon0(c1) and isnon0(c2):
305 E = _ellipsoidal(datum, cosineForsytheAndoyerLambert_)
306 if E.f: # ellipsoidal
307 sr, cr, s2r, _ = sincos2_(r, r * 2)
308 if isnon0(sr) and fabs(cr) < EPS1:
309 s = (s1 + s2)**2 / (1 + cr)
310 t = (s1 - s2)**2 / (1 - cr)
311 x = s + t
312 y = s - t
314 s = 8 * r**2 / sr
315 a = 64 * r + s * cr * 2 # 16 * r**2 / tan(r)
316 d = 48 * sr + s # 8 * r**2 / tan(r)
317 b = -2 * d
318 e = 30 * s2r
319 c = fsumf_(30 * r, e * _0_5, s * cr) # 8 * r**2 / tan(r)
320 t = fsumf_( a * x, e * y**2, b * y, -c * x**2, d * x * y)
322 r += fsumf_(-r * x, 3 * y * sr, t * E.f / _32_0) * E.f * _0_25
323 return r
326def cosineLaw(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
327 '''Compute the distance between two points using the U{spherical Law of
328 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
329 formula.
331 @arg lat1: Start latitude (C{degrees}).
332 @arg lon1: Start longitude (C{degrees}).
333 @arg lat2: End latitude (C{degrees}).
334 @arg lon2: End longitude (C{degrees}).
335 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
336 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
337 L{a_f2Tuple}) to use.
338 @kwarg wrap: If C{True}, wrap or I{normalize} and B{C{lat2}}
339 and B{C{lon2}} (C{bool}).
341 @return: Distance (C{meter}, same units as B{C{radius}} or the
342 ellipsoid or datum axes).
344 @raise TypeError: Invalid B{C{radius}}.
346 @see: Functions L{cosineLaw_}, L{cosineAndoyerLambert},
347 L{cosineForsytheAndoyerLambert}, L{equirectangular}, L{euclidean},
348 L{flatLocal}/L{hubeny}, L{flatPolar}, L{haversine}, L{thomas} and
349 L{vincentys} and method L{Ellipsoid.distance2}.
351 @note: See note at function L{vincentys_}.
352 '''
353 return _dS(cosineLaw_, radius, wrap, lat1, lon1, lat2, lon2)
356def cosineLaw_(phi2, phi1, lam21):
357 '''Compute the I{angular} distance between two points using the U{spherical
358 Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
359 formula.
361 @arg phi2: End latitude (C{radians}).
362 @arg phi1: Start latitude (C{radians}).
363 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
365 @return: Angular distance (C{radians}).
367 @see: Functions L{cosineLaw}, L{cosineAndoyerLambert_},
368 L{cosineForsytheAndoyerLambert_}, L{equirectangular_},
369 L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},
370 L{haversine_}, L{thomas_} and L{vincentys_}.
372 @note: See note at function L{vincentys_}.
373 '''
374 return _sincosa6(phi2, phi1, lam21)[4]
377def _d3(wrap, lat1, lon1, lat2, lon2):
378 '''(INTERNAL) Helper for _dE, _dS and _eA.
379 '''
380 if wrap:
381 d_lon, lat2, _ = _Wrap.latlon3(lon1, lat2, lon2, wrap)
382 return radians(lat2), Phi_(lat1=lat1), radians(d_lon)
383 else: # for backward compaibility
384 return Phi_(lat2=lat2), Phi_(lat1=lat1), Phi_(d_lon=lon2 - lon1)
387def _dE(func_, earth, *wrap_lls):
388 '''(INTERNAL) Helper for ellipsoidal distances.
389 '''
390 E = _ellipsoidal(earth, func_)
391 r = func_(*_d3(*wrap_lls), datum=E)
392 return r * E.a
395def _dS(func_, radius, *wrap_lls, **adjust):
396 '''(INTERNAL) Helper for spherical distances.
397 '''
398 r = func_(*_d3(*wrap_lls), **adjust)
399 if radius is not R_M:
400 _, lat1, _, lat2, _ = wrap_lls
401 radius = _mean_radius(radius, lat1, lat2)
402 return r * radius
405def _eA(excess_, radius, *wrap_lls):
406 '''(INTERNAL) Helper for spherical excess or area.
407 '''
408 r = excess_(*_d3(*wrap_lls))
409 if radius:
410 _, lat1, _, lat2, _ = wrap_lls
411 r *= _mean_radius(radius, lat1, lat2)**2
412 return r
415def _ellipsoidal(earth, where):
416 '''(INTERNAL) Helper for distances.
417 '''
418 return _EWGS84 if earth in (_WGS84, _EWGS84) else (
419 earth if isinstance(earth, Ellipsoid) else
420 (earth if isinstance(earth, Datum) else # PYCHOK indent
421 _ellipsoidal_datum(earth, name=where.__name__)).ellipsoid)
424def equirectangular(lat1, lon1, lat2, lon2, radius=R_M, **adjust_limit_wrap):
425 '''Compute the distance between two points using
426 the U{Equirectangular Approximation / Projection
427 <https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.
429 @arg lat1: Start latitude (C{degrees}).
430 @arg lon1: Start longitude (C{degrees}).
431 @arg lat2: End latitude (C{degrees}).
432 @arg lon2: End longitude (C{degrees}).
433 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
434 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
435 L{a_f2Tuple}).
436 @kwarg adjust_limit_wrap: Optional keyword arguments for
437 function L{equirectangular_}.
439 @return: Distance (C{meter}, same units as B{C{radius}} or
440 the ellipsoid or datum axes).
442 @raise TypeError: Invalid B{C{radius}}.
444 @see: Function L{equirectangular_} for more details, the
445 available B{C{options}}, errors, restrictions and other,
446 approximate or accurate distance functions.
447 '''
448 d = sqrt(equirectangular_(Lat(lat1=lat1), Lon(lon1=lon1),
449 Lat(lat2=lat2), Lon(lon2=lon2),
450 **adjust_limit_wrap).distance2) # PYCHOK 4 vs 2-3
451 return degrees2m(d, radius=_mean_radius(radius, lat1, lat2))
454def _equirectangular(lat1, lon1, lat2, lon2, **adjust_limit_wrap):
455 '''(INTERNAL) Helper for the L{frechet._FrecherMeterRadians}
456 and L{hausdorff._HausdorffMeterRedians} classes.
457 '''
458 return equirectangular_(lat1, lon1, lat2, lon2, **adjust_limit_wrap).distance2 * _D2_R2
461def equirectangular_(lat1, lon1, lat2, lon2, adjust=True, limit=45, wrap=False):
462 '''Compute the distance between two points using the U{Equirectangular
463 Approximation / Projection
464 <https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.
466 This approximation is valid for short distance of several hundred Km
467 or Miles, see the B{C{limit}} keyword argument and L{LimitError}.
469 @arg lat1: Start latitude (C{degrees}).
470 @arg lon1: Start longitude (C{degrees}).
471 @arg lat2: End latitude (C{degrees}).
472 @arg lon2: End longitude (C{degrees}).
473 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
474 by the cosine of the mean latitude (C{bool}).
475 @kwarg limit: Optional limit for lat- and longitudinal deltas
476 (C{degrees}) or C{None} or C{0} for unlimited.
477 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
478 and B{C{lon2}} (C{bool}).
480 @return: A L{Distance4Tuple}C{(distance2, delta_lat, delta_lon,
481 unroll_lon2)} in C{degrees squared}.
483 @raise LimitError: If the lat- and/or longitudinal delta exceeds the
484 B{C{-limit..limit}} range and L{pygeodesy.limiterrors}
485 set to C{True}.
487 @see: U{Local, flat earth approximation
488 <https://www.EdWilliams.org/avform.htm#flat>}, functions
489 L{equirectangular}, L{cosineAndoyerLambert},
490 L{cosineForsytheAndoyerLambert}, L{cosineLaw}, L{euclidean},
491 L{flatLocal}/L{hubeny}, L{flatPolar}, L{haversine}, L{thomas}
492 and L{vincentys} and methods L{Ellipsoid.distance2},
493 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.
494 '''
495 d_lon, lat2, ulon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)
496 d_lat = lat2 - lat1
498 if limit and limit > 0 and limiterrors():
499 d = max(fabs(d_lat), fabs(d_lon))
500 if d > limit:
501 t = _SPACE_(_delta_, Fmt.PAREN_g(d), Fmt.exceeds_limit(limit))
502 s = unstr(equirectangular_, lat1, lon1, lat2, lon2,
503 limit=limit, wrap=wrap)
504 raise LimitError(s, txt=t)
506 if adjust: # scale delta lon
507 d_lon *= _scale_deg(lat1, lat2)
509 d2 = hypot2(d_lat, d_lon) # degrees squared!
510 return Distance4Tuple(d2, d_lat, d_lon, ulon2 - lon2)
513def euclidean(lat1, lon1, lat2, lon2, radius=R_M, adjust=True, wrap=False):
514 '''Approximate the C{Euclidean} distance between two (spherical) points.
516 @arg lat1: Start latitude (C{degrees}).
517 @arg lon1: Start longitude (C{degrees}).
518 @arg lat2: End latitude (C{degrees}).
519 @arg lon2: End longitude (C{degrees}).
520 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
521 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
522 L{a_f2Tuple}) to use.
523 @kwarg adjust: Adjust the longitudinal delta by the cosine of
524 the mean latitude (C{bool}).
525 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
526 and B{C{lon2}} (C{bool}).
528 @return: Distance (C{meter}, same units as B{C{radius}} or the
529 ellipsoid or datum axes).
531 @raise TypeError: Invalid B{C{radius}}.
533 @see: U{Distance between two (spherical) points
534 <https://www.EdWilliams.org/avform.htm#Dist>}, functions L{euclid},
535 L{euclidean_}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},
536 L{cosineLaw}, L{equirectangular}, L{flatLocal}/L{hubeny}, L{flatPolar},
537 L{haversine}, L{thomas} and L{vincentys} and methods L{Ellipsoid.distance2},
538 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.
539 '''
540 return _dS(euclidean_, radius, wrap, lat1, lon1, lat2, lon2, adjust=adjust)
543def euclidean_(phi2, phi1, lam21, adjust=True):
544 '''Approximate the I{angular} C{Euclidean} distance between two
545 (spherical) points.
547 @arg phi2: End latitude (C{radians}).
548 @arg phi1: Start latitude (C{radians}).
549 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
550 @kwarg adjust: Adjust the longitudinal delta by the cosine
551 of the mean latitude (C{bool}).
553 @return: Angular distance (C{radians}).
555 @see: Functions L{euclid}, L{euclidean}, L{cosineAndoyerLambert_},
556 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_}, L{equirectangular_},
557 L{flatLocal_}/L{hubeny_}, L{flatPolar_}, L{haversine_}, L{thomas_}
558 and L{vincentys_}.
559 '''
560 if adjust:
561 lam21 *= _scale_rad(phi2, phi1)
562 return euclid(phi2 - phi1, lam21)
565def excessAbc_(A, b, c):
566 '''Compute the I{spherical excess} C{E} of a (spherical) triangle
567 from two sides and the included (small) angle.
569 @arg A: An interior triangle angle (C{radians}).
570 @arg b: Frist adjacent triangle side (C{radians}).
571 @arg c: Second adjacent triangle side (C{radians}).
573 @return: Spherical excess (C{radians}).
575 @raise UnitError: Invalid B{C{A}}, B{C{b}} or B{C{c}}.
577 @see: Functions L{excessGirard_}, L{excessLHuilier_} and U{Spherical
578 trigonometry<https://WikiPedia.org/wiki/Spherical_trigonometry>}.
579 '''
580 A = Radians_(A=A)
581 b = Radians_(b=b) * _0_5
582 c = Radians_(c=c) * _0_5
584 sA, cA, sb, cb, sc, cc = sincos2_(A, b, c)
585 return atan2(sA * sb * sc, cb * cc + cA * sb * sc) * _2_0
588def excessCagnoli_(a, b, c):
589 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using
590 U{Cagnoli's<https://Zenodo.org/record/35392>} (D.34) formula.
592 @arg a: First triangle side (C{radians}).
593 @arg b: Second triangle side (C{radians}).
594 @arg c: Third triangle side (C{radians}).
596 @return: Spherical excess (C{radians}).
598 @raise UnitError: Invalid B{C{a}}, B{C{b}} or B{C{c}}.
600 @see: Function L{excessLHuilier_} and U{Spherical trigonometry
601 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
602 '''
603 a = Radians_(a=a)
604 b = Radians_(b=b)
605 c = Radians_(c=c)
607 s = fsumf_(a, b, c) * _0_5
608 r = sin(s) * sin(s - a) * sin(s - b) * sin(s - c)
609 c = cos(a * _0_5) * cos(b * _0_5) * cos(c * _0_5)
610 r = asin(sqrt(r) * _0_5 / c) if c and r > 0 else _0_0
611 return Radians(Cagnoli=r * _2_0)
614def excessGirard_(A, B, C):
615 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using
616 U{Girard's<https://MathWorld.Wolfram.com/GirardsSphericalExcessFormula.html>}
617 formula.
619 @arg A: First interior triangle angle (C{radians}).
620 @arg B: Second interior triangle angle (C{radians}).
621 @arg C: Third interior triangle angle (C{radians}).
623 @return: Spherical excess (C{radians}).
625 @raise UnitError: Invalid B{C{A}}, B{C{B}} or B{C{C}}.
627 @see: Function L{excessLHuilier_} and U{Spherical trigonometry
628 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
629 '''
630 return Radians(Girard=fsumf_(Radians_(A=A),
631 Radians_(B=B),
632 Radians_(C=C), -PI))
635def excessLHuilier_(a, b, c):
636 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using
637 U{L'Huilier's<https://MathWorld.Wolfram.com/LHuiliersTheorem.html>}
638 Theorem.
640 @arg a: First triangle side (C{radians}).
641 @arg b: Second triangle side (C{radians}).
642 @arg c: Third triangle side (C{radians}).
644 @return: Spherical excess (C{radians}).
646 @raise UnitError: Invalid B{C{a}}, B{C{b}} or B{C{c}}.
648 @see: Function L{excessCagnoli_}, L{excessGirard_} and U{Spherical
649 trigonometry<https://WikiPedia.org/wiki/Spherical_trigonometry>}.
650 '''
651 a = Radians_(a=a)
652 b = Radians_(b=b)
653 c = Radians_(c=c)
655 s = fsumf_(a, b, c) * _0_5
656 r = tan_2(s) * tan_2(s - a) * tan_2(s - b) * tan_2(s - c)
657 r = atan(sqrt(r)) if r > 0 else _0_0
658 return Radians(LHuilier=r * _4_0)
661def excessKarney(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
662 '''Compute the surface area of a (spherical) quadrilateral bounded by a
663 segment of a great circle, two meridians and the equator using U{Karney's
664 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}
665 method.
667 @arg lat1: Start latitude (C{degrees}).
668 @arg lon1: Start longitude (C{degrees}).
669 @arg lat2: End latitude (C{degrees}).
670 @arg lon2: End longitude (C{degrees}).
671 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
672 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
673 L{a_f2Tuple}) or C{None}.
674 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
675 B{C{lat2}} and B{C{lon2}} (C{bool}).
677 @return: Surface area, I{signed} (I{square} C{meter} or the same units as
678 B{C{radius}} I{squared}) or the I{spherical excess} (C{radians})
679 if C{B{radius}=0} or C{None}.
681 @raise TypeError: Invalid B{C{radius}}.
683 @raise UnitError: Invalid B{C{lat2}} or B{C{lat1}}.
685 @raise ValueError: Semi-circular longitudinal delta.
687 @see: Functions L{excessKarney_} and L{excessQuad}.
688 '''
689 return _eA(excessKarney_, radius, wrap, lat1, lon1, lat2, lon2)
692def excessKarney_(phi2, phi1, lam21):
693 '''Compute the I{spherical excess} C{E} of a (spherical) quadrilateral bounded
694 by a segment of a great circle, two meridians and the equator using U{Karney's
695 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}
696 method.
698 @arg phi2: End latitude (C{radians}).
699 @arg phi1: Start latitude (C{radians}).
700 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
702 @return: Spherical excess, I{signed} (C{radians}).
704 @raise ValueError: Semi-circular longitudinal delta B{C{lam21}}.
706 @see: Function L{excessKarney} and U{Area of a spherical polygon
707 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}.
708 '''
709 # from: Veness <https://www.Movable-Type.co.UK/scripts/latlong.html> Area
710 # method due to Karney: for each edge of the polygon,
711 #
712 # tan(Δλ / 2) · (tan(φ1 / 2) + tan(φ2 / 2))
713 # tan(E / 2) = -----------------------------------------
714 # 1 + tan(φ1 / 2) · tan(φ2 / 2)
715 #
716 # where E is the spherical excess of the trapezium obtained by extending
717 # the edge to the equator-circle vector for each edge (see also ***).
718 t2 = tan_2(phi2)
719 t1 = tan_2(phi1)
720 t = tan_2(lam21, lam21=None)
721 return Radians(Karney=atan2(t * (t1 + t2),
722 _1_0 + (t1 * t2)) * _2_0)
725# ***) Original post no longer available, following is a copy of the main part
726# <http://OSGeo-org.1560.x6.Nabble.com/Area-of-a-spherical-polygon-td3841625.html>
727#
728# The area of a polygon on a (unit) sphere is given by the spherical excess
729#
730# A = 2 * pi - sum(exterior angles)
731#
732# However this is badly conditioned if the polygon is small. In this case, use
733#
734# A = sum(S12{i, i+1}) over the edges of the polygon
735#
736# where S12 is the area of the quadrilateral bounded by an edge of the polygon,
737# two meridians and the equator, i.e. with vertices (phi1, lambda1), (phi2,
738# lambda2), (0, lambda1) and (0, lambda2). S12 is given by
739#
740# tan(S12 / 2) = tan(lambda21 / 2) * (tan(phi1 / 2) + tan(phi2 / 2)) /
741# (tan(phi1 / 2) * tan(phi2 / 2) + 1)
742#
743# = tan(lambda21 / 2) * tanh((Lamb(phi1) + Lamb(phi2)) / 2)
744#
745# where lambda21 = lambda2 - lambda1 and Lamb(x) is the Lambertian (or the
746# inverse Gudermannian) function
747#
748# Lambertian(x) = asinh(tan(x)) = atanh(sin(x)) = 2 * atanh(tan(x / 2))
749#
750# Notes: The formula for S12 is exact, except that...
751# - it is indeterminate if an edge is a semi-circle
752# - the formula for A applies only if the polygon does not include a pole
753# (if it does, then add +/- 2 * pi to the result)
754# - in the limit of small phi and lambda, S12 reduces to the trapezoidal
755# formula, S12 = (lambda2 - lambda1) * (phi1 + phi2) / 2
756# - I derived this result from the equation for the area of a spherical
757# triangle in terms of two edges and the included angle given by, e.g.
758# U{Todhunter, I. - Spherical Trigonometry (1871), Sec. 103, Eq. (2)
759# <http://Books.Google.com/books?id=3uBHAAAAIAAJ&pg=PA71>}
760# - I would be interested to know if this formula for S12 is already known
761# - Charles Karney
764def excessQuad(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
765 '''Compute the surface area of a (spherical) quadrilateral bounded by a segment
766 of a great circle, two meridians and the equator.
768 @arg lat1: Start latitude (C{degrees}).
769 @arg lon1: Start longitude (C{degrees}).
770 @arg lat2: End latitude (C{degrees}).
771 @arg lon2: End longitude (C{degrees}).
772 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
773 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
774 L{a_f2Tuple}) or C{None}.
775 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
776 B{C{lat2}} and B{C{lon2}} (C{bool}).
778 @return: Surface area, I{signed} (I{square} C{meter} or the same units as
779 B{C{radius}} I{squared}) or the I{spherical excess} (C{radians})
780 if C{B{radius}=0} or C{None}.
782 @raise TypeError: Invalid B{C{radius}}.
784 @raise UnitError: Invalid B{C{lat2}} or B{C{lat1}}.
786 @see: Function L{excessQuad_} and L{excessKarney}.
787 '''
788 return _eA(excessQuad_, radius, wrap, lat1, lon1, lat2, lon2)
791def excessQuad_(phi2, phi1, lam21):
792 '''Compute the I{spherical excess} C{E} of a (spherical) quadrilateral bounded
793 by a segment of a great circle, two meridians and the equator.
795 @arg phi2: End latitude (C{radians}).
796 @arg phi1: Start latitude (C{radians}).
797 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
799 @return: Spherical excess, I{signed} (C{radians}).
801 @see: Function L{excessQuad} and U{Spherical trigonometry
802 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
803 '''
804 s = sin((phi2 + phi1) * _0_5)
805 c = cos((phi2 - phi1) * _0_5)
806 return Radians(Quad=atan2(tan_2(lam21) * s, c) * _2_0)
809def flatLocal(lat1, lon1, lat2, lon2, datum=_WGS84, scaled=True, wrap=False):
810 '''Compute the distance between two (ellipsoidal) points using
811 the U{ellipsoidal Earth to plane projection<https://WikiPedia.org/
812 wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}
813 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.
815 @arg lat1: Start latitude (C{degrees}).
816 @arg lon1: Start longitude (C{degrees}).
817 @arg lat2: End latitude (C{degrees}).
818 @arg lon2: End longitude (C{degrees}).
819 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
820 L{Ellipsoid2} or L{a_f2Tuple}) to use.
821 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}),
822 see method L{pygeodesy.Ellipsoid.roc2_}.
823 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
824 B{C{lat2}} and B{C{lon2}} (C{bool}).
826 @return: Distance (C{meter}, same units as the B{C{datum}}'s
827 ellipsoid axes).
829 @raise TypeError: Invalid B{C{datum}}.
831 @note: The meridional and prime_vertical radii of curvature
832 are taken and scaled at the mean of both latitude.
834 @see: Functions L{flatLocal_} or L{hubeny_}, L{cosineLaw}, L{flatPolar},
835 L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},
836 L{equirectangular}, L{euclidean}, L{haversine}, L{thomas},
837 L{vincentys}, method L{Ellipsoid.distance2} and U{local, flat
838 earth approximation<https://www.EdWilliams.org/avform.htm#flat>}.
839 '''
840 E = _ellipsoidal(datum, flatLocal)
841 return E._hubeny_2(*_d3(wrap, lat1, lon1, lat2, lon2),
842 scaled=scaled, squared=False) * E.a
844hubeny = flatLocal # PYCHOK for Karl Hubeny
847def flatLocal_(phi2, phi1, lam21, datum=_WGS84, scaled=True):
848 '''Compute the I{angular} distance between two (ellipsoidal) points using
849 the U{ellipsoidal Earth to plane projection<https://WikiPedia.org/
850 wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}
851 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.
853 @arg phi2: End latitude (C{radians}).
854 @arg phi1: Start latitude (C{radians}).
855 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
856 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
857 L{Ellipsoid2} or L{a_f2Tuple}) to use.
858 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}),
859 see method L{pygeodesy.Ellipsoid.roc2_}.
861 @return: Angular distance (C{radians}).
863 @raise TypeError: Invalid B{C{datum}}.
865 @note: The meridional and prime_vertical radii of curvature
866 are taken and scaled I{at the mean of both latitude}.
868 @see: Functions L{flatLocal} or L{hubeny}, L{cosineAndoyerLambert_},
869 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_}, L{flatPolar_},
870 L{equirectangular_}, L{euclidean_}, L{haversine_}, L{thomas_}
871 and L{vincentys_} and U{local, flat earth approximation
872 <https://www.EdWilliams.org/avform.htm#flat>}.
873 '''
874 E = _ellipsoidal(datum, flatLocal_)
875 return E._hubeny_2(phi2, phi1, lam21, scaled=scaled, squared=False)
877hubeny_ = flatLocal_ # PYCHOK for Karl Hubeny
880def flatPolar(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
881 '''Compute the distance between two (spherical) points using
882 the U{polar coordinate flat-Earth <https://WikiPedia.org/wiki/
883 Geographical_distance#Polar_coordinate_flat-Earth_formula>}
884 formula.
886 @arg lat1: Start latitude (C{degrees}).
887 @arg lon1: Start longitude (C{degrees}).
888 @arg lat2: End latitude (C{degrees}).
889 @arg lon2: End longitude (C{degrees}).
890 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
891 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
892 L{a_f2Tuple}) to use.
893 @kwarg wrap: If C{True}, wrap or I{normalize} and B{C{lat2}}
894 and B{C{lon2}} (C{bool}).
896 @return: Distance (C{meter}, same units as B{C{radius}} or the
897 ellipsoid or datum axes).
899 @raise TypeError: Invalid B{C{radius}}.
901 @see: Functions L{flatPolar_}, L{cosineAndoyerLambert},
902 L{cosineForsytheAndoyerLambert},L{cosineLaw},
903 L{flatLocal}/L{hubeny}, L{equirectangular},
904 L{euclidean}, L{haversine}, L{thomas} and
905 L{vincentys}.
906 '''
907 return _dS(flatPolar_, radius, wrap, lat1, lon1, lat2, lon2)
910def flatPolar_(phi2, phi1, lam21):
911 '''Compute the I{angular} distance between two (spherical) points
912 using the U{polar coordinate flat-Earth<https://WikiPedia.org/wiki/
913 Geographical_distance#Polar_coordinate_flat-Earth_formula>}
914 formula.
916 @arg phi2: End latitude (C{radians}).
917 @arg phi1: Start latitude (C{radians}).
918 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
920 @return: Angular distance (C{radians}).
922 @see: Functions L{flatPolar}, L{cosineAndoyerLambert_},
923 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
924 L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
925 L{haversine_}, L{thomas_} and L{vincentys_}.
926 '''
927 a = fabs(PI_2 - phi1) # co-latitude
928 b = fabs(PI_2 - phi2) # co-latitude
929 if a < b:
930 a, b = b, a
931 if a < EPS0:
932 a = _0_0
933 elif b > 0:
934 b = b / a # /= chokes PyChecker
935 c = b * cos(lam21) * _2_0
936 c = fsumf_(_1_0, b**2, -fabs(c))
937 a *= sqrt0(c)
938 return a
941def hartzell(pov, los=None, earth=_WGS84, name=NN, **LatLon_and_kwds):
942 '''Compute the intersection of the earth's surface and a Line-Of-Sight
943 from a Point-Of-View in space.
945 @arg pov: Point-Of-View outside the earth (C{Cartesian}, L{Ecef9Tuple}
946 or L{Vector3d}).
947 @kwarg los: Line-Of-Sight, I{direction} to earth (L{Vector3d}) or
948 C{None} to point to the earth' center.
949 @kwarg earth: The earth model (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},
950 L{a_f2Tuple} or C{scalar} radius in C{meter}).
951 @kwarg name: Optional name (C{str}).
952 @kwarg LatLon_and_kwds: Optional C{LatLon} class for the intersection
953 point plus C{LatLon} keyword arguments, include
954 B{C{datum}} if different from B{C{earth}}.
956 @return: The earth intersection (L{Vector3d}, C{Cartesian type} of
957 B{C{pov}} or B{C{LatLon}}).
959 @raise IntersectionError: Null B{C{pov}} or B{C{los}} vector, B{C{pov}}
960 is inside the earth or B{C{los}} points outside
961 the earth or points in an opposite direction.
963 @raise TypeError: Invalid B{C{pov}}, B{C{los}} or B{C{earth}}.
965 @see: Function L{pygeodesy.hartzell4}, L{pygeodesy.tyr3d} for B{C{los}},
966 method L{Ellipsoid.hartzell4} and U{I{Satellite Line-of-Sight
967 Intersection with Earth}<https://StephenHartzell.Medium.com/
968 satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}.
969 '''
970 D = earth if isinstance(earth, Datum) else \
971 _spherical_datum(earth, name=hartzell.__name__)
972 try:
973 r, _ = _MODS.triaxials._hartzell3d2(pov, los, D.ellipsoid._triaxial)
974 except Exception as x:
975 raise IntersectionError(pov=pov, los=los, earth=earth, cause=x)
977# else:
978# E = D.ellipsoid
979# # Triaxial(a, b, c) == (E.a, E.a, E.b)
980#
981# def _Error(txt):
982# return IntersectionError(pov=pov, los=los, earth=earth, txt=txt)
983#
984# a2 = b2 = E.a2 # earth' x, y, ...
985# c2 = E.b2 # ... z semi-axis squared
986# q2 = E.b2_a2 # == c2 / a2
987# bc = E.a * E.b # == b * c
988#
989# V3 = _MODS.vector3d._otherV3d
990# p3 = V3(pov=pov)
991# u3 = V3(los=los) if los else p3.negate()
992# u3 = u3.unit() # unit vector, opposing signs
993#
994# x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz
995# ux, vy, wz = u3.times_(p3).xyz
996# u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz
997#
998# t = c2, c2, b2
999# m = fdot(t, u2, v2, w2) # a2 factored out
1000# if m < EPS0: # zero or near-null LOS vector
1001# raise _Error(_near_(_null_))
1002#
1003# # a2 and b2 factored out, b2 == a2 and b2 / a2 == 1
1004# r = fsumf_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2,
1005# c2 * u2, -u2 * z2, -w2 * x2, ux * wz * 2,
1006# -w2 * y2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2)
1007# if r > 0:
1008# r = sqrt(r) * bc # == a * a * b * c / a2
1009# elif r < 0: # LOS pointing away from or missing the earth
1010# raise _Error(_opposite_ if max(ux, vy, wz) > 0 else _outside_)
1011#
1012# d = Fdot(t, ux, vy, wz).fadd_(r).fover(m) # -r for antipode, a2 factored out
1013# if d > 0: # POV inside or LOS missing, outside the earth
1014# s = fsumf_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0) # like _sideOf
1015# raise _Error(_outside_ if s > 0 else _inside_)
1016# elif fsumf_(x2, y2, z2) < d**2: # d past earth center
1017# raise _Error(_too_(_distant_))
1018#
1019# r = p3.minus(u3.times(d))
1020# # h = p3.minus(r).length # distance to ellipsoid
1022 r = _xnamed(r, name or hartzell.__name__)
1023 if LatLon_and_kwds:
1024 c = _MODS.cartesianBase.CartesianBase(r, datum=D, name=r.name)
1025 r = c.toLatLon(**LatLon_and_kwds)
1026 return r
1029def haversine(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
1030 '''Compute the distance between two (spherical) points using the
1031 U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}
1032 formula.
1034 @arg lat1: Start latitude (C{degrees}).
1035 @arg lon1: Start longitude (C{degrees}).
1036 @arg lat2: End latitude (C{degrees}).
1037 @arg lon2: End longitude (C{degrees}).
1038 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
1039 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
1040 L{a_f2Tuple}) to use.
1041 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1042 B{C{lat2}} and B{C{lon2}} (C{bool}).
1044 @return: Distance (C{meter}, same units as B{C{radius}}).
1046 @raise TypeError: Invalid B{C{radius}}.
1048 @see: U{Distance between two (spherical) points
1049 <https://www.EdWilliams.org/avform.htm#Dist>}, functions
1050 L{cosineLaw}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},
1051 L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny}, L{flatPolar},
1052 L{thomas} and L{vincentys} and methods L{Ellipsoid.distance2},
1053 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.
1055 @note: See note at function L{vincentys_}.
1056 '''
1057 return _dS(haversine_, radius, wrap, lat1, lon1, lat2, lon2)
1060def haversine_(phi2, phi1, lam21):
1061 '''Compute the I{angular} distance between two (spherical) points
1062 using the U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}
1063 formula.
1065 @arg phi2: End latitude (C{radians}).
1066 @arg phi1: Start latitude (C{radians}).
1067 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
1069 @return: Angular distance (C{radians}).
1071 @see: Functions L{haversine}, L{cosineAndoyerLambert_},
1072 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
1073 L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
1074 L{flatPolar_}, L{thomas_} and L{vincentys_}.
1076 @note: See note at function L{vincentys_}.
1077 '''
1078 def _hsin(rad):
1079 return sin(rad * _0_5)**2
1081 h = _hsin(phi2 - phi1) + cos(phi1) * cos(phi2) * _hsin(lam21) # haversine
1082 return atan2(sqrt0(h), sqrt0(_1_0 - h)) * _2_0 # == asin(sqrt(h)) * 2
1085def heightOf(angle, distance, radius=R_M):
1086 '''Determine the height above the (spherical) earth' surface after
1087 traveling along a straight line at a given tilt.
1089 @arg angle: Tilt angle above horizontal (C{degrees}).
1090 @arg distance: Distance along the line (C{meter} or same units as
1091 B{C{radius}}).
1092 @kwarg radius: Optional mean earth radius (C{meter}).
1094 @return: Height (C{meter}, same units as B{C{distance}} and B{C{radius}}).
1096 @raise ValueError: Invalid B{C{angle}}, B{C{distance}} or B{C{radius}}.
1098 @see: U{MultiDop geog_lib.GeogBeamHt<https://GitHub.com/NASA/MultiDop>}
1099 (U{Shapiro et al. 2009, JTECH
1100 <https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}
1101 and U{Potvin et al. 2012, JTECH
1102 <https://Journals.AMetSoc.org/doi/abs/10.1175/JTECH-D-11-00019.1>}).
1103 '''
1104 r = h = Radius(radius)
1105 d = fabs(Distance(distance))
1106 if d > h:
1107 d, h = h, d
1109 if d > EPS0: # and h > EPS0
1110 d = d / h # /= h chokes PyChecker
1111 s = sin(Phi_(angle=angle, clip=_180_0))
1112 s = fsumf_(_1_0, _2_0 * s * d, d**2)
1113 if s > 0:
1114 return h * sqrt(s) - r
1116 raise _ValueError(angle=angle, distance=distance, radius=radius)
1119def heightOrthometric(h_ll, N):
1120 '''Get the I{orthometric} height B{H}, the height above the geoid, earth surface.
1122 @arg h_ll: The height above the ellipsoid (C{meter}) or an I{ellipsoidal}
1123 location (C{LatLon} with a C{height} or C{h} attribute).
1124 @arg N: The I{geoid} height (C{meter}), the height of the geoid above the
1125 ellipsoid at the same B{C{h_ll}} location.
1127 @return: I{Orthometric} height C{B{H} = B{h} - B{N}} (C{meter}, same units
1128 as B{C{h}} and B{C{N}}).
1130 @see: U{Ellipsoid, Geoid, and Othometric Heights<https://www.NGS.NOAA.gov/
1131 GEOID/PRESENTATIONS/2007_02_24_CCPS/Roman_A_PLSC2007notes.pdf>}, page
1132 6 and module L{pygeodesy.geoids}.
1133 '''
1134 h = h_ll if isscalar(h_ll) else _xattr(h_ll, height=_xattr(h_ll, h=0))
1135 return Height(H=Height(h=h) - Height(N=N))
1138def horizon(height, radius=R_M, refraction=False):
1139 '''Determine the distance to the horizon from a given altitude
1140 above the (spherical) earth.
1142 @arg height: Altitude (C{meter} or same units as B{C{radius}}).
1143 @kwarg radius: Optional mean earth radius (C{meter}).
1144 @kwarg refraction: Consider atmospheric refraction (C{bool}).
1146 @return: Distance (C{meter}, same units as B{C{height}} and B{C{radius}}).
1148 @raise ValueError: Invalid B{C{height}} or B{C{radius}}.
1150 @see: U{Distance to horizon<https://www.EdWilliams.org/avform.htm#Horizon>}.
1151 '''
1152 h, r = Height(height), Radius(radius)
1153 if min(h, r) < 0:
1154 raise _ValueError(height=height, radius=radius)
1156 if refraction:
1157 d2 = 2.415750694528 * h * r # 2.0 / 0.8279
1158 else:
1159 d2 = h * fsumf_(r, r, h)
1160 return sqrt0(d2)
1163class _idllmn6(object): # see also .geodesicw._wargs, .vector2d._numpy
1164 '''(INTERNAL) Helper for C{intersection2} and C{intersections2}.
1165 '''
1166 @contextmanager # <https://www.python.org/dev/peps/pep-0343/> Examples
1167 def __call__(self, datum, lat1, lon1, lat2, lon2, small, wrap, s, **kwds):
1168 try:
1169 if wrap:
1170 _, lat2, lon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)
1171 kwds = _xkwds(kwds, wrap=wrap) # for _xError
1172 m = small if small is _100km else Meter_(small=small)
1173 n = (intersections2 if s else intersection2).__name__
1174 if datum is None or euclidean(lat1, lon1, lat2, lon2) < m:
1175 d, m = None, _MODS.vector3d
1176 _i = m._intersects2 if s else m._intersect3d3
1177 elif isscalar(datum) and datum < 0 and not s:
1178 d = _spherical_datum(-datum, name=n)
1179 m = _MODS.sphericalNvector
1180 _i = m.intersection
1181 else:
1182 d = _spherical_datum(datum, name=n)
1183 if d.isSpherical:
1184 m = _MODS.sphericalTrigonometry
1185 _i = m._intersects2 if s else m._intersect
1186 elif d.isEllipsoidal:
1187 try:
1188 if d.ellipsoid.geodesic:
1189 pass
1190 m = _MODS.ellipsoidalKarney
1191 except ImportError:
1192 m = _MODS.ellipsoidalExact
1193 _i = m._intersections2 if s else m._intersection3 # ellispoidalBaseDI
1194 else:
1195 raise _TypeError(datum=datum)
1196 yield _i, d, lat2, lon2, m, n
1198 except (TypeError, ValueError) as x:
1199 raise _xError(x, lat1=lat1, lon1=lon1, datum=datum,
1200 lat2=lat2, lon2=lon2, small=small, **kwds)
1202_idllmn6 = _idllmn6() # PYCHOK singleton
1205def intersection2(lat1, lon1, bearing1,
1206 lat2, lon2, bearing2, datum=None, wrap=False, small=_100km): # was=True
1207 '''I{Conveniently} compute the intersection of two lines each defined
1208 by a (geodetic) point and a bearing from North, using either ...
1210 1) L{vector3d.intersection3d3} for B{C{small}} distances (below 100 Km
1211 or about 0.88 degrees) or if I{no} B{C{datum}} is specified, or ...
1213 2) L{sphericalTrigonometry.intersection} for a spherical B{C{datum}}
1214 or a C{scalar B{datum}} representing the earth radius, conventionally
1215 in C{meter} or ...
1217 3) L{sphericalNvector.intersection} if B{C{datum}} is a I{negative}
1218 C{scalar}, (negative) earth radius, conventionally in C{meter} or ...
1220 4) L{ellipsoidalKarney.intersection3} for an ellipsoidal B{C{datum}}
1221 and if I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
1222 is installed, otherwise ...
1224 5) L{ellipsoidalExact.intersection3}, provided B{C{datum}} is ellipsoidal.
1226 @arg lat1: Latitude of the first point (C{degrees}).
1227 @arg lon1: Longitude of the first point (C{degrees}).
1228 @arg bearing1: Bearing at the first point (compass C{degrees}).
1229 @arg lat2: Latitude of the second point (C{degrees}).
1230 @arg lon2: Longitude of the second point (C{degrees}).
1231 @arg bearing2: Bearing at the second point (compass C{degrees}).
1232 @kwarg datum: Optional datum (L{Datum}) or ellipsoid (L{Ellipsoid},
1233 L{Ellipsoid2} or L{a_f2Tuple}) or C{scalar} earth
1234 radius (C{meter}, same units as B{C{radius1}} and
1235 B{C{radius2}}) or C{None}.
1236 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
1237 and B{C{lon2}} (C{bool}).
1238 @kwarg small: Upper limit for small distances (C{meter}).
1240 @return: A L{LatLon2Tuple}C{(lat, lon)} with the lat- and
1241 longitude of the intersection point.
1243 @raise IntersectionError: Ambiguous or infinite intersection
1244 or colinear, parallel or otherwise
1245 non-intersecting lines.
1247 @raise TypeError: Invalid B{C{datum}}.
1249 @raise UnitError: Invalid B{C{lat1}}, B{C{lon1}}, B{C{bearing1}},
1250 B{C{lat2}}, B{C{lon2}} or B{C{bearing2}}.
1252 @see: Method L{RhumbLine.intersection2}.
1254 @note: The returned intersections may be near-antipodal.
1255 '''
1256 b1 = Bearing(bearing1=bearing1)
1257 b2 = Bearing(bearing2=bearing2)
1258 with _idllmn6(datum, lat1, lon1, lat2, lon2,
1259 small, wrap, False, bearing1=b1, bearing2=b2) as t:
1260 _i, d, lat2, lon2, m, n = t
1261 if d is None:
1262 t, _, _ = _i(m.Vector3d(lon1, lat1, 0), b1,
1263 m.Vector3d(lon2, lat2, 0), b2, useZ=False)
1264 t = LatLon2Tuple(t.y, t.x, name=n)
1266 else:
1267 t = _i(m.LatLon(lat1, lon1, datum=d), b1,
1268 m.LatLon(lat2, lon2, datum=d), b2, height=0, wrap=False)
1269 if isinstance(t, Intersection3Tuple): # ellipsoidal
1270 t, _, _ = t
1271 t = LatLon2Tuple(t.lat, t.lon, name=n)
1272 return t
1275def intersections2(lat1, lon1, radius1,
1276 lat2, lon2, radius2, datum=None, wrap=False, small=_100km): # was=True
1277 '''I{Conveniently} compute the intersections of two circles each defined
1278 by a (geodetic) center point and a radius, using either ...
1280 1) L{vector3d.intersections2} for B{C{small}} distances (below 100 Km
1281 or about 0.88 degrees) or if I{no} B{C{datum}} is specified, or ...
1283 2) L{sphericalTrigonometry.intersections2} for a spherical B{C{datum}}
1284 or a C{scalar B{datum}} representing the earth radius, conventionally
1285 in C{meter} or ...
1287 3) L{ellipsoidalKarney.intersections2} for an ellipsoidal B{C{datum}}
1288 and if I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
1289 is installed, otherwise ...
1291 4) L{ellipsoidalExact.intersections2}, provided B{C{datum}} is ellipsoidal.
1293 @arg lat1: Latitude of the first circle center (C{degrees}).
1294 @arg lon1: Longitude of the first circle center (C{degrees}).
1295 @arg radius1: Radius of the first circle (C{meter}, conventionally).
1296 @arg lat2: Latitude of the second circle center (C{degrees}).
1297 @arg lon2: Longitude of the second circle center (C{degrees}).
1298 @arg radius2: Radius of the second circle (C{meter}, same units as B{C{radius1}}).
1299 @kwarg datum: Optional datum (L{Datum}) or ellipsoid (L{Ellipsoid},
1300 L{Ellipsoid2} or L{a_f2Tuple}) or C{scalar} earth
1301 radius (C{meter}, same units as B{C{radius1}} and
1302 B{C{radius2}}) or C{None}.
1303 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}
1304 and B{C{lon2}} (C{bool}).
1305 @kwarg small: Upper limit for small distances (C{meter}).
1307 @return: 2-Tuple of the intersection points, each a
1308 L{LatLon2Tuple}C{(lat, lon)}. For abutting circles, the
1309 points are the same instance, aka the I{radical center}.
1311 @raise IntersectionError: Concentric, antipodal, invalid or
1312 non-intersecting circles or no
1313 convergence.
1315 @raise TypeError: Invalid B{C{datum}}.
1317 @raise UnitError: Invalid B{C{lat1}}, B{C{lon1}}, B{C{radius1}},
1318 B{C{lat2}}, B{C{lon2}} or B{C{radius2}}.
1319 '''
1320 r1 = Radius_(radius1=radius1)
1321 r2 = Radius_(radius2=radius2)
1322 with _idllmn6(datum, lat1, lon1, lat2, lon2,
1323 small, wrap, True, radius1=r1, radius2=r2) as t:
1324 _i, d, lat2, lon2, m, n = t
1325 if d is None:
1326 r1 = m2degrees(r1, radius=R_M, lat=lat1)
1327 r2 = m2degrees(r2, radius=R_M, lat=lat2)
1329 def _V2T(x, y, _, **unused): # _ == z unused
1330 return LatLon2Tuple(y, x, name=n)
1332 t = _i(m.Vector3d(lon1, lat1, 0), r1,
1333 m.Vector3d(lon2, lat2, 0), r2, sphere=False,
1334 Vector=_V2T)
1335 else:
1336 def _LL2T(lat, lon, **unused):
1337 return LatLon2Tuple(lat, lon, name=n)
1339 t = _i(m.LatLon(lat1, lon1, datum=d), r1,
1340 m.LatLon(lat2, lon2, datum=d), r2,
1341 LatLon=_LL2T, height=0, wrap=False)
1342 return t
1345def isantipode(lat1, lon1, lat2, lon2, eps=EPS):
1346 '''Check whether two points are I{antipodal}, on diametrically
1347 opposite sides of the earth.
1349 @arg lat1: Latitude of one point (C{degrees}).
1350 @arg lon1: Longitude of one point (C{degrees}).
1351 @arg lat2: Latitude of the other point (C{degrees}).
1352 @arg lon2: Longitude of the other point (C{degrees}).
1353 @kwarg eps: Tolerance for near-equality (C{degrees}).
1355 @return: C{True} if points are antipodal within the
1356 B{C{eps}} tolerance, C{False} otherwise.
1358 @see: Functions L{isantipode_} and L{antipode}.
1359 '''
1360 return (fabs(lat1 + lat2) <= eps and
1361 fabs(lon1 + lon2) <= eps) or _isequalTo(
1362 normal(lat1, lon1), antipode(lat2, lon2), eps)
1365def isantipode_(phi1, lam1, phi2, lam2, eps=EPS):
1366 '''Check whether two points are I{antipodal}, on diametrically
1367 opposite sides of the earth.
1369 @arg phi1: Latitude of one point (C{radians}).
1370 @arg lam1: Longitude of one point (C{radians}).
1371 @arg phi2: Latitude of the other point (C{radians}).
1372 @arg lam2: Longitude of the other point (C{radians}).
1373 @kwarg eps: Tolerance for near-equality (C{radians}).
1375 @return: C{True} if points are antipodal within the
1376 B{C{eps}} tolerance, C{False} otherwise.
1378 @see: Functions L{isantipode} and L{antipode_}.
1379 '''
1380 return (fabs(phi1 + phi2) <= eps and
1381 fabs(lam1 + lam2) <= eps) or _isequalTo_(
1382 normal_(phi1, lam1), antipode_(phi2, lam2), eps)
1385def _isequalTo(p1, p2, eps=EPS):
1386 '''Compare 2 point lat-/lons ignoring C{class}.
1387 '''
1388 return (fabs(p1.lat - p2.lat) <= eps and
1389 fabs(p1.lon - p2.lon) <= eps) if eps else (p1.latlon == p2.latlon)
1392def _isequalTo_(p1, p2, eps=EPS):
1393 '''(INTERNAL) Compare 2 point phi-/lams ignoring C{class}.
1394 '''
1395 return (fabs(p1.phi - p2.phi) <= eps and
1396 fabs(p1.lam - p2.lam) <= eps) if eps else (p1.philam == p2.philam)
1399def isnormal(lat, lon, eps=0):
1400 '''Check whether B{C{lat}} I{and} B{C{lon}} are within their
1401 respective I{normal} range in C{degrees}.
1403 @arg lat: Latitude (C{degrees}).
1404 @arg lon: Longitude (C{degrees}).
1405 @kwarg eps: Optional tolerance C{degrees}).
1407 @return: C{True} if C{(abs(B{lat}) + B{eps}) <= 90} and
1408 C{(abs(B{lon}) + B{eps}) <= 180}, C{False} othwerwise.
1410 @see: Functions L{isnormal_} and L{normal}.
1411 '''
1412 return (_90_0 - fabs(lat)) >= eps and (_180_0 - fabs(lon)) >= eps
1415def isnormal_(phi, lam, eps=0):
1416 '''Check whether B{C{phi}} I{and} B{C{lam}} are within their
1417 respective I{normal} range in C{radians}.
1419 @arg phi: Latitude (C{radians}).
1420 @arg lam: Longitude (C{radians}).
1421 @kwarg eps: Optional tolerance C{radians}).
1423 @return: C{True} if C{(abs(B{phi}) + B{eps}) <= PI/2} and
1424 C{(abs(B{lam}) + B{eps}) <= PI}, C{False} othwerwise.
1426 @see: Functions L{isnormal} and L{normal_}.
1427 '''
1428 return (PI_2 - fabs(phi)) >= eps and (PI - fabs(lam)) >= eps
1431def latlon2n_xyz(lat, lon, name=NN):
1432 '''Convert lat-, longitude to C{n-vector} (I{normal} to the
1433 earth's surface) X, Y and Z components.
1435 @arg lat: Latitude (C{degrees}).
1436 @arg lon: Longitude (C{degrees}).
1437 @kwarg name: Optional name (C{str}).
1439 @return: A L{Vector3Tuple}C{(x, y, z)}.
1441 @see: Function L{philam2n_xyz}.
1443 @note: These are C{n-vector} x, y and z components,
1444 I{NOT} geocentric ECEF x, y and z coordinates!
1445 '''
1446 return _2n_xyz(name, *sincos2d_(lat, lon))
1449def _normal2(a, b, n_2, n, n2):
1450 '''(INTERNAL) Helper for C{normal} and C{normal_}.
1451 '''
1452 if fabs(b) > n:
1453 b = remainder(b, n2)
1454 if fabs(a) > n_2:
1455 r = remainder(a, n)
1456 if r != a:
1457 a = -r
1458 b -= n if b > 0 else -n
1459 return float0_(a, b)
1462def normal(lat, lon, name=NN):
1463 '''Normalize a lat- I{and} longitude pair in C{degrees}.
1465 @arg lat: Latitude (C{degrees}).
1466 @arg lon: Longitude (C{degrees}).
1467 @kwarg name: Optional name (C{str}).
1469 @return: L{LatLon2Tuple}C{(lat, lon)} with C{abs(lat) <= 90}
1470 and C{abs(lon) <= 180}.
1472 @see: Functions L{normal_} and L{isnormal}.
1473 '''
1474 return LatLon2Tuple(*_normal2(lat, lon, _90_0, _180_0, _360_0),
1475 name=name or normal.__name__)
1478def normal_(phi, lam, name=NN):
1479 '''Normalize a lat- I{and} longitude pair in C{radians}.
1481 @arg phi: Latitude (C{radians}).
1482 @arg lam: Longitude (C{radians}).
1483 @kwarg name: Optional name (C{str}).
1485 @return: L{PhiLam2Tuple}C{(phi, lam)} with C{abs(phi) <= PI/2}
1486 and C{abs(lam) <= PI}.
1488 @see: Functions L{normal} and L{isnormal_}.
1489 '''
1490 return PhiLam2Tuple(*_normal2(phi, lam, PI_2, PI, PI2),
1491 name=name or normal_.__name__)
1494def _2n_xyz(name, sa, ca, sb, cb):
1495 '''(INTERNAL) Helper for C{latlon2n_xyz} and C{philam2n_xyz}.
1496 '''
1497 # Kenneth Gade eqn 3, but using right-handed
1498 # vector x -> 0°E,0°N, y -> 90°E,0°N, z -> 90°N
1499 return Vector3Tuple(ca * cb, ca * sb, sa, name=name)
1502def n_xyz2latlon(x, y, z, name=NN):
1503 '''Convert C{n-vector} components to lat- and longitude in C{degrees}.
1505 @arg x: X component (C{scalar}).
1506 @arg y: Y component (C{scalar}).
1507 @arg z: Z component (C{scalar}).
1508 @kwarg name: Optional name (C{str}).
1510 @return: A L{LatLon2Tuple}C{(lat, lon)}.
1512 @see: Function L{n_xyz2philam}.
1513 '''
1514 return LatLon2Tuple(atan2d(z, hypot(x, y)), atan2d(y, x), name=name)
1517def n_xyz2philam(x, y, z, name=NN):
1518 '''Convert C{n-vector} components to lat- and longitude in C{radians}.
1520 @arg x: X component (C{scalar}).
1521 @arg y: Y component (C{scalar}).
1522 @arg z: Z component (C{scalar}).
1523 @kwarg name: Optional name (C{str}).
1525 @return: A L{PhiLam2Tuple}C{(phi, lam)}.
1527 @see: Function L{n_xyz2latlon}.
1528 '''
1529 return PhiLam2Tuple(atan2(z, hypot(x, y)), atan2(y, x), name=name)
1532def _opposes(d, m, n, n2):
1533 '''(INTERNAL) Helper for C{opposing} and C{opposing_}.
1534 '''
1535 d = d % n2 # -20 % 360 == 340, -1 % PI2 == PI2 - 1
1536 return False if d < m or d > (n2 - m) else (
1537 True if (n - m) < d < (n + m) else None)
1540def opposing(bearing1, bearing2, margin=_90_0):
1541 '''Compare the direction of two bearings given in C{degrees}.
1543 @arg bearing1: First bearing (compass C{degrees}).
1544 @arg bearing2: Second bearing (compass C{degrees}).
1545 @kwarg margin: Optional, interior angle bracket (C{degrees}).
1547 @return: C{True} if both bearings point in opposite, C{False} if
1548 in similar or C{None} if in perpendicular directions.
1550 @see: Function L{opposing_}.
1551 '''
1552 m = Degrees_(margin=margin, low=EPS0, high=_90_0)
1553 return _opposes(bearing2 - bearing1, m, _180_0, _360_0)
1556def opposing_(radians1, radians2, margin=PI_2):
1557 '''Compare the direction of two bearings given in C{radians}.
1559 @arg radians1: First bearing (C{radians}).
1560 @arg radians2: Second bearing (C{radians}).
1561 @kwarg margin: Optional, interior angle bracket (C{radians}).
1563 @return: C{True} if both bearings point in opposite, C{False} if
1564 in similar or C{None} if in perpendicular directions.
1566 @see: Function L{opposing}.
1567 '''
1568 m = Radians_(margin=margin, low=EPS0, high=PI_2)
1569 return _opposes(radians2 - radians1, m, PI, PI2)
1572def philam2n_xyz(phi, lam, name=NN):
1573 '''Convert lat-, longitude to C{n-vector} (I{normal} to the
1574 earth's surface) X, Y and Z components.
1576 @arg phi: Latitude (C{radians}).
1577 @arg lam: Longitude (C{radians}).
1578 @kwarg name: Optional name (C{str}).
1580 @return: A L{Vector3Tuple}C{(x, y, z)}.
1582 @see: Function L{latlon2n_xyz}.
1584 @note: These are C{n-vector} x, y and z components,
1585 I{NOT} geocentric ECEF x, y and z coordinates!
1586 '''
1587 return _2n_xyz(name, *sincos2_(phi, lam))
1590def _radical2(d, r1, r2): # in .ellipsoidalBaseDI, .sphericalTrigonometry, .vector3d
1591 # (INTERNAL) See C{radical2} below
1592 # assert d > EPS0
1593 r = fsumf_(_1_0, (r1 / d)**2, -(r2 / d)**2) * _0_5
1594 return Radical2Tuple(max(_0_0, min(_1_0, r)), r * d)
1597def radical2(distance, radius1, radius2):
1598 '''Compute the I{radical ratio} and I{radical line} of two
1599 U{intersecting circles<https://MathWorld.Wolfram.com/
1600 Circle-CircleIntersection.html>}.
1602 The I{radical line} is perpendicular to the axis thru the
1603 centers of the circles at C{(0, 0)} and C{(B{distance}, 0)}.
1605 @arg distance: Distance between the circle centers (C{scalar}).
1606 @arg radius1: Radius of the first circle (C{scalar}).
1607 @arg radius2: Radius of the second circle (C{scalar}).
1609 @return: A L{Radical2Tuple}C{(ratio, xline)} where C{0.0 <=
1610 ratio <= 1.0} and C{xline} is along the B{C{distance}}.
1612 @raise IntersectionError: The B{C{distance}} exceeds the sum
1613 of B{C{radius1}} and B{C{radius2}}.
1615 @raise UnitError: Invalid B{C{distance}}, B{C{radius1}} or
1616 B{C{radius2}}.
1618 @see: U{Circle-Circle Intersection
1619 <https://MathWorld.Wolfram.com/Circle-CircleIntersection.html>}.
1620 '''
1621 d = Distance_(distance, low=_0_0)
1622 r1 = Radius_(radius1=radius1)
1623 r2 = Radius_(radius2=radius2)
1624 if d > (r1 + r2):
1625 raise IntersectionError(distance=d, radius1=r1, radius2=r2,
1626 txt=_too_(_distant_))
1627 return _radical2(d, r1, r2) if d > EPS0 else \
1628 Radical2Tuple(_0_5, _0_0)
1631class Radical2Tuple(_NamedTuple):
1632 '''2-Tuple C{(ratio, xline)} of the I{radical} C{ratio} and
1633 I{radical} C{xline}, both C{scalar} and C{0.0 <= ratio <= 1.0}
1634 '''
1635 _Names_ = (_ratio_, _xline_)
1636 _Units_ = ( Scalar, Scalar)
1639def _radistance(inst):
1640 '''(INTERNAL) Helper for the L{frechet._FrecherMeterRadians}
1641 and L{hausdorff._HausdorffMeterRedians} classes.
1642 '''
1643 kwds_ = _xkwds(inst._kwds, wrap=False)
1644 wrap_ = _xkwds_pop(kwds_, wrap=False)
1645 func_ = inst._func_
1646 try: # calling lower-overhead C{func_}
1647 func_(0, _0_25, _0_5, **kwds_)
1648 wrap_ = _Wrap._philamop(wrap_)
1649 except TypeError:
1650 return inst.distance
1652 def _philam(p):
1653 try:
1654 return p.phi, p.lam
1655 except AttributeError: # no .phi or .lam
1656 return radians(p.lat), radians(p.lon)
1658 def _func_wrap(point1, point2):
1659 phi1, lam1 = wrap_(*_philam(point1))
1660 phi2, lam2 = wrap_(*_philam(point2))
1661 return func_(phi2, phi1, lam2 - lam1, **kwds_)
1663 inst._units = inst._units_
1664 return _func_wrap
1667def _scale_deg(lat1, lat2): # degrees
1668 # scale factor cos(mean of lats) for delta lon
1669 m = fabs(lat1 + lat2) * _0_5
1670 return cos(radians(m)) if m < 90 else _0_0
1673def _scale_rad(phi1, phi2): # radians, by .frechet, .hausdorff, .heights
1674 # scale factor cos(mean of phis) for delta lam
1675 m = fabs(phi1 + phi2) * _0_5
1676 return cos(m) if m < PI_2 else _0_0
1679def _sincosa6(phi2, phi1, lam21): # [4] in cosineLaw
1680 '''(INTERNAL) C{sin}es, C{cos}ines and C{acos}ine.
1681 '''
1682 s2, c2, s1, c1, _, c21 = sincos2_(phi2, phi1, lam21)
1683 return s2, c2, s1, c1, acos1(s1 * s2 + c1 * c2 * c21), c21
1686def thomas(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):
1687 '''Compute the distance between two (ellipsoidal) points using
1688 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}
1689 formula.
1691 @arg lat1: Start latitude (C{degrees}).
1692 @arg lon1: Start longitude (C{degrees}).
1693 @arg lat2: End latitude (C{degrees}).
1694 @arg lon2: End longitude (C{degrees}).
1695 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},
1696 L{Ellipsoid2} or L{a_f2Tuple}) to use.
1697 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1698 B{C{lat2}} and B{C{lon2}} (C{bool}).
1700 @return: Distance (C{meter}, same units as the B{C{datum}}'s
1701 ellipsoid axes).
1703 @raise TypeError: Invalid B{C{datum}}.
1705 @see: Functions L{thomas_}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},
1706 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},
1707 L{flatPolar}, L{haversine}, L{vincentys} and method L{Ellipsoid.distance2}.
1708 '''
1709 return _dE(thomas_, datum, wrap, lat1, lon1, lat2, lon2)
1712def thomas_(phi2, phi1, lam21, datum=_WGS84):
1713 '''Compute the I{angular} distance between two (ellipsoidal) points using
1714 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}
1715 formula.
1717 @arg phi2: End latitude (C{radians}).
1718 @arg phi1: Start latitude (C{radians}).
1719 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
1720 @kwarg datum: Datum or ellipsoid to use (L{Datum}, L{Ellipsoid},
1721 L{Ellipsoid2} or L{a_f2Tuple}).
1723 @return: Angular distance (C{radians}).
1725 @raise TypeError: Invalid B{C{datum}}.
1727 @see: Functions L{thomas}, L{cosineAndoyerLambert_},
1728 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
1729 L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
1730 L{flatPolar_}, L{haversine_} and L{vincentys_} and U{Geodesy-PHP
1731 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/
1732 Distance/ThomasFormula.php>}.
1733 '''
1734 s2, c2, s1, c1, r, _ = _sincosa6(phi2, phi1, lam21)
1735 if r and isnon0(c1) and isnon0(c2):
1736 E = _ellipsoidal(datum, thomas_)
1737 if E.f:
1738 r1 = atan2(E.b_a * s1, c1)
1739 r2 = atan2(E.b_a * s2, c2)
1741 j = (r2 + r1) * _0_5
1742 k = (r2 - r1) * _0_5
1743 sj, cj, sk, ck, h, _ = sincos2_(j, k, lam21 * _0_5)
1745 h = fsumf_(sk**2, (ck * h)**2, -(sj * h)**2)
1746 u = _1_0 - h
1747 if isnon0(u) and isnon0(h):
1748 r = atan(sqrt0(h / u)) * 2 # == acos(1 - 2 * h)
1749 sr, cr = sincos2(r)
1750 if isnon0(sr):
1751 u = 2 * (sj * ck)**2 / u
1752 h = 2 * (sk * cj)**2 / h
1753 x = u + h
1754 y = u - h
1756 s = r / sr
1757 e = 4 * s**2
1758 d = 2 * cr
1759 a = e * d
1760 b = 2 * r
1761 c = s - (a - d) * _0_5
1762 f = E.f * _0_25
1764 t = fsumf_(a * x, -b * y, c * x**2, -d * y**2, e * x * y)
1765 r -= fsumf_(s * x, -y, -t * f * _0_25) * f * sr
1766 return r
1769def vincentys(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):
1770 '''Compute the distance between two (spherical) points using
1771 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}
1772 spherical formula.
1774 @arg lat1: Start latitude (C{degrees}).
1775 @arg lon1: Start longitude (C{degrees}).
1776 @arg lat2: End latitude (C{degrees}).
1777 @arg lon2: End longitude (C{degrees}).
1778 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})
1779 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or
1780 L{a_f2Tuple}) to use.
1781 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1782 B{C{lat2}} and B{C{lon2}} (C{bool}).
1784 @return: Distance (C{meter}, same units as B{C{radius}}).
1786 @raise UnitError: Invalid B{C{radius}}.
1788 @see: Functions L{vincentys_}, L{cosineAndoyerLambert},
1789 L{cosineForsytheAndoyerLambert},L{cosineLaw}, L{equirectangular},
1790 L{euclidean}, L{flatLocal}/L{hubeny}, L{flatPolar},
1791 L{haversine} and L{thomas} and methods L{Ellipsoid.distance2},
1792 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.
1794 @note: See note at function L{vincentys_}.
1795 '''
1796 return _dS(vincentys_, radius, wrap, lat1, lon1, lat2, lon2)
1799def vincentys_(phi2, phi1, lam21):
1800 '''Compute the I{angular} distance between two (spherical) points using
1801 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}
1802 spherical formula.
1804 @arg phi2: End latitude (C{radians}).
1805 @arg phi1: Start latitude (C{radians}).
1806 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).
1808 @return: Angular distance (C{radians}).
1810 @see: Functions L{vincentys}, L{cosineAndoyerLambert_},
1811 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
1812 L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
1813 L{flatPolar_}, L{haversine_} and L{thomas_}.
1815 @note: Functions L{vincentys_}, L{haversine_} and L{cosineLaw_}
1816 produce equivalent results, but L{vincentys_} is suitable
1817 for antipodal points and slightly more expensive (M{3 cos,
1818 3 sin, 1 hypot, 1 atan2, 6 mul, 2 add}) than L{haversine_}
1819 (M{2 cos, 2 sin, 2 sqrt, 1 atan2, 5 mul, 1 add}) and
1820 L{cosineLaw_} (M{3 cos, 3 sin, 1 acos, 3 mul, 1 add}).
1821 '''
1822 s1, c1, s2, c2, s21, c21 = sincos2_(phi1, phi2, lam21)
1824 c = c2 * c21
1825 x = s1 * s2 + c1 * c
1826 y = c1 * s2 - s1 * c
1827 return atan2(hypot(c2 * s21, y), x)
1829# **) MIT License
1830#
1831# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1832#
1833# Permission is hereby granted, free of charge, to any person obtaining a
1834# copy of this software and associated documentation files (the "Software"),
1835# to deal in the Software without restriction, including without limitation
1836# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1837# and/or sell copies of the Software, and to permit persons to whom the
1838# Software is furnished to do so, subject to the following conditions:
1839#
1840# The above copyright notice and this permission notice shall be included
1841# in all copies or substantial portions of the Software.
1842#
1843# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1844# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1845# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1846# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1847# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1848# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1849# OTHER DEALINGS IN THE SOFTWARE.