1

2# -*- coding: utf-8 -*-

3

4 u'''Utilities for point lists, tuples, etc.

5

6 Functions to handle collections and sequences of C{LatLon} points

7 specified as 2-d U{NumPy<https://www.NumPy.org>}, C{arrays} or tuples

8 as C{LatLon} or as C{pseudo-x/-y} pairs.

9

10 C{NumPy} arrays are assumed to contain rows of points with a lat-, a

11 longitude -and possibly other- values in different columns. While

12 iterating over the array rows, create an instance of a given C{LatLon}

13 class "on-the-fly" for each row with the row's lat- and longitude.

14

15 The original C{NumPy} array is read-accessed only and never duplicated,

16 except to return a I{subset} of the original array.

17

18 For example, to process a C{NumPy} array, wrap the array by instantiating

19 class L{Numpy2LatLon} and specifying the column index for the lat- and

20 longitude in each row. Then, pass the L{Numpy2LatLon} instance to any

21 L{pygeodesy} function or method accepting a I{points} argument.

22

23 Similarly, class L{Tuple2LatLon} is used to instantiate a C{LatLon} from

24 each 2+tuple in a sequence of such 2+tuples using the C{ilat} lat- and

25 C{ilon} longitude index in each 2+tuple.

26'''

27

28 from pygeodesy.basics import isclass, isint, isscalar, issequence, \

29 issubclassof, _Sequence, _xcopy, _xdup, \

30 _xinstanceof

31 from 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

35 from pygeodesy.dms import F_D, parseDMS

36 from pygeodesy.errors import CrossError, crosserrors, _IndexError, \

37 _IsnotError, _TypeError, _ValueError, \

38 _xattr, _xkwds, _xkwds_pop

39 from pygeodesy.fmath import favg, fdot, hypot, Fsum, fsum

40# from pygeodesy.fsums import Fsum, fsum # from .fmath

41 from pygeodesy.formy import _bearingTo2, equirectangular_, _spherical_datum

42 from pygeodesy.interns import NN, _colinear_, _COMMASPACE_, _composite_, \

43 _DEQUALSPACED_, _ELLIPSIS_, _EW_, _immutable_, \

44 _near_, _no_, _not_, _NS_, _point_, _SPACE_, \

45 _UNDER_, _valid_ # _lat_, _lon_

46 from pygeodesy.iters import LatLon2PsxyIter, PointsIter, points2

47 from pygeodesy.latlonBase import LatLonBase, _latlonheight3, \

48 _ALL_DOCS, _ALL_LAZY, _MODS

49# from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS

50 from pygeodesy.named import classname, nameof, notImplemented, notOverloaded, \

51 _NamedTuple

52 from pygeodesy.namedTuples import Bounds2Tuple, Bounds4Tuple, LatLon2Tuple, \

53 NearestOn3Tuple, NearestOn5Tuple, \

54 Point3Tuple, Vector3Tuple, \

55 PhiLam2Tuple # PYCHOK shared

56 from pygeodesy.props import Property_RO, property_doc_, property_RO

57 from pygeodesy.streprs import Fmt, instr

58 from pygeodesy.units import Number_, Radius, Scalar, Scalar_

59 from pygeodesy.utily import atan2b, degrees90, degrees180, degrees2m, \

60 unroll180, _unrollon, unrollPI, _Wrap, wrap180

61

62 from math import cos, fabs, fmod, radians, sin

63

64__all__ = _ALL_LAZY.points

65__version__ = '23.10.07'

66

67_ilat_ = 'ilat'

68_ilon_ = 'ilon'

69_ncols_ = 'ncols'

70_nrows_ = 'nrows'

71

72

73 class LatLon_(LatLonBase): # XXX in heights._HeightBase.height

74 '''Low-overhead C{LatLon} class, mainly for L{Numpy2LatLon} and L{Tuple2LatLon}.

75 '''

76 # __slots__ efficiency is voided if the __slots__ class attribute is

77 # used in a subclass of a class with the traditional __dict__, @see

78 # <https://docs.Python.org/2/reference/datamodel.html#slots> plus ...

79 #

80 # __slots__ must be repeated in sub-classes, @see Luciano Ramalho,

81 # "Fluent Python", O'Reilly, 2016 p. 276+ "Problems with __slots__",

82 # 2nd Ed, 2022 p. 390 "Summarizing the Issues with __slots__".

83 #

84 # __slots__ = (_lat_, _lon_, _height_, _datum_, _name_)

85 # Property_RO = property_RO # no __dict__ with __slots__!

86 #

87 # In addition, both size and overhead have shrunk in recent Python:

88 #

89 # sys.getsizeof(LatLon_(1, 2)) is 72-88 I{with} __slots__, but

90 # only 48-56 bytes I{without in Python 2.7.18+ and Python 3+}.

91 #

92 # python3 -m timeit -s "from pygeodesy... import LatLonBase as LL" "LL(0, 0)" 2.14 usec

93 # python3 -m timeit -s "from pygeodesy import LatLon_" "LatLon_(0, 0)" 216 nsec

94

95 def __init__(self, latlonh, lon=None, height=0, wrap=False, name=NN, datum=None):

96 '''New L{LatLon_}.

97

98 @note: The lat- and longitude values are taken I{as-given,

99 un-clipped and un-validated}.

100

101 @see: L{latlonBase.LatLonBase} for further details.

102 '''

103 if name:

104 self.name = name

105

106 if lon is None: # PYCHOK no cover

107 lat, lon, height = _latlonheight3(latlonh, height, wrap)

108 elif wrap: # PYCHOK no cover

109 lat, lon = _Wrap.latlonDMS2(latlonh, lon)

110 else: # must be latNS, lonEW

111 try:

112 lat, lon = float(latlonh), float(lon)

113 except (TypeError, ValueError):

114 lat = parseDMS(latlonh, suffix=_NS_)

115 lon = parseDMS(lon, suffix=_EW_)

116

117 # get the minimal __dict__, see _isLatLon_ below

118 self._lat = lat # un-clipped and ...

119 self._lon = lon # ... un-validated

120 self._datum = None if datum is None else \

121 _spherical_datum(datum, name=self.name)

122 self._height = height

123

124 def __eq__(self, other):

125 return isinstance(other, LatLon_) and \

126 other.lat == self.lat and \

127 other.lon == self.lon

128

129 def __ne__(self, other):

130 return not self.__eq__(other)

131

132 @Property_RO

133 def datum(self):

134 '''Get the C{datum} (L{Datum}) or C{None}.

135 '''

136 return self._datum

137

138 def intermediateTo(self, other, fraction, height=None, wrap=False):

139 '''Locate the point at a given fraction, I{linearly} between

140 (or along) this and an other point.

141

142 @arg other: The other point (C{LatLon}).

143 @arg fraction: Fraction between both points (C{float},

144 0.0 for this and 1.0 for the other point).

145 @kwarg height: Optional height (C{meter}), overriding the

146 intermediate height.

147 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

148 the B{C{other}} point (C{bool}).

149

150 @return: Intermediate point (this C{LatLon}).

151

152 @raise TypeError: Incompatible B{C{other}} C{type}.

153 '''

154 f = Scalar(fraction=fraction)

155 if isnear0(f):

156 r = self

157 else:

158 r = self.others(other)

159 if wrap or not isnear1(f):

160 _, lat, lon = _Wrap.latlon3(self.lon, r.lat, r.lon, wrap)

161 lat = favg(self.lat, lat, f=f)

162 lon = favg(self.lon, lon, f=f)

163 h = height if height is not None else \

164 favg(self.height, r.height, f=f)

165 # = self._havg(r, f=f, h=height)

166 r = self.classof(lat, lon, height=h, datum=r.datum,

167 name=r.intermediateTo.__name__)

168 return r

169

170 def toRepr(self, **kwds):

171 '''This L{LatLon_} as a string "class(<degrees>, ...)",

172 ignoring keyword argument C{B{std}=N/A}.

173

174 @see: L{latlonBase.LatLonBase.toRepr} for further details.

175 '''

176 _ = _xkwds_pop(kwds, std=NotImplemented)

177 return LatLonBase.toRepr(self, **kwds)

178

179 def toStr(self, form=F_D, joined=_COMMASPACE_, **m_prec_sep_s_D_M_S): # PYCHOK expected

180 '''Convert this point to a "lat, lon[, height][, name][, ...]"

181 string, formatted in the given C{B{form}at}.

182

183 @see: L{latlonBase.LatLonBase.toStr} for further details.

184 '''

