Coverage for pygeodesy/points.py: 92%

610 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-09-01 13:41 -0400

1 

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

3 

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

5 

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

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

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

9 

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

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

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

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

14 

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

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

17 

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

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

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

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

22 

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

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

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

26''' 

27 

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

29 issubclassof, _Sequence, _xcopy, _xdup, \ 

30 _xinstanceof 

31from pygeodesy.constants import EPS, EPS1, PI_2, R_M, isnear0, isnear1, \ 

32 _umod_360, _0_0, _0_5, _1_0, _2_0, _6_0, \ 

33 _90_0, _N_90_0, _180_0, _360_0 

34# from pygeodesy.datums import _spherical_datum # from .formy 

35from pygeodesy.dms import F_D, latDMS, lonDMS, parseDMS, S_DEG, S_DMS, \ 

36 S_MIN, S_SEC, S_SEP 

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

38 _IsnotError, _TypeError, _ValueError, \ 

39 _xattr, _xkwds, _xkwds_pop 

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

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

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

43 isnormal, latlon2n_xyz, normal, _spherical_datum 

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

45 _DEQUALSPACED_, _ELLIPSIS_, _EW_, _immutable_, \ 

46 _lat_, _lon_, _LatLon_, _near_, _no_, _not_, \ 

47 _NS_, _point_, _SPACE_, _UNDER_, _valid_ 

48from pygeodesy.iters import LatLon2PsxyIter, PointsIter, points2 

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

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

51 _NamedTuple, _xnamed, _xother3, _xotherError 

52from pygeodesy.namedTuples import Bounds2Tuple, Bounds4Tuple, \ 

53 LatLon2Tuple, NearestOn3Tuple, \ 

54 NearestOn5Tuple, PhiLam2Tuple, \ 

55 Point3Tuple, Vector3Tuple, Vector4Tuple 

56from pygeodesy.nvectorBase import NvectorBase, _N_vector_ 

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

58 property_doc_, property_RO, _update_all 

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

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

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

62 unroll180, _unrollon, unrollPI, _Wrap, wrap180 

63 

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

65 

66__all__ = _ALL_LAZY.points 

67__version__ = '23.08.23' 

68 

69_ilat_ = 'ilat' 

70_ilon_ = 'ilon' 

71_ncols_ = 'ncols' 

72_nrows_ = 'nrows' 

73 

74 

75class LatLon_(object): # XXX in heights._HeightBase.height 

76 '''Low-overhead C{LatLon} class for L{Numpy2LatLon} and L{Tuple2LatLon}. 

77 ''' 

78 # __slots__ efficiency is voided if the __slots__ class attribute 

79 # is used in a subclass of a class with the traditional __dict__, 

80 # see <https://docs.Python.org/2/reference/datamodel.html#slots> 

81 # and __slots__ must be repeated in sub-classes, see Luciano 

82 # Ramalho, "Fluent Python", O'Reilly, 2016 p. 276+ "Problems 

83 # with __slots__" (also at <https://Books.Google.ie/books? 

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

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

86 # 2022 p. 390 "Summarizing the Issues with __slots__". 

87 # 

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

89 # Property_RO = property_RO # no __dict__ with __slots__! 

90 # 

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

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

93 

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

95 datum=None, wrap=False): 

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

97 

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

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

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

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

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

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

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

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

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

107 (C{bool}). 

108 

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

110 

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

112 un-clipped and un-validated . 

113 ''' 

114 if lon is None: # PYCHOK no cover 

115 try: 

116 ll = latlonh.lat, latlonh.lon 

117 height = _xattr(latlonh, height=height) 

118 except AttributeError: 

119 raise _IsnotError(_LatLon_, latlonh=latlonh) 

120 if wrap: 

121 ll = _Wrap.latlon(*ll) 

122 elif wrap: # PYCHOK no cover 

123 ll = _Wrap.latlonDMS2(latlonh, lon) 

124 else: # must be latNS, lonEW 

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

126 

127 self.lat, self.lon = ll 

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

129 self.height = height 

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

131 _spherical_datum(datum, name=self.name) 

132 

133 def __eq__(self, other): 

134 return isinstance(other, LatLon_) and \ 

135 other.lat == self.lat and \ 

136 other.lon == self.lon 

137 

138 def __ne__(self, other): 

139 return not self.__eq__(other) 

140 

141 def __repr__(self): 

142 return self.toRepr() 

143 

144 def __str__(self): 

145 return self.toStr() 

146 

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

148 '''Instantiate this very class. 

149 

150 @arg args: Optional, positional arguments. 

151 @kwarg kwds: Optional, keyword arguments. 

152 

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

154 ''' 

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

156 

157 def copy(self, deep=False): 

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

