Coverage for pygeodesy/karney.py: 92%

1 

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

3 

4u'''Wrapper around several C{geomath.Math} functions from I{Karney}'s Python package U{geographiclib 

5<https://PyPI.org/project/geographiclib>}, provided that package is installed. 

6 

7Methods of the I{wrapped} L{Geodesic<geodesicw.Geodesic>} and L{GeodesicLine<geodesicw.GeodesicLine>} 

8classes return a L{GDict} instance offering access to the C{dict} items either by C{key} or by 

9C{attribute} name. 

10 

11With env variable C{PYGEODESY_GEOGRAPHICLIB} left undefined or set to C{"2"}, modules L{geodesicw}, 

12L{geodesicx} and this module will use U{GeographicLib 2.0+<https://GeographicLib.SourceForge.io/C++/doc/>} 

13and newer transcoding, otherwise C{1.52} or older. Set C{PYGEODESY_GEOGRAPHICLIB=2.4} to default to the 

14C{Jacobi amplitude} instead of C{Bulirsch}' function in methods L{ExactTransverseMercator.forward 

15<pygeodesy.ExactTransverseMercator.forward>} and L{reverse <pygeodesy.ExactTransverseMercator.reverse>}. 

16 

17Karney-based functionality 

18========================== 

19 

201. The following classes and functions in C{pygeodesy} 

21 

22 - L{AlbersEqualArea}, L{AlbersEqualArea2}, L{AlbersEqualArea4}, L{AlbersEqualAreaCylindrical}, 

23 L{AlbersEqualAreaNorth}, L{AlbersEqualAreaSouth} -- U{AlbersEqualArea 

24 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AlbersEqualArea.html>} 

25 

26 - L{Ang}, L{Deg}, L{Rad} -- U{AngleT<https://GeographicLib.SourceForge.io/C++/doc/ 

27 classGeographicLib_1_1AngleT.html>} in I{GeographicLib 2.7}. 

28 

29 - L{AuxAngle}, L{AuxDST}, L{AuxDLat}, L{AuxLat} -- U{AuxAngle<https://GeographicLib.SourceForge.io/C++/doc/ 

30 classGeographicLib_1_1AuxAngle.html>}, U{DST<https://GeographicLib.SourceForge.io/C++/doc/ 

31 classGeographicLib_1_1DST.html>}, U{DAuxLatitude<https://GeographicLib.SourceForge.io/C++/doc/ 

32 classGeographicLib_1_1DAuxLatitude.html>}, U{AuxLatitude<https://GeographicLib.SourceForge.io/C++/doc/ 

33 classGeographicLib_1_1AuxLatitude.html>} in I{GeographicLib 2.2+} 

34 

35 - L{CassiniSoldner} -- U{CassiniSoldner<https://GeographicLib.SourceForge.io/C++/doc/ 

36 classGeographicLib_1_1CassiniSoldner.html>} 

37 

38 - L{Conformal3}, L{Conformal3B}, L{Conformal3Sphere} -- U{Conformal3<https://GeographicLib.SourceForge.io/C++/doc/ 

39 classGeographicLib_1_1Triaxial_1_1Conformal3.html>} in I{GeographicLib 2.7} 

40 

41 - L{EcefKarney} -- U{Geocentric<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>} 

42 

43 - L{Elliptic} -- U{EllipticFunction<https://GeographicLib.SourceForge.io/C++/doc/ 

44 classGeographicLib_1_1EllipticFunction.html>} 

45 

46 - L{EquidistantExact}, L{EquidistantGeodSolve}, L{EquidistantKarney} -- U{AzimuthalEquidistant 

47 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AzimuthalEquidistant.html>} 

48 

49 - L{Etm}, L{ExactTransverseMercator} -- U{TransverseMercatorExact<https://GeographicLib.SourceForge.io/C++/doc/ 

50 classGeographicLib_1_1TransverseMercatorExact.html>} 

51 

52 - L{Geodesic}, L{GeodesicLine} -- I{wrapped} U{geodesic.Geodesic<https://PyPI.org/project/geographiclib>}, 

53 I{wrapped} U{geodesicline.GeodesicLine<https://PyPI.org/project/geographiclib>} 

54 

55 - L{GeodesicAreaExact}, L{PolygonArea} -- U{PolygonArea<https://GeographicLib.SourceForge.io/C++/doc/ 

56 classGeographicLib_1_1PolygonAreaT.html>} 

57 

58 - L{GeodesicExact}, L{GeodesicLineExact} -- U{GeodesicExact<https://GeographicLib.SourceForge.io/C++/doc/ 

59 classGeographicLib_1_1GeodesicExact.html>}, U{GeodesicLineExact<https://GeographicLib.SourceForge.io/C++/doc/ 

60 classGeographicLib_1_1GeodesicLineExact.html>} 

61 

62 - L{GeoidKarney} -- U{Geoid<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>} 

63 

64 - L{Georef} -- U{Georef<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Georef.html>} 

65 

66 - L{GnomonicExact}, L{GnomonicGeodSolve}, L{GnomonicKarney} -- U{Gnomonic 

67 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Gnomonic.html>} 

68 

69 - L{Intersector} -- U{Intersect<https://GeographicLib.SourceForge.io/C++/doc/ 

70 classGeographicLib_1_1Intersect.html>} from I{GeographicLib 2.3+} 

71 

72 - L{JacobiConformal} -- U{JacobiConformal 

73 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1experimental_1_1JacobiConformal.html>} 

74 

75 - L{KTransverseMercator} - U{TransverseMercator<https://GeographicLib.SourceForge.io/C++/doc/ 

76 classGeographicLib_1_1TransverseMercator.html>} 

77 

78 - L{LocalCartesian}, L{Ltp} -- U{LocalCartesian<https://GeographicLib.SourceForge.io/C++/doc/ 

79 classGeographicLib_1_1LocalCartesian.html>} 

80 

81 - L{Osgr} -- U{OSGB<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1OSGB.html>} 

82 

83 - L{rhumb.aux_}, L{RhumbAux}, L{RhumbLineAux} -- U{Rhumb<https://GeographicLib.SourceForge.io/C++/doc/ 

84 classGeographicLib_1_1Rhumb.html>} and U{RhumbLine<https://GeographicLib.SourceForge.io/C++/doc/ 

85 classGeographicLib_1_1RhumbLine.html>} from I{GeographicLib 2.2+} 

86 

87 - L{rhumb.ekx}, L{Rhumb}, L{RhumbLine} -- U{Rhumb<https://GeographicLib.SourceForge.io/C++/doc/ 

88 classGeographicLib_1_1Rhumb.html>}, U{RhumbLine<https://GeographicLib.SourceForge.io/C++/doc/ 

89 classGeographicLib_1_1RhumbLine.html>}, U{TransverseMercator<https://GeographicLib.SourceForge.io/C++/doc/ 

90 classGeographicLib_1_1TransverseMercator.html>} from I{GeographicLib 2.0} 

91 

92 - L{Triaxial3}, L{Triaxial3B} -- U{Cartesian3<https://GeographicLib.SourceForge.io/C++/doc/ 

93 classGeographicLib_1_1Triaxial_1_1Cartesian3.html>}, U{Ellipsoidal3<https://GeographicLib.SourceForge.io/C++/doc/ 

94 classGeographicLib_1_1Triaxial_1_1Ellipsoidal3.html>} in I{GeographicLib 2.7} 

95 

96 - L{Ups} -- U{PolarStereographic<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1PolarStereographic.html>} 

97 

98 - L{Utm} -- U{TransverseMercator<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1TransverseMercator.html>} 

99 

100 - L{UtmUps}, L{Epsg} -- U{UTMUPS<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1UTMUPS.html>} 

101 

102 - L{atan1d}, L{atan2d}, L{sincos2}, L{sincos2d}, L{tand} -- U{geomath.Math<https://GeographicLib.SourceForge.io/C++/doc/ 

103 classGeographicLib_1_1Math.html>} 

104 

105are I{transcoded} from C++ classes in I{Karney}'s U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>}. 

106 

1072. These C{pygeodesy} modules and classes 

108 

109 - L{ellipsoidalGeodSolve}, L{ellipsoidalKarney}, L{geodesici}, L{geodsolve}, L{karney}, L{rhumb.solve} 

110 - L{EquidistantKarney}, L{FrechetKarney}, L{GeodesicSolve}, L{GeodesicLineSolve}, L{Geodesic3Solve}, 

111 L{GeodesicLine3Solve}, L{GnomonicGeodSolve}, L{GnomonicKarney}, L{HeightIDWkarney}, L{Intersectool} 

112 

113are or use I{wrappers} around I{Karney}'s Python U{geographiclib<https://PyPI.org/project/geographiclib>} or 

114C++ utility U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>}, U{Geod3Solve 

115<https://GeographicLib.SourceForge.io/C++/doc/Geod3Solve.1.html>}, U{IntersectTool 

116<https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} or U{RhumbSolve 

117<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>}. 

118 

1193. All C{pygeodesy} functions and methods to compute I{ellipsoidal} intersections, nearest points and trilaterations 

120 

121 - L{ellipsoidalExact.intersection3}, L{ellipsoidalExact.intersections2}, L{ellipsoidalExact.nearestOn}, 

122 L{ellipsoidalExact.LatLon.intersection3}, L{ellipsoidalExact.LatLon.intersections2}, 

123 L{ellipsoidalExact.LatLon.nearestOn}, L{ellipsoidalExact.LatLon.trilaterate5} 

124 

125 - L{ellipsoidalKarney.intersection3}, L{ellipsoidalKarney.intersections2}, L{ellipsoidalKarney.nearestOn}, 

126 L{ellipsoidalKarney.LatLon.intersection3}, L{ellipsoidalKarney.LatLon.intersections2}, 

127 L{ellipsoidalKarney.LatLon.nearestOn}, L{ellipsoidalKarney.LatLon.trilaterate5} 

128 

129 - L{ellipsoidalVincenty.intersection3}, L{ellipsoidalVincenty.intersections2}, L{ellipsoidalVincenty.nearestOn}, 

130 L{ellipsoidalVincenty.LatLon.intersection3}, L{ellipsoidalVincenty.LatLon.intersections2}, 

131 L{ellipsoidalVincenty.LatLon.nearestOn}, L{ellipsoidalVincenty.LatLon.trilaterate5} 

132 

133 - L{RhumbLineAux.Intersection} and L{RhumbLine.Intersection} 

134 

135are implementations of I{Karney}'s iterative solution posted under U{The B{ellipsoidal} case 

136<https://GIS.StackExchange.com/questions/48937/calculating-intersection-of-two-circles>} and in paper U{Geodesics 

137on an ellipsoid of revolution<https://ArXiv.org/pdf/1102.1215.pdf>} (pp 20-21, section B{14. MARITIME BOUNDARIES}). 

138 

1394. The C{pygeodesy} methods to compute I{ellipsoidal} intersections and nearest points 

140 

141 - L{RhumbLineAux.Intersecant2}, L{RhumbLineAux.PlumbTo}, L{RhumbLine.Intersecant2} and L{RhumbLine.PlumbTo} 

142 

143are I{transcoded} of I{Karney}'s iterative C++ function U{rhumb-intercept 

144<https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}. 

145 

1465. Spherical functions 

147 

148 - L{pygeodesy.excessKarney_}, L{sphericalTrigonometry.areaOf} 

149 

150in C{pygeodesy} are based on I{Karney}'s post U{Area of a spherical polygon 

151<https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}, 3rd Answer. 

152''' 

