Coverage for pygeodesy/points.py: 92%
610 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-23 12:10 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-23 12:10 -0400
2# -*- coding: utf-8 -*-
4u'''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, latDMS, lonDMS, parseDMS, S_DEG, S_DMS, \
36 S_MIN, S_SEC, S_SEP
37from pygeodesy.errors import CrossError, crosserrors, _IndexError, \
38 _IsnotError, _TypeError, _ValueError, \
39 _xattr, _xkwds, _xkwds_pop
40from pygeodesy.fmath import favg, fdot, hypot, Fsum, fsum
41# from pygeodesy.fsums import Fsum, fsum # from .fmath
42from pygeodesy.formy import _bearingTo2, equirectangular_, _isequalTo, \
43 isnormal, latlon2n_xyz, normal, _spherical_datum
44from pygeodesy.interns import NN, _colinear_, _COMMASPACE_, _composite_, \
45 _DEQUALSPACED_, _ELLIPSIS_, _EW_, _immutable_, \
46 _lat_, _lon_, _LatLon_, _near_, _no_, _not_, \
47 _NS_, _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, Property_RO, \
58 property_doc_, property_RO, _update_all
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, _unrollon, unrollPI, _Wrap, wrap180
64from math import cos, fabs, fmod, radians, sin
66__all__ = _ALL_LAZY.points
67__version__ = '23.08.23'
69_ilat_ = 'ilat'
70_ilon_ = 'ilon'
71_ncols_ = 'ncols'
72_nrows_ = 'nrows'
75class LatLon_(object): # XXX in heights._HeightBase.height
76 '''Low-overhead C{LatLon} class for L{Numpy2LatLon} and L{Tuple2LatLon}.
77 '''
78 # __slots__ efficiency is voided if the __slots__ class attribute
79 # is used in a subclass of a class with the traditional __dict__,
80 # see <https://docs.Python.org/2/reference/datamodel.html#slots>
81 # and __slots__ must be repeated in sub-classes, see Luciano
82 # Ramalho, "Fluent Python", O'Reilly, 2016 p. 276+ "Problems
83 # with __slots__" (also at <https://Books.Google.ie/books?
84 # id=bIZHCgAAQBAJ&lpg=PP1&dq=fluent%20python&pg=PT364#
85 # v=onepage&q=“Problems%20with%20__slots__”&f=false>),
86 # 2022 p. 390 "Summarizing the Issues with __slots__".
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-56 bytes without in Python 2.7.18+ and Python 3+.
94 def __init__(self, latlonh, lon=None, name=NN, height=0,
95 datum=None, wrap=False):
96 '''Creat a new, mininal, low-overhead L{LatLon_} instance.
98 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or
99 a previous C{LatLon} instance provided C{B{lon}=None}.
100 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or
101 C(None), indicating B{C{latlonh}} is a C{LatLon}.
102 @kwarg name: Optional name (C{str}).
103 @kwarg height: Optional height (C{meter}, conventionally).
104 @kwarg datum: Optional datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},
105 L{a_f2Tuple} or I{scalar} radius) or C{None}.
106 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}
107 (C{bool}).
109 @raise TypeError: Invalid B{C{datum}} or B{C{latlonh}} not a C{LatLon}.
111 @note: The lat- and longitude are taken as-given,
112 un-clipped and un-validated .
113 '''
114 if lon is None: # PYCHOK no cover
115 try:
116 ll = latlonh.lat, latlonh.lon
117 height = _xattr(latlonh, height=height)
118 except AttributeError:
119 raise _IsnotError(_LatLon_, latlonh=latlonh)
120 if wrap:
121 ll = _Wrap.latlon(*ll)
122 elif wrap: # PYCHOK no cover
123 ll = _Wrap.latlonDMS2(latlonh, lon)
124 else: # must be latNS, lonEW
125 ll = parseDMS(latlonh, suffix=_NS_), parseDMS(lon, suffix=_EW_)
127 self.lat, self.lon = ll
128 self.name = str(name) if name else NN
129 self.height = height
130 self.datum = datum if datum is None else \
131 _spherical_datum(datum, name=self.name)
133 def __eq__(self, other):
134 return isinstance(other, LatLon_) and \
135 other.lat == self.lat and \
136 other.lon == self.lon
138 def __ne__(self, other):
139 return not self.__eq__(other)
141 def __repr__(self):
142 return self.toRepr()
144 def __str__(self):
145 return self.toStr()
147 def classof(self, *args, **kwds):
148 '''Instantiate this very class.
150 @arg args: Optional, positional arguments.
151 @kwarg kwds: Optional, keyword arguments.
153 @return: New instance (C{self.__class__}).
154 '''
155 return _xnamed(self.__class__(*args, **kwds), self.name)
157 def copy(self, deep=False):
158 '''Make a shallow or deep copy of this instance.
160 @kwarg deep: If C{True} make a deep, otherwise a
161 shallow copy (C{bool}).
163 @return: The copy (C{This} (sub-)class).
164 '''
165 return _xcopy(self, deep=deep)
167 def dup(self, **items):
168 '''Duplicate this instance, replacing some items.
170 @kwarg items: Attributes to be changed (C{any}).
172 @return: The duplicate (C{This} (sub-)class).
174 @raise AttributeError: Some B{C{items}} invalid.
175 '''
176 d = _xdup(self, **items)
177 if items:
178 _update_all(d)
179 return d
181 def heightStr(self, prec=-2):
182 '''Return a string for the height B{C{height}}.
184 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
186 @see: Function L{pygeodesy.hstr}.
187 '''
188 return hstr(self.height, prec=prec)
190 def intermediateTo(self, other, fraction, height=None, wrap=False):
191 '''Locate the point at a given fraction between (or along) this
192 and an other point.
194 @arg other: The other point (C{LatLon}).
195 @arg fraction: Fraction between both points (C{float},
196 0.0 for this and 1.0 for the other point).
197 @kwarg height: Optional height (C{meter}), overriding the
198 intermediate height.
199 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
200 the B{C{other}} point (C{bool}).
202 @return: Intermediate point (this C{LatLon}).
204 @raise TypeError: Incompatible B{C{other}} C{type}.
205 '''
206 f = Scalar(fraction=fraction)
207 if isnear0(f):
208 r = self
209 elif isnear1(f) and not wrap:
210 r = self.others(other)
211 else:
212 p = self.others(other)
213 h = favg(self.height, p.height, f=f) if height is None else height
214 _, lat, lon = _Wrap.latlon3(self.lon, p.lat, p.lon, wrap=wrap)
215 r = self.classof(favg(self.lat, lat, f=f),
216 favg(self.lon, lon, f=f),
217 height=h, datum=self.datum,
218 name=self.intermediateTo.__name__)
219 return r
221 @Property_RO # PYCHOK no cover
222 def isEllipsoidal(self):
223 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).
224 '''
225 return self.datum.isEllipsoidal if self.datum else None
227 @Property_RO # PYCHOK no cover
228 def isEllipsoidalLatLon(self):
229 '''Get C{LatLon} base.
230 '''
231 return False
233 def isequalTo(self, other, eps=None):
234 '''Compare this point with an other point, I{ignoring} height.
236 @arg other: The other point (C{LatLon}).
237 @kwarg eps: Tolerance for equality (C{degrees}).
239 @return: C{True} if both points are identical,
240 I{ignoring} height, C{False} otherwise.
242 @raise UnitError: Invalid B{C{eps}}.
243 '''
244 return _isequalTo(self, self.others(other), eps=eps)
246 @Property_RO
247 def isnormal(self):
248 '''Return C{True} if this point is normal (C{bool}),
249 meaning C{abs(lat) <= 90} and C{abs(lon) <= 180}.
251 @see: Methods L{normal}, L{toNormal} and functions
252 L{pygeodesy.isnormal} and L{pygeodesy.normal}.
253 '''
254 return isnormal(self.lat, self.lon, eps=0)
256 @Property_RO
257 def isSpherical(self): # PYCHOK no cover
258 '''Check whether this point is spherical (C{bool} or C{None} if unknown).
259 '''
260 return self.datum.isSpherical if self.datum else None
262 @Property_RO
263 def lam(self):
264 '''Get the longitude (B{C{radians}}).
265 '''
266 return radians(self.lon)
268 @Property
269 def latlon(self):
270 '''Get the lat- and longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
271 '''
272 return LatLon2Tuple(self.lat, self.lon, name=self.name)
274 @latlon.setter # PYCHOK setter!
275 def latlon(self, latlon):
276 '''Set the lat- and longitude.
278 @arg latlon: New lat- and longitude in C{degrees} (C{2-tuple} or {-list}).
279 '''
280 if self.latlon != latlon[:2]:
281 _update_all(self)
282 self.lat, self.lon = latlon[:2]
284 @Property_RO
285 def latlonheight(self):
286 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).
287 '''
288 return self.latlon.to3Tuple(self.height)
290 @Property_RO
291 def _N_vector(self):
292 '''(INTERNAL) Get the minimal, low-overhead (C{nvectorBase._N_vector_})
293 '''
294 return _N_vector_(*latlon2n_xyz(self.lat, self.lon),
295 h=self.height, name=self.name)
297 def others(self, *other, **name_other_up): # see .named._namedBase.others
298 '''Refined class comparison.
300 @arg other: The other instance (any C{type}).
301 @kwarg name_other_up: Overriding C{name=other} and C{up=1}
302 keyword arguments.
304 @return: The B{C{other}} if compatible.
306 @raise TypeError: Incompatible B{C{other}} C{type}.
307 '''
308 other, name, up = _xother3(self, other, **name_other_up)
309 if isinstance(other, self.__class__) or (hasattr(other, _lat_)
310 and hasattr(other, _lon_)):
311 return other
312 raise _xotherError(self, other, name=name, up=up + 1)
314 def normal(self):
315 '''Normalize this point I{in-place} to C{abs(lat) <= 90} and
316 C{abs(lon) <= 180}.
318 @return: C{True} if this point was I{normal}, C{False} if it
319 wasn't (but is now).
321 @see: Property L{isnormal} and method L{toNormal}.
322 '''
323 n = self.isnormal
324 if not n:
325 self.latlon = normal(*self.latlon)
326 return n
328 @Property_RO
329 def phi(self):
330 '''Get the latitude (B{C{radians}}).
331 '''
332 return radians(self.lat)
334 @Property_RO
335 def philam(self):
336 '''Get the lat- and longitude (L{PhiLam2Tuple}C{(phi, lam)}).
337 '''
338 return PhiLam2Tuple(self.phi, self.lam, name=self.name)
340 @Property_RO
341 def philamheight(self):
342 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
343 '''
344 return self.philam.to3Tuple(self.height)
346 @deprecated_method
347 def points(self, points, closed=False, base=None): # PYCHOK no cover
348 '''DEPRECATED, use method C{points2}.'''
349 return points2(points, closed=closed, base=base)
351 def points2(self, points, closed=False, base=None):
352 '''Check a path or polygon represented by points.
354 @arg points: The path or polygon points (C{LatLon}[])
355 @kwarg closed: Optionally, consider the polygon closed,
356 ignoring any duplicate or closing final
357 B{C{points}} (C{bool}).
358 @kwarg base: Optionally, check all B{C{points}} against
359 this base class, if C{None} don't check.
361 @return: A L{Points2Tuple}C{(number, points)} with the number
362 of points and the points C{list} or C{tuple}.
364 @raise PointsError: Insufficient number of B{C{points}}.
366 @raise TypeError: Some B{C{points}} are not B{C{base}}.
367 '''
368 return points2(points, closed=closed, base=base)
370 def PointsIter(self, points, loop=0, dedup=False):
371 '''Return a points iterator.
373 @arg points: The path or polygon points (C{LatLon}[])
374 @kwarg loop: Number of loop-back points (non-negative C{int}).
375 @kwarg dedup: Skip duplicate points (C{bool}).
377 @return: A new C{PointsIter} iterator.
379 @raise PointsError: Insufficient number of B{C{points}}.
380 '''
381 return PointsIter(points, base=self, loop=loop, dedup=dedup)
383 @deprecated_method
384 def to2ab(self): # PYCHOK no cover
385 '''DEPRECATED, use property L{philam}.'''
386 return self.philam
388 def toNormal(self, deep=False, name=NN):
389 '''Get a copy of this point normalized to C{abs(lat) <= 90}
390 and C{abs(lon) <= 180}.
392 @kwarg deep: If C{True} make a deep, otherwise a
393 shallow copy (C{bool}).
394 @kwarg name: Optional name of the copy (C{str}).
396 @return: A copy of this point, I{normalized} and
397 optionally renamed (C{LatLon}).
399 @see: Property L{isnormal}, method L{normal} and
400 function L{pygeodesy.normal}.
401 '''
402 ll = self.copy(deep=deep)
403 _ = ll.normal()
404 if name:
405 ll.name = str(name)
406 return ll
408 def toNvector(self, h=None, Nvector=NvectorBase, **Nvector_kwds):
409 '''Convert this point to C{n-vector} (normal to the earth's
410 surface) components, I{including height}.
412 @kwarg h: Optional height, overriding this point's height
413 (C{meter}).
414 @kwarg Nvector: Optional class to return the C{n-vector}
415 components (C{Nvector}) or C{None}.
416 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} keyword
417 arguments, ignored if C{B{Nvector} is None}.
419 @return: The C{n-vector} components B{C{Nvector}} or if
420 B{C{Nvector}} is C{None}, a L{Vector4Tuple}C{(x,
421 y, z, h)}.
423 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}
424 argument.
425 '''
426 x, y, z = latlon2n_xyz(self.lat, self.lon)
427 r = Vector4Tuple(x, y, z, self.height if h is None else h)
428 if Nvector is not None:
429 r = Nvector(x, y, z, **_xkwds(Nvector_kwds, h=r.h, ll=self))
430 return _xnamed(r, self.name)
432 def toRepr(self, **kwds):
433 '''This L{LatLon_} as a string "class(<degrees>, ...)".
435 @kwarg kwds: Optional, keyword arguments.
437 @return: Class instance (C{str}).
438 '''
439 _ = _xkwds_pop(kwds, std=None) # PYCHOK std unused
440 return Fmt.PAREN(classname(self), self.toStr(**kwds))
442 def toStr(self, form=F_D, joined=_COMMASPACE_, **prec_sep_s_D_M_S_kwds):
443 '''Convert this point to a "lat, lon[, height][, name][, ...]" string,
444 formatted in the given C{B{form}at}.
446 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see
447 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.
448 @kwarg joined: Separator to join the lat-, longitude, heigth, name
449 and other strings (C{str} or C{None} or C{NN} for
450 non-joined).
451 @kwarg prec_sep_s_D_M_S_kwds: Optional C{B{prec}ision}, C{B{sep}arator},
452 B{C{s_D}}, B{C{s_M}}, B{C{s_S}}, B{C{s_DMS}} and possibly
453 other keyword arguments, see functions L{pygeodesy.latDMS}
454 or L{pygeodesy.lonDMS}.
456 @return: This point in the specified C{B{form}at}, etc. (C{str} or
457 a 2- or 3+tuple C{(lat_str, lon_str[, height_str][, name_str][,
458 ...])} if C{B{joined}=NN} or C{B{joined}=None} and with the
459 C{height_str} and C{name_str} only included if non-zero
460 respectively non-empty).
462 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more
463 details about keyword arguments C{B{form}at}, C{B{prec}ision},
464 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.
465 '''
466 def _t3(prec=None, sep=S_SEP, s_D=S_DEG, s_M=S_MIN, s_S=S_SEC, s_DMS=S_DMS, **kwds):
467 return dict(prec=prec, sep=sep, s_D=s_D, s_M=s_M, s_S=s_S, s_DMS=s_DMS), kwds, prec
469 prec_sep_s_D_M_S, kwds, prec = _t3(**prec_sep_s_D_M_S_kwds)
470 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),
471 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))
472 if self.height:
473 t += (self.heightStr(),)
474 if self.name:
475 t += (repr(self.name),)
476 if kwds:
477 t += pairs(kwds) if prec is None else pairs(kwds, prec=prec)
478 return joined.join(t) if joined else t
480 @deprecated_method
481 def toStr2(self, **kwds): # PYCHOK no cover
482 '''DEPRECATED, used method L{toRepr}.'''
483 return self.toRepr(**kwds)
486def _isLatLon(inst):
487 '''(INTERNAL) Check a C{LatLon} or C{LatLon_} instance.
488 '''
489 return isinstance(inst, (LatLon_, _MODS.latlonBase.LatLonBase))
492def _isLatLon_(LL):
493 '''(INTERNAL) Check a (sub-)class of C{LatLon_}.
494 '''
495 return issubclassof(LL, LatLon_) or (isclass(LL) and
496 all(hasattr(LL, a) for a in _ALL_ATTRS_))
499# get all pseudo-slots for class C{LatLon_}
500_ALL_ATTRS_ = tuple(LatLon_(0, 0).__dict__.keys())
503class _Basequence(_Sequence): # immutable, on purpose
504 '''(INTERNAL) Base class.
505 '''
506 _array = []
507 _epsilon = EPS
508 _itemname = _point_
510 def _contains(self, point):
511 '''(INTERNAL) Check for a matching point.
512 '''
513 return any(self._findall(point, ()))
515 def copy(self, deep=False): # PYCHOK no cover
516 '''Make a shallow or deep copy of this instance.
518 @kwarg deep: If C{True} make a deep, otherwise a
519 shallow copy (C{bool}).
521 @return: The copy (C{This class} or subclass thereof).
522 '''
523 return _xcopy(self, deep=deep)
525 def _count(self, point):
526 '''(INTERNAL) Count the number of matching points.
527 '''
528 return sum(1 for _ in self._findall(point, ())) # NOT len()!
530 def dup(self, **items): # PYCHOK no cover
531 '''Duplicate this instance, I{without replacing items}.
533 @kwarg items: No attributes (I{not allowed}).
535 @return: The duplicate (C{This} (sub-)class).
537 @raise TypeError: Any B{C{items}} invalid.
538 '''
539 if items:
540 t = _SPACE_(classname(self), _immutable_)
541 raise _TypeError(txt=t, this=self, **items)
542 return _xdup(self)
544 @property_doc_(''' the equality tolerance (C{float}).''')
545 def epsilon(self):
546 '''Get the tolerance for equality tests (C{float}).
547 '''
548 return self._epsilon
550 @epsilon.setter # PYCHOK setter!
551 def epsilon(self, tol):
552 '''Set the tolerance for equality tests (C{scalar}).
554 @raise UnitError: Non-scalar or invalid B{C{tol}}.
555 '''
556 self._epsilon = Scalar_(tolerance=tol)
558 def _find(self, point, start_end):
559 '''(INTERNAL) Find the first matching point index.
560 '''
561 for i in self._findall(point, start_end):
562 return i
563 return -1
565 def _findall(self, point, start_end): # PYCHOK no cover
566 '''(INTERNAL) I{Must be implemented/overloaded}.
567 '''
568 notImplemented(self, point, start_end)
570 def _getitem(self, index):
571 '''(INTERNAL) Return point [index] or return a slice.
572 '''
573 # Luciano Ramalho, "Fluent Python", O'Reilly, 2016 p. 290+, 2022 p. 405+
574 if isinstance(index, slice):
575 # XXX an numpy.[nd]array slice is a view, not a copy
576 return self.__class__(self._array[index], **self._slicekwds())
577 else:
578 return self.point(self._array[index])
580 def _index(self, point, start_end):
581 '''(INTERNAL) Find the first matching point index.
582 '''
583 for i in self._findall(point, start_end):
584 return i
585 raise _IndexError(self._itemname, point, txt=_not_('found'))
587 @property_RO
588 def isNumpy2(self): # PYCHOK no cover
589 '''Is this a Numpy2 wrapper?
590 '''
591 return False # isinstance(self, (Numpy2LatLon, ...))
593 @property_RO
594 def isPoints2(self): # PYCHOK no cover
595 '''Is this a LatLon2 wrapper/converter?
596 '''
597 return False # isinstance(self, (LatLon2psxy, ...))
599 @property_RO
600 def isTuple2(self): # PYCHOK no cover
601 '''Is this a Tuple2 wrapper?
602 '''
603 return False # isinstance(self, (Tuple2LatLon, ...))
605 def _iter(self):
606 '''(INTERNAL) Yield all points.
607 '''
608 _array, _point = self._array, self.point
609 for i in range(len(self)):
610 yield _point(_array[i])
612 def point(self, *attrs): # PYCHOK no cover
613 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded in}.
615 @arg attrs: Optional arguments.
616 '''
617 notOverloaded(self, *attrs)
619 def _range(self, start=None, end=None, step=1):
620 '''(INTERNAL) Return the range.
621 '''
622 if step > 0:
623 if start is None:
624 start = 0
625 if end is None:
626 end = len(self)
627 elif step < 0:
628 if start is None:
629 start = len(self) - 1
630 if end is None:
631 end = -1
632 else:
633 raise _ValueError(step=step)
634 return range(start, end, step)
636 def _repr(self):
637 '''(INTERNAL) Return a string representation.
638 '''
639 # XXX use Python 3+ reprlib.repr
640 t = repr(self._array[:1]) # only first row
641 t = _SPACE_(t[:-1], _ELLIPSIS_, Fmt.SQUARE(t[-1:], len(self)))
642 t = _SPACE_.join(t.split()) # coalesce spaces
643 return instr(self, t, **self._slicekwds())
645 def _reversed(self): # PYCHOK false
646 '''(INTERNAL) Yield all points in reverse order.
647 '''
648 _array, point = self._array, self.point
649 for i in range(len(self) - 1, -1, -1):
650 yield point(_array[i])
652 def _rfind(self, point, start_end):
653 '''(INTERNAL) Find the last matching point index.
654 '''
655 def _r3(start=None, end=None, step=-1):
656 return (start, end, step) # PYCHOK returns
658 for i in self._findall(point, _r3(*start_end)):
659 return i
660 return -1
662 def _slicekwds(self): # PYCHOK no cover
663 '''(INTERNAL) I{Should be overloaded}.
664 '''
665 return {}
668class _Array2LatLon(_Basequence): # immutable, on purpose
669 '''(INTERNAL) Base class for Numpy2LatLon or Tuple2LatLon.
670 '''
671 _array = ()
672 _ilat = 0 # row column index
673 _ilon = 0 # row column index
674 _LatLon = LatLon_ # default
675 _shape = ()
677 def __init__(self, array, ilat=0, ilon=1, LatLon=None, shape=()):
678 '''Handle a C{NumPy} or C{Tuple} array as a sequence of C{LatLon} points.
679 '''
680 ais = (_ilat_, ilat), (_ilon_, ilon)
682 if len(shape) != 2 or shape[0] < 1 or shape[1] < len(ais):
683 raise _IndexError('array.shape', shape)
685 self._array = array
686 self._shape = Shape2Tuple(shape) # *shape
688 if LatLon: # check the point class
689 if not _isLatLon_(LatLon):
690 raise _IsnotError(_valid_, LatLon=LatLon)
691 self._LatLon = LatLon
693 # check the attr indices
694 for n, (ai, i) in enumerate(ais):
695 if not isint(i):
696 raise _IsnotError(int.__name__, **{ai: i})
697 i = int(i)
698 if not 0 <= i < shape[1]:
699 raise _ValueError(ai, i)
700 for aj, j in ais[:n]:
701 if int(j) == i:
702 raise _ValueError(_DEQUALSPACED_(ai, aj, i))
703 setattr(self, NN(_UNDER_, ai), i)
705 def __contains__(self, latlon):
706 '''Check for a specific lat-/longitude.
708 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
709 C{(lat, lon)}).
711 @return: C{True} if B{C{latlon}} is present, C{False} otherwise.
713 @raise TypeError: Invalid B{C{latlon}}.
714 '''
715 return self._contains(latlon)
717 def __getitem__(self, index):
718 '''Return row[index] as C{LatLon} or return a L{Numpy2LatLon} slice.
719 '''
720 return self._getitem(index)
722 def __iter__(self):
723 '''Yield rows as C{LatLon}.
724 '''
725 return self._iter()
727 def __len__(self):
728 '''Return the number of rows.
729 '''
730 return self._shape[0]
732 def __repr__(self):
733 '''Return a string representation.
734 '''
735 return self._repr()
737 def __reversed__(self): # PYCHOK false
738 '''Yield rows as C{LatLon} in reverse order.
739 '''
740 return self._reversed()
742 __str__ = __repr__
744 def count(self, latlon):
745 '''Count the number of rows with a specific lat-/longitude.
747 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
748 C{(lat, lon)}).
750 @return: Count (C{int}).
752 @raise TypeError: Invalid B{C{latlon}}.
753 '''
754 return self._count(latlon)
756 def find(self, latlon, *start_end):
757 '''Find the first row with a specific lat-/longitude.
759 @arg latlon: Point (C{LatLon}) or 2-tuple (lat, lon).
760 @arg start_end: Optional C{[start[, end]]} index (integers).
762 @return: Index or -1 if not found (C{int}).
764 @raise TypeError: Invalid B{C{latlon}}.
765 '''
766 return self._find(latlon, start_end)
768 def _findall(self, latlon, start_end):
769 '''(INTERNAL) Yield indices of all matching rows.
770 '''
771 try:
772 lat, lon = latlon.lat, latlon.lon
773 except AttributeError:
774 try:
775 lat, lon = latlon
776 except (TypeError, ValueError):
777 raise _IsnotError(_valid_, latlon=latlon)
779 _ilat, _ilon = self._ilat, self._ilon
780 _array, _eps = self._array, self._epsilon
781 for i in self._range(*start_end):
782 row = _array[i]
783 if fabs(row[_ilat] - lat) <= _eps and \
784 fabs(row[_ilon] - lon) <= _eps:
785 yield i
787 def findall(self, latlon, *start_end):
788 '''Yield indices of all rows with a specific lat-/longitude.
790 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
791 C{(lat, lon)}).
792 @arg start_end: Optional C{[start[, end]]} index (C{int}).
794 @return: Indices (C{iterable}).
796 @raise TypeError: Invalid B{C{latlon}}.
797 '''
798 return self._findall(latlon, start_end)
800 def index(self, latlon, *start_end): # PYCHOK Python 2- issue
801 '''Find index of the first row with a specific lat-/longitude.
803 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
804 C{(lat, lon)}).
805 @arg start_end: Optional C{[start[, end]]} index (C{int}).
807 @return: Index (C{int}).
809 @raise IndexError: Point not found.
811 @raise TypeError: Invalid B{C{latlon}}.
812 '''
813 return self._index(latlon, start_end)
815 @Property_RO
816 def ilat(self):
817 '''Get the latitudes column index (C{int}).
818 '''
819 return self._ilat
821 @Property_RO
822 def ilon(self):
823 '''Get the longitudes column index (C{int}).
824 '''
825 return self._ilon
827# next = __iter__
829 def point(self, row): # PYCHOK *attrs
830 '''Instantiate a point C{LatLon}.
832 @arg row: Array row (numpy.array).
834 @return: Point (C{LatLon}).
835 '''
836 return self._LatLon(row[self._ilat], row[self._ilon])
838 def rfind(self, latlon, *start_end):
839 '''Find the last row with a specific lat-/longitude.
841 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
842 C{(lat, lon)}).
843 @arg start_end: Optional C{[start[, end]]} index (C{int}).
845 @note: Keyword order, first stop, then start.
847 @return: Index or -1 if not found (C{int}).
849 @raise TypeError: Invalid B{C{latlon}}.
850 '''
851 return self._rfind(latlon, start_end)
853 def _slicekwds(self):
854 '''(INTERNAL) Slice kwds.
855 '''
856 return dict(ilat=self._ilat, ilon=self._ilon)
858 @Property_RO
859 def shape(self):
860 '''Get the shape of the C{NumPy} array or the C{Tuples} as
861 L{Shape2Tuple}C{(nrows, ncols)}.
862 '''
863 return self._shape
865 def _subset(self, indices): # PYCHOK no cover
866 '''(INTERNAL) I{Must be implemented/overloaded}.
867 '''
868 notImplemented(self, indices)
870 def subset(self, indices):
871 '''Return a subset of the C{NumPy} array.
873 @arg indices: Row indices (C{range} or C{int}[]).
875 @note: A C{subset} is different from a C{slice} in 2 ways:
876 (a) the C{subset} is typically specified as a list of
877 (un-)ordered indices and (b) the C{subset} allocates
878 a new, separate C{NumPy} array while a C{slice} is
879 just an other C{view} of the original C{NumPy} array.
881 @return: Sub-array (C{numpy.array}).
883 @raise IndexError: Out-of-range B{C{indices}} value.
885 @raise TypeError: If B{C{indices}} is not a C{range}
886 nor an C{int}[].
887 '''
888 if not issequence(indices, tuple): # NO tuple, only list
889 # and range work properly to get Numpy array sub-sets
890 raise _IsnotError(_valid_, indices=type(indices))
892 n = len(self)
893 for i, v in enumerate(indices):
894 if not isint(v):
895 raise _TypeError(Fmt.SQUARE(indices=i), v)
896 elif not 0 <= v < n:
897 raise _IndexError(Fmt.SQUARE(indices=i), v)
899 return self._subset(indices)
902class LatLon2psxy(_Basequence):
903 '''Wrapper for C{LatLon} points as "on-the-fly" pseudo-xy coordinates.
904 '''
905 _closed = False
906 _len = 0
907 _deg2m = None # default, keep degrees
908 _radius = None
909 _wrap = True
911 def __init__(self, latlons, closed=False, radius=None, wrap=True):
912 '''Handle C{LatLon} points as pseudo-xy coordinates.
914 @note: The C{LatLon} latitude is considered the I{pseudo-y}
915 and longitude the I{pseudo-x} coordinate, likewise
916 for L{LatLon2Tuple}. However, 2-tuples C{(x, y)} are
917 considered as I{(longitude, latitude)}.
919 @arg latlons: Points C{list}, C{sequence}, C{set}, C{tuple},
920 etc. (C{LatLon[]}).
921 @kwarg closed: Optionally, close the polygon (C{bool}).
922 @kwarg radius: Mean earth radius (C{meter}).
923 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
924 the B{C{latlons}} points (C{bool}).
926 @raise PointsError: Insufficient number of B{C{latlons}}.
928 @raise TypeError: Some B{C{points}} are not B{C{base}}.
929 '''
930 self._closed = closed
931 self._len, self._array = points2(latlons, closed=closed)
932 if radius:
933 self._radius = r = Radius(radius)
934 self._deg2m = degrees2m(_1_0, r)
935 if not wrap:
936 self._wrap = False
938 def __contains__(self, xy):
939 '''Check for a matching point.
941 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
942 C{(x, y)}) in (C{degrees}.
944 @return: C{True} if B{C{xy}} is present, C{False} otherwise.
946 @raise TypeError: Invalid B{C{xy}}.
947 '''
948 return self._contains(xy)
950 def __getitem__(self, index):
951 '''Return the pseudo-xy or return a L{LatLon2psxy} slice.
952 '''
953 return self._getitem(index)
955 def __iter__(self):
956 '''Yield all pseudo-xy's.
957 '''
958 return self._iter()
960 def __len__(self):
961 '''Return the number of pseudo-xy's.
962 '''
963 return self._len
965 def __repr__(self):
966 '''Return a string representation.
967 '''
968 return self._repr()
970 def __reversed__(self): # PYCHOK false
971 '''Yield all pseudo-xy's in reverse order.
972 '''
973 return self._reversed()
975 __str__ = __repr__
977 def count(self, xy):
978 '''Count the number of matching points.
980 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
981 C{(x, y)}) in (C{degrees}.
983 @return: Count (C{int}).
985 @raise TypeError: Invalid B{C{xy}}.
986 '''
987 return self._count(xy)
989 def find(self, xy, *start_end):
990 '''Find the first matching point.
992 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
993 C{(x, y)}) in (C{degrees}.
994 @arg start_end: Optional C{[start[, end]]} index (C{int}).
996 @return: Index or -1 if not found (C{int}).
998 @raise TypeError: Invalid B{C{xy}}.
999 '''
1000 return self._find(xy, start_end)
1002 def _findall(self, xy, start_end):
1003 '''(INTERNAL) Yield indices of all matching points.
1004 '''
1005 try:
1006 x, y = xy.lon, xy.lat
1008 def _x_y_ll3(ll): # match LatLon
1009 return ll.lon, ll.lat, ll
1011 except AttributeError:
1012 try:
1013 x, y = xy[:2]
1014 except (IndexError, TypeError, ValueError):
1015 raise _IsnotError(_valid_, xy=xy)
1017 _x_y_ll3 = self.point # PYCHOK expected
1019 _array, _eps = self._array, self._epsilon
1020 for i in self._range(*start_end):
1021 xi, yi, _ = _x_y_ll3(_array[i])
1022 if fabs(xi - x) <= _eps and \
1023 fabs(yi - y) <= _eps:
1024 yield i
1026 def findall(self, xy, *start_end):
1027 '''Yield indices of all matching points.
1029 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
1030 C{(x, y)}) in (C{degrees}.
1031 @arg start_end: Optional C{[start[, end]]} index (C{int}).
1033 @return: Indices (C{iterator}).
1035 @raise TypeError: Invalid B{C{xy}}.
1036 '''
1037 return self._findall(xy, start_end)
1039 def index(self, xy, *start_end): # PYCHOK Python 2- issue
1040 '''Find the first matching point.
1042 @arg xy: Point (C{LatLon}) or 2-tuple (x, y) in (C{degrees}).
1043 @arg start_end: Optional C{[start[, end]]} index (C{int}).
1045 @return: Index (C{int}).
1047 @raise IndexError: Point not found.
1049 @raise TypeError: Invalid B{C{xy}}.
1050 '''
1051 return self._index(xy, start_end)
1053 @property_RO
1054 def isPoints2(self):
1055 '''Is this a LatLon2 wrapper/converter?
1056 '''
1057 return True # isinstance(self, (LatLon2psxy, ...))
1059 def point(self, ll): # PYCHOK *attrs
1060 '''Create a pseudo-xy.
1062 @arg ll: Point (C{LatLon}).
1064 @return: An L{Point3Tuple}C{(x, y, ll)}.
1065 '''
1066 x, y = ll.lon, ll.lat # note, x, y = lon, lat
1067 if self._wrap:
1068 y, x = _Wrap.latlon(y, x)
1069 d = self._deg2m
1070 if d: # convert degrees to meter (or radians)
1071 x *= d
1072 y *= d
1073 return Point3Tuple(x, y, ll)
1075 def rfind(self, xy, *start_end):
1076 '''Find the last matching point.
1078 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple
1079 C{(x, y)}) in (C{degrees}.
1080 @arg start_end: Optional C{[start[, end]]} index (C{int}).
1082 @return: Index or -1 if not found (C{int}).
1084 @raise TypeError: Invalid B{C{xy}}.
1085 '''
1086 return self._rfind(xy, start_end)
1088 def _slicekwds(self):
1089 '''(INTERNAL) Slice kwds.
1090 '''
1091 return dict(closed=self._closed, radius=self._radius, wrap=self._wrap)
1094class Numpy2LatLon(_Array2LatLon): # immutable, on purpose
1095 '''Wrapper for C{NumPy} arrays as "on-the-fly" C{LatLon} points.
1096 '''
1097 def __init__(self, array, ilat=0, ilon=1, LatLon=None):
1098 '''Handle a C{NumPy} array as a sequence of C{LatLon} points.
1100 @arg array: C{NumPy} array (C{numpy.array}).
1101 @kwarg ilat: Optional index of the latitudes column (C{int}).
1102 @kwarg ilon: Optional index of the longitudes column (C{int}).
1103 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
1105 @raise IndexError: If B{C{array.shape}} is not (1+, 2+).
1107 @raise TypeError: If B{C{array}} is not a C{NumPy} array or
1108 C{LatLon} is not a class with C{lat}
1109 and C{lon} attributes.
1111 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
1112 are the same or out of range.
1114 @example:
1116 >>> type(array)
1117 <type 'numpy.ndarray'> # <class ...> in Python 3+
1118 >>> points = Numpy2LatLon(array, lat=0, lon=1)
1119 >>> simply = simplifyRDP(points, ...)
1120 >>> type(simply)
1121 <type 'numpy.ndarray'> # <class ...> in Python 3+
1122 >>> sliced = points[1:-1]
1123 >>> type(sliced)
1124 <class '...Numpy2LatLon'>
1125 '''
1126 try: # get shape and check some other numpy.array attrs
1127 s, _, _ = array.shape, array.nbytes, array.ndim # PYCHOK expected
1128 except AttributeError:
1129 raise _IsnotError('NumPy', array=type(array))
1131 _Array2LatLon.__init__(self, array, ilat=ilat, ilon=ilon,
1132 LatLon=LatLon, shape=s)
1134 @property_RO
1135 def isNumpy2(self):
1136 '''Is this a Numpy2 wrapper?
1137 '''
1138 return True # isinstance(self, (Numpy2LatLon, ...))
1140 def _subset(self, indices):
1141 return self._array[indices] # NumPy special
1144class Shape2Tuple(_NamedTuple):
1145 '''2-Tuple C{(nrows, ncols)}, the number of rows and columns,
1146 both C{int}.
1147 '''
1148 _Names_ = (_nrows_, _ncols_)
1149 _Units_ = ( Number_, Number_)
1152class Tuple2LatLon(_Array2LatLon):
1153 '''Wrapper for tuple sequences as "on-the-fly" C{LatLon} points.
1154 '''
1155 def __init__(self, tuples, ilat=0, ilon=1, LatLon=None):
1156 '''Handle a list of tuples, each containing a lat- and longitude
1157 and perhaps other values as a sequence of C{LatLon} points.
1159 @arg tuples: The C{list}, C{tuple} or C{sequence} of tuples (C{tuple}[]).
1160 @kwarg ilat: Optional index of the latitudes value (C{int}).
1161 @kwarg ilon: Optional index of the longitudes value (C{int}).
1162 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).
1164 @raise IndexError: If C{(len(B{tuples}), min(len(t) for t
1165 in B{tuples}))} is not (1+, 2+).
1167 @raise TypeError: If B{C{tuples}} is not a C{list}, C{tuple}
1168 or C{sequence} or if B{C{LatLon}} is not a
1169 C{LatLon} with C{lat}, C{lon} and C{name}
1170 attributes.
1172 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values
1173 are the same or out of range.
1175 @example:
1177 >>> tuples = [(0, 1), (2, 3), (4, 5)]
1178 >>> type(tuples)
1179 <type 'list'> # <class ...> in Python 3+
1180 >>> points = Tuple2LatLon(tuples, lat=0, lon=1)
1181 >>> simply = simplifyRW(points, 0.5, ...)
1182 >>> type(simply)
1183 <type 'list'> # <class ...> in Python 3+
1184 >>> simply
1185 [(0, 1), (4, 5)]
1186 >>> sliced = points[1:-1]
1187 >>> type(sliced)
1188 <class '...Tuple2LatLon'>
1189 >>> sliced
1190 ...Tuple2LatLon([(2, 3), ...][1], ilat=0, ilon=1)
1192 >>> closest, _ = nearestOn2(LatLon_(2, 1), points, adjust=False)
1193 >>> closest
1194 LatLon_(lat=1.0, lon=2.0)
1196 >>> closest, _ = nearestOn2(LatLon_(3, 2), points)
1197 >>> closest
1198 LatLon_(lat=2.001162, lon=3.001162)
1199 '''
1200 _xinstanceof(list, tuple, tuples=tuples)
1201 s = len(tuples), min(len(_) for _ in tuples)
1202 _Array2LatLon.__init__(self, tuples, ilat=ilat, ilon=ilon,
1203 LatLon=LatLon, shape=s)
1205 @property_RO
1206 def isTuple2(self):
1207 '''Is this a Tuple2 wrapper?
1208 '''
1209 return True # isinstance(self, (Tuple2LatLon, ...))
1211 def _subset(self, indices):
1212 return type(self._array)(self._array[i] for i in indices)
1215def _area2(points, adjust, wrap):
1216 '''(INTERNAL) Approximate the area in radians squared, I{signed}.
1217 '''
1218 if adjust:
1219 # approximate trapezoid by a rectangle, adjusting
1220 # the top width by the cosine of the latitudinal
1221 # average and bottom width by some fudge factor
1222 def _adjust(w, h):
1223 c = cos(h) if fabs(h) < PI_2 else _0_0
1224 return w * h * (c + 1.2876) * _0_5
1225 else:
1226 def _adjust(w, h): # PYCHOK expected
1227 return w * h
1229 # setting radius=1 converts degrees to radians
1230 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)
1231 x1, y1, ll = Ps[0]
1232 pts = [ll] # for _areaError
1234 A2 = Fsum() # trapezoidal area in radians**2
1235 for p in Ps.iterate(closed=True):
1236 x2, y2, ll = p
1237 if len(pts) < 4:
1238 pts.append(ll)
1239 w, x2 = unrollPI(x1, x2, wrap=wrap and not Ps.looped)
1240 A2 += _adjust(w, (y2 + y1) * _0_5)
1241 x1, y1 = x2, y2
1243 return A2.fsum(), tuple(pts)
1246def _areaError(pts, near_=NN): # imported by .ellipsoidalKarney
1247 '''(INTERNAL) Area issue.
1248 '''
1249 t = _ELLIPSIS_(pts[:3], NN)
1250 return _ValueError(NN(near_, 'zero or polar area'), txt=t)
1253def areaOf(points, adjust=True, radius=R_M, wrap=True):
1254 '''Approximate the area of a polygon or composite.
1256 @arg points: The polygon points or clips (C{LatLon}[],
1257 L{BooleanFHP} or L{BooleanGH}).
1258 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1259 by the cosine of the mean latitude (C{bool}).
1260 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1261 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1262 the B{C{points}} (C{bool}).
1264 @return: Approximate area (I{square} C{meter}, same units as
1265 B{C{radius}} or C{radians} I{squared} if B{C{radius}}
1266 is C{None}).
1268 @raise PointsError: Insufficient number of B{C{points}}
1270 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1272 @raise ValueError: Invalid B{C{radius}}.
1274 @note: This area approximation has limited accuracy and is
1275 ill-suited for regions exceeding several hundred Km
1276 or Miles or with near-polar latitudes.
1278 @see: L{sphericalNvector.areaOf}, L{sphericalTrigonometry.areaOf},
1279 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.
1280 '''
1281 if _MODS.booleans.isBoolean(points):
1282 a = points._sum1(areaOf, adjust=adjust, radius=None, wrap=wrap)
1283 else:
1284 a, _ = _area2(points, adjust, wrap)
1285 return fabs(a if radius is None else (Radius(radius)**2 * a))
1288def boundsOf(points, wrap=False, LatLon=None): # was=True
1289 '''Determine the bottom-left SW and top-right NE corners of a
1290 path or polygon.
1292 @arg points: The path or polygon points (C{LatLon}[]).
1293 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1294 the B{C{points}} (C{bool}).
1295 @kwarg LatLon: Optional class to return the C{bounds}
1296 corners (C{LatLon}) or C{None}.
1298 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)} as
1299 B{C{LatLon}}s if B{C{LatLon}} is C{None} a
1300 L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}.
1302 @raise PointsError: Insufficient number of B{C{points}}
1304 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1306 @see: Function L{quadOf}.
1308 @example:
1310 >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1311 >>> boundsOf(b) # False
1312 >>> 45.0, 1.0, 46.0, 2.0
1313 '''
1314 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)
1315 w, s, _ = e, n, _ = Ps[0]
1317 v = w
1318 for x, y, _ in Ps.iterate(closed=False): # [1:]
1319 if wrap:
1320 _, x = unroll180(v, x, wrap=True)
1321 v = x
1323 if w > x:
1324 w = x
1325 elif e < x:
1326 e = x
1328 if s > y:
1329 s = y
1330 elif n < y:
1331 n = y
1333 return Bounds4Tuple(s, w, n, e) if LatLon is None else \
1334 Bounds2Tuple(LatLon(s, w), LatLon(n, e)) # PYCHOK inconsistent
1337def _distanceTo(Error, **name_points): # .frechet, .hausdorff, .heights
1338 '''(INTERNAL) Chech callable C{distanceTo} methods.
1339 '''
1340 name, ps = name_points.popitem()
1341 for i, p in enumerate(ps):
1342 if not callable(_xattr(p, distanceTo=None)):
1343 n = _distanceTo.__name__[1:]
1344 t = _SPACE_(_no_, callable.__name__, n)
1345 raise Error(Fmt.SQUARE(name, i), p, txt=t)
1346 return ps
1349def centroidOf(points, wrap=False, LatLon=None): # was=True
1350 '''Determine the centroid of a polygon.
1352 @arg points: The polygon points (C{LatLon}[]).
1353 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1354 B{C{points}} (C{bool}).
1355 @kwarg LatLon: Optional class to return the centroid (C{LatLon})
1356 or C{None}.
1358 @return: Centroid (B{C{LatLon}}) or a L{LatLon2Tuple}C{(lat, lon)}
1359 if C{B{LatLon} is None}.
1361 @raise PointsError: Insufficient number of B{C{points}}
1363 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1365 @raise ValueError: The B{C{points}} enclose a pole or
1366 near-zero area.
1368 @see: U{Centroid<https://WikiPedia.org/wiki/Centroid#Of_a_polygon>} and
1369 Paul Bourke's U{Calculating The Area And Centroid Of A Polygon
1370 <https://www.SEAS.UPenn.edu/~ese502/lab-content/extra_materials/
1371 Polygon%20Area%20and%20Centroid.pdf>}, 1988.
1372 '''
1373 A, X, Y = Fsum(), Fsum(), Fsum()
1375 # setting radius=1 converts degrees to radians
1376 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)
1377 x1, y1, ll = Ps[0]
1378 pts = [ll] # for _areaError
1379 for p in Ps.iterate(closed=True):
1380 x2, y2, ll = p
1381 if len(pts) < 4:
1382 pts.append(ll)
1383 if wrap and not Ps.looped:
1384 _, x2 = unrollPI(x1, x2, wrap=True)
1385 t = x1 * y2 - x2 * y1
1386 A += t
1387 X += t * (x1 + x2)
1388 Y += t * (y1 + y2)
1389 # XXX more elaborately:
1390 # t1, t2 = x1 * y2, -(x2 * y1)
1391 # A.fadd_(t1, t2)
1392 # X.fadd_(t1 * x1, t1 * x2, t2 * x1, t2 * x2)
1393 # Y.fadd_(t1 * y1, t1 * y2, t2 * y1, t2 * y2)
1394 x1, y1 = x2, y2
1396 a = A.fmul(_6_0).fover(_2_0)
1397 if isnear0(a):
1398 raise _areaError(pts, near_=_near_)
1399 y, x = degrees90(Y.fover(a)), degrees180(X.fover(a))
1400 return LatLon2Tuple(y, x) if LatLon is None else LatLon(y, x)
1403def fractional(points, fi, j=None, wrap=None, LatLon=None, Vector=None, **kwds):
1404 '''Return the point at a given I{fractional} index.
1406 @arg points: The points (C{LatLon}[], L{Numpy2LatLon}[],
1407 L{Tuple2LatLon}[], C{Cartesian}[], C{Vector3d}[],
1408 L{Vector3Tuple}[]).
1409 @arg fi: The fractional index (L{FIx}, C{float} or C{int}).
1410 @kwarg j: Optionally, index of the other point (C{int}).
1411 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1412 B{{points}} (C{bool}) or C{None} for a backward
1413 compatible L{LatLon2Tuple} or B{C{LatLon}} with
1414 averaged lat- and longitudes. Use C{True} or
1415 C{False} to get the I{fractional} point computed
1416 by method C{B{points}[fi].intermediateTo}.
1417 @kwarg LatLon: Optional class to return the I{intermediate},
1418 I{fractional} point (C{LatLon}) or C{None}.
1419 @kwarg Vector: Optional class to return the I{intermediate},
1420 I{fractional} point (C{Cartesian}, C{Vector3d})
1421 or C{None}.
1422 @kwarg kwds: Optional, additional B{C{LatLon}} I{or} B{C{Vector}}
1423 keyword arguments, ignored if both C{B{LatLon}} and
1424 C{B{Vector}} are C{None}.
1426 @return: A L{LatLon2Tuple}C{(lat, lon)} if B{C{wrap}}, B{C{LatLon}}
1427 and B{C{Vector}} all are C{None}, the defaults.
1429 An instance of B{C{LatLon}} if not C{None} I{or} an instance
1430 of B{C{Vector}} if not C{None}.
1432 Otherwise with B{C{wrap}} either C{True} or C{False} and
1433 B{C{LatLon}} and B{C{Vector}} both C{None}, an instance of
1434 B{C{points}}' (sub-)class C{intermediateTo} I{fractional}.
1436 Summarized as follows:
1438 >>> wrap | LatLon | Vector | returned type/value
1439 # -------+--------+--------+--------------+------
1440 # | | | LatLon2Tuple | favg
1441 # None | None | None | or** |
1442 # | | | Vector3Tuple | favg
1443 # None | LatLon | None | LatLon | favg
1444 # None | None | Vector | Vector | favg
1445 # -------+--------+--------+--------------+------
1446 # True | None | None | points' | .iTo
1447 # True | LatLon | None | LatLon | .iTo
1448 # True | None | Vector | Vector | .iTo
1449 # -------+--------+--------+--------------+------
1450 # False | None | None | points' | .iTo
1451 # False | LatLon | None | LatLon | .iTo
1452 # False | None | Vector | Vector | .iTo
1453 # _____
1454 # favg) averaged lat, lon or x, y, z values
1455 # .iTo) value from points[fi].intermediateTo
1456 # **) depends on base class of points[fi]
1458 @raise IndexError: Fractional index B{C{fi}} invalid or B{C{points}}
1459 not subscriptable or not closed.
1461 @raise TypeError: Invalid B{C{LatLon}}, B{C{Vector}} or B{C{kwds}}
1462 argument.
1464 @see: Class L{FIx} and method L{FIx.fractional}.
1465 '''
1466 if LatLon and Vector: # PYCHOK no cover
1467 kwds = _xkwds(kwds, fi=fi, LatLon=LatLon, Vector=Vector)
1468 raise _TypeError(txt=fractional.__name__, **kwds)
1469 w = wrap if LatLon else False # intermediateTo
1470 try:
1471 if not isscalar(fi) or fi < 0:
1472 raise IndexError
1473 n = _xattr(fi, fin=0)
1474 p = _fractional(points, fi, j, fin=n, wrap=w) # see .units.FIx
1475 if LatLon:
1476 p = LatLon(p.lat, p.lon, **kwds)
1477 elif Vector:
1478 p = Vector(p.x, p.y, p.z, **kwds)
1479 except (IndexError, TypeError):
1480 raise _IndexError(fi=fi, points=points, wrap=w, txt=fractional.__name__)
1481 return p
1484def _fractional(points, fi, j, fin=None, wrap=None): # in .frechet.py
1485 '''(INTERNAL) Compute point at L{fractional} index C{fi} and C{j}.
1486 '''
1487 i = int(fi)
1488 p = points[i]
1489 r = fi - float(i)
1490 if r > EPS: # EPS0?
1491 if j is None: # in .frechet.py
1492 j = i + 1
1493 if fin:
1494 j %= fin
1495 q = points[j]
1496 if r >= EPS1: # PYCHOK no cover
1497 p = q
1498 elif wrap is not None: # in (True, False)
1499 p = p.intermediateTo(q, r, wrap=wrap)
1500 elif _isLatLon(p): # backward compatible default
1501 p = LatLon2Tuple(favg(p.lat, q.lat, f=r),
1502 favg(p.lon, q.lon, f=r),
1503 name=fractional.__name__)
1504 else: # assume p and q are cartesian or vectorial
1505 z = p.z if p.z is q.z else favg(p.z, q.z, f=r)
1506 p = Vector3Tuple(favg(p.x, q.x, f=r),
1507 favg(p.y, q.y, f=r), z,
1508 name=fractional.__name__)
1509 return p
1512def isclockwise(points, adjust=False, wrap=True):
1513 '''Determine the direction of a path or polygon.
1515 @arg points: The path or polygon points (C{LatLon}[]).
1516 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1517 by the cosine of the mean latitude (C{bool}).
1518 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1519 B{C{points}} (C{bool}).
1521 @return: C{True} if B{C{points}} are clockwise, C{False} otherwise.
1523 @raise PointsError: Insufficient number of B{C{points}}
1525 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1527 @raise ValueError: The B{C{points}} enclose a pole or zero area.
1529 @example:
1531 >>> f = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
1532 >>> isclockwise(f) # False
1533 >>> isclockwise(reversed(f)) # True
1534 '''
1535 a, pts = _area2(points, adjust, wrap)
1536 if a > 0: # opposite of ellipsoidalExact and -Karney
1537 return True
1538 elif a < 0:
1539 return False
1540 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1541 raise _areaError(pts)
1544def isconvex(points, adjust=False, wrap=False): # was=True
1545 '''Determine whether a polygon is convex.
1547 @arg points: The polygon points (C{LatLon}[]).
1548 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1549 by the cosine of the mean latitude (C{bool}).
1550 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1551 B{C{points}} (C{bool}).
1553 @return: C{True} if B{C{points}} are convex, C{False} otherwise.
1555 @raise CrossError: Some B{C{points}} are colinear.
1557 @raise PointsError: Insufficient number of B{C{points}}
1559 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1561 @example:
1563 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1564 >>> isconvex(t) # True
1566 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1567 >>> isconvex(f) # False
1568 '''
1569 return bool(isconvex_(points, adjust=adjust, wrap=wrap))
1572def isconvex_(points, adjust=False, wrap=False): # was=True
1573 '''Determine whether a polygon is convex I{and clockwise}.
1575 @arg points: The polygon points (C{LatLon}[]).
1576 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1577 by the cosine of the mean latitude (C{bool}).
1578 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1579 B{C{points}} (C{bool}).
1581 @return: C{+1} if B{C{points}} are convex clockwise, C{-1} for
1582 convex counter-clockwise B{C{points}}, C{0} otherwise.
1584 @raise CrossError: Some B{C{points}} are colinear.
1586 @raise PointsError: Insufficient number of B{C{points}}
1588 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1590 @example:
1592 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)
1593 >>> isconvex_(t) # +1
1595 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)
1596 >>> isconvex_(f) # 0
1597 '''
1598 if adjust:
1599 def _unroll2(x1, x2, w, y1, y2):
1600 x21, x2 = unroll180(x1, x2, wrap=w)
1601 y = radians(y1 + y2) * _0_5
1602 x21 *= cos(y) if fabs(y) < PI_2 else _0_0
1603 return x21, x2
1604 else:
1605 def _unroll2(x1, x2, w, *unused): # PYCHOK expected
1606 return unroll180(x1, x2, wrap=w)
1608 c, s = crosserrors(), 0
1610 Ps = LatLon2PsxyIter(points, loop=2, wrap=wrap)
1611 x1, y1, _ = Ps[0]
1612 x2, y2, _ = Ps[1]
1614 x21, x2 = _unroll2(x1, x2, False, y1, y2)
1615 for i, p in Ps.enumerate(closed=True):
1616 x3, y3, ll = p
1617 x32, x3 = _unroll2(x2, x3, bool(wrap and not Ps.looped), y2, y3)
1619 # get the sign of the distance from point
1620 # x3, y3 to the line from x1, y1 to x2, y2
1621 # <https://WikiPedia.org/wiki/Distance_from_a_point_to_a_line>
1622 s3 = fdot((x3, y3, x1, y1), y2 - y1, -x21, -y2, x2)
1623 if s3 > 0: # x3, y3 on the right
1624 if s < 0: # non-convex
1625 return 0
1626 s = +1
1628 elif s3 < 0: # x3, y3 on the left
1629 if s > 0: # non-convex
1630 return 0
1631 s = -1
1633 elif c and fdot((x32, y1 - y2), y3 - y2, -x21) < 0: # PYCHOK no cover
1634 # colinear u-turn: x3, y3 not on the
1635 # opposite side of x2, y2 as x1, y1
1636 t = Fmt.SQUARE(points=i)
1637 raise CrossError(t, ll, txt=_colinear_)
1639 x1, y1, x2, y2, x21 = x2, y2, x3, y3, x32
1641 return s # all points on the same side
1644def isenclosedBy(point, points, wrap=False): # MCCABE 15
1645 '''Determine whether a point is enclosed by a polygon or composite.
1647 @arg point: The point (C{LatLon} or 2-tuple C{(lat, lon)}).
1648 @arg points: The polygon points or clips (C{LatLon}[], L{BooleanFHP}
1649 or L{BooleanGH}).
1650 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1651 B{C{points}} (C{bool}).
1653 @return: C{True} if the B{C{point}} is inside the polygon or
1654 composite, C{False} otherwise.
1656 @raise PointsError: Insufficient number of B{C{points}}
1658 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1660 @raise ValueError: Invalid B{C{point}}, lat- or longitude.
1662 @see: Functions L{pygeodesy.isconvex} and L{pygeodesy.ispolar} especially
1663 if the B{C{points}} may enclose a pole or wrap around the earth
1664 I{longitudinally}, methods L{sphericalNvector.LatLon.isenclosedBy},
1665 L{sphericalTrigonometry.LatLon.isenclosedBy} and U{MultiDop
1666 GeogContainPt<https://GitHub.com/NASA/MultiDop>} (U{Shapiro et.al. 2009,
1667 JTECH<https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}
1668 and U{Potvin et al. 2012, JTECH <https://Journals.AMetSoc.org/doi/abs/
1669 10.1175/JTECH-D-11-00019.1>}).
1670 '''
1671 try:
1672 y0, x0 = point.lat, point.lon
1673 except AttributeError:
1674 try:
1675 y0, x0 = map(float, point[:2])
1676 except (IndexError, TypeError, ValueError) as x:
1677 raise _ValueError(point=point, cause=x)
1679 if wrap:
1680 y0, x0 = _Wrap.latlon(y0, x0)
1682 def _dxy3(x, x2, y2, Ps):
1683 dx, x2 = unroll180(x, x2, wrap=not Ps.looped)
1684 return dx, x2, y2
1686 else:
1687 x0 = fmod(x0, _360_0) # not x0 % 360!
1688 x0_180_ = x0 - _180_0
1689 x0_180 = x0 + _180_0
1691 def _dxy3(x1, x, y, unused): # PYCHOK expected
1692 x = _umod_360(float(x))
1693 if x < x0_180_:
1694 x += _360_0
1695 elif x >= x0_180:
1696 x -= _360_0
1697 return (x - x1), x, y
1699 if _MODS.booleans.isBoolean(points):
1700 return points._encloses(y0, x0, wrap=wrap)
1702 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)
1703 p = Ps[0]
1704 e = m = False
1705 S = Fsum()
1707 _, x1, y1 = _dxy3(x0, p.x, p.y, False)
1708 for p in Ps.iterate(closed=True):
1709 dx, x2, y2 = _dxy3(x1, p.x, p.y, Ps)
1710 # ignore duplicate and near-duplicate pts
1711 if fabs(dx) > EPS or fabs(y2 - y1) > EPS:
1712 # determine if polygon edge (x1, y1)..(x2, y2) straddles
1713 # point (lat, lon) or is on boundary, but do not count
1714 # edges on boundary as more than one crossing
1715 if fabs(dx) < 180 and (x1 < x0 <= x2 or x2 < x0 <= x1):
1716 m = not m
1717 dy = (x0 - x1) * (y2 - y1) - (y0 - y1) * dx
1718 if (dy > 0 and dx >= 0) or (dy < 0 and dx <= 0):
1719 e = not e
1721 S += sin(radians(y2))
1722 x1, y1 = x2, y2
1724 # An odd number of meridian crossings means, the polygon
1725 # contains a pole. Assume it is the pole on the hemisphere
1726 # containing the polygon mean point and if the polygon does
1727 # contain the North Pole, flip the result.
1728 if m and S.fsum() > 0:
1729 e = not e
1730 return e
1733def ispolar(points, wrap=False):
1734 '''Check whether a polygon encloses a pole.
1736 @arg points: The polygon points (C{LatLon}[]).
1737 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1738 the B{C{points}} (C{bool}).
1740 @return: C{True} if the polygon encloses a pole, C{False}
1741 otherwise.
1743 @raise PointsError: Insufficient number of B{C{points}}
1745 @raise TypeError: Some B{C{points}} are not C{LatLon} or don't
1746 have C{bearingTo2}, C{initialBearingTo}
1747 and C{finalBearingTo} methods.
1748 '''
1749 def _cds(ps, w): # iterate over course deltas
1750 Ps = PointsIter(ps, loop=2, wrap=w)
1751 p2, p1 = Ps[0:2]
1752 b1, _ = _bearingTo2(p2, p1, wrap=False)
1753 for p2 in Ps.iterate(closed=True):
1754 if not p2.isequalTo(p1, EPS):
1755 if w and not Ps.looped:
1756 p2 = _unrollon(p1, p2)
1757 b, b2 = _bearingTo2(p1, p2, wrap=False)
1758 yield wrap180(b - b1) # (b - b1 + 540) % 360 - 180
1759 yield wrap180(b2 - b) # (b2 - b + 540) % 360 - 180
1760 p1, b1 = p2, b2
1762 # summation of course deltas around pole is 0° rather than normally ±360°
1763 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
1764 s = fsum(_cds(points, wrap), floats=True)
1765 # XXX fix (intermittant) edge crossing pole - eg (85,90), (85,0), (85,-90)
1766 return fabs(s) < 90 # "zero-ish"
1769def luneOf(lon1, lon2, closed=False, LatLon=LatLon_, **LatLon_kwds):
1770 '''Generate an ellipsoidal or spherical U{lune
1771 <https://WikiPedia.org/wiki/Spherical_lune>}-shaped path or polygon.
1773 @arg lon1: Left longitude (C{degrees90}).
1774 @arg lon2: Right longitude (C{degrees90}).
1775 @kwarg closed: Optionally, close the path (C{bool}).
1776 @kwarg LatLon: Class to use (L{LatLon_}).
1777 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1778 keyword arguments.
1780 @return: A tuple of 4 or 5 B{C{LatLon}} instances outlining
1781 the lune shape.
1783 @see: U{Latitude-longitude quadrangle
1784 <https://www.MathWorks.com/help/map/ref/areaquad.html>}.
1785 '''
1786 t = (LatLon( _0_0, lon1, **LatLon_kwds),
1787 LatLon( _90_0, lon1, **LatLon_kwds),
1788 LatLon( _0_0, lon2, **LatLon_kwds),
1789 LatLon(_N_90_0, lon2, **LatLon_kwds))
1790 if closed:
1791 t += t[:1]
1792 return t
1795def nearestOn5(point, points, closed=False, wrap=False, adjust=True,
1796 limit=9, **LatLon_and_kwds):
1797 '''Locate the point on a path or polygon closest to a reference point.
1799 The closest point on each polygon edge is either the nearest of that
1800 edge's end points or a point in between.
1802 @arg point: The reference point (C{LatLon}).
1803 @arg points: The path or polygon points (C{LatLon}[]).
1804 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1805 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1806 B{C{points}} (C{bool}).
1807 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1808 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1809 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1810 spherical earth radius L{R_KM}).
1811 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use for
1812 the closest point and additional B{C{LatLon}} keyword
1813 arguments, ignored if C{B{LatLon}=None} or not given.
1815 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1816 {closest} point (B{C{LatLon}}) or if C{B{LatLon} is None},
1817 a L{NearestOn5Tuple}C{(lat, lon, distance, angle, height)}.
1818 The C{distance} is the L{pygeodesy.equirectangular} distance
1819 between the C{closest} and reference B{C{point}} in C{degrees}.
1820 The C{angle} from the B{C{point}} to the C{closest} is in
1821 compass C{degrees}, like function L{pygeodesy.compassAngle}.
1823 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1824 see function L{pygeodesy.equirectangular_}.
1826 @raise PointsError: Insufficient number of B{C{points}}
1828 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1830 @note: Distances are I{approximated} by function L{pygeodesy.equirectangular_}.
1831 For more accuracy use one of the C{LatLon.nearestOn6} methods.
1833 @see: Function L{pygeodesy.degrees2m}.
1834 '''
1835 def _d2yx4(p2, p1, u, alw):
1836 # w = wrap if (i < (n - 1) or not closed) else False
1837 # equirectangular_ returns a Distance4Tuple(distance
1838 # in degrees squared, delta lat, delta lon, p2.lon
1839 # unroll/wrap'd); the previous p2.lon unroll/wrap'd
1840 # is also applied to the next edge's p1.lon
1841 return equirectangular_(p1.lat, p1.lon + u,
1842 p2.lat, p2.lon, **alw)
1844 def _h(p): # get height or default 0
1845 return _xattr(p, height=0) or 0
1847 # 3-D version used in .vector3d._nearestOn2
1848 #
1849 # point (x, y) on axis rotated ccw by angle a:
1850 # x' = y * sin(a) + x * cos(a)
1851 # y' = y * cos(a) - x * sin(a)
1852 #
1853 # distance (w) along and (h) perpendicular to
1854 # a line thru point (dx, dy) and the origin:
1855 # w = (y * dy + x * dx) / hypot(dx, dy)
1856 # h = (y * dx - x * dy) / hypot(dx, dy)
1857 #
1858 # closest point on that line thru (dx, dy):
1859 # xc = dx * w / hypot(dx, dy)
1860 # yc = dy * w / hypot(dx, dy)
1861 # or
1862 # xc = dx * f
1863 # yc = dy * f
1864 # with
1865 # f = w / hypot(dx, dy)
1866 # or
1867 # f = (y * dy + x * dx) / hypot2(dx, dy)
1868 #
1869 # i.e. no need for sqrt or hypot
1871 Ps = PointsIter(points, loop=1, wrap=wrap)
1872 p1 = c = Ps[0]
1873 u1 = u = _0_0
1874 kw = dict(adjust=adjust, limit=limit, wrap=False)
1875 d, dy, dx, _ = _d2yx4(p1, point, u1, kw)
1876 for p2 in Ps.iterate(closed=closed):
1877 # iff wrapped, unroll lon1 (actually previous
1878 # lon2) like function unroll180/-PI would've
1879 if wrap:
1880 kw.update(wrap=not (closed and Ps.looped))
1881 d21, y21, x21, u2 = _d2yx4(p2, p1, u1, kw)
1882 if d21 > EPS:
1883 # distance point to p1, y01 and x01 negated
1884 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)
1885 if d2 > EPS:
1886 w2 = y01 * y21 + x01 * x21
1887 if w2 > 0:
1888 if w2 < d21:
1889 # closest is between p1 and p2, use
1890 # original delta's, not y21 and x21
1891 f = w2 / d21
1892 p1 = LatLon_(favg(p1.lat, p2.lat, f=f),
1893 favg(p1.lon, p2.lon + u2, f=f),
1894 height=favg(_h(p1), _h(p2), f=f))
1895 u1 = _0_0
1896 else: # p2 is closest
1897 p1, u1 = p2, u2
1898 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)
1899 if d2 < d: # p1 is closer, y01 and x01 negated
1900 c, u, d, dy, dx = p1, u1, d2, -y01, -x01
1901 p1, u1 = p2, u2
1903 a = atan2b(dx, dy)
1904 d = hypot(dx, dy)
1905 h = _h(c)
1906 n = nameof(point) or nearestOn5.__name__
1907 if LatLon_and_kwds:
1908 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1909 LL = _xkwds_pop(kwds, LatLon=None)
1910 if LL is not None:
1911 r = LL(c.lat, c.lon + u, **kwds)
1912 return NearestOn3Tuple(r, d, a, name=n)
1913 return NearestOn5Tuple(c.lat, c.lon + u, d, a, h, name=n) # PYCHOK expected
1916def perimeterOf(points, closed=False, adjust=True, radius=R_M, wrap=True):
1917 '''I{Approximate} the perimeter of a path, polygon. or composite.
1919 @arg points: The path or polygon points or clips (C{LatLon}[],
1920 L{BooleanFHP} or L{BooleanGH}).
1921 @kwarg closed: Optionally, close the path or polygon (C{bool}).
1922 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta
1923 by the cosine of the mean latitude (C{bool}).
1924 @kwarg radius: Mean earth radius (C{meter}).
1925 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1926 B{C{points}} (C{bool}).
1928 @return: Approximate perimeter (C{meter}, same units as
1929 B{C{radius}}).
1931 @raise PointsError: Insufficient number of B{C{points}}
1933 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1935 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1936 C{B{points}} a composite.
1938 @note: This perimeter is based on the L{pygeodesy.equirectangular_}
1939 distance approximation and is ill-suited for regions exceeding
1940 several hundred Km or Miles or with near-polar latitudes.
1942 @see: Functions L{sphericalTrigonometry.perimeterOf} and
1943 L{ellipsoidalKarney.perimeterOf}.
1944 '''
1945 def _degs(ps, c, a, w): # angular edge lengths in degrees
1946 Ps = LatLon2PsxyIter(ps, loop=1) # wrap=w
1947 p1, u = Ps[0], _0_0 # previous x2's unroll/wrap
1948 for p2 in Ps.iterate(closed=c):
1949 if w and c:
1950 w = not Ps.looped
1951 # apply previous x2's unroll/wrap'd to new x1
1952 _, dy, dx, u = equirectangular_(p1.y, p1.x + u,
1953 p2.y, p2.x,
1954 adjust=a, limit=None,
1955 wrap=w) # PYCHOK non-seq
1956 yield hypot(dx, dy)
1957 p1 = p2
1959 if _MODS.booleans.isBoolean(points):
1960 if not closed:
1961 notImplemented(None, closed=closed, points=_composite_)
1962 d = points._sum1(perimeterOf, closed=True, adjust=adjust,
1963 radius=radius, wrap=wrap)
1964 else:
1965 d = fsum(_degs(points, closed, adjust, wrap), floats=True)
1966 return degrees2m(d, radius=radius)
1969def quadOf(latS, lonW, latN, lonE, closed=False, LatLon=LatLon_, **LatLon_kwds):
1970 '''Generate a quadrilateral path or polygon from two points.
1972 @arg latS: Souther-nmost latitude (C{degrees90}).
1973 @arg lonW: Western-most longitude (C{degrees180}).
1974 @arg latN: Norther-nmost latitude (C{degrees90}).
1975 @arg lonE: Eastern-most longitude (C{degrees180}).
1976 @kwarg closed: Optionally, close the path (C{bool}).
1977 @kwarg LatLon: Class to use (L{LatLon_}).
1978 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}
1979 keyword arguments.
1981 @return: Return a tuple of 4 or 5 B{C{LatLon}} instances
1982 outlining the quadrilateral.
1984 @see: Function L{boundsOf}.
1985 '''
1986 t = (LatLon(latS, lonW, **LatLon_kwds),
1987 LatLon(latN, lonW, **LatLon_kwds),
1988 LatLon(latN, lonE, **LatLon_kwds),
1989 LatLon(latS, lonE, **LatLon_kwds))
1990 if closed:
1991 t += t[:1]
1992 return t
1995__all__ += _ALL_DOCS(_Array2LatLon, _Basequence)
1997# **) MIT License
1998#
1999# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
2000#
2001# Permission is hereby granted, free of charge, to any person obtaining a
2002# copy of this software and associated documentation files (the "Software"),
2003# to deal in the Software without restriction, including without limitation
2004# the rights to use, copy, modify, merge, publish, distribute, sublicense,
2005# and/or sell copies of the Software, and to permit persons to whom the
2006# Software is furnished to do so, subject to the following conditions:
2007#
2008# The above copyright notice and this permission notice shall be included
2009# in all copies or substantial portions of the Software.
2010#
2011# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
2012# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
2013# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
2014# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
2015# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
2016# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
2017# OTHER DEALINGS IN THE SOFTWARE.