159 

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

161 shallow copy (C{bool}). 

162 

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

164 ''' 

165 return _xcopy(self, deep=deep) 

166 

167 def dup(self, **items): 

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

169 

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

171 

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

173 

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

175 ''' 

176 d = _xdup(self, **items) 

177 if items: 

178 _update_all(d) 

179 return d 

180 

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

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

183 

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

185 

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

187 ''' 

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

189 

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

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

192 and an other point. 

193 

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

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

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

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

198 intermediate height. 

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

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

201 

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

203 

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

205 ''' 

206 f = Scalar(fraction=fraction) 

207 if isnear0(f): 

208 r = self 

209 elif isnear1(f) and not wrap: 

210 r = self.others(other) 

211 else: 

212 p = self.others(other) 

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

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

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

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

217 height=h, datum=self.datum, 

218 name=self.intermediateTo.__name__) 

219 return r 

220 

221 @Property_RO # PYCHOK no cover 

222 def isEllipsoidal(self): 

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

224 ''' 

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

226 

227 @Property_RO # PYCHOK no cover 

228 def isEllipsoidalLatLon(self): 

229 '''Get C{LatLon} base. 

230 ''' 

231 return False 

232 

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

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

235 

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

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

238 

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

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

241 

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

243 ''' 

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

245 

246 @Property_RO 

247 def isnormal(self): 

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

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

250 

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

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

253 ''' 

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

255 

256 @Property_RO 

257 def isSpherical(self): # PYCHOK no cover 

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

259 ''' 

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

261 

262 @Property_RO 

263 def lam(self): 

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

265 ''' 

266 return radians(self.lon) 

267 

268 @Property 

269 def latlon(self): 

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

271 ''' 

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

273 

274 @latlon.setter # PYCHOK setter! 

275 def latlon(self, latlon): 

276 '''Set the lat- and longitude. 

277 

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

279 ''' 

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

281 _update_all(self) 

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

283 

284 @Property_RO 

285 def latlonheight(self): 

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

287 ''' 

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

289 

290 @Property_RO 

291 def _N_vector(self): 

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

293 ''' 

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

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

296 

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

298 '''Refined class comparison. 

299 

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

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

302 keyword arguments. 

303 

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

305 

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

307 ''' 

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

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

310 and hasattr(other, _lon_)): 

311 return other 

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

313 

314 def normal(self): 

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

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

317 

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

319 wasn't (but is now). 

320 

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

322 ''' 

323 n = self.isnormal 

324 if not n: 

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

326 return n 

327 

328 @Property_RO 

329 def phi(self): 

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

331 ''' 

332 return radians(self.lat) 

333 

334 @Property_RO 

335 def philam(self): 

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

337 ''' 

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

339 

340 @Property_RO 

341 def philamheight(self): 

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

343 ''' 

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

345 

346 @deprecated_method 

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

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

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

350 

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

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

353 

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

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

356 ignoring any duplicate or closing final 

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

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

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

360 

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

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

363 

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

365 

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

367 ''' 

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

369 

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

371 '''Return a points iterator. 

372 

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

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

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

376 

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

378 

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

380 ''' 

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

382 

383 @deprecated_method 

384 def to2ab(self): # PYCHOK no cover 

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

386 return self.philam 

387 

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

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

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

391 

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

393 shallow copy (C{bool}). 

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

395 

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

397 optionally renamed (C{LatLon}). 

398 

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

400 function L{pygeodesy.normal}. 

401 ''' 

402 ll = self.copy(deep=deep) 

403 _ = ll.normal() 

404 if name: 

405 ll.name = str(name) 

406 return ll 

407 

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

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

410 surface) components, I{including height}. 

411 

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

413 (C{meter}). 

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

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

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

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

418 

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

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

421 y, z, h)}. 

422 

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

424 argument. 

425 ''' 

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

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

428 if Nvector is not None: 

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

430 return _xnamed(r, self.name) 

431 

432 def toRepr(self, **kwds): 

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

434 

435 @kwarg kwds: Optional, keyword arguments. 

436 

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

438 ''' 

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

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

441 

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

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

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

445 

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

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

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

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

450 non-joined). 

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

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

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

454 or L{pygeodesy.lonDMS}. 

455 

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

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

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

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

460 respectively non-empty). 

461 

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

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

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

465 ''' 

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

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

468 

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

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

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

472 if self.height: 

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

474 if self.name: 

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

476 if kwds: 

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

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

479 

480 @deprecated_method 

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

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

483 return self.toRepr(**kwds) 

484 

485 

486def _isLatLon(inst): 

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

488 ''' 

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

490 

491 

492def _isLatLon_(LL): 

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