153# make sure int/int division yields float quotient, see .basics 

154from __future__ import division as _; del _ # noqa: E702 ; 

155 

156from pygeodesy.basics import _copysign, _isin, isint, neg, _xgeographiclib, _zip 

157from pygeodesy.constants import NAN, _isfinite as _math_isfinite, _over, \ 

158 _0_0, _1_0, _2_0, _180_0, _N_180_0, _360_0 

159# from pygeodesy.deprecated.classes import Rhumb7Tuple # _MODS 

160# from pygeodesy.ellipsoids import _Taus # _MODS 

161from pygeodesy.errors import GeodesicError, _ValueError, _xkwds 

162# from pygeodesy.geod3Solve import Geod3Solve8Tuple # _MODS 

163from pygeodesy.fmath import cbrt, fhorner, fremainder, ftruncate, norm2 

164from pygeodesy.fsums import _Ksum 

165from pygeodesy.internals import _getenv, _popen2, _PYGEODESY_ENV, typename, \ 

166 _version_info 

167from pygeodesy.interns import NN, _a12_, _area_, _azi2_, _azi12_, _composite_, \ 

168 _lat1_, _lat2_, _lon1_, _lon2_, _m12_, _M12_, \ 

169 _M21_, _number_, _s12_, _S12_, _SPACE_, \ 