185 t = LatLonBase.toStr(self, form=form, joined=NN,

186 **_xkwds(m_prec_sep_s_D_M_S, m=NN))

187 if self.name:

188 t += (repr(self.name),)

189 return joined.join(t) if joined else t

190

191

192 def _isLatLon(inst):

193 '''(INTERNAL) Check a C{LatLon} or C{LatLon_} instance.

194 '''

195 return isinstance(inst, (LatLon_, _MODS.latlonBase.LatLonBase))

196

197

198 def _isLatLon_(LL):

199 '''(INTERNAL) Check a (sub-)class of C{LatLon_}.

200 '''

201 return issubclassof(LL, LatLon_) or (isclass(LL) and

202 all(hasattr(LL, _) for _ in LatLon_(0, 0).__dict__.keys()))

203

204

205 class _Basequence(_Sequence): # immutable, on purpose

206 '''(INTERNAL) Base class.

207 '''

208 _array = []

209 _epsilon = EPS

210 _itemname = _point_

211

212 def _contains(self, point):

213 '''(INTERNAL) Check for a matching point.

214 '''

215 return any(self._findall(point, ()))

216

217 def copy(self, deep=False): # PYCHOK no cover

218 '''Make a shallow or deep copy of this instance.

219

220 @kwarg deep: If C{True} make a deep, otherwise a

221 shallow copy (C{bool}).

222

223 @return: The copy (C{This class} or subclass thereof).

224 '''

225 return _xcopy(self, deep=deep)

226

227 def _count(self, point):

228 '''(INTERNAL) Count the number of matching points.

229 '''

230 return sum(1 for _ in self._findall(point, ())) # NOT len()!

231

232 def dup(self, **items): # PYCHOK no cover

233 '''Duplicate this instance, I{without replacing items}.

234

235 @kwarg items: No attributes (I{not allowed}).

236

237 @return: The duplicate (C{This} (sub-)class).

238

239 @raise TypeError: Any B{C{items}} invalid.

240 '''

241 if items:

242 t = _SPACE_(classname(self), _immutable_)

243 raise _TypeError(txt=t, this=self, **items)

244 return _xdup(self)

245

246 @property_doc_(''' the equality tolerance (C{float}).''')

247 def epsilon(self):

248 '''Get the tolerance for equality tests (C{float}).

249 '''

250 return self._epsilon

251

252 @epsilon.setter # PYCHOK setter!

253 def epsilon(self, tol):

254 '''Set the tolerance for equality tests (C{scalar}).

255

256 @raise UnitError: Non-scalar or invalid B{C{tol}}.

257 '''

258 self._epsilon = Scalar_(tolerance=tol)

259

260 def _find(self, point, start_end):

261 '''(INTERNAL) Find the first matching point index.

262 '''

263 for i in self._findall(point, start_end):

264 return i

265 return -1

266

267 def _findall(self, point, start_end): # PYCHOK no cover

268 '''(INTERNAL) I{Must be implemented/overloaded}.

269 '''

270 notImplemented(self, point, start_end)

271

272 def _getitem(self, index):

273 '''(INTERNAL) Return point [index] or return a slice.

274 '''

275 # Luciano Ramalho, "Fluent Python", O'Reilly, 2016 p. 290+, 2022 p. 405+

276 if isinstance(index, slice):

277 # XXX an numpy.[nd]array slice is a view, not a copy

278 return self.__class__(self._array[index], **self._slicekwds())

279 else:

280 return self.point(self._array[index])

281

282 def _index(self, point, start_end):

283 '''(INTERNAL) Find the first matching point index.

284 '''

285 for i in self._findall(point, start_end):

286 return i

287 raise _IndexError(self._itemname, point, txt=_not_('found'))

288

289 @property_RO

290 def isNumpy2(self): # PYCHOK no cover

291 '''Is this a Numpy2 wrapper?

292 '''

293 return False # isinstance(self, (Numpy2LatLon, ...))

294

295 @property_RO

296 def isPoints2(self): # PYCHOK no cover

297 '''Is this a LatLon2 wrapper/converter?

298 '''

299 return False # isinstance(self, (LatLon2psxy, ...))

300

301 @property_RO

302 def isTuple2(self): # PYCHOK no cover

303 '''Is this a Tuple2 wrapper?

304 '''

305 return False # isinstance(self, (Tuple2LatLon, ...))

306

307 def _iter(self):

308 '''(INTERNAL) Yield all points.

309 '''

310 _array, _point = self._array, self.point

311 for i in range(len(self)):

312 yield _point(_array[i])

313

314 def point(self, *attrs): # PYCHOK no cover

315 '''I{Must be overloaded}.

316

317 @arg attrs: Optional arguments.

318 '''

319 notOverloaded(self, *attrs)

320

321 def _range(self, start=None, end=None, step=1):

322 '''(INTERNAL) Return the range.

323 '''

324 if step > 0:

325 if start is None:

326 start = 0

327 if end is None:

328 end = len(self)

329 elif step < 0:

330 if start is None:

331 start = len(self) - 1

332 if end is None:

333 end = -1

334 else:

335 raise _ValueError(step=step)

336 return range(start, end, step)

337

338 def _repr(self):

339 '''(INTERNAL) Return a string representation.

340 '''

341 # XXX use Python 3+ reprlib.repr

342 t = repr(self._array[:1]) # only first row

343 t = _SPACE_(t[:-1], _ELLIPSIS_, Fmt.SQUARE(t[-1:], len(self)))

344 t = _SPACE_.join(t.split()) # coalesce spaces

345 return instr(self, t, **self._slicekwds())

346

347 def _reversed(self): # PYCHOK false

348 '''(INTERNAL) Yield all points in reverse order.

349 '''

350 _array, point = self._array, self.point

351 for i in range(len(self) - 1, -1, -1):

352 yield point(_array[i])

353

354 def _rfind(self, point, start_end):

355 '''(INTERNAL) Find the last matching point index.

356 '''

357 def _r3(start=None, end=None, step=-1):

358 return (start, end, step) # PYCHOK returns

359

360 for i in self._findall(point, _r3(*start_end)):

361 return i

362 return -1

363

364 def _slicekwds(self): # PYCHOK no cover

365 '''(INTERNAL) I{Should be overloaded}.

366 '''

367 return {}

368

369

370 class _Array2LatLon(_Basequence): # immutable, on purpose

371 '''(INTERNAL) Base class for Numpy2LatLon or Tuple2LatLon.

372 '''

373 _array = ()

374 _ilat = 0 # row column index

375 _ilon = 0 # row column index

376 _LatLon = LatLon_ # default

377 _shape = ()

378

379 def __init__(self, array, ilat=0, ilon=1, LatLon=None, shape=()):

380 '''Handle a C{NumPy} or C{Tuple} array as a sequence of C{LatLon} points.

381 '''

382 ais = (_ilat_, ilat), (_ilon_, ilon)

383

384 if len(shape) != 2 or shape[0] < 1 or shape[1] < len(ais):

385 raise _IndexError('array.shape', shape)

386

387 self._array = array

388 self._shape = Shape2Tuple(shape) # *shape

389

390 if LatLon: # check the point class

391 if not _isLatLon_(LatLon):

392 raise _IsnotError(_valid_, LatLon=LatLon)

393 self._LatLon = LatLon

394

395 # check the attr indices

396 for n, (ai, i) in enumerate(ais):

397 if not isint(i):

398 raise _IsnotError(int.__name__, **{ai: i})

399 i = int(i)

400 if not 0 <= i < shape[1]:

401 raise _ValueError(ai, i)

402 for aj, j in ais[:n]:

403 if int(j) == i:

404 raise _ValueError(_DEQUALSPACED_(ai, aj, i))

405 setattr(self, NN(_UNDER_, ai), i)

406

407 def __contains__(self, latlon):

408 '''Check for a specific lat-/longitude.

409

410 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

411 C{(lat, lon)}).

412

413 @return: C{True} if B{C{latlon}} is present, C{False} otherwise.

414

415 @raise TypeError: Invalid B{C{latlon}}.

416 '''

417 return self._contains(latlon)

418

419 def __getitem__(self, index):

420 '''Return row[index] as C{LatLon} or return a L{Numpy2LatLon} slice.

421 '''

422 return self._getitem(index)

423

424 def __iter__(self):

425 '''Yield rows as C{LatLon}.

426 '''

427 return self._iter()

428

429 def __len__(self):

430 '''Return the number of rows.

431 '''

432 return self._shape[0]

433

434 def __repr__(self):

