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

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

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>}.

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}.

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

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 _ # noqa: E702 ;

20from pygeodesy.basics import _copysign, isodd, unsigned0

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

22from pygeodesy.fmath import _fma

23from pygeodesy.internals import printf, typename

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

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

26 _norm180, _remainder, _sum2

27from pygeodesy.lazily import _ALL_DOCS, _COMMASPACE_

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

29from pygeodesy.props import Property, Property_RO, property_RO

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

32from math import fabs, fmod as _fmod

34__all__ = ()

35__version__ = '25.12.23'

38class GeodesicAreaExact(_NamedBase):

39 '''Area and perimeter of a geodesic polygon, an enhanced version of I{Karney}'s

40 Python class U{PolygonArea<https://GeographicLib.SourceForge.io/html/python/

41 code.html#module-geographiclib.polygonarea>} using the more accurate surface area.

43 @note: The name of this class C{*Exact} is a misnomer, see I{Karney}'s comments at

44 C++ attribute U{GeodesicExact._c2<https://GeographicLib.SourceForge.io/C++/doc/

45 GeodesicExact_8cpp_source.html>}.

46 '''

47 _Area = None

48 _g_gX = None # Exact or not

49 _lat0 = _lon0 = \

50 _lat1 = _lon1 = NAN

51 _mask = 0

52 _num = 0

53 _Peri = None

54 _verbose = False

55 _xings = 0

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

58 '''New L{GeodesicAreaExact} instance.

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

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

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

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

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

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

67 '''

68 try: # results returned as L{GDict}

69 if not (callable(geodesic._GDictDirect) and

70 callable(geodesic._GDictInverse)):

71 raise TypeError()

72 except (AttributeError, TypeError):

73 raise GeodesicError(geodesic=geodesic)

75 self._g_gX = g = geodesic

76 # use the class-level Caps since the values

77 # differ between GeodesicExact and Geodesic

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

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

80 if not polyline: # perimeter and area

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

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

83 if g.debug: # PYCHOK no cover

84 self.verbose = True # debug as verbosity

85 if name:

86 self.name = name

88 def AddEdge(self, azi, s):

89 '''Add another polygon edge.

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

92 C{degrees360}).

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

94 '''

95 if self.num < 1:

96 raise GeodesicError(num=self.num)

97 r = self._Direct(azi, s)

98 p = self._Peri.Add(s)

99 if self._Area:

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

101 self._xings += r.xing

102 else:

103 a = NAN

104 self._lat1 = r.lat2

105 self._lon1 = r.lon2

106 self._num += 1

107 if self.verbose: # PYCHOK no cover

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

109 azi=azi, s=s)

110 return self.num

111

112 def AddPoint(self, lat, lon):

113 '''Add another polygon point.

114

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

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

117 '''

118 if self.num > 0:

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

120 s = r.s12

121 p = self._Peri.Add(s)

122 if self._Area:

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

124 self._xings += r.xing

125 else:

126 a = NAN

127 else:

128 self._lat0 = lat

129 self._lon0 = lon

130 a = p = s = _0_0

131 r = None

132 self._lat1 = lat

133 self._lon1 = lon

134 self._num += 1

135 if self.verbose: # PYCHOK no cover

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

137 return self.num

138

139 @Property_RO

140 def area0x(self):

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

142 for very I{oblate} ellipsoids.

143 '''

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

145

146 area0 = area0x # for C{geographiclib} compatibility

147

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

149 '''Compute the accumulated perimeter and area.

150

151 @kwarg reverse: If C{True}, clockwise traversal counts as a positive area instead

152 of counter-clockwise (C{bool}).

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

154 traversed in the "wrong" direction instead of returning the area for

155 the rest of the earth.

156 @kwarg polar: Use C{B{polar}=True} if the polygon encloses a pole (C{bool}), see

