Coverage for pygeodesy/latlonBase.py: 96%
409 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private base class L{LatLonBase} for elliposiodal, spherical
5and N-vectorial C{LatLon}s.
7After I{(C) Chris Veness 2011-2015} published under the same MIT Licence**,
8see U{https://www.Movable-Type.co.UK/scripts/latlong.html},
9U{<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}
10and U{https://www.Movable-Type.co.UK/scripts/latlong-vectors.html}.
11'''
13from pygeodesy.basics import isscalar, isstr, _xinstanceof
14from pygeodesy.constants import EPS, EPS0, EPS1, EPS4, INT0, R_M, \
15 _0_0, _0_5, _1_0
16# from pygeodesy.datums import _spherical_datum # from .formy
17from pygeodesy.dms import F_D, F_DMS, latDMS, lonDMS, parse3llh
18from pygeodesy.errors import _incompatible, IntersectionError, _TypeError, \
19 _ValueError, _xdatum, _xError, _xkwds, _xkwds_not
20from pygeodesy.formy import antipode, compassAngle, cosineAndoyerLambert_, \
21 cosineForsytheAndoyerLambert_, cosineLaw, \
22 equirectangular, euclidean, flatLocal_, \
23 flatPolar, hartzell, haversine, isantipode, \
24 isnormal, normal, philam2n_xyz, \
25 _spherical_datum, thomas_, vincentys
26from pygeodesy.interns import NN, _COMMASPACE_, _concentric_, _height_, \
27 _intersection_, _m_, _no_, _overlap_, _point_
28from pygeodesy.iters import PointsIter, points2
29from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
30from pygeodesy.named import _NamedBase, notOverloaded
31from pygeodesy.namedTuples import Bounds2Tuple, LatLon2Tuple, PhiLam2Tuple, \
32 Trilaterate5Tuple, Vector3Tuple
33from pygeodesy.props import deprecated_method, Property, Property_RO, \
34 property_RO, _update_all
35from pygeodesy.streprs import Fmt, hstr
36from pygeodesy.units import Distance_, Lat, Lon, Height, Radius, Radius_, \
37 Scalar, Scalar_
38from pygeodesy.utily import _unrollon, unrollPI
39from pygeodesy.vector2d import _circin6, Circin6Tuple, _circum3, Circum3Tuple, \
40 circum4_, Circum4Tuple, _radii11ABC
41from pygeodesy.vector3d import nearestOn6, Vector3d
43from math import asin, cos, degrees, fabs, radians
45__all__ = _ALL_LAZY.latlonBase
46__version__ = '23.03.30'
49class LatLonBase(_NamedBase):
50 '''(INTERNAL) Base class for C{LatLon} points on spherical or
51 ellipsoidal earth models.
52 '''
53 _clipid = INT0 # polygonal clip, see .booleans
54 _datum = None # L{Datum}, to be overriden
55 _height = INT0 # height (C{meter}), default
56 _lat = 0 # latitude (C{degrees})
57 _lon = 0 # longitude (C{degrees})
59 def __init__(self, lat, lon, height=0, name=NN):
60 '''New C{LatLon}.
62 @arg lat: Latitude (C{degrees} or DMS C{str} with N or S suffix).
63 @arg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix).
64 @kwarg height: Optional height (C{meter} above or below the earth surface).
65 @kwarg name: Optional name (C{str}).
67 @return: New instance (C{LatLon}).
69 @raise RangeError: Value of B{C{lat}} or B{C{lon}} outside the valid
70 range and C{rangerrors} set to C{True}.
72 @raise UnitError: Invalid B{C{lat}}, B{C{lon}} or B{C{height}}.
74 @example:
76 >>> p = LatLon(50.06632, -5.71475)
77 >>> q = LatLon('50°03′59″N', """005°42'53"W""")
78 '''
79 if name:
80 self.name = name
82 self._lat = Lat(lat) # parseDMS2(lat, lon)
83 self._lon = Lon(lon) # PYCHOK LatLon2Tuple
84 if height: # elevation
85 self._height = Height(height)
87 def __eq__(self, other):
88 return self.isequalTo(other)
90 def __ne__(self, other):
91 return not self.isequalTo(other)
93 def __str__(self):
94 return self.toStr(form=F_D, prec=6)
96 def antipode(self, height=None):
97 '''Return the antipode, the point diametrically opposite
98 to this point.
100 @kwarg height: Optional height of the antipode (C{meter}),
101 this point's height otherwise.
103 @return: The antipodal point (C{LatLon}).
104 '''
105 a, b = antipode(self.lat, self.lon) # PYCHOK LatLon2Tuple
106 h = self.height if height is None else Height(height)
107 return self.classof(a, b, height=h)
109 @deprecated_method
110 def bounds(self, wide, tall, radius=R_M): # PYCHOK no cover
111 '''DEPRECATED, use method C{boundsOf}.'''
112 return self.boundsOf(wide, tall, radius=radius)
114 def boundsOf(self, wide, tall, radius=R_M, height=None):
115 '''Return the SW and NE lat-/longitude of a great circle
116 bounding box centered at this location.
118 @arg wide: Longitudinal box width (C{meter}, same units as
119 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).
120 @arg tall: Latitudinal box size (C{meter}, same units as
121 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).
122 @kwarg radius: Mean earth radius (C{meter}).
123 @kwarg height: Height for C{latlonSW} and C{latlonNE} (C{meter}),
124 overriding the point's height.
126 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)}, the
127 lower-left and upper-right corner (C{LatLon}).
129 @see: U{https://www.Movable-Type.co.UK/scripts/latlong-db.html}
130 '''
131 x = Scalar_(wide=wide) * _0_5
132 y = Scalar_(tall=tall) * _0_5
133 if radius is not None:
134 r = Radius_(radius)
135 c = cos(self.phi)
136 x = degrees(asin(x / r) / c) if fabs(c) > EPS0 else _0_0 # XXX
137 y = degrees(y / r)
138 x, y = fabs(x), fabs(y)
140 h = self.height if height is None else Height(height)
141 sw = self.classof(self.lat - y, self.lon - x, height=h)
142 ne = self.classof(self.lat + y, self.lon + x, height=h)
143 return Bounds2Tuple(sw, ne, name=self.name)
145 def chordTo(self, other, height=None):
146 '''Compute the length of the chord through the earth between
147 this and an other point.
149 @arg other: The other point (C{LatLon}).
150 @kwarg height: Overriding height for both points (C{meter})
151 or C{None} for each point's height.
153 @return: The chord length (conventionally C{meter}).
155 @raise TypeError: The B{C{other}} point is not C{LatLon}.
156 '''
157 def _v3d(ll):
158 t = ll.toEcef(height=height) # .toVector(Vector=Vector3d)
159 return Vector3d(t.x, t.y, t.z)
161 self.others(other)
162 return _v3d(self).minus(_v3d(other)).length
164 def circin6(self, point2, point3, eps=EPS4):
165 '''Return the radius and center of the I{inscribed} aka I{In-}circle
166 of the (planar) triangle formed by this and two other points.
168 @arg point2: Second point (C{LatLon}).
169 @arg point3: Third point (C{LatLon}).
170 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.
172 @return: L{Circin6Tuple}C{(radius, center, deltas, cA, cB, cC)}. The
173 C{center} and contact points C{cA}, C{cB} and C{cC}, each an
174 instance of this (sub-)class, are co-planar with this and the
175 two given points, see the B{Note} below.
177 @raise ImportError: Package C{numpy} not found, not installed or older
178 than version 1.10.
180 @raise IntersectionError: Near-coincident or -colinear points or
181 a trilateration or C{numpy} issue.
183 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
185 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted
186 back to geodetic lat-, longitude and height. The latter, conventionally
187 in C{meter} indicates whether the C{center} is above, below or on the
188 surface of the earth model. If C{deltas} is C{None}, the C{center} is
189 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,
190 height)} representing the differences between both results from
191 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
193 @see: Function L{pygeodesy.circin6}, method L{circum3}, U{Incircle
194 <https://MathWorld.Wolfram.com/Incircle.html>} and U{Contact Triangle
195 <https://MathWorld.Wolfram.com/ContactTriangle.html>}.
196 '''
197 try:
198 cs = self._toCartesian3(point2, point3)
199 r, c, d, cA, cB, cC = _circin6(*cs, eps=eps, useZ=True, dLL3=True,
200 datum=self.datum) # PYCHOK unpack
201 return Circin6Tuple(r, c.toLatLon(), d, cA.toLatLon(), cB.toLatLon(), cC.toLatLon())
202 except (AssertionError, TypeError, ValueError) as x:
203 raise _xError(x, point=self, point2=point2, point3=point3)
205 def circum3(self, point2, point3, circum=True, eps=EPS4):
206 '''Return the radius and center of the smallest circle I{through} or I{containing}
207 this and two other points.
209 @arg point2: Second point (C{LatLon}).
210 @arg point3: Third point (C{LatLon}).
211 @kwarg circum: If C{True} return the C{circumradius} and C{circumcenter},
212 always, ignoring the I{Meeus}' Type I case (C{bool}).
213 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.
215 @return: A L{Circum3Tuple}C{(radius, center, deltas)}. The C{center}, an
216 instance of this (sub-)class, is co-planar with this and the two
217 given points. If C{deltas} is C{None}, the C{center} is
218 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat,
219 lon, height)} representing the difference between both results
220 from L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
222 @raise ImportError: Package C{numpy} not found, not installed or older than
223 version 1.10.
225 @raise IntersectionError: Near-concentric, -coincident or -colinear points,
226 incompatible C{Ecef} classes or a trilateration
227 or C{numpy} issue.
229 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
231 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted
232 back to geodetic lat-, longitude and height. The latter, conventionally
233 in C{meter} indicates whether the C{center} is above, below or on the
234 surface of the earth model. If C{deltas} is C{None}, the C{center} is
235 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,
236 height)} representing the difference between both results from
237 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
239 @see: Function L{pygeodesy.circum3} and methods L{circin6} and L{circum4_}.
240 '''
241 try:
242 cs = self._toCartesian3(point2, point3)
243 r, c, d = _circum3(*cs, circum=circum, eps=eps, useZ=True, dLL3=True, # XXX -3d2
244 clas=cs[0].classof, datum=self.datum) # PYCHOK unpack
245 return Circum3Tuple(r, c.toLatLon(), d)
246 except (AssertionError, TypeError, ValueError) as x:
247 raise _xError(x, point=self, point2=point2, point3=point3, circum=circum)
249 def circum4_(self, *points):
250 '''Best-fit a sphere through this and two or more other points.
252 @arg points: The other points (each a C{LatLon}).
254 @return: L{Circum4Tuple}C{(radius, center, rank, residuals)} with C{center}
255 an instance of this (sub-)class.
257 @raise ImportError: Package C{numpy} not found, not installed or older than
258 version 1.10.
260 @raise NumPyError: Some C{numpy} issue.
262 @raise TypeError: One of the B{C{points}} invalid.
264 @raise ValueError: Too few B{C{points}}.
266 @see: Function L{pygeodesy.circum4_} and L{circum3}.
267 '''
268 C = self._toCartesianEcef
269 c = C(point=self)
270 t = circum4_(c, Vector=c.classof, *(C(i=i, points=p) for i, p in enumerate(points)))
271 c = t.center.toLatLon(LatLon=self.classof, name=t.name)
272 return Circum4Tuple(t.radius, c, t.rank, t.residuals, name=c.name)
274 @property
275 def clipid(self):
276 '''Get the (polygonal) clip (C{int}).
277 '''
278 return self._clipid
280 @clipid.setter # PYCHOK setter!
281 def clipid(self, clipid):
282 '''Get the (polygonal) clip (C{int}).
283 '''
284 self._clipid = int(clipid)
286 @deprecated_method
287 def compassAngle(self, other, adjust=True, wrap=False): # PYCHOK no cover
288 '''DEPRECATED, use method L{compassAngleTo}.'''
289 return self.compassAngleTo(other, adjust=adjust, wrap=wrap)
291 def compassAngleTo(self, other, adjust=True, wrap=False):
292 '''Return the angle from North for the direction vector between
293 this and an other point.
295 Suitable only for short, non-near-polar vectors up to a few
296 hundred Km or Miles. Use method C{initialBearingTo} for
297 larger distances.
299 @arg other: The other point (C{LatLon}).
300 @kwarg adjust: Adjust the longitudinal delta by the
301 cosine of the mean latitude (C{bool}).
302 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes and
303 longitudinal delta (C{bool}).
305 @return: Compass angle from North (C{degrees360}).
307 @raise TypeError: The B{C{other}} point is not C{LatLon}.
309 @note: Courtesy of Martin Schultz.
311 @see: U{Local, flat earth approximation
312 <https://www.EdWilliams.org/avform.htm#flat>}.
313 '''
314 self.others(other)
315 return compassAngle(self.lat, self.lon, other.lat, other.lon,
316 adjust=adjust, wrap=wrap)
318 def cosineAndoyerLambertTo(self, other, wrap=False):
319 '''Compute the distance between this and an other point using
320 the U{Andoyer-Lambert correction<https://navlib.net/wp-content/uploads/
321 2013/10/admiralty-manual-of-navigation-vol-1-1964-english501c.pdf>} of
322 the U{Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
323 formula.
325 @arg other: The other point (C{LatLon}).
326 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
328 @return: Distance (C{meter}, same units as the axes of
329 this point's datum ellipsoid).
331 @raise TypeError: The B{C{other}} point is not C{LatLon}.
333 @see: Function L{pygeodesy.cosineAndoyerLambert} and methods
334 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo},
335 C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},
336 L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo}, L{haversineTo},
337 L{thomasTo} and L{vincentysTo}.
338 '''
339 return self._distanceTo_(cosineAndoyerLambert_, other, wrap=wrap)
341 def cosineForsytheAndoyerLambertTo(self, other, wrap=False):
342 '''Compute the distance between this and an other point using
343 the U{Forsythe-Andoyer-Lambert correction
344 <https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} of the U{Law of Cosines
345 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
346 formula.
348 @arg other: The other point (C{LatLon}).
349 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
351 @return: Distance (C{meter}, same units as the axes of
352 this point's datum ellipsoid).
354 @raise TypeError: The B{C{other}} point is not C{LatLon}.
356 @see: Function L{pygeodesy.cosineForsytheAndoyerLambert} and methods
357 L{cosineAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
358 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
359 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo}.
360 '''
361 return self._distanceTo_(cosineForsytheAndoyerLambert_, other, wrap=wrap)
363 def cosineLawTo(self, other, radius=None, wrap=False):
364 '''Compute the distance between this and an other point using the
365 U{spherical Law of Cosines
366 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
367 formula.
369 @arg other: The other point (C{LatLon}).
370 @kwarg radius: Mean earth radius (C{meter}) or C{None}
371 for the mean radius of this point's datum
372 ellipsoid.
373 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
375 @return: Distance (C{meter}, same units as B{C{radius}}).
377 @raise TypeError: The B{C{other}} point is not C{LatLon}.
379 @see: Function L{pygeodesy.cosineLaw} and methods L{cosineAndoyerLambertTo},
380 L{cosineForsytheAndoyerLambertTo}, C{distanceTo*}, L{equirectangularTo},
381 L{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
382 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
383 '''
384 return self._distanceTo(cosineLaw, other, radius, wrap=wrap)
386 @property_RO
387 def datum(self): # PYCHOK no cover
388 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
389 '''
390 notOverloaded(self)
392 def destinationXyz(self, delta, LatLon=None, **LatLon_kwds):
393 '''Calculate the destination using a I{local} delta from this point.
395 @arg delta: Local delta to the destination (L{XyzLocal}, L{Enu},
396 L{Ned} or L{Local9Tuple}).
397 @kwarg LatLon: Optional (geodetic) class to return the destination
398 or C{None}.
399 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
400 arguments, ignored if C{B{LatLon} is None}.
402 @return: Destination as a C{B{LatLon}(lat, lon, **B{LatLon_kwds})}
403 instance or if C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat,
404 lon, height)} respectively L{LatLon4Tuple}C{(lat, lon,
405 height, datum)} depending on whether a C{datum} keyword
406 is un-/specified.
408 @raise TypeError: Invalid B{C{delta}}, B{C{LatLon}} or B{C{LatLon_kwds}}.
409 '''
410 t = self._ltp._local2ecef(delta, nine=True)
411 return t.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, name=self.name))
413 def _distanceTo(self, func, other, radius, **options):
414 '''(INTERNAL) Helper for methods C{<func>To}.
415 '''
416 self.others(other) # up=2
417 if radius is None:
418 radius = self._datum.ellipsoid.R1 if self._datum else R_M
419 return func(self.lat, self.lon, other.lat, other.lon,
420 radius=radius, **options)
422 def _distanceTo_(self, func_, other, wrap=False):
423 '''(INTERNAL) Helper for (ellipsoidal) methods C{<func>To}.
424 '''
425 self.others(other) # up=2
426 r, _ = unrollPI(self.lam, other.lam, wrap=wrap)
427 r = func_(other.phi, self.phi, r, datum=self.datum)
428 return r * self.datum.ellipsoid.a
430 @Property_RO
431 def Ecef(self):
432 '''Get the ECEF I{class} (L{EcefKarney}), I{lazily}.
433 '''
434 return _MODS.ecef.EcefKarney # default
436 @Property_RO
437 def _Ecef_forward(self):
438 '''(INTERNAL) Helper for L{_ecef9} and L{toEcef} (C{callable}).
439 '''
440 return self.Ecef(self.datum, name=self.name).forward
442 @Property_RO
443 def _ecef9(self):
444 '''(INTERNAL) Helper for L{toCartesian}, L{toEcef} and L{toCartesian} (L{Ecef9Tuple}).
445 '''
446 return self._Ecef_forward(self, M=True)
448 @deprecated_method
449 def equals(self, other, eps=None): # PYCHOK no cover
450 '''DEPRECATED, use method L{isequalTo}.'''
451 return self.isequalTo(other, eps=eps)
453 @deprecated_method
454 def equals3(self, other, eps=None): # PYCHOK no cover
455 '''DEPRECATED, use method L{isequalTo3}.'''
456 return self.isequalTo3(other, eps=eps)
458 def equirectangularTo(self, other, radius=None, **options):
459 '''Compute the distance between this and an other point
460 using the U{Equirectangular Approximation / Projection
461 <https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.
463 Suitable only for short, non-near-polar distances up to a
464 few hundred Km or Miles. Use method L{haversineTo} or
465 C{distanceTo*} for more accurate and/or larger distances.
467 See function L{pygeodesy.equirectangular_} for more details,
468 the available B{C{options}} and errors raised.
470 @arg other: The other point (C{LatLon}).
471 @kwarg radius: Mean earth radius (C{meter}) or C{None} for
472 the mean radius of this point's datum ellipsoid.
473 @kwarg options: Optional keyword arguments for function
474 L{pygeodesy.equirectangular}.
476 @return: Distance (C{meter}, same units as B{C{radius}}).
478 @raise TypeError: The B{C{other}} point is not C{LatLon}.
480 @see: Function L{pygeodesy.equirectangular} and methods L{cosineAndoyerLambertTo},
481 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
482 C{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
483 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
484 '''
485 return self._distanceTo(equirectangular, other, radius, **options)
487 def euclideanTo(self, other, radius=None, **options):
488 '''Approximate the C{Euclidian} distance between this and
489 an other point.
491 See function L{pygeodesy.euclidean} for the available B{C{options}}.
493 @arg other: The other point (C{LatLon}).
494 @kwarg radius: Mean earth radius (C{meter}) or C{None} for
495 the mean radius of this point's datum ellipsoid.
496 @kwarg options: Optional keyword arguments for function
497 L{pygeodesy.euclidean}.
499 @return: Distance (C{meter}, same units as B{C{radius}}).
501 @raise TypeError: The B{C{other}} point is not C{LatLon}.
503 @see: Function L{pygeodesy.euclidean} and methods L{cosineAndoyerLambertTo},
504 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
505 L{equirectangularTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
506 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
507 '''
508 return self._distanceTo(euclidean, other, radius, **options)
510 def flatLocalTo(self, other, radius=None, wrap=False):
511 '''Compute the distance between this and an other point using the
512 U{ellipsoidal Earth to plane projection
513 <https://WikiPedia.org/wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}
514 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.
516 @arg other: The other point (C{LatLon}).
517 @kwarg radius: Mean earth radius (C{meter}) or C{None} for the
518 equatorial radius of this point's datum/ellipsoid.
519 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
521 @return: Distance (C{meter}, same units as B{C{radius}}).
523 @raise TypeError: The B{C{other}} point is not C{LatLon}.
525 @raise ValueError: Invalid B{C{radius}}.
527 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}, methods
528 L{cosineAndoyerLambertTo}, L{cosineForsytheAndoyerLambertTo},
529 L{cosineLawTo}, C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},
530 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo} and
531 U{local, flat Earth approximation<https://www.edwilliams.org/avform.htm#flat>}.
532 '''
533 E = self.datum.ellipsoid
534 r = self._distanceTo_(flatLocal_, other, wrap=wrap) * E.a2_
535 a = E.a if radius in (None, 1, _1_0) else Radius(radius)
536 return r * a
538 hubenyTo = flatLocalTo # for Karl Hubeny
540 def flatPolarTo(self, other, radius=None, wrap=False):
541 '''Compute the distance between this and an other point using
542 the U{polar coordinate flat-Earth<https://WikiPedia.org/wiki/
543 Geographical_distance#Polar_coordinate_flat-Earth_formula>}formula.
545 @arg other: The other point (C{LatLon}).
546 @kwarg radius: Mean earth radius (C{meter}) or C{None} for the
547 mean radius of this point's datum ellipsoid.
548 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
550 @return: Distance (C{meter}, same units as B{C{radius}}).
552 @raise TypeError: The B{C{other}} point is not C{LatLon}.
554 @see: Function L{pygeodesy.flatPolar} and methods L{cosineAndoyerLambertTo},
555 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
556 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
557 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
558 '''
559 return self._distanceTo(flatPolar, other, radius, wrap=wrap)
561 def hartzell(self, los=None, earth=None):
562 '''Compute the intersection of a Line-Of-Sight (los) from this Point-Of-View
563 (pov) with this point's ellipsoid surface.
565 @kwarg los: Line-Of-Sight, I{direction} to earth (L{Vector3d}) or
566 C{None} to point to the ellipsoid's center.
567 @kwarg earth: The earth model (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},
568 L{a_f2Tuple} or C{scalar} radius in C{meter}) overriding
569 this point's C{datum} ellipsoid.
571 @return: The ellipsoid intersection (C{LatLon}) or this very instance
572 if this C{pov's height} is C{0}.
574 @raise IntersectionError: Null C{pov} or B{C{los}} vector, this
575 C{pov's height} is negative or B{C{los}}
576 points outside the ellipsoid or in an
577 opposite direction.
579 @raise TypeError: Invalid B{C{los}}.
581 @see: Function C{hartzell} for further details.
582 '''
583 h = self.height
584 if not h:
585 r = self
586 elif h < 0:
587 raise IntersectionError(pov=self, los=los, height=h, txt=_no_(_height_))
588 elif los is None:
589 d = self.datum if earth is None else _spherical_datum(earth)
590 r = self.dup(datum=d, height=0, name=self.hartzell.__name__)
591 else:
592 c = self.toCartesian()
593 r = hartzell(c, los=los, earth=earth or self.datum, LatLon=self.classof)
594 return r
596 def haversineTo(self, other, radius=None, wrap=False):
597 '''Compute the distance between this and an other point using the
598 U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}
599 formula.
601 @arg other: The other point (C{LatLon}).
602 @kwarg radius: Mean earth radius (C{meter}) or C{None} for
603 the mean radius of this point's datum ellipsoid.
604 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
606 @return: Distance (C{meter}, same units as B{C{radius}}).
608 @raise TypeError: The B{C{other}} point is not C{LatLon}.
610 @see: Function L{pygeodesy.haversine} and methods L{cosineAndoyerLambertTo},
611 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
612 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
613 L{flatPolarTo}, L{thomasTo} and L{vincentysTo}.
614 '''
615 return self._distanceTo(haversine, other, radius, wrap=wrap)
617 def _havg(self, other, f=_0_5):
618 '''(INTERNAL) Weighted, average height.
620 @arg other: An other point (C{LatLon}).
621 @kwarg f: Optional fraction (C{float}).
623 @return: Average, fractional height (C{float}).
624 '''
625 return _MODS.fmath.favg(self.height, other.height, f=f)
627 @Property
628 def height(self):
629 '''Get the height (C{meter}).
630 '''
631 return self._height
633 @height.setter # PYCHOK setter!
634 def height(self, height):
635 '''Set the height (C{meter}).
637 @raise TypeError: Invalid B{C{height}} C{type}.
639 @raise ValueError: Invalid B{C{height}}.
640 '''
641 h = Height(height)
642 if self._height != h:
643 _update_all(self)
644 self._height = h
646 def height4(self, earth=None, normal=True, LatLon=None, **LatLon_kwds):
647 '''Compute the height above or below and the projection on this datum's
648 ellipsoid surface.
650 @kwarg earth: A datum, ellipsoid, triaxial ellipsoid or earth radius
651 I{overriding} this datum (L{Datum}, L{Ellipsoid},
652 L{Ellipsoid2}, L{a_f2Tuple}, L{Triaxial}, L{Triaxial_},
653 L{JacobiConformal} or C{meter}, conventionally).
654 @kwarg normal: If C{True} the projection is the nearest point on the
655 ellipsoid's surface, otherwise the intersection of the
656 radial line to the center and the ellipsoid's surface.
657 @kwarg LatLon: Optional class to return the height and projection
658 (C{LatLon}) or C{None}.
659 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
660 ignored if C{B{LatLon} is None}.
662 @note: Use keyword argument C{height=0} to override C{B{LatLon}.height}
663 to {0} or any other C{scalar}, conventionally in C{meter}.
665 @return: An instance of B{C{LatLon}} or if C{B{LatLon} is None}, a
666 L{Vector4Tuple}C{(x, y, z, h)} with the I{projection} C{x}, C{y}
667 and C{z} coordinates and height C{h} in C{meter}, conventionally.
669 @raise TriaxialError: No convergence in triaxial root finding.
671 @raise TypeError: Invalid B{C{earth}}.
673 @see: L{Ellipsoid.height4} and L{Triaxial_.height4} for more information.
674 '''
675 c = self.toCartesian()
676 if LatLon is None:
677 r = c.height4(earth=earth, normal=normal)
678 else:
679 r = c.height4(earth=earth, normal=normal, Cartesian=c.classof, height=0)
680 r = r.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, height=r.height))
681 return r
683 def heightStr(self, prec=-2, m=_m_):
684 '''Return this B{C{height}} as C{str}ing.
686 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
687 @kwarg m: Optional unit of the height (C{str}).
689 @see: Function L{pygeodesy.hstr}.
690 '''
691 return hstr(self.height, prec=prec, m=m)
693 @deprecated_method
694 def isantipode(self, other, eps=EPS): # PYCHOK no cover
695 '''DEPRECATED, use method L{isantipodeTo}.'''
696 return self.isantipodeTo(other, eps=eps)
698 def isantipodeTo(self, other, eps=EPS):
699 '''Check whether this and an other point are antipodal,
700 on diametrically opposite sides of the earth.
702 @arg other: The other point (C{LatLon}).
703 @kwarg eps: Tolerance for near-equality (C{degrees}).
705 @return: C{True} if points are antipodal within the given
706 tolerance, C{False} otherwise.
707 '''
708 return isantipode(self.lat, self.lon,
709 other.lat, other.lon, eps=eps)
711 @Property_RO
712 def isEllipsoidal(self):
713 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).
714 '''
715 return self.datum.isEllipsoidal if self._datum else None
717 @Property_RO
718 def isEllipsoidalLatLon(self):
719 '''Get C{LatLon} base.
720 '''
721 return False
723 def isequalTo(self, other, eps=None):
724 '''Compare this point with an other point, I{ignoring} height.
726 @arg other: The other point (C{LatLon}).
727 @kwarg eps: Tolerance for equality (C{degrees}).
729 @return: C{True} if both points are identical,
730 I{ignoring} height, C{False} otherwise.
732 @raise TypeError: The B{C{other}} point is not C{LatLon}
733 or mismatch of the B{C{other}} and
734 this C{class} or C{type}.
736 @raise UnitError: Invalid B{C{eps}}.
738 @see: Method L{isequalTo3}.
739 '''
740 self.others(other)
742 return _isequalTo(self, other, eps=Scalar_(eps=eps)) if eps else \
743 (self.lat == other.lat and self.lon == other.lon)
745 def isequalTo3(self, other, eps=None):
746 '''Compare this point with an other point, I{including} height.
748 @arg other: The other point (C{LatLon}).
749 @kwarg eps: Tolerance for equality (C{degrees}).
751 @return: C{True} if both points are identical
752 I{including} height, C{False} otherwise.
754 @raise TypeError: The B{C{other}} point is not C{LatLon}
755 or mismatch of the B{C{other}} and
756 this C{class} or C{type}.
758 @see: Method L{isequalTo}.
759 '''
760 return self.height == other.height and self.isequalTo(other, eps=eps)
762 @Property_RO
763 def isnormal(self):
764 '''Return C{True} if this point is normal (C{bool}),
765 meaning C{abs(lat) <= 90} and C{abs(lon) <= 180}.
767 @see: Functions L{pygeodesy.isnormal} and
768 L{pygeodesy.normal}.
769 '''
770 return isnormal(self.lat, self.lon, eps=0)
772 @Property_RO
773 def isSpherical(self):
774 '''Check whether this point is spherical (C{bool} or C{None} if unknown).
775 '''
776 return self.datum.isSpherical if self._datum else None
778 @Property_RO
779 def lam(self):
780 '''Get the longitude (B{C{radians}}).
781 '''
782 return radians(self.lon)
784 @Property
785 def lat(self):
786 '''Get the latitude (C{degrees90}).
787 '''
788 return self._lat
790 @lat.setter # PYCHOK setter!
791 def lat(self, lat):
792 '''Set the latitude (C{str[N|S]} or C{degrees}).
794 @raise ValueError: Invalid B{C{lat}}.
795 '''
796 lat = Lat(lat) # parseDMS(lat, suffix=_NS_, clip=90)
797 if self._lat != lat:
798 _update_all(self)
799 self._lat = lat
801 @Property
802 def latlon(self):
803 '''Get the lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).
804 '''
805 return LatLon2Tuple(self._lat, self._lon, name=self.name)
807 @latlon.setter # PYCHOK setter!
808 def latlon(self, latlonh):
809 '''Set the lat- and longitude and optionally the height
810 (2- or 3-tuple or comma- or space-separated C{str}
811 of C{degrees90}, C{degrees180} and C{meter}).
813 @raise TypeError: Height of B{C{latlonh}} not C{scalar} or
814 B{C{latlonh}} not C{list} or C{tuple}.
816 @raise ValueError: Invalid B{C{latlonh}} or M{len(latlonh)}.
818 @see: Function L{pygeodesy.parse3llh} to parse a B{C{latlonh}}
819 string into a 3-tuple (lat, lon, h).
820 '''
821 if isstr(latlonh):
822 latlonh = parse3llh(latlonh, height=self.height)
823 else:
824 _xinstanceof(list, tuple, latlonh=latlonh)
825 if len(latlonh) == 3:
826 h = Height(latlonh[2], name=Fmt.SQUARE(latlonh=2))
827 elif len(latlonh) != 2:
828 raise _ValueError(latlonh=latlonh)
829 else:
830 h = self.height
832 llh = Lat(latlonh[0]), Lon(latlonh[1]), h # parseDMS2(latlonh[0], latlonh[1])
833 if (self._lat, self._lon, self._height) != llh:
834 _update_all(self)
835 self._lat, self._lon, self._height = llh
837 def latlon2(self, ndigits=0):
838 '''Return this point's lat- and longitude in C{degrees}, rounded.
840 @kwarg ndigits: Number of (decimal) digits (C{int}).
842 @return: A L{LatLon2Tuple}C{(lat, lon)}, both C{float}
843 and rounded away from zero.
845 @note: The C{round}ed values are always C{float}, also
846 if B{C{ndigits}} is omitted.
847 '''
848 return LatLon2Tuple(round(self.lat, ndigits),
849 round(self.lon, ndigits), name=self.name)
851 @deprecated_method
852 def latlon_(self, ndigits=0): # PYCHOK no cover
853 '''DEPRECATED, use method L{latlon2}.'''
854 return self.latlon2(ndigits=ndigits)
856 latlon2round = latlon_ # PYCHOK no cover
858 @Property_RO
859 def latlonheight(self):
860 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).
861 '''
862 return self.latlon.to3Tuple(self.height)
864 @Property
865 def lon(self):
866 '''Get the longitude (C{degrees180}).
867 '''
868 return self._lon
870 @lon.setter # PYCHOK setter!
871 def lon(self, lon):
872 '''Set the longitude (C{str[E|W]} or C{degrees}).
874 @raise ValueError: Invalid B{C{lon}}.
875 '''
876 lon = Lon(lon) # parseDMS(lon, suffix=_EW_, clip=180)
877 if self._lon != lon:
878 _update_all(self)
879 self._lon = lon
881 @Property_RO
882 def _ltp(self):
883 '''(INTERNAL) Cache for L{toLtp}.
884 '''
885 return _MODS.ltp.Ltp(self, ecef=self.Ecef(self.datum), name=self.name)
887 def nearestOn6(self, points, closed=False, height=None, wrap=False):
888 '''Locate the point on a path or polygon closest to this point.
890 Points are converted to and distances are computed in
891 I{geocentric}, cartesian space.
893 @arg points: The path or polygon points (C{LatLon}[]).
894 @kwarg closed: Optionally, close the polygon (C{bool}).
895 @kwarg height: Optional height, overriding the height of
896 this and all other points (C{meter}). If
897 C{None}, take the height of points into
898 account for distances.
899 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes
900 (C{bool}).
902 @return: A L{NearestOn6Tuple}C{(closest, distance, fi, j,
903 start, end)} with the C{closest}, the C{start}
904 and the C{end} point each an instance of this
905 C{LatLon} and C{distance} in C{meter}, same
906 units as the cartesian axes.
908 @raise PointsError: Insufficient number of B{C{points}}.
910 @raise TypeError: Some B{C{points}} or some B{C{points}}'
911 C{Ecef} invalid.
913 @raise ValueError: Some B{C{points}}' C{Ecef} is incompatible.
915 @see: Function L{pygeodesy.nearestOn6}.
916 '''
917 def _cs(Ps, h, w, C):
918 p = None # not used
919 for i, q in Ps.enumerate():
920 if w and i != 0:
921 q = _unrollon(p, q)
922 yield C(height=h, i=i, up=3, points=q)
923 p = q
925 C = self._toCartesianEcef # to verify datum and Ecef
926 Ps = self.PointsIter(points)
928 c = C(height=height, this=self) # this Cartesian
929 t = nearestOn6(c, _cs(Ps, height, wrap, C), closed=closed)
930 c, s, e = t.closest, t.start, t.end
932 kwds = _xkwds_not(None, LatLon=self.classof, # this LatLon
933 height=height)
934 r = self.Ecef(self.datum).reverse
935 p = r(c).toLatLon(**kwds)
936 s = r(s).toLatLon(**kwds) if s is not c else p
937 e = r(e).toLatLon(**kwds) if e is not c else p
938 return t.dup(closest=p, start=s, end=e)
940 def normal(self):
941 '''Normalize this point to C{abs(lat) <= 90} and C{abs(lon) <= 180}.
943 @return: C{True} if this point was I{normal}, C{False}
944 if it wasn't (but is now).
946 @see: Property L{isnormal} and function L{pygeodesy.normal}.
947 '''
948 n = self.isnormal
949 if not n:
950 self.latlon = normal(self.lat, self.lon)
951 return n
953 @Property_RO
954 def _N_vector(self):
955 '''(INTERNAL) Get the (C{nvectorBase._N_vector_})
956 '''
957 return _MODS.nvectorBase._N_vector_(*self.xyzh)
959 @Property_RO
960 def phi(self):
961 '''Get the latitude (B{C{radians}}).
962 '''
963 return radians(self.lat)
965 @Property_RO
966 def philam(self):
967 '''Get the lat- and longitude (L{PhiLam2Tuple}C{(phi, lam)}).
968 '''
969 return PhiLam2Tuple(self.phi, self.lam, name=self.name)
971 def philam2(self, ndigits=0):
972 '''Return this point's lat- and longitude in C{radians}, rounded.
974 @kwarg ndigits: Number of (decimal) digits (C{int}).
976 @return: A L{PhiLam2Tuple}C{(phi, lam)}, both C{float}
977 and rounded away from zero.
979 @note: The C{round}ed values are always C{float}, also
980 if B{C{ndigits}} is omitted.
981 '''
982 return PhiLam2Tuple(round(self.phi, ndigits),
983 round(self.lam, ndigits), name=self.name)
985 @Property_RO
986 def philamheight(self):
987 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
988 '''
989 return self.philam.to3Tuple(self.height)
991 @deprecated_method
992 def points(self, points, closed=True): # PYCHOK no cover
993 '''DEPRECATED, use method L{points2}.'''
994 return self.points2(points, closed=closed)
996 def points2(self, points, closed=True):
997 '''Check a path or polygon represented by points.
999 @arg points: The path or polygon points (C{LatLon}[])
1000 @kwarg closed: Optionally, consider the polygon closed,
1001 ignoring any duplicate or closing final
1002 B{C{points}} (C{bool}).
1004 @return: A L{Points2Tuple}C{(number, points)}, an C{int}
1005 and C{list} or C{tuple}.
1007 @raise PointsError: Insufficient number of B{C{points}}.
1009 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1010 '''
1011 return points2(points, closed=closed, base=self)
1013 def PointsIter(self, points, loop=0, dedup=False):
1014 '''Return a C{PointsIter} iterator.
1016 @arg points: The path or polygon points (C{LatLon}[])
1017 @kwarg loop: Number of loop-back points (non-negative C{int}).
1018 @kwarg dedup: Skip duplicate points (C{bool}).
1020 @return: A new C{PointsIter} iterator.
1022 @raise PointsError: Insufficient number of B{C{points}}.
1023 '''
1024 return PointsIter(points, base=self, loop=loop, dedup=dedup)
1026 def radii11(self, point2, point3):
1027 '''Return the radii of the C{Circum-}, C{In-}, I{Soddy} and C{Tangent}
1028 circles of a (planar) triangle formed by this and two other points.
1030 @arg point2: Second point (C{LatLon}).
1031 @arg point3: Third point (C{LatLon}).
1033 @return: L{Radii11Tuple}C{(rA, rB, rC, cR, rIn, riS, roS, a, b, c, s)}.
1035 @raise IntersectionError: Near-coincident or -colinear points.
1037 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
1039 @see: Function L{pygeodesy.radii11}, U{Incircle
1040 <https://MathWorld.Wolfram.com/Incircle.html>}, U{Soddy Circles
1041 <https://MathWorld.Wolfram.com/SoddyCircles.html>} and U{Tangent
1042 Circles<https://MathWorld.Wolfram.com/TangentCircles.html>}.
1043 '''
1044 try:
1045 cs = self._toCartesian3(point2, point3)
1046 return _radii11ABC(*cs, useZ=True)[0]
1047 except (TypeError, ValueError) as x:
1048 raise _xError(x, point=self, point2=point2, point3=point3)
1050 def _rhumb2(self, exact, radius):
1051 '''(INTERNAL) Get the C{rhumb} for this point's datum or for
1052 the earth model or earth B{C{radius}} if not C{None}.
1053 '''
1054 D = self.datum if radius is None else _spherical_datum(radius) # ellipsoidal OK
1055 r = D.ellipsoid.rhumbx if exact else \
1056 _MODS.rhumbx.Rhumb(D, exact=False, name=D.name)
1057 return r, D
1059 def rhumbAzimuthTo(self, other, exact=False, radius=None):
1060 '''Return the azimuth (bearing) of a rhumb line (loxodrome)
1061 between this and an other (ellipsoidal) point.
1063 @arg other: The other point (C{LatLon}).
1064 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1065 default C{False}.
1066 @kwarg radius: Optional earth radius (C{meter}) or earth model
1067 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1068 L{a_f2Tuple}), overriding this point's datum.
1070 @return: Rhumb azimuth (compass C{degrees360}).
1072 @raise TypeError: The B{C{other}} point is incompatible or
1073 B{C{radius}} is invalid.
1074 '''
1075 self.others(other)
1076 r, _ = self._rhumb2(exact, radius)
1077 C = _MODS.rhumbx.Caps
1078 return r.Inverse(self.lat, self.lon, other.lat, other.lon,
1079 outmask=C.AZIMUTH).azi12
1081 def rhumbDestination(self, distance, azimuth, exact=False, radius=None, height=None):
1082 '''Return the destination point having travelled the given distance
1083 from this point along a rhumb line (loxodrome) at the given azimuth.
1085 @arg distance: Distance travelled (C{meter}, same units as this
1086 point's datum (ellipsoid) axes or B{C{radius}},
1087 may be negative.
1088 @arg azimuth: Azimuth (bearing) at this point (compass C{degrees}).
1089 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1090 default C{False}.
1091 @kwarg radius: Optional earth radius (C{meter}) or earth model
1092 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1093 L{a_f2Tuple}), overriding this point's datum.
1094 @kwarg height: Optional height, overriding the default height
1095 (C{meter}).
1097 @return: The destination point (ellipsoidal C{LatLon}).
1099 @raise TypeError: Invalid B{C{radius}}.
1101 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}},
1102 B{C{radius}} or B{C{height}}.
1103 '''
1104 r, D = self._rhumb2(exact, radius)
1105 d = r.Direct(self.lat, self.lon, azimuth, distance)
1106 h = self.height if height is None else Height(height)
1107 return self.classof(d.lat2, d.lon2, datum=D, height=h)
1109 def rhumbDistanceTo(self, other, exact=False, radius=None):
1110 '''Return the distance from this to an other point along
1111 a rhumb line (loxodrome).
1113 @arg other: The other point (C{LatLon}).
1114 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1115 default C{False}.
1116 @kwarg radius: Optional earth radius (C{meter}) or earth model
1117 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1118 L{a_f2Tuple}), overriding this point's datum.
1120 @return: Distance (C{meter}, the same units as this point's
1121 datum (ellipsoid) axes or B{C{radius}}.
1123 @raise TypeError: The B{C{other}} point is incompatible or
1124 B{C{radius}} is invalid.
1126 @raise ValueError: Invalid B{C{radius}}.
1127 '''
1128 self.others(other)
1129 r, _ = self._rhumb2(exact, radius)
1130 C = _MODS.rhumbx.Caps
1131 return r.Inverse(self.lat, self.lon, other.lat, other.lon,
1132 outmask=C.DISTANCE).s12
1134 def rhumbLine(self, azimuth_other, exact=False, radius=None, name=NN, **caps):
1135 '''Get a rhumb line through this point at a given azimuth or
1136 through this and an other point.
1138 @arg azimuth_other: The azimuth of the rhumb line (compass)
1139 C{degrees} or the other point (C{LatLon}).
1140 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1141 default C{False}.
1142 @kwarg radius: Optional earth radius (C{meter}) or earth model
1143 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1144 L{a_f2Tuple}), overriding this point's datum.
1145 @kwarg name: Optional name (C{str}).
1146 @kwarg caps: Optional C{caps}, see L{RhumbLine} C{B{caps}}.
1148 @return: A L{RhumbLine} instance.
1150 @raise TypeError: Invalid B{C{radius}} or BC{C{azimuth_other}}
1151 not a C{scalar} nor a C{LatLon}.
1153 @see: Classes L{RhumbLine} and L{Rhumb}, property L{Rhumb.exact}
1154 and methods L{Rhumb.DirectLine} and L{Rhumb.InverseLine}.
1155 '''
1156 r, _ = self._rhumb2(exact, radius)
1157 a = azimuth_other
1158 if isscalar(a):
1159 r = r.DirectLine(self.lat, self.lon, a,
1160 name=name or self.name, **caps)
1161 elif isinstance(a, LatLonBase):
1162 self.others(a)
1163 r = r.InverseLine(self.lat, self.lon, a.lat, a.lon,
1164 name=name or self.name, **caps)
1165 else:
1166 raise _TypeError(azimuth_other=a)
1167 return r
1169 def rhumbMidpointTo(self, other, exact=False, radius=None,
1170 height=None, fraction=_0_5):
1171 '''Return the (loxodromic) midpoint on the rhumb line between
1172 this and an other point.
1174 @arg other: The other point (C{LatLon}).
1175 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1176 default C{False}.
1177 @kwarg radius: Optional earth radius (C{meter}) or earth model
1178 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1179 L{a_f2Tuple}), overriding this point's datum.
1180 @kwarg height: Optional height, overriding the mean height
1181 (C{meter}).
1182 @kwarg fraction: Midpoint location from this point (C{scalar}),
1183 may be negative or greater than 1.0.
1185 @return: The midpoint at the given B{C{fraction}} along the
1186 rhumb line (C{LatLon}).
1188 @raise TypeError: The B{C{other}} point is incompatible or
1189 B{C{radius}} is invalid.
1191 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}.
1192 '''
1193 self.others(other)
1194 r, D = self._rhumb2(exact, radius)
1195 f = Scalar(fraction=fraction)
1196 d = r.Inverse(self.lat, self.lon, other.lat, other.lon)
1197 d = r.Direct( self.lat, self.lon, d.azi12, d.s12 * f)
1198 h = self._havg(other, f=f) if height is None else Height(height)
1199 return self.classof(d.lat2, d.lon2, datum=D, height=h)
1201 def thomasTo(self, other, wrap=False):
1202 '''Compute the distance between this and an other point using
1203 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}
1204 formula.
1206 @arg other: The other point (C{LatLon}).
1207 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
1209 @return: Distance (C{meter}, same units as the axes of
1210 this point's datum ellipsoid).
1212 @raise TypeError: The B{C{other}} point is not C{LatLon}.
1214 @see: Function L{pygeodesy.thomas} and methods L{cosineAndoyerLambertTo},
1215 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
1216 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
1217 L{flatPolarTo}, L{haversineTo} and L{vincentysTo}.
1218 '''
1219 return self._distanceTo_(thomas_, other, wrap=wrap)
1221 @deprecated_method
1222 def to2ab(self): # PYCHOK no cover
1223 '''DEPRECATED, use property L{philam}.'''
1224 return self.philam
1226 def toCartesian(self, height=None, Cartesian=None, **Cartesian_kwds):
1227 '''Convert this point to cartesian, I{geocentric} coordinates,
1228 also known as I{Earth-Centered, Earth-Fixed} (ECEF).
1230 @kwarg height: Optional height, overriding this point's height
1231 (C{meter}, conventionally).
1232 @kwarg Cartesian: Optional class to return the geocentric
1233 coordinates (C{Cartesian}) or C{None}.
1234 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}}
1235 keyword arguments, ignored if
1236 C{B{Cartesian} is None}.
1238 @return: A B{C{Cartesian}} or if B{C{Cartesian}} is C{None},
1239 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
1240 datum)} with C{C=0} and C{M} if available.
1242 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}.
1243 '''
1244 r = self._ecef9 if height is None else self.toEcef(height=height)
1245 if Cartesian is not None: # class or .classof
1246 r = self._xnamed(Cartesian(r, **Cartesian_kwds))
1247 _xdatum(r.datum, self.datum)
1248 return r
1250 def _toCartesian3(self, point2, point3):
1251 '''(INTERNAL) Convert this and 2 other points.
1252 '''
1253 return (self. toCartesian().copy(name=_point_), # copy to rename
1254 self._toCartesianEcef(up=3, point2=point2),
1255 self._toCartesianEcef(up=3, point3=point3))
1257 def _toCartesianEcef(self, height=None, i=None, up=2, **name_point):
1258 '''(INTERNAL) Convert to cartesian and check Ecef's before and after.
1259 '''
1260 p = self.others(up=up, **name_point)
1261 c = p.toCartesian(height=height)
1262 E = self.Ecef
1263 if E:
1264 for p in (p, c):
1265 e = getattr(p, LatLonBase.Ecef.name, None)
1266 if e not in (None, E): # PYCHOK no cover
1267 n, _ = name_point.popitem()
1268 if i is not None:
1269 Fmt.SQUARE(n, i)
1270 raise _ValueError(n, e, txt=_incompatible(E.__name__))
1271 return c
1273 def toEcef(self, height=None, M=False):
1274 '''Convert this point to I{geocentric} coordinates, also known as
1275 I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).
1277 @kwarg height: Optional height, overriding this point's height
1278 (C{meter}, conventionally).
1279 @kwarg M: Optionally, include the rotation L{EcefMatrix} (C{bool}).
1281 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
1282 with C{C=0} and C{M} if available.
1284 @raise EcefError: A C{.datum} or an ECEF issue.
1285 '''
1286 return self._ecef9 if height in (None, self.height) else \
1287 self._Ecef_forward(self.lat, self.lon, height=height, M=M)
1289 @deprecated_method
1290 def to3llh(self, height=None): # PYCHOK no cover
1291 '''DEPRECATED, use property L{latlonheight} or C{latlon.to3Tuple(B{height})}.'''
1292 return self.latlonheight if height in (None, self.height) else \
1293 self.latlon.to3Tuple(height)
1295 def toLocal(self, Xyz=None, ltp=None, **Xyz_kwds):
1296 '''Convert this I{geodetic} point to I{local} C{X}, C{Y} and C{Z}.
1298 @kwarg Xyz: Optional class to return C{X}, C{Y} and C{Z}
1299 (L{XyzLocal}, L{Enu}, L{Ned}) or C{None}.
1300 @kwarg ltp: The I{local tangent plane} (LTP) to use,
1301 overriding this point's LTP (L{Ltp}).
1302 @kwarg Xyz_kwds: Optional, additional B{C{Xyz}} keyword
1303 arguments, ignored if C{B{Xyz} is None}.
1305 @return: An B{C{Xyz}} instance or if C{B{Xyz} is None},
1306 a L{Local9Tuple}C{(x, y, z, lat, lon, height,
1307 ltp, ecef, M)} with C{M=None}, always.
1309 @raise TypeError: Invalid B{C{ltp}}.
1310 '''
1311 p = _MODS.ltp._xLtp(ltp, self._ltp)
1312 return p._ecef2local(self._ecef9, Xyz, Xyz_kwds)
1314 def toLtp(self, Ecef=None):
1315 '''Return the I{local tangent plane} (LTP) for this point.
1317 @kwarg Ecef: Optional ECEF I{class} (L{EcefKarney}, ...
1318 L{EcefYou}), overriding this point's C{Ecef}.
1319 '''
1320 return self._ltp if Ecef in (None, self.Ecef) else _MODS.ltp.Ltp(
1321 self, ecef=Ecef(self.datum), name=self.name)
1323 def toNvector(self, h=None, Nvector=None, **Nvector_kwds):
1324 '''Convert this point to C{n-vector} (normal to the earth's
1325 surface) components, I{including height}.
1327 @kwarg h: Optional height, overriding this point's
1328 height (C{meter}).
1329 @kwarg Nvector: Optional class to return the C{n-vector}
1330 components (C{Nvector}) or C{None}.
1331 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}
1332 keyword arguments, ignored if
1333 C{B{Nvector} is None}.
1335 @return: A B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}
1336 if B{C{Nvector}} is C{None}.
1338 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.
1339 '''
1340 return self.toVector(Vector=Nvector, h=self.height if h is None else h,
1341 ll=self, **Nvector_kwds)
1343 def toStr(self, form=F_DMS, joined=_COMMASPACE_, m=_m_, **prec_sep_s_D_M_S): # PYCHOK expected
1344 '''Convert this point to a "lat, lon[, +/-height]" string, formatted
1345 in the given C{B{form}at}.
1347 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see
1348 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.
1349 @kwarg joined: Separator to join the lat-, longitude and heigth
1350 strings (C{str} or C{None} or C{NN} for non-joined).
1351 @kwarg m: Optional unit of the height (C{str}), use C{None} to
1352 exclude height from the returned string.
1353 @kwarg prec_sep_s_D_M_S: Optional C{B{prec}ision}, C{B{sep}arator},
1354 B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}} keyword
1355 arguments, see function L{pygeodesy.latDMS} or
1356 L{pygeodesy.lonDMS}.
1358 @return: This point in the specified C{B{form}at}, etc. (C{str} or
1359 a 2- or 3-tuple C{(lat_str, lon_str[, height_str])} if
1360 C{B{joined}=NN} or C{B{joined}=None}).
1362 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more
1363 details about keyword arguments C{B{form}at}, C{B{prec}ision},
1364 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.
1366 @example:
1368 >>> LatLon(51.4778, -0.0016).toStr() # 51°28′40″N, 000°00′06″W
1369 >>> LatLon(51.4778, -0.0016).toStr(F_D) # 51.4778°N, 000.0016°W
1370 >>> LatLon(51.4778, -0.0016, 42).toStr() # 51°28′40″N, 000°00′06″W, +42.00m
1371 '''
1372 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),
1373 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))
1374 if self.height and m is not None:
1375 t += (self.heightStr(m=m),)
1376 return joined.join(t) if joined else t
1378 def toVector(self, Vector=None, **Vector_kwds):
1379 '''Convert this point to C{n-vector} (normal to the earth's
1380 surface) components, I{ignoring height}.
1382 @kwarg Vector: Optional class to return the C{n-vector}
1383 components (L{Vector3d}) or C{None}.
1384 @kwarg Vector_kwds: Optional, additional B{C{Vector}}
1385 keyword arguments, ignored if
1386 C{B{Vector} is None}.
1388 @return: A B{C{Vector}} or a L{Vector3Tuple}C{(x, y, z)}
1389 if B{C{Vector}} is C{None}.
1391 @raise TypeError: Invalid B{C{Vector}} or B{C{kwds}}.
1393 @note: These are C{n-vector} x, y and z components,
1394 I{NOT} geocentric (ECEF) x, y and z coordinates!
1395 '''
1396 r = self._vector3tuple
1397 if Vector is not None:
1398 r = Vector(*r, **_xkwds(Vector_kwds, name=self.name))
1399 return r
1401 def toVector3d(self):
1402 '''Convert this point to C{n-vector} (normal to the earth's
1403 surface) components, I{ignoring height}.
1405 @return: Unit vector (L{Vector3d}).
1407 @note: These are C{n-vector} x, y and z components,
1408 I{NOT} geocentric (ECEF) x, y and z coordinates!
1409 '''
1410 return self._vector3d # XXX .unit()
1412 @deprecated_method
1413 def to3xyz(self): # PYCHOK no cover
1414 '''DEPRECATED, use property L{xyz} or method L{toNvector}, L{toVector},
1415 L{toVector3d} or perhaps (geocentric) L{toEcef}.'''
1416 return self.xyz # self.toVector()
1418 @Property_RO
1419 def _vector3d(self):
1420 '''(INTERNAL) Cache for L{toVector3d}.
1421 '''
1422 return self.toVector(Vector=Vector3d) # XXX .unit()
1424 @Property_RO
1425 def _vector3tuple(self):
1426 '''(INTERNAL) Cache for L{toVector}.
1427 '''
1428 return philam2n_xyz(self.phi, self.lam, name=self.name)
1430 def vincentysTo(self, other, radius=None, wrap=False):
1431 '''Compute the distance between this and an other point using
1432 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}
1433 spherical formula.
1435 @arg other: The other point (C{LatLon}).
1436 @kwarg radius: Mean earth radius (C{meter}) or C{None}
1437 for the mean radius of this point's datum
1438 ellipsoid.
1439 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
1441 @return: Distance (C{meter}, same units as B{C{radius}}).
1443 @raise TypeError: The B{C{other}} point is not C{LatLon}.
1445 @see: Function L{pygeodesy.vincentys} and methods L{cosineAndoyerLambertTo},
1446 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
1447 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
1448 L{flatPolarTo}, L{haversineTo} and L{thomasTo}.
1449 '''
1450 return self._distanceTo(vincentys, other, radius, wrap=wrap)
1452 @Property_RO
1453 def xyz(self):
1454 '''Get the C{n-vector} X, Y and Z components (L{Vector3Tuple}C{(x, y, z)})
1456 @note: These are C{n-vector} x, y and z components, I{NOT}
1457 geocentric (ECEF) x, y and z coordinates!
1458 '''
1459 return self.toVector(Vector=Vector3Tuple)
1461 @Property_RO
1462 def xyzh(self):
1463 '''Get the C{n-vector} X, Y, Z and H components (L{Vector4Tuple}C{(x, y, z, h)})
1465 @note: These are C{n-vector} x, y and z components, I{NOT}
1466 geocentric (ECEF) x, y and z coordinates!
1467 '''
1468 return self.xyz.to4Tuple(self.height)
1471def _isequalTo(point1, point2, eps=EPS): # in .ellipsoidalBaseDI._intersect3._on, .formy
1472 '''(INTERNAL) Compare point lat-/lon without type.
1473 '''
1474 return fabs(point1.lat - point2.lat) <= eps and \
1475 fabs(point1.lon - point2.lon) <= eps
1478def _isequalTo_(point1, point2, eps=EPS): # PYCHOK in .formy
1479 '''(INTERNAL) Compare point phi-/lam without type.
1480 '''
1481 return fabs(point1.phi - point2.phi) <= eps and \
1482 fabs(point1.lam - point2.lam) <= eps
1485def _trilaterate5(p1, d1, p2, d2, p3, d3, area=True, eps=EPS1,
1486 radius=R_M, wrap=False):
1487 '''(INTERNAL) Trilaterate three points by area overlap or by
1488 perimeter intersection of three circles.
1490 @note: The B{C{radius}} is only needed for both the n-vectorial
1491 and C{sphericalTrigonometry.LatLon.distanceTo} methods and
1492 silently ignored by the C{ellipsoidalExact}, C{-GeodSolve},
1493 C{-Karney} and C{-Vincenty.LatLon.distanceTo} methods.
1494 '''
1495 r1 = Distance_(distance1=d1)
1496 r2 = Distance_(distance2=d2)
1497 r3 = Distance_(distance3=d3)
1499 m = 0 if area else (r1 + r2 + r3)
1500 pc = 0
1501 t = []
1502 for _ in range(3):
1503 try: # intersection of circle (p1, r1) and (p2, r2)
1504 c1, c2 = p1.intersections2(r1, p2, r2, wrap=wrap)
1506 if area: # check overlap
1507 if c1 is c2: # abutting
1508 c = c1
1509 else: # nearest point on radical
1510 c = p3.nearestOn(c1, c2, within=True, wrap=wrap)
1511 d = r3 - p3.distanceTo(c, radius=radius, wrap=wrap)
1512 if d > eps: # sufficient overlap
1513 t.append((d, c))
1514 m = max(m, d)
1516 else: # check intersection
1517 for c in ((c1,) if c1 is c2 else (c1, c2)):
1518 d = fabs(r3 - p3.distanceTo(c, radius=radius, wrap=wrap))
1519 if d < eps: # below margin
1520 t.append((d, c))
1521 m = min(m, d)
1523 except IntersectionError as x:
1524 if _concentric_ in str(x): # XXX ConcentricError?
1525 pc += 1
1527 p1, r1, p2, r2, p3, r3 = p2, r2, p3, r3, p1, r1 # rotate
1529 if t: # get min, max, points and count ...
1530 t = tuple(sorted(t))
1531 n = (len(t),) # as 1-tuple
1532 # ... or for a single trilaterated result,
1533 # min *is* max, min- *is* maxPoint and n=1
1534 return Trilaterate5Tuple(t[0] + t[-1] + n) # *(t[0] + ...)
1536 elif area and pc == 3: # all pairwise concentric ...
1537 r, p = min((r1, p1), (r2, p2), (r3, p3))
1538 m = max(r1, r2, r3)
1539 # ... return "smallest" point twice, the smallest
1540 # and largest distance and n=0 for concentric
1541 return Trilaterate5Tuple(float(r), p, float(m), p, 0)
1543 n, f = (_overlap_, max) if area else (_intersection_, min)
1544 t = '%s (%s %.3f)' % (_no_(n), f.__name__, m)
1545 raise IntersectionError(area=area, eps=eps, wrap=wrap, txt=t)
1548__all__ += _ALL_DOCS(LatLonBase)
1550# **) MIT License
1551#
1552# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1553#
1554# Permission is hereby granted, free of charge, to any person obtaining a
1555# copy of this software and associated documentation files (the "Software"),
1556# to deal in the Software without restriction, including without limitation
1557# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1558# and/or sell copies of the Software, and to permit persons to whom the
1559# Software is furnished to do so, subject to the following conditions:
1560#
1561# The above copyright notice and this permission notice shall be included
1562# in all copies or substantial portions of the Software.
1563#
1564# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1565# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1566# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1567# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1568# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1569# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1570# OTHER DEALINGS IN THE SOFTWARE.