Coverage for pygeodesy/points.py: 93%
533 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-11 16:04 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-11 16:04 -0400
2# -*- coding: utf-8 -*-
4u'''Utilities for point lists, tuples, etc.
6Functions to handle collections and sequences of C{LatLon} points
7specified as 2-d U{NumPy<https://www.NumPy.org>}, C{arrays} or tuples
8as C{LatLon} or as C{pseudo-x/-y} pairs.
10C{NumPy} arrays are assumed to contain rows of points with a lat-, a
11longitude -and possibly other- values in different columns. While
12iterating over the array rows, create an instance of a given C{LatLon}
13class "on-the-fly" for each row with the row's lat- and longitude.
15The original C{NumPy} array is read-accessed only and never duplicated,
16except to return a I{subset} of the original array.
18For example, to process a C{NumPy} array, wrap the array by instantiating
19class L{Numpy2LatLon} and specifying the column index for the lat- and
20longitude in each row. Then, pass the L{Numpy2LatLon} instance to any
21L{pygeodesy} function or method accepting a I{points} argument.
23Similarly, class L{Tuple2LatLon} is used to instantiate a C{LatLon} from
24each 2+tuple in a sequence of such 2+tuples using the C{ilat} lat- and
25C{ilon} longitude index in each 2+tuple.
26'''
28from pygeodesy.basics import isclass, isint, isscalar, issequence, \
29 issubclassof, _Sequence, _xcopy, _xdup, \
30 _xinstanceof
31from pygeodesy.constants import EPS, EPS1, PI_2, R_M, isnear0, isnear1, \
32 _umod_360, _0_0, _0_5, _1_0, _2_0, _6_0, \
33 _90_0, _N_90_0, _180_0, _360_0
34# from pygeodesy.datums import _spherical_datum # from .formy
35from pygeodesy.dms import F_D, parseDMS
36from pygeodesy.errors import CrossError, crosserrors, _IndexError, \
37 _IsnotError, _TypeError, _ValueError, \
38 _xattr, _xkwds, _xkwds_pop
39from pygeodesy.fmath import favg, fdot, hypot, Fsum, fsum
40# from pygeodesy.fsums import Fsum, fsum # from .fmath
41from pygeodesy.formy import _bearingTo2, equirectangular_, _spherical_datum
42from pygeodesy.interns import NN, _colinear_, _COMMASPACE_, _composite_, \
43 _DEQUALSPACED_, _ELLIPSIS_, _EW_, _immutable_, \
44 _near_, _no_, _not_, _NS_, _point_, _SPACE_, \
45 _UNDER_, _valid_ # _lat_, _lon_
46from pygeodesy.iters import LatLon2PsxyIter, PointsIter, points2
47from pygeodesy.latlonBase import LatLonBase, _latlonheight3, \
48 _ALL_DOCS, _ALL_LAZY, _MODS
49# from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
50from pygeodesy.named import classname, nameof, notImplemented, notOverloaded, \
51 _NamedTuple
52from pygeodesy.namedTuples import Bounds2Tuple, Bounds4Tuple, LatLon2Tuple, \
53 NearestOn3Tuple, NearestOn5Tuple, \
54 Point3Tuple, Vector3Tuple, \
55 PhiLam2Tuple # PYCHOK shared
56from pygeodesy.props import Property_RO, property_doc_, property_RO
57from pygeodesy.streprs import Fmt, instr
58from pygeodesy.units import Number_, Radius, Scalar, Scalar_
59from pygeodesy.utily import atan2b, degrees90, degrees180, degrees2m, \
60 unroll180, _unrollon, unrollPI, _Wrap, wrap180
62from math import cos, fabs, fmod, radians, sin
64__all__ = _ALL_LAZY.points
65__version__ = '23.10.07'
67_ilat_ = 'ilat'
68_ilon_ = 'ilon'
69_ncols_ = 'ncols'
70_nrows_ = 'nrows'
73class LatLon_(LatLonBase): # XXX in heights._HeightBase.height
74 '''Low-overhead C{LatLon} class, mainly for L{Numpy2LatLon} and L{Tuple2LatLon}.
75 '''
76 # __slots__ efficiency is voided if the __slots__ class attribute is
77 # used in a subclass of a class with the traditional __dict__, @see
78 # <https://docs.Python.org/2/reference/datamodel.html#slots> plus ...
79 #
80 # __slots__ must be repeated in sub-classes, @see Luciano Ramalho,
81 # "Fluent Python", O'Reilly, 2016 p. 276+ "Problems with __slots__",
82 # 2nd Ed, 2022 p. 390 "Summarizing the Issues with __slots__".
83 #
84 # __slots__ = (_lat_, _lon_, _height_, _datum_, _name_)
85 # Property_RO = property_RO # no __dict__ with __slots__!
86 #
87 # In addition, both size and overhead have shrunk in recent Python:
88 #
89 # sys.getsizeof(LatLon_(1, 2)) is 72-88 I{with} __slots__, but
90 # only 48-56 bytes I{without in Python 2.7.18+ and Python 3+}.
91 #
92 # python3 -m timeit -s "from pygeodesy... import LatLonBase as LL" "LL(0, 0)" 2.14 usec
93 # python3 -m timeit -s "from pygeodesy import LatLon_" "LatLon_(0, 0)" 216 nsec
95 def __init__(self, latlonh, lon=None, height=0, wrap=False, name=NN, datum=None):
96 '''New L{LatLon_}.
98 @note: The lat- and longitude values are taken I{as-given,
99 un-clipped and un-validated}.
101 @see: L{latlonBase.LatLonBase} for further details.
102 '''
103 if name:
104 self.name = name
106 if lon is None: # PYCHOK no cover
107 lat, lon, height = _latlonheight3(latlonh, height, wrap)
108 elif wrap: # PYCHOK no cover
109 lat, lon = _Wrap.latlonDMS2(latlonh, lon)
110 else: # must be latNS, lonEW
111 try:
112 lat, lon = float(latlonh), float(lon)
113 except (TypeError, ValueError):
114 lat = parseDMS(latlonh, suffix=_NS_)
115 lon = parseDMS(lon, suffix=_EW_)
117 # get the minimal __dict__, see _isLatLon_ below
118 self._lat = lat # un-clipped and ...
119 self._lon = lon # ... un-validated
120 self._datum = None if datum is None else \
121 _spherical_datum(datum, name=self.name)
122 self._height = height
124 def __eq__(self, other):
125 return isinstance(other, LatLon_) and \
126 other.lat == self.lat and \
127 other.lon == self.lon
129 def __ne__(self, other):
130 return not self.__eq__(other)
132 @Property_RO
133 def datum(self):
134 '''Get the C{datum} (L{Datum}) or C{None}.
135 '''
136 return self._datum
138 def intermediateTo(self, other, fraction, height=None, wrap=False):
139 '''Locate the point at a given fraction, I{linearly} between
140 (or along) this and an other point.
142 @arg other: The other point (C{LatLon}).
143 @arg fraction: Fraction between both points (C{float},
144 0.0 for this and 1.0 for the other point).
145 @kwarg height: Optional height (C{meter}), overriding the
146 intermediate height.
147 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
148 the B{C{other}} point (C{bool}).
150 @return: Intermediate point (this C{LatLon}).
152 @raise TypeError: Incompatible B{C{other}} C{type}.
153 '''
154 f = Scalar(fraction=fraction)
155 if isnear0(f):
156 r = self
157 else:
158 r = self.others(other)
159 if wrap or not isnear1(f):
160 _, lat, lon = _Wrap.latlon3(self.lon, r.lat, r.lon, wrap)
161 lat = favg(self.lat, lat, f=f)
162 lon = favg(self.lon, lon, f=f)
163 h = height if height is not None else \
164 favg(self.height, r.height, f=f)
165 # = self._havg(r, f=f, h=height)
166 r = self.classof(lat, lon, height=h, datum=r.datum,
167 name=r.intermediateTo.__name__)
168 return r
170 def toRepr(self, **kwds):
171 '''This L{LatLon_} as a string "class(<degrees>, ...)",
172 ignoring keyword argument C{B{std}=N/A}.
174 @see: L{latlonBase.LatLonBase.toRepr} for further details.
175 '''
176 _ = _xkwds_pop(kwds, std=NotImplemented)
177 return LatLonBase.toRepr(self, **kwds)
179 def toStr(self, form=F_D, joined=_COMMASPACE_, **m_prec_sep_s_D_M_S): # PYCHOK expected
180 '''Convert this point to a "lat, lon[, height][, name][, ...]"
181 string, formatted in the given C{B{form}at}.
183 @see: L{latlonBase.LatLonBase.toStr} for further details.
184 '''
185 t = LatLonBase.toStr(self, form=form, joined=NN,
186 **_xkwds(m_prec_sep_s_D_M_S, m=NN))
187 if self.name:
188 t += (repr(self.name),)
189 return joined.join(t) if joined else t
192def _isLatLon(inst):
193 '''(INTERNAL) Check a C{LatLon} or C{LatLon_} instance.
194 '''
195 return isinstance(inst, (LatLon_, _MODS.latlonBase.LatLonBase))
198def _isLatLon_(LL):
199 '''(INTERNAL) Check a (sub-)class of C{LatLon_}.
200 '''
201 return issubclassof(LL, LatLon_) or (isclass(LL) and
202 all(hasattr(LL, _) for _ in LatLon_(0, 0).__dict__.keys()))
205class _Basequence(_Sequence): # immutable, on purpose
206 '''(INTERNAL) Base class.
207 '''
208 _array = []
209 _epsilon = EPS
210 _itemname = _point_
212 def _contains(self, point):
213 '''(INTERNAL) Check for a matching point.
214 '''
215 return any(self._findall(point, ()))
217 def copy(self, deep=False): # PYCHOK no cover
218 '''Make a shallow or deep copy of this instance.
220 @kwarg deep: If C{True} make a deep, otherwise a
221 shallow copy (C{bool}).
223 @return: The copy (C{This class} or subclass thereof).
224 '''
225 return _xcopy(self, deep=deep)
227 def _count(self, point):
228 '''(INTERNAL) Count the number of matching points.
229 '''
230 return sum(1 for _ in self._findall(point, ())) # NOT len()!
232 def dup(self, **items): # PYCHOK no cover
233 '''Duplicate this instance, I{without replacing items}.
235 @kwarg items: No attributes (I{not allowed}).
237 @return: The duplicate (C{This} (sub-)class).
239 @raise TypeError: Any B{C{items}} invalid.
240 '''
241 if items:
242 t = _SPACE_(classname(self), _immutable_)
243 raise _TypeError(txt=t, this=self, **items)
244 return _xdup(self)
246 @property_doc_(''' the equality tolerance (C{float}).''')
247 def epsilon(self):
248 '''Get the tolerance for equality tests (C{float}).
249 '''
250 return self._epsilon
252 @epsilon.setter # PYCHOK setter!
253 def epsilon(self, tol):
254 '''Set the tolerance for equality tests (C{scalar}).
256 @raise UnitError: Non-scalar or invalid B{C{tol}}.
257 '''
258 self._epsilon = Scalar_(tolerance=tol)
260 def _find(self, point, start_end):
261 '''(INTERNAL) Find the first matching point index.
262 '''
263 for i in self._findall(point, start_end):
264 return i
265 return -1
267 def _findall(self, point, start_end): # PYCHOK no cover
268 '''(INTERNAL) I{Must be implemented/overloaded}.
269 '''
270 notImplemented(self, point, start_end)
272 def _getitem(self, index):
273 '''(INTERNAL) Return point [index] or return a slice.
274 '''
275 # Luciano Ramalho, "Fluent Python", O'Reilly, 2016 p. 290+, 2022 p. 405+
276 if isinstance(index, slice):
277 # XXX an numpy.[nd]array slice is a view, not a copy
278 return self.__class__(self._array[index], **self._slicekwds())
279 else:
280 return self.point(self._array[index])
282 def _index(self, point, start_end):
283 '''(INTERNAL) Find the first matching point index.
284 '''
285 for i in self._findall(point, start_end):
286 return i
287 raise _IndexError(self._itemname, point, txt=_not_('found'))
289 @property_RO
290 def isNumpy2(self): # PYCHOK no cover
291 '''Is this a Numpy2 wrapper?
292 '''
293 return False # isinstance(self, (Numpy2LatLon, ...))
295 @property_RO
296 def isPoints2(self): # PYCHOK no cover
297 '''Is this a LatLon2 wrapper/converter?
298 '''
299 return False # isinstance(self, (LatLon2psxy, ...))
301 @property_RO
302 def isTuple2(self): # PYCHOK no cover
303 '''Is this a Tuple2 wrapper?
304 '''
305 return False # isinstance(self, (Tuple2LatLon, ...))
307 def _iter(self):
308 '''(INTERNAL) Yield all points.
309 '''
310 _array, _point = self._array, self.point
311 for i in range(len(self)):
312 yield _point(_array[i])
314 def point(self, *attrs): # PYCHOK no cover
315 '''I{Must be overloaded}.
317 @arg attrs: Optional arguments.
318 '''
319 notOverloaded(self, *attrs)
321 def _range(self, start=None, end=None, step=1):
322 '''(INTERNAL) Return the range.
323 '''
324 if step > 0:
325 if start is None:
326 start = 0
327 if end is None:
328 end = len(self)
329 elif step < 0:
330 if start is None:
331 start = len(self) - 1
332 if end is None:
333 end = -1
334 else:
335 raise _ValueError(step=step)
336 return range(start, end, step)
338 def _repr(self):
339 '''(INTERNAL) Return a string representation.
340 '''
341 # XXX use Python 3+ reprlib.repr
342 t = repr(self._array[:1]) # only first row
343 t = _SPACE_(t[:-1], _ELLIPSIS_, Fmt.SQUARE(t[-1:], len(self)))
344 t = _SPACE_.join(t.split()) # coalesce spaces
345 return instr(self, t, **self._slicekwds())
347 def _reversed(self): # PYCHOK false
348 '''(INTERNAL) Yield all points in reverse order.
349 '''
350 _array, point = self._array, self.point
351 for i in range(len(self) - 1, -1, -1):
352 yield point(_array[i])
354 def _rfind(self, point, start_end):
355 '''(INTERNAL) Find the last matching point index.
356 '''
357 def _r3(start=None, end=None, step=-1):
358 return (start, end, step) # PYCHOK returns
360 for i in self._findall(point, _r3(*start_end)):
361 return i
362 return -1
364 def _slicekwds(self): # PYCHOK no cover
365 '''(INTERNAL) I{Should be overloaded}.
366 '''
367 return {}
370class _Array2LatLon(_Basequence): # immutable, on purpose
371 '''(INTERNAL) Base class for Numpy2LatLon or Tuple2LatLon.
372 '''
373 _array = ()
374 _ilat = 0 # row column index
375 _ilon = 0 # row column index
376 _LatLon = LatLon_ # default
377 _shape = ()
379 def __init__(self, array, ilat=0, ilon=1, LatLon=None, shape=()):
380 '''Handle a C{NumPy} or C{Tuple} array as a sequence of C{LatLon} points.
381 '''
382 ais = (_ilat_, ilat), (_ilon_, ilon)
384 if len(shape) != 2 or shape[0] < 1 or shape[1] < len(ais):
385 raise _IndexError('array.shape', shape)
387 self._array = array
388 self._shape = Shape2Tuple(shape) # *shape
390 if LatLon: # check the point class
391 if not _isLatLon_(LatLon):
392 raise _IsnotError(_valid_, LatLon=LatLon)
393 self._LatLon = LatLon
395 # check the attr indices
396 for n, (ai, i) in enumerate(ais):
397 if not isint(i):
398 raise _IsnotError(int.__name__, **{ai: i})
399 i = int(i)
400 if not 0 <= i < shape[1]:
401 raise _ValueError(ai, i)
402 for aj, j in ais[:n]:
403 if int(j) == i:
404 raise _ValueError(_DEQUALSPACED_(ai, aj, i))
405 setattr(self, NN(_UNDER_, ai), i)
407 def __contains__(self, latlon):
408 '''Check for a specific lat-/longitude.
410 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
411 C{(lat, lon)}).
413 @return: C{True} if B{C{latlon}} is present, C{False} otherwise.
415 @raise TypeError: Invalid B{C{latlon}}.
416 '''
417 return self._contains(latlon)
419 def __getitem__(self, index):
420 '''Return row[index] as C{LatLon} or return a L{Numpy2LatLon} slice.
421 '''
422 return self._getitem(index)
424 def __iter__(self):
425 '''Yield rows as C{LatLon}.
426 '''
427 return self._iter()
429 def __len__(self):
430 '''Return the number of rows.
431 '''
432 return self._shape[0]
434 def __repr__(self):
435 '''Return a string representation.
436 '''
437 return self._repr()
439 def __reversed__(self): # PYCHOK false
440 '''Yield rows as C{LatLon} in reverse order.
441 '''
442 return self._reversed()
444 __str__ = __repr__
446 def count(self, latlon):
447 '''Count the number of rows with a specific lat-/longitude.
449 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
450 C{(lat, lon)}).
452 @return: Count (C{int}).
454 @raise TypeError: Invalid B{C{latlon}}.
455 '''
456 return self._count(latlon)
458 def find(self, latlon, *start_end):
459 '''Find the first row with a specific lat-/longitude.
461 @arg latlon: Point (C{LatLon}) or 2-tuple (lat, lon).
462 @arg start_end: Optional C{[start[, end]]} index (integers).
464 @return: Index or -1 if not found (C{int}).
466 @raise TypeError: Invalid B{C{latlon}}.
467 '''
468 return self._find(latlon, start_end)
470 def _findall(self, latlon, start_end):
471 '''(INTERNAL) Yield indices of all matching rows.
472 '''
473 try:
474 lat, lon = latlon.lat, latlon.lon
475 except AttributeError:
476 try:
477 lat, lon = latlon
478 except (TypeError, ValueError):
479 raise _IsnotError(_valid_, latlon=latlon)
481 _ilat, _ilon = self._ilat, self._ilon
482 _array, _eps = self._array, self._epsilon
483 for i in self._range(*start_end):
484 row = _array[i]
485 if fabs(row[_ilat] - lat) <= _eps and \
486 fabs(row[_ilon] - lon) <= _eps:
487 yield i
489 def findall(self, latlon, *start_end):
490 '''Yield indices of all rows with a specific lat-/longitude.
492 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
493 C{(lat, lon)}).
494 @arg start_end: Optional C{[start[, end]]} index (C{int}).
496 @return: Indices (C{iterable}).
498 @raise TypeError: Invalid B{C{latlon}}.
499 '''
500 return self._findall(latlon, start_end)
502 def index(self, latlon, *start_end): # PYCHOK Python 2- issue
503 '''Find index of the first row with a specific lat-/longitude.
505 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
506 C{(lat, lon)}).
507 @arg start_end: Optional C{[start[, end]]} index (C{int}).
509 @return: Index (C{int}).
511 @raise IndexError: Point not found.
513 @raise TypeError: Invalid B{C{latlon}}.
514 '''
515 return self._index(latlon, start_end)
517 @Property_RO
518 def ilat(self):
519 '''Get the latitudes column index (C{int}).
520 '''
521 return self._ilat
523 @Property_RO
524 def ilon(self):
525 '''Get the longitudes column index (C{int}).
526 '''
527 return self._ilon
529# next = __iter__
531 def point(self, row): # PYCHOK *attrs
532 '''Instantiate a point C{LatLon}.
534 @arg row: Array row (numpy.array).
536 @return: Point (C{LatLon}).
537 '''
538 return self._LatLon(row[self._ilat], row[self._ilon])
540 def rfind(self, latlon, *start_end):
541 '''Find the last row with a specific lat-/longitude.
543 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
544 C{(lat, lon)}).
545 @arg start_end: Optional C{[start[, end]]} index (C{int}).
547 @note: Keyword order, first stop, then start.
549 @return: Index or -1 if not found (C{int}).
551 @raise TypeError: Invalid B{C{latlon}}.
552 '''
553 return self._rfind(latlon, start_end)
555 def _slicekwds(self):
556 '''(INTERNAL) Slice kwds.
557 '''
558 return dict(ilat=self._ilat, ilon=self._ilon)
560 @Property_RO
561 def shape(self):
562 '''Get the shape of the C{NumPy} array or the C{Tuples} as
563 L{Shape2Tuple}C{(nrows, ncols)}.
564 '''
565 return self._shape
567 def _subset(self, indices): # PYCHOK no cover
568 '''(INTERNAL) I{Must be implemented/overloaded}.
569 '''
570 notImplemented(self, indices)
572 def subset(self, indices):
573 '''Return a subset of the C{NumPy} array.
575 @arg indices: Row indices (C{range} or C{int}[]).
577 @note: A C{subset} is different from a C{slice} in 2 ways:
578 (a) the C{subset} is typically specified as a list of
579 (un-)ordered indices and (b) the C{subset} allocates
580 a new, separate C{NumPy} array while a C{slice} is
581 just an other C{view} of the original C{NumPy} array.
583 @return: Sub-array (C{numpy.array}).
585 @raise IndexError: Out-of-range B{C{indices}} value.
587 @raise TypeError: If B{C{indices}} is not a C{range}
588 nor an C{int}[].
589 '''
590 if not issequence(indices, tuple): # NO tuple, only list
591 # and range work properly to get Numpy array sub-sets
592 raise _IsnotError(_valid_, indices=type(indices))
594 n = len(self)
595 for i, v in enumerate(indices):
596 if not isint(v):
597 raise _TypeError(Fmt.SQUARE(indices=i), v)
598 elif not 0 <= v < n:
599 raise _IndexError(Fmt.SQUARE(indices=i), v)
601 return self._subset(indices)
604class LatLon2psxy(_Basequence):
605 '''Wrapper for C{LatLon} points as "on-the-fly" pseudo-xy coordinates.
606 '''
607 _closed = False
608 _len = 0
609 _deg2m = None # default, keep degrees
610 _radius = None
611 _wrap = True
613 def __init__(self, latlons, closed=False, radius=None, wrap=True):
614 '''Handle C{LatLon} points as pseudo-xy coordinates.
616 @note: The C{LatLon} latitude is considered the I{pseudo-y}
617 and longitude the I{pseudo-x} coordinate, likewise
618 for L{LatLon2Tuple}. However, 2-tuples C{(x, y)} are
619 considered as I{(longitude, latitude)}.
621 @arg latlons: Points C{list}, C{sequence}, C{set}, C{tuple},
622 etc. (C{LatLon[]}).
623 @kwarg closed: Optionally, close the polygon (C{bool}).
624 @kwarg radius: Mean earth radius (C{meter}).
625 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
626 the B{C{latlons}} points (C{bool}).
628 @raise PointsError: Insufficient number of B{C{latlons}}.
630 @raise TypeError: Some B{C{points}} are not B{C{base}}.
631 '''
632 self._closed = closed
633 self._len, self._array = points2(latlons, closed=closed)
634 if radius:
635 self._radius = r = Radius(radius)
636 self._deg2m = degrees2m(_1_0, r)
637 if not wrap:
638 self._wrap = False
640 def __contains__(self, xy):
641 '''Check for a matching point.
643 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
644 C{(x, y)}) in (C{degrees}.
646 @return: C{True} if B{C{xy}} is present, C{False} otherwise.
648 @raise TypeError: Invalid B{C{xy}}.
649 '''
650 return self._contains(xy)
652 def __getitem__(self, index):
653 '''Return the pseudo-xy or return a L{LatLon2psxy} slice.
654 '''
655 return self._getitem(index)
657 def __iter__(self):
658 '''Yield all pseudo-xy's.
659 '''
660 return self._iter()
662 def __len__(self):
663 '''Return the number of pseudo-xy's.
664 '''
665 return self._len
667 def __repr__(self):
668 '''Return a string representation.
669 '''
670 return self._repr()
672 def __reversed__(self): # PYCHOK false
673 '''Yield all pseudo-xy's in reverse order.
674 '''
675 return self._reversed()
677 __str__ = __repr__
679 def count(self, xy):
680 '''Count the number of matching points.
682 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
683 C{(x, y)}) in (C{degrees}.
685 @return: Count (C{int}).
687 @raise TypeError: Invalid B{C{xy}}.
688 '''
689 return self._count(xy)
691 def find(self, xy, *start_end):
692 '''Find the first matching point.
694 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
695 C{(x, y)}) in (C{degrees}.
696 @arg start_end: Optional C{[start[, end]]} index (C{int}).
698 @return: Index or -1 if not found (C{int}).
700 @raise TypeError: Invalid B{C{xy}}.
701 '''
702 return self._find(xy, start_end)
704 def _findall(self, xy, start_end):
705 '''(INTERNAL) Yield indices of all matching points.
706 '''
707 try:
708 x, y = xy.lon, xy.lat
710 def _x_y_ll3(ll): # match LatLon
711 return ll.lon, ll.lat, ll
713 except AttributeError:
714 try:
715 x, y = xy[:2]
716 except (IndexError, TypeError, ValueError):
717 raise _IsnotError(_valid_, xy=xy)
719 _x_y_ll3 = self.point # PYCHOK expected
721 _array, _eps = self._array, self._epsilon
722 for i in self._range(*start_end):
723 xi, yi, _ = _x_y_ll3(_array[i])
724 if fabs(xi - x) <= _eps and \
725 fabs(yi - y) <= _eps:
726 yield i
728 def findall(self, xy, *start_end):
729 '''Yield indices of all matching points.
731 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
732 C{(x, y)}) in (C{degrees}.
733 @arg start_end: Optional C{[start[, end]]} index (C{int}).
735 @return: Indices (C{iterator}).
737 @raise TypeError: Invalid B{C{xy}}.
738 '''
739 return self._findall(xy, start_end)
741 def index(self, xy, *start_end): # PYCHOK Python 2- issue
742 '''Find the first matching point.
744 @arg xy: Point (C{LatLon}) or 2-tuple (x, y) in (C{degrees}).
745 @arg start_end: Optional C{[start[, end]]} index (C{int}).
747 @return: Index (C{int}).
749 @raise IndexError: Point not found.
751 @raise TypeError: Invalid B{C{xy}}.
752 '''
753 return self._index(xy, start_end)
755 @property_RO
756 def isPoints2(self):
757 '''Is this a LatLon2 wrapper/converter?
758 '''
759 return True # isinstance(self, (LatLon2psxy, ...))
761 def point(self, ll): # PYCHOK *attrs
762 '''Create a pseudo-xy.
764 @arg ll: Point (C{LatLon}).
766 @return: An L{Point3Tuple}C{(x, y, ll)}.
767 '''
768 x, y = ll.lon, ll.lat # note, x, y = lon, lat
769 if self._wrap:
770 y, x = _Wrap.latlon(y, x)
771 d = self._deg2m
772 if d: # convert degrees to meter (or radians)
773 x *= d
774 y *= d
775 return Point3Tuple(x, y, ll)
777 def rfind(self, xy, *start_end):
778 '''Find the last matching point.
780 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
781 C{(x, y)}) in (C{degrees}.
782 @arg start_end: Optional C{[start[, end]]} index (C{int}).
784 @return: Index or -1 if not found (C{int}).
786 @raise TypeError: Invalid B{C{xy}}.
787 '''
788 return self._rfind(xy, start_end)
790 def _slicekwds(self):
791 '''(INTERNAL) Slice kwds.
792 '''
793 return dict(closed=self._closed, radius=self._radius, wrap=self._wrap)
796class Numpy2LatLon(_Array2LatLon): # immutable, on purpose
797 '''Wrapper for C{NumPy} arrays as "on-the-fly" C{LatLon} points.
798 '''
799 def __init__(self, array, ilat=0, ilon=1, LatLon=None):
800 '''Handle a C{NumPy} array as a sequence of C{LatLon} points.
802 @arg array: C{NumPy} array (C{numpy.array}).
803 @kwarg ilat: Optional index of the latitudes column (C{int}).
804 @kwarg ilon: Optional index of the longitudes column (C{int}).
805 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
807 @raise IndexError: If B{C{array.shape}} is not (1+, 2+).
809 @raise TypeError: If B{C{array}} is not a C{NumPy} array or
810 C{LatLon} is not a class with C{lat}
811 and C{lon} attributes.
813 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
814 are the same or out of range.
816 @example:
818 >>> type(array)
819 <type 'numpy.ndarray'> # <class ...> in Python 3+
820 >>> points = Numpy2LatLon(array, lat=0, lon=1)
821 >>> simply = simplifyRDP(points, ...)
822 >>> type(simply)
823 <type 'numpy.ndarray'> # <class ...> in Python 3+
824 >>> sliced = points[1:-1]
825 >>> type(sliced)
826 <class '...Numpy2LatLon'>
827 '''
828 try: # get shape and check some other numpy.array attrs
829 s, _, _ = array.shape, array.nbytes, array.ndim # PYCHOK expected
830 except AttributeError:
831 raise _IsnotError('NumPy', array=type(array))
833 _Array2LatLon.__init__(self, array, ilat=ilat, ilon=ilon,
834 LatLon=LatLon, shape=s)
836 @property_RO
837 def isNumpy2(self):
838 '''Is this a Numpy2 wrapper?
839 '''
840 return True # isinstance(self, (Numpy2LatLon, ...))
842 def _subset(self, indices):
843 return self._array[indices] # NumPy special
846class Shape2Tuple(_NamedTuple):
847 '''2-Tuple C{(nrows, ncols)}, the number of rows and columns,
848 both C{int}.
849 '''
850 _Names_ = (_nrows_, _ncols_)
851 _Units_ = ( Number_, Number_)
854class Tuple2LatLon(_Array2LatLon):
855 '''Wrapper for tuple sequences as "on-the-fly" C{LatLon} points.
856 '''
857 def __init__(self, tuples, ilat=0, ilon=1, LatLon=None):
858 '''Handle a list of tuples, each containing a lat- and longitude
859 and perhaps other values as a sequence of C{LatLon} points.
861 @arg tuples: The C{list}, C{tuple} or C{sequence} of tuples (C{tuple}[]).
862 @kwarg ilat: Optional index of the latitudes value (C{int}).
863 @kwarg ilon: Optional index of the longitudes value (C{int}).
864 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
866 @raise IndexError: If C{(len(B{tuples}), min(len(t) for t
867 in B{tuples}))} is not (1+, 2+).
869 @raise TypeError: If B{C{tuples}} is not a C{list}, C{tuple}
870 or C{sequence} or if B{C{LatLon}} is not a
871 C{LatLon} with C{lat}, C{lon} and C{name}
872 attributes.
874 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
875 are the same or out of range.
877 @example:
879 >>> tuples = [(0, 1), (2, 3), (4, 5)]
880 >>> type(tuples)
881 <type 'list'> # <class ...> in Python 3+
882 >>> points = Tuple2LatLon(tuples, lat=0, lon=1)
883 >>> simply = simplifyRW(points, 0.5, ...)
884 >>> type(simply)
885 <type 'list'> # <class ...> in Python 3+
886 >>> simply
887 [(0, 1), (4, 5)]
888 >>> sliced = points[1:-1]
889 >>> type(sliced)
890 <class '...Tuple2LatLon'>
891 >>> sliced
892 ...Tuple2LatLon([(2, 3), ...][1], ilat=0, ilon=1)
894 >>> closest, _ = nearestOn2(LatLon_(2, 1), points, adjust=False)
895 >>> closest
896 LatLon_(lat=1.0, lon=2.0)
898 >>> closest, _ = nearestOn2(LatLon_(3, 2), points)
899 >>> closest
900 LatLon_(lat=2.001162, lon=3.001162)
901 '''
902 _xinstanceof(list, tuple, tuples=tuples)
903 s = len(tuples), min(len(_) for _ in tuples)
904 _Array2LatLon.__init__(self, tuples, ilat=ilat, ilon=ilon,
905 LatLon=LatLon, shape=s)
907 @property_RO
908 def isTuple2(self):
909 '''Is this a Tuple2 wrapper?
910 '''
911 return True # isinstance(self, (Tuple2LatLon, ...))
913 def _subset(self, indices):
914 return type(self._array)(self._array[i] for i in indices)
917def _area2(points, adjust, wrap):
918 '''(INTERNAL) Approximate the area in radians squared, I{signed}.
919 '''
920 if adjust:
921 # approximate trapezoid by a rectangle, adjusting
922 # the top width by the cosine of the latitudinal
923 # average and bottom width by some fudge factor
924 def _adjust(w, h):
925 c = cos(h) if fabs(h) < PI_2 else _0_0
926 return w * h * (c + 1.2876) * _0_5
927 else:
928 def _adjust(w, h): # PYCHOK expected
929 return w * h
931 # setting radius=1 converts degrees to radians
932 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)
933 x1, y1, ll = Ps[0]
934 pts = [ll] # for _areaError
936 A2 = Fsum() # trapezoidal area in radians**2
937 for p in Ps.iterate(closed=True):
938 x2, y2, ll = p
939 if len(pts) < 4:
940 pts.append(ll)
941 w, x2 = unrollPI(x1, x2, wrap=wrap and not Ps.looped)
942 A2 += _adjust(w, (y2 + y1) * _0_5)
943 x1, y1 = x2, y2
945 return A2.fsum(), tuple(pts)
948def _areaError(pts, near_=NN): # imported by .ellipsoidalKarney
949 '''(INTERNAL) Area issue.
950 '''
951 t = _ELLIPSIS_(pts[:3], NN)
952 return _ValueError(NN(near_, 'zero or polar area'), txt=t)
955def areaOf(points, adjust=True, radius=R_M, wrap=True):
956 '''Approximate the area of a polygon or composite.
958 @arg points: The polygon points or clips (C{LatLon}[],
959 L{BooleanFHP} or L{BooleanGH}).
960 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
961 by the cosine of the mean latitude (C{bool}).
962 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
963 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
964 the B{C{points}} (C{bool}).
966 @return: Approximate area (I{square} C{meter}, same units as
967 B{C{radius}} or C{radians} I{squared} if B{C{radius}}
968 is C{None}).
970 @raise PointsError: Insufficient number of B{C{points}}
972 @raise TypeError: Some B{C{points}} are not C{LatLon}.
974 @raise ValueError: Invalid B{C{radius}}.
976 @note: This area approximation has limited accuracy and is
977 ill-suited for regions exceeding several hundred Km
978 or Miles or with near-polar latitudes.
980 @see: L{sphericalNvector.areaOf}, L{sphericalTrigonometry.areaOf},
981 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.
982 '''
983 if _MODS.booleans.isBoolean(points):
984 a = points._sum1(areaOf, adjust=adjust, radius=None, wrap=wrap)
985 else:
986 a, _ = _area2(points, adjust, wrap)
987 return fabs(a if radius is None else (Radius(radius)**2 * a))
990def boundsOf(points, wrap=False, LatLon=None): # was=True
991 '''Determine the bottom-left SW and top-right NE corners of a
992 path or polygon.
994 @arg points: The path or polygon points (C{LatLon}[]).
995 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
996 the B{C{points}} (C{bool}).
997 @kwarg LatLon: Optional class to return the C{bounds}
998 corners (C{LatLon}) or C{None}.
1000 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)} as
1001 B{C{LatLon}}s if B{C{LatLon}} is C{None} a
1002 L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}.
1004 @raise PointsError: Insufficient number of B{C{points}}
1006 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1008 @see: Function L{quadOf}.
1010 @example:
1012 >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1013 >>> boundsOf(b) # False
1014 >>> 45.0, 1.0, 46.0, 2.0
1015 '''
1016 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)
1017 w, s, _ = e, n, _ = Ps[0]
1019 v = w
1020 for x, y, _ in Ps.iterate(closed=False): # [1:]
1021 if wrap:
1022 _, x = unroll180(v, x, wrap=True)
1023 v = x
1025 if w > x:
1026 w = x
1027 elif e < x:
1028 e = x
1030 if s > y:
1031 s = y
1032 elif n < y:
1033 n = y
1035 return Bounds4Tuple(s, w, n, e) if LatLon is None else \
1036 Bounds2Tuple(LatLon(s, w), LatLon(n, e)) # PYCHOK inconsistent
1039def _distanceTo(Error, **name_points): # .frechet, .hausdorff, .heights
1040 '''(INTERNAL) Chech callable C{distanceTo} methods.
1041 '''
1042 name, ps = name_points.popitem()
1043 for i, p in enumerate(ps):
1044 if not callable(_xattr(p, distanceTo=None)):
1045 n = _distanceTo.__name__[1:]
1046 t = _SPACE_(_no_, callable.__name__, n)
1047 raise Error(Fmt.SQUARE(name, i), p, txt=t)
1048 return ps
1051def centroidOf(points, wrap=False, LatLon=None): # was=True
1052 '''Determine the centroid of a polygon.
1054 @arg points: The polygon points (C{LatLon}[]).
1055 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1056 B{C{points}} (C{bool}).
1057 @kwarg LatLon: Optional class to return the centroid (C{LatLon})
1058 or C{None}.
1060 @return: Centroid (B{C{LatLon}}) or a L{LatLon2Tuple}C{(lat, lon)}
1061 if C{B{LatLon} is None}.
1063 @raise PointsError: Insufficient number of B{C{points}}
1065 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1067 @raise ValueError: The B{C{points}} enclose a pole or
1068 near-zero area.
1070 @see: U{Centroid<https://WikiPedia.org/wiki/Centroid#Of_a_polygon>} and
1071 Paul Bourke's U{Calculating The Area And Centroid Of A Polygon
1072 <https://www.SEAS.UPenn.edu/~ese502/lab-content/extra_materials/
1073 Polygon%20Area%20and%20Centroid.pdf>}, 1988.
1074 '''
1075 A, X, Y = Fsum(), Fsum(), Fsum()
1077 # setting radius=1 converts degrees to radians
1078 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)
1079 x1, y1, ll = Ps[0]
1080 pts = [ll] # for _areaError
1081 for p in Ps.iterate(closed=True):
1082 x2, y2, ll = p
1083 if len(pts) < 4:
1084 pts.append(ll)
1085 if wrap and not Ps.looped:
1086 _, x2 = unrollPI(x1, x2, wrap=True)
1087 t = x1 * y2 - x2 * y1
1088 A += t
1089 X += t * (x1 + x2)
1090 Y += t * (y1 + y2)
1091 # XXX more elaborately:
1092 # t1, t2 = x1 * y2, -(x2 * y1)
1093 # A.fadd_(t1, t2)
1094 # X.fadd_(t1 * x1, t1 * x2, t2 * x1, t2 * x2)
1095 # Y.fadd_(t1 * y1, t1 * y2, t2 * y1, t2 * y2)
1096 x1, y1 = x2, y2
1098 a = A.fmul(_6_0).fover(_2_0)
1099 if isnear0(a):
1100 raise _areaError(pts, near_=_near_)
1101 y, x = degrees90(Y.fover(a)), degrees180(X.fover(a))
1102 return LatLon2Tuple(y, x) if LatLon is None else LatLon(y, x)
1105def fractional(points, fi, j=None, wrap=None, LatLon=None, Vector=None, **kwds):
1106 '''Return the point at a given I{fractional} index.
1108 @arg points: The points (C{LatLon}[], L{Numpy2LatLon}[],
1109 L{Tuple2LatLon}[], C{Cartesian}[], C{Vector3d}[],
1110 L{Vector3Tuple}[]).
1111 @arg fi: The fractional index (L{FIx}, C{float} or C{int}).
1112 @kwarg j: Optionally, index of the other point (C{int}).
1113 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1114 B{{points}} (C{bool}) or C{None} for a backward
1115 compatible L{LatLon2Tuple} or B{C{LatLon}} with
1116 averaged lat- and longitudes. Use C{True} or
1117 C{False} to get the I{fractional} point computed
1118 by method C{B{points}[fi].intermediateTo}.
1119 @kwarg LatLon: Optional class to return the I{intermediate},
1120 I{fractional} point (C{LatLon}) or C{None}.
1121 @kwarg Vector: Optional class to return the I{intermediate},
1122 I{fractional} point (C{Cartesian}, C{Vector3d})
1123 or C{None}.
1124 @kwarg kwds: Optional, additional B{C{LatLon}} I{or} B{C{Vector}}
1125 keyword arguments, ignored if both C{B{LatLon}} and
1126 C{B{Vector}} are C{None}.
1128 @return: A L{LatLon2Tuple}C{(lat, lon)} if B{C{wrap}}, B{C{LatLon}}
1129 and B{C{Vector}} all are C{None}, the defaults.
1131 An instance of B{C{LatLon}} if not C{None} I{or} an instance
1132 of B{C{Vector}} if not C{None}.
1134 Otherwise with B{C{wrap}} either C{True} or C{False} and
1135 B{C{LatLon}} and B{C{Vector}} both C{None}, an instance of
1136 B{C{points}}' (sub-)class C{intermediateTo} I{fractional}.
1138 Summarized as follows:
1140 >>> wrap | LatLon | Vector | returned type/value
1141 # -------+--------+--------+--------------+------
1142 # | | | LatLon2Tuple | favg
1143 # None | None | None | or** |
1144 # | | | Vector3Tuple | favg
1145 # None | LatLon | None | LatLon | favg
1146 # None | None | Vector | Vector | favg
1147 # -------+--------+--------+--------------+------
1148 # True | None | None | points' | .iTo
1149 # True | LatLon | None | LatLon | .iTo
1150 # True | None | Vector | Vector | .iTo
1151 # -------+--------+--------+--------------+------
1152 # False | None | None | points' | .iTo
1153 # False | LatLon | None | LatLon | .iTo
1154 # False | None | Vector | Vector | .iTo
1155 # _____
1156 # favg) averaged lat, lon or x, y, z values
1157 # .iTo) value from points[fi].intermediateTo
1158 # **) depends on base class of points[fi]
1160 @raise IndexError: Fractional index B{C{fi}} invalid or B{C{points}}
1161 not subscriptable or not closed.
1163 @raise TypeError: Invalid B{C{LatLon}}, B{C{Vector}} or B{C{kwds}}
1164 argument.
1166 @see: Class L{FIx} and method L{FIx.fractional}.
1167 '''
1168 if LatLon and Vector: # PYCHOK no cover
1169 kwds = _xkwds(kwds, fi=fi, LatLon=LatLon, Vector=Vector)
1170 raise _TypeError(txt=fractional.__name__, **kwds)
1171 w = wrap if LatLon else False # intermediateTo
1172 try:
1173 if not isscalar(fi) or fi < 0:
1174 raise IndexError
1175 n = _xattr(fi, fin=0)
1176 p = _fractional(points, fi, j, fin=n, wrap=w) # see .units.FIx
1177 if LatLon:
1178 p = LatLon(p.lat, p.lon, **kwds)
1179 elif Vector:
1180 p = Vector(p.x, p.y, p.z, **kwds)
1181 except (IndexError, TypeError):
1182 raise _IndexError(fi=fi, points=points, wrap=w, txt=fractional.__name__)
1183 return p
1186def _fractional(points, fi, j, fin=None, wrap=None): # in .frechet.py
1187 '''(INTERNAL) Compute point at L{fractional} index C{fi} and C{j}.
1188 '''
1189 i = int(fi)
1190 p = points[i]
1191 r = fi - float(i)
1192 if r > EPS: # EPS0?
1193 if j is None: # in .frechet.py
1194 j = i + 1
1195 if fin:
1196 j %= fin
1197 q = points[j]
1198 if r >= EPS1: # PYCHOK no cover
1199 p = q
1200 elif wrap is not None: # in (True, False)
1201 p = p.intermediateTo(q, r, wrap=wrap)
1202 elif _isLatLon(p): # backward compatible default
1203 p = LatLon2Tuple(favg(p.lat, q.lat, f=r),
1204 favg(p.lon, q.lon, f=r),
1205 name=fractional.__name__)
1206 else: # assume p and q are cartesian or vectorial
1207 z = p.z if p.z is q.z else favg(p.z, q.z, f=r)
1208 p = Vector3Tuple(favg(p.x, q.x, f=r),
1209 favg(p.y, q.y, f=r), z,
1210 name=fractional.__name__)
1211 return p
1214def isclockwise(points, adjust=False, wrap=True):
1215 '''Determine the direction of a path or polygon.
1217 @arg points: The path or polygon points (C{LatLon}[]).
1218 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1219 by the cosine of the mean latitude (C{bool}).
1220 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1221 B{C{points}} (C{bool}).
1223 @return: C{True} if B{C{points}} are clockwise, C{False} otherwise.
1225 @raise PointsError: Insufficient number of B{C{points}}
1227 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1229 @raise ValueError: The B{C{points}} enclose a pole or zero area.
1231 @example:
1233 >>> f = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1234 >>> isclockwise(f) # False
1235 >>> isclockwise(reversed(f)) # True
1236 '''
1237 a, pts = _area2(points, adjust, wrap)
1238 if a > 0: # opposite of ellipsoidalExact and -Karney
1239 return True
1240 elif a < 0:
1241 return False
1242 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1243 raise _areaError(pts)
1246def isconvex(points, adjust=False, wrap=False): # was=True
1247 '''Determine whether a polygon is convex.
1249 @arg points: The polygon points (C{LatLon}[]).
1250 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1251 by the cosine of the mean latitude (C{bool}).
1252 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1253 B{C{points}} (C{bool}).
1255 @return: C{True} if B{C{points}} are convex, C{False} otherwise.
1257 @raise CrossError: Some B{C{points}} are colinear.
1259 @raise PointsError: Insufficient number of B{C{points}}
1261 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1263 @example:
1265 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1266 >>> isconvex(t) # True
1268 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1269 >>> isconvex(f) # False
1270 '''
1271 return bool(isconvex_(points, adjust=adjust, wrap=wrap))
1274def isconvex_(points, adjust=False, wrap=False): # was=True
1275 '''Determine whether a polygon is convex I{and clockwise}.
1277 @arg points: The polygon points (C{LatLon}[]).
1278 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1279 by the cosine of the mean latitude (C{bool}).
1280 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1281 B{C{points}} (C{bool}).
1283 @return: C{+1} if B{C{points}} are convex clockwise, C{-1} for
1284 convex counter-clockwise B{C{points}}, C{0} otherwise.
1286 @raise CrossError: Some B{C{points}} are colinear.
1288 @raise PointsError: Insufficient number of B{C{points}}
1290 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1292 @example:
1294 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1295 >>> isconvex_(t) # +1
1297 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1298 >>> isconvex_(f) # 0
1299 '''
1300 if adjust:
1301 def _unroll2(x1, x2, w, y1, y2):
1302 x21, x2 = unroll180(x1, x2, wrap=w)
1303 y = radians(y1 + y2) * _0_5
1304 x21 *= cos(y) if fabs(y) < PI_2 else _0_0
1305 return x21, x2
1306 else:
1307 def _unroll2(x1, x2, w, *unused): # PYCHOK expected
1308 return unroll180(x1, x2, wrap=w)
1310 c, s = crosserrors(), 0
1312 Ps = LatLon2PsxyIter(points, loop=2, wrap=wrap)
1313 x1, y1, _ = Ps[0]
1314 x2, y2, _ = Ps[1]
1316 x21, x2 = _unroll2(x1, x2, False, y1, y2)
1317 for i, p in Ps.enumerate(closed=True):
1318 x3, y3, ll = p
1319 x32, x3 = _unroll2(x2, x3, bool(wrap and not Ps.looped), y2, y3)
1321 # get the sign of the distance from point
1322 # x3, y3 to the line from x1, y1 to x2, y2
1323 # <https://WikiPedia.org/wiki/Distance_from_a_point_to_a_line>
1324 s3 = fdot((x3, y3, x1, y1), y2 - y1, -x21, -y2, x2)
1325 if s3 > 0: # x3, y3 on the right
1326 if s < 0: # non-convex
1327 return 0
1328 s = +1
1330 elif s3 < 0: # x3, y3 on the left
1331 if s > 0: # non-convex
1332 return 0
1333 s = -1
1335 elif c and fdot((x32, y1 - y2), y3 - y2, -x21) < 0: # PYCHOK no cover
1336 # colinear u-turn: x3, y3 not on the
1337 # opposite side of x2, y2 as x1, y1
1338 t = Fmt.SQUARE(points=i)
1339 raise CrossError(t, ll, txt=_colinear_)
1341 x1, y1, x2, y2, x21 = x2, y2, x3, y3, x32
1343 return s # all points on the same side
1346def isenclosedBy(point, points, wrap=False): # MCCABE 15
1347 '''Determine whether a point is enclosed by a polygon or composite.
1349 @arg point: The point (C{LatLon} or 2-tuple C{(lat, lon)}).
1350 @arg points: The polygon points or clips (C{LatLon}[], L{BooleanFHP}
1351 or L{BooleanGH}).
1352 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1353 B{C{points}} (C{bool}).
1355 @return: C{True} if the B{C{point}} is inside the polygon or
1356 composite, C{False} otherwise.
1358 @raise PointsError: Insufficient number of B{C{points}}
1360 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1362 @raise ValueError: Invalid B{C{point}}, lat- or longitude.
1364 @see: Functions L{pygeodesy.isconvex} and L{pygeodesy.ispolar} especially
1365 if the B{C{points}} may enclose a pole or wrap around the earth
1366 I{longitudinally}, methods L{sphericalNvector.LatLon.isenclosedBy},
1367 L{sphericalTrigonometry.LatLon.isenclosedBy} and U{MultiDop
1368 GeogContainPt<https://GitHub.com/NASA/MultiDop>} (U{Shapiro et.al. 2009,
1369 JTECH<https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}
1370 and U{Potvin et al. 2012, JTECH <https://Journals.AMetSoc.org/doi/abs/
1371 10.1175/JTECH-D-11-00019.1>}).
1372 '''
1373 try:
1374 y0, x0 = point.lat, point.lon
1375 except AttributeError:
1376 try:
1377 y0, x0 = map(float, point[:2])
1378 except (IndexError, TypeError, ValueError) as x:
1379 raise _ValueError(point=point, cause=x)
1381 if wrap:
1382 y0, x0 = _Wrap.latlon(y0, x0)
1384 def _dxy3(x, x2, y2, Ps):
1385 dx, x2 = unroll180(x, x2, wrap=not Ps.looped)
1386 return dx, x2, y2
1388 else:
1389 x0 = fmod(x0, _360_0) # not x0 % 360!
1390 x0_180_ = x0 - _180_0
1391 x0_180 = x0 + _180_0
1393 def _dxy3(x1, x, y, unused): # PYCHOK expected
1394 x = _umod_360(float(x))
1395 if x < x0_180_:
1396 x += _360_0
1397 elif x >= x0_180:
1398 x -= _360_0
1399 return (x - x1), x, y
1401 if _MODS.booleans.isBoolean(points):
1402 return points._encloses(y0, x0, wrap=wrap)
1404 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)
1405 p = Ps[0]
1406 e = m = False
1407 S = Fsum()
1409 _, x1, y1 = _dxy3(x0, p.x, p.y, False)
1410 for p in Ps.iterate(closed=True):
1411 dx, x2, y2 = _dxy3(x1, p.x, p.y, Ps)
1412 # ignore duplicate and near-duplicate pts
1413 if fabs(dx) > EPS or fabs(y2 - y1) > EPS:
1414 # determine if polygon edge (x1, y1)..(x2, y2) straddles
1415 # point (lat, lon) or is on boundary, but do not count
1416 # edges on boundary as more than one crossing
1417 if fabs(dx) < 180 and (x1 < x0 <= x2 or x2 < x0 <= x1):
1418 m = not m
1419 dy = (x0 - x1) * (y2 - y1) - (y0 - y1) * dx
1420 if (dy > 0 and dx >= 0) or (dy < 0 and dx <= 0):
1421 e = not e
1423 S += sin(radians(y2))
1424 x1, y1 = x2, y2
1426 # An odd number of meridian crossings means, the polygon
1427 # contains a pole. Assume it is the pole on the hemisphere
1428 # containing the polygon mean point and if the polygon does
1429 # contain the North Pole, flip the result.
1430 if m and S.fsum() > 0:
1431 e = not e
1432 return e
1435def ispolar(points, wrap=False):
1436 '''Check whether a polygon encloses a pole.
1438 @arg points: The polygon points (C{LatLon}[]).
1439 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1440 the B{C{points}} (C{bool}).
1442 @return: C{True} if the polygon encloses a pole, C{False}
1443 otherwise.
1445 @raise PointsError: Insufficient number of B{C{points}}
1447 @raise TypeError: Some B{C{points}} are not C{LatLon} or don't
1448 have C{bearingTo2}, C{initialBearingTo}
1449 and C{finalBearingTo} methods.
1450 '''
1451 def _cds(ps, w): # iterate over course deltas
1452 Ps = PointsIter(ps, loop=2, wrap=w)
1453 p2, p1 = Ps[0:2]
1454 b1, _ = _bearingTo2(p2, p1, wrap=False)
1455 for p2 in Ps.iterate(closed=True):
1456 if not p2.isequalTo(p1, EPS):
1457 if w and not Ps.looped:
1458 p2 = _unrollon(p1, p2)
1459 b, b2 = _bearingTo2(p1, p2, wrap=False)
1460 yield wrap180(b - b1) # (b - b1 + 540) % 360 - 180
1461 yield wrap180(b2 - b) # (b2 - b + 540) % 360 - 180
1462 p1, b1 = p2, b2
1464 # summation of course deltas around pole is 0° rather than normally ±360°
1465 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1466 s = fsum(_cds(points, wrap), floats=True)
1467 # XXX fix (intermittant) edge crossing pole - eg (85,90), (85,0), (85,-90)
1468 return fabs(s) < 90 # "zero-ish"
1471def luneOf(lon1, lon2, closed=False, LatLon=LatLon_, **LatLon_kwds):
1472 '''Generate an ellipsoidal or spherical U{lune
1473 <https://WikiPedia.org/wiki/Spherical_lune>}-shaped path or polygon.
1475 @arg lon1: Left longitude (C{degrees90}).
1476 @arg lon2: Right longitude (C{degrees90}).
1477 @kwarg closed: Optionally, close the path (C{bool}).
1478 @kwarg LatLon: Class to use (L{LatLon_}).
1479 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1480 keyword arguments.
1482 @return: A tuple of 4 or 5 B{C{LatLon}} instances outlining
1483 the lune shape.
1485 @see: U{Latitude-longitude quadrangle
1486 <https://www.MathWorks.com/help/map/ref/areaquad.html>}.
1487 '''
1488 t = (LatLon( _0_0, lon1, **LatLon_kwds),
1489 LatLon( _90_0, lon1, **LatLon_kwds),
1490 LatLon( _0_0, lon2, **LatLon_kwds),
1491 LatLon(_N_90_0, lon2, **LatLon_kwds))
1492 if closed:
1493 t += t[:1]
1494 return t
1497def nearestOn5(point, points, closed=False, wrap=False, adjust=True,
1498 limit=9, **LatLon_and_kwds):
1499 '''Locate the point on a path or polygon closest to a reference point.
1501 The closest point on each polygon edge is either the nearest of that
1502 edge's end points or a point in between.
1504 @arg point: The reference point (C{LatLon}).
1505 @arg points: The path or polygon points (C{LatLon}[]).
1506 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1507 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1508 B{C{points}} (C{bool}).
1509 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1510 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1511 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1512 spherical earth radius L{R_KM}).
1513 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use for
1514 the closest point and additional B{C{LatLon}} keyword
1515 arguments, ignored if C{B{LatLon}=None} or not given.
1517 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1518 {closest} point (B{C{LatLon}}) or if C{B{LatLon} is None},
1519 a L{NearestOn5Tuple}C{(lat, lon, distance, angle, height)}.
1520 The C{distance} is the L{pygeodesy.equirectangular} distance
1521 between the C{closest} and reference B{C{point}} in C{degrees}.
1522 The C{angle} from the B{C{point}} to the C{closest} is in
1523 compass C{degrees}, like function L{pygeodesy.compassAngle}.
1525 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1526 see function L{pygeodesy.equirectangular_}.
1528 @raise PointsError: Insufficient number of B{C{points}}
1530 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1532 @note: Distances are I{approximated} by function L{pygeodesy.equirectangular_}.
1533 For more accuracy use one of the C{LatLon.nearestOn6} methods.
1535 @see: Function L{pygeodesy.degrees2m}.
1536 '''
1537 def _d2yx4(p2, p1, u, alw):
1538 # w = wrap if (i < (n - 1) or not closed) else False
1539 # equirectangular_ returns a Distance4Tuple(distance
1540 # in degrees squared, delta lat, delta lon, p2.lon
1541 # unroll/wrap'd); the previous p2.lon unroll/wrap'd
1542 # is also applied to the next edge's p1.lon
1543 return equirectangular_(p1.lat, p1.lon + u,
1544 p2.lat, p2.lon, **alw)
1546 def _h(p): # get height or default 0
1547 return _xattr(p, height=0) or 0
1549 # 3-D version used in .vector3d._nearestOn2
1550 #
1551 # point (x, y) on axis rotated ccw by angle a:
1552 # x' = x * cos(a) + y * sin(a)
1553 # y' = y * cos(a) - x * sin(a)
1554 #
1555 # distance (w) along and (h) perpendicular to
1556 # a line thru point (dx, dy) and the origin:
1557 # d = hypot(dx, dy)
1558 # w = (x * dx + y * dy) / d
1559 # h = (y * dx - x * dy) / d
1560 #
1561 # closest point on that line thru (dx, dy):
1562 # xc = dx * w / d
1563 # yc = dy * w / d
1564 # or
1565 # xc = dx * f
1566 # yc = dy * f
1567 # with
1568 # f = w / d
1569 # or
1570 # f = (y * dy + x * dx) / hypot2(dx, dy)
1571 #
1572 # i.e. no need for sqrt or hypot
1574 Ps = PointsIter(points, loop=1, wrap=wrap)
1575 p1 = c = Ps[0]
1576 u1 = u = _0_0
1577 kw = dict(adjust=adjust, limit=limit, wrap=False)
1578 d, dy, dx, _ = _d2yx4(p1, point, u1, kw)
1579 for p2 in Ps.iterate(closed=closed):
1580 # iff wrapped, unroll lon1 (actually previous
1581 # lon2) like function unroll180/-PI would've
1582 if wrap:
1583 kw.update(wrap=not (closed and Ps.looped))
1584 d21, y21, x21, u2 = _d2yx4(p2, p1, u1, kw)
1585 if d21 > EPS:
1586 # distance point to p1, y01 and x01 negated
1587 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)
1588 if d2 > EPS:
1589 w2 = y01 * y21 + x01 * x21
1590 if w2 > 0:
1591 if w2 < d21:
1592 # closest is between p1 and p2, use
1593 # original delta's, not y21 and x21
1594 f = w2 / d21
1595 p1 = LatLon_(favg(p1.lat, p2.lat, f=f),
1596 favg(p1.lon, p2.lon + u2, f=f),
1597 height=favg(_h(p1), _h(p2), f=f))
1598 u1 = _0_0
1599 else: # p2 is closest
1600 p1, u1 = p2, u2
1601 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)
1602 if d2 < d: # p1 is closer, y01 and x01 negated
1603 c, u, d, dy, dx = p1, u1, d2, -y01, -x01
1604 p1, u1 = p2, u2
1606 a = atan2b(dx, dy) # azimuth
1607 d = hypot( dx, dy)
1608 h = _h(c)
1609 n = nameof(point) or nearestOn5.__name__
1610 if LatLon_and_kwds:
1611 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1612 LL = _xkwds_pop(kwds, LatLon=None)
1613 if LL is not None:
1614 r = LL(c.lat, c.lon + u, **kwds)
1615 return NearestOn3Tuple(r, d, a, name=n)
1616 return NearestOn5Tuple(c.lat, c.lon + u, d, a, h, name=n) # PYCHOK expected
1619def perimeterOf(points, closed=False, adjust=True, radius=R_M, wrap=True):
1620 '''I{Approximate} the perimeter of a path, polygon. or composite.
1622 @arg points: The path or polygon points or clips (C{LatLon}[],
1623 L{BooleanFHP} or L{BooleanGH}).
1624 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1625 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1626 by the cosine of the mean latitude (C{bool}).
1627 @kwarg radius: Mean earth radius (C{meter}).
1628 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1629 B{C{points}} (C{bool}).
1631 @return: Approximate perimeter (C{meter}, same units as
1632 B{C{radius}}).
1634 @raise PointsError: Insufficient number of B{C{points}}
1636 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1638 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1639 C{B{points}} a composite.
1641 @note: This perimeter is based on the L{pygeodesy.equirectangular_}
1642 distance approximation and is ill-suited for regions exceeding
1643 several hundred Km or Miles or with near-polar latitudes.
1645 @see: Functions L{sphericalTrigonometry.perimeterOf} and
1646 L{ellipsoidalKarney.perimeterOf}.
1647 '''
1648 def _degs(ps, c, a, w): # angular edge lengths in degrees
1649 Ps = LatLon2PsxyIter(ps, loop=1) # wrap=w
1650 p1, u = Ps[0], _0_0 # previous x2's unroll/wrap
1651 for p2 in Ps.iterate(closed=c):
1652 if w and c:
1653 w = not Ps.looped
1654 # apply previous x2's unroll/wrap'd to new x1
1655 _, dy, dx, u = equirectangular_(p1.y, p1.x + u,
1656 p2.y, p2.x,
1657 adjust=a, limit=None,
1658 wrap=w) # PYCHOK non-seq
1659 yield hypot(dx, dy)
1660 p1 = p2
1662 if _MODS.booleans.isBoolean(points):
1663 if not closed:
1664 notImplemented(None, closed=closed, points=_composite_)
1665 d = points._sum1(perimeterOf, closed=True, adjust=adjust,
1666 radius=radius, wrap=wrap)
1667 else:
1668 d = fsum(_degs(points, closed, adjust, wrap), floats=True)
1669 return degrees2m(d, radius=radius)
1672def quadOf(latS, lonW, latN, lonE, closed=False, LatLon=LatLon_, **LatLon_kwds):
1673 '''Generate a quadrilateral path or polygon from two points.
1675 @arg latS: Souther-nmost latitude (C{degrees90}).
1676 @arg lonW: Western-most longitude (C{degrees180}).
1677 @arg latN: Norther-nmost latitude (C{degrees90}).
1678 @arg lonE: Eastern-most longitude (C{degrees180}).
1679 @kwarg closed: Optionally, close the path (C{bool}).
1680 @kwarg LatLon: Class to use (L{LatLon_}).
1681 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1682 keyword arguments.
1684 @return: Return a tuple of 4 or 5 B{C{LatLon}} instances
1685 outlining the quadrilateral.
1687 @see: Function L{boundsOf}.
1688 '''
1689 t = (LatLon(latS, lonW, **LatLon_kwds),
1690 LatLon(latN, lonW, **LatLon_kwds),
1691 LatLon(latN, lonE, **LatLon_kwds),
1692 LatLon(latS, lonE, **LatLon_kwds))
1693 if closed:
1694 t += t[:1]
1695 return t
1698__all__ += _ALL_DOCS(_Array2LatLon, _Basequence)
1700# **) MIT License
1701#
1702# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1703#
1704# Permission is hereby granted, free of charge, to any person obtaining a
1705# copy of this software and associated documentation files (the "Software"),
1706# to deal in the Software without restriction, including without limitation
1707# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1708# and/or sell copies of the Software, and to permit persons to whom the
1709# Software is furnished to do so, subject to the following conditions:
1710#
1711# The above copyright notice and this permission notice shall be included
1712# in all copies or substantial portions of the Software.
1713#
1714# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1715# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1716# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1717# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1718# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1719# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1720# OTHER DEALINGS IN THE SOFTWARE.