Coverage for pygeodesy/geodesicx/gxarea.py: 95%

213 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-25 12:04 -0400

1# -*- coding: utf-8 -*- 

2 

3u'''Slightly enhanced versions of classes U{PolygonArea 

4<https://GeographicLib.SourceForge.io/1.52/python/code.html# 

5module-geographiclib.polygonarea>} and C{Accumulator} from 

6I{Karney}'s Python U{geographiclib 

7<https://GeographicLib.SourceForge.io/1.52/python/index.html>}. 

8 

9Class L{GeodesicAreaExact} is intended to work with instances 

10of class L{GeodesicExact} and of I{wrapped} class C{Geodesic}, 

11see module L{pygeodesy.karney}. 

12 

13Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023) 

14and licensed under the MIT/X11 License. For more information, see the 

15U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

16''' 

17# make sure int/int division yields float quotient 

18from __future__ import division as _; del _ # PYCHOK semicolon 

19 

20from pygeodesy.basics import isodd, unsigned0 

21from pygeodesy.constants import NAN, _0_0, _0_5, _720_0 

22# from pygeodesy.interns import _COMMASPACE_ # from .lazily 

23from pygeodesy.karney import Area3Tuple, _diff182, GeodesicError, \ 

24 _norm180, _remainder, _sum2_ 

25from pygeodesy.lazily import _ALL_DOCS, printf, _COMMASPACE_ 

26from pygeodesy.named import ADict, callername, _NamedBase, pairs 

27from pygeodesy.props import Property, Property_RO, property_RO 

28# from pygeodesy.streprs import pairs # from .named 

29 

30from math import fmod as _fmod 

31 

32__all__ = () 

33__version__ = '24.05.19' 

34 

35 

36class GeodesicAreaExact(_NamedBase): 

37 '''Area and perimeter of a geodesic polygon, an enhanced 

38 version of I{Karney}'s Python class U{PolygonArea 

39 <https://GeographicLib.SourceForge.io/html/python/ 

40 code.html#module-geographiclib.polygonarea>} using 

41 the more accurate surface area. 

42 

43 @note: The name of this class C{*Exact} is a misnomer, see 

44 I{Karney}'s comments at C++ attribute U{GeodesicExact._c2 

45 <https://GeographicLib.SourceForge.io/C++/doc/ 

46 GeodesicExact_8cpp_source.html>}. 

47 ''' 

48 _Area = None 

49 _g_gX = None # Exact or not 

50 _lat0 = _lon0 = \ 

51 _lat1 = _lon1 = NAN 

52 _mask = 0 

53 _num = 0 

54 _Peri = None 

55 _verbose = False 

56 _xings = 0 

57 

58 def __init__(self, geodesic, polyline=False, **name): 

59 '''New L{GeodesicAreaExact} instance. 

60 

61 @arg geodesic: A geodesic (L{GeodesicExact}, I{wrapped} 

62 C{Geodesic} or L{GeodesicSolve}). 

63 @kwarg polyline: If C{True}, compute the perimeter only, 

64 otherwise area and perimeter (C{bool}). 

65 @kwarg name: Optional C{B{name}=NN} (C{str}). 

66 

67 @raise GeodesicError: Invalid B{C{geodesic}}. 

68 ''' 

69 try: # results returned as L{GDict} 

70 if not (callable(geodesic._GDictDirect) and 

71 callable(geodesic._GDictInverse)): 

72 raise TypeError() 

73 except (AttributeError, TypeError): 

74 raise GeodesicError(geodesic=geodesic) 

75 

76 self._g_gX = g = geodesic 

77 # use the class-level Caps since the values 

78 # differ between GeodesicExact and Geodesic 

79 self._mask = g.DISTANCE | g.LATITUDE | g.LONGITUDE 

80 self._Peri = _Accumulator(name='_Peri') 

81 if not polyline: # perimeter and area 

82 self._mask |= g.AREA | g.LONG_UNROLL 

83 self._Area = _Accumulator(name='_Area') 

84 if g.debug: # PYCHOK no cover 

85 self.verbose = True # debug as verbosity 

86 if name: 

87 self.name = name 

88 

89 def AddEdge(self, azi, s): 