170 _UNDER_, _X_, _1_, _2_, _BAR_ # PYCHOK used! 

171from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS 

172from pygeodesy.named import ADict, _NamedBase, _NamedTuple, notImplemented, _Pass 

173from pygeodesy.props import deprecated_method, property_RO, property_ROnce 

174from pygeodesy.units import Azimuth as _Azi, Degrees as _Deg, Lat, Lon, \ 

175 Meter as _M, Meter2 as _M2, Number_ 

176from pygeodesy.utily import atan2d, sincos2d, tand, _unrollon, fabs 

177 

178# from math import fabs # from .utily 

179 

180__all__ = _ALL_LAZY.karney 

181__version__ = '26.10.30' 

182 

183_1_16th = _1_0 / 16 

184_2_4_ = '2.4' 

185_K_2_0 = _getenv(_PYGEODESY_ENV(typename(_xgeographiclib)[2:]), _2_) 

186_K_2_4 = _K_2_0 == _2_4_ 

187_K_2_0 = _K_2_4 or (_K_2_0 == _2_) 

188_perimeter_ = 'perimeter' 

189 

190 

191class _GTuple(_NamedTuple): # in .testNamedTuples 

192 '''(INTERNAL) Helper. 

193 ''' 

194 def _dup_iteration(self, r): 

195 if self._iteration is not None: 

196 r._iteration = self._iteration 

197 return r 

198 

199 def toGDict(self, **updates): # NO name=NN 

200 '''Convert this C{*Tuple} to a L{GDict}. 

201 

202 @kwarg updates: Optional items to apply (C{name=value} pairs) 

203 ''' 

204 r = GDict(_zip(self._Names_, self)) # strict=True 

205 if updates: 

206 r.update(updates) 

207 return self._dup_iteration(r) 

208 

209 def _toTuple(self, nTuple, dflt, updates={}): # PYLINT W0102 

210 '''(INTERNAL) Convert this C{_GTuple} to an B{C{nTuple}}. 

211 ''' 

212 _g = self.toGDict(**updates).get 

213 t = tuple(_g(n, dflt) for n in nTuple._Names_) 

214 t = nTuple(t, name=self.name) 

215 return self._dup_iteration(t) 

216 

217 

218class _Lat(Lat): 

219 '''(INTERNAL) Latitude B{C{lat}}. 

220 ''' 

221 def __init__(self, *lat, **Error_name): 

222 kwds = _xkwds(Error_name, clip=0, Error=GeodesicError) 

223 Lat.__new__(_Lat, *lat, **kwds) 

224 

225 

226class _Lon(Lon): 

227 '''(INTERNAL) Longitude B{C{lon}}. 

228 ''' 

229 def __init__(self, *lon, **Error_name): 

230 kwds = _xkwds(Error_name, clip=0, Error=GeodesicError) 

231 Lon.__new__(_Lon, *lon, **kwds) 

232 

233 

234class Area3Tuple(_NamedTuple): # in .geodesicx.gxarea 

235 '''3-Tuple C{(number, perimeter, area)} with the C{number} 

236 of points of the polygon or polyline, the C{perimeter} in 

237 C{meter} and the C{area} in C{meter} I{squared}. 

238 ''' 

239 _Names_ = (_number_, _perimeter_, _area_) 

240 _Units_ = ( Number_, _M, _M2) 

241 

242 

243class Caps(object): 

244 '''I{Enum}-style masks to be bit-C{or}'ed to specify C{geodesic}, C{rhumb} 

245 or C{-lines} capabilities (C{caps}) and method results (C{outmask}). 

246 

247 C{AREA} - compute the area C{S12}, 

248 

249 C{AZIMUTH} - include azimuths C{azi1} and C{azi2}, 

250 

251 C{DISTANCE} - compute distance C{s12} and angular distance C{a12}, 

252 

253 C{DISTANCE_IN} - allow distance C{s12} in method C{.Direct}, 

254 

255 C{EMPTY} - nothing, formerly aka C{NONE}, 

256 

257 C{GEODESICSCALE} - compute geodesic scales C{M12} and C{M21}, 

258 

259 C{LATITUDE} - compute latitude C{lat2}, 

260 

261 C{LINE_OFF} - Line without updates from parent C{geodesic} or C{rhumb}. 

262 

263 C{LONGITUDE} - compute longitude C{lon2}, 

264 

265 C{LONG_UNROLL} - unroll C{lon2} in method C{.Direct}, 

266 

267 C{NONFINITONAN} - set all C{GDict} items to C{NAN} iff any argument is C{non-finite}. 

268 

269 C{REDUCEDLENGTH} - compute the reduced length C{m12}, 

270 

271 C{REVERSE2} - reverse azimuth C{azi2} by 180 degrees, 

272 

273 and C{ALL} - all of the above. 

274 

275 C{STANDARD} = C{AZIMUTH | DISTANCE | LATITUDE | LONGITUDE} 

276 

277 C{LINE_CAPS} = C{STANDARD | DISTANCE_IN | REDUCEDLENGTH | GEODESICSCALE} 

278 ''' 

279 EMPTY = 0 # formerly aka NONE 

280 _CAP_1 = 1 << 0 # for goedesici/-w 

281 _CAP_1p = 1 << 1 # for goedesici/-w 

