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, latDMS, lonDMS, parseDMS, S_DEG, S_DMS, \

36 S_MIN, S_SEC, S_SEP

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

38 _IsnotError, _TypeError, _ValueError, \

39 _xattr, _xkwds, _xkwds_pop

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

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

42 from pygeodesy.formy import _bearingTo2, equirectangular_, _isequalTo, \

43 isnormal, latlon2n_xyz, normal, _spherical_datum

44 from 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_

48 from pygeodesy.iters import LatLon2PsxyIter, PointsIter, points2

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

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

51 _NamedTuple, _xnamed, _xother3, _xotherError

52 from pygeodesy.namedTuples import Bounds2Tuple, Bounds4Tuple, \

53 LatLon2Tuple, NearestOn3Tuple, \

54 NearestOn5Tuple, PhiLam2Tuple, \

55 Point3Tuple, Vector3Tuple, Vector4Tuple

56 from pygeodesy.nvectorBase import NvectorBase, _N_vector_

57 from pygeodesy.props import deprecated_method, Property, Property_RO, \

58 property_doc_, property_RO, _update_all

59 from pygeodesy.streprs import Fmt, hstr, instr, pairs

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

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

62 unroll180, _unrollon, unrollPI, _Wrap, wrap180

63

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

65

66__all__ = _ALL_LAZY.points

67__version__ = '23.05.26'

68

69_ilat_ = 'ilat'

70_ilon_ = 'ilon'

71_ncols_ = 'ncols'

72_nrows_ = 'nrows'

73

74

75 class 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 "Problems

82 # with __slots__" in Luciano Ramalho, "Fluent Python", page

83 # 276+, O'Reilly, 2016, also at <https://Books.Google.ie/

84 # books?id=bIZHCgAAQBAJ&lpg=PP1&dq=fluent%20python&pg=

85 # PT364#v=onepage&q=“Problems%20with%20__slots__”&f=false>

86 #

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

88 # Property_RO = property_RO # no __dict__ with __slots__!

89 #

90 # However, sys.getsizeof(LatLon_(1, 2)) is 72-88 with __slots__

91 # but only 48-56 bytes without in Python 2.7.18+ and Python 3+.

92

93 def __init__(self, latlonh, lon=None, name=NN, height=0,

94 datum=None, wrap=False):

95 '''Creat a new, mininal, low-overhead L{LatLon_} instance.

96

97 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or

98 a previous C{LatLon} instance provided C{B{lon}=None}.

99 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or

100 C(None), indicating B{C{latlonh}} is a C{LatLon}.

101 @kwarg name: Optional name (C{str}).

102 @kwarg height: Optional height (C{meter}, conventionally).

103 @kwarg datum: Optional datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},

104 L{a_f2Tuple} or I{scalar} radius) or C{None}.

105 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}

106 (C{bool}).

107

108 @raise TypeError: Invalid B{C{datum}} or B{C{latlonh}} not a C{LatLon}.

109

110 @note: The lat- and longitude are taken as-given,

111 un-clipped and un-validated .

112 '''

113 if lon is None:

114 try:

115 ll = latlonh.lat, latlonh.lon

116 height = _xattr(latlonh, height=height)

117 except AttributeError:

118 raise _IsnotError(_LatLon_, latlonh=latlonh)

119 if wrap:

120 ll = _Wrap.latlon(*ll)

121 elif wrap:

122 ll = _Wrap.latlonDMS2(latlonh, lon)

123 else: # must be latNS, lonEW

124 ll = parseDMS(latlonh, suffix=_NS_), parseDMS(lon, suffix=_EW_)

125

126 self.lat, self.lon = ll

127 self.name = str(name) if name else NN

128 self.height = height

129 self.datum = datum if datum is None else \

130 _spherical_datum(datum, name=self.name)

131

132 def __eq__(self, other):

133 return isinstance(other, LatLon_) and \

134 other.lat == self.lat and \

135 other.lon == self.lon

136

137 def __ne__(self, other):

138 return not self.__eq__(other)

139

140 def __repr__(self):

141 return self.toRepr()

142

143 def __str__(self):

144 return self.toStr()

145

146 def classof(self, *args, **kwds):

147 '''Instantiate this very class.

148

149 @arg args: Optional, positional arguments.

150 @kwarg kwds: Optional, keyword arguments.

151

152 @return: New instance (C{self.__class__}).

153 '''

154 return _xnamed(self.__class__(*args, **kwds), self.name)

155

156 def copy(self, deep=False):

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

158

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

160 shallow copy (C{bool}).

161

162 @return: The copy (C{This} (sub-)class).

163 '''

164 return _xcopy(self, deep=deep)

165

166 def dup(self, **items):

167 '''Duplicate this instance, replacing some items.

168

169 @kwarg items: Attributes to be changed (C{any}).

170

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

172

173 @raise AttributeError: Some B{C{items}} invalid.