494 ''' 

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

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

497 

498 

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

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

501 

502 

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

504 '''(INTERNAL) Base class. 

505 ''' 

506 _array = [] 

507 _epsilon = EPS 

508 _itemname = _point_ 

509 

510 def _contains(self, point): 

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

512 ''' 

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

514 

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

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

517 

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

519 shallow copy (C{bool}). 

520 

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

522 ''' 

523 return _xcopy(self, deep=deep) 

524 

525 def _count(self, point): 

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

527 ''' 

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

529 

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

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

532 

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

534 

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

536 

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

538 ''' 

539 if items: 

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

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

542 return _xdup(self) 

543 

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

545 def epsilon(self): 

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

547 ''' 

548 return self._epsilon 

549 

550 @epsilon.setter # PYCHOK setter! 

551 def epsilon(self, tol): 

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

553 

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

555 ''' 

556 self._epsilon = Scalar_(tolerance=tol) 

557 

558 def _find(self, point, start_end): 

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

560 ''' 

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

562 return i 

563 return -1 

564 

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

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

567 ''' 

568 notImplemented(self, point, start_end) 

569 

570 def _getitem(self, index): 

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

572 ''' 

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

574 if isinstance(index, slice): 

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

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

577 else: 

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

579 

580 def _index(self, point, start_end): 

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

582 ''' 

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

584 return i 

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

586 

587 @property_RO 

588 def isNumpy2(self): # PYCHOK no cover 

589 '''Is this a Numpy2 wrapper? 

590 ''' 

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

592 

593 @property_RO 

594 def isPoints2(self): # PYCHOK no cover 

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

596 ''' 

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

598 

599 @property_RO 

600 def isTuple2(self): # PYCHOK no cover 

601 '''Is this a Tuple2 wrapper? 

602 ''' 

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

604 

605 def _iter(self): 

606 '''(INTERNAL) Yield all points. 

607 ''' 

608 _array, _point = self._array, self.point 

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

610 yield _point(_array[i]) 

611 

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

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

614 

615 @arg attrs: Optional arguments. 

616 ''' 

617 notOverloaded(self, *attrs) 

618 

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

620 '''(INTERNAL) Return the range. 

621 ''' 

622 if step > 0: 

623 if start is None: 

624 start = 0 

625 if end is None: 

626 end = len(self) 

627 elif step < 0: 

628 if start is None: 

629 start = len(self) - 1 

630 if end is None: 

631 end = -1 

632 else: 

633 raise _ValueError(step=step) 

634 return range(start, end, step) 

635 

636 def _repr(self): 

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

638 ''' 

639 # XXX use Python 3+ reprlib.repr 

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

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

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

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

644 

645 def _reversed(self): # PYCHOK false 

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

647 ''' 

648 _array, point = self._array, self.point 

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

650 yield point(_array[i]) 

651 

652 def _rfind(self, point, start_end): 

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

654 ''' 

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

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

657 

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

659 return i 

660 return -1 

661 

662 def _slicekwds(self): # PYCHOK no cover 

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

664 ''' 

665 return {} 

666 

667 

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

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

670 ''' 

671 _array = () 

672 _ilat = 0 # row column index 

673 _ilon = 0 # row column index 

674 _LatLon = LatLon_ # default 

675 _shape = () 

676 

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

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

679 ''' 

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

681 

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

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

684 

685 self._array = array 

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

687 

688 if LatLon: # check the point class 

689 if not _isLatLon_(LatLon): 

690 raise _IsnotError(_valid_, LatLon=LatLon) 

691 self._LatLon = LatLon 

692 

693 # check the attr indices 

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

695 if not isint(i): 

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

697 i = int(i) 

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

699 raise _ValueError(ai, i) 

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

701 if int(j) == i: 

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

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

704 

705 def __contains__(self, latlon): 

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

707 

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

709 C{(lat, lon)}). 

710 

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

712 

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

714 ''' 

715 return self._contains(latlon) 

716 

717 def __getitem__(self, index): 

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

719 ''' 

720 return self._getitem(index) 

721 

722 def __iter__(self): 

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

724 ''' 

725 return self._iter() 

726 

727 def __len__(self): 

728 '''Return the number of rows. 

729 ''' 

730 return self._shape[0] 

731 

732 def __repr__(self): 

733 '''Return a string representation. 

734 ''' 

735 return self._repr() 

736 

737 def __reversed__(self): # PYCHOK false 

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

739 ''' 

740 return self._reversed() 

741 

742 __str__ = __repr__ 

743 

744 def count(self, latlon): 

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

746 

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

748 C{(lat, lon)}). 

749 

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

751 

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

753 ''' 

754 return self._count(latlon) 

755 

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

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