282 _CAP_2 = 1 << 2 

283 _CAP_3 = 1 << 3 # for goedesici/-w 

284 _CAP_4 = 1 << 4 # for goedesicw 

285# = 1 << 5 # unused 

286# = 1 << 6 # unused 

287 _CAP_ALL = 0x1F 

288# _CAP_MASK = _CAP_ALL 

289 LATITUDE = 1 << 7 # compute latitude C{lat2} 

290 LONGITUDE = 1 << 8 | _CAP_3 # compute longitude C{lon2} 

291 AZIMUTH = 1 << 9 # azimuths C{azi1} and C{azi2} 

292 DISTANCE = 1 << 10 | _CAP_1 # compute distance C{s12} 

293 DISTANCE_IN = 1 << 11 | _CAP_1 | _CAP_1p # allow distance C{s12} in Direct 

294 REDUCEDLENGTH = 1 << 12 | _CAP_1 | _CAP_2 # compute reduced length C{m12} 

295 GEODESICSCALE = 1 << 13 | _CAP_1 | _CAP_2 # compute geodesic scales C{M12} and C{M21} 

296 AREA = 1 << 14 | _CAP_4 # compute area C{S12} 

297 ALL = 0x7F80 | _CAP_ALL # without LONG_UNROLL, LINE_OFF, NONFINITONAN, REVERSE2 and _DEBUG_* 

298 

299 STANDARD = AZIMUTH | DISTANCE | LATITUDE | LONGITUDE 

300 STANDARD_LINE = STANDARD | DISTANCE_IN # for goedesici/-w 

301 

302 LINE_CAPS = STANDARD_LINE | REDUCEDLENGTH | GEODESICSCALE # .geodesici only 

303 LONG_UNROLL = 1 << 15 # unroll C{lon2} in .Direct and .Position 

304 LINE_OFF = 1 << 16 # Line without updates from parent geodesic or rhumb 

305 NONFINITONAN = 1 << 17 # see method GDict._toNAN 

306 REVERSE2 = 1 << 18 # reverse C{azi2} 

307 

308 AZIMUTH_DISTANCE = AZIMUTH | DISTANCE 

309 AZIMUTH_DISTANCE_AREA = AZIMUTH | DISTANCE | AREA 

310 

311 LATITUDE_LONGITUDE = LATITUDE | LONGITUDE 

312 LATITUDE_LONGITUDE_AREA = LATITUDE | LONGITUDE | AREA 

313 

314 _SALP_CALPs_ = 1 << 19 # (INTERNAL) GeodesicExact._GenInverse 

315 

316 _DEBUG_AREA = 1 << 20 # (INTERNAL) include Line details 

317 _DEBUG_DIRECT = 1 << 21 # (INTERNAL) include Direct details 

318 _DEBUG_INVERSE = 1 << 22 # (INTERNAL) include Inverse details 

319 _DEBUG_LINE = 1 << 23 # (INTERNAL) include Line details 

320 _DEBUG_ALL = _DEBUG_AREA | _DEBUG_DIRECT | _DEBUG_INVERSE | \ 

321 _DEBUG_LINE | _SALP_CALPs_ 

322 

323 _DIRECT3 = AZIMUTH | LATITUDE | LONGITUDE # for goedesicw only 

324 _INVERSE3 = AZIMUTH | DISTANCE # for goedesicw only 

325 

326 _OUT_ALL = ALL # see geographiclib.geodesiccapabilities.py 

327 _OUT_MASK = ALL | LONG_UNROLL | NONFINITONAN | REVERSE2 | _DEBUG_ALL 

328 

329 _AZIMUTH_LATITUDE_LONGITUDE = AZIMUTH | LATITUDE | LONGITUDE 

330 _AZIMUTH_LATITUDE_LONG_UNROLL = AZIMUTH | LATITUDE | LONG_UNROLL 

331 _DEBUG_DIRECT_LINE = _DEBUG_DIRECT | _DEBUG_LINE 

332# _DISTANCE_IN_OUT = DISTANCE_IN & _OUT_MASK # == DISTANCE_IN in .gx, .gxline 

333 _REDUCEDLENGTH_GEODESICSCALE = REDUCEDLENGTH | GEODESICSCALE 

334# _REDUCEDLENGTH_GEODESICSCALE_DISTANCE = REDUCEDLENGTH | GEODESICSCALE | DISTANCE 

335 

336 def items(self): 

337 '''Yield all I{public} C{Caps} as 2-tuple C{(NAME, mask)}. 

338 ''' 

339 for n, C in type(self).__dict__.items(): 

340 if isint(C) and not n.startswith(_UNDER_) \ 

341 and n.replace(_UNDER_, NN).isupper(): 

342 yield n, C 

343 

344 def toStr(self, Cask, sep=_BAR_): 

345 '''Return C{Caps} or an C{outmask} as C{str} or tuple of C{str}s. 

346 ''' 

347 s = (n for n, C in sorted(Caps.items()) # PYLINT E1120 

348 if C and (Cask & C) == C) # and int1s(C) <= 3 

349 return sep.join(s) if sep else tuple(s) 

350 

351if _FOR_DOCS: # PYCHOK force ... 

352 pass # ... CAPS.__doc__ into epydoc 

353else: 

354 Caps = Caps() # PYCHOK singleton 

355 

356_key2Caps = dict(a12 =Caps.DISTANCE, # in GDict._toNAN, -._unCaps 

357 azi1=Caps.AZIMUTH, 

358 azi2=Caps.AZIMUTH, 

359 lat1=Caps.LATITUDE, 

360 lat2=Caps.LATITUDE, 

361 lon1=Caps.LONGITUDE, 

362 lon2=Caps.LONGITUDE, 

363 m12 =Caps.REDUCEDLENGTH, 

364 M12 =Caps.GEODESICSCALE, 

365 M21 =Caps.GEODESICSCALE, 

366 s12 =Caps.DISTANCE, 

367 S12 =Caps.AREA) 

368 

369 

370class _CapsBase(_NamedBase): # in .auxilats, .geodesicx.gxbases 