90 '''Add another polygon edge. 

91 

92 @arg azi: Azimuth at the current point (compass 

93 C{degrees360}). 

94 @arg s: Length of the edge (C{meter}). 

95 ''' 

96 if self.num < 1: 

97 raise GeodesicError(num=self.num) 

98 r = self._Direct(azi, s) 

99 p = self._Peri.Add(s) 

100 if self._Area: 

101 a = self._Area.Add(r.S12) 

102 self._xings += r.xing 

103 else: 

104 a = NAN 

105 self._lat1 = r.lat2 

106 self._lon1 = r.lon2 

107 self._num += 1 

108 if self.verbose: # PYCHOK no cover 

109 self._print(self.num, p, a, r, lat1=r.lat2, lon1=r.lon2, 

110 azi=azi, s=s) 

111 return self.num 

112 

113 def AddPoint(self, lat, lon): 

114 '''Add another polygon point. 

115 

116 @arg lat: Latitude of the point (C{degrees}). 

117 @arg lon: Longitude of the point (C{degrees}). 

118 ''' 

119 if self.num > 0: 

120 r = self._Inverse(self.lat1, self.lon1, lat, lon) 

121 s = r.s12 

122 p = self._Peri.Add(s) 

123 if self._Area: 

124 a = self._Area.Add(r.S12) 

125 self._xings += r.xing 

126 else: 

127 a = NAN 

128 else: 

129 self._lat0 = lat 

130 self._lon0 = lon 

131 a = p = s = _0_0 

132 r = None 

133 self._lat1 = lat 

134 self._lon1 = lon 

135 self._num += 1 

136 if self.verbose: # PYCHOK no cover 

137 self._print(self.num, p, a, r, lat1=lat, lon1=lon, s=s) 

138 return self.num 

139 

140 @Property_RO 

141 def area0x(self): 

142 '''Get the ellipsoid's surface area (C{meter} I{squared}), 

143 more accurate for very I{oblate} ellipsoids. 

144 ''' 

145 return self.ellipsoid.areax # not .area! 

146 

147 area0 = area0x # for C{geographiclib} compatibility 

148 

149 def Compute(self, reverse=False, sign=True): 

150 '''Compute the accumulated perimeter and area. 

151 

152 @kwarg reverse: If C{True}, clockwise traversal counts as a 

153 positive area instead of counter-clockwise 

154 (C{bool}). 

155 @kwarg sign: If C{True}, return a signed result for the area if 

156 the polygon is traversed in the "wrong" direction 

157 instead of returning the area for the rest of the 

158 earth. 

159 

160 @return: L{Area3Tuple}C{(number, perimeter, area)} with the 

161 number of points, the perimeter in C{meter} and the 

162 area in C{meter**2}. The perimeter includes the 

163 length of a final edge, connecting the current to 

164 the initial point, if this polygon was initialized 

165 with C{polyline=False}. For perimeter only, i.e. 

166 C{polyline=True}, area is C{NAN}. 

167 

168 @note: Arbitrarily complex polygons are allowed. In the case 

169 of self-intersecting polygons, the area is accumulated 

170 "algebraically". E.g., the areas of the 2 loops in a 

171 I{figure-8} polygon will partially cancel. 

172 

173 @note: More points and edges can be added after this call. 

174 ''' 

175 r, n = None, self.num 

176 if n < 2: 

177 p = _0_0 

178 a = NAN if self.polyline else p 

179 elif self._Area: 

180 r = self._Inverse(self.lat1, self.lon1, self.lat0, self.lon0) 

181 a = self._reduced(r.S12, reverse, sign, r.xing) 

182 p = self._Peri.Sum(r.s12) 

183 else: 

184 p = self._Peri.Sum() 

185 a = NAN 

186 if self.verbose: # PYCHOK no cover 

187 self._print(n, p, a, r, lat0=self.lat0, lon0=self.lon0) 

188 return Area3Tuple(n, p, a) 

189 

190 def _Direct(self, azi, s): 

191 '''(INTERNAL) Edge helper. 

192 ''' 

193 lon1 = self.lon1 

194 r = self._g_gX._GDictDirect(self.lat1, lon1, azi, False, s, self._mask) 

195 if self._Area: # aka transitDirect 

196 # Count crossings of prime meridian exactly as 