435 '''Return a string representation.

436 '''

437 return self._repr()

438

439 def __reversed__(self): # PYCHOK false

440 '''Yield rows as C{LatLon} in reverse order.

441 '''

442 return self._reversed()

443

444 __str__ = __repr__

445

446 def count(self, latlon):

447 '''Count the number of rows with a specific lat-/longitude.

448

449 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

450 C{(lat, lon)}).

451

452 @return: Count (C{int}).

453

454 @raise TypeError: Invalid B{C{latlon}}.

455 '''

456 return self._count(latlon)

457

458 def find(self, latlon, *start_end):

459 '''Find the first row with a specific lat-/longitude.

460

461 @arg latlon: Point (C{LatLon}) or 2-tuple (lat, lon).

462 @arg start_end: Optional C{[start[, end]]} index (integers).

463

464 @return: Index or -1 if not found (C{int}).

465

466 @raise TypeError: Invalid B{C{latlon}}.

467 '''

468 return self._find(latlon, start_end)

469

470 def _findall(self, latlon, start_end):

471 '''(INTERNAL) Yield indices of all matching rows.

472 '''

473 try:

474 lat, lon = latlon.lat, latlon.lon

475 except AttributeError:

476 try:

477 lat, lon = latlon

478 except (TypeError, ValueError):

479 raise _IsnotError(_valid_, latlon=latlon)

480

481 _ilat, _ilon = self._ilat, self._ilon

482 _array, _eps = self._array, self._epsilon

483 for i in self._range(*start_end):

484 row = _array[i]

485 if fabs(row[_ilat] - lat) <= _eps and \

486 fabs(row[_ilon] - lon) <= _eps:

487 yield i

488

489 def findall(self, latlon, *start_end):

490 '''Yield indices of all rows with a specific lat-/longitude.

491

492 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

493 C{(lat, lon)}).

494 @arg start_end: Optional C{[start[, end]]} index (C{int}).

495

496 @return: Indices (C{iterable}).

497

498 @raise TypeError: Invalid B{C{latlon}}.

499 '''

500 return self._findall(latlon, start_end)

501

502 def index(self, latlon, *start_end): # PYCHOK Python 2- issue

503 '''Find index of the first row with a specific lat-/longitude.

504

505 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

506 C{(lat, lon)}).

507 @arg start_end: Optional C{[start[, end]]} index (C{int}).

508

509 @return: Index (C{int}).

510

511 @raise IndexError: Point not found.

512

513 @raise TypeError: Invalid B{C{latlon}}.

514 '''

515 return self._index(latlon, start_end)

516

517 @Property_RO

518 def ilat(self):

519 '''Get the latitudes column index (C{int}).

520 '''

521 return self._ilat

522

523 @Property_RO

524 def ilon(self):

525 '''Get the longitudes column index (C{int}).

526 '''

527 return self._ilon

528

529# next = __iter__

530

531 def point(self, row): # PYCHOK *attrs

532 '''Instantiate a point C{LatLon}.

533

534 @arg row: Array row (numpy.array).

535

536 @return: Point (C{LatLon}).

537 '''

538 return self._LatLon(row[self._ilat], row[self._ilon])

539

540 def rfind(self, latlon, *start_end):

541 '''Find the last row with a specific lat-/longitude.

542

543 @arg latlon: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

544 C{(lat, lon)}).

545 @arg start_end: Optional C{[start[, end]]} index (C{int}).

546

547 @note: Keyword order, first stop, then start.

548

549 @return: Index or -1 if not found (C{int}).

550

551 @raise TypeError: Invalid B{C{latlon}}.

552 '''

553 return self._rfind(latlon, start_end)

554

555 def _slicekwds(self):

556 '''(INTERNAL) Slice kwds.

557 '''

558 return dict(ilat=self._ilat, ilon=self._ilon)

559

560 @Property_RO

561 def shape(self):

562 '''Get the shape of the C{NumPy} array or the C{Tuples} as

563 L{Shape2Tuple}C{(nrows, ncols)}.

564 '''

565 return self._shape

566

567 def _subset(self, indices): # PYCHOK no cover

568 '''(INTERNAL) I{Must be implemented/overloaded}.

569 '''

570 notImplemented(self, indices)

571

572 def subset(self, indices):

573 '''Return a subset of the C{NumPy} array.

574

575 @arg indices: Row indices (C{range} or C{int}[]).

576

577 @note: A C{subset} is different from a C{slice} in 2 ways:

578 (a) the C{subset} is typically specified as a list of

579 (un-)ordered indices and (b) the C{subset} allocates

580 a new, separate C{NumPy} array while a C{slice} is

581 just an other C{view} of the original C{NumPy} array.

582

583 @return: Sub-array (C{numpy.array}).

584

585 @raise IndexError: Out-of-range B{C{indices}} value.

586

587 @raise TypeError: If B{C{indices}} is not a C{range}

588 nor an C{int}[].

589 '''

590 if not issequence(indices, tuple): # NO tuple, only list

591 # and range work properly to get Numpy array sub-sets

592 raise _IsnotError(_valid_, indices=type(indices))

593

594 n = len(self)

595 for i, v in enumerate(indices):

596 if not isint(v):

597 raise _TypeError(Fmt.SQUARE(indices=i), v)

598 elif not 0 <= v < n:

599 raise _IndexError(Fmt.SQUARE(indices=i), v)

600

601 return self._subset(indices)

602

603

604 class LatLon2psxy(_Basequence):

605 '''Wrapper for C{LatLon} points as "on-the-fly" pseudo-xy coordinates.

606 '''

607 _closed = False

608 _len = 0

609 _deg2m = None # default, keep degrees

610 _radius = None

611 _wrap = True

612

613 def __init__(self, latlons, closed=False, radius=None, wrap=True):

614 '''Handle C{LatLon} points as pseudo-xy coordinates.

615

616 @note: The C{LatLon} latitude is considered the I{pseudo-y}

617 and longitude the I{pseudo-x} coordinate, likewise

618 for L{LatLon2Tuple}. However, 2-tuples C{(x, y)} are

619 considered as I{(longitude, latitude)}.

620

621 @arg latlons: Points C{list}, C{sequence}, C{set}, C{tuple},

622 etc. (C{LatLon[]}).

623 @kwarg closed: Optionally, close the polygon (C{bool}).

624 @kwarg radius: Mean earth radius (C{meter}).

625 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

626 the B{C{latlons}} points (C{bool}).

627

628 @raise PointsError: Insufficient number of B{C{latlons}}.

629

630 @raise TypeError: Some B{C{points}} are not B{C{base}}.

631 '''

632 self._closed = closed

633 self._len, self._array = points2(latlons, closed=closed)

634 if radius:

635 self._radius = r = Radius(radius)

636 self._deg2m = degrees2m(_1_0, r)

637 if not wrap:

638 self._wrap = False

639

640 def __contains__(self, xy):

641 '''Check for a matching point.

642

643 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

644 C{(x, y)}) in (C{degrees}.

645

646 @return: C{True} if B{C{xy}} is present, C{False} otherwise.

647

648 @raise TypeError: Invalid B{C{xy}}.

649 '''

650 return self._contains(xy)

651

652 def __getitem__(self, index):

653 '''Return the pseudo-xy or return a L{LatLon2psxy} slice.

654 '''

655 return self._getitem(index)

656

657 def __iter__(self):

658 '''Yield all pseudo-xy's.

659 '''

660 return self._iter()

661

662 def __len__(self):

663 '''Return the number of pseudo-xy's.

664 '''

665 return self._len

666

667 def __repr__(self):

668 '''Return a string representation.

669 '''

670 return self._repr()

671

672 def __reversed__(self): # PYCHOK false

673 '''Yield all pseudo-xy's in reverse order.

674 '''

675 return self._reversed()

676

677 __str__ = __repr__

678

679 def count(self, xy):

680 '''Count the number of matching points.

681

682 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

683 C{(x, y)}) in (C{degrees}.

684

685 @return: Count (C{int}).

686

687 @raise TypeError: Invalid B{C{xy}}.

688 '''

689 return self._count(xy)

690

691 def find(self, xy, *start_end):

692 '''Find the first matching point.

693

694 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

695 C{(x, y)}) in (C{degrees}.

696 @arg start_end: Optional C{[start[, end]]} index (C{int}).

697

698 @return: Index or -1 if not found (C{int}).

699

700 @raise TypeError: Invalid B{C{xy}}.

701 '''

702 return self._find(xy, start_end)

703