371 '''(INTERNAL) Base class for C{Geodesic*}, C{Geodesic*Exact}, C{Intersectool} and C{Rhumb*Base}. 

372 ''' 

373 ALL = Caps.ALL 

374 AREA = Caps.AREA 

375 AZIMUTH = Caps.AZIMUTH 

376 DISTANCE = Caps.DISTANCE 

377 DISTANCE_IN = Caps.DISTANCE_IN 

378 EMPTY = Caps.EMPTY # aka NONE 

379 GEODESICSCALE = Caps.GEODESICSCALE 

380 LATITUDE = Caps.LATITUDE 

381 LINE_CAPS = Caps.LINE_CAPS 

382 LINE_OFF = Caps.LINE_OFF 

383 LONGITUDE = Caps.LONGITUDE 

384 LONG_UNROLL = Caps.LONG_UNROLL 

385 NONFINITONAN = Caps.NONFINITONAN 

386 REDUCEDLENGTH = Caps.REDUCEDLENGTH 

387 STANDARD = Caps.STANDARD 

388 STANDARD_LINE = Caps.STANDARD_LINE # for geodesici 

389 

390 _caps = 0 # None 

391 _debug = 0 # or Caps._DEBUG_... 

392 

393 @property_RO 

394 def caps(self): 

395 '''Get the capabilities (bit-or'ed C{Caps}). 

396 ''' 

397 return self._caps 

398 

399 def caps_(self, caps): 

400 '''Check the available capabilities. 

401 

402 @arg caps: Bit-or'ed combination of L{Caps<pygeodesy.karney.Caps>} 

403 values for all capabilities to be checked. 

404 

405 @return: C{True} if I{all} B{C{caps}} are available, C{False} 

406 otherwise (C{bool}). 

407 ''' 

408 caps &= Caps._OUT_ALL 

409 return (self.caps & caps) == caps 

410 

411 @property 

412 def debug(self): 

413 '''Get the C{debug} option (C{bool}). 

414 ''' 

415 return bool(self._debug) 

416 

417 @debug.setter # PYCHOK setter! 

418 def debug(self, debug): 

419 '''Set the C{debug} option (C{bool}) to include 

420 more details in L{GDict} results. 

421 ''' 

422 self._debug = Caps._DEBUG_ALL if debug else 0 

423 

424 def _iter2tion(self, r, iter=None, **unused): # PYLINT W0622 

425 '''(INTERNAL) Copy C{C{s}.iter} into C{B{r}._iteration}. 

426 ''' 

427 if iter is not None: 

428 self._iteration = r._iteration = iter 

429 return r 

430 

431 

432class Direct9Tuple(_GTuple): 

433 '''9-Tuple C{(a12, lat2, lon2, azi2, s12, m12, M12, M21, S12)} with arc 

434 length C{a12}, angles C{lat2}, C{lon2} and azimuth C{azi2} in C{degrees}, 

435 distance C{s12} and reduced length C{m12} in C{meter} and area C{S12} in 

436 C{meter} I{squared}. 

437 ''' 

438 _Names_ = (_a12_, _lat2_, _lon2_, _azi2_, _s12_, _m12_, _M12_, _M21_, _S12_) 

439 _Units_ = (_Azi, _Lat, _Lon, _Azi, _M, _Pass, _Pass, _Pass, _M2) 

440 

441 

442class GDict(ADict): # XXX _NamedDict 

443 '''A C{dict} with both C{key} I{and} C{attribute} access to the C{dict} items. 

444 

445 Results of all C{geodesic} and C{rhumb} methods (with capitalized named) are 

446 returned as L{GDict} instances, see for example L{GeodesicExact} and L{RhumbAux}. 

447 ''' 

448 def toDirect9Tuple(self, dflt=NAN): 

449 '''Convert this L{GDict} result to a 9-tuple, like I{Karney}'s method 

450 C{geographiclib.geodesic.Geodesic._GenDirect}. 

451 

452 @kwarg dflt: Default value for missing items (C{any}). 

453 

454 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, 

455 s12, m12, M12, M21, S12)} 

456 ''' 

457 return self._toTuple(Direct9Tuple, dflt) 

458 

459 def toGeodSolve12Tuple(self, dflt=NAN): # PYCHOK 12 args 

460 '''Convert this L{GDict} result to a 12-Tuple, compatible with I{Karney}'s 

461 U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} 

462 result. 

463 

464 @kwarg dflt: Default value for missing items (C{any}). 

465 

466 @return: L{GeodSolve12Tuple}C{(lat1, lon1, azi1, lat2, lon2, azi2, 

467 s12, a12, m12, M12, M21, S12)}. 

468 ''' 

469 return self._toTuple(_MODS.geodsolve.GeodSolve12Tuple, dflt) 

470 

471 def toGeod3Solve8Tuple(self, dflt=NAN): # PYCHOK 12 args 

472 '''Convert this L{GDict} result to an 8-Tuple, compatible with I{Karney}'s 

473 U{Geod3Solve<https://GeographicLib.SourceForge.io/C++/doc/Geod3Solve.1.html>} 

474 result. 

475 

476 @kwarg dflt: Default value for missing items (C{any}). 

477 

478 @return: L{Geod3Solve8Tuple}C{(bet1, omg1, alp1, bet2, omg2, alp2, s12, a12)}. 

479 ''' 

480 return self._toTuple(_MODS.geod3solve.Geod3Solve8Tuple, dflt) 

481 

482 def toInverse10Tuple(self, dflt=NAN): 

483 '''Convert this L{GDict} result to a 10-tuple, like I{Karney}'s 

484 method C{geographiclib.geodesic.Geodesic._GenInverse}. 

485 

486 @kwarg dflt: Default value for missing items (C{any}). 

487 

488 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, 

489 m12, M12, M21, S12)}. 

490 ''' 

491 return self._toTuple(Inverse10Tuple, dflt) 

492 

493 def _toNAN(self, outmask, **specs): # .GeodesicExact._GDistInverse, .GeodesicLineExact._GenPosition 

494 '''(INTERNAL) Convert this C{GDict} to all C{NAN}s. 

495 ''' 

