Coverage for pygeodesy/points.py: 95%
561 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
2# -*- coding: utf-8 -*-
4u'''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, map1, _Sequence, _xcopy, \
30 _xdup, _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, latDMS, lonDMS, parseDMS2, S_DEG, S_DMS, \
36 S_MIN, S_SEC, S_SEP
37from pygeodesy.errors import CrossError, crosserrors, _IndexError, \
38 _IsnotError, _TypeError, _ValueError, \
39 _xkwds, _xkwds_pop
40from pygeodesy.fmath import favg, fdot, Fsum, fsum, hypot
41# from pygeodesy.fsums import Fsum, fsum # from .fmath
42from pygeodesy.formy import _bearingTo2, equirectangular_, latlon2n_xyz, \
43 _spherical_datum
44from pygeodesy.interns import NN, _colinear_, _COMMASPACE_, _composite_, \
45 _DEQUALSPACED_, _ELLIPSIS_, _height_, \
46 _immutable_, _lat_, _lon_, _near_, _not_, \
47 _point_, _SPACE_, _UNDER_, _valid_
48from pygeodesy.iters import LatLon2PsxyIter, PointsIter, points2
49from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
50from pygeodesy.named import classname, nameof, notImplemented, notOverloaded, \
51 _NamedTuple, _xnamed, _xother3, _xotherError
52from pygeodesy.namedTuples import Bounds2Tuple, Bounds4Tuple, \
53 LatLon2Tuple, NearestOn3Tuple, \
54 NearestOn5Tuple, PhiLam2Tuple, \
55 Point3Tuple, Vector3Tuple, Vector4Tuple
56from pygeodesy.nvectorBase import NvectorBase, _N_vector_
57from pygeodesy.props import deprecated_method, Property_RO, property_doc_, \
58 property_RO
59from pygeodesy.streprs import Fmt, hstr, instr, pairs
60from pygeodesy.units import Number_, Radius, Scalar, Scalar_
61from pygeodesy.utily import atan2b, degrees90, degrees180, degrees2m, \
62 unroll180, unrollPI, wrap90, wrap180
64from math import cos, fabs, fmod, radians, sin
66__all__ = _ALL_LAZY.points
67__version__ = '23.03.30'
69_fin_ = 'fin'
70_ilat_ = 'ilat'
71_ilon_ = 'ilon'
72_ncols_ = 'ncols'
73_nrows_ = 'nrows'
76class LatLon_(object): # XXX in heights._HeightBase.height
77 '''Low-overhead C{LatLon} class for L{Numpy2LatLon} and L{Tuple2LatLon}.
78 '''
79 # __slots__ efficiency is voided if the __slots__ class attribute
80 # is used in a subclass of a class with the traditional __dict__,
81 # see <https://docs.Python.org/2/reference/datamodel.html#slots>
82 # and __slots__ must be repeated in sub-classes, see "Problems
83 # with __slots__" in Luciano Ramalho, "Fluent Python", page
84 # 276+, O'Reilly, 2016, also at <https://Books.Google.ie/
85 # books?id=bIZHCgAAQBAJ&lpg=PP1&dq=fluent%20python&pg=
86 # PT364#v=onepage&q=“Problems%20with%20__slots__”&f=false>
87 #
88 # __slots__ = (_lat_, _lon_, _height_, _datum_, _name_)
89 # Property_RO = property_RO # no __dict__ with __slots__!
90 #
91 # However, sys.getsizeof(LatLon_(1, 2)) is 72-88 with __slots__
92 # but only 48-64 bytes without in Python 2.7.18+ and Python 3+.
94 def __init__(self, lat, lon, name=NN, height=0, datum=None):
95 '''Creat a new, mininal, low-overhead L{LatLon_} instance,
96 without height and datum.
98 @arg lat: Latitude (C{degrees}).
99 @arg lon: Longitude (C{degrees}).
100 @kwarg name: Optional name (C{str}).
101 @kwarg height: Optional height (C{float} or C{int}).
102 @kwarg datum: Optional datum (L{Datum}, L{Ellipsoid},
103 L{Ellipsoid2}, L{a_f2Tuple} or I{scalar}
104 radius) or C{None}.
106 @raise TypeError: Invalid B{C{datum}}.
108 @note: The lat- and longitude are taken as-given,
109 un-clipped and un-validated .
110 '''
111 try: # most common use case
112 self.lat, self.lon = float(lat), float(lon) # Lat(lat), Lon(lon)
113 except (TypeError, ValueError):
114 self.lat, self.lon = parseDMS2(lat, lon, clipLat=0, clipLon=0) # PYCHOK LatLon2Tuple
115 self.name = str(name) if name else NN
116 self.height = height
117 self.datum = datum if datum is None else \
118 _spherical_datum(datum, name=self.name)
120 def __eq__(self, other):
121 return isinstance(other, LatLon_) and \
122 other.lat == self.lat and \
123 other.lon == self.lon
125 def __ne__(self, other):
126 return not self.__eq__(other)
128 def __repr__(self):
129 return self.toRepr()
131 def __str__(self):
132 return self.toStr()
134 def classof(self, *args, **kwds):
135 '''Instantiate this very class.
137 @arg args: Optional, positional arguments.
138 @kwarg kwds: Optional, keyword arguments.
140 @return: New instance (C{self.__class__}).
141 '''
142 return _xnamed(self.__class__(*args, **kwds), self.name)
144 def copy(self, deep=False):
145 '''Make a shallow or deep copy of this instance.
147 @kwarg deep: If C{True} make a deep, otherwise a
148 shallow copy (C{bool}).
150 @return: The copy (C{This} (sub-)class).
151 '''
152 return _xcopy(self, deep=deep)
154 def dup(self, **items):
155 '''Duplicate this instance, replacing some items.
157 @kwarg items: Attributes to be changed (C{any}).
159 @return: The duplicate (C{This} (sub-)class).
161 @raise AttributeError: Some B{C{items}} invalid.
162 '''
163 return _xdup(self, **items)
165 def heightStr(self, prec=-2):
166 '''Return a string for the height B{C{height}}.
168 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
170 @see: Function L{pygeodesy.hstr}.
171 '''
172 return hstr(self.height, prec=prec)
174 def intermediateTo(self, other, fraction, height=None, wrap=False):
175 '''Locate the point at a given fraction between (or along) this
176 and an other point.
178 @arg other: The other point (C{LatLon}).
179 @arg fraction: Fraction between both points (C{float},
180 0.0 for this and 1.0 for the other point).
181 @kwarg height: Optional height (C{meter}), overriding the
182 intermediate height.
183 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
185 @return: Intermediate point (this C{LatLon}).
187 @raise TypeError: Incompatible B{C{other}} C{type}.
188 '''
189 f = Scalar(fraction=fraction)
190 if isnear0(f):
191 r = self
192 elif isnear1(f) and not wrap:
193 r = self.others(other)
194 else:
195 r = self.others(other)
196 h = favg(self.height, r.height, f=f) if height is None else height
197 _, lon = unroll180(self.lon, r.lon, wrap=wrap)
198 r = self.classof(favg(self.lat, r.lat, f=f),
199 favg(self.lon, lon, f=f),
200 height=h, datum=self.datum,
201 name=self.intermediateTo.__name__)
202 return r
204 @Property_RO # PYCHOK no cover
205 def isEllipsoidal(self):
206 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).
207 '''
208 return self.datum.isEllipsoidal if self.datum else None
210 @Property_RO # PYCHOK no cover
211 def isEllipsoidalLatLon(self):
212 '''Get C{LatLon} base.
213 '''
214 return False
216 def isequalTo(self, other, eps=None):
217 '''Compare this point with an other point, I{ignoring} height.
219 @arg other: The other point (C{LatLon}).
220 @kwarg eps: Tolerance for equality (C{degrees}).
222 @return: C{True} if both points are identical,
223 I{ignoring} height, C{False} otherwise.
225 @raise UnitError: Invalid B{C{eps}}.
226 '''
227 self.others(other)
229 if eps:
230 return max(fabs(self.lat - other.lat),
231 fabs(self.lon - other.lon)) < Scalar_(eps=eps)
232 else:
233 return self.lat == other.lat and \
234 self.lon == other.lon
236 @Property_RO
237 def isSpherical(self): # PYCHOK no cover
238 '''Check whether this point is spherical (C{bool} or C{None} if unknown).
239 '''
240 return self.datum.isSpherical if self.datum else None
242 @Property_RO
243 def latlon(self):
244 '''Get the lat- and longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
245 '''
246 return LatLon2Tuple(self.lat, self.lon, name=self.name)
248 @Property_RO
249 def latlonheight(self):
250 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).
251 '''
252 return self.latlon.to3Tuple(self.height)
254 @Property_RO
255 def _N_vector(self):
256 '''(INTERNAL) Get the minimal, low-overhead (C{nvectorBase._N_vector_})
257 '''
258 return _N_vector_(*latlon2n_xyz(self.lat, self.lon),
259 h=self.height, name=self.name)
261 def others(self, *other, **name_other_up): # see .named._namedBase.others
262 '''Refined class comparison.
264 @arg other: The other instance (any C{type}).
265 @kwarg name_other_up: Overriding C{name=other} and C{up=1}
266 keyword arguments.
268 @return: The B{C{other}} if compatible.
270 @raise TypeError: Incompatible B{C{other}} C{type}.
271 '''
272 other, name, up = _xother3(self, other, **name_other_up)
273 if isinstance(other, self.__class__) or (hasattr(other, _lat_)
274 and hasattr(other, _lon_)):
275 return other
276 raise _xotherError(self, other, name=name, up=up + 1)
278 @Property_RO
279 def philam(self):
280 '''Get the lat- and longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
281 '''
282 return PhiLam2Tuple(radians(self.lat), radians(self.lon), name=self.name)
284 @Property_RO
285 def philamheight(self):
286 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
287 '''
288 return self.philam.to3Tuple(self.height)
290 @deprecated_method
291 def points(self, points, closed=False, base=None): # PYCHOK no cover
292 '''DEPRECATED, use method C{points2}.'''
293 return points2(points, closed=closed, base=base)
295 def points2(self, points, closed=False, base=None):
296 '''Check a path or polygon represented by points.
298 @arg points: The path or polygon points (C{LatLon}[])
299 @kwarg closed: Optionally, consider the polygon closed,
300 ignoring any duplicate or closing final
301 B{C{points}} (C{bool}).
302 @kwarg base: Optionally, check all B{C{points}} against
303 this base class, if C{None} don't check.
305 @return: A L{Points2Tuple}C{(number, points)} with the number
306 of points and the points C{list} or C{tuple}.
308 @raise PointsError: Insufficient number of B{C{points}}.
310 @raise TypeError: Some B{C{points}} are not B{C{base}}.
311 '''
312 return points2(points, closed=closed, base=base)
314 def PointsIter(self, points, loop=0, dedup=False):
315 '''Return a points iterator.
317 @arg points: The path or polygon points (C{LatLon}[])
318 @kwarg loop: Number of loop-back points (non-negative C{int}).
319 @kwarg dedup: Skip duplicate points (C{bool}).
321 @return: A new C{PointsIter} iterator.
323 @raise PointsError: Insufficient number of B{C{points}}.
324 '''
325 return PointsIter(points, loop=loop, base=self, dedup=dedup)
327 @deprecated_method
328 def to2ab(self): # PYCHOK no cover
329 '''DEPRECATED, use property L{philam}.'''
330 return self.philam
332 def toNvector(self, h=None, Nvector=NvectorBase, **Nvector_kwds):
333 '''Convert this point to C{n-vector} (normal to the earth's
334 surface) components, I{including height}.
336 @kwarg h: Optional height, overriding this point's height
337 (C{meter}).
338 @kwarg Nvector: Optional class to return the C{n-vector}
339 components (C{Nvector}) or C{None}.
340 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} keyword
341 arguments, ignored if C{B{Nvector} is None}.
343 @return: The C{n-vector} components B{C{Nvector}} or if
344 B{C{Nvector}} is C{None}, a L{Vector4Tuple}C{(x,
345 y, z, h)}.
347 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}
348 argument.
349 '''
350 x, y, z = latlon2n_xyz(self.lat, self.lon)
351 r = Vector4Tuple(x, y, z, self.height if h is None else h)
352 if Nvector is not None:
353 r = Nvector(x, y, z, **_xkwds(Nvector_kwds, h=r.h, ll=self))
354 return _xnamed(r, self.name)
356 def toRepr(self, **kwds):
357 '''This L{LatLon_} as a string "class(<degrees>, ...)".
359 @kwarg kwds: Optional, keyword arguments.
361 @return: Class instance (C{str}).
362 '''
363 _ = _xkwds_pop(kwds, std=None) # PYCHOK std unused
364 return Fmt.PAREN(classname(self), self.toStr(**kwds))
366 def toStr(self, form=F_D, joined=_COMMASPACE_, **prec_sep_s_D_M_S_kwds):
367 '''Convert this point to a "lat, lon[, height][, name][, ...]" string,
368 formatted in the given C{B{form}at}.
370 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see
371 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.
372 @kwarg joined: Separator to join the lat-, longitude, heigth, name
373 and other strings (C{str} or C{None} or C{NN} for
374 non-joined).
375 @kwarg prec_sep_s_D_M_S_kwds: Optional C{B{prec}ision}, C{B{sep}arator},
376 B{C{s_D}}, B{C{s_M}}, B{C{s_S}}, B{C{s_DMS}} and possibly
377 other keyword arguments, see functions L{pygeodesy.latDMS}
378 or L{pygeodesy.lonDMS}.
380 @return: This point in the specified C{B{form}at}, etc. (C{str} or
381 a 2- or 3+tuple C{(lat_str, lon_str[, height_str][, name_str][,
382 ...])} if C{B{joined}=NN} or C{B{joined}=None} and with the
383 C{height_str} and C{name_str} only included if non-zero
384 respectively non-empty).
386 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more
387 details about keyword arguments C{B{form}at}, C{B{prec}ision},
388 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.
389 '''
390 def _t3(prec=None, sep=S_SEP, s_D=S_DEG, s_M=S_MIN, s_S=S_SEC, s_DMS=S_DMS, **kwds):
391 return dict(prec=prec, sep=sep, s_D=s_D, s_M=s_M, s_S=s_S, s_DMS=s_DMS), kwds, prec
393 prec_sep_s_D_M_S, kwds, prec = _t3(**prec_sep_s_D_M_S_kwds)
394 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),
395 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))
396 if self.height:
397 t += (self.heightStr(),)
398 if self.name:
399 t += (repr(self.name),)
400 if kwds:
401 t += pairs(kwds) if prec is None else pairs(kwds, prec=prec)
402 return joined.join(t) if joined else t
404 @deprecated_method
405 def toStr2(self, **kwds): # PYCHOK no cover
406 '''DEPRECATED, used method L{toRepr}.'''
407 return self.toRepr(**kwds)
410def _isLatLon(inst):
411 '''(INTERNAL) Check a C{LatLon} or C{LatLon_} instance.
412 '''
413 return isinstance(inst, (LatLon_, _MODS.latlonBase.LatLonBase))
416def _isLatLon_(LL):
417 '''(INTERNAL) Check a (sub-)class of C{LatLon_}.
418 '''
419 return issubclassof(LL, LatLon_) or (isclass(LL) and
420 all(hasattr(LL, a) for a in _ALL_ATTRS_))
423# get all pseudo-slots for class C{LatLon_}
424_ALL_ATTRS_ = tuple(LatLon_(0, 0).__dict__.keys())
427class _Basequence(_Sequence): # immutable, on purpose
428 '''(INTERNAL) Base class.
429 '''
430 _array = []
431 _epsilon = EPS
432 _itemname = _point_
434 def _contains(self, point):
435 '''(INTERNAL) Check for a matching point.
436 '''
437 return any(self._findall(point, ()))
439 def copy(self, deep=False): # PYCHOK no cover
440 '''Make a shallow or deep copy of this instance.
442 @kwarg deep: If C{True} make a deep, otherwise a
443 shallow copy (C{bool}).
445 @return: The copy (C{This class} or subclass thereof).
446 '''
447 return _xcopy(self, deep=deep)
449 def _count(self, point):
450 '''(INTERNAL) Count the number of matching points.
451 '''
452 return sum(1 for _ in self._findall(point, ())) # NOT len()!
454 def dup(self, **items): # PYCHOK no cover
455 '''Duplicate this instance, I{without replacing items}.
457 @kwarg items: No attributes (I{not allowed}).
459 @return: The duplicate (C{This} (sub-)class).
461 @raise TypeError: Any B{C{items}} invalid.
462 '''
463 if items:
464 t = _SPACE_(classname(self), _immutable_)
465 raise _TypeError(txt=t, this=self, **items)
466 return _xdup(self)
468 @property_doc_(''' the equality tolerance (C{float}).''')
469 def epsilon(self):
470 '''Get the tolerance for equality tests (C{float}).
471 '''
472 return self._epsilon
474 @epsilon.setter # PYCHOK setter!
475 def epsilon(self, tol):
476 '''Set the tolerance for equality tests (C{scalar}).
478 @raise UnitError: Non-scalar or invalid B{C{tol}}.
479 '''
480 self._epsilon = Scalar_(tolerance=tol)
482 def _find(self, point, start_end):
483 '''(INTERNAL) Find the first matching point index.
484 '''
485 for i in self._findall(point, start_end):
486 return i
487 return -1
489 def _findall(self, point, start_end): # PYCHOK no cover
490 '''(INTERNAL) I{Must be implemented/overloaded}.
491 '''
492 notImplemented(self, point, start_end)
494 def _getitem(self, index):
495 '''(INTERNAL) Return point [index] or return a slice.
496 '''
497 # Luciano Ramalho, "Fluent Python", page 290+, O'Reilly, 2016
498 if isinstance(index, slice):
499 # XXX an numpy.[nd]array slice is a view, not a copy
500 return self.__class__(self._array[index], **self._slicekwds())
501 else:
502 return self.point(self._array[index])
504 def _index(self, point, start_end):
505 '''(INTERNAL) Find the first matching point index.
506 '''
507 for i in self._findall(point, start_end):
508 return i
509 raise _IndexError(self._itemname, point, txt=_not_('found'))
511 @property_RO
512 def isNumpy2(self): # PYCHOK no cover
513 '''Is this a Numpy2 wrapper?
514 '''
515 return False # isinstance(self, (Numpy2LatLon, ...))
517 @property_RO
518 def isPoints2(self): # PYCHOK no cover
519 '''Is this a LatLon2 wrapper/converter?
520 '''
521 return False # isinstance(self, (LatLon2psxy, ...))
523 @property_RO
524 def isTuple2(self): # PYCHOK no cover
525 '''Is this a Tuple2 wrapper?
526 '''
527 return False # isinstance(self, (Tuple2LatLon, ...))
529 def _iter(self):
530 '''(INTERNAL) Yield all points.
531 '''
532 _array, point = self._array, self.point
533 for i in range(len(self)):
534 yield point(_array[i])
536 def point(self, *attrs): # PYCHOK no cover
537 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded in}.
539 @arg attrs: Optional arguments.
540 '''
541 notOverloaded(self, *attrs)
543 def _range(self, start=None, end=None, step=1):
544 '''(INTERNAL) Return the range.
545 '''
546 if step > 0:
547 if start is None:
548 start = 0
549 if end is None:
550 end = len(self)
551 elif step < 0:
552 if start is None:
553 start = len(self) - 1
554 if end is None:
555 end = -1
556 else:
557 raise _ValueError(step=step)
558 return range(start, end, step)
560 def _repr(self):
561 '''(INTERNAL) Return a string representation.
562 '''
563 # XXX use Python 3+ reprlib.repr
564 t = repr(self._array[:1]) # only first row
565 t = _SPACE_(t[:-1], _ELLIPSIS_, Fmt.SQUARE(t[-1:], len(self)))
566 t = _SPACE_.join(t.split()) # coalesce spaces
567 return instr(self, t, **self._slicekwds())
569 def _reversed(self): # PYCHOK false
570 '''(INTERNAL) Yield all points in reverse order.
571 '''
572 _array, point = self._array, self.point
573 for i in range(len(self) - 1, -1, -1):
574 yield point(_array[i])
576 def _rfind(self, point, start_end):
577 '''(INTERNAL) Find the last matching point index.
578 '''
579 def _r3(start=None, end=None, step=-1):
580 return (start, end, step) # PYCHOK returns
582 for i in self._findall(point, _r3(*start_end)):
583 return i
584 return -1
586 def _slicekwds(self): # PYCHOK no cover
587 '''(INTERNAL) I{Should be overloaded}.
588 '''
589 return {}
592class _Array2LatLon(_Basequence): # immutable, on purpose
593 '''(INTERNAL) Base class for Numpy2LatLon or Tuple2LatLon.
594 '''
595 _array = ()
596 _ilat = 0 # row column index
597 _ilon = 0 # row column index
598 _LatLon = LatLon_ # default
599 _shape = ()
601 def __init__(self, array, ilat=0, ilon=1, LatLon=None, shape=()):
602 '''Handle a C{NumPy} or C{Tuple} array as a sequence of C{LatLon} points.
603 '''
604 ais = (_ilat_, ilat), (_ilon_, ilon)
606 if len(shape) != 2 or shape[0] < 1 or shape[1] < len(ais):
607 raise _IndexError('array.shape', shape)
609 self._array = array
610 self._shape = Shape2Tuple(shape) # *shape
612 if LatLon: # check the point class
613 if not _isLatLon_(LatLon):
614 raise _IsnotError(_valid_, LatLon=LatLon)
615 self._LatLon = LatLon
617 # check the attr indices
618 for n, (ai, i) in enumerate(ais):
619 if not isint(i):
620 raise _IsnotError(int.__name__, **{ai: i})
621 i = int(i)
622 if not 0 <= i < shape[1]:
623 raise _ValueError(ai, i)
624 for aj, j in ais[:n]:
625 if int(j) == i:
626 raise _ValueError(_DEQUALSPACED_(ai, aj, i))
627 setattr(self, NN(_UNDER_, ai), i)
629 def __contains__(self, latlon):
630 '''Check for a specific lat-/longitude.
632 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
633 C{(lat, lon)}).
635 @return: C{True} if B{C{latlon}} is present, C{False} otherwise.
637 @raise TypeError: Invalid B{C{latlon}}.
638 '''
639 return self._contains(latlon)
641 def __getitem__(self, index):
642 '''Return row[index] as C{LatLon} or return a L{Numpy2LatLon} slice.
643 '''
644 return self._getitem(index)
646 def __iter__(self):
647 '''Yield rows as C{LatLon}.
648 '''
649 return self._iter()
651 def __len__(self):
652 '''Return the number of rows.
653 '''
654 return self._shape[0]
656 def __repr__(self):
657 '''Return a string representation.
658 '''
659 return self._repr()
661 def __reversed__(self): # PYCHOK false
662 '''Yield rows as C{LatLon} in reverse order.
663 '''
664 return self._reversed()
666 __str__ = __repr__
668 def count(self, latlon):
669 '''Count the number of rows with a specific lat-/longitude.
671 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
672 C{(lat, lon)}).
674 @return: Count (C{int}).
676 @raise TypeError: Invalid B{C{latlon}}.
677 '''
678 return self._count(latlon)
680 def find(self, latlon, *start_end):
681 '''Find the first row with a specific lat-/longitude.
683 @arg latlon: Point (C{LatLon}) or 2-tuple (lat, lon).
684 @arg start_end: Optional C{[start[, end]]} index (integers).
686 @return: Index or -1 if not found (C{int}).
688 @raise TypeError: Invalid B{C{latlon}}.
689 '''
690 return self._find(latlon, start_end)
692 def _findall(self, latlon, start_end):
693 '''(INTERNAL) Yield indices of all matching rows.
694 '''
695 try:
696 lat, lon = latlon.lat, latlon.lon
697 except AttributeError:
698 try:
699 lat, lon = latlon
700 except (TypeError, ValueError):
701 raise _IsnotError(_valid_, latlon=latlon)
703 _ilat, _ilon = self._ilat, self._ilon
704 _array, _eps = self._array, self._epsilon
705 for i in self._range(*start_end):
706 row = _array[i]
707 if fabs(row[_ilat] - lat) <= _eps and \
708 fabs(row[_ilon] - lon) <= _eps:
709 yield i
711 def findall(self, latlon, *start_end):
712 '''Yield indices of all rows with a specific lat-/longitude.
714 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
715 C{(lat, lon)}).
716 @arg start_end: Optional C{[start[, end]]} index (C{int}).
718 @return: Indices (C{iterable}).
720 @raise TypeError: Invalid B{C{latlon}}.
721 '''
722 return self._findall(latlon, start_end)
724 def index(self, latlon, *start_end): # PYCHOK Python 2- issue
725 '''Find index of the first row with a specific lat-/longitude.
727 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
728 C{(lat, lon)}).
729 @arg start_end: Optional C{[start[, end]]} index (C{int}).
731 @return: Index (C{int}).
733 @raise IndexError: Point not found.
735 @raise TypeError: Invalid B{C{latlon}}.
736 '''
737 return self._index(latlon, start_end)
739 @Property_RO
740 def ilat(self):
741 '''Get the latitudes column index (C{int}).
742 '''
743 return self._ilat
745 @Property_RO
746 def ilon(self):
747 '''Get the longitudes column index (C{int}).
748 '''
749 return self._ilon
751# next = __iter__
753 def point(self, row): # PYCHOK *attrs
754 '''Instantiate a point C{LatLon}.
756 @arg row: Array row (numpy.array).
758 @return: Point (C{LatLon}).
759 '''
760 return self._LatLon(row[self._ilat], row[self._ilon])
762 def rfind(self, latlon, *start_end):
763 '''Find the last row with a specific lat-/longitude.
765 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
766 C{(lat, lon)}).
767 @arg start_end: Optional C{[start[, end]]} index (C{int}).
769 @note: Keyword order, first stop, then start.
771 @return: Index or -1 if not found (C{int}).
773 @raise TypeError: Invalid B{C{latlon}}.
774 '''
775 return self._rfind(latlon, start_end)
777 def _slicekwds(self):
778 '''(INTERNAL) Slice kwds.
779 '''
780 return dict(ilat=self._ilat, ilon=self._ilon)
782 @Property_RO
783 def shape(self):
784 '''Get the shape of the C{NumPy} array or the C{Tuples} as
785 L{Shape2Tuple}C{(nrows, ncols)}.
786 '''
787 return self._shape
789 def _subset(self, indices): # PYCHOK no cover
790 '''(INTERNAL) I{Must be implemented/overloaded}.
791 '''
792 notImplemented(self, indices)
794 def subset(self, indices):
795 '''Return a subset of the C{NumPy} array.
797 @arg indices: Row indices (C{range} or C{int}[]).
799 @note: A C{subset} is different from a C{slice} in 2 ways:
800 (a) the C{subset} is typically specified as a list of
801 (un-)ordered indices and (b) the C{subset} allocates
802 a new, separate C{NumPy} array while a C{slice} is
803 just an other C{view} of the original C{NumPy} array.
805 @return: Sub-array (C{numpy.array}).
807 @raise IndexError: Out-of-range B{C{indices}} value.
809 @raise TypeError: If B{C{indices}} is not a C{range}
810 nor an C{int}[].
811 '''
812 if not issequence(indices, tuple): # NO tuple, only list
813 # and range work properly to get Numpy array sub-sets
814 raise _IsnotError(_valid_, indices=type(indices))
816 n = len(self)
817 for i, v in enumerate(indices):
818 if not isint(v):
819 raise _TypeError(Fmt.SQUARE(indices=i), v)
820 elif not 0 <= v < n:
821 raise _IndexError(Fmt.SQUARE(indices=i), v)
823 return self._subset(indices)
826class LatLon2psxy(_Basequence):
827 '''Wrapper for C{LatLon} points as "on-the-fly" pseudo-xy coordinates.
828 '''
829 _closed = False
830 _len = 0
831 _deg2m = None # default, keep degrees
832 _radius = None
833 _wrap = True
835 def __init__(self, latlons, closed=False, radius=None, wrap=True):
836 '''Handle C{LatLon} points as pseudo-xy coordinates.
838 @note: The C{LatLon} latitude is considered the I{pseudo-y}
839 and longitude the I{pseudo-x} coordinate, likewise
840 for L{LatLon2Tuple}. However, 2-tuples C{(x, y)} are
841 considered as I{(longitude, latitude)}.
843 @arg latlons: Points C{list}, C{sequence}, C{set}, C{tuple},
844 etc. (C{LatLon[]}).
845 @kwarg closed: Optionally, close the polygon (C{bool}).
846 @kwarg radius: Mean earth radius (C{meter}).
847 @kwarg wrap: Wrap lat- and longitudes (C{bool}).
849 @raise PointsError: Insufficient number of B{C{latlons}}.
851 @raise TypeError: Some B{C{points}} are not B{C{base}}.
852 '''
853 self._closed = closed
854 self._len, self._array = points2(latlons, closed=closed)
855 if radius:
856 self._radius = r = Radius(radius)
857 self._deg2m = degrees2m(_1_0, r)
858 self._wrap = wrap
860 def __contains__(self, xy):
861 '''Check for a matching point.
863 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
864 C{(x, y)}) in (C{degrees}.
866 @return: C{True} if B{C{xy}} is present, C{False} otherwise.
868 @raise TypeError: Invalid B{C{xy}}.
869 '''
870 return self._contains(xy)
872 def __getitem__(self, index):
873 '''Return the pseudo-xy or return a L{LatLon2psxy} slice.
874 '''
875 return self._getitem(index)
877 def __iter__(self):
878 '''Yield all pseudo-xy's.
879 '''
880 return self._iter()
882 def __len__(self):
883 '''Return the number of pseudo-xy's.
884 '''
885 return self._len
887 def __repr__(self):
888 '''Return a string representation.
889 '''
890 return self._repr()
892 def __reversed__(self): # PYCHOK false
893 '''Yield all pseudo-xy's in reverse order.
894 '''
895 return self._reversed()
897 __str__ = __repr__
899 def count(self, xy):
900 '''Count the number of matching points.
902 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
903 C{(x, y)}) in (C{degrees}.
905 @return: Count (C{int}).
907 @raise TypeError: Invalid B{C{xy}}.
908 '''
909 return self._count(xy)
911 def find(self, xy, *start_end):
912 '''Find the first matching point.
914 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
915 C{(x, y)}) in (C{degrees}.
916 @arg start_end: Optional C{[start[, end]]} index (C{int}).
918 @return: Index or -1 if not found (C{int}).
920 @raise TypeError: Invalid B{C{xy}}.
921 '''
922 return self._find(xy, start_end)
924 def _findall(self, xy, start_end):
925 '''(INTERNAL) Yield indices of all matching points.
926 '''
927 try:
928 x, y = xy.lon, xy.lat
930 def _x_y_ll3(ll): # match LatLon
931 return ll.lon, ll.lat, ll
933 except AttributeError:
934 try:
935 x, y = xy[:2]
936 except (IndexError, TypeError, ValueError):
937 raise _IsnotError(_valid_, xy=xy)
939 _x_y_ll3 = self.point # PYCHOK expected
941 _array, _eps = self._array, self._epsilon
942 for i in self._range(*start_end):
943 xi, yi, _ = _x_y_ll3(_array[i])
944 if fabs(xi - x) <= _eps and \
945 fabs(yi - y) <= _eps:
946 yield i
948 def findall(self, xy, *start_end):
949 '''Yield indices of all matching points.
951 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
952 C{(x, y)}) in (C{degrees}.
953 @arg start_end: Optional C{[start[, end]]} index (C{int}).
955 @return: Indices (C{iterator}).
957 @raise TypeError: Invalid B{C{xy}}.
958 '''
959 return self._findall(xy, start_end)
961 def index(self, xy, *start_end): # PYCHOK Python 2- issue
962 '''Find the first matching point.
964 @arg xy: Point (C{LatLon}) or 2-tuple (x, y) in (C{degrees}).
965 @arg start_end: Optional C{[start[, end]]} index (C{int}).
967 @return: Index (C{int}).
969 @raise IndexError: Point not found.
971 @raise TypeError: Invalid B{C{xy}}.
972 '''
973 return self._index(xy, start_end)
975 @property_RO
976 def isPoints2(self):
977 '''Is this a LatLon2 wrapper/converter?
978 '''
979 return True # isinstance(self, (LatLon2psxy, ...))
981 def point(self, ll): # PYCHOK *attrs
982 '''Create a pseudo-xy.
984 @arg ll: Point (C{LatLon}).
986 @return: An L{Point3Tuple}C{(x, y, ll)}.
987 '''
988 x, y = ll.lon, ll.lat # note, x, y = lon, lat
989 if self._wrap:
990 x, y = wrap180(x), wrap90(y)
991 d = self._deg2m
992 if d: # convert degrees to meter (or radians)
993 x *= d
994 y *= d
995 return Point3Tuple(x, y, ll)
997 def rfind(self, xy, *start_end):
998 '''Find the last matching point.
1000 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
1001 C{(x, y)}) in (C{degrees}.
1002 @arg start_end: Optional C{[start[, end]]} index (C{int}).
1004 @return: Index or -1 if not found (C{int}).
1006 @raise TypeError: Invalid B{C{xy}}.
1007 '''
1008 return self._rfind(xy, start_end)
1010 def _slicekwds(self):
1011 '''(INTERNAL) Slice kwds.
1012 '''
1013 return dict(closed=self._closed, radius=self._radius, wrap=self._wrap)
1016class Numpy2LatLon(_Array2LatLon): # immutable, on purpose
1017 '''Wrapper for C{NumPy} arrays as "on-the-fly" C{LatLon} points.
1018 '''
1019 def __init__(self, array, ilat=0, ilon=1, LatLon=None):
1020 '''Handle a C{NumPy} array as a sequence of C{LatLon} points.
1022 @arg array: C{NumPy} array (C{numpy.array}).
1023 @kwarg ilat: Optional index of the latitudes column (C{int}).
1024 @kwarg ilon: Optional index of the longitudes column (C{int}).
1025 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
1027 @raise IndexError: If B{C{array.shape}} is not (1+, 2+).
1029 @raise TypeError: If B{C{array}} is not a C{NumPy} array or
1030 C{LatLon} is not a class with C{lat}
1031 and C{lon} attributes.
1033 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
1034 are the same or out of range.
1036 @example:
1038 >>> type(array)
1039 <type 'numpy.ndarray'> # <class ...> in Python 3+
1040 >>> points = Numpy2LatLon(array, lat=0, lon=1)
1041 >>> simply = simplifyRDP(points, ...)
1042 >>> type(simply)
1043 <type 'numpy.ndarray'> # <class ...> in Python 3+
1044 >>> sliced = points[1:-1]
1045 >>> type(sliced)
1046 <class '...Numpy2LatLon'>
1047 '''
1048 try: # get shape and check some other numpy.array attrs
1049 s, _, _ = array.shape, array.nbytes, array.ndim # PYCHOK expected
1050 except AttributeError:
1051 raise _IsnotError('NumPy', array=type(array))
1053 _Array2LatLon.__init__(self, array, ilat=ilat, ilon=ilon,
1054 LatLon=LatLon, shape=s)
1056 @property_RO
1057 def isNumpy2(self):
1058 '''Is this a Numpy2 wrapper?
1059 '''
1060 return True # isinstance(self, (Numpy2LatLon, ...))
1062 def _subset(self, indices):
1063 return self._array[indices] # NumPy special
1066class Shape2Tuple(_NamedTuple):
1067 '''2-Tuple C{(nrows, ncols)}, the number of rows and columns,
1068 both C{int}.
1069 '''
1070 _Names_ = (_nrows_, _ncols_)
1071 _Units_ = ( Number_, Number_)
1074class Tuple2LatLon(_Array2LatLon):
1075 '''Wrapper for tuple sequences as "on-the-fly" C{LatLon} points.
1076 '''
1077 def __init__(self, tuples, ilat=0, ilon=1, LatLon=None):
1078 '''Handle a list of tuples, each containing a lat- and longitude
1079 and perhaps other values as a sequence of C{LatLon} points.
1081 @arg tuples: The C{list}, C{tuple} or C{sequence} of tuples (C{tuple}[]).
1082 @kwarg ilat: Optional index of the latitudes value (C{int}).
1083 @kwarg ilon: Optional index of the longitudes value (C{int}).
1084 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
1086 @raise IndexError: If C{(len(B{tuples}), min(len(t) for t
1087 in B{tuples}))} is not (1+, 2+).
1089 @raise TypeError: If B{C{tuples}} is not a C{list}, C{tuple}
1090 or C{sequence} or if B{C{LatLon}} is not a
1091 C{LatLon} with C{lat}, C{lon} and C{name}
1092 attributes.
1094 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
1095 are the same or out of range.
1097 @example:
1099 >>> tuples = [(0, 1), (2, 3), (4, 5)]
1100 >>> type(tuples)
1101 <type 'list'> # <class ...> in Python 3+
1102 >>> points = Tuple2LatLon(tuples, lat=0, lon=1)
1103 >>> simply = simplifyRW(points, 0.5, ...)
1104 >>> type(simply)
1105 <type 'list'> # <class ...> in Python 3+
1106 >>> simply
1107 [(0, 1), (4, 5)]
1108 >>> sliced = points[1:-1]
1109 >>> type(sliced)
1110 <class '...Tuple2LatLon'>
1111 >>> sliced
1112 ...Tuple2LatLon([(2, 3), ...][1], ilat=0, ilon=1)
1114 >>> closest, _ = nearestOn2(LatLon_(2, 1), points, adjust=False)
1115 >>> closest
1116 LatLon_(lat=1.0, lon=2.0)
1118 >>> closest, _ = nearestOn2(LatLon_(3, 2), points)
1119 >>> closest
1120 LatLon_(lat=2.001162, lon=3.001162)
1121 '''
1122 _xinstanceof(list, tuple, tuples=tuples)
1123 s = len(tuples), min(len(_) for _ in tuples)
1124 _Array2LatLon.__init__(self, tuples, ilat=ilat, ilon=ilon,
1125 LatLon=LatLon, shape=s)
1127 @property_RO
1128 def isTuple2(self):
1129 '''Is this a Tuple2 wrapper?
1130 '''
1131 return True # isinstance(self, (Tuple2LatLon, ...))
1133 def _subset(self, indices):
1134 return type(self._array)(self._array[i] for i in indices)
1137def _area2(points, adjust, wrap):
1138 '''(INTERNAL) Approximate the area in radians squared, I{signed}.
1139 '''
1140 if adjust:
1141 # approximate trapezoid by a rectangle, adjusting
1142 # the top width by the cosine of the latitudinal
1143 # average and bottom width by some fudge factor
1144 def _adjust(w, h):
1145 c = cos(h) if fabs(h) < PI_2 else _0_0
1146 return w * h * (c + 1.2876) * _0_5
1147 else:
1148 def _adjust(w, h): # PYCHOK expected
1149 return w * h
1151 # setting radius=1 converts degrees to radians
1152 Ps = LatLon2PsxyIter(points, wrap=wrap, radius=_1_0, loop=1)
1153 x1, y1, ll = Ps[0]
1154 pts = [ll]
1156 A2 = Fsum() # trapezoidal area in radians**2
1157 for i, p in Ps.enumerate(closed=True):
1158 x2, y2, ll = p
1159 if 0 < i < 4:
1160 pts.append(ll)
1161 w, x2 = unrollPI(x1, x2, wrap=wrap if i else False)
1162 A2 += _adjust(w, (y2 + y1) * _0_5)
1163 x1, y1 = x2, y2
1165 return A2.fsum(), tuple(pts)
1168def _areaError(pts, near_=NN): # imported by .ellipsoidalKarney
1169 '''(INTERNAL) Area issue.
1170 '''
1171 t = _ELLIPSIS_(pts[:3], NN)
1172 return _ValueError(NN(near_, 'zero or polar area'), txt=t)
1175def areaOf(points, adjust=True, radius=R_M, wrap=True):
1176 '''Approximate the area of a polygon or composite.
1178 @arg points: The polygon points or clips (C{LatLon}[],
1179 L{BooleanFHP} or L{BooleanGH}).
1180 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1181 by the cosine of the mean latitude (C{bool}).
1182 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1183 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1185 @return: Approximate area (I{square} C{meter}, same units as
1186 B{C{radius}} or C{radians} I{squared} if B{C{radius}}
1187 is C{None}).
1189 @raise PointsError: Insufficient number of B{C{points}}
1191 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1193 @raise ValueError: Invalid B{C{radius}}.
1195 @note: This area approximation has limited accuracy and is
1196 ill-suited for regions exceeding several hundred Km
1197 or Miles or with near-polar latitudes.
1199 @see: L{sphericalNvector.areaOf}, L{sphericalTrigonometry.areaOf},
1200 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.
1201 '''
1202 if _MODS.booleans.isBoolean(points):
1203 a = points._sum1(areaOf, adjust=adjust, radius=None, wrap=wrap)
1204 else:
1205 a, _ = _area2(points, adjust, wrap)
1206 return fabs(a if radius is None else (Radius(radius)**2 * a))
1209def boundsOf(points, wrap=True, LatLon=None):
1210 '''Determine the bottom-left SW and top-right NE corners of a
1211 path or polygon.
1213 @arg points: The path or polygon points (C{LatLon}[]).
1214 @kwarg wrap: Wrap lat- and longitudes (C{bool}).
1215 @kwarg LatLon: Optional class to return the C{bounds}
1216 corners (C{LatLon}) or C{None}.
1218 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)} as
1219 B{C{LatLon}}s if B{C{LatLon}} is C{None} a
1220 L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}.
1222 @raise PointsError: Insufficient number of B{C{points}}
1224 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1226 @see: Function L{quadOf}.
1228 @example:
1230 >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1231 >>> boundsOf(b) # False
1232 >>> 45.0, 1.0, 46.0, 2.0
1233 '''
1234 Ps = LatLon2PsxyIter(points, wrap=wrap, loop=1)
1235 w, s, _ = e, n, _ = Ps[0]
1237 for x, y, _ in Ps.iterate(closed=False): # [1:]
1238 if w > x:
1239 w = x
1240 elif e < x:
1241 e = x
1243 if s > y:
1244 s = y
1245 elif n < y:
1246 n = y
1248 return Bounds4Tuple(s, w, n, e) if LatLon is None else \
1249 Bounds2Tuple(LatLon(s, w), LatLon(n, e)) # PYCHOK inconsistent
1252def centroidOf(points, wrap=True, LatLon=None):
1253 '''Determine the centroid of a polygon.
1255 @arg points: The polygon points (C{LatLon}[]).
1256 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1257 @kwarg LatLon: Optional class to return the centroid (C{LatLon})
1258 or C{None}.
1260 @return: Centroid (B{C{LatLon}}) or a L{LatLon2Tuple}C{(lat, lon)}
1261 if C{B{LatLon} is None}.
1263 @raise PointsError: Insufficient number of B{C{points}}
1265 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1267 @raise ValueError: The B{C{points}} enclose a pole or
1268 near-zero area.
1270 @see: U{Centroid<https://WikiPedia.org/wiki/Centroid#Of_a_polygon>} and
1271 Paul Bourke's U{Calculating The Area And Centroid Of A Polygon
1272 <https://www.SEAS.UPenn.edu/~ese502/lab-content/extra_materials/
1273 Polygon%20Area%20and%20Centroid.pdf>}, 1988.
1274 '''
1275 A, X, Y = Fsum(), Fsum(), Fsum()
1277 # setting radius=1 converts degrees to radians
1278 Ps = LatLon2PsxyIter(points, wrap=wrap, radius=_1_0, loop=1)
1279 x1, y1, ll = Ps[0]
1280 pts = [ll] # for _areaError
1281 for i, p in Ps.enumerate(closed=True):
1282 x2, y2, ll = p
1283 if 0 < i < 4:
1284 pts.append(ll)
1285 if wrap and i != 0:
1286 _, x2 = unrollPI(x1, x2, wrap=wrap)
1287 t = x1 * y2 - x2 * y1
1288 A += t
1289 X += t * (x1 + x2)
1290 Y += t * (y1 + y2)
1291 # XXX more elaborately:
1292 # t1, t2 = x1 * y2, -(x2 * y1)
1293 # A.fadd_(t1, t2)
1294 # X.fadd_(t1 * x1, t1 * x2, t2 * x1, t2 * x2)
1295 # Y.fadd_(t1 * y1, t1 * y2, t2 * y1, t2 * y2)
1296 x1, y1 = x2, y2
1298 a = A.fmul(_6_0).fover(_2_0)
1299 if isnear0(a):
1300 raise _areaError(pts, near_=_near_)
1301 y, x = degrees90(Y.fover(a)), degrees180(X.fover(a))
1302 return LatLon2Tuple(y, x) if LatLon is None else LatLon(y, x)
1305def fractional(points, fi, j=None, wrap=None, LatLon=None, Vector=None, **kwds):
1306 '''Return the point at a given I{fractional} index.
1308 @arg points: The points (C{LatLon}[], L{Numpy2LatLon}[],
1309 L{Tuple2LatLon}[], C{Cartesian}[], C{Vector3d}[],
1310 L{Vector3Tuple}[]).
1311 @arg fi: The fractional index (L{FIx}, C{float} or C{int}).
1312 @kwarg j: Optionally, index of the other point (C{int}).
1313 @kwarg wrap: Wrap and unroll longitudes (C{bool}) or C{None} for
1314 backward compatible L{LatLon2Tuple} or B{C{LatLon}}
1315 with averaged lat- and longitudes. Use C{True} or
1316 C{False} to get the I{fractional} point computed
1317 method C{points[fi].intermediateTo}.
1318 @kwarg LatLon: Optional class to return the I{intermediate},
1319 I{fractional} point (C{LatLon}) or C{None}.
1320 @kwarg Vector: Optional class to return the I{intermediate},
1321 I{fractional} point (C{Cartesian}, C{Vector3d})
1322 or C{None}.
1323 @kwarg kwds: Optional, additional B{C{LatLon}} I{or} B{C{Vector}}
1324 keyword arguments, ignored if both C{B{LatLon}} and
1325 C{B{Vector}} are C{None}.
1327 @return: A L{LatLon2Tuple}C{(lat, lon)} if B{C{wrap}}, B{C{LatLon}}
1328 and B{C{Vector}} all are C{None}, the defaults.
1330 An instance of B{C{LatLon}} if not C{None} I{or} an instance
1331 of B{C{Vector}} if not C{None}.
1333 Otherwise with B{C{wrap}} either C{True} or C{False} and
1334 B{C{LatLon}} and B{C{Vector}} both C{None}, an instance of
1335 B{C{points}}' (sub-)class C{intermediateTo} I{fractional}.
1337 Summarized as follows:
1339 >>> wrap | LatLon | Vector | returned type/value
1340 # -------+--------+--------+--------------+------
1341 # | | | LatLon2Tuple | favg
1342 # None | None | None | or** |
1343 # | | | Vector3Tuple | favg
1344 # None | LatLon | None | LatLon | favg
1345 # None | None | Vector | Vector | favg
1346 # -------+--------+--------+--------------+------
1347 # True | None | None | points' | .iTo
1348 # True | LatLon | None | LatLon | .iTo
1349 # True | None | Vector | Vector | .iTo
1350 # -------+--------+--------+--------------+------
1351 # False | None | None | points' | .iTo
1352 # False | LatLon | None | LatLon | .iTo
1353 # False | None | Vector | Vector | .iTo
1354 # _____
1355 # favg) averaged lat, lon or x, y, z values
1356 # .iTo) value from points[fi].intermediateTo
1357 # **) depends on base class of points[fi]
1359 @raise IndexError: Fractional index B{C{fi}} invalid or B{C{points}}
1360 not subscriptable or not closed.
1362 @raise TypeError: Invalid B{C{LatLon}}, B{C{Vector}} or B{C{kwds}}
1363 argument.
1365 @see: Class L{FIx} and method L{FIx.fractional}.
1366 '''
1367 if LatLon and Vector: # PYCHOK no cover
1368 kwds = _xkwds(kwds, fi=fi, LatLon=LatLon, Vector=Vector)
1369 raise _TypeError(txt=fractional.__name__, **kwds)
1370 try:
1371 if not isscalar(fi) or fi < 0:
1372 raise IndexError
1373 n = getattr(fi, _fin_, 0)
1374 w = wrap if Vector is None else False # intermediateTo
1375 p = _fractional(points, fi, j, fin=n, wrap=w) # see .units.FIx
1376 if LatLon:
1377 p = LatLon(p.lat, p.lon, **kwds)
1378 elif Vector:
1379 p = Vector(p.x, p.y, p.z, **kwds)
1380 except (IndexError, TypeError):
1381 raise _IndexError(fi=fi, points=points, wrap=wrap, txt=fractional.__name__)
1382 return p
1385def _fractional(points, fi, j, fin=None, wrap=None): # in .frechet.py
1386 '''(INTERNAL) Compute point at L{fractional} index C{fi} and C{j}.
1387 '''
1388 i = int(fi)
1389 p = points[i]
1390 r = fi - float(i)
1391 if r > EPS: # EPS0?
1392 if j is None:
1393 j = i + 1
1394 if fin:
1395 j %= fin
1396 q = points[j]
1397 if r >= EPS1: # PYCHOK no cover
1398 p = q
1399 elif wrap is not None: # in (True, False)
1400 p = p.intermediateTo(q, r, wrap=wrap)
1401 elif _isLatLon(p): # backward compatible default
1402 p = LatLon2Tuple(favg(p.lat, q.lat, f=r),
1403 favg(p.lon, q.lon, f=r),
1404 name=fractional.__name__)
1405 else: # assume p and q are cartesian or vectorial
1406 z = p.z if p.z is q.z else favg(p.z, q.z, f=r)
1407 p = Vector3Tuple(favg(p.x, q.x, f=r),
1408 favg(p.y, q.y, f=r), z,
1409 name=fractional.__name__)
1410 return p
1413def isclockwise(points, adjust=False, wrap=True):
1414 '''Determine the direction of a path or polygon.
1416 @arg points: The path or polygon points (C{LatLon}[]).
1417 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1418 by the cosine of the mean latitude (C{bool}).
1419 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1421 @return: C{True} if B{C{points}} are clockwise, C{False} otherwise.
1423 @raise PointsError: Insufficient number of B{C{points}}
1425 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1427 @raise ValueError: The B{C{points}} enclose a pole or zero area.
1429 @example:
1431 >>> f = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1432 >>> isclockwise(f) # False
1433 >>> isclockwise(reversed(f)) # True
1434 '''
1435 a, pts = _area2(points, adjust, wrap)
1436 if a > 0: # opposite of ellipsoidalExact and -Karney
1437 return True
1438 elif a < 0:
1439 return False
1440 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1441 raise _areaError(pts)
1444def isconvex(points, adjust=False, wrap=True):
1445 '''Determine whether a polygon is convex.
1447 @arg points: The polygon points (C{LatLon}[]).
1448 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1449 by the cosine of the mean latitude (C{bool}).
1450 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1452 @return: C{True} if B{C{points}} are convex, C{False} otherwise.
1454 @raise CrossError: Some B{C{points}} are colinear.
1456 @raise PointsError: Insufficient number of B{C{points}}
1458 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1460 @example:
1462 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1463 >>> isconvex(t) # True
1465 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1466 >>> isconvex(f) # False
1467 '''
1468 return bool(isconvex_(points, adjust=adjust, wrap=wrap))
1471def isconvex_(points, adjust=False, wrap=True):
1472 '''Determine whether a polygon is convex I{and clockwise}.
1474 @arg points: The polygon points (C{LatLon}[]).
1475 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1476 by the cosine of the mean latitude (C{bool}).
1477 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1479 @return: C{+1} if B{C{points}} are convex clockwise, C{-1} for
1480 convex counter-clockwise B{C{points}}, C{0} otherwise.
1482 @raise CrossError: Some B{C{points}} are colinear.
1484 @raise PointsError: Insufficient number of B{C{points}}
1486 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1488 @example:
1490 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1491 >>> isconvex_(t) # +1
1493 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1494 >>> isconvex_(f) # 0
1495 '''
1496 def _unroll2(x1, y1, x2, y2, a, w):
1497 x21, x2 = unroll180(x1, x2, wrap=w)
1498 if a:
1499 y = radians(y1 + y2) * _0_5
1500 x21 *= cos(y) if fabs(y) < PI_2 else _0_0
1501 return x21, x2
1503 c, s = crosserrors(), 0
1505 Ps = LatLon2PsxyIter(points, wrap=wrap, loop=2)
1506 x1, y1, _ = Ps[0]
1507 x2, y2, _ = Ps[1]
1508 x21, x2 = _unroll2(x1, y1, x2, y2, adjust, False)
1510 for i, p in Ps.enumerate(closed=True):
1511 x3, y3, ll = p
1512 x32, x3 = _unroll2(x2, y2, x3, y3, adjust, (wrap if i > 1 else False))
1514 # get the sign of the distance from point
1515 # x3, y3 to the line from x1, y1 to x2, y2
1516 # <https://WikiPedia.org/wiki/Distance_from_a_point_to_a_line>
1517 s3 = fdot((x3, y3, x1, y1), y2 - y1, -x21, -y2, x2)
1518 if s3 > 0: # x3, y3 on the right
1519 if s < 0: # non-convex
1520 return 0
1521 s = +1
1523 elif s3 < 0: # x3, y3 on the left
1524 if s > 0: # non-convex
1525 return 0
1526 s = -1
1528 elif c and fdot((x32, y1 - y2), y3 - y2, -x21) < 0: # PYCHOK no cover
1529 # colinear u-turn: x3, y3 not on the
1530 # opposite side of x2, y2 as x1, y1
1531 t = Fmt.SQUARE(points=i)
1532 raise CrossError(t, ll, txt=_colinear_)
1534 x1, y1, x2, y2, x21 = x2, y2, x3, y3, x32
1536 return s # all points on the same side
1539def isenclosedBy(point, points, wrap=False): # MCCABE 15
1540 '''Determine whether a point is enclosed by a polygon or composite.
1542 @arg point: The point (C{LatLon} or 2-tuple C{(lat, lon)}).
1543 @arg points: The polygon points or clips (C{LatLon}[], L{BooleanFHP}
1544 or L{BooleanGH}).
1545 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1547 @return: C{True} if the B{C{point}} is inside the polygon or
1548 composite, C{False} otherwise.
1550 @raise PointsError: Insufficient number of B{C{points}}
1552 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1554 @raise ValueError: Invalid B{C{point}}, lat- or longitude.
1556 @see: Functions L{pygeodesy.isconvex} and L{pygeodesy.ispolar} especially
1557 if the B{C{points}} may enclose a pole or wrap around the earth
1558 I{longitudinally}, methods L{sphericalNvector.LatLon.isenclosedBy},
1559 L{sphericalTrigonometry.LatLon.isenclosedBy} and U{MultiDop
1560 GeogContainPt<https://GitHub.com/NASA/MultiDop>} (U{Shapiro et.al. 2009,
1561 JTECH<https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}
1562 and U{Potvin et al. 2012, JTECH <https://Journals.AMetSoc.org/doi/abs/
1563 10.1175/JTECH-D-11-00019.1>}).
1564 '''
1565 try:
1566 y0, x0 = point.lat, point.lon
1567 except AttributeError:
1568 try:
1569 y0, x0 = map1(float, *point[:2])
1570 except (IndexError, TypeError, ValueError) as x:
1571 raise _ValueError(point=point, cause=x)
1573 if wrap:
1574 x0, y0 = wrap180(x0), wrap90(y0)
1576 def _dxy3(x1, x2, y2, w):
1577 dx, x2 = unroll180(x1, x2, wrap=w)
1578 return dx, x2, y2
1580 else:
1581 x0 = fmod(x0, _360_0) # not x0 % 360!
1582 x0_180_ = x0 - _180_0
1583 x0_180 = x0 + _180_0
1585 def _dxy3(x1, x, y, unused): # PYCHOK expected
1586 x = _umod_360(float(x))
1587 if x < x0_180_:
1588 x += _360_0
1589 elif x >= x0_180:
1590 x -= _360_0
1591 return (x - x1), x, y
1593 if _MODS.booleans.isBoolean(points):
1594 return points._encloses(y0, x0)
1596 Ps = LatLon2PsxyIter(points, wrap=wrap, loop=1)
1597 p = Ps[0]
1598 e = m = False
1599 S = Fsum()
1601 _, x1, y1 = _dxy3(x0, p.x, p.y, False)
1602 for i, p in Ps.enumerate(closed=True):
1603 dx, x2, y2 = _dxy3(x1, p.x, p.y, (wrap if i else False))
1604 # ignore duplicate and near-duplicate pts
1605 if max(fabs(dx), fabs(y2 - y1)) > EPS:
1606 # determine if polygon edge (x1, y1)..(x2, y2) straddles
1607 # point (lat, lon) or is on boundary, but do not count
1608 # edges on boundary as more than one crossing
1609 if fabs(dx) < 180 and (x1 < x0 <= x2 or x2 < x0 <= x1):
1610 m = not m
1611 dy = (x0 - x1) * (y2 - y1) - (y0 - y1) * dx
1612 if (dy > 0 and dx >= 0) or (dy < 0 and dx <= 0):
1613 e = not e
1615 S += sin(radians(y2))
1616 x1, y1 = x2, y2
1618 # An odd number of meridian crossings means, the polygon
1619 # contains a pole. Assume it is the pole on the hemisphere
1620 # containing the polygon mean point and if the polygon does
1621 # contain the North Pole, flip the result.
1622 if m and S.fsum() > 0:
1623 e = not e
1624 return e
1627def ispolar(points, wrap=False):
1628 '''Check whether a polygon encloses a pole.
1630 @arg points: The polygon points (C{LatLon}[]).
1631 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1633 @return: C{True} if the polygon encloses a pole, C{False}
1634 otherwise.
1636 @raise PointsError: Insufficient number of B{C{points}}
1638 @raise TypeError: Some B{C{points}} are not C{LatLon} or don't
1639 have C{bearingTo2}, C{initialBearingTo}
1640 and C{finalBearingTo} methods.
1641 '''
1642 def _cds(points, wrap): # iterate over course deltas
1643 Ps = PointsIter(points, loop=2)
1644 p2, p1 = Ps[0:2]
1645 b1, _ = _bearingTo2(p2, p1, wrap=wrap)
1647 for p2 in Ps.iterate(closed=True):
1648 if not p2.isequalTo(p1, EPS):
1649 b, b2 = _bearingTo2(p1, p2, wrap=wrap)
1650 yield wrap180(b - b1) # (b - b1 + 540) % 360 - 180
1651 yield wrap180(b2 - b) # (b2 - b + 540) % 360 - 180
1652 p1, b1 = p2, b2
1654 # summation of course deltas around pole is 0° rather than normally ±360°
1655 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1656 s = fsum(_cds(points, wrap))
1657 # XXX fix (intermittant) edge crossing pole - eg (85,90), (85,0), (85,-90)
1658 return fabs(s) < 90 # "zero-ish"
1661def luneOf(lon1, lon2, closed=False, LatLon=LatLon_, **LatLon_kwds):
1662 '''Generate an ellipsoidal or spherical U{lune
1663 <https://WikiPedia.org/wiki/Spherical_lune>}-shaped path or polygon.
1665 @arg lon1: Left longitude (C{degrees90}).
1666 @arg lon2: Right longitude (C{degrees90}).
1667 @kwarg closed: Optionally, close the path (C{bool}).
1668 @kwarg LatLon: Class to use (L{LatLon_}).
1669 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1670 keyword arguments.
1672 @return: A tuple of 4 or 5 B{C{LatLon}} instances outlining
1673 the lune shape.
1675 @see: U{Latitude-longitude quadrangle
1676 <https://www.MathWorks.com/help/map/ref/areaquad.html>}.
1677 '''
1678 t = (LatLon( _0_0, lon1, **LatLon_kwds),
1679 LatLon( _90_0, lon1, **LatLon_kwds),
1680 LatLon( _0_0, lon2, **LatLon_kwds),
1681 LatLon(_N_90_0, lon2, **LatLon_kwds))
1682 if closed:
1683 t += t[:1]
1684 return t
1687def nearestOn5(point, points, closed=False, wrap=False, LatLon=None, **options):
1688 '''Locate the point on a path or polygon closest to a reference point.
1690 The closest point is either on and within the extent of a polygon edge or
1691 the nearest of that edge's end points.
1693 @arg point: Reference point (C{LatLon}).
1694 @arg points: The path or polygon points (C{LatLon}[]).
1695 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1696 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes and longitudinal
1697 delta (C{bool}) in function L{pygeodesy.equirectangular_}.
1698 @kwarg LatLon: Optional class to return the closest point (C{LatLon})
1699 or C{None}.
1700 @kwarg options: Other keyword arguments for function
1701 L{pygeodesy.equirectangular_}.
1703 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1704 {closest} point (B{C{LatLon}}) or if C{B{LatLon} is None},
1705 a L{NearestOn5Tuple}C{(lat, lon, distance, angle, height)}.
1706 The C{distance} is the L{pygeodesy.equirectangular} distance
1707 between the C{closest} and reference B{C{point}} in C{degrees}.
1708 The C{angle} from the reference B{C{point}} to the C{closest}
1709 is in compass C{degrees360}, like function
1710 L{pygeodesy.compassAngle}.
1712 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1713 see function L{pygeodesy.equirectangular_}.
1715 @raise PointsError: Insufficient number of B{C{points}}
1717 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1719 @note: Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1720 subject to the supplied B{C{options}}. Method C{LatLon.nearestOn6}
1721 measures distances more accurately.
1723 @see: Function L{pygeodesy.nearestOn6} for cartesian points. Use function
1724 L{pygeodesy.degrees2m} to convert C{degrees} to C{meter}.
1725 '''
1726 def _d2yx(p2, p1, u, w):
1727 # w = wrap if (not closed or w < (n - 1)) else False
1728 # equirectangular_ returns a Distance4Tuple(distance
1729 # in degrees squared, delta lat, delta lon, p2.lon
1730 # unroll/wrap); the previous p2.lon unroll/wrap
1731 # is also applied to the next edge's p1.lon
1732 return equirectangular_(p1.lat, p1.lon + u,
1733 p2.lat, p2.lon, wrap=w, **options)
1735 def _h(p): # get height or default 0
1736 return getattr(p, _height_, 0) or 0
1738 # 3-D version used in .vector3d._nearestOn2
1739 #
1740 # point (x, y) on axis rotated ccw by angle a:
1741 # x' = y * sin(a) + x * cos(a)
1742 # y' = y * cos(a) - x * sin(a)
1743 #
1744 # distance (w) along and (h) perpendicular to
1745 # a line thru point (dx, dy) and the origin:
1746 # w = (y * dy + x * dx) / hypot(dx, dy)
1747 # h = (y * dx - x * dy) / hypot(dx, dy)
1748 #
1749 # closest point on that line thru (dx, dy):
1750 # xc = dx * w / hypot(dx, dy)
1751 # yc = dy * w / hypot(dx, dy)
1752 # or
1753 # xc = dx * f
1754 # yc = dy * f
1755 # with
1756 # f = w / hypot(dx, dy)
1757 # or
1758 # f = (y * dy + x * dx) / hypot2(dx, dy)
1759 #
1760 # i.e. no need for sqrt or hypot
1762 Ps = PointsIter(points, loop=1)
1763 p1 = c = Ps[0]
1764 u1 = u = _0_0
1765 d, dy, dx, _ = _d2yx(p1, point, u1, False)
1766 for i, p2 in Ps.enumerate(closed=closed):
1767 # iff wrapped, unroll lon1 (actually previous
1768 # lon2) like function unroll180/-PI would've
1769 w = False if closed and i == 0 else wrap
1770 d21, y21, x21, u2 = _d2yx(p2, p1, u1, w)
1771 if d21 > EPS:
1772 # distance point to p1, y01 and x01 inverted
1773 d2, y01, x01, _ = _d2yx(point, p1, u1, closed)
1774 if d2 > EPS:
1775 w2 = y01 * y21 + x01 * x21
1776 if w2 > 0:
1777 if w2 < d21:
1778 # closest is between p1 and p2, use
1779 # original delta's, not y21 and x21
1780 f = w2 / d21
1781 p1 = LatLon_(favg(p1.lat, p2.lat, f=f),
1782 favg(p1.lon, p2.lon + u2, f=f),
1783 height=favg(_h(p1), _h(p2), f=f))
1784 u1 = _0_0
1785 else: # p2 is closest
1786 p1, u1 = p2, u2
1787 d2, y01, x01, _ = _d2yx(point, p1, u1, closed)
1788 if d2 < d: # p1 is closer, y01 and x01 negated
1789 c, u, d, dy, dx = p1, u1, d2, -y01, -x01
1790 p1, u1 = p2, u2
1792 d, a, h = hypot(dx, dy), atan2b(dx, dy), _h(c)
1793 if LatLon is None:
1794 r = NearestOn5Tuple(c.lat, c.lon + u, d, a, h)
1795 else:
1796 r = LatLon(c.lat, c.lon + u, height=h)
1797 r = NearestOn3Tuple(r, d, a)
1798 return _xnamed(r, nameof(point))
1801def perimeterOf(points, closed=False, adjust=True, radius=R_M, wrap=True):
1802 '''Approximate the perimeter of a path, polygon. or composite
1804 @arg points: The path or polygon points or clips (C{LatLon}[],
1805 L{BooleanFHP} or L{BooleanGH}).
1806 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1807 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1808 by the cosine of the mean latitude (C{bool}).
1809 @kwarg radius: Mean earth radius (C{meter}).
1810 @kwarg wrap: Wrap lat-, wrap and unroll longitudes (C{bool}).
1812 @return: Approximate perimeter (C{meter}, same units as
1813 B{C{radius}}).
1815 @raise PointsError: Insufficient number of B{C{points}}
1817 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1819 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1820 C{B{points}} a composite.
1822 @note: This perimeter is based on the L{pygeodesy.equirectangular_}
1823 distance approximation and is ill-suited for regions exceeding
1824 several hundred Km or Miles or with near-polar latitudes.
1826 @see: Functions L{sphericalTrigonometry.perimeterOf} and
1827 L{ellipsoidalKarney.perimeterOf}.
1828 '''
1829 def _degs(points, closed, wrap): # angular edge lengths in degrees
1830 Ps = LatLon2PsxyIter(points, wrap=wrap, loop=1)
1831 p1, u = Ps[0], _0_0 # previous x2's unroll/wrap
1832 for i, p2 in Ps.enumerate(closed=closed):
1833 w = False if closed and i == 0 else wrap
1834 # apply previous x2's unroll/wrap to new x1
1835 _, dy, dx, u = equirectangular_(p1.y, p1.x + u, p2.y, p2.x,
1836 adjust=adjust,
1837 limit=None,
1838 wrap=w) # PYCHOK non-sequence
1839 yield hypot(dx, dy)
1840 p1 = p2
1842 if _MODS.booleans.isBoolean(points):
1843 if not closed:
1844 raise _ValueError(closed=closed, points=_composite_)
1845 d = points._sum1(perimeterOf, closed=True, adjust=adjust,
1846 radius=radius, wrap=wrap)
1847 else:
1848 d = fsum(_degs(points, closed, wrap))
1849 return degrees2m(d, radius=radius)
1852def quadOf(latS, lonW, latN, lonE, closed=False, LatLon=LatLon_, **LatLon_kwds):
1853 '''Generate a quadrilateral path or polygon from two points.
1855 @arg latS: Southernmost latitude (C{degrees90}).
1856 @arg lonW: Westernmost longitude (C{degrees180}).
1857 @arg latN: Northernmost latitude (C{degrees90}).
1858 @arg lonE: Easternmost longitude (C{degrees180}).
1859 @kwarg closed: Optionally, close the path (C{bool}).
1860 @kwarg LatLon: Class to use (L{LatLon_}).
1861 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1862 keyword arguments.
1864 @return: Return a tuple of 4 or 5 B{C{LatLon}} instances
1865 outlining the quadrilateral.
1867 @see: Function L{boundsOf}.
1868 '''
1869 t = (LatLon(latS, lonW, **LatLon_kwds),
1870 LatLon(latN, lonW, **LatLon_kwds),
1871 LatLon(latN, lonE, **LatLon_kwds),
1872 LatLon(latS, lonE, **LatLon_kwds))
1873 if closed:
1874 t += t[:1]
1875 return t
1878__all__ += _ALL_DOCS(_Array2LatLon, _Basequence)
1880# **) MIT License
1881#
1882# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1883#
1884# Permission is hereby granted, free of charge, to any person obtaining a
1885# copy of this software and associated documentation files (the "Software"),
1886# to deal in the Software without restriction, including without limitation
1887# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1888# and/or sell copies of the Software, and to permit persons to whom the
1889# Software is furnished to do so, subject to the following conditions:
1890#
1891# The above copyright notice and this permission notice shall be included
1892# in all copies or substantial portions of the Software.
1893#
1894# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1895# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1896# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1897# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1898# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1899# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1900# OTHER DEALINGS IN THE SOFTWARE.