157 function L{ispolar<pygeodesy.points.ispolar>} and U{area of a polygon

158 enclosing a pole<https://GeographicLib.SourceForge.io/C++/doc/

159 classGeographicLib_1_1GeodesicExact.html#a3d7a9155e838a09a48dc14d0c3fac525>}.

160

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

162 perimeter in C{meter} and the (signed) area in C{meter**2}. The perimeter

163 includes the length of a final edge, connecting the current to the initial

164 point, if this polygon was initialized with C{polyline=False}. For perimeter

165 only, i.e. C{polyline=True}, area is C{NAN}.

166

167 @note: Arbitrarily complex polygons are allowed. In the case of self-intersecting

168 polygons, the area is accumulated "algebraically". E.g., the areas of both

169 loops in a I{figure-8} polygon will partially cancel.

170

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

172 '''

173 r, n = None, self.num

174 if n < 2:

175 p = _0_0

176 a = NAN if n > 0 and self.polyline else p

177 elif self._Area:

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

179 a = self._reduced(r.S12, r.xing, n, reverse=reverse, sign=sign, polar=polar)

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

181 else:

182 p = self._Peri.Sum()

183 a = NAN

184 if self.verbose: # PYCHOK no cover

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

186 return Area3Tuple(n, p, a)

187

188 def _Direct(self, azi, s):

189 '''(INTERNAL) Edge helper.

190 '''

191 lon1 = self.lon1

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

193 if self._Area: # aka transitDirect

194 # Count crossings of prime meridian exactly as

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

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

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

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

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

200 lon2 = _fmod(r.lon2, _720_0)

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

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

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

204 return r

205

206 @Property_RO

207 def ellipsoid(self):

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

209 '''

210 return self._g_gX.ellipsoid

211

212 @Property_RO

213 def geodesic(self):

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

215 '''

216 return self._g_gX

217

218 earth = geodesic # for C{geographiclib} compatibility

219

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

221 '''(INTERNAL) Point helper.

222 '''

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

224 if self._Area: # aka transit

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

226 # if in east or west direction, otherwise 0

227 lon1 = _norm180(lon1)

228 lon2 = _norm180(lon2)

229 lon12, _ = _diff182(lon1, lon2)

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

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

232 return r

233

234 @property_RO

235 def lat0(self):

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

237 '''

238 return self._lat0

239

240 @property_RO

241 def lat1(self):

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

243 '''

244 return self._lat1

245

246 @property_RO

247 def lon0(self):

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

249 '''

250 return self._lon0

251

252 @property_RO

253 def lon1(self):

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

255 '''

256 return self._lon1

257

258 @property_RO

259 def num(self):

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

261 '''

262 return self._num

263

264 @Property_RO

265 def polyline(self):

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

267 '''

268 return self._Area is None

269

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

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

272 '''

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

274 if self._Area:

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

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

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

278

279 def _reduced(self, S12, xing, n, reverse=False, sign=True, polar=False):

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

281 '''

282 a0 = self.area0x

283 A = _Accumulator(self._Area)

284 _ = A.Add(S12)

285 a = A.Remainder(a0) # clockwise

286 if isodd(self._xings + xing):

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

288 if not reverse:

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

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

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

292 if a > a0_:

293 a = A.Add(-a0)

294 elif a <= -a0_:

295 a = A.Add( a0)

296 if polar: # see .geodesicw._gwrapped.Geodesic.Area

297 a = A.Add(_copysign(a0 * _0_5 * n, a)) # - if reverse or sign?

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, typename(self.Reset))

311 return n

312

313 Clear = Reset

314

315 def TestEdge(self, azi, s, **reverse_sign_polar):

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_sign_polar: Optional C{B{reverse}=False}, C{B{sign}=True} and

321 C{B{polar}=False} keyword arguments, see method L{Compute}.

322

323 @return: L{Area3Tuple}C{(number, perimeter, area)}, with C{perimeter} and

324 C{area} both C{NAN} for insuffcient C{number} of points.

325 '''

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

327 if n < 2: # raise GeodesicError(num=self.num)

328 a = p = NAN # like .test_Planimeter19

329 else:

330 p = self._Peri.Sum(s)

331 if self.polyline:

332 a = NAN

333 else:

334 d = self._Direct(azi, s)

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

336 a = self._reduced(d.S12 + r.S12, d.xing + r.xing, n, **reverse_sign_polar)

337 p += r.s12

338 if self.verbose: # PYCHOK no cover

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

340 return Area3Tuple(n, p, a)

341

342 def TestPoint(self, lat, lon, **reverse_sign_polar):

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

344

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

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

347 @kwarg reverse_sign_polar: Optional C{B{reverse}=False}, C{B{sign}=True} and

348 C{B{polar}=False} keyword arguments, see method L{Compute}.

349

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

351 '''

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

353 if n < 2:

354 p = _0_0

355 a = NAN if self.polyline else p

356 else:

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

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

359 if self._Area:

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

361 a = self._reduced(i.S12 + r.S12, i.xing + r.xing, n, **reverse_sign_polar)

362 p += r.s12

363 else:

364 a = NAN

365 if self.verbose: # PYCHOK no cover

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

367 return Area3Tuple(n, p, a)

368

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

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

371

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

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

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

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

376

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

378 '''

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

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

381 perimeter=p, polyline=self.polyline)

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