704 def _findall(self, xy, start_end):

705 '''(INTERNAL) Yield indices of all matching points.

706 '''

707 try:

708 x, y = xy.lon, xy.lat

709

710 def _x_y_ll3(ll): # match LatLon

711 return ll.lon, ll.lat, ll

712

713 except AttributeError:

714 try:

715 x, y = xy[:2]

716 except (IndexError, TypeError, ValueError):

717 raise _IsnotError(_valid_, xy=xy)

718

719 _x_y_ll3 = self.point # PYCHOK expected

720

721 _array, _eps = self._array, self._epsilon

722 for i in self._range(*start_end):

723 xi, yi, _ = _x_y_ll3(_array[i])

724 if fabs(xi - x) <= _eps and \

725 fabs(yi - y) <= _eps:

726 yield i

727

728 def findall(self, xy, *start_end):

729 '''Yield indices of all matching points.

730

731 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

732 C{(x, y)}) in (C{degrees}.

733 @arg start_end: Optional C{[start[, end]]} index (C{int}).

734

735 @return: Indices (C{iterator}).

736

737 @raise TypeError: Invalid B{C{xy}}.

738 '''

739 return self._findall(xy, start_end)

740

741 def index(self, xy, *start_end): # PYCHOK Python 2- issue

742 '''Find the first matching point.

743

744 @arg xy: Point (C{LatLon}) or 2-tuple (x, y) in (C{degrees}).

745 @arg start_end: Optional C{[start[, end]]} index (C{int}).

746

747 @return: Index (C{int}).

748

749 @raise IndexError: Point not found.

750

751 @raise TypeError: Invalid B{C{xy}}.

752 '''

753 return self._index(xy, start_end)

754

755 @property_RO

756 def isPoints2(self):

757 '''Is this a LatLon2 wrapper/converter?

758 '''

759 return True # isinstance(self, (LatLon2psxy, ...))

760

761 def point(self, ll): # PYCHOK *attrs

762 '''Create a pseudo-xy.

763

764 @arg ll: Point (C{LatLon}).

765

766 @return: An L{Point3Tuple}C{(x, y, ll)}.

767 '''

768 x, y = ll.lon, ll.lat # note, x, y = lon, lat

769 if self._wrap:

770 y, x = _Wrap.latlon(y, x)

771 d = self._deg2m

772 if d: # convert degrees to meter (or radians)

773 x *= d

774 y *= d

775 return Point3Tuple(x, y, ll)

776

777 def rfind(self, xy, *start_end):

778 '''Find the last matching point.

779

780 @arg xy: Point (C{LatLon}, L{LatLon2Tuple} or 2-tuple

781 C{(x, y)}) in (C{degrees}.

782 @arg start_end: Optional C{[start[, end]]} index (C{int}).

783

784 @return: Index or -1 if not found (C{int}).

785

786 @raise TypeError: Invalid B{C{xy}}.

787 '''

788 return self._rfind(xy, start_end)

789

790 def _slicekwds(self):

791 '''(INTERNAL) Slice kwds.

792 '''

793 return dict(closed=self._closed, radius=self._radius, wrap=self._wrap)

794

795

796 class Numpy2LatLon(_Array2LatLon): # immutable, on purpose

797 '''Wrapper for C{NumPy} arrays as "on-the-fly" C{LatLon} points.

798 '''

799 def __init__(self, array, ilat=0, ilon=1, LatLon=None):

800 '''Handle a C{NumPy} array as a sequence of C{LatLon} points.

801

802 @arg array: C{NumPy} array (C{numpy.array}).

803 @kwarg ilat: Optional index of the latitudes column (C{int}).

804 @kwarg ilon: Optional index of the longitudes column (C{int}).

805 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).

806

807 @raise IndexError: If B{C{array.shape}} is not (1+, 2+).

808

809 @raise TypeError: If B{C{array}} is not a C{NumPy} array or

810 C{LatLon} is not a class with C{lat}

811 and C{lon} attributes.

812

813 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values

814 are the same or out of range.

815

816 @example:

817

818 >>> type(array)

819 <type 'numpy.ndarray'> # <class ...> in Python 3+

820 >>> points = Numpy2LatLon(array, lat=0, lon=1)

821 >>> simply = simplifyRDP(points, ...)

822 >>> type(simply)

823 <type 'numpy.ndarray'> # <class ...> in Python 3+

824 >>> sliced = points[1:-1]

825 >>> type(sliced)

826 <class '...Numpy2LatLon'>

827 '''

828 try: # get shape and check some other numpy.array attrs

829 s, _, _ = array.shape, array.nbytes, array.ndim # PYCHOK expected

830 except AttributeError:

831 raise _IsnotError('NumPy', array=type(array))

832

833 _Array2LatLon.__init__(self, array, ilat=ilat, ilon=ilon,

834 LatLon=LatLon, shape=s)

835

836 @property_RO

837 def isNumpy2(self):

838 '''Is this a Numpy2 wrapper?

839 '''

840 return True # isinstance(self, (Numpy2LatLon, ...))

841

842 def _subset(self, indices):

843 return self._array[indices] # NumPy special

844

845

846 class Shape2Tuple(_NamedTuple):

847 '''2-Tuple C{(nrows, ncols)}, the number of rows and columns,

848 both C{int}.

849 '''

850 _Names_ = (_nrows_, _ncols_)

851 _Units_ = ( Number_, Number_)

852

853

854 class Tuple2LatLon(_Array2LatLon):

855 '''Wrapper for tuple sequences as "on-the-fly" C{LatLon} points.

856 '''

857 def __init__(self, tuples, ilat=0, ilon=1, LatLon=None):

858 '''Handle a list of tuples, each containing a lat- and longitude

859 and perhaps other values as a sequence of C{LatLon} points.

860

861 @arg tuples: The C{list}, C{tuple} or C{sequence} of tuples (C{tuple}[]).

862 @kwarg ilat: Optional index of the latitudes value (C{int}).

863 @kwarg ilon: Optional index of the longitudes value (C{int}).

864 @kwarg LatLon: Optional C{LatLon} class to use (L{LatLon_}).

865

866 @raise IndexError: If C{(len(B{tuples}), min(len(t) for t

867 in B{tuples}))} is not (1+, 2+).

868

869 @raise TypeError: If B{C{tuples}} is not a C{list}, C{tuple}

870 or C{sequence} or if B{C{LatLon}} is not a

871 C{LatLon} with C{lat}, C{lon} and C{name}

872 attributes.

873

874 @raise ValueError: If the B{C{ilat}} and/or B{C{ilon}} values

875 are the same or out of range.

876

877 @example:

878

879 >>> tuples = [(0, 1), (2, 3), (4, 5)]

880 >>> type(tuples)

881 <type 'list'> # <class ...> in Python 3+

882 >>> points = Tuple2LatLon(tuples, lat=0, lon=1)

883 >>> simply = simplifyRW(points, 0.5, ...)

884 >>> type(simply)

885 <type 'list'> # <class ...> in Python 3+

886 >>> simply

887 [(0, 1), (4, 5)]

888 >>> sliced = points[1:-1]

889 >>> type(sliced)

890 <class '...Tuple2LatLon'>

891 >>> sliced

892 ...Tuple2LatLon([(2, 3), ...][1], ilat=0, ilon=1)

893

894 >>> closest, _ = nearestOn2(LatLon_(2, 1), points, adjust=False)

895 >>> closest

896 LatLon_(lat=1.0, lon=2.0)

897

898 >>> closest, _ = nearestOn2(LatLon_(3, 2), points)

899 >>> closest

900 LatLon_(lat=2.001162, lon=3.001162)

901 '''

902 _xinstanceof(list, tuple, tuples=tuples)

903 s = len(tuples), min(len(_) for _ in tuples)

904 _Array2LatLon.__init__(self, tuples, ilat=ilat, ilon=ilon,

905 LatLon=LatLon, shape=s)

906

907 @property_RO

908 def isTuple2(self):

909 '''Is this a Tuple2 wrapper?

910 '''

911 return True # isinstance(self, (Tuple2LatLon, ...))

912

913 def _subset(self, indices):

914 return type(self._array)(self._array[i] for i in indices)

915

916

917 def _area2(points, adjust, wrap):

918 '''(INTERNAL) Approximate the area in radians squared, I{signed}.

919 '''

920 if adjust:

921 # approximate trapezoid by a rectangle, adjusting

922 # the top width by the cosine of the latitudinal

923 # average and bottom width by some fudge factor

924 def _adjust(w, h):