174 '''

175 d = _xdup(self, **items)

176 if items:

177 _update_all(d)

178 return d

179

180 def heightStr(self, prec=-2):

181 '''Return a string for the height B{C{height}}.

182

183 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

184

185 @see: Function L{pygeodesy.hstr}.

186 '''

187 return hstr(self.height, prec=prec)

188

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

190 '''Locate the point at a given fraction between (or along) this

191 and an other point.

192

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

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

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

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

197 intermediate height.

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

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

200

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

202

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

204 '''

205 f = Scalar(fraction=fraction)

206 if isnear0(f):

207 r = self

208 elif isnear1(f) and not wrap:

209 r = self.others(other)

210 else:

211 p = self.others(other)

212 h = favg(self.height, p.height, f=f) if height is None else height

213 _, lat, lon = _Wrap.latlon3(self.lon, p.lat, p.lon, wrap=wrap)

214 r = self.classof(favg(self.lat, lat, f=f),

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

216 height=h, datum=self.datum,

217 name=self.intermediateTo.__name__)

218 return r

219

220 @Property_RO # PYCHOK no cover

221 def isEllipsoidal(self):

222 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).

223 '''

224 return self.datum.isEllipsoidal if self.datum else None

225

226 @Property_RO # PYCHOK no cover

227 def isEllipsoidalLatLon(self):

228 '''Get C{LatLon} base.

229 '''

230 return False

231

232 def isequalTo(self, other, eps=None):

233 '''Compare this point with an other point, I{ignoring} height.

234

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

236 @kwarg eps: Tolerance for equality (C{degrees}).

237

238 @return: C{True} if both points are identical,

239 I{ignoring} height, C{False} otherwise.

240

241 @raise UnitError: Invalid B{C{eps}}.

242 '''

243 return _isequalTo(self, self.others(other), eps=eps)

244

245 @Property_RO

246 def isnormal(self):

247 '''Return C{True} if this point is normal (C{bool}),

248 meaning C{abs(lat) <= 90} and C{abs(lon) <= 180}.

249

250 @see: Methods L{normal}, L{toNormal} and functions

251 L{pygeodesy.isnormal} and L{pygeodesy.normal}.

252 '''

253 return isnormal(self.lat, self.lon, eps=0)

254

255 @Property_RO

256 def isSpherical(self): # PYCHOK no cover

257 '''Check whether this point is spherical (C{bool} or C{None} if unknown).

258 '''

259 return self.datum.isSpherical if self.datum else None

260

261 @Property_RO

262 def lam(self):

263 '''Get the longitude (B{C{radians}}).

264 '''

265 return radians(self.lon)

266

267 @Property

268 def latlon(self):

269 '''Get the lat- and longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).

270 '''

271 return LatLon2Tuple(self.lat, self.lon, name=self.name)

272

273 @latlon.setter # PYCHOK setter!

274 def latlon(self, latlon):

275 '''Set the lat- and longitude.

276

277 @arg latlon: New lat- and longitude in C{degrees} (C{2-tuple} or {-list}).

278 '''

279 if self.latlon != latlon[:2]:

280 _update_all(self)

281 self.lat, self.lon = latlon[:2]

282

283 @Property_RO

284 def latlonheight(self):

285 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).

286 '''

287 return self.latlon.to3Tuple(self.height)

288

289 @Property_RO

290 def _N_vector(self):

291 '''(INTERNAL) Get the minimal, low-overhead (C{nvectorBase._N_vector_})

292 '''

293 return _N_vector_(*latlon2n_xyz(self.lat, self.lon),

294 h=self.height, name=self.name)

295

296 def others(self, *other, **name_other_up): # see .named._namedBase.others

297 '''Refined class comparison.

298

299 @arg other: The other instance (any C{type}).

300 @kwarg name_other_up: Overriding C{name=other} and C{up=1}

301 keyword arguments.

302

303 @return: The B{C{other}} if compatible.

304

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

306 '''

307 other, name, up = _xother3(self, other, **name_other_up)

308 if isinstance(other, self.__class__) or (hasattr(other, _lat_)

309 and hasattr(other, _lon_)):

310 return other

311 raise _xotherError(self, other, name=name, up=up + 1)

312

313 def normal(self):

314 '''Normalize this point I{in-place} to C{abs(lat) <= 90} and

315 C{abs(lon) <= 180}.

316

317 @return: C{True} if this point was I{normal}, C{False} if it

318 wasn't (but is now).

319

320 @see: Property L{isnormal} and method L{toNormal}.

321 '''

322 n = self.isnormal

323 if not n:

324 self.latlon = normal(*self.latlon)

325 return n

326

327 @Property_RO

328 def phi(self):

329 '''Get the latitude (B{C{radians}}).

330 '''

331 return radians(self.lat)

332

333 @Property_RO

334 def philam(self):

335 '''Get the lat- and longitude (L{PhiLam2Tuple}C{(phi, lam)}).