197 # int(ceil(lon2 / 360)) - int(ceil(lon1 / 360)) 

198 # Since we only need the parity of the result we 

199 # can use std::remquo but this is buggy with g++ 

200 # 4.8.3 and requires C++11. So instead we do: 

201 lon1 = _fmod( lon1, _720_0) # r.lon1 

202 lon2 = _fmod(r.lon2, _720_0) 

203 # int(True) == 1, int(False) == 0 

204 r.set_(xing=int(lon2 > 360 or -360 < lon2 <= 0) - 

205 int(lon1 > 360 or -360 < lon1 <= 0)) 

206 return r 

207 

208 @Property_RO 

209 def ellipsoid(self): 

210 '''Get this area's ellipsoid (C{Ellipsoid[2]}). 

211 ''' 

212 return self._g_gX.ellipsoid 

213 

214 @Property_RO 

215 def geodesic(self): 

216 '''Get this area's geodesic object (C{Geodesic[Exact]}). 

217 ''' 

218 return self._g_gX 

219 

220 earth = geodesic # for C{geographiclib} compatibility 

221 

222 def _Inverse(self, lat1, lon1, lat2, lon2): 

223 '''(INTERNAL) Point helper. 

224 ''' 

225 r = self._g_gX._GDictInverse(lat1, lon1, lat2, lon2, self._mask) 

226 if self._Area: # aka transit 

227 # count crossings of prime meridian as +1 or -1 

228 # if in east or west direction, otherwise 0 

229 lon1 = _norm180(lon1) 

230 lon2 = _norm180(lon2) 

231 lon12, _ = _diff182(lon1, lon2) 

232 r.set_(xing=int(lon12 > 0 and lon1 <= 0 and lon2 > 0) or 

233 -int(lon12 < 0 and lon2 <= 0 and lon1 > 0)) 

234 return r 

235 

236 @property_RO 

237 def lat0(self): 

238 '''Get the first point's latitude (C{degrees}). 

239 ''' 

240 return self._lat0 

241 

242 @property_RO 

243 def lat1(self): 

244 '''Get the most recent point's latitude (C{degrees}). 

245 ''' 

246 return self._lat1 

247 

248 @property_RO 

249 def lon0(self): 

250 '''Get the first point's longitude (C{degrees}). 

251 ''' 

252 return self._lon0 

253 

254 @property_RO 

255 def lon1(self): 

256 '''Get the most recent point's longitude (C{degrees}). 

257 ''' 

258 return self._lon1 

259 

260 @property_RO 

261 def num(self): 

262 '''Get the current number of points (C{int}). 

263 ''' 

264 return self._num 

265 

266 @Property_RO 

267 def polyline(self): 

268 '''Is this perimeter only (C{bool}), area NAN? 

269 ''' 

270 return self._Area is None 

271 

272 def _print(self, n, p, a, r, **kwds): # PYCHOK no cover 

273 '''(INTERNAL) Print a verbose line. 

274 ''' 

275 d = ADict(p=p, s12=r.s12 if r else NAN, **kwds) 

276 if self._Area: 

277 d.set_(a=a, S12=r.S12 if r else NAN) 

278 t = _COMMASPACE_.join(pairs(d, prec=10)) 

279 printf('%s %s: %s (%s)', self.named2, n, t, callername(up=2)) 

280 

281 def _reduced(self, S12, reverse, sign, xing): 

282 '''(INTERNAL) Accumulate and reduce area to allowed range. 

283 ''' 

284 a0 = self.area0x 

285 A = _Accumulator(self._Area) 

286 _ = A.Add(S12) 

287 a = A.Remainder(a0) # clockwise 

288 if isodd(self._xings + xing): 

289 a = A.Add((a0 if a < 0 else -a0) * _0_5) 

290 if not reverse: 

291 a = A.Negate() # counter-clockwise 

292 # (-area0x/2, area0x/2] if sign else [0, area0x) 

293 a0_ = a0 if sign else (a0 * _0_5) 

294 if a > a0_: 

295 a = A.Add(-a0) 

296 elif a <= -a0_: 

297 a = A.Add( a0) 

298 return unsigned0(a) 

299 

300 def Reset(self): 

301 '''Reset this polygon to empty. 

302 ''' 