758 

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

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

761 

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

763 

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

765 ''' 

766 return self._find(latlon, start_end) 

767 

768 def _findall(self, latlon, start_end): 

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

770 ''' 

771 try: 

772 lat, lon = latlon.lat, latlon.lon 

773 except AttributeError: 

774 try: 

775 lat, lon = latlon 

776 except (TypeError, ValueError): 

777 raise _IsnotError(_valid_, latlon=latlon) 

778 

779 _ilat, _ilon = self._ilat, self._ilon 

780 _array, _eps = self._array, self._epsilon 

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

782 row = _array[i] 

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

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

785 yield i 

786 

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

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

789 

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

791 C{(lat, lon)}). 

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

793 

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

795 

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

797 ''' 

798 return self._findall(latlon, start_end) 

799 

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

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

802 

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

804 C{(lat, lon)}). 

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

806 

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

808 

809 @raise IndexError: Point not found. 

810 

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

812 ''' 

813 return self._index(latlon, start_end) 

814 

815 @Property_RO 

816 def ilat(self): 

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

818 ''' 

819 return self._ilat 

820 

821 @Property_RO 

822 def ilon(self): 

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

824 ''' 

825 return self._ilon 

826 

827# next = __iter__ 

828 

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

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

831 

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

833 

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

835 ''' 

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

837 

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

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

840 

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

842 C{(lat, lon)}). 

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

844 

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

846 

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

848 

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

850 ''' 

851 return self._rfind(latlon, start_end) 

852 

853 def _slicekwds(self): 

854 '''(INTERNAL) Slice kwds. 

855 ''' 

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

857 

858 @Property_RO 

859 def shape(self): 

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

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

862 ''' 

863 return self._shape 

864 

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

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

867 ''' 

868 notImplemented(self, indices) 

869 

870 def subset(self, indices): 

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

872 

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

874 

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

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

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

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

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

880 

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

882 

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

884 

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

886 nor an C{int}[]. 

887 ''' 

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

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

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

891 

892 n = len(self) 

893 for i, v in enumerate(indices): 

894 if not isint(v): 

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

896 elif not 0 <= v < n: 

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

898 

899 return self._subset(indices) 

900 

901 

902class LatLon2psxy(_Basequence): 

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

904 ''' 

905 _closed = False 

906 _len = 0 

907 _deg2m = None # default, keep degrees 

908 _radius = None 

909 _wrap = True 

910 

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

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

913 

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

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

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

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

918 

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

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

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

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

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

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

925 

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

927 

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

929 ''' 

930 self._closed = closed 

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

932 if radius: 

933 self._radius = r = Radius(radius) 

934 self._deg2m = degrees2m(_1_0, r) 

935 if not wrap: 

936 self._wrap = False 

937 

938 def __contains__(self, xy): 

939 '''Check for a matching point. 

940 

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

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

943 

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

945 

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

947 ''' 

948 return self._contains(xy) 

949 

950 def __getitem__(self, index): 

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

952 ''' 

953 return self._getitem(index) 

954 

955 def __iter__(self): 

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

957 ''' 

958 return self._iter() 

959 

960 def __len__(self): 

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

962 ''' 

963 return self._len 

964 

965 def __repr__(self): 

966 '''Return a string representation. 

967 ''' 

968 return self._repr() 

969 

970 def __reversed__(self): # PYCHOK false 

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

972 ''' 

973 return self._reversed() 

974 

975 __str__ = __repr__ 

976 

977 def count(self, xy): 

978 '''Count the number of matching points. 

979 

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

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

982 

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

984 

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

986 ''' 

987 return self._count(xy) 

988 

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

990 '''Find the first matching point. 

991 

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

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

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

995 

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

997 

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

999 ''' 

1000 return self._find(xy, start_end) 

1001 

1002 def _findall(self, xy, start_end): 

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

1004 ''' 

1005 try: 

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

1007 

1008 def _x_y_ll3(ll): # match LatLon 

1009 return ll.lon, ll.lat, ll 

1010 

1011 except AttributeError: 

1012 try: 

1013 x, y = xy[:2] 

1014 except (IndexError, TypeError, ValueError): 

1015 raise _IsnotError(_valid_, xy=xy) 

1016 

1017 _x_y_ll3 = self.point # PYCHOK expected 

1018 

1019 _array, _eps = self._array, self._epsilon 

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

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

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

1023 fabs(yi - y) <= _eps: 

1024 yield i 

1025 

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

1027 '''Yield indices of all matching points. 

1028 

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

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

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

1032 

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

1034 

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

1036 ''' 

1037 return self._findall(xy, start_end) 

1038 

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