336 '''

337 return PhiLam2Tuple(self.phi, self.lam, name=self.name)

338

339 @Property_RO

340 def philamheight(self):

341 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).

342 '''

343 return self.philam.to3Tuple(self.height)

344

345 @deprecated_method

346 def points(self, points, closed=False, base=None): # PYCHOK no cover

347 '''DEPRECATED, use method C{points2}.'''

348 return points2(points, closed=closed, base=base)

349

350 def points2(self, points, closed=False, base=None):

351 '''Check a path or polygon represented by points.

352

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

354 @kwarg closed: Optionally, consider the polygon closed,

355 ignoring any duplicate or closing final

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

357 @kwarg base: Optionally, check all B{C{points}} against

358 this base class, if C{None} don't check.

359

360 @return: A L{Points2Tuple}C{(number, points)} with the number

361 of points and the points C{list} or C{tuple}.

362

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

364

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

366 '''

367 return points2(points, closed=closed, base=base)

368

369 def PointsIter(self, points, loop=0, dedup=False):

370 '''Return a points iterator.

371

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

373 @kwarg loop: Number of loop-back points (non-negative C{int}).

374 @kwarg dedup: Skip duplicate points (C{bool}).

375

376 @return: A new C{PointsIter} iterator.

377

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

379 '''

380 return PointsIter(points, base=self, loop=loop, dedup=dedup)

381

382 @deprecated_method

383 def to2ab(self): # PYCHOK no cover

384 '''DEPRECATED, use property L{philam}.'''

385 return self.philam

386

387 def toNormal(self, deep=False, name=NN):

388 '''Get a copy of this point normalized to C{abs(lat) <= 90}

389 and C{abs(lon) <= 180}.

390

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

392 shallow copy (C{bool}).

393 @kwarg name: Optional name of the copy (C{str}).

394

395 @return: A copy of this point, I{normalized} and

396 optionally renamed (C{LatLon}).

397

398 @see: Property L{isnormal}, method L{normal} and

399 function L{pygeodesy.normal}.

400 '''

401 ll = self.copy(deep=deep)

402 _ = ll.normal()

403 if name:

404 ll.name = str(name)

405 return ll

406

407 def toNvector(self, h=None, Nvector=NvectorBase, **Nvector_kwds):

408 '''Convert this point to C{n-vector} (normal to the earth's

409 surface) components, I{including height}.

410

411 @kwarg h: Optional height, overriding this point's height

412 (C{meter}).

413 @kwarg Nvector: Optional class to return the C{n-vector}

414 components (C{Nvector}) or C{None}.

415 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} keyword

416 arguments, ignored if C{B{Nvector} is None}.

417

418 @return: The C{n-vector} components B{C{Nvector}} or if

419 B{C{Nvector}} is C{None}, a L{Vector4Tuple}C{(x,

420 y, z, h)}.

421

422 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}

423 argument.

424 '''

425 x, y, z = latlon2n_xyz(self.lat, self.lon)

426 r = Vector4Tuple(x, y, z, self.height if h is None else h)

427 if Nvector is not None:

428 r = Nvector(x, y, z, **_xkwds(Nvector_kwds, h=r.h, ll=self))

429 return _xnamed(r, self.name)

430

431 def toRepr(self, **kwds):

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

433

434 @kwarg kwds: Optional, keyword arguments.

435

436 @return: Class instance (C{str}).

437 '''

438 _ = _xkwds_pop(kwds, std=None) # PYCHOK std unused

439 return Fmt.PAREN(classname(self), self.toStr(**kwds))

440

441 def toStr(self, form=F_D, joined=_COMMASPACE_, **prec_sep_s_D_M_S_kwds):

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

443 formatted in the given C{B{form}at}.

444

445 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see

446 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.

447 @kwarg joined: Separator to join the lat-, longitude, heigth, name

448 and other strings (C{str} or C{None} or C{NN} for

449 non-joined).

450 @kwarg prec_sep_s_D_M_S_kwds: Optional C{B{prec}ision}, C{B{sep}arator},

451 B{C{s_D}}, B{C{s_M}}, B{C{s_S}}, B{C{s_DMS}} and possibly

452 other keyword arguments, see functions L{pygeodesy.latDMS}

453 or L{pygeodesy.lonDMS}.

454

455 @return: This point in the specified C{B{form}at}, etc. (C{str} or

456 a 2- or 3+tuple C{(lat_str, lon_str[, height_str][, name_str][,

457 ...])} if C{B{joined}=NN} or C{B{joined}=None} and with the

458 C{height_str} and C{name_str} only included if non-zero

459 respectively non-empty).

460

461 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more

462 details about keyword arguments C{B{form}at}, C{B{prec}ision},

463 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.

464 '''