925 c = cos(h) if fabs(h) < PI_2 else _0_0

926 return w * h * (c + 1.2876) * _0_5

927 else:

928 def _adjust(w, h): # PYCHOK expected

929 return w * h

930

931 # setting radius=1 converts degrees to radians

932 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)

933 x1, y1, ll = Ps[0]

934 pts = [ll] # for _areaError

935

936 A2 = Fsum() # trapezoidal area in radians**2

937 for p in Ps.iterate(closed=True):

938 x2, y2, ll = p

939 if len(pts) < 4:

940 pts.append(ll)

941 w, x2 = unrollPI(x1, x2, wrap=wrap and not Ps.looped)

942 A2 += _adjust(w, (y2 + y1) * _0_5)

943 x1, y1 = x2, y2

944

945 return A2.fsum(), tuple(pts)

946

947

948 def _areaError(pts, near_=NN): # imported by .ellipsoidalKarney

949 '''(INTERNAL) Area issue.

950 '''

951 t = _ELLIPSIS_(pts[:3], NN)

952 return _ValueError(NN(near_, 'zero or polar area'), txt=t)

953

954

955 def areaOf(points, adjust=True, radius=R_M, wrap=True):

956 '''Approximate the area of a polygon or composite.

957

958 @arg points: The polygon points or clips (C{LatLon}[],

959 L{BooleanFHP} or L{BooleanGH}).

960 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta

961 by the cosine of the mean latitude (C{bool}).

962 @kwarg radius: Mean earth radius (C{meter}) or C{None}.

963 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

964 the B{C{points}} (C{bool}).

965

966 @return: Approximate area (I{square} C{meter}, same units as

967 B{C{radius}} or C{radians} I{squared} if B{C{radius}}

968 is C{None}).

969

970 @raise PointsError: Insufficient number of B{C{points}}

971

972 @raise TypeError: Some B{C{points}} are not C{LatLon}.

973

974 @raise ValueError: Invalid B{C{radius}}.

975

976 @note: This area approximation has limited accuracy and is

977 ill-suited for regions exceeding several hundred Km

978 or Miles or with near-polar latitudes.

979

980 @see: L{sphericalNvector.areaOf}, L{sphericalTrigonometry.areaOf},

981 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.

982 '''

983 if _MODS.booleans.isBoolean(points):

984 a = points._sum1(areaOf, adjust=adjust, radius=None, wrap=wrap)

985 else:

986 a, _ = _area2(points, adjust, wrap)

987 return fabs(a if radius is None else (Radius(radius)**2 * a))

988

989

990 def boundsOf(points, wrap=False, LatLon=None): # was=True

991 '''Determine the bottom-left SW and top-right NE corners of a

992 path or polygon.

993

994 @arg points: The path or polygon points (C{LatLon}[]).

995 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

996 the B{C{points}} (C{bool}).

997 @kwarg LatLon: Optional class to return the C{bounds}

998 corners (C{LatLon}) or C{None}.

999

1000 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)} as

1001 B{C{LatLon}}s if B{C{LatLon}} is C{None} a

1002 L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}.

1003

1004 @raise PointsError: Insufficient number of B{C{points}}

1005

1006 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1007

1008 @see: Function L{quadOf}.

1009

1010 @example:

1011

1012 >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)

1013 >>> boundsOf(b) # False

1014 >>> 45.0, 1.0, 46.0, 2.0

1015 '''

1016 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)

1017 w, s, _ = e, n, _ = Ps[0]

1018

1019 v = w

1020 for x, y, _ in Ps.iterate(closed=False): # [1:]

1021 if wrap:

1022 _, x = unroll180(v, x, wrap=True)

1023 v = x

1024

1025 if w > x:

1026 w = x

1027 elif e < x:

1028 e = x

1029

1030 if s > y:

1031 s = y

1032 elif n < y:

1033 n = y

1034

1035 return Bounds4Tuple(s, w, n, e) if LatLon is None else \

1036 Bounds2Tuple(LatLon(s, w), LatLon(n, e)) # PYCHOK inconsistent

1037

1038

1039 def _distanceTo(Error, **name_points): # .frechet, .hausdorff, .heights

1040 '''(INTERNAL) Chech callable C{distanceTo} methods.

1041 '''

1042 name, ps = name_points.popitem()

1043 for i, p in enumerate(ps):

1044 if not callable(_xattr(p, distanceTo=None)):

1045 n = _distanceTo.__name__[1:]

1046 t = _SPACE_(_no_, callable.__name__, n)

1047 raise Error(Fmt.SQUARE(name, i), p, txt=t)

1048 return ps

1049

1050

1051 def centroidOf(points, wrap=False, LatLon=None): # was=True

1052 '''Determine the centroid of a polygon.

1053

1054 @arg points: The polygon points (C{LatLon}[]).

1055 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1056 B{C{points}} (C{bool}).

1057 @kwarg LatLon: Optional class to return the centroid (C{LatLon})

1058 or C{None}.

1059

1060 @return: Centroid (B{C{LatLon}}) or a L{LatLon2Tuple}C{(lat, lon)}

1061 if C{B{LatLon} is None}.

1062

1063 @raise PointsError: Insufficient number of B{C{points}}

1064

1065 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1066

1067 @raise ValueError: The B{C{points}} enclose a pole or

1068 near-zero area.

1069

1070 @see: U{Centroid<https://WikiPedia.org/wiki/Centroid#Of_a_polygon>} and

1071 Paul Bourke's U{Calculating The Area And Centroid Of A Polygon

1072 <https://www.SEAS.UPenn.edu/~ese502/lab-content/extra_materials/

1073 Polygon%20Area%20and%20Centroid.pdf>}, 1988.

1074 '''

1075 A, X, Y = Fsum(), Fsum(), Fsum()

1076

1077 # setting radius=1 converts degrees to radians

1078 Ps = LatLon2PsxyIter(points, loop=1, radius=_1_0, wrap=wrap)

1079 x1, y1, ll = Ps[0]

1080 pts = [ll] # for _areaError

1081 for p in Ps.iterate(closed=True):

1082 x2, y2, ll = p

1083 if len(pts) < 4:

1084 pts.append(ll)

1085 if wrap and not Ps.looped:

1086 _, x2 = unrollPI(x1, x2, wrap=True)

1087 t = x1 * y2 - x2 * y1

1088 A += t

1089 X += t * (x1 + x2)

1090 Y += t * (y1 + y2)

1091 # XXX more elaborately:

1092 # t1, t2 = x1 * y2, -(x2 * y1)

1093 # A.fadd_(t1, t2)

1094 # X.fadd_(t1 * x1, t1 * x2, t2 * x1, t2 * x2)

1095 # Y.fadd_(t1 * y1, t1 * y2, t2 * y1, t2 * y2)

1096 x1, y1 = x2, y2

1097

1098 a = A.fmul(_6_0).fover(_2_0)

1099 if isnear0(a):

1100 raise _areaError(pts, near_=_near_)

1101 y, x = degrees90(Y.fover(a)), degrees180(X.fover(a))

1102 return LatLon2Tuple(y, x) if LatLon is None else LatLon(y, x)

1103

1104

1105 def fractional(points, fi, j=None, wrap=None, LatLon=None, Vector=None, **kwds):

