Coverage for pygeodesy/latlonBase.py: 93%
430 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -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, map1, _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, _IsnotError, \
19 _TypeError, _ValueError, _xdatum, _xError, \
20 _xkwds, _xkwds_not
21from pygeodesy.formy import antipode, compassAngle, cosineAndoyerLambert_, \
22 cosineForsytheAndoyerLambert_, cosineLaw, \
23 equirectangular, euclidean, flatLocal_, \
24 flatPolar, hartzell, haversine, isantipode, \
25 _isequalTo, isnormal, normal, philam2n_xyz, \
26 _spherical_datum, thomas_, vincentys
27from pygeodesy.interns import NN, _COMMASPACE_, _concentric_, _height_, \
28 _intersection_, _m_, _LatLon_, _no_, \
29 _overlap_, _point_ # PYCHOK used!
30from pygeodesy.iters import PointsIter, points2
31from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
32from pygeodesy.named import _NamedBase, notOverloaded
33from pygeodesy.namedTuples import Bounds2Tuple, LatLon2Tuple, PhiLam2Tuple, \
34 Trilaterate5Tuple, Vector3Tuple
35from pygeodesy.props import deprecated_method, Property, Property_RO, \
36 property_RO, _update_all
37from pygeodesy.streprs import Fmt, hstr
38from pygeodesy.units import Distance_, Lat, Lon, Height, Radius, Radius_, \
39 Scalar, Scalar_
40from pygeodesy.utily import _unrollon, _unrollon3, _Wrap
41from pygeodesy.vector2d import _circin6, Circin6Tuple, _circum3, circum4_, \
42 Circum3Tuple, _radii11ABC
43from pygeodesy.vector3d import nearestOn6, Vector3d
45from contextlib import contextmanager
46from math import asin, cos, degrees, fabs, radians
48__all__ = _ALL_LAZY.latlonBase
49__version__ = '23.05.12'
52class LatLonBase(_NamedBase):
53 '''(INTERNAL) Base class for C{LatLon} points on spherical or
54 ellipsoidal earth models.
55 '''
56 _clipid = INT0 # polygonal clip, see .booleans
57 _datum = None # L{Datum}, to be overriden
58 _height = INT0 # height (C{meter}), default
59 _lat = 0 # latitude (C{degrees})
60 _lon = 0 # longitude (C{degrees})
62 def __init__(self, latlonh, lon=None, height=0, wrap=False, name=NN):
63 '''New C{LatLon}.
65 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or
66 a previous C{LatLon} instance provided C{B{lon}=None}.
67 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or
68 C(None), indicating B{C{latlonh}} is a C{LatLon}.
69 @kwarg height: Optional height above (or below) the earth surface
70 (C{meter}, conventionally).
71 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}
72 (C{bool}).
73 @kwarg name: Optional name (C{str}).
75 @return: New instance (C{LatLon}).
77 @raise RangeError: A B{C{lon}} or C{lat} value outside the valid
78 range and L{rangerrors} set to C{True}.
80 @raise TypeError: If B{C{latlonh}} is not a C{LatLon}.
82 @raise UnitError: Invalid B{C{lat}}, B{C{lon}} or B{C{height}}.
84 @example:
86 >>> p = LatLon(50.06632, -5.71475)
87 >>> q = LatLon('50°03′59″N', """005°42'53"W""")
88 >>> r = LatLon(p)
89 '''
90 if name:
91 self.name = name
93 if lon is None:
94 try:
95 lat, lon = latlonh.lat, latlonh.lon
96 height = latlonh.get(_height_, height)
97 except AttributeError:
98 raise _IsnotError(_LatLon_, latlonh=latlonh)
99 if wrap:
100 lat, lon = _Wrap.latlon(lat, lon)
101 elif wrap:
102 lat, lon = _Wrap.latlonDMS2(latlonh, lon)
103 else:
104 lat = latlonh
106 self._lat = Lat(lat) # parseDMS2(lat, lon)
107 self._lon = Lon(lon) # PYCHOK LatLon2Tuple
108 if height: # elevation
109 self._height = Height(height)
111 def __eq__(self, other):
112 return self.isequalTo(other)
114 def __ne__(self, other):
115 return not self.isequalTo(other)
117 def __str__(self):
118 return self.toStr(form=F_D, prec=6)
120 def antipode(self, height=None):
121 '''Return the antipode, the point diametrically opposite
122 to this point.
124 @kwarg height: Optional height of the antipode (C{meter}),
125 this point's height otherwise.
127 @return: The antipodal point (C{LatLon}).
128 '''
129 h = self._heigHt(height)
130 return self.classof(*antipode(*self.latlon), height=h)
132 @deprecated_method
133 def bounds(self, wide, tall, radius=R_M): # PYCHOK no cover
134 '''DEPRECATED, use method C{boundsOf}.'''
135 return self.boundsOf(wide, tall, radius=radius)
137 def boundsOf(self, wide, tall, radius=R_M, height=None):
138 '''Return the SW and NE lat-/longitude of a great circle
139 bounding box centered at this location.
141 @arg wide: Longitudinal box width (C{meter}, same units as
142 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).
143 @arg tall: Latitudinal box size (C{meter}, same units as
144 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).
145 @kwarg radius: Mean earth radius (C{meter}) or C{None} if I{both}
146 B{C{wide}} and B{C{tall}} are in C{degrees}.
147 @kwarg height: Height for C{latlonSW} and C{latlonNE} (C{meter}),
148 overriding the point's height.
150 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)}, the
151 lower-left and upper-right corner (C{LatLon}).
153 @see: U{https://www.Movable-Type.co.UK/scripts/latlong-db.html}
154 '''
155 w = Scalar_(wide=wide) * _0_5
156 t = Scalar_(tall=tall) * _0_5
157 if radius is not None:
158 r = Radius_(radius)
159 c = cos(self.phi)
160 w = degrees(asin(w / r) / c) if fabs(c) > EPS0 else _0_0 # XXX
161 t = degrees(t / r)
162 y, t = self.lat, fabs(t)
163 x, w = self.lon, fabs(w)
165 h = self._heigHt(height)
166 sw = self.classof(y - t, x - w, height=h)
167 ne = self.classof(y + t, x + w, height=h)
168 return Bounds2Tuple(sw, ne, name=self.name)
170 def chordTo(self, other, height=None, wrap=False):
171 '''Compute the length of the chord through the earth between
172 this and an other point.
174 @arg other: The other point (C{LatLon}).
175 @kwarg height: Overriding height for both points (C{meter})
176 or C{None} for each point's height.
177 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{other}}
178 point (C{bool}).
180 @return: The chord length (conventionally C{meter}).
182 @raise TypeError: The B{C{other}} point is not C{LatLon}.
183 '''
184 def _v3d(ll):
185 t = ll.toEcef(height=height) # .toVector(Vector=Vector3d)
186 return Vector3d(t.x, t.y, t.z)
188 p = self.others(other)
189 if wrap:
190 p = _Wrap.point(p)
191 return _v3d(self).minus(_v3d(p)).length
193 def circin6(self, point2, point3, eps=EPS4, wrap=False):
194 '''Return the radius and center of the I{inscribed} aka I{In-}circle
195 of the (planar) triangle formed by this and two other points.
197 @arg point2: Second point (C{LatLon}).
198 @arg point3: Third point (C{LatLon}).
199 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.
200 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{point2}} and
201 B{C{point3}} (C{bool}).
203 @return: L{Circin6Tuple}C{(radius, center, deltas, cA, cB, cC)}. The
204 C{center} and contact points C{cA}, C{cB} and C{cC}, each an
205 instance of this (sub-)class, are co-planar with this and the
206 two given points, see the B{Note} below.
208 @raise ImportError: Package C{numpy} not found, not installed or older
209 than version 1.10.
211 @raise IntersectionError: Near-coincident or -colinear points or
212 a trilateration or C{numpy} issue.
214 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
216 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted
217 back to geodetic lat-, longitude and height. The latter, conventionally
218 in C{meter} indicates whether the C{center} is above, below or on the
219 surface of the earth model. If C{deltas} is C{None}, the C{center} is
220 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,
221 height)} representing the differences between both results from
222 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
224 @see: Function L{pygeodesy.circin6}, method L{circum3}, U{Incircle
225 <https://MathWorld.Wolfram.com/Incircle.html>} and U{Contact Triangle
226 <https://MathWorld.Wolfram.com/ContactTriangle.html>}.
227 '''
228 with _toCartesian3(self, point2, point3, wrap) as cs:
229 r, c, d, cA, cB, cC = _circin6(*cs, eps=eps, useZ=True, dLL3=True,
230 datum=self.datum) # PYCHOK unpack
231 return Circin6Tuple(r, c.toLatLon(), d, cA.toLatLon(), cB.toLatLon(), cC.toLatLon())
233 def circum3(self, point2, point3, circum=True, eps=EPS4, wrap=False):
234 '''Return the radius and center of the smallest circle I{through} or I{containing}
235 this and two other points.
237 @arg point2: Second point (C{LatLon}).
238 @arg point3: Third point (C{LatLon}).
239 @kwarg circum: If C{True} return the C{circumradius} and C{circumcenter},
240 always, ignoring the I{Meeus}' Type I case (C{bool}).
241 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.
242 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{point2}} and
243 B{C{point3}} (C{bool}).
245 @return: A L{Circum3Tuple}C{(radius, center, deltas)}. The C{center}, an
246 instance of this (sub-)class, is co-planar with this and the two
247 given points. If C{deltas} is C{None}, the C{center} is
248 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat,
249 lon, height)} representing the difference between both results
250 from L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
252 @raise ImportError: Package C{numpy} not found, not installed or older than
253 version 1.10.
255 @raise IntersectionError: Near-concentric, -coincident or -colinear points,
256 incompatible C{Ecef} classes or a trilateration
257 or C{numpy} issue.
259 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
261 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted
262 back to geodetic lat-, longitude and height. The latter, conventionally
263 in C{meter} indicates whether the C{center} is above, below or on the
264 surface of the earth model. If C{deltas} is C{None}, the C{center} is
265 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,
266 height)} representing the difference between both results from
267 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.
269 @see: Function L{pygeodesy.circum3} and methods L{circin6} and L{circum4_}.
270 '''
271 with _toCartesian3(self, point2, point3, wrap, circum=circum) as cs:
272 r, c, d = _circum3(*cs, circum=circum, eps=eps, useZ=True, dLL3=True, # XXX -3d2
273 clas=cs[0].classof, datum=self.datum) # PYCHOK unpack
274 return Circum3Tuple(r, c.toLatLon(), d)
276 def circum4_(self, *points, **wrap):
277 '''Best-fit a sphere through this and two or more other points.
279 @arg points: The other points (each a C{LatLon}).
280 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}}
281 (C{bool}), default C{False}.
283 @return: L{Circum4Tuple}C{(radius, center, rank, residuals)} with C{center}
284 an instance of this (sub-)class.
286 @raise ImportError: Package C{numpy} not found, not installed or older than
287 version 1.10.
289 @raise NumPyError: Some C{numpy} issue.
291 @raise TypeError: One of the B{C{points}} invalid.
293 @raise ValueError: Too few B{C{points}}.
295 @see: Function L{pygeodesy.circum4_} and L{circum3}.
296 '''
297 def _cs(ps, C, wrap=False):
298 _wp = _Wrap.point if wrap else (lambda p: p)
299 for i, p in enumerate(ps):
300 yield C(i=i, points=_wp(p))
302 C = self._toCartesianEcef
303 c = C(point=self)
304 t = circum4_(c, Vector=c.classof, *_cs(points, C, **wrap))
305 c = t.center.toLatLon(LatLon=self.classof)
306 return t.dup(center=c)
308 @property
309 def clipid(self):
310 '''Get the (polygonal) clip (C{int}).
311 '''
312 return self._clipid
314 @clipid.setter # PYCHOK setter!
315 def clipid(self, clipid):
316 '''Get the (polygonal) clip (C{int}).
317 '''
318 self._clipid = int(clipid)
320 @deprecated_method
321 def compassAngle(self, other, **adjust_wrap): # PYCHOK no cover
322 '''DEPRECATED, use method L{compassAngleTo}.'''
323 return self.compassAngleTo(other, **adjust_wrap)
325 def compassAngleTo(self, other, **adjust_wrap):
326 '''Return the angle from North for the direction vector between
327 this and an other point.
329 Suitable only for short, non-near-polar vectors up to a few
330 hundred Km or Miles. Use method C{initialBearingTo} for
331 larger distances.
333 @arg other: The other point (C{LatLon}).
334 @kwarg adjust_wrap: Optional keyword arguments for function
335 L{pygeodesy.compassAngle}.
337 @return: Compass angle from North (C{degrees360}).
339 @raise TypeError: The B{C{other}} point is not C{LatLon}.
341 @note: Courtesy of Martin Schultz.
343 @see: U{Local, flat earth approximation
344 <https://www.EdWilliams.org/avform.htm#flat>}.
345 '''
346 p = self.others(other)
347 return compassAngle(self.lat, self.lon, p.lat, p.lon, **adjust_wrap)
349 def cosineAndoyerLambertTo(self, other, wrap=False):
350 '''Compute the distance between this and an other point using the U{Andoyer-Lambert correction<https://
351 navlib.net/wp-content/uploads/2013/10/admiralty-manual-of-navigation-vol-1-1964-english501c.pdf>}
352 of the U{Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.
354 @arg other: The other point (C{LatLon}).
355 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
356 the B{C{other}} point (C{bool}).
358 @return: Distance (C{meter}, same units as the axes of this
359 point's datum ellipsoid).
361 @raise TypeError: The B{C{other}} point is not C{LatLon}.
363 @see: Function L{pygeodesy.cosineAndoyerLambert} and methods
364 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo},
365 C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},
366 L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo}, L{haversineTo},
367 L{thomasTo} and L{vincentysTo}.
368 '''
369 return self._distanceTo_(cosineAndoyerLambert_, other, wrap=wrap)
371 def cosineForsytheAndoyerLambertTo(self, other, wrap=False):
372 '''Compute the distance between this and an other point using
373 the U{Forsythe-Andoyer-Lambert correction
374 <https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} of the U{Law of Cosines
375 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
376 formula.
378 @arg other: The other point (C{LatLon}).
379 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
380 the B{C{other}} point (C{bool}).
382 @return: Distance (C{meter}, same units as the axes of
383 this point's datum ellipsoid).
385 @raise TypeError: The B{C{other}} point is not C{LatLon}.
387 @see: Function L{pygeodesy.cosineForsytheAndoyerLambert} and methods
388 L{cosineAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
389 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
390 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo}.
391 '''
392 return self._distanceTo_(cosineForsytheAndoyerLambert_, other, wrap=wrap)
394 def cosineLawTo(self, other, radius=None, wrap=False):
395 '''Compute the distance between this and an other point using the
396 U{spherical Law of Cosines
397 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}
398 formula.
400 @arg other: The other point (C{LatLon}).
401 @kwarg radius: Mean earth radius (C{meter}) or C{None}
402 for the mean radius of this point's datum
403 ellipsoid.
404 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
405 the B{C{other}} point (C{bool}).
407 @return: Distance (C{meter}, same units as B{C{radius}}).
409 @raise TypeError: The B{C{other}} point is not C{LatLon}.
411 @see: Function L{pygeodesy.cosineLaw} and methods L{cosineAndoyerLambertTo},
412 L{cosineForsytheAndoyerLambertTo}, C{distanceTo*}, L{equirectangularTo},
413 L{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
414 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
415 '''
416 return self._distanceTo(cosineLaw, other, radius, wrap=wrap)
418 @property_RO
419 def datum(self): # PYCHOK no cover
420 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
421 '''
422 notOverloaded(self)
424 def destinationXyz(self, delta, LatLon=None, **LatLon_kwds):
425 '''Calculate the destination using a I{local} delta from this point.
427 @arg delta: Local delta to the destination (L{XyzLocal}, L{Enu},
428 L{Ned} or L{Local9Tuple}).
429 @kwarg LatLon: Optional (geodetic) class to return the destination
430 or C{None}.
431 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
432 arguments, ignored if C{B{LatLon} is None}.
434 @return: Destination as a C{B{LatLon}(lat, lon, **B{LatLon_kwds})}
435 instance or if C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat,
436 lon, height)} respectively L{LatLon4Tuple}C{(lat, lon,
437 height, datum)} depending on whether a C{datum} keyword
438 is un-/specified.
440 @raise TypeError: Invalid B{C{delta}}, B{C{LatLon}} or B{C{LatLon_kwds}}.
441 '''
442 t = self._ltp._local2ecef(delta, nine=True)
443 return t.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, name=self.name))
445 def _distanceTo(self, func, other, radius=None, **kwds):
446 '''(INTERNAL) Helper for distance methods C{<func>To}.
447 '''
448 p, r = self.others(other, up=2), radius
449 if r is None:
450 r = self._datum.ellipsoid.R1 if self._datum else R_M
451 return func(self.lat, self.lon, p.lat, p.lon, radius=r, **kwds)
453 def _distanceTo_(self, func_, other, wrap=False, radius=None):
454 '''(INTERNAL) Helper for (ellipsoidal) methods C{<func>To}.
455 '''
456 p = self.others(other, up=2)
457 D = self.datum
458 lam21, phi2, _ = _Wrap.philam3(self.lam, p.phi, p.lam, wrap)
459 r = func_(phi2, self.phi, lam21, datum=D)
460 return r * (D.ellipsoid.a if radius is None else radius)
462 @Property_RO
463 def Ecef(self):
464 '''Get the ECEF I{class} (L{EcefKarney}), I{lazily}.
465 '''
466 return _MODS.ecef.EcefKarney # default
468 @Property_RO
469 def _Ecef_forward(self):
470 '''(INTERNAL) Helper for L{_ecef9} and L{toEcef} (C{callable}).
471 '''
472 return self.Ecef(self.datum, name=self.name).forward
474 @Property_RO
475 def _ecef9(self):
476 '''(INTERNAL) Helper for L{toCartesian}, L{toEcef} and L{toCartesian} (L{Ecef9Tuple}).
477 '''
478 return self._Ecef_forward(self, M=True)
480 @deprecated_method
481 def equals(self, other, eps=None): # PYCHOK no cover
482 '''DEPRECATED, use method L{isequalTo}.'''
483 return self.isequalTo(other, eps=eps)
485 @deprecated_method
486 def equals3(self, other, eps=None): # PYCHOK no cover
487 '''DEPRECATED, use method L{isequalTo3}.'''
488 return self.isequalTo3(other, eps=eps)
490 def equirectangularTo(self, other, **radius_adjust_limit_wrap):
491 '''Compute the distance between this and an other point
492 using the U{Equirectangular Approximation / Projection
493 <https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.
495 Suitable only for short, non-near-polar distances up to a
496 few hundred Km or Miles. Use method L{haversineTo} or
497 C{distanceTo*} for more accurate and/or larger distances.
499 @arg other: The other point (C{LatLon}).
500 @kwarg radius_adjust_limit_wrap: Optional keyword arguments
501 for function L{pygeodesy.equirectangular},
502 overriding the default mean C{radius} of this
503 point's datum ellipsoid.
505 @return: Distance (C{meter}, same units as B{C{radius}}).
507 @raise TypeError: The B{C{other}} point is not C{LatLon}.
509 @see: Function L{pygeodesy.equirectangular} and methods L{cosineAndoyerLambertTo},
510 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
511 C{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
512 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
513 '''
514 return self._distanceTo(equirectangular, other, **radius_adjust_limit_wrap)
516 def euclideanTo(self, other, **radius_adjust_wrap):
517 '''Approximate the C{Euclidian} distance between this and
518 an other point.
520 See function L{pygeodesy.euclidean} for the available B{C{options}}.
522 @arg other: The other point (C{LatLon}).
523 @kwarg radius_adjust_wrap: Optional keyword arguments for function
524 L{pygeodesy.euclidean}, overriding the default mean
525 C{radius} of this point's datum ellipsoid.
527 @return: Distance (C{meter}, same units as B{C{radius}}).
529 @raise TypeError: The B{C{other}} point is not C{LatLon}.
531 @see: Function L{pygeodesy.euclidean} and methods L{cosineAndoyerLambertTo},
532 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
533 L{equirectangularTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},
534 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
535 '''
536 return self._distanceTo(euclidean, other, **radius_adjust_wrap)
538 def flatLocalTo(self, other, radius=None, wrap=False):
539 '''Compute the distance between this and an other point using the
540 U{ellipsoidal Earth to plane projection
541 <https://WikiPedia.org/wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}
542 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.
544 @arg other: The other point (C{LatLon}).
545 @kwarg radius: Mean earth radius (C{meter}) or C{None} for
546 the I{equatorial radius} of this point's
547 datum ellipsoid.
548 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
549 the B{C{other}} point (C{bool}).
551 @return: Distance (C{meter}, same units as B{C{radius}}).
553 @raise TypeError: The B{C{other}} point is not C{LatLon}.
555 @raise ValueError: Invalid B{C{radius}}.
557 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}, methods
558 L{cosineAndoyerLambertTo}, L{cosineForsytheAndoyerLambertTo},
559 L{cosineLawTo}, C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},
560 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo} and
561 U{local, flat Earth approximation<https://www.edwilliams.org/avform.htm#flat>}.
562 '''
563 return self._distanceTo_(flatLocal_, other, wrap=wrap, radius=
564 radius if radius in (None, R_M, _1_0, 1) else Radius(radius)) # PYCHOK kwargs
566 hubenyTo = flatLocalTo # for Karl Hubeny
568 def flatPolarTo(self, other, **radius_wrap):
569 '''Compute the distance between this and an other point using
570 the U{polar coordinate flat-Earth<https://WikiPedia.org/wiki/
571 Geographical_distance#Polar_coordinate_flat-Earth_formula>}formula.
573 @arg other: The other point (C{LatLon}).
574 @kwarg radius_wrap: Optional keyword arguments for function
575 L{pygeodesy.flatPolar}, overriding the
576 default mean C{radius} of this point's
577 datum ellipsoid.
579 @return: Distance (C{meter}, same units as B{C{radius}}).
581 @raise TypeError: The B{C{other}} point is not C{LatLon}.
583 @see: Function L{pygeodesy.flatPolar} and methods L{cosineAndoyerLambertTo},
584 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
585 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
586 L{haversineTo}, L{thomasTo} and L{vincentysTo}.
587 '''
588 return self._distanceTo(flatPolar, other, **radius_wrap)
590 def hartzell(self, los=None, earth=None):
591 '''Compute the intersection of a Line-Of-Sight (los) from this Point-Of-View
592 (pov) with this point's ellipsoid surface.
594 @kwarg los: Line-Of-Sight, I{direction} to earth (L{Vector3d}) or
595 C{None} to point to the ellipsoid's center.
596 @kwarg earth: The earth model (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},
597 L{a_f2Tuple} or C{scalar} radius in C{meter}) overriding
598 this point's C{datum} ellipsoid.
600 @return: The ellipsoid intersection (C{LatLon}) or this very instance
601 if this C{pov's height} is C{0}.
603 @raise IntersectionError: Null C{pov} or B{C{los}} vector, this
604 C{pov's height} is negative or B{C{los}}
605 points outside the ellipsoid or in an
606 opposite direction.
608 @raise TypeError: Invalid B{C{los}}.
610 @see: Function C{hartzell} for further details.
611 '''
612 h = self.height
613 if not h:
614 r = self
615 elif h < 0:
616 raise IntersectionError(pov=self, los=los, height=h, txt=_no_(_height_))
617 elif los is None:
618 d = self.datum if earth is None else _spherical_datum(earth)
619 r = self.dup(datum=d, height=0, name=self.hartzell.__name__)
620 else:
621 c = self.toCartesian()
622 r = hartzell(c, los=los, earth=earth or self.datum, LatLon=self.classof)
623 return r
625 def haversineTo(self, other, **radius_wrap):
626 '''Compute the distance between this and an other point using the
627 U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}
628 formula.
630 @arg other: The other point (C{LatLon}).
631 @kwarg radius_wrap: Optional keyword arguments for function
632 L{pygeodesy.haversine}, overriding the
633 default mean C{radius} of this point's
634 datum ellipsoid.
636 @return: Distance (C{meter}, same units as B{C{radius}}).
638 @raise TypeError: The B{C{other}} point is not C{LatLon}.
640 @see: Function L{pygeodesy.haversine} and methods L{cosineAndoyerLambertTo},
641 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
642 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
643 L{flatPolarTo}, L{thomasTo} and L{vincentysTo}.
644 '''
645 return self._distanceTo(haversine, other, **radius_wrap)
647 def _havg(self, other, f=_0_5, h=None):
648 '''(INTERNAL) Weighted, average height.
650 @arg other: An other point (C{LatLon}).
651 @kwarg f: Optional fraction (C{float}).
652 @kwarg h: Overriding height (C{meter}).
654 @return: Average, fractional height (C{float}) or
655 the overriding B{C{height}} (C{Height}).
656 '''
657 return Height(h) if h is not None else \
658 _MODS.fmath.favg(self.height, other.height, f=f)
660 def _heigHt(self, height):
661 '''(INTERNAL) Overriding C{height}.
662 '''
663 return self.height if height is None else Height(height)
665 @Property
666 def height(self):
667 '''Get the height (C{meter}).
668 '''
669 return self._height
671 @height.setter # PYCHOK setter!
672 def height(self, height):
673 '''Set the height (C{meter}).
675 @raise TypeError: Invalid B{C{height}} C{type}.
677 @raise ValueError: Invalid B{C{height}}.
678 '''
679 h = Height(height)
680 if self._height != h:
681 _update_all(self)
682 self._height = h
684 def height4(self, earth=None, normal=True, LatLon=None, **LatLon_kwds):
685 '''Compute the height above or below and the projection of this point
686 on this datum's or on an other earth's ellipsoid surface.
688 @kwarg earth: A datum, ellipsoid, triaxial ellipsoid or earth radius
689 I{overriding} this datum (L{Datum}, L{Ellipsoid},
690 L{Ellipsoid2}, L{a_f2Tuple}, L{Triaxial}, L{Triaxial_},
691 L{JacobiConformal} or C{meter}, conventionally).
692 @kwarg normal: If C{True} the projection is the nearest point on the
693 ellipsoid's surface, otherwise the intersection of the
694 radial line to the center and the ellipsoid's surface.
695 @kwarg LatLon: Optional class to return the height and projection
696 (C{LatLon}) or C{None}.
697 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
698 ignored if C{B{LatLon} is None}.
700 @note: Use keyword argument C{height=0} to override C{B{LatLon}.height}
701 to {0} or any other C{scalar}, conventionally in C{meter}.
703 @return: An instance of B{C{LatLon}} or if C{B{LatLon} is None}, a
704 L{Vector4Tuple}C{(x, y, z, h)} with the I{projection} C{x}, C{y}
705 and C{z} coordinates and height C{h} in C{meter}, conventionally.
707 @raise TriaxialError: No convergence in triaxial root finding.
709 @raise TypeError: Invalid B{C{earth}}.
711 @see: L{Ellipsoid.height4} and L{Triaxial_.height4} for more information.
712 '''
713 c = self.toCartesian()
714 if LatLon is None:
715 r = c.height4(earth=earth, normal=normal)
716 else:
717 r = c.height4(earth=earth, normal=normal, Cartesian=c.classof, height=0)
718 r = r.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, height=r.height))
719 return r
721 def heightStr(self, prec=-2, m=_m_):
722 '''Return this point's B{C{height}} as C{str}ing.
724 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
725 @kwarg m: Optional unit of the height (C{str}).
727 @see: Function L{pygeodesy.hstr}.
728 '''
729 return hstr(self.height, prec=prec, m=m)
731 @deprecated_method
732 def isantipode(self, other, eps=EPS): # PYCHOK no cover
733 '''DEPRECATED, use method L{isantipodeTo}.'''
734 return self.isantipodeTo(other, eps=eps)
736 def isantipodeTo(self, other, eps=EPS):
737 '''Check whether this and an other point are antipodal,
738 on diametrically opposite sides of the earth.
740 @arg other: The other point (C{LatLon}).
741 @kwarg eps: Tolerance for near-equality (C{degrees}).
743 @return: C{True} if points are antipodal within the given
744 tolerance, C{False} otherwise.
745 '''
746 p = self.others(other)
747 return isantipode(*(self.latlon + p.latlon), eps=eps)
749 @Property_RO
750 def isEllipsoidal(self):
751 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).
752 '''
753 return self.datum.isEllipsoidal if self._datum else None
755 @Property_RO
756 def isEllipsoidalLatLon(self):
757 '''Get C{LatLon} base.
758 '''
759 return False
761 def isequalTo(self, other, eps=None):
762 '''Compare this point with an other point, I{ignoring} height.
764 @arg other: The other point (C{LatLon}).
765 @kwarg eps: Tolerance for equality (C{degrees}).
767 @return: C{True} if both points are identical,
768 I{ignoring} height, C{False} otherwise.
770 @raise TypeError: The B{C{other}} point is not C{LatLon}
771 or mismatch of the B{C{other}} and
772 this C{class} or C{type}.
774 @raise UnitError: Invalid B{C{eps}}.
776 @see: Method L{isequalTo3}.
777 '''
778 return _isequalTo(self, self.others(other), eps=eps)
780 def isequalTo3(self, other, eps=None):
781 '''Compare this point with an other point, I{including} height.
783 @arg other: The other point (C{LatLon}).
784 @kwarg eps: Tolerance for equality (C{degrees}).
786 @return: C{True} if both points are identical
787 I{including} height, C{False} otherwise.
789 @raise TypeError: The B{C{other}} point is not C{LatLon}
790 or mismatch of the B{C{other}} and
791 this C{class} or C{type}.
793 @see: Method L{isequalTo}.
794 '''
795 return self.height == self.others(other).height and \
796 _isequalTo(self, other, eps=eps)
798 @Property_RO
799 def isnormal(self):
800 '''Return C{True} if this point is normal (C{bool}),
801 meaning C{abs(lat) <= 90} and C{abs(lon) <= 180}.
803 @see: Methods L{normal}, L{toNormal} and functions
804 L{pygeodesy.isnormal} and L{pygeodesy.normal}.
805 '''
806 return isnormal(self.lat, self.lon, eps=0)
808 @Property_RO
809 def isSpherical(self):
810 '''Check whether this point is spherical (C{bool} or C{None} if unknown).
811 '''
812 return self.datum.isSpherical if self._datum else None
814 @Property_RO
815 def lam(self):
816 '''Get the longitude (B{C{radians}}).
817 '''
818 return radians(self.lon)
820 @Property
821 def lat(self):
822 '''Get the latitude (C{degrees90}).
823 '''
824 return self._lat
826 @lat.setter # PYCHOK setter!
827 def lat(self, lat):
828 '''Set the latitude (C{str[N|S]} or C{degrees}).
830 @raise ValueError: Invalid B{C{lat}}.
831 '''
832 lat = Lat(lat) # parseDMS(lat, suffix=_NS_, clip=90)
833 if self._lat != lat:
834 _update_all(self)
835 self._lat = lat
837 @Property
838 def latlon(self):
839 '''Get the lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).
840 '''
841 return LatLon2Tuple(self._lat, self._lon, name=self.name)
843 @latlon.setter # PYCHOK setter!
844 def latlon(self, latlonh):
845 '''Set the lat- and longitude and optionally the height
846 (2- or 3-tuple or comma- or space-separated C{str}
847 of C{degrees90}, C{degrees180} and C{meter}).
849 @raise TypeError: Height of B{C{latlonh}} not C{scalar} or
850 B{C{latlonh}} not C{list} or C{tuple}.
852 @raise ValueError: Invalid B{C{latlonh}} or M{len(latlonh)}.
854 @see: Function L{pygeodesy.parse3llh} to parse a B{C{latlonh}}
855 string into a 3-tuple C{(lat, lon, h)}.
856 '''
857 if isstr(latlonh):
858 latlonh = parse3llh(latlonh, height=self.height)
859 else:
860 _xinstanceof(list, tuple, latlonh=latlonh)
861 if len(latlonh) == 3:
862 h = Height(latlonh[2], name=Fmt.SQUARE(latlonh=2))
863 elif len(latlonh) != 2:
864 raise _ValueError(latlonh=latlonh)
865 else:
866 h = self.height
868 llh = Lat(latlonh[0]), Lon(latlonh[1]), h # parseDMS2(latlonh[0], latlonh[1])
869 if (self._lat, self._lon, self._height) != llh:
870 _update_all(self)
871 self._lat, self._lon, self._height = llh
873 def latlon2(self, ndigits=0):
874 '''Return this point's lat- and longitude in C{degrees}, rounded.
876 @kwarg ndigits: Number of (decimal) digits (C{int}).
878 @return: A L{LatLon2Tuple}C{(lat, lon)}, both C{float}
879 and rounded away from zero.
881 @note: The C{round}ed values are always C{float}, also
882 if B{C{ndigits}} is omitted.
883 '''
884 return LatLon2Tuple(round(self.lat, ndigits),
885 round(self.lon, ndigits), name=self.name)
887 @deprecated_method
888 def latlon_(self, ndigits=0): # PYCHOK no cover
889 '''DEPRECATED, use method L{latlon2}.'''
890 return self.latlon2(ndigits=ndigits)
892 latlon2round = latlon_ # PYCHOK no cover
894 @Property
895 def latlonheight(self):
896 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).
897 '''
898 return self.latlon.to3Tuple(self.height)
900 @latlonheight.setter # PYCHOK setter!
901 def latlonheight(self, latlonh):
902 '''Set the lat- and longitude and optionally the height
903 (2- or 3-tuple or comma- or space-separated C{str}
904 of C{degrees90}, C{degrees180} and C{meter}).
906 @see: Property L{latlon} for more details.
907 '''
908 self.latlon = latlonh
910 @Property
911 def lon(self):
912 '''Get the longitude (C{degrees180}).
913 '''
914 return self._lon
916 @lon.setter # PYCHOK setter!
917 def lon(self, lon):
918 '''Set the longitude (C{str[E|W]} or C{degrees}).
920 @raise ValueError: Invalid B{C{lon}}.
921 '''
922 lon = Lon(lon) # parseDMS(lon, suffix=_EW_, clip=180)
923 if self._lon != lon:
924 _update_all(self)
925 self._lon = lon
927 @Property_RO
928 def _ltp(self):
929 '''(INTERNAL) Cache for L{toLtp}.
930 '''
931 return _MODS.ltp.Ltp(self, ecef=self.Ecef(self.datum), name=self.name)
933 def nearestOn6(self, points, closed=False, height=None, wrap=False):
934 '''Locate the point on a path or polygon closest to this point.
936 Points are converted to and distances are computed in
937 I{geocentric}, cartesian space.
939 @arg points: The path or polygon points (C{LatLon}[]).
940 @kwarg closed: Optionally, close the polygon (C{bool}).
941 @kwarg height: Optional height, overriding the height of
942 this and all other points (C{meter}). If
943 C{None}, take the height of points into
944 account for distances.
945 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
946 the B{C{points}} (C{bool}).
948 @return: A L{NearestOn6Tuple}C{(closest, distance, fi, j,
949 start, end)} with the C{closest}, the C{start}
950 and the C{end} point each an instance of this
951 C{LatLon} and C{distance} in C{meter}, same
952 units as the cartesian axes.
954 @raise PointsError: Insufficient number of B{C{points}}.
956 @raise TypeError: Some B{C{points}} or some B{C{points}}'
957 C{Ecef} invalid.
959 @raise ValueError: Some B{C{points}}' C{Ecef} is incompatible.
961 @see: Function L{pygeodesy.nearestOn6}.
962 '''
963 def _cs(Ps, h, w, C):
964 p = None # not used
965 for i, q in Ps.enumerate():
966 if w and i:
967 q = _unrollon(p, q)
968 yield C(height=h, i=i, up=3, points=q)
969 p = q
971 C = self._toCartesianEcef # to verify datum and Ecef
972 Ps = self.PointsIter(points, wrap=wrap)
974 c = C(height=height, this=self) # this Cartesian
975 t = nearestOn6(c, _cs(Ps, height, wrap, C), closed=closed)
976 c, s, e = t.closest, t.start, t.end
978 kwds = _xkwds_not(None, LatLon=self.classof, # this LatLon
979 height=height)
980 r = self.Ecef(self.datum).reverse
981 p = r(c).toLatLon(**kwds)
982 s = r(s).toLatLon(**kwds) if s is not c else p
983 e = r(e).toLatLon(**kwds) if e is not c else p
984 return t.dup(closest=p, start=s, end=e)
986 def normal(self):
987 '''Normalize this point I{in-place} to C{abs(lat) <= 90} and
988 C{abs(lon) <= 180}.
990 @return: C{True} if this point was I{normal}, C{False} if it
991 wasn't (but is now).
993 @see: Property L{isnormal} and method L{toNormal}.
994 '''
995 n = self.isnormal
996 if not n:
997 self.latlon = normal(*self.latlon)
998 return n
1000 @Property_RO
1001 def _N_vector(self):
1002 '''(INTERNAL) Get the (C{nvectorBase._N_vector_})
1003 '''
1004 return _MODS.nvectorBase._N_vector_(*self.xyzh)
1006 @Property_RO
1007 def phi(self):
1008 '''Get the latitude (B{C{radians}}).
1009 '''
1010 return radians(self.lat)
1012 @Property_RO
1013 def philam(self):
1014 '''Get the lat- and longitude (L{PhiLam2Tuple}C{(phi, lam)}).
1015 '''
1016 return PhiLam2Tuple(self.phi, self.lam, name=self.name)
1018 def philam2(self, ndigits=0):
1019 '''Return this point's lat- and longitude in C{radians}, rounded.
1021 @kwarg ndigits: Number of (decimal) digits (C{int}).
1023 @return: A L{PhiLam2Tuple}C{(phi, lam)}, both C{float}
1024 and rounded away from zero.
1026 @note: The C{round}ed values are always C{float}, also
1027 if B{C{ndigits}} is omitted.
1028 '''
1029 return PhiLam2Tuple(round(self.phi, ndigits),
1030 round(self.lam, ndigits), name=self.name)
1032 @Property_RO
1033 def philamheight(self):
1034 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
1035 '''
1036 return self.philam.to3Tuple(self.height)
1038 @deprecated_method
1039 def points(self, points, closed=True): # PYCHOK no cover
1040 '''DEPRECATED, use method L{points2}.'''
1041 return self.points2(points, closed=closed)
1043 def points2(self, points, closed=True):
1044 '''Check a path or polygon represented by points.
1046 @arg points: The path or polygon points (C{LatLon}[])
1047 @kwarg closed: Optionally, consider the polygon closed,
1048 ignoring any duplicate or closing final
1049 B{C{points}} (C{bool}).
1051 @return: A L{Points2Tuple}C{(number, points)}, an C{int}
1052 and C{list} or C{tuple}.
1054 @raise PointsError: Insufficient number of B{C{points}}.
1056 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1057 '''
1058 return points2(points, closed=closed, base=self)
1060 def PointsIter(self, points, loop=0, dedup=False, wrap=False):
1061 '''Return a C{PointsIter} iterator.
1063 @arg points: The path or polygon points (C{LatLon}[])
1064 @kwarg loop: Number of loop-back points (non-negative C{int}).
1065 @kwarg dedup: Skip duplicate points (C{bool}).
1066 @kwarg wrap: If C{True}, wrap or I{normalize} the
1067 enum-/iterated B{C{points}} (C{bool}).
1069 @return: A new C{PointsIter} iterator.
1071 @raise PointsError: Insufficient number of B{C{points}}.
1072 '''
1073 return PointsIter(points, base=self, loop=loop, dedup=dedup, wrap=wrap)
1075 def radii11(self, point2, point3, wrap=False):
1076 '''Return the radii of the C{Circum-}, C{In-}, I{Soddy} and C{Tangent}
1077 circles of a (planar) triangle formed by this and two other points.
1079 @arg point2: Second point (C{LatLon}).
1080 @arg point3: Third point (C{LatLon}).
1081 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{point2}} and
1082 B{C{point3}} (C{bool}).
1084 @return: L{Radii11Tuple}C{(rA, rB, rC, cR, rIn, riS, roS, a, b, c, s)}.
1086 @raise IntersectionError: Near-coincident or -colinear points.
1088 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
1090 @see: Function L{pygeodesy.radii11}, U{Incircle
1091 <https://MathWorld.Wolfram.com/Incircle.html>}, U{Soddy Circles
1092 <https://MathWorld.Wolfram.com/SoddyCircles.html>} and U{Tangent
1093 Circles<https://MathWorld.Wolfram.com/TangentCircles.html>}.
1094 '''
1095 with _toCartesian3(self, point2, point3, wrap) as cs:
1096 return _radii11ABC(*cs, useZ=True)[0]
1098 def _rhumbx3(self, exact, radius): # != .sphericalBase._rhumbs3
1099 '''(INTERNAL) Get the C{rhumb} for this point's datum or for
1100 the earth model or earth B{C{radius}} if not C{None}.
1101 '''
1102 D = self.datum if radius is None else _spherical_datum(radius) # ellipsoidal OK
1103 x = _MODS.rhumbx # XXX Property_RO?
1104 r = D.ellipsoid.rhumbx if exact else \
1105 x.Rhumb(D, exact=False, name=D.name)
1106 return r, D, x.Caps
1108 def rhumbAzimuthTo(self, other, exact=False, radius=None, wrap=False):
1109 '''Return the azimuth (bearing) of a rhumb line (loxodrome)
1110 between this and an other (ellipsoidal) point.
1112 @arg other: The other point (C{LatLon}).
1113 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1114 default C{False}.
1115 @kwarg radius: Optional earth radius (C{meter}) or earth model
1116 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1117 L{a_f2Tuple}), overriding this point's datum.
1118 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1119 B{C{other}} point (C{bool}).
1121 @return: Rhumb azimuth (compass C{degrees360}).
1123 @raise TypeError: The B{C{other}} point is incompatible or
1124 B{C{radius}} is invalid.
1125 '''
1126 r, _, C = self._rhumbx3(exact, radius)
1127 return r._Inverse(self, other, wrap, outmask=C.AZIMUTH).azi12
1129 def rhumbDestination(self, distance, azimuth, exact=False, radius=None, height=None):
1130 '''Return the destination point having travelled the given distance
1131 from this point along a rhumb line (loxodrome) at the given azimuth.
1133 @arg distance: Distance travelled (C{meter}, same units as this
1134 point's datum (ellipsoid) axes or B{C{radius}},
1135 may be negative.
1136 @arg azimuth: Azimuth (bearing) at this point (compass C{degrees}).
1137 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1138 default C{False}.
1139 @kwarg radius: Optional earth radius (C{meter}) or earth model
1140 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1141 L{a_f2Tuple}), overriding this point's datum.
1142 @kwarg height: Optional height, overriding the default height
1143 (C{meter}).
1145 @return: The destination point (ellipsoidal C{LatLon}).
1147 @raise TypeError: Invalid B{C{radius}}.
1149 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}},
1150 B{C{radius}} or B{C{height}}.
1151 '''
1152 r, D, _ = self._rhumbx3(exact, radius)
1153 d = r._Direct(self, azimuth, distance)
1154 h = self._heigHt(height)
1155 return self.classof(d.lat2, d.lon2, datum=D, height=h)
1157 def rhumbDistanceTo(self, other, exact=False, radius=None, wrap=False):
1158 '''Return the distance from this to an other point along
1159 a rhumb line (loxodrome).
1161 @arg other: The other point (C{LatLon}).
1162 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1163 default C{False}.
1164 @kwarg radius: Optional earth radius (C{meter}) or earth model
1165 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1166 L{a_f2Tuple}), overriding this point's datum.
1167 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1168 B{C{other}} point (C{bool}).
1170 @return: Distance (C{meter}, the same units as this point's
1171 datum (ellipsoid) axes or B{C{radius}}.
1173 @raise TypeError: The B{C{other}} point is incompatible or
1174 B{C{radius}} is invalid.
1176 @raise ValueError: Invalid B{C{radius}}.
1177 '''
1178 r, _, C = self._rhumbx3(exact, radius)
1179 return r._Inverse(self, other, wrap, outmask=C.DISTANCE).s12
1181 def rhumbLine(self, azimuth_other, exact=False, radius=None, wrap=False,
1182 **name_caps):
1183 '''Get a rhumb line through this point at a given azimuth or
1184 through this and an other point.
1186 @arg azimuth_other: The azimuth of the rhumb line (compass
1187 C{degrees}) or the other point (C{LatLon}).
1188 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1189 default C{False}.
1190 @kwarg radius: Optional earth radius (C{meter}) or earth model
1191 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1192 L{a_f2Tuple}), overriding this point's datum.
1193 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1194 C{azimuth_B{other}} point (C{bool}).
1195 @kwarg name_caps: Optional C{B{name}=str} and C{caps}, see
1196 L{RhumbLine} C{B{caps}}.
1198 @return: A L{RhumbLine} instance.
1200 @raise TypeError: Invalid B{C{radius}} or BC{C{azimuth_other}}
1201 not a C{scalar} nor a C{LatLon}.
1203 @see: Classes L{RhumbLine} and L{Rhumb}, property L{Rhumb.exact}
1204 and methods L{Rhumb.DirectLine} and L{Rhumb.InverseLine}.
1205 '''
1206 r, _, _ = self._rhumbx3(exact, radius)
1207 a, kwds = azimuth_other, _xkwds(name_caps, name=self.name)
1208 if isscalar(a):
1209 r = r._DirectLine(self, a, **kwds)
1210 elif isinstance(a, LatLonBase):
1211 r = r._InverseLine(self, a, wrap, **kwds)
1212 else:
1213 raise _TypeError(azimuth_other=a)
1214 return r
1216 def rhumbMidpointTo(self, other, exact=False, radius=None,
1217 height=None, fraction=_0_5, wrap=False):
1218 '''Return the (loxodromic) midpoint on the rhumb line between
1219 this and an other point.
1221 @arg other: The other point (C{LatLon}).
1222 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),
1223 default C{False}.
1224 @kwarg radius: Optional earth radius (C{meter}) or earth model
1225 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1226 L{a_f2Tuple}), overriding this point's datum.
1227 @kwarg height: Optional height, overriding the mean height
1228 (C{meter}).
1229 @kwarg fraction: Midpoint location from this point (C{scalar}),
1230 may be negative or greater than 1.0.
1231 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1232 B{C{other}} point (C{bool}).
1234 @return: The midpoint at the given B{C{fraction}} along the
1235 rhumb line (C{LatLon}).
1237 @raise TypeError: The B{C{other}} point is incompatible or
1238 B{C{radius}} is invalid.
1240 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}.
1241 '''
1242 r, D, _ = self._rhumbx3(exact, radius)
1243 f = Scalar(fraction=fraction)
1244 d = r._Inverse(self, other, wrap) # C.AZIMUTH_DISTANCE
1245 d = r._Direct( self, d.azi12, d.s12 * f)
1246 h = self._havg(other, f=f, h=height)
1247 return self.classof(d.lat2, d.lon2, datum=D, height=h)
1249 def thomasTo(self, other, wrap=False):
1250 '''Compute the distance between this and an other point using
1251 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}
1252 formula.
1254 @arg other: The other point (C{LatLon}).
1255 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1256 the B{C{other}} point (C{bool}).
1258 @return: Distance (C{meter}, same units as the axes of
1259 this point's datum ellipsoid).
1261 @raise TypeError: The B{C{other}} point is not C{LatLon}.
1263 @see: Function L{pygeodesy.thomas} and methods L{cosineAndoyerLambertTo},
1264 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
1265 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
1266 L{flatPolarTo}, L{haversineTo} and L{vincentysTo}.
1267 '''
1268 return self._distanceTo_(thomas_, other, wrap=wrap)
1270 @deprecated_method
1271 def to2ab(self): # PYCHOK no cover
1272 '''DEPRECATED, use property L{philam}.'''
1273 return self.philam
1275 def toCartesian(self, height=None, Cartesian=None, **Cartesian_kwds):
1276 '''Convert this point to cartesian, I{geocentric} coordinates,
1277 also known as I{Earth-Centered, Earth-Fixed} (ECEF).
1279 @kwarg height: Optional height, overriding this point's height
1280 (C{meter}, conventionally).
1281 @kwarg Cartesian: Optional class to return the geocentric
1282 coordinates (C{Cartesian}) or C{None}.
1283 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}}
1284 keyword arguments, ignored if
1285 C{B{Cartesian} is None}.
1287 @return: A B{C{Cartesian}} or if B{C{Cartesian}} is C{None},
1288 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
1289 datum)} with C{C=0} and C{M} if available.
1291 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}.
1292 '''
1293 r = self._ecef9 if height is None else self.toEcef(height=height)
1294 if Cartesian is not None: # class or .classof
1295 r = self._xnamed(Cartesian(r, **Cartesian_kwds))
1296 _xdatum(r.datum, self.datum)
1297 return r
1299 def _toCartesianEcef(self, height=None, i=None, up=2, **name_point):
1300 '''(INTERNAL) Convert to cartesian and check Ecef's before and after.
1301 '''
1302 p = self.others(up=up, **name_point)
1303 c = p.toCartesian(height=height)
1304 E = self.Ecef
1305 if E:
1306 for p in (p, c):
1307 e = getattr(p, LatLonBase.Ecef.name, None)
1308 if e not in (None, E): # PYCHOK no cover
1309 n, _ = name_point.popitem()
1310 if i is not None:
1311 Fmt.SQUARE(n, i)
1312 raise _ValueError(n, e, txt=_incompatible(E.__name__))
1313 return c
1315 def toEcef(self, height=None, M=False):
1316 '''Convert this point to I{geocentric} coordinates, also known as
1317 I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).
1319 @kwarg height: Optional height, overriding this point's height
1320 (C{meter}, conventionally).
1321 @kwarg M: Optionally, include the rotation L{EcefMatrix} (C{bool}).
1323 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
1324 with C{C=0} and C{M} if available.
1326 @raise EcefError: A C{.datum} or an ECEF issue.
1327 '''
1328 return self._ecef9 if height in (None, self.height) else \
1329 self._Ecef_forward(self.lat, self.lon, height=height, M=M)
1331 @deprecated_method
1332 def to3llh(self, height=None): # PYCHOK no cover
1333 '''DEPRECATED, use property L{latlonheight} or C{latlon.to3Tuple(B{height})}.'''
1334 return self.latlonheight if height in (None, self.height) else \
1335 self.latlon.to3Tuple(height)
1337 def toLocal(self, Xyz=None, ltp=None, **Xyz_kwds):
1338 '''Convert this I{geodetic} point to I{local} C{X}, C{Y} and C{Z}.
1340 @kwarg Xyz: Optional class to return C{X}, C{Y} and C{Z}
1341 (L{XyzLocal}, L{Enu}, L{Ned}) or C{None}.
1342 @kwarg ltp: The I{local tangent plane} (LTP) to use,
1343 overriding this point's LTP (L{Ltp}).
1344 @kwarg Xyz_kwds: Optional, additional B{C{Xyz}} keyword
1345 arguments, ignored if C{B{Xyz} is None}.
1347 @return: An B{C{Xyz}} instance or if C{B{Xyz} is None},
1348 a L{Local9Tuple}C{(x, y, z, lat, lon, height,
1349 ltp, ecef, M)} with C{M=None}, always.
1351 @raise TypeError: Invalid B{C{ltp}}.
1352 '''
1353 p = _MODS.ltp._xLtp(ltp, self._ltp)
1354 return p._ecef2local(self._ecef9, Xyz, Xyz_kwds)
1356 def toLtp(self, Ecef=None):
1357 '''Return the I{local tangent plane} (LTP) for this point.
1359 @kwarg Ecef: Optional ECEF I{class} (L{EcefKarney}, ...
1360 L{EcefYou}), overriding this point's C{Ecef}.
1361 '''
1362 return self._ltp if Ecef in (None, self.Ecef) else _MODS.ltp.Ltp(
1363 self, ecef=Ecef(self.datum), name=self.name)
1365 def toNormal(self, deep=False, name=NN):
1366 '''Get this point I{normalized} to C{abs(lat) <= 90}
1367 and C{abs(lon) <= 180}.
1369 @kwarg deep: If C{True} make a deep, otherwise a
1370 shallow copy (C{bool}).
1371 @kwarg name: Optional name of the copy (C{str}).
1373 @return: A copy of this point, I{normalized} and
1374 optionally renamed (C{LatLon}).
1376 @see: Property L{isnormal}, method L{normal} and function
1377 L{pygeodesy.normal}.
1378 '''
1379 ll = self.copy(deep=deep)
1380 _ = ll.normal()
1381 if name:
1382 ll.rename(name)
1383 return ll
1385 def toNvector(self, h=None, Nvector=None, **Nvector_kwds):
1386 '''Convert this point to C{n-vector} (normal to the earth's surface)
1387 components, I{including height}.
1389 @kwarg h: Optional height, overriding this point's
1390 height (C{meter}).
1391 @kwarg Nvector: Optional class to return the C{n-vector}
1392 components (C{Nvector}) or C{None}.
1393 @kwarg Nvector_kwds_wrap: Optional, additional B{C{Nvector}}
1394 keyword arguments, ignored if C{B{Nvector}
1395 is None}.
1397 @return: A B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}
1398 if B{C{Nvector}} is C{None}.
1400 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.
1401 '''
1402 return self.toVector(Vector=Nvector, h=self.height if h is None else h,
1403 ll=self, **Nvector_kwds)
1405 def toStr(self, form=F_DMS, joined=_COMMASPACE_, m=_m_, **prec_sep_s_D_M_S): # PYCHOK expected
1406 '''Convert this point to a "lat, lon[, +/-height]" string, formatted
1407 in the given C{B{form}at}.
1409 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see
1410 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.
1411 @kwarg joined: Separator to join the lat-, longitude and heigth
1412 strings (C{str} or C{None} or C{NN} for non-joined).
1413 @kwarg m: Optional unit of the height (C{str}), use C{None} to
1414 exclude height from the returned string.
1415 @kwarg prec_sep_s_D_M_S: Optional C{B{prec}ision}, C{B{sep}arator},
1416 B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}} keyword
1417 arguments, see function L{pygeodesy.latDMS} or
1418 L{pygeodesy.lonDMS}.
1420 @return: This point in the specified C{B{form}at}, etc. (C{str} or
1421 a 2- or 3-tuple C{(lat_str, lon_str[, height_str])} if
1422 C{B{joined}=NN} or C{B{joined}=None}).
1424 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more
1425 details about keyword arguments C{B{form}at}, C{B{prec}ision},
1426 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.
1428 @example:
1430 >>> LatLon(51.4778, -0.0016).toStr() # 51°28′40″N, 000°00′06″W
1431 >>> LatLon(51.4778, -0.0016).toStr(F_D) # 51.4778°N, 000.0016°W
1432 >>> LatLon(51.4778, -0.0016, 42).toStr() # 51°28′40″N, 000°00′06″W, +42.00m
1433 '''
1434 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),
1435 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))
1436 if self.height and m is not None:
1437 t += (self.heightStr(m=m),)
1438 return joined.join(t) if joined else t
1440 def toVector(self, Vector=None, **Vector_kwds):
1441 '''Convert this point to C{n-vector} (normal to the earth's
1442 surface) components, I{ignoring height}.
1444 @kwarg Vector: Optional class to return the C{n-vector}
1445 components (L{Vector3d}) or C{None}.
1446 @kwarg Vector_kwds: Optional, additional B{C{Vector}}
1447 keyword arguments, ignored if
1448 C{B{Vector} is None}.
1450 @return: A B{C{Vector}} or a L{Vector3Tuple}C{(x, y, z)}
1451 if B{C{Vector}} is C{None}.
1453 @raise TypeError: Invalid B{C{Vector}} or B{C{kwds}}.
1455 @note: These are C{n-vector} x, y and z components,
1456 I{NOT} geocentric (ECEF) x, y and z coordinates!
1457 '''
1458 r = self._vector3tuple
1459 if Vector is not None:
1460 r = Vector(*r, **_xkwds(Vector_kwds, name=self.name))
1461 return r
1463 def toVector3d(self):
1464 '''Convert this point to C{n-vector} (normal to the earth's
1465 surface) components, I{ignoring height}.
1467 @return: Unit vector (L{Vector3d}).
1469 @note: These are C{n-vector} x, y and z components,
1470 I{NOT} geocentric (ECEF) x, y and z coordinates!
1471 '''
1472 return self._vector3d # XXX .unit()
1474 @deprecated_method
1475 def to3xyz(self): # PYCHOK no cover
1476 '''DEPRECATED, use property L{xyz} or method L{toNvector}, L{toVector},
1477 L{toVector3d} or perhaps (geocentric) L{toEcef}.'''
1478 return self.xyz # self.toVector()
1480 @Property_RO
1481 def _vector3d(self):
1482 '''(INTERNAL) Cache for L{toVector3d}.
1483 '''
1484 return self.toVector(Vector=Vector3d) # XXX .unit()
1486 @Property_RO
1487 def _vector3tuple(self):
1488 '''(INTERNAL) Cache for L{toVector}.
1489 '''
1490 return philam2n_xyz(self.phi, self.lam, name=self.name)
1492 def vincentysTo(self, other, **radius_wrap):
1493 '''Compute the distance between this and an other point using
1494 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}
1495 spherical formula.
1497 @arg other: The other point (C{LatLon}).
1498 @kwarg radius_wrap: Optional keyword arguments for function
1499 L{pygeodesy.vincentys}, overriding the
1500 default mean C{radius} of this point's
1501 datum ellipsoid.
1503 @return: Distance (C{meter}, same units as B{C{radius}}).
1505 @raise TypeError: The B{C{other}} point is not C{LatLon}.
1507 @see: Function L{pygeodesy.vincentys} and methods L{cosineAndoyerLambertTo},
1508 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},
1509 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},
1510 L{flatPolarTo}, L{haversineTo} and L{thomasTo}.
1511 '''
1512 return self._distanceTo(vincentys, other, **_xkwds(radius_wrap, radius=None))
1514 @Property_RO
1515 def xyz(self):
1516 '''Get the C{n-vector} X, Y and Z components (L{Vector3Tuple}C{(x, y, z)})
1518 @note: These are C{n-vector} x, y and z components, I{NOT}
1519 geocentric (ECEF) x, y and z coordinates!
1520 '''
1521 return self.toVector(Vector=Vector3Tuple)
1523 @Property_RO
1524 def xyzh(self):
1525 '''Get the C{n-vector} X, Y, Z and H components (L{Vector4Tuple}C{(x, y, z, h)})
1527 @note: These are C{n-vector} x, y and z components, I{NOT}
1528 geocentric (ECEF) x, y and z coordinates!
1529 '''
1530 return self.xyz.to4Tuple(self.height)
1533class _toCartesian3(object): # see also .geodesicw._wargs, .vector2d._numpy
1534 '''(INTERNAL) Wrapper convert 2 other points.
1535 '''
1536 @contextmanager # <https://www.python.org/dev/peps/pep-0343/> Examples
1537 def __call__(self, p, p2, p3, wrap, **kwds):
1538 try:
1539 if wrap:
1540 p2, p3 = map1(_Wrap.point, p2, p3)
1541 kwds = _xkwds(kwds, wrap=wrap)
1542 yield (p. toCartesian().copy(name=_point_), # copy to rename
1543 p._toCartesianEcef(up=4, point2=p2),
1544 p._toCartesianEcef(up=4, point3=p3))
1545 except (AssertionError, TypeError, ValueError) as x:
1546 raise _xError(x, point=p, point2=p2, point3=p3, **kwds)
1548_toCartesian3 = _toCartesian3() # PYCHOK singleton
1551def _trilaterate5(p1, d1, p2, d2, p3, d3, area=True, eps=EPS1, # MCCABE 13
1552 radius=R_M, wrap=False):
1553 '''(INTERNAL) Trilaterate three points by area overlap or by
1554 perimeter intersection of three circles.
1556 @note: The B{C{radius}} is only needed for both the n-vectorial
1557 and C{sphericalTrigonometry.LatLon.distanceTo} methods and
1558 silently ignored by the C{ellipsoidalExact}, C{-GeodSolve},
1559 C{-Karney} and C{-Vincenty.LatLon.distanceTo} methods.
1560 '''
1561 p2, p3, w = _unrollon3(p1, p2, p3, wrap)
1563 r1 = Distance_(distance1=d1)
1564 r2 = Distance_(distance2=d2)
1565 r3 = Distance_(distance3=d3)
1566 m = 0 if area else (r1 + r2 + r3)
1567 pc = 0
1568 t = []
1569 for _ in range(3):
1570 try: # intersection of circle (p1, r1) and (p2, r2)
1571 c1, c2 = p1.intersections2(r1, p2, r2, wrap=w)
1573 if area: # check overlap
1574 if c1 is c2: # abutting
1575 c = c1
1576 else: # nearest point on radical
1577 c = p3.nearestOn(c1, c2, within=True, wrap=w)
1578 d = r3 - p3.distanceTo(c, radius=radius, wrap=w)
1579 if d > eps: # sufficient overlap
1580 t.append((d, c))
1581 m = max(m, d)
1583 else: # check intersection
1584 for c in ((c1,) if c1 is c2 else (c1, c2)):
1585 d = fabs(r3 - p3.distanceTo(c, radius=radius, wrap=w))
1586 if d < eps: # below margin
1587 t.append((d, c))
1588 m = min(m, d)
1590 except IntersectionError as x:
1591 if _concentric_ in str(x): # XXX ConcentricError?
1592 pc += 1
1594 p1, r1, p2, r2, p3, r3 = p2, r2, p3, r3, p1, r1 # rotate
1596 if t: # get min, max, points and count ...
1597 t = tuple(sorted(t))
1598 n = (len(t),) # as 1-tuple
1599 # ... or for a single trilaterated result,
1600 # min *is* max, min- *is* maxPoint and n=1
1601 return Trilaterate5Tuple(t[0] + t[-1] + n) # *(t[0] + ...)
1603 elif area and pc == 3: # all pairwise concentric ...
1604 r, p = min((r1, p1), (r2, p2), (r3, p3))
1605 m = max(r1, r2, r3)
1606 # ... return "smallest" point twice, the smallest
1607 # and largest distance and n=0 for concentric
1608 return Trilaterate5Tuple(float(r), p, float(m), p, 0)
1610 n, f = (_overlap_, max) if area else (_intersection_, min)
1611 t = '%s (%s %.3f)' % (_no_(n), f.__name__, m)
1612 raise IntersectionError(area=area, eps=eps, wrap=wrap, txt=t)
1615__all__ += _ALL_DOCS(LatLonBase)
1617# **) MIT License
1618#
1619# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1620#
1621# Permission is hereby granted, free of charge, to any person obtaining a
1622# copy of this software and associated documentation files (the "Software"),
1623# to deal in the Software without restriction, including without limitation
1624# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1625# and/or sell copies of the Software, and to permit persons to whom the
1626# Software is furnished to do so, subject to the following conditions:
1627#
1628# The above copyright notice and this permission notice shall be included
1629# in all copies or substantial portions of the Software.
1630#
1631# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1632# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1633# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1634# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1635# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1636# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1637# OTHER DEALINGS IN THE SOFTWARE.