1040 '''Find the first matching point. 

1041 

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

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

1044 

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

1046 

1047 @raise IndexError: Point not found. 

1048 

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

1050 ''' 

1051 return self._index(xy, start_end) 

1052 

1053 @property_RO 

1054 def isPoints2(self): 

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

1056 ''' 

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

1058 

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

1060 '''Create a pseudo-xy. 

1061 

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

1063 

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

1065 ''' 

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

1067 if self._wrap: 

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

1069 d = self._deg2m 

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

1071 x *= d 

1072 y *= d 

1073 return Point3Tuple(x, y, ll) 

1074 

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

1076 '''Find the last matching point. 

1077 

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

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

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

1081 

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

1083 

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

1085 ''' 

1086 return self._rfind(xy, start_end) 

1087 

1088 def _slicekwds(self): 

1089 '''(INTERNAL) Slice kwds. 

1090 ''' 

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

1092 

1093 

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

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

1096 ''' 

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

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

1099 

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

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

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

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

1104 

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

1106 

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

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

1109 and C{lon} attributes. 

1110 

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

1112 are the same or out of range. 

1113 

1114 @example: 

1115 

1116 >>> type(array) 

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

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

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

1120 >>> type(simply) 

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

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

1123 >>> type(sliced) 

1124 <class '...Numpy2LatLon'> 

1125 ''' 

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

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

1128 except AttributeError: 

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

1130 

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

1132 LatLon=LatLon, shape=s) 

1133 

1134 @property_RO 

1135 def isNumpy2(self): 

1136 '''Is this a Numpy2 wrapper? 

1137 ''' 

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

1139 

1140 def _subset(self, indices): 

1141 return self._array[indices] # NumPy special 

1142 

1143 

1144class Shape2Tuple(_NamedTuple): 

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

1146 both C{int}. 

1147 ''' 

1148 _Names_ = (_nrows_, _ncols_) 

1149 _Units_ = ( Number_, Number_) 

1150 

1151 

1152class Tuple2LatLon(_Array2LatLon): 

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

1154 ''' 

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

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

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

1158 

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

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

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

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

1163 

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

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

1166 

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

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

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

1170 attributes. 

1171 

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

1173 are the same or out of range. 

1174 

1175 @example: 

1176 

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

1178 >>> type(tuples) 

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

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

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

1182 >>> type(simply) 

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

1184 >>> simply 

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

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

1187 >>> type(sliced) 

1188 <class '...Tuple2LatLon'> 

1189 >>> sliced 

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

1191 

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

1193 >>> closest 

1194 LatLon_(lat=1.0, lon=2.0) 

1195 

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

1197 >>> closest 

1198 LatLon_(lat=2.001162, lon=3.001162) 

1199 ''' 

1200 _xinstanceof(list, tuple, tuples=tuples) 

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

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

1203 LatLon=LatLon, shape=s) 

1204 

1205 @property_RO 

1206 def isTuple2(self): 

1207 '''Is this a Tuple2 wrapper? 

1208 ''' 

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

1210 

1211 def _subset(self, indices): 

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

1213 

1214 

1215def _area2(points, adjust, wrap): 

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

1217 ''' 

1218 if adjust: 

1219 # approximate trapezoid by a rectangle, adjusting 

1220 # the top width by the cosine of the latitudinal 

1221 # average and bottom width by some fudge factor 

1222 def _adjust(w, h): 

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

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

1225 else: 

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

1227 return w * h 

1228 

1229 # setting radius=1 converts degrees to radians 

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

1231 x1, y1, ll = Ps[0] 

1232 pts = [ll] # for _areaError 

1233 

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

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

1236 x2, y2, ll = p 

1237 if len(pts) < 4: 

1238 pts.append(ll) 

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

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

1241 x1, y1 = x2, y2 

1242 

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

1244 

1245 

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

1247 '''(INTERNAL) Area issue. 

1248 ''' 

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

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

1251 

1252 

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

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

1255 

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

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

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

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

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

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

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

1263 

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

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

1266 is C{None}). 

1267 

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

1269 

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

1271 

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

1273 

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

1275 ill-suited for regions exceeding several hundred Km 

1276 or Miles or with near-polar latitudes. 

1277 

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

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

1280 ''' 

1281 if _MODS.booleans.isBoolean(points): 

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

1283 else: 

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

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

1286 

1287 

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

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

1290 path or polygon. 

1291 

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

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

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

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

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

1297 

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

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

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

1301 

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

1303 

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

1305 

1306 @see: Function L{quadOf}. 

1307 

1308 @example: 

1309 

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

1311 >>> boundsOf(b) # False 

1312 >>> 45.0, 1.0, 46.0, 2.0 