1106 '''Return the point at a given I{fractional} index.

1107

1108 @arg points: The points (C{LatLon}[], L{Numpy2LatLon}[],

1109 L{Tuple2LatLon}[], C{Cartesian}[], C{Vector3d}[],

1110 L{Vector3Tuple}[]).

1111 @arg fi: The fractional index (L{FIx}, C{float} or C{int}).

1112 @kwarg j: Optionally, index of the other point (C{int}).

1113 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1114 B{{points}} (C{bool}) or C{None} for a backward

1115 compatible L{LatLon2Tuple} or B{C{LatLon}} with

1116 averaged lat- and longitudes. Use C{True} or

1117 C{False} to get the I{fractional} point computed

1118 by method C{B{points}[fi].intermediateTo}.

1119 @kwarg LatLon: Optional class to return the I{intermediate},

1120 I{fractional} point (C{LatLon}) or C{None}.

1121 @kwarg Vector: Optional class to return the I{intermediate},

1122 I{fractional} point (C{Cartesian}, C{Vector3d})

1123 or C{None}.

1124 @kwarg kwds: Optional, additional B{C{LatLon}} I{or} B{C{Vector}}

1125 keyword arguments, ignored if both C{B{LatLon}} and

1126 C{B{Vector}} are C{None}.

1127

1128 @return: A L{LatLon2Tuple}C{(lat, lon)} if B{C{wrap}}, B{C{LatLon}}

1129 and B{C{Vector}} all are C{None}, the defaults.

1130

1131 An instance of B{C{LatLon}} if not C{None} I{or} an instance

1132 of B{C{Vector}} if not C{None}.

1133

1134 Otherwise with B{C{wrap}} either C{True} or C{False} and

1135 B{C{LatLon}} and B{C{Vector}} both C{None}, an instance of

1136 B{C{points}}' (sub-)class C{intermediateTo} I{fractional}.

1137

1138 Summarized as follows:

1139

1140 >>> wrap | LatLon | Vector | returned type/value

1141 # -------+--------+--------+--------------+------

1142 # | | | LatLon2Tuple | favg

1143 # None | None | None | or** |

1144 # | | | Vector3Tuple | favg

1145 # None | LatLon | None | LatLon | favg

1146 # None | None | Vector | Vector | favg

1147 # -------+--------+--------+--------------+------

1148 # True | None | None | points' | .iTo

1149 # True | LatLon | None | LatLon | .iTo

1150 # True | None | Vector | Vector | .iTo

1151 # -------+--------+--------+--------------+------

1152 # False | None | None | points' | .iTo

1153 # False | LatLon | None | LatLon | .iTo

1154 # False | None | Vector | Vector | .iTo

1155 # _____

1156 # favg) averaged lat, lon or x, y, z values

1157 # .iTo) value from points[fi].intermediateTo

1158 # **) depends on base class of points[fi]

1159

1160 @raise IndexError: Fractional index B{C{fi}} invalid or B{C{points}}

1161 not subscriptable or not closed.

1162

1163 @raise TypeError: Invalid B{C{LatLon}}, B{C{Vector}} or B{C{kwds}}

1164 argument.

1165

1166 @see: Class L{FIx} and method L{FIx.fractional}.

1167 '''

1168 if LatLon and Vector: # PYCHOK no cover

1169 kwds = _xkwds(kwds, fi=fi, LatLon=LatLon, Vector=Vector)

1170 raise _TypeError(txt=fractional.__name__, **kwds)

1171 w = wrap if LatLon else False # intermediateTo

1172 try:

1173 if not isscalar(fi) or fi < 0:

1174 raise IndexError

1175 n = _xattr(fi, fin=0)

1176 p = _fractional(points, fi, j, fin=n, wrap=w) # see .units.FIx

1177 if LatLon:

1178 p = LatLon(p.lat, p.lon, **kwds)

1179 elif Vector:

1180 p = Vector(p.x, p.y, p.z, **kwds)

1181 except (IndexError, TypeError):

1182 raise _IndexError(fi=fi, points=points, wrap=w, txt=fractional.__name__)

1183 return p

1184

1185

1186 def _fractional(points, fi, j, fin=None, wrap=None): # in .frechet.py

1187 '''(INTERNAL) Compute point at L{fractional} index C{fi} and C{j}.

1188 '''

1189 i = int(fi)

1190 p = points[i]

1191 r = fi - float(i)

1192 if r > EPS: # EPS0?

1193 if j is None: # in .frechet.py

1194 j = i + 1

1195 if fin:

1196 j %= fin

1197 q = points[j]

1198 if r >= EPS1: # PYCHOK no cover

1199 p = q

1200 elif wrap is not None: # in (True, False)

1201 p = p.intermediateTo(q, r, wrap=wrap)

1202 elif _isLatLon(p): # backward compatible default

1203 p = LatLon2Tuple(favg(p.lat, q.lat, f=r),

1204 favg(p.lon, q.lon, f=r),

1205 name=fractional.__name__)

1206 else: # assume p and q are cartesian or vectorial

1207 z = p.z if p.z is q.z else favg(p.z, q.z, f=r)

1208 p = Vector3Tuple(favg(p.x, q.x, f=r),

1209 favg(p.y, q.y, f=r), z,

1210 name=fractional.__name__)

1211 return p

1212

1213

1214 def isclockwise(points, adjust=False, wrap=True):

1215 '''Determine the direction of a path or polygon.

1216

1217 @arg points: The path or polygon points (C{LatLon}[]).

1218 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta

1219 by the cosine of the mean latitude (C{bool}).

1220 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1221 B{C{points}} (C{bool}).

1222

1223 @return: C{True} if B{C{points}} are clockwise, C{False} otherwise.

1224

1225 @raise PointsError: Insufficient number of B{C{points}}

1226

1227 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1228

1229 @raise ValueError: The B{C{points}} enclose a pole or zero area.

1230

1231 @example:

1232

1233 >>> f = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)

1234 >>> isclockwise(f) # False

1235 >>> isclockwise(reversed(f)) # True

1236 '''

1237 a, pts = _area2(points, adjust, wrap)

1238 if a > 0: # opposite of ellipsoidalExact and -Karney

1239 return True

1240 elif a < 0:

1241 return False

1242 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>

1243 raise _areaError(pts)

1244

1245

1246 def isconvex(points, adjust=False, wrap=False): # was=True

1247 '''Determine whether a polygon is convex.

1248

1249 @arg points: The polygon points (C{LatLon}[]).

1250 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta

1251 by the cosine of the mean latitude (C{bool}).

1252 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1253 B{C{points}} (C{bool}).

1254

1255 @return: C{True} if B{C{points}} are convex, C{False} otherwise.

1256

1257 @raise CrossError: Some B{C{points}} are colinear.

1258

1259 @raise PointsError: Insufficient number of B{C{points}}

1260

1261 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1262

1263 @example:

1264

1265 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)

1266 >>> isconvex(t) # True

1267

1268 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)

1269 >>> isconvex(f) # False

1270 '''

1271 return bool(isconvex_(points, adjust=adjust, wrap=wrap))

1272

1273

1274 def isconvex_(points, adjust=False, wrap=False): # was=True

1275 '''Determine whether a polygon is convex I{and clockwise}.

1276

1277 @arg points: The polygon points (C{LatLon}[]).

1278 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta

1279 by the cosine of the mean latitude (C{bool}).

1280 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1281 B{C{points}} (C{bool}).

1282

1283 @return: C{+1} if B{C{points}} are convex clockwise, C{-1} for

1284 convex counter-clockwise B{C{points}}, C{0} otherwise.

1285

1286 @raise CrossError: Some B{C{points}} are colinear.

1287

1288 @raise PointsError: Insufficient number of B{C{points}}

1289

1290 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1291

1292 @example:

1293

1294 >>> t = LatLon(45,1), LatLon(46,1), LatLon(46,2)

1295 >>> isconvex_(t) # +1

1296

1297 >>> f = LatLon(45,1), LatLon(46,2), LatLon(45,2), LatLon(46,1)

1298 >>> isconvex_(f) # 0

1299 '''

1300 if adjust:

1301 def _unroll2(x1, x2, w, y1, y2):

1302 x21, x2 = unroll180(x1, x2, wrap=w)

1303 y = radians(y1 + y2) * _0_5

1304 x21 *= cos(y) if fabs(y) < PI_2 else _0_0

1305 return x21, x2

1306 else:

1307 def _unroll2(x1, x2, w, *unused): # PYCHOK expected

1308 return unroll180(x1, x2, wrap=w)

1309

1310 c, s = crosserrors(), 0

1311

1312 Ps = LatLon2PsxyIter(points, loop=2, wrap=wrap)

1313 x1, y1, _ = Ps[0]

1314 x2, y2, _ = Ps[1]

1315

1316 x21, x2 = _unroll2(x1, x2, False, y1, y2)

1317 for i, p in Ps.enumerate(closed=True):

1318 x3, y3, ll = p

1319 x32, x3 = _unroll2(x2, x3, bool(wrap and not Ps.looped), y2, y3)

1320

1321 # get the sign of the distance from point

1322 # x3, y3 to the line from x1, y1 to x2, y2

1323 # <https://WikiPedia.org/wiki/Distance_from_a_point_to_a_line>

1324 s3 = fdot((x3, y3, x1, y1), y2 - y1, -x21, -y2, x2)

1325 if s3 > 0: # x3, y3 on the right

1326 if s < 0: # non-convex

1327 return 0

1328 s = +1

1329

1330 elif s3 < 0: # x3, y3 on the left

1331 if s > 0: # non-convex

1332 return 0

1333 s = -1

1334

1335 elif c and fdot((x32, y1 - y2), y3 - y2, -x21) < 0: # PYCHOK no cover