465 def _t3(prec=None, sep=S_SEP, s_D=S_DEG, s_M=S_MIN, s_S=S_SEC, s_DMS=S_DMS, **kwds):

466 return dict(prec=prec, sep=sep, s_D=s_D, s_M=s_M, s_S=s_S, s_DMS=s_DMS), kwds, prec

467

468 prec_sep_s_D_M_S, kwds, prec = _t3(**prec_sep_s_D_M_S_kwds)

469 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),

470 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))

471 if self.height:

472 t += (self.heightStr(),)

473 if self.name:

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

475 if kwds:

476 t += pairs(kwds) if prec is None else pairs(kwds, prec=prec)

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

478

479 @deprecated_method

480 def toStr2(self, **kwds): # PYCHOK no cover

481 '''DEPRECATED, used method L{toRepr}.'''

482 return self.toRepr(**kwds)

483

484

485 def _isLatLon(inst):

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

487 '''

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

489

490

491 def _isLatLon_(LL):

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

493 '''

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

495 all(hasattr(LL, a) for a in _ALL_ATTRS_))

496

497

498# get all pseudo-slots for class C{LatLon_}

499_ALL_ATTRS_ = tuple(LatLon_(0, 0).__dict__.keys())

500

501

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

503 '''(INTERNAL) Base class.

504 '''

505 _array = []

506 _epsilon = EPS

507 _itemname = _point_

508

509 def _contains(self, point):

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

511 '''

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

513

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

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

516

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

518 shallow copy (C{bool}).

519

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

521 '''

522 return _xcopy(self, deep=deep)

523

524 def _count(self, point):

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

526 '''

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

528

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

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

531

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

533

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

535

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

537 '''

538 if items:

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

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

541 return _xdup(self)

542

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

544 def epsilon(self):

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

546 '''

547 return self._epsilon

548

549 @epsilon.setter # PYCHOK setter!

550 def epsilon(self, tol):

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

552

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

554 '''

555 self._epsilon = Scalar_(tolerance=tol)

556

557 def _find(self, point, start_end):

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

559 '''

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

561 return i

562 return -1

563

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

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

566 '''

567 notImplemented(self, point, start_end)

568

569 def _getitem(self, index):

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

571 '''

572 # Luciano Ramalho, "Fluent Python", page 290+, O'Reilly, 2016

573 if isinstance(index, slice):

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

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

576 else:

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

578

579 def _index(self, point, start_end):

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

581 '''

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

583 return i

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

585

586 @property_RO

587 def isNumpy2(self): # PYCHOK no cover

588 '''Is this a Numpy2 wrapper?

589 '''

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

591

592 @property_RO

593 def isPoints2(self): # PYCHOK no cover

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

595 '''

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

597

598 @property_RO

599 def isTuple2(self): # PYCHOK no cover

600 '''Is this a Tuple2 wrapper?

601 '''

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

603

604 def _iter(self):

605 '''(INTERNAL) Yield all points.

606 '''

607 _array, _point = self._array, self.point

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

609 yield _point(_array[i])

610

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

612 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded in}.

613

614 @arg attrs: Optional arguments.

615 '''

616 notOverloaded(self, *attrs)

617

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

619 '''(INTERNAL) Return the range.

620 '''

621 if step > 0:

622 if start is None:

623 start = 0

624 if end is None:

625 end = len(self)

626 elif step < 0:

627 if start is None:

628 start = len(self) - 1

629 if end is None:

630 end = -1

631 else:

632 raise _ValueError(step=step)

633 return range(start, end, step)

634

635 def _repr(self):

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

637 '''

638 # XXX use Python 3+ reprlib.repr

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

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

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

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

643

644 def _reversed(self): # PYCHOK false

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

646 '''

647 _array, point = self._array, self.point

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

649 yield point(_array[i])

650

651 def _rfind(self, point, start_end):

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

653 '''

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

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

656

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

658 return i

659 return -1

660

661 def _slicekwds(self): # PYCHOK no cover

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

663 '''

664 return {}

665

666

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

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

669 '''

670 _array = ()

671 _ilat = 0 # row column index

672 _ilon = 0 # row column index

673 _LatLon = LatLon_ # default

674 _shape = ()

675

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

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

678 '''

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

680

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

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

683

684 self._array = array

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

686

687 if LatLon: # check the point class

688 if not _isLatLon_(LatLon):

689 raise _IsnotError(_valid_, LatLon=LatLon)

690 self._LatLon = LatLon

691

692 # check the attr indices

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

694 if not isint(i):

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

696 i = int(i)

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

698 raise _ValueError(ai, i)

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

700 if int(j) == i:

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

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

703

704 def __contains__(self, latlon):

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

706

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

708 C{(lat, lon)}).

709

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

711

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

713 '''

714 return self._contains(latlon)

715

716 def __getitem__(self, index):

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

718 '''

719 return self._getitem(index)

720

721 def __iter__(self):

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