1313 ''' 

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

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

1316 

1317 v = w 

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

1319 if wrap: 

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

1321 v = x 

1322 

1323 if w > x: 

1324 w = x 

1325 elif e < x: 

1326 e = x 

1327 

1328 if s > y: 

1329 s = y 

1330 elif n < y: 

1331 n = y 

1332 

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

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

1335 

1336 

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

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

1339 ''' 

1340 name, ps = name_points.popitem() 

1341 for i, p in enumerate(ps): 

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

1343 n = _distanceTo.__name__[1:] 

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

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

1346 return ps 

1347 

1348 

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

1350 '''Determine the centroid of a polygon. 

1351 

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

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

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

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

1356 or C{None}. 

1357 

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

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

1360 

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

1362 

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

1364 

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

1366 near-zero area. 

1367 

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

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

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

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

1372 ''' 

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

1374 

1375 # setting radius=1 converts degrees to radians 

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

1377 x1, y1, ll = Ps[0] 

1378 pts = [ll] # for _areaError 

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

1380 x2, y2, ll = p 

1381 if len(pts) < 4: 

1382 pts.append(ll) 

1383 if wrap and not Ps.looped: 

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

1385 t = x1 * y2 - x2 * y1 

1386 A += t 

1387 X += t * (x1 + x2) 

1388 Y += t * (y1 + y2) 

1389 # XXX more elaborately: 

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

1391 # A.fadd_(t1, t2) 

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

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

1394 x1, y1 = x2, y2 

1395 

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

1397 if isnear0(a): 

1398 raise _areaError(pts, near_=_near_) 

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

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

1401 

1402 

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

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

1405 

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

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

1408 L{Vector3Tuple}[]). 

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

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

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

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

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

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

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

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

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

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

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

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

1421 or C{None}. 

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

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

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

1425 

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

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

1428 

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

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

1431 

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

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

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

1435 

1436 Summarized as follows: 

1437 

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

1439 # -------+--------+--------+--------------+------ 

1440 # | | | LatLon2Tuple | favg 

1441 # None | None | None | or** | 

1442 # | | | Vector3Tuple | favg 

1443 # None | LatLon | None | LatLon | favg 

1444 # None | None | Vector | Vector | favg 

1445 # -------+--------+--------+--------------+------ 

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

1447 # True | LatLon | None | LatLon | .iTo 

1448 # True | None | Vector | Vector | .iTo 

1449 # -------+--------+--------+--------------+------ 

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

1451 # False | LatLon | None | LatLon | .iTo 

1452 # False | None | Vector | Vector | .iTo 

1453 # _____ 

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

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

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

1457 

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

1459 not subscriptable or not closed. 

1460 

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

1462 argument. 

1463 

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

1465 ''' 

1466 if LatLon and Vector: # PYCHOK no cover 

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

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

1469 w = wrap if LatLon else False # intermediateTo 

1470 try: 

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

1472 raise IndexError 

1473 n = _xattr(fi, fin=0) 

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

1475 if LatLon: 

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

1477 elif Vector: 

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

1479 except (IndexError, TypeError): 

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

1481 return p 

1482 

1483 

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

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

1486 ''' 

1487 i = int(fi) 

1488 p = points[i] 

1489 r = fi - float(i) 

1490 if r > EPS: # EPS0? 

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

1492 j = i + 1 

1493 if fin: 

1494 j %= fin 

1495 q = points[j] 

1496 if r >= EPS1: # PYCHOK no cover 

1497 p = q 

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

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

1500 elif _isLatLon(p): # backward compatible default 

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

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

1503 name=fractional.__name__) 

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

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

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

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

1508 name=fractional.__name__) 

1509 return p 

1510 

1511 

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

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

1514 

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

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

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

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

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

1520 

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

1522 

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

1524 

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

1526 

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

1528 

1529 @example: 

1530 

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

1532 >>> isclockwise(f) # False 

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

1534 ''' 

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

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

1537 return True 

1538 elif a < 0: 

1539 return False 

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

1541 raise _areaError(pts) 

1542 

1543 

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

1545 '''Determine whether a polygon is convex. 

1546 

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

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

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

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

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

1552 

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

1554 

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

1556 

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

1558 

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

1560 

1561 @example: 

1562 

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

1564 >>> isconvex(t) # True 

1565 

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

1567 >>> isconvex(f) # False 

1568 ''' 

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

1570 

1571 

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

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

1574 

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

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

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

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

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

1580 

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

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

1583 

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

1585 

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

1587 

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

1589 

1590 @example: 

1591 

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

1593 >>> isconvex_(t) # +1 

1594 

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

1596 >>> isconvex_(f) # 0 

1597 ''' 

1598 if adjust: 

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

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

1601 y = radians(y1 + y2) * _0_5 

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