496 if (outmask & Caps.NONFINITONAN): 

497 def _t2(k): 

498 return k, specs.get(k, NAN) 

499 

500 d = dict(_t2(k) for k, C in _key2Caps.items() 

501 if (outmask & C) == C) 

502 self.set_(**d) # self.update(d) 

503 return self 

504 

505 @deprecated_method 

506 def toRhumb7Tuple(self, dflt=NAN): # PYCHOK no cover 

507 '''DEPRECATED on 23.12.07, use method C{toRhumb8Tuple}. 

508 

509 @return: A I{DEPRECATED} L{Rhumb7Tuple}. 

510 ''' 

511 return self._toTuple(_MODS.deprecated.classes.Rhumb7Tuple, dflt) 

512 

513 def toRhumb8Tuple(self, dflt=NAN): 

514 '''Convert this L{GDict} result to a 8-tuple. 

515 

516 @kwarg dflt: Default value for missing items (C{any}). 

517 

518 @return: L{Rhumb8Tuple}C{(lat1, lon1, lat2, lon2, azi12, 

519 s12, S12, a12)}. 

520 ''' 

521 return self._toTuple(Rhumb8Tuple, dflt) 

522 

523 def toRhumbSolve7Tuple(self, dflt=NAN): 

524 '''Convert this L{GDict} result to a 8-tuple. 

525 

526 @kwarg dflt: Default value for missing items (C{any}). 

527 

528 @return: L{RhumbSolve7Tuple}C{(lat1, lon1, lat2, lon2, 

529 azi12, s12, S12)}. 

530 ''' 

531 return self._toTuple(_MODS.rhumb.solve.RhumbSolve7Tuple, dflt) 

532 

533 def _toTuple(self, nTuple, dflt): 

534 '''(INTERNAL) Convert this C{GDict} to an B{C{nTuple}}. 

535 ''' 

536 t = tuple(self.get(n, dflt) for n in nTuple._Names_) 

537 return nTuple(t, iteration=self._iteration, name=self.name) # PYCHOK name 

538 

539 def _2X(self, gl, _2X=_X_): # .Intersectool, .Intersector 

540 '''(INTERNAL) Rename C{-2} attr to C{-X} or C{-M}. 

541 ''' 

542 X = GDict(self) 

543 for n in (_lat2_, _lon2_, _azi2_, _s12_, _a12_): 

544 if n in X: # X._X = X._2 

545 X[n[:-1] + _2X] = X.pop(n) 

546 v = getattr(gl, n, X) 

547 if v is not X: # X._2 = gl._2 

548 X[n] = v 

549 return X 

550 

551 def _unCaps(self, outmask): # in .geodsolve 

552 '''(INTERNAL) Remove superfluous items. 

553 ''' 

554 for k, C in _key2Caps.items(): 

555 if k in self and (outmask & C) != C: 

556 self.pop(k) # delattr(self, k) 

557 return self 

558 

559 

560class Inverse10Tuple(_GTuple): 

561 '''10-Tuple C{(a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12)} with arc length 

562 C{a12} in C{degrees}, distance C{s12} and reduced length C{m12} in C{meter}, area 

563 C{S12} in C{meter} I{squared} and the sines C{salp1}, C{salp2} and cosines C{calp1}, 

564 C{calp2} of the initial C{1} and final C{2} (forward) azimuths. 

565 ''' 

566 _Names_ = (_a12_, _s12_, 'salp1', 'calp1', 'salp2', 'calp2', _m12_, _M12_, _M21_, _S12_) 

567 _Units_ = (_Azi, _M, _Pass, _Pass, _Pass, _Pass, _Pass, _Pass, _Pass, _M2) 

568 

569 def toGDict(self, **updates): 

570 '''Convert this C{Inverse10Tuple} to a L{GDict}. 

571 

572 @kwarg updates: Optional items to apply (C{nam=value} pairs) 

573 ''' 

574 return _GTuple.toGDict(self, azi1=atan2d(self.salp1, self.calp1), # PYCHOK namedTuple 

575 azi2=atan2d(self.salp2, self.calp2), # PYCHOK namedTuple 

576 **updates) # PYCHOK indent 

577 

578 

579class Rhumb8Tuple(_GTuple): 

580 '''8-Tuple C{(lat1, lon1, lat2, lon2, azi12, s12, S12, a12)} with lat- C{lat1}, 

581 C{lat2} and longitudes C{lon1}, C{lon2} of both points, the azimuth of the 

582 rhumb line C{azi12}, the distance C{s12}, the area C{S12} under the rhumb 

583 line and the angular distance C{a12} between both points. 

584 ''' 

585 _Names_ = (_lat1_, _lon1_, _lat2_, _lon2_, _azi12_, _s12_, _S12_, _a12_) 

586 _Units_ = ( Lat, Lon, Lat, Lon, _Azi, _M, _M2, _Deg) 

587 

588 def toDirect9Tuple(self, dflt=NAN, **a12_azi1_azi2_m12_M12_M21): 

589 '''Convert this L{Rhumb8Tuple} result to a 9-tuple, like I{Karney}'s 

590 method C{geographiclib.geodesic.Geodesic._GenDirect}. 

591 

592 @kwarg dflt: Default value for missing items (C{any}). 

593 @kwarg a12_azi1_azi2_m12_M12_M21: Optional keyword arguments 

594 to specify or override L{Inverse10Tuple} items. 

595 

596 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, s12, 

597 m12, M12, M21, S12)} 

598 ''' 

599 d = dict(azi1=self.azi12, M12=_1_0, m12=self.s12, # PYCHOK attr 

600 azi2=self.azi12, M21=_1_0) # PYCHOK attr 

601 if a12_azi1_azi2_m12_M12_M21: 

602 d.update(a12_azi1_azi2_m12_M12_M21) 

603 return self._toTuple(Direct9Tuple, dflt, d) 

604 

605 def toInverse10Tuple(self, dflt=NAN, **a12_m12_M12_M21_salp1_calp1_salp2_calp2): 