723 '''

724 return self._iter()

725

726 def __len__(self):

727 '''Return the number of rows.

728 '''

729 return self._shape[0]

730

731 def __repr__(self):

732 '''Return a string representation.

733 '''

734 return self._repr()

735

736 def __reversed__(self): # PYCHOK false

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

738 '''

739 return self._reversed()

740

741 __str__ = __repr__

742

743 def count(self, latlon):

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

745

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

747 C{(lat, lon)}).

748

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

750

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

752 '''

753 return self._count(latlon)

754

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

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

757

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

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

760

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

762

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

764 '''

765 return self._find(latlon, start_end)

766

767 def _findall(self, latlon, start_end):

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

769 '''

770 try:

771 lat, lon = latlon.lat, latlon.lon

772 except AttributeError:

773 try:

774 lat, lon = latlon

775 except (TypeError, ValueError):

776 raise _IsnotError(_valid_, latlon=latlon)

777

778 _ilat, _ilon = self._ilat, self._ilon

779 _array, _eps = self._array, self._epsilon

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

781 row = _array[i]

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

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

784 yield i

785

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

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

788

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

790 C{(lat, lon)}).

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

792

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

794

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

796 '''

797 return self._findall(latlon, start_end)

798

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

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

801

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

803 C{(lat, lon)}).

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

805

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

807

808 @raise IndexError: Point not found.

809

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

811 '''

812 return self._index(latlon, start_end)

813

814 @Property_RO

815 def ilat(self):

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

817 '''

818 return self._ilat

819

820 @Property_RO

821 def ilon(self):

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

823 '''

824 return self._ilon

825

826# next = __iter__

827

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

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

830

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

832

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

834 '''

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

836

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

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

839

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

841 C{(lat, lon)}).

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

843

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

845

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

847

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

849 '''

850 return self._rfind(latlon, start_end)

851

852 def _slicekwds(self):

853 '''(INTERNAL) Slice kwds.

854 '''

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

856

857 @Property_RO

858 def shape(self):

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

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

861 '''

862 return self._shape

863

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

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

866 '''

867 notImplemented(self, indices)

868

869 def subset(self, indices):

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

871

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

873

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

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

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

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

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

879

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

881

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

883

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

885 nor an C{int}[].

886 '''

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

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

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

890

891 n = len(self)

892 for i, v in enumerate(indices):

893 if not isint(v):

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

895 elif not 0 <= v < n:

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

897

898 return self._subset(indices)

899

900

901 class LatLon2psxy(_Basequence):

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

903 '''

904 _closed = False

905 _len = 0

906 _deg2m = None # default, keep degrees

907 _radius = None

908 _wrap = True

909

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

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

912

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

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

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

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

917

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

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

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

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

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

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

924

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

926

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

928 '''

929 self._closed = closed

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

931 if radius:

932 self._radius = r = Radius(radius)

933 self._deg2m = degrees2m(_1_0, r)

934 if not wrap:

935 self._wrap = False

936

937 def __contains__(self, xy):

938 '''Check for a matching point.

939

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

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

942

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

944

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

946 '''

947 return self._contains(xy)

948

949 def __getitem__(self, index):

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

951 '''

952 return self._getitem(index)

953

954 def __iter__(self):

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

956 '''

957 return self._iter()

958

959 def __len__(self):

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

961 '''

962 return self._len

963

964 def __repr__(self):

965 '''Return a string representation.

966 '''

967 return self._repr()

968

969 def __reversed__(self): # PYCHOK false

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

971 '''

972 return self._reversed()

973

974 __str__ = __repr__

975

976 def count(self, xy):

977 '''Count the number of matching points.

978

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

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

981

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

983

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

985 '''

986 return self._count(xy)

987

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

989 '''Find the first matching point.

990

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

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

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

994

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

996

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

998 '''

999 return self._find(xy, start_end)

1000

1001 def _findall(self, xy, start_end):

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

1003 '''

1004 try:

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

1006

1007 def _x_y_ll3(ll): # match LatLon

1008 return ll.lon, ll.lat, ll

1009

1010 except AttributeError:

1011 try:

1012 x, y = xy[:2]

1013 except (IndexError, TypeError, ValueError):

1014 raise _IsnotError(_valid_, xy=xy)

1015

1016 _x_y_ll3 = self.point # PYCHOK expected

1017

1018 _array, _eps = self._array, self._epsilon

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

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

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

1022 fabs(yi - y) <= _eps:

1023 yield i

1024

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

1026 '''Yield indices of all matching points.

1027

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

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

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

1031

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

1033

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

1035 '''

1036 return self._findall(xy, start_end)

1037

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

1039 '''Find the first matching point.

1040

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

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

1043

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

1045

1046 @raise IndexError: Point not found.

1047

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

1049 '''

1050 return self._index(xy, start_end)

1051

1052 @property_RO

1053 def isPoints2(self):

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

1055 '''

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