1603 return x21, x2 

1604 else: 

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

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

1607 

1608 c, s = crosserrors(), 0 

1609 

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

1611 x1, y1, _ = Ps[0] 

1612 x2, y2, _ = Ps[1] 

1613 

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

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

1616 x3, y3, ll = p 

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

1618 

1619 # get the sign of the distance from point 

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

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

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

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

1624 if s < 0: # non-convex 

1625 return 0 

1626 s = +1 

1627 

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

1629 if s > 0: # non-convex 

1630 return 0 

1631 s = -1 

1632 

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

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

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

1636 t = Fmt.SQUARE(points=i) 

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

1638 

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

1640 

1641 return s # all points on the same side 

1642 

1643 

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

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

1646 

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

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

1649 or L{BooleanGH}). 

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

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

1652 

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

1654 composite, C{False} otherwise. 

1655 

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

1657 

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

1659 

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

1661 

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

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

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

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

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

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

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

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

1670 ''' 

1671 try: 

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

1673 except AttributeError: 

1674 try: 

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

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

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

1678 

1679 if wrap: 

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

1681 

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

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

1684 return dx, x2, y2 

1685 

1686 else: 

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

1688 x0_180_ = x0 - _180_0 

1689 x0_180 = x0 + _180_0 

1690 

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

1692 x = _umod_360(float(x)) 

1693 if x < x0_180_: 

1694 x += _360_0 

1695 elif x >= x0_180: 

1696 x -= _360_0 

1697 return (x - x1), x, y 

1698 

1699 if _MODS.booleans.isBoolean(points): 

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

1701 

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

1703 p = Ps[0] 

1704 e = m = False 

1705 S = Fsum() 

1706 

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

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

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

1710 # ignore duplicate and near-duplicate pts 

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

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

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

1714 # edges on boundary as more than one crossing 

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

1716 m = not m 

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

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

1719 e = not e 

1720 

1721 S += sin(radians(y2)) 

1722 x1, y1 = x2, y2 

1723 

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

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

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

1727 # contain the North Pole, flip the result. 

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

1729 e = not e 

1730 return e 

1731 

1732 

1733def ispolar(points, wrap=False): 

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

1735 

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

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

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

1739 

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

1741 otherwise. 

1742 

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

1744 

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

1746 have C{bearingTo2}, C{initialBearingTo} 

1747 and C{finalBearingTo} methods. 

1748 ''' 

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

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

1751 p2, p1 = Ps[0:2] 

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

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

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

1755 if w and not Ps.looped: 

1756 p2 = _unrollon(p1, p2) 

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

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

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

1760 p1, b1 = p2, b2 

1761 

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

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

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

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

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

1767 

1768 

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

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

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

1772 

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

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

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

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

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

1778 keyword arguments. 

1779 

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

1781 the lune shape. 

1782 

