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:Charles@Karney.com>} (2008-2022)

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 # from .karney

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

22 from pygeodesy.interns import NN, _COMMASPACE_

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

24 _norm180, _remainder, _sum2_, unsigned0

25 from pygeodesy.lazily import _ALL_DOCS, printf

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

27 from pygeodesy.props import Property, Property_RO, property_RO

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

29

30 from math import fmod

31

32__all__ = ()

33__version__ = '23.04.05'

34

35

36 class 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/C++/doc/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=NN):

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 name (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 = _Dict(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

408 class PolygonArea(GeodesicAreaExact):

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

410 '''

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

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 name (C{str}).

419

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

421 '''

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

423

424

425 class _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=NN):

436 '''New L{_Accumulator}.

437 '''

438 if isinstance(y, _Accumulator):

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

440 elif y:

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

442 if name:

443 self.name = name

444

445 def Add(self, y):

446 '''Add a value.

447

448 @return: Current C{sum}.

449 '''

450 self._n += 1

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

452 return self._s # current .Sum()

453

454 def Negate(self):

455 '''Negate sum.

456

457 @return: Current C{sum}.

458 '''

459 self._s = s = -self._s

460 self._t = -self._t

461 return s # current .Sum()

462

463 @property_RO

464 def num(self):

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

466 '''

467 return self._n

468

469 def Remainder(self, y):

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

471

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

473 '''

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

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

476 self._n = -1

477 return self.Add(_0_0)

478

479 def Reset(self, y=0):

480 '''Set value from argument.

481 '''

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

483

484 Set = Reset

485

486 def Sum(self, y=0):

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

488

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

490 result, but I{not} accumulated.

491 '''

492 if y:

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

494 s.Add(y)

495 else:

496 s = self

497 return s._s

498

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

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

501

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

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

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

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

506

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

508 '''

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

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

511

512

513__all__ += _ALL_DOCS(GeodesicAreaExact, PolygonArea)

514

515# **) MIT License

516#

517# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

518#

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

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

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

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

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

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

525#

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

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

528#

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

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

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

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

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

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

535# OTHER DEALINGS IN THE SOFTWARE.