1057

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

1059 '''Create a pseudo-xy.

1060

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

1062

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

1064 '''

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

1066 if self._wrap:

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

1068 d = self._deg2m

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

1070 x *= d

1071 y *= d

1072 return Point3Tuple(x, y, ll)

1073

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

1075 '''Find the last matching point.

1076

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

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

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

1080

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

1082

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

1084 '''

1085 return self._rfind(xy, start_end)

1086

1087 def _slicekwds(self):

1088 '''(INTERNAL) Slice kwds.

1089 '''

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

1091

1092

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

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

1095 '''

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

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

1098

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

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

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

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

1103

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

1105

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

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

1108 and C{lon} attributes.

1109

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

1111 are the same or out of range.

1112

1113 @example:

1114

1115 >>> type(array)

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

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

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

1119 >>> type(simply)

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

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

1122 >>> type(sliced)

1123 <class '...Numpy2LatLon'>

1124 '''

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

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

1127 except AttributeError:

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

1129

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

1131 LatLon=LatLon, shape=s)

1132

1133 @property_RO

1134 def isNumpy2(self):

1135 '''Is this a Numpy2 wrapper?

1136 '''

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

1138

1139 def _subset(self, indices):

1140 return self._array[indices] # NumPy special

1141

1142

1143 class Shape2Tuple(_NamedTuple):

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

1145 both C{int}.

1146 '''

1147 _Names_ = (_nrows_, _ncols_)

1148 _Units_ = ( Number_, Number_)

1149

1150

1151 class Tuple2LatLon(_Array2LatLon):

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

1153 '''

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

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

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

1157

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

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

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

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

1162

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

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

1165

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

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

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

1169 attributes.

1170

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

1172 are the same or out of range.

1173

1174 @example:

1175

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

1177 >>> type(tuples)

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

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

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

1181 >>> type(simply)

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

1183 >>> simply

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

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

1186 >>> type(sliced)

1187 <class '...Tuple2LatLon'>

1188 >>> sliced

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

1190

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

1192 >>> closest

1193 LatLon_(lat=1.0, lon=2.0)

1194

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

1196 >>> closest

1197 LatLon_(lat=2.001162, lon=3.001162)

1198 '''

1199 _xinstanceof(list, tuple, tuples=tuples)

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

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

1202 LatLon=LatLon, shape=s)

1203

1204 @property_RO

1205 def isTuple2(self):

1206 '''Is this a Tuple2 wrapper?

1207 '''

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

1209

1210 def _subset(self, indices):

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

1212

1213

1214 def _area2(points, adjust, wrap):

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

1216 '''

1217 if adjust:

1218 # approximate trapezoid by a rectangle, adjusting

1219 # the top width by the cosine of the latitudinal

1220 # average and bottom width by some fudge factor

1221 def _adjust(w, h):

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

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

1224 else:

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

1226 return w * h

1227

1228 # setting radius=1 converts degrees to radians

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

1230 x1, y1, ll = Ps[0]

1231 pts = [ll] # for _areaError

1232

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

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

1235 x2, y2, ll = p

1236 if len(pts) < 4:

1237 pts.append(ll)

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

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

1240 x1, y1 = x2, y2

1241

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

1243

1244

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

1246 '''(INTERNAL) Area issue.

1247 '''

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

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

1250

1251

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

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

1254

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

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

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

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

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

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

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

1262

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

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

1265 is C{None}).

1266

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

1268

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

1270

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

1272

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

1274 ill-suited for regions exceeding several hundred Km

1275 or Miles or with near-polar latitudes.

1276

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

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

1279 '''

1280 if _MODS.booleans.isBoolean(points):

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

1282 else:

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

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

1285

1286

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

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

1289 path or polygon.

1290

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

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

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

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

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

1296

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

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

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

1300

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

1302

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

1304

1305 @see: Function L{quadOf}.

1306

1307 @example:

1308

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

1310 >>> boundsOf(b) # False

1311 >>> 45.0, 1.0, 46.0, 2.0

1312 '''

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

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

1315

1316 v = w

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

1318 if wrap:

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

1320 v = x

1321

1322 if w > x:

1323 w = x

1324 elif e < x:

1325 e = x

1326

1327 if s > y:

1328 s = y

1329 elif n < y:

1330 n = y

1331

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

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

1334

1335

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

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

1338 '''

1339 name, ps = name_points.popitem()

1340 for i, p in enumerate(ps):

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

1342 n = _distanceTo.__name__[1:]

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

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

1345 return ps

1346

1347

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

1349 '''Determine the centroid of a polygon.

1350

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

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

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

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

1355 or C{None}.

1356

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

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

1359

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

1361

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

1363

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

1365 near-zero area.

1366

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

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

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

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