1783 @see: U{Latitude-longitude quadrangle 

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

1785 ''' 

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

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

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

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

1790 if closed: 

1791 t += t[:1] 

1792 return t 

1793 

1794 

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

1796 limit=9, **LatLon_and_kwds): 

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

1798 

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

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

1801 

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

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

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

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

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

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

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

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

1810 spherical earth radius L{R_KM}). 

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

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

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

1814 

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

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

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

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

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

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

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

1822 

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

1824 see function L{pygeodesy.equirectangular_}. 

1825 

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

1827 

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

1829 

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

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

1832 

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

1834 ''' 

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

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

1837 # equirectangular_ returns a Distance4Tuple(distance 

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

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

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

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

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

1843 

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

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

1846 

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

1848 # 

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

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

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

1852 # 

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

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

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

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

1857 # 

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

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

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

1861 # or 

1862 # xc = dx * f 

1863 # yc = dy * f 

1864 # with 

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

1866 # or 

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

1868 # 

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

1870 

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

1872 p1 = c = Ps[0] 

1873 u1 = u = _0_0 

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

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

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

1877 # iff wrapped, unroll lon1 (actually previous 

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

1879 if wrap: 

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

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

1882 if d21 > EPS: 

1883 # distance point to p1, y01 and x01 negated 

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

1885 if d2 > EPS: 

1886 w2 = y01 * y21 + x01 * x21 

1887 if w2 > 0: 

1888 if w2 < d21: 

1889 # closest is between p1 and p2, use 

1890 # original delta's, not y21 and x21 

1891 f = w2 / d21 

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

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

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

1895 u1 = _0_0 

1896 else: # p2 is closest 

1897 p1, u1 = p2, u2 

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

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

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

1901 p1, u1 = p2, u2 

1902 

1903 a = atan2b(dx, dy) 

1904 d = hypot(dx, dy) 

1905 h = _h(c) 

1906 n = nameof(point) or nearestOn5.__name__ 

1907 if LatLon_and_kwds: 

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

1909 LL = _xkwds_pop(kwds, LatLon=None) 

1910 if LL is not None: 

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

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

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

1914 

1915 

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

1917 '''I{Approximate} the perimeter of a path, polygon. or composite. 

1918 

1919 @arg points: The path or polygon points or clips (C{LatLon}[], 

1920 L{BooleanFHP} or L{BooleanGH}). 

1921 @kwarg closed: Optionally, close the path or polygon (C{bool}). 

1922 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta 

1923 by the cosine of the mean latitude (C{bool}). 

1924 @kwarg radius: Mean earth radius (C{meter}). 

1925 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

1926 B{C{points}} (C{bool}). 

1927 

1928 @return: Approximate perimeter (C{meter}, same units as 

1929 B{C{radius}}). 

1930 

1931 @raise PointsError: Insufficient number of B{C{points}} 

1932 

1933 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

1934 

1935 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with 

1936 C{B{points}} a composite. 

1937 

1938 @note: This perimeter is based on the L{pygeodesy.equirectangular_} 

1939 distance approximation and is ill-suited for regions exceeding 

1940 several hundred Km or Miles or with near-polar latitudes. 

1941 

1942 @see: Functions L{sphericalTrigonometry.perimeterOf} and 

1943 L{ellipsoidalKarney.perimeterOf}. 

1944 ''' 

1945 def _degs(ps, c, a, w): # angular edge lengths in degrees 

1946 Ps = LatLon2PsxyIter(ps, loop=1) # wrap=w 

1947 p1, u = Ps[0], _0_0 # previous x2's unroll/wrap 

1948 for p2 in Ps.iterate(closed=c): 

1949 if w and c: 

1950 w = not Ps.looped 

1951 # apply previous x2's unroll/wrap'd to new x1 

1952 _, dy, dx, u = equirectangular_(p1.y, p1.x + u, 

1953 p2.y, p2.x, 

1954 adjust=a, limit=None, 

1955 wrap=w) # PYCHOK non-seq 

1956 yield hypot(dx, dy) 

1957 p1 = p2 

1958 

1959 if _MODS.booleans.isBoolean(points): 

1960 if not closed: 

1961 notImplemented(None, closed=closed, points=_composite_) 

1962 d = points._sum1(perimeterOf, closed=True, adjust=adjust, 

1963 radius=radius, wrap=wrap) 

1964 else: 

1965 d = fsum(_degs(points, closed, adjust, wrap), floats=True) 

1966 return degrees2m(d, radius=radius) 

1967 

1968 

1969def quadOf(latS, lonW, latN, lonE, closed=False, LatLon=LatLon_, **LatLon_kwds): 

1970 '''Generate a quadrilateral path or polygon from two points. 

1971 

1972 @arg latS: Souther-nmost latitude (C{degrees90}). 

1973 @arg lonW: Western-most longitude (C{degrees180}). 

1974 @arg latN: Norther-nmost latitude (C{degrees90}). 

1975 @arg lonE: Eastern-most longitude (C{degrees180}). 

1976 @kwarg closed: Optionally, close the path (C{bool}). 

1977 @kwarg LatLon: Class to use (L{LatLon_}). 

1978 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} 

1979 keyword arguments. 

1980 

1981 @return: Return a tuple of 4 or 5 B{C{LatLon}} instances 

1982 outlining the quadrilateral. 

1983 

1984 @see: Function L{boundsOf}. 

1985 ''' 

1986 t = (LatLon(latS, lonW, **LatLon_kwds), 

1987 LatLon(latN, lonW, **LatLon_kwds), 

1988 LatLon(latN, lonE, **LatLon_kwds), 

1989 LatLon(latS, lonE, **LatLon_kwds)) 

1990 if closed: 

1991 t += t[:1] 

1992 return t 

1993 

1994 

1995__all__ += _ALL_DOCS(_Array2LatLon, _Basequence) 

1996 

1997# **) MIT License 

1998# 

1999# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

2000# 

2001# Permission is hereby granted, free of charge, to any person obtaining a 

2002# copy of this software and associated documentation files (the "Software"), 

2003# to deal in the Software without restriction, including without limitation 

2004# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

2005# and/or sell copies of the Software, and to permit persons to whom the 

2006# Software is furnished to do so, subject to the following conditions: 

2007# 

2008# The above copyright notice and this permission notice shall be included 

2009# in all copies or substantial portions of the Software. 

2010# 

2011# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

2012# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

2013# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

2014# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

2015# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

2016# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

2017# OTHER DEALINGS IN THE SOFTWARE.