383

384 @Property

385 def verbose(self):

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

387 '''

388 return self._verbose

389

390 @verbose.setter # PYCHOK setter!

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

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

393 a message after each method invokation.

394 '''

395 self._verbose = bool(verbose)

396

397

398class PolygonArea(GeodesicAreaExact):

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

400 '''

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

402 '''New L{PolygonArea} instance.

403

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

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

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

407 perimeter and area (C{bool}).

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

409

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

411 '''

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

413

414

415class _Accumulator(_NamedBase):

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

417

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

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

420 enhanced to return the current sum by most methods.

421 '''

422 _n = 0 # len()

423 _s = _t = _0_0

424

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

426 '''New L{_Accumulator}.

427

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

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

430 '''

431 self._s_t(*_s_t2(y))

432 self._n = 1

433 if name:

434 self.name = name

435

436 def Add(self, *ys):

437 '''Add one or more scalars or L{_Accumulator}s.

438

439 @return: Current C{sum}.

440

441 @see: C++ U{Accumulator.Add<https://GeographicLib.sourceforge.io/C++/doc/Accumulator_8hpp_source.html>}

442 for more details about Karney's and Shewchuk's addition.

443 '''

444 # _Accumulator().Add(1, 1e20, 2, 100, 5000, -1e20) ... 5103.0

445 s, t = self._s, self._t

446 for y in ys:

447 for y in _s_t2(y):

448 if y:

449 t, u = _sum2(t, y)

450 s, t = _sum2(s, t)

451 if s: # accumlate u in t

452 t += u

453 else: # s == 0 implies t == 0

454 s = u

455 self._n += 1

456 return self._s_t(s, t)

457

458 def Negate(self):

459 '''Negate sum.

460

461 @return: Current C{sum}.

462 '''

463 return self._s_t(-self._s, -self._t)

464

465 @property_RO

466 def num(self):

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

468 '''

469 return self._n

470

471 def Remainder(self, y):

472 '''Remainder of division by B{C{y}}.

473

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

475 '''

476 return self._s_t(_remainder(self._s, y),

477 _remainder(self._t, y))

478

479 def Reset(self, y=0):

480 '''Reset from scalar or L{_Accumulator}.

481 '''

482 self._s_t(*_s_t2(y))

483 self._n = 0

484

485 Set = Reset

486

487 def _s_t(self, s, t=0):

488 if t and fabs(s) < fabs(t):

489 s, t = t, s

490 self._s, self._t = s, t

491 return s

492

493 def Sum(self, y=0):

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

495

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

497 result, but I{not} accumulated.

498 '''

499 if y:

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

501 s.Add(y)

502 else:

503 s = self

504 return s._s

505

506 def Times(self, y):

507 '''Multiply by a scalar.

508

509 @return: Current C{sum}.

510 '''

511 s = d = self._s

512 s *= y

513 d = _fma(y, d, -s)

514 t = _fma(y, self._t, d)

515 return self._s_t(s, t) # current .Sum()

516

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

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

519

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

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

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

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

524

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

526 '''

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

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

529

530

531def _s_t2(y):

532 return (y._s, y._t) if isinstance(y, _Accumulator) else (float(y),) # PYCHOK OK

533

534

535__all__ += _ALL_DOCS(GeodesicAreaExact, PolygonArea)

536

537# **) MIT License

538#

540#

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

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

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

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

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

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

547#

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

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

550#

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

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

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

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

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

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

557# OTHER DEALINGS IN THE SOFTWARE.