1371 '''

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

1373

1374 # setting radius=1 converts degrees to radians

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

1376 x1, y1, ll = Ps[0]

1377 pts = [ll] # for _areaError

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

1379 x2, y2, ll = p

1380 if len(pts) < 4:

1381 pts.append(ll)

1382 if wrap and not Ps.looped:

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

1384 t = x1 * y2 - x2 * y1

1385 A += t

1386 X += t * (x1 + x2)

1387 Y += t * (y1 + y2)

1388 # XXX more elaborately:

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

1390 # A.fadd_(t1, t2)

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

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

1393 x1, y1 = x2, y2

1394

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

1396 if isnear0(a):

1397 raise _areaError(pts, near_=_near_)

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

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

1400

1401

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

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

1404

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

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

1407 L{Vector3Tuple}[]).

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

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

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

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

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

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

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

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

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

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

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

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

1420 or C{None}.

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

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

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

1424

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

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

1427

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

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

1430

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

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

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

1434

1435 Summarized as follows:

1436

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

1438 # -------+--------+--------+--------------+------

1439 # | | | LatLon2Tuple | favg

1440 # None | None | None | or** |

1441 # | | | Vector3Tuple | favg

1442 # None | LatLon | None | LatLon | favg

1443 # None | None | Vector | Vector | favg

1444 # -------+--------+--------+--------------+------

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

1446 # True | LatLon | None | LatLon | .iTo

1447 # True | None | Vector | Vector | .iTo

1448 # -------+--------+--------+--------------+------

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

1450 # False | LatLon | None | LatLon | .iTo

1451 # False | None | Vector | Vector | .iTo

1452 # _____

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

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

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

1456

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

1458 not subscriptable or not closed.

1459

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

1461 argument.

1462

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

1464 '''

1465 if LatLon and Vector: # PYCHOK no cover

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

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

1468 w = wrap if LatLon else False # intermediateTo

1469 try:

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

1471 raise IndexError

1472 n = _xattr(fi, fin=0)

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

1474 if LatLon:

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

1476 elif Vector:

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

1478 except (IndexError, TypeError):

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

1480 return p

1481

1482

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

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

1485 '''

1486 i = int(fi)

1487 p = points[i]

1488 r = fi - float(i)

1489 if r > EPS: # EPS0?

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

1491 j = i + 1

1492 if fin:

1493 j %= fin

1494 q = points[j]

1495 if r >= EPS1: # PYCHOK no cover

1496 p = q

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

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

1499 elif _isLatLon(p): # backward compatible default

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

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

1502 name=fractional.__name__)

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

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

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

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

1507 name=fractional.__name__)

1508 return p

1509

1510

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

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

1513

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

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

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

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

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

1519

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

1521

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

1523

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

1525

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

1527

1528 @example:

1529

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

1531 >>> isclockwise(f) # False

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

1533 '''

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

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

1536 return True

1537 elif a < 0:

1538 return False

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

1540 raise _areaError(pts)

1541

1542

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

1544 '''Determine whether a polygon is convex.

1545

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

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

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

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

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

1551

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

1553

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

1555

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

1557

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

1559

1560 @example:

1561

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

1563 >>> isconvex(t) # True

1564

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

1566 >>> isconvex(f) # False

1567 '''

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

1569

1570

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

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

1573

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

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

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

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

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

1579

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

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

1582

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

1584

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

1586

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

1588

1589 @example:

1590

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

1592 >>> isconvex_(t) # +1

1593

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

1595 >>> isconvex_(f) # 0

1596 '''

1597 if adjust:

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

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

1600 y = radians(y1 + y2) * _0_5

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

1602 return x21, x2

1603 else:

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

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

1606

1607 c, s = crosserrors(), 0

1608

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

1610 x1, y1, _ = Ps[0]

1611 x2, y2, _ = Ps[1]

1612

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

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

1615 x3, y3, ll = p

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

1617

1618 # get the sign of the distance from point

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

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

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

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

1623 if s < 0: # non-convex

1624 return 0

1625 s = +1

1626

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

1628 if s > 0: # non-convex

1629 return 0

1630 s = -1

1631

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

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

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

1635 t = Fmt.SQUARE(points=i)

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

1637

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

1639

1640 return s # all points on the same side

1641

1642

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

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

1645

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

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

1648 or L{BooleanGH}).

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

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

1651

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

1653 composite, C{False} otherwise.

1654

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

1656

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

1658

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

1660

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

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

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

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

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

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

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

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

1669 '''

1670 try:

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

1672 except AttributeError:

1673 try:

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

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

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

1677

1678 if wrap:

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

1680

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

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

1683 return dx, x2, y2

1684

1685 else:

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

1687 x0_180_ = x0 - _180_0

1688 x0_180 = x0 + _180_0

1689

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

1691 x = _umod_360(float(x))

1692 if x < x0_180_:

1693 x += _360_0

1694 elif x >= x0_180:

1695 x -= _360_0

1696 return (x - x1), x, y

1697