303 if self._Area: 

304 self._Area.Reset() 

305 self._Peri.Reset() 

306 self._lat0 = self._lon0 = \ 

307 self._lat1 = self._lon1 = NAN 

308 self._num = self._xings = n = 0 

309 if self.verbose: # PYCHOK no cover 

310 printf('%s %s: (%s)', self.named2, n, self.Reset.__name__) 

311 return n 

312 

313 Clear = Reset 

314 

315 def TestEdge(self, azi, s, reverse=False, sign=True): 

316 '''Compute the properties for a tentative, additional edge 

317 

318 @arg azi: Azimuth at the current the point (compass C{degrees}). 

319 @arg s: Length of the edge (C{meter}). 

320 @kwarg reverse: If C{True}, clockwise traversal counts as a 

321 positive area instead of counter-clockwise 

322 (C{bool}). 

323 @kwarg sign: If C{True}, return a signed result for the area if 

324 the polygon is traversed in the "wrong" direction 

325 instead of returning the area for the rest of the 

326 earth. 

327 

328 @return: L{Area3Tuple}C{(number, perimeter, area)}. 

329 

330 @raise GeodesicError: No points. 

331 ''' 

332 n = self.num + 1 

333 p = self._Peri.Sum(s) 

334 if self.polyline: 

335 a, r = NAN, None 

336 elif n < 2: 

337 raise GeodesicError(num=self.num) 

338 else: 

339 d = self._Direct(azi, s) 

340 r = self._Inverse(d.lat2, d.lon2, self.lat0, self.lon0) 

341 a = self._reduced(d.S12 + r.S12, reverse, sign, d.xing + r.xing) 

342 p += r.s12 

343 if self.verbose: # PYCHOK no cover 

344 self._print(n, p, a, r, azi=azi, s=s) 

345 return Area3Tuple(n, p, a) 

346 

347 def TestPoint(self, lat, lon, reverse=False, sign=True): 

348 '''Compute the properties for a tentative, additional vertex 

349 

350 @arg lat: Latitude of the point (C{degrees}). 

351 @arg lon: Longitude of the point (C{degrees}). 

352 @kwarg reverse: If C{True}, clockwise traversal counts as a 

353 positive area instead of counter-clockwise 

354 (C{bool}). 

355 @kwarg sign: If C{True}, return a signed result for the area if 

356 the polygon is traversed in the "wrong" direction 

357 instead of returning the area for the rest of the 

358 earth. 

359 

360 @return: L{Area3Tuple}C{(number, perimeter, area)}. 

361 ''' 

362 r, n = None, self.num + 1 

363 if n < 2: 

364 p = _0_0 

365 a = NAN if self.polyline else p 

366 else: 

367 i = self._Inverse(self.lat1, self.lon1, lat, lon) 

368 p = self._Peri.Sum(i.s12) 

369 if self._Area: 

370 r = self._Inverse(lat, lon, self.lat0, self.lon0) 

371 a = self._reduced(i.S12 + r.S12, reverse, sign, i.xing + r.xing) 

372 p += r.s12 

373 else: 

374 a = NAN 

375 if self.verbose: # PYCHOK no cover 

376 self._print(n, p, a, r, lat=lat, lon=lon) 

377 return Area3Tuple(n, p, a) 

378 

379 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature 

380 '''Return this C{GeodesicExactArea} as string. 

381 

382 @kwarg prec: The C{float} precision, number of decimal digits (0..9). 

383 Trailing zero decimals are stripped for B{C{prec}} values 

384 of 1 and above, but kept for negative B{C{prec}} values. 

385 @kwarg sep: Separator to join (C{str}). 

386 

387 @return: Area items (C{str}). 

388 ''' 

389 n, p, a = self.Compute() 

390 d = dict(geodesic=self.geodesic, num=n, area=a, 

391 perimeter=p, polyline=self.polyline) 

392 return sep.join(pairs(d, prec=prec)) 

393 

394 @Property 

395 def verbose(self): 

396 '''Get the C{verbose} option (C{bool}). 

397 ''' 

398 return self._verbose 

399 

400 @verbose.setter # PYCHOK setter! 

401 def verbose(self, verbose): # PYCHOK no cover 

