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

2

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

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

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

6 I{Karney}'s Python U{geographiclib

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

8

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

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

11 see module L{pygeodesy.karney}.

12

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

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

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

16'''

17# make sure int/int division yields float quotient

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

19

20 from pygeodesy.basics import isodd, unsigned0

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

22 from pygeodesy.internals import printf

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

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

25 _norm180, _remainder, _sum3

26 from pygeodesy.lazily import _ALL_DOCS, _COMMASPACE_

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

28 from pygeodesy.props import Property, Property_RO, property_RO

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

30

31 from math import fmod as _fmod

32

33__all__ = ()

34__version__ = '24.10.14'

35

36

37 class GeodesicAreaExact(_NamedBase):

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

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

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

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

42 the more accurate surface area.

43

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

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

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

47 GeodesicExact_8cpp_source.html>}.

48 '''

49 _Area = None

50 _g_gX = None # Exact or not

51 _lat0 = _lon0 = \

52 _lat1 = _lon1 = NAN

53 _mask = 0

54 _num = 0

55 _Peri = None

56 _verbose = False

57 _xings = 0

58

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

60 '''New L{GeodesicAreaExact} instance.

61

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

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

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

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

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

67

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

69 '''

70 try: # results returned as L{GDict}

71 if not (callable(geodesic._GDictDirect) and

72 callable(geodesic._GDictInverse)):

73 raise TypeError()

74 except (AttributeError, TypeError):

75 raise GeodesicError(geodesic=geodesic)

76

77 self._g_gX = g = geodesic

78 # use the class-level Caps since the values

79 # differ between GeodesicExact and Geodesic

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

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

82 if not polyline: # perimeter and area

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

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

85 if g.debug: # PYCHOK no cover

86 self.verbose = True # debug as verbosity

87 if name:

88 self.name = name

89

90 def AddEdge(self, azi, s):

91 '''Add another polygon edge.

92

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

94 C{degrees360}).

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

96 '''

97 if self.num < 1:

98 raise GeodesicError(num=self.num)

99 r = self._Direct(azi, s)

100 p = self._Peri.Add(s)

101 if self._Area:

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

103 self._xings += r.xing

104 else:

105 a = NAN

106 self._lat1 = r.lat2

107 self._lon1 = r.lon2

108 self._num += 1

109 if self.verbose: # PYCHOK no cover

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

111 azi=azi, s=s)

112 return self.num

113

114 def AddPoint(self, lat, lon):

115 '''Add another polygon point.

116

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

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

119 '''

120 if self.num > 0:

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

122 s = r.s12

123 p = self._Peri.Add(s)

124 if self._Area:

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

126 self._xings += r.xing

127 else:

128 a = NAN

129 else:

130 self._lat0 = lat

131 self._lon0 = lon

132 a = p = s = _0_0

133 r = None

134 self._lat1 = lat

135 self._lon1 = lon

136 self._num += 1

137 if self.verbose: # PYCHOK no cover

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

139 return self.num

140

141 @Property_RO

142 def area0x(self):

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

144 more accurate for very I{oblate} ellipsoids.

145 '''

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

147

148 area0 = area0x # for C{geographiclib} compatibility

149

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

151 '''Compute the accumulated perimeter and area.

152

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

154 positive area instead of counter-clockwise

155 (C{bool}).

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

157 the polygon is traversed in the "wrong" direction

158 instead of returning the area for the rest of the

159 earth.

160

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

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

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

164 length of a final edge, connecting the current to

165 the initial point, if this polygon was initialized

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

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

168

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

170 of self-intersecting polygons, the area is accumulated

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

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

173

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

175 '''

176 r, n = None, self.num

177 if n < 2:

178 p = _0_0

179 a = NAN if self.polyline else p

180 elif self._Area:

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

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

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

184 else:

185 p = self._Peri.Sum()

186 a = NAN

187 if self.verbose: # PYCHOK no cover

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

189 return Area3Tuple(n, p, a)

190

191 def _Direct(self, azi, s):

192 '''(INTERNAL) Edge helper.

193 '''

194 lon1 = self.lon1

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

196 if self._Area: # aka transitDirect

197 # Count crossings of prime meridian exactly as

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

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

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

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

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

203 lon2 = _fmod(r.lon2, _720_0)

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

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

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

207 return r

208

209 @Property_RO

210 def ellipsoid(self):

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

212 '''

213 return self._g_gX.ellipsoid

214

215 @Property_RO

216 def geodesic(self):

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

218 '''

219 return self._g_gX

220

221 earth = geodesic # for C{geographiclib} compatibility

222

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

224 '''(INTERNAL) Point helper.

225 '''

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

227 if self._Area: # aka transit

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

229 # if in east or west direction, otherwise 0

230 lon1 = _norm180(lon1)

231 lon2 = _norm180(lon2)

232 lon12, _ = _diff182(lon1, lon2)

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

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

235 return r

236

237 @property_RO

238 def lat0(self):

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

240 '''