1698 if _MODS.booleans.isBoolean(points):

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

1700

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

1702 p = Ps[0]

1703 e = m = False

1704 S = Fsum()

1705

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

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

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

1709 # ignore duplicate and near-duplicate pts

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

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

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

1713 # edges on boundary as more than one crossing

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

1715 m = not m

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

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

1718 e = not e

1719

1720 S += sin(radians(y2))

1721 x1, y1 = x2, y2

1722

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

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

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

1726 # contain the North Pole, flip the result.

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

1728 e = not e

1729 return e

1730

1731

1732 def ispolar(points, wrap=False):

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

1734

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

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

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

1738

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

1740 otherwise.

1741

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

1743

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

1745 have C{bearingTo2}, C{initialBearingTo}

1746 and C{finalBearingTo} methods.

1747 '''

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

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

1750 p2, p1 = Ps[0:2]

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

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

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

1754 if w and not Ps.looped:

1755 p2 = _unrollon(p1, p2)

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

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

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

1759 p1, b1 = p2, b2

1760

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

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

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

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

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

1766

1767

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

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

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

1771

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

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

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

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

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

1777 keyword arguments.

1778

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

1780 the lune shape.

1781

1782 @see: U{Latitude-longitude quadrangle

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

1784 '''

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

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

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

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

1789 if closed:

1790 t += t[:1]

1791 return t

1792

1793

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

1795 limit=9, **LatLon_and_kwds):

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

1797

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

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

1800

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

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

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

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

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

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

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

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

1809 spherical earth radius L{R_KM}).

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

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

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

1813

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

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

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

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

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

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

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

1821

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

1823 see function L{pygeodesy.equirectangular_}.

1824

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

1826

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

1828

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

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

1831

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

1833 '''

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

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

1836 # equirectangular_ returns a Distance4Tuple(distance

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

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

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

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

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

1842

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

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

1845

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

1847 #

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

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

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

1851 #

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

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

1854 # w = (y * dy + x * dx) / hypot(dx, dy)

1855 # h = (y * dx - x * dy) / hypot(dx, dy)

1856 #

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

1858 # xc = dx * w / hypot(dx, dy)

1859 # yc = dy * w / hypot(dx, dy)

1860 # or

1861 # xc = dx * f

1862 # yc = dy * f

1863 # with

1864 # f = w / hypot(dx, dy)

1865 # or

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

1867 #

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

1869

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

1871 p1 = c = Ps[0]

1872 u1 = u = _0_0

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

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

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

1876 # iff wrapped, unroll lon1 (actually previous

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

1878 if wrap:

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

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

1881 if d21 > EPS:

1882 # distance point to p1, y01 and x01 negated

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

1884 if d2 > EPS:

1885 w2 = y01 * y21 + x01 * x21

1886 if w2 > 0:

1887 if w2 < d21:

1888 # closest is between p1 and p2, use

1889 # original delta's, not y21 and x21

1890 f = w2 / d21

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

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

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

1894 u1 = _0_0

1895 else: # p2 is closest

1896 p1, u1 = p2, u2

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

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

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

1900 p1, u1 = p2, u2

1901

1902 a = atan2b(dx, dy)

1903 d = hypot(dx, dy)

1904 h = _h(c)

1905 n = nameof(point) or nearestOn5.__name__

1906 if LatLon_and_kwds:

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

1908 LL = _xkwds_pop(kwds, LatLon=None)

1909 if LL is not None:

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

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

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

1913

1914

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

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

1917

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

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

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

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

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

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

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

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

1926

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

1928 B{C{radius}}).

1929

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

1931

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

1933

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

1935 C{B{points}} a composite.

1936

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

1938 distance approximation and is ill-suited for regions exceeding

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

1940

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

1942 L{ellipsoidalKarney.perimeterOf}.

1943 '''

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

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

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

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

1948 if w and c:

1949 w = not Ps.looped

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

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

1952 p2.y, p2.x,

1953 adjust=a, limit=None,

1954 wrap=w) # PYCHOK non-seq

1955 yield hypot(dx, dy)

1956 p1 = p2

1957

1958 if _MODS.booleans.isBoolean(points):

1959 if not closed:

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

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

1962 radius=radius, wrap=wrap)

1963 else:

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

1965 return degrees2m(d, radius=radius)

1966

1967

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

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

1970

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

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

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

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

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

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

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

1978 keyword arguments.

1979

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

1981 outlining the quadrilateral.

1982

1983 @see: Function L{boundsOf}.

1984 '''

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

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

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

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

1989 if closed:

1990 t += t[:1]

1991 return t

1992

1993

1994__all__ += _ALL_DOCS(_Array2LatLon, _Basequence)

1995

1996# **) MIT License

1997#

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

1999#

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

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

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

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

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

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

2006#

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

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

2009#

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

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

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

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

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

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

2016# OTHER DEALINGS IN THE SOFTWARE.