606 '''Convert this L{Rhumb8Tuple} to a 10-tuple, like I{Karney}'s 

607 method C{geographiclib.geodesic.Geodesic._GenInverse}. 

608 

609 @kwarg dflt: Default value for missing items (C{any}). 

610 @kwarg a12_m12_M12_M21_salp1_calp1_salp2_calp2: Optional keyword 

611 arguments to specify or override L{Inverse10Tuple} items. 

612 

613 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, 

614 m12, M12, M21, S12)}. 

615 ''' 

616 s, c = _sincos2d(self.azi12) # PYCHOK attr 

617 d = dict(salp1=s, calp1=c, M12=_1_0, m12=self.s12, # PYCHOK attr 

618 salp2=s, calp2=c, M21=_1_0) 

619 if a12_m12_M12_M21_salp1_calp1_salp2_calp2: 

620 d.update(a12_m12_M12_M21_salp1_calp1_salp2_calp2) 

621 return self._toTuple(Inverse10Tuple, dflt, d) 

622 

623 @deprecated_method 

624 def _to7Tuple(self): # in DEPRECATED Rhumb7Tuple 

625 '''DEPRECATED, I{do not use!}''' 

626 return _MODS.deprecated.classes.Rhumb7Tuple(self[:-1]) 

627 

628 

629class _kWrapped(_CapsBase): # in .geodesicw 

630 '''(INTERNAL) Wrapper for some of I{Karney}'s U{geographiclib 

631 <https://PyPI.org/project/geographiclib>} classes. 

632 ''' 

633 

634 @property_ROnce 

635 def geographiclib(self): 

636 '''Lazily import C{geographiclib}, provided the U{geographiclib 

637 <https://PyPI.org/project/geographiclib>} package is installed, 

638 otherwise raise a C{LazyImportError}. 

639 ''' 

640 g = _xgeographiclib(type(self), 1, 49) 

641 from geographiclib.geodesic import Geodesic 

642 g.Geodesic = Geodesic 

643 from geographiclib.geodesicline import GeodesicLine 

644 g.GeodesicLine = GeodesicLine 

645 from geographiclib.geomath import Math 

646 g.Math = Math 

647# from geographiclib.polygonarea import PolygonArea 

648# g.PolygonArea = PolygonArea # see below 

649 return g 

650 

651 @property_ROnce 

652 def Math(self): 

653 '''Wrap the C{geomath.Math} class, provided the U{geographiclib 

654 <https://PyPI.org/project/geographiclib>} package is installed, 

655 otherwise C{None}. 

656 ''' 

657 try: 

658 g = self.geographiclib 

659 M = g.Math 

660 if _version_info(g) < (2,): 

661 if _K_2_0: 

662 M = None 

663# elif not _K_2_0: # XXX set 2.0? 

664# _K_2_0 = True 

665 except (AttributeError, ImportError): 

666 M = None 

667 return M 

668 

669 @property_RO 

670 def Math_K_2(self): 

671 return (_2_4_ if _K_2_4 else 

672 (_2_ if _K_2_0 else _1_)) if self.Math else NN 

673 

674 @property_ROnce 

675 def _PolygonArea(self): # lazy import 

676 from geographiclib.polygonarea import PolygonArea 

677 return PolygonArea 

678 

679 @property_ROnce 

680 def _Taus(self): # lazy import 

681 return _MODS.ellipsoids._Taus 

682 

683_wrapped = _kWrapped() # PYCHOK singleton, .datum, .test/base.py 

684 

685 

686class _Xables(object): 

687 '''(INTERNAL) Get I{Karney}'s executable paths from/and env vars. 

688 ''' 

689 bin_ = '/opt/local/bin/' # '/opt/local/Cellar/geographiclib/2.X/bin/' # HomeBrew on macOS 

690 ENV = NN 

691 

692 def GeoConvert(self, *dir_): 

693 return self._path(self.GeoConvert, *dir_) 

694 

695 def GeodSolve(self, *dir_): 

696 return self._path(self.GeodSolve, *dir_) 

697 

698 def Geod3Solve(self, *dir_): 

699 return self._path(self.Geod3Solve, *dir_) 

700 

701 def IntersectTool(self, *dir_): 

702 return self._path(self.IntersectTool, *dir_) 

703 

704 def RhumbSolve(self, *dir_): 

705 return self._path(self.RhumbSolve, *dir_) 

706 

707 def name_version(self, path, base=True): 

708 # return C{(path + ' ' + version)} of an executable 

709 if path: 

710 try: 

711 r, s = _popen2((path, '--version')) 

712 if base: 

713 path = _MODS.os.path.basename(path) # PYLINT E1101 

714 r = _SPACE_(path, r.split()[-1]) 

715 else: 

716 r = _MODS.streprs.Fmt.PARENSPACED(r, s) # PYLINT E1101 

717 return r 

718 except (IndexError, IOError, OSError): 

719 pass 

720 return NN 

721 

722 def _path(self, which, *dir_): 

723 n = typename(which, which) 

724 self.ENV = E = _PYGEODESY_ENV(n) 

725 return _getenv(E, NN) or (NN(dir_[0], n) if dir_ else E) 

726 

727 def X_not(self, path): 

728 return 'env %s=%r not executable' % (self.ENV, path) 

729 

730 def X_OK(self, path): # is C{path} an executable? 

731 os = _MODS.os # import os 

732 return os.access(path, os.X_OK) if path else False # PYLINT E1101 

733 

734_Xables = _Xables() # PYCHOK singleton 

735 

736 

737def _around(x): # in .utily.sincos2d 

738 '''I{Coarsen} a scalar by rounding small values to underflow to C{0.0}. 

739 

740 @return: Coarsened value (C{float}). 

741 

742 @see: I{Karney}'s U{geomath.Math.AngRound<https://SourceForge.net/p/ 

743 geographiclib/code/ci/release/tree/python/geographiclib/geomath.py>} 

744 ''' 

745 try: 

746 return _wrapped.Math.AngRound(x) 

747 except AttributeError: 

748 z = _1_16th 

749 w = z - fabs(x) 

750 if w > 0: # don't "simplify" z - (z - x) to x 