241 return self._lat0

242

243 @property_RO

244 def lat1(self):

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

246 '''

247 return self._lat1

248

249 @property_RO

250 def lon0(self):

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

252 '''

253 return self._lon0

254

255 @property_RO

256 def lon1(self):

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

258 '''

259 return self._lon1

260

261 @property_RO

262 def num(self):

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

264 '''

265 return self._num

266

267 @Property_RO

268 def polyline(self):

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

270 '''

271 return self._Area is None

272

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

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

275 '''

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

277 if self._Area:

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

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

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

281

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

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

284 '''

285 a0 = self.area0x

286 A = _Accumulator(self._Area)

287 _ = A.Add(S12)

288 a = A.Remainder(a0) # clockwise

289 if isodd(self._xings + xing):

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

291 if not reverse:

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

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

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

295 if a > a0_:

296 a = A.Add(-a0)

297 elif a <= -a0_:

298 a = A.Add( a0)

299 return unsigned0(a)

300

301 def Reset(self):

302 '''Reset this polygon to empty.

303 '''

304 if self._Area:

305 self._Area.Reset()

306 self._Peri.Reset()

307 self._lat0 = self._lon0 = \

308 self._lat1 = self._lon1 = NAN

309 self._num = self._xings = n = 0

310 if self.verbose: # PYCHOK no cover

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

312 return n

313

314 Clear = Reset

315

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

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

318

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

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

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

322 positive area instead of counter-clockwise

323 (C{bool}).

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

325 the polygon is traversed in the "wrong" direction

326 instead of returning the area for the rest of the

327 earth.

328

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

330

331 @raise GeodesicError: No points.

332 '''

333 n = self.num + 1

334 p = self._Peri.Sum(s)

335 if self.polyline:

336 a, r = NAN, None

337 elif n < 2:

338 raise GeodesicError(num=self.num)

339 else:

340 d = self._Direct(azi, s)

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

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

343 p += r.s12

344 if self.verbose: # PYCHOK no cover

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

346 return Area3Tuple(n, p, a)

347

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

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

350

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

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

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

354 positive area instead of counter-clockwise

355 (C{bool}).

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

357 the polygon is traversed in the "wrong" direction

358 instead of returning the area for the rest of the

359 earth.

360

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

362 '''

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

364 if n < 2:

365 p = _0_0

366 a = NAN if self.polyline else p

367 else:

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

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

370 if self._Area:

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

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

373 p += r.s12

374 else:

375 a = NAN

376 if self.verbose: # PYCHOK no cover

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

378 return Area3Tuple(n, p, a)

379

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

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

382

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

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

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

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

387

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

389 '''

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

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

392 perimeter=p, polyline=self.polyline)

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

394

395 @Property

396 def verbose(self):

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

398 '''

399 return self._verbose

400

401 @verbose.setter # PYCHOK setter!

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

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

404 a message after each method invokation.

405 '''

406 self._verbose = bool(verbose)

407

408

409 class PolygonArea(GeodesicAreaExact):

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

411 '''

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

413 '''New L{PolygonArea} instance.

414

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

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

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

418 perimeter and area (C{bool}).

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

420

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

422 '''

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

424

425

426 class _Accumulator(_NamedBase):

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

428

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

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

431 enhanced to return the current sum by most methods.

432 '''

433 _n = 0 # len()

434 _s = _t = _0_0

435

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

437 '''New L{_Accumulator}.

438

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

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

441 '''

442 if isinstance(y, _Accumulator):

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

444 elif y:

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

446 if name:

447 self.name = name

448

449 def Add(self, y):

450 '''Add a value.

451

452 @return: Current C{sum}.

453 '''

454 self._n += 1

455 self._s, self._t, _ = _sum3(self._s, self._t, y)

456 return self._s # current .Sum()

457

458 def Negate(self):

459 '''Negate sum.

460

461 @return: Current C{sum}.

462 '''

463 self._s = s = -self._s

464 self._t = -self._t

465 return s # current .Sum()

466

467 @property_RO

468 def num(self):

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

470 '''

471 return self._n

472

473 def Remainder(self, y):

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

475

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

477 '''

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

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

480 self._n = -1

481 return self.Add(_0_0)

482

483 def Reset(self, y=0):

484 '''Set value from argument.

485 '''

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

487

488 Set = Reset

489

490 def Sum(self, y=0):

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

492

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

494 result, but I{not} accumulated.

495 '''

496 if y:

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

498 s.Add(y)

499 else:

500 s = self

501 return s._s

502

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

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

505

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

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

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

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

510

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

512 '''

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

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

515

516

517__all__ += _ALL_DOCS(GeodesicAreaExact, PolygonArea)

518

519# **) MIT License

520#

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

522#

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

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

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

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

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

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

529#

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

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

532#

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

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

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

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

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

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

539# OTHER DEALINGS IN THE SOFTWARE.