1336 # colinear u-turn: x3, y3 not on the

1337 # opposite side of x2, y2 as x1, y1

1338 t = Fmt.SQUARE(points=i)

1339 raise CrossError(t, ll, txt=_colinear_)

1340

1341 x1, y1, x2, y2, x21 = x2, y2, x3, y3, x32

1342

1343 return s # all points on the same side

1344

1345

1346 def isenclosedBy(point, points, wrap=False): # MCCABE 15

1347 '''Determine whether a point is enclosed by a polygon or composite.

1348

1349 @arg point: The point (C{LatLon} or 2-tuple C{(lat, lon)}).

1350 @arg points: The polygon points or clips (C{LatLon}[], L{BooleanFHP}

1351 or L{BooleanGH}).

1352 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1353 B{C{points}} (C{bool}).

1354

1355 @return: C{True} if the B{C{point}} is inside the polygon or

1356 composite, C{False} otherwise.

1357

1358 @raise PointsError: Insufficient number of B{C{points}}

1359

1360 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1361

1362 @raise ValueError: Invalid B{C{point}}, lat- or longitude.

1363

1364 @see: Functions L{pygeodesy.isconvex} and L{pygeodesy.ispolar} especially

1365 if the B{C{points}} may enclose a pole or wrap around the earth

1366 I{longitudinally}, methods L{sphericalNvector.LatLon.isenclosedBy},

1367 L{sphericalTrigonometry.LatLon.isenclosedBy} and U{MultiDop

1368 GeogContainPt<https://GitHub.com/NASA/MultiDop>} (U{Shapiro et.al. 2009,

1369 JTECH<https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}

1370 and U{Potvin et al. 2012, JTECH <https://Journals.AMetSoc.org/doi/abs/

1371 10.1175/JTECH-D-11-00019.1>}).

1372 '''

1373 try:

1374 y0, x0 = point.lat, point.lon

1375 except AttributeError:

1376 try:

1377 y0, x0 = map(float, point[:2])

1378 except (IndexError, TypeError, ValueError) as x:

1379 raise _ValueError(point=point, cause=x)

1380

1381 if wrap:

1382 y0, x0 = _Wrap.latlon(y0, x0)

1383

1384 def _dxy3(x, x2, y2, Ps):

1385 dx, x2 = unroll180(x, x2, wrap=not Ps.looped)

1386 return dx, x2, y2

1387

1388 else:

1389 x0 = fmod(x0, _360_0) # not x0 % 360!

1390 x0_180_ = x0 - _180_0

1391 x0_180 = x0 + _180_0

1392

1393 def _dxy3(x1, x, y, unused): # PYCHOK expected

1394 x = _umod_360(float(x))

1395 if x < x0_180_:

1396 x += _360_0

1397 elif x >= x0_180:

1398 x -= _360_0

1399 return (x - x1), x, y

1400

1401 if _MODS.booleans.isBoolean(points):

1402 return points._encloses(y0, x0, wrap=wrap)

1403

1404 Ps = LatLon2PsxyIter(points, loop=1, wrap=wrap)

1405 p = Ps[0]

1406 e = m = False

1407 S = Fsum()

1408

1409 _, x1, y1 = _dxy3(x0, p.x, p.y, False)

1410 for p in Ps.iterate(closed=True):

1411 dx, x2, y2 = _dxy3(x1, p.x, p.y, Ps)

1412 # ignore duplicate and near-duplicate pts

1413 if fabs(dx) > EPS or fabs(y2 - y1) > EPS:

1414 # determine if polygon edge (x1, y1)..(x2, y2) straddles

1415 # point (lat, lon) or is on boundary, but do not count

1416 # edges on boundary as more than one crossing

1417 if fabs(dx) < 180 and (x1 < x0 <= x2 or x2 < x0 <= x1):

1418 m = not m

1419 dy = (x0 - x1) * (y2 - y1) - (y0 - y1) * dx

1420 if (dy > 0 and dx >= 0) or (dy < 0 and dx <= 0):

1421 e = not e

1422

1423 S += sin(radians(y2))

1424 x1, y1 = x2, y2

1425

1426 # An odd number of meridian crossings means, the polygon

1427 # contains a pole. Assume it is the pole on the hemisphere

1428 # containing the polygon mean point and if the polygon does

1429 # contain the North Pole, flip the result.

1430 if m and S.fsum() > 0:

1431 e = not e

1432 return e

1433

1434

1435 def ispolar(points, wrap=False):

1436 '''Check whether a polygon encloses a pole.

1437

1438 @arg points: The polygon points (C{LatLon}[]).

1439 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

1440 the B{C{points}} (C{bool}).

1441

1442 @return: C{True} if the polygon encloses a pole, C{False}

1443 otherwise.

1444

1445 @raise PointsError: Insufficient number of B{C{points}}

1446

1447 @raise TypeError: Some B{C{points}} are not C{LatLon} or don't

1448 have C{bearingTo2}, C{initialBearingTo}

1449 and C{finalBearingTo} methods.

1450 '''

1451 def _cds(ps, w): # iterate over course deltas

1452 Ps = PointsIter(ps, loop=2, wrap=w)

1453 p2, p1 = Ps[0:2]

1454 b1, _ = _bearingTo2(p2, p1, wrap=False)

1455 for p2 in Ps.iterate(closed=True):

1456 if not p2.isequalTo(p1, EPS):

1457 if w and not Ps.looped:

1458 p2 = _unrollon(p1, p2)

1459 b, b2 = _bearingTo2(p1, p2, wrap=False)

1460 yield wrap180(b - b1) # (b - b1 + 540) % 360 - 180

1461 yield wrap180(b2 - b) # (b2 - b + 540) % 360 - 180

1462 p1, b1 = p2, b2

1463

1464 # summation of course deltas around pole is 0° rather than normally ±360°

1465 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>

1466 s = fsum(_cds(points, wrap), floats=True)

1467 # XXX fix (intermittant) edge crossing pole - eg (85,90), (85,0), (85,-90)

1468 return fabs(s) < 90 # "zero-ish"

1469

1470

1471 def luneOf(lon1, lon2, closed=False, LatLon=LatLon_, **LatLon_kwds):

1472 '''Generate an ellipsoidal or spherical U{lune

1473 <https://WikiPedia.org/wiki/Spherical_lune>}-shaped path or polygon.

1474

1475 @arg lon1: Left longitude (C{degrees90}).

1476 @arg lon2: Right longitude (C{degrees90}).

1477 @kwarg closed: Optionally, close the path (C{bool}).

1478 @kwarg LatLon: Class to use (L{LatLon_}).

1479 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}

1480 keyword arguments.

1481

1482 @return: A tuple of 4 or 5 B{C{LatLon}} instances outlining

1483 the lune shape.

1484

1485 @see: U{Latitude-longitude quadrangle

1486 <https://www.MathWorks.com/help/map/ref/areaquad.html>}.

1487 '''

1488 t = (LatLon( _0_0, lon1, **LatLon_kwds),

1489 LatLon( _90_0, lon1, **LatLon_kwds),

1490 LatLon( _0_0, lon2, **LatLon_kwds),

1491 LatLon(_N_90_0, lon2, **LatLon_kwds))

1492 if closed:

1493 t += t[:1]

1494 return t

1495

1496

1497 def nearestOn5(point, points, closed=False, wrap=False, adjust=True,

1498 limit=9, **LatLon_and_kwds):

1499 '''Locate the point on a path or polygon closest to a reference point.

1500

1501 The closest point on each polygon edge is either the nearest of that

1502 edge's end points or a point in between.

1503

1504 @arg point: The reference point (C{LatLon}).

1505 @arg points: The path or polygon points (C{LatLon}[]).

1506 @kwarg closed: Optionally, close the path or polygon (C{bool}).

1507 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1508 B{C{points}} (C{bool}).

1509 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).

1510 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),

1511 default C{9 degrees} is about C{1,000 Kmeter} (for mean

1512 spherical earth radius L{R_KM}).

1513 @kwarg LatLon_and_kwds: Optional, C{B{LatLon}=None} class to use for

1514 the closest point and additional B{C{LatLon}} keyword

1515 arguments, ignored if C{B{LatLon}=None} or not given.

1516

1517 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the

1518 {closest} point (B{C{LatLon}}) or if C{B{LatLon} is None},

1519 a L{NearestOn5Tuple}C{(lat, lon, distance, angle, height)}.

1520 The C{distance} is the L{pygeodesy.equirectangular} distance

1521 between the C{closest} and reference B{C{point}} in C{degrees}.

1522 The C{angle} from the B{C{point}} to the C{closest} is in

1523 compass C{degrees}, like function L{pygeodesy.compassAngle}.

1524

1525 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},

1526 see function L{pygeodesy.equirectangular_}.

1527

1528 @raise PointsError: Insufficient number of B{C{points}}

1529

1530 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1531

1532 @note: Distances are I{approximated} by function L{pygeodesy.equirectangular_}.

1533 For more accuracy use one of the C{LatLon.nearestOn6} methods.

1534

1535 @see: Function L{pygeodesy.degrees2m}.

1536 '''