751 x = _copysign(z - w, x) 

752 return x 

753 

754 

755# def _atan2(y, x): 

756# # equivalent to math.atan2 

757# r = fabs(y) 

758# if r > fabs(x): # y, x = x, y 

759# r = atan2(x, r) 

760# # q = 3 if y < 0 else 2 

761# r = (r - PI_2) if y < 0 else (PI_2 - r) 

762# elif x < 0: # q = 1 

763# r = (-PI if y < 0 else PI) - atan2(y, -x) 

764# else: # q = 0 

765# r = atan2(y, x) 

766# return r 

767 

768 

769def _atan2d(y, x): 

770 '''Return C{atan2(B{y}, B{x})} in C{degrees}. 

771 ''' 

772 try: 

773 return _wrapped.Math.atan2d(y, x) 

774 except AttributeError: 

775 return atan2d(y, x) 

776 

777 

778def _cbrt(x): 

779 '''Return C{cubic root(B{x})}. 

780 ''' 

781 try: 

782 return _wrapped.Math.cbrt(x) 

783 except AttributeError: 

784 return cbrt(x) 

785 

786 

787def _copyBit(x, y): 

788 '''Like C{copysign0(B{x}, B{y})}, with C{B{x} > 0}. 

789 ''' 

790 return (-x) if _signBit(y) else x 

791 

792 

793def _2cos2x(cx, sx): # in .auxDST, .auxLat, .gxbases 

794 '''Return M{2 * cos(2 * x)} from cos(x) and sin(x). 

795 ''' 

796 r = cx - sx 

797 if r: 

798 r *= (cx + sx) * _2_0 

799 return r 

800 

801 

802def _diff182(deg0, deg, K_2_0=False): 

803 '''Compute C{deg - deg0}, reduced to C{[-180,180]} accurately. 

804 

805 @return: 2-Tuple C{(delta_angle, residual)} in C{degrees}. 

806 ''' 

807 try: 

808 return _wrapped.Math.AngDiff(deg0, deg) 

809 except AttributeError: 

810 if K_2_0 or _K_2_0: # geographiclib 2.0+ 

811 _r = fremainder 

812 d, t = _sum2(_r(-deg0, _360_0), 

813 _r( deg, _360_0)) 

814 d, t = _sum2(_r( d, _360_0), t) 

815 if _isin(d, _0_0, _180_0, _N_180_0): 

816 d = _copysign(d, (-t) if t else (deg - deg0)) 

817 else: 

818 _n = _norm180 

819 d, t = _sum2(_n(-deg0), _n(deg)) 

820 d = _n(d) 

821 if t > 0 and d == _180_0: 

822 d = _N_180_0 

823 d, t = _sum2(d, t) 

824 return d, t 

825 

826 

827def _eatanhe(x, es): 

828 '''Compute GeographicLib 2.7+ C{Math.eatanhe(x, es)}. 

829 ''' 

830 return _wrapped._Taus(es).eatanhe(x) 

831 

832 

833def _fix90(deg): # mimick Math.LatFix 

834 '''Replace angle in C{degrees} outside [-90,90] by NAN. 

835 

836 @return: Angle C{degrees} or NAN. 

837 ''' 

838 try: 

839 return _wrapped.Math.LatFix(deg) 

840 except AttributeError: 

841 return NAN if fabs(deg) > 90 else deg 

842 

843 

844def _isfinite(x): # mimick geomath.Math.isfinite 

845 '''Check finiteness of C{x}. 

846 

847 @return: C{True} if finite. 

848 ''' 

849 try: 

850 return _wrapped.Math.isfinite(x) 

851 except AttributeError: 

852 return _math_isfinite(x) # and fabs(x) <= _MAX 

853 

854 

855def _llz2gl(gl, **llz2): # see .geodesici._llz2G 

856 '''(INTERNAL) Set C{gl.lat2, .lon2, .azi2} from C{llz2}. 

857 ''' 

858 if llz2: 

859 for n in (_lat2_, _lon2_, _azi2_): # _lat1_, _lon1_, _azi1_ 

860 v = llz2.get(n, None) 

861 if v is not None: 

862 setattr(gl, n, v) 

863 return gl 

864 

865 

866def _norm2(x, y): # mimick geomath.Math.norm 

867 '''Normalize C{B{x}} and C{B{y}}. 

868 

869 @return: 2-Tuple of C{(B{x}, B{y})}, normalized. 

870 ''' 

871 try: 

872 return _wrapped.Math.norm(x, y) 

873 except AttributeError: 

874 return norm2(x, y) 

875 

876 

877def _norm180(deg): # mimick geomath.Math.AngNormalize 

878 '''Reduce angle in C{degrees} to (-180,180]. 

879 

880 @return: Reduced angle C{degrees}. 

881 ''' 

882 try: 

883 return _wrapped.Math.AngNormalize(deg) 

884 except AttributeError: 

885 d = fremainder(deg, _360_0) 

886 if _isin(d, _180_0, _N_180_0): 

887 d = _copysign(_180_0, deg) if _K_2_0 else _180_0 

888 return d 

889 

890 

891def _polygon(geodesic, points, closed, line, wrap, polar): 

892 '''(INTERNAL) Compute the area or perimeter of a polygon, 

893 using a L{GeodesicExact}, L{GeodesicSolve} or (if the 

894 C{geographiclib} package is installed) a C{Geodesic} 

895 or C{geodesicw.Geodesic} instance. 

896 ''' 

897 if not wrap: # capability LONG_UNROLL can't be off 

898 notImplemented(None, wrap=wrap, up=3) 

899 

900 if _MODS.booleans.isBoolean(points): 

901 # recursive call for each boolean clip 

902 

903 def _a_p(clip, *args, **unused): 

904 return _polygon(geodesic, clip, *args) 

905 

906 if not closed: # closed only 

907 raise _ValueError(closed=closed, points=_composite_) 

908 

909 return points._sum1(_a_p, closed, line, wrap, polar) 

910 

911 gP = geodesic.Polygon(line) 

912 _A = gP.AddPoint 

913