402 '''Set the C{verbose} option (C{bool}) to print 

403 a message after each method invokation. 

404 ''' 

405 self._verbose = bool(verbose) 

406 

407 

408class PolygonArea(GeodesicAreaExact): 

409 '''For C{geographiclib} compatibility, sub-class of L{GeodesicAreaExact}. 

410 ''' 

411 def __init__(self, earth, polyline=False, **name): 

412 '''New L{PolygonArea} instance. 

413 

414 @arg earth: A geodesic (L{GeodesicExact}, I{wrapped} 

415 C{Geodesic} or L{GeodesicSolve}). 

416 @kwarg polyline: If C{True} perimeter only, otherwise 

417 area and perimeter (C{bool}). 

418 @kwarg name: Optional C{B{name}=NN} (C{str}). 

419 

420 @raise GeodesicError: Invalid B{C{earth}}. 

421 ''' 

422 GeodesicAreaExact.__init__(self, earth, polyline=polyline, **name) 

423 

424 

425class _Accumulator(_NamedBase): 

426 '''Like C{math.fsum}, but allowing a running sum. 

427 

428 Original from I{Karney}'s U{geographiclib 

429 <https://PyPI.org/project/geographiclib>}C{.accumulator}, 

430 enhanced to return the current sum by most methods. 

431 ''' 

432 _n = 0 # len() 

433 _s = _t = _0_0 

434 

435 def __init__(self, y=0, **name): 

436 '''New L{_Accumulator}. 

437 

438 @kwarg y: Initial value (C{scalar}). 

439 @kwarg name: Optional C{B{name}=NN} (C{str}). 

440 ''' 

441 if isinstance(y, _Accumulator): 

442 self._s, self._t, self._n = y._s, y._t, 1 

443 elif y: 

444 self._s, self._n = float(y), 1 

445 if name: 

446 self.name = name 

447 

448 def Add(self, y): 

449 '''Add a value. 

450 

451 @return: Current C{sum}. 

452 ''' 

453 self._n += 1 

454 self._s, self._t = _sum2_(self._s, self._t, y) 

455 return self._s # current .Sum() 

456 

457 def Negate(self): 

458 '''Negate sum. 

459 

460 @return: Current C{sum}. 

461 ''' 

462 self._s = s = -self._s 

463 self._t = -self._t 

464 return s # current .Sum() 

465 

466 @property_RO 

467 def num(self): 

468 '''Get the current number of C{Add}itions (C{int}). 

469 ''' 

470 return self._n 

471 

472 def Remainder(self, y): 

473 '''Remainder on division by B{C{y}}. 

474 

475 @return: Remainder of C{sum} / B{C{y}}. 

476 ''' 

477 self._s = _remainder(self._s, y) 

478# self._t = _remainder(self._t, y) 

479 self._n = -1 

480 return self.Add(_0_0) 

481 

482 def Reset(self, y=0): 

483 '''Set value from argument. 

484 ''' 

485 self._n, self._s, self._t = 0, float(y), _0_0 

486 

487 Set = Reset 

488 

489 def Sum(self, y=0): 

490 '''Return C{sum + B{y}}. 

491 

492 @note: B{C{y}} is included in the returned 

493 result, but I{not} accumulated. 

494 ''' 

495 if y: 

496 s = _Accumulator(self, name='_Sum') 

497 s.Add(y) 

498 else: 

499 s = self 

500 return s._s 

501 

502 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature 

503 '''Return this C{_Accumulator} as string. 

504 

505 @kwarg prec: The C{float} precision, number of decimal digits (0..9). 

506 Trailing zero decimals are stripped for B{C{prec}} values 

507 of 1 and above, but kept for negative B{C{prec}} values. 

508 @kwarg sep: Separator to join (C{str}). 

509 

510 @return: Accumulator (C{str}). 

511 ''' 

512 d = dict(n=self.num, s=self._s, t=self._t) 

513 return sep.join(pairs(d, prec=prec)) 

514 

515 

516__all__ += _ALL_DOCS(GeodesicAreaExact, PolygonArea) 

517 

518# **) MIT License 

519# 

520# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

521# 

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

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

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

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

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

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

528# 

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

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

531# 

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

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

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

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

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

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

538# OTHER DEALINGS IN THE SOFTWARE.