1537 def _d2yx4(p2, p1, u, alw):

1538 # w = wrap if (i < (n - 1) or not closed) else False

1539 # equirectangular_ returns a Distance4Tuple(distance

1540 # in degrees squared, delta lat, delta lon, p2.lon

1541 # unroll/wrap'd); the previous p2.lon unroll/wrap'd

1542 # is also applied to the next edge's p1.lon

1543 return equirectangular_(p1.lat, p1.lon + u,

1544 p2.lat, p2.lon, **alw)

1545

1546 def _h(p): # get height or default 0

1547 return _xattr(p, height=0) or 0

1548

1549 # 3-D version used in .vector3d._nearestOn2

1550 #

1551 # point (x, y) on axis rotated ccw by angle a:

1552 # x' = x * cos(a) + y * sin(a)

1553 # y' = y * cos(a) - x * sin(a)

1554 #

1555 # distance (w) along and (h) perpendicular to

1556 # a line thru point (dx, dy) and the origin:

1557 # d = hypot(dx, dy)

1558 # w = (x * dx + y * dy) / d

1559 # h = (y * dx - x * dy) / d

1560 #

1561 # closest point on that line thru (dx, dy):

1562 # xc = dx * w / d

1563 # yc = dy * w / d

1564 # or

1565 # xc = dx * f

1566 # yc = dy * f

1567 # with

1568 # f = w / d

1569 # or

1570 # f = (y * dy + x * dx) / hypot2(dx, dy)

1571 #

1572 # i.e. no need for sqrt or hypot

1573

1574 Ps = PointsIter(points, loop=1, wrap=wrap)

1575 p1 = c = Ps[0]

1576 u1 = u = _0_0

1577 kw = dict(adjust=adjust, limit=limit, wrap=False)

1578 d, dy, dx, _ = _d2yx4(p1, point, u1, kw)

1579 for p2 in Ps.iterate(closed=closed):

1580 # iff wrapped, unroll lon1 (actually previous

1581 # lon2) like function unroll180/-PI would've

1582 if wrap:

1583 kw.update(wrap=not (closed and Ps.looped))

1584 d21, y21, x21, u2 = _d2yx4(p2, p1, u1, kw)

1585 if d21 > EPS:

1586 # distance point to p1, y01 and x01 negated

1587 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)

1588 if d2 > EPS:

1589 w2 = y01 * y21 + x01 * x21

1590 if w2 > 0:

1591 if w2 < d21:

1592 # closest is between p1 and p2, use

1593 # original delta's, not y21 and x21

1594 f = w2 / d21

1595 p1 = LatLon_(favg(p1.lat, p2.lat, f=f),

1596 favg(p1.lon, p2.lon + u2, f=f),

1597 height=favg(_h(p1), _h(p2), f=f))

1598 u1 = _0_0

1599 else: # p2 is closest

1600 p1, u1 = p2, u2

1601 d2, y01, x01, _ = _d2yx4(point, p1, u1, kw)

1602 if d2 < d: # p1 is closer, y01 and x01 negated

1603 c, u, d, dy, dx = p1, u1, d2, -y01, -x01

1604 p1, u1 = p2, u2

1605

1606 a = atan2b(dx, dy) # azimuth

1607 d = hypot( dx, dy)

1608 h = _h(c)

1609 n = nameof(point) or nearestOn5.__name__

1610 if LatLon_and_kwds:

1611 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)

1612 LL = _xkwds_pop(kwds, LatLon=None)

1613 if LL is not None:

1614 r = LL(c.lat, c.lon + u, **kwds)

1615 return NearestOn3Tuple(r, d, a, name=n)

1616 return NearestOn5Tuple(c.lat, c.lon + u, d, a, h, name=n) # PYCHOK expected

1617

1618

1619 def perimeterOf(points, closed=False, adjust=True, radius=R_M, wrap=True):

1620 '''I{Approximate} the perimeter of a path, polygon. or composite.

1621

1622 @arg points: The path or polygon points or clips (C{LatLon}[],

1623 L{BooleanFHP} or L{BooleanGH}).

1624 @kwarg closed: Optionally, close the path or polygon (C{bool}).

1625 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta

1626 by the cosine of the mean latitude (C{bool}).

1627 @kwarg radius: Mean earth radius (C{meter}).

1628 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the

1629 B{C{points}} (C{bool}).

1630

1631 @return: Approximate perimeter (C{meter}, same units as

1632 B{C{radius}}).

1633

1634 @raise PointsError: Insufficient number of B{C{points}}

1635

1636 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1637

1638 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with

1639 C{B{points}} a composite.

1640

1641 @note: This perimeter is based on the L{pygeodesy.equirectangular_}

1642 distance approximation and is ill-suited for regions exceeding

1643 several hundred Km or Miles or with near-polar latitudes.

1644

1645 @see: Functions L{sphericalTrigonometry.perimeterOf} and

1646 L{ellipsoidalKarney.perimeterOf}.

1647 '''

1648 def _degs(ps, c, a, w): # angular edge lengths in degrees

1649 Ps = LatLon2PsxyIter(ps, loop=1) # wrap=w

1650 p1, u = Ps[0], _0_0 # previous x2's unroll/wrap

1651 for p2 in Ps.iterate(closed=c):

1652 if w and c:

1653 w = not Ps.looped

1654 # apply previous x2's unroll/wrap'd to new x1

1655 _, dy, dx, u = equirectangular_(p1.y, p1.x + u,

1656 p2.y, p2.x,

1657 adjust=a, limit=None,

1658 wrap=w) # PYCHOK non-seq

1659 yield hypot(dx, dy)

1660 p1 = p2

1661

1662 if _MODS.booleans.isBoolean(points):

1663 if not closed:

1664 notImplemented(None, closed=closed, points=_composite_)

1665 d = points._sum1(perimeterOf, closed=True, adjust=adjust,

1666 radius=radius, wrap=wrap)

1667 else:

1668 d = fsum(_degs(points, closed, adjust, wrap), floats=True)

1669 return degrees2m(d, radius=radius)

1670

1671

1672 def quadOf(latS, lonW, latN, lonE, closed=False, LatLon=LatLon_, **LatLon_kwds):

1673 '''Generate a quadrilateral path or polygon from two points.

1674

1675 @arg latS: Souther-nmost latitude (C{degrees90}).

1676 @arg lonW: Western-most longitude (C{degrees180}).

1677 @arg latN: Norther-nmost latitude (C{degrees90}).

1678 @arg lonE: Eastern-most longitude (C{degrees180}).

1679 @kwarg closed: Optionally, close the path (C{bool}).

1680 @kwarg LatLon: Class to use (L{LatLon_}).

1681 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}}

1682 keyword arguments.

1683

1684 @return: Return a tuple of 4 or 5 B{C{LatLon}} instances

1685 outlining the quadrilateral.

1686

1687 @see: Function L{boundsOf}.

1688 '''

1689 t = (LatLon(latS, lonW, **LatLon_kwds),

1690 LatLon(latN, lonW, **LatLon_kwds),

1691 LatLon(latN, lonE, **LatLon_kwds),

1692 LatLon(latS, lonE, **LatLon_kwds))

1693 if closed:

1694 t += t[:1]

1695 return t

1696

1697

1698__all__ += _ALL_DOCS(_Array2LatLon, _Basequence)

1699

1700# **) MIT License

1701#

1702# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

1703#

1704# Permission is hereby granted, free of charge, to any person obtaining a

1705# copy of this software and associated documentation files (the "Software"),

1706# to deal in the Software without restriction, including without limitation

1707# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1708# and/or sell copies of the Software, and to permit persons to whom the

1709# Software is furnished to do so, subject to the following conditions:

1710#

1711# The above copyright notice and this permission notice shall be included

1712# in all copies or substantial portions of the Software.

1713#

1714# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1715# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1716# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1717# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1718# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1719# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1720# OTHER DEALINGS IN THE SOFTWARE.