1

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

3

4 u'''Wrapper to invoke I{Karney}'s U{GeodSolve

5<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} utility

6 as an (exact) geodesic, but intended I{for testing purposes only}.

7

8 Set env variable C{PYGEODESY_GEODSOLVE} to the (fully qualified) path

9 of the C{GeodSolve} executable.

10'''

11

12# from pygeodesy.basics import _xinstanceof # from .karney

13 from pygeodesy.interns import NN, _a12_, _azi1_, _azi2_, \

14 _lat1_, _lat2_, _lon1_, _lon2_, _m12_, \

15 _M12_, _M21_, _s12_, _S12_, _UNDER_

16 from pygeodesy.interns import _UNUSED_, _not_ # PYCHOK used!

17 from pygeodesy.karney import _Azi, Caps, _Deg, GeodesicError, _GTuple, \

18 _Pass, _Lat, _Lon, _M, _M2, _sincos2d, \

19 _xinstanceof

20 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, \

21 _PYGEODESY_GEODSOLVE_, _getenv, printf

22 from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple

23 from pygeodesy.props import Property, Property_RO

24 from pygeodesy.solveBase import _SolveBase, _SolveLineBase, wrap360

25# from pygeodesy.utily import wrap360 # from .solveBase

26

27__all__ = _ALL_LAZY.geodsolve

28__version__ = '23.09.14'

29

30

31 class GeodSolve12Tuple(_GTuple):

32 '''12-Tuple C{(lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12)} with

33 angles C{lat1}, C{lon1}, C{azi1}, C{lat2}, C{lon2} and C{azi2} and arc C{a12} all in

34 C{degrees}, initial C{azi1} and final C{azi2} forward azimuths, distance C{s12} and

35 reduced length C{m12} in C{meter}, area C{S12} in C{meter} I{squared} and geodesic

36 scale factors C{M12} and C{M21}, both C{scalar}, see U{GeodSolve

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

38 '''

39 # from GeodSolve --help option -f ... lat1 lon1 azi1 lat2 lon2 azi2 s12 a12 m12 M12 M21 S12

40 _Names_ = (_lat1_, _lon1_, _azi1_, _lat2_, _lon2_, _azi2_, _s12_, _a12_, _m12_, _M12_, _M21_, _S12_)

41 _Units_ = (_Lat, _Lon, _Azi, _Lat, _Lon, _Azi, _M, _Deg, _Pass, _Pass, _Pass, _M2)

42

43

44 class _GeodesicSolveBase(_SolveBase):

45 '''(INTERNAL) Base class for L{GeodesicSolve} and L{GeodesicLineSolve}.

46 '''

47 _Error = GeodesicError

48 _Names_Direct = \

49 _Names_Inverse = GeodSolve12Tuple._Names_

50 _Solve_name = 'GeodSolve'

51 _Solve_path = _getenv(_PYGEODESY_GEODSOLVE_, _PYGEODESY_GEODSOLVE_)

52

53 @Property_RO

54 def _b_option(self):

55 return ('-b',) if self.reverse2 else ()

56

57 @Property_RO

58 def _cmdBasic(self):

59 '''(INTERNAL) Get the basic C{GeodSolve} cmd (C{tuple}).

60 '''

61 return (self.GeodSolve,) + self._b_option \

62 + self._e_option \

63 + self._E_option \

64 + self._p_option \

65 + self._u_option + ('-f',)

66

67 @Property_RO

68 def _E_option(self):

69 return ('-E',) if self.Exact else ()

70

71 @Property

72 def GeodSolve(self):

73 '''Get the U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>}

74 executable (C{filename}).

75 '''

76 return self._Solve_path

77

78 @GeodSolve.setter # PYCHOK setter!

79 def GeodSolve(self, path):

80 '''Set the U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>}

81 executable (C{filename}), the (fully qualified) path to the C{GeodSolve} executable.

82

83 @raise GeodesicError: Invalid B{C{path}}, B{C{path}} doesn't exist or

84 isn't the C{GeodSolve} executable.

85 '''

86 self._setSolve(path)

87

88 def toStr(self, **prec_sep): # PYCHOK signature

89 '''Return this C{GeodesicSolve} as string.

90

91 @kwarg prec_sep: Keyword argumens C{B{prec}=6} and C{B{sep}=', '}

92 for the C{float} C{prec}ision, number of decimal digits

93 (0..9) and the C{sep}arator string to join. Trailing

94 zero decimals are stripped for B{C{prec}} values of

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

96

97 @return: GeodesicSolve items (C{str}).

98 '''

99 return _SolveBase._toStr(self, GeodSolve=self.GeodSolve, **prec_sep)

100

101 @Property_RO

102 def _u_option(self):

103 return ('-u',) if self.unroll else ()

104

105

106 class GeodesicSolve(_GeodesicSolveBase):

107 '''Wrapper to invoke I{Karney}'s U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>}

108 as an C{Exact} version of I{Karney}'s Python class U{Geodesic<https://GeographicLib.SourceForge.io/C++/doc/

109 python/code.html#geographiclib.geodesic.Geodesic>}.

110

111 @note: Use property C{GeodSolve} or env variable C{PYGEODESY_GEODSOLVE} to specify the (fully

112 qualified) path to the C{GeodSolve} executable.

113

114 @note: This C{geodesic} is intended I{for testing purposes only}, it invokes the C{GeodSolve}

115 executable for I{every} method call.

116 '''

117

118 def Area(self, polyline=False, name=NN):

119 '''Set up a L{GeodesicAreaExact} to compute area and

120 perimeter of a polygon.

121

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

123 area and perimeter (C{bool}).

124 @kwarg name: Optional name (C{str}).

125

126 @return: A L{GeodesicAreaExact} instance.

127

128 @note: The B{C{debug}} setting is passed as C{verbose}

129 to the returned L{GeodesicAreaExact} instance.

130 '''

131 gaX = _MODS.geodesicx.GeodesicAreaExact(self, polyline=polyline,

132 name=name or self.name)

133 if self.verbose or self.debug: # PYCHOK no cover

134 gaX.verbose = True

135 return gaX

136

137 Polygon = Area # for C{geographiclib} compatibility

138

139 def Direct3(self, lat1, lon1, azi1, s12): # PYCHOK outmask

140 '''Return the destination lat, lon and reverse azimuth

141 (final bearing) in C{degrees}.

142

143 @return: L{Destination3Tuple}C{(lat, lon, final)}.

144 '''

145 r = self._GDictDirect(lat1, lon1, azi1, False, s12, floats=False)

146 return Destination3Tuple(float(r.lat2), float(r.lon2), wrap360(r.azi2),

147 iteration=r._iteration)

148

149 def Inverse3(self, lat1, lon1, lat2, lon2): # PYCHOK outmask

150 '''Return the distance in C{meter} and the forward and

151 reverse azimuths (initial and final bearing) in C{degrees}.

152

153 @return: L{Distance3Tuple}C{(distance, initial, final)}.

154 '''

155 r = self._GDictInverse(lat1, lon1, lat2, lon2, floats=False)

156 return Distance3Tuple(float(r.s12), wrap360(r.azi1), wrap360(r.azi2),

157 iteration=r._iteration)

158

159 def Line(self, lat1, lon1, azi1, caps=Caps.ALL):

160 '''Set up a L{GeodesicLineSolve} to compute several points

161 on a single geodesic.

162

163 @arg lat1: Latitude of the first point (C{degrees}).

164 @arg lon1: Longitude of the first point (C{degrees}).

165 @arg azi1: Azimuth at the first point (compass C{degrees}).

166 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying

167 the capabilities the L{GeodesicLineSolve} instance

168 should possess, always C{Caps.ALL}.

169

170 @return: A L{GeodesicLineSolve} instance.

171

172 @note: If the point is at a pole, the azimuth is defined by keeping

173 B{C{lon1}} fixed, writing C{B{lat1} = ±(90 − ε)}, and taking

174 the limit C{ε → 0+}.

175

176 @see: C++ U{GeodesicExact.Line

177 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}

178 and Python U{Geodesic.Line<https://GeographicLib.SourceForge.io/Python/doc/code.html>}.

179 '''

180 return GeodesicLineSolve(self, lat1, lon1, azi1, caps=caps, name=self.name)

181

182

183 class GeodesicLineSolve(_GeodesicSolveBase, _SolveLineBase):

184 '''Wrapper to invoke I{Karney}'s U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>}

185 as an C{Exact} version of I{Karney}'s Python class U{GeodesicLine<https://GeographicLib.SourceForge.io/C++/doc/

186 python/code.html#geographiclib.geodesicline.GeodesicLine>}.

187

188 @note: Use property C{GeodSolve} or env variable C{PYGEODESY_GEODSOLVE} to specify the (fully

189 qualified) path to the C{GeodSolve} executable.

190

191 @note: This C{geodesic} is intended I{for testing purposes only}, it invokes the C{GeodSolve}

192 executable for I{every} method call.

193 '''

194

195 def __init__(self, geodesic, lat1, lon1, azi1, caps=Caps.ALL, name=NN):

196 '''New L{GeodesicLineSolve} instance, allowing points to be found along

197 a geodesic starting at C{(B{lat1}, B{lon1})} with azimuth B{C{azi1}}.

198

199 @arg geodesic: The geodesic to use (L{GeodesicSolve}).

200 @arg lat1: Latitude of the first point (C{degrees}).

201 @arg lon1: Longitude of the first point (C{degrees}).

202 @arg azi1: Azimuth at the first points (compass C{degrees}).

203 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying

204 the capabilities the L{GeodesicLineSolve} instance

205 should possess, always C{Caps.ALL}. Use C{Caps.LINE_OFF}

206 if updates to the B{C{geodesic}} should I{not} be

207 reflected in this L{GeodesicLineSolve} instance.

208 @kwarg name: Optional name (C{str}).

209

210 @raise GeodesicError: Invalid path for the C{GeodSolve} executable or

211 or isn't the C{GeodSolve} executable, see

212 property C{geodesic.GeodSolve}.

213

214 @raise TypeError: Invalid B{C{geodesic}}.

215 '''

216 _xinstanceof(GeodesicSolve, geodesic=geodesic)

217 if (caps & Caps.LINE_OFF): # copy to avoid updates

218 geodesic = geodesic.copy(deep=False, name=NN(_UNDER_, geodesic.name))

219 _SolveLineBase.__init__(self, geodesic, lat1, lon1, caps, name, azi1=azi1)

220 try:

221 self.GeodSolve = geodesic.GeodSolve # geodesic or copy of geodesic

222 except GeodesicError:

223 pass

224

225 def ArcPosition(self, a12, outmask=_UNUSED_): # PYCHOK unused

226 '''Find the position on the line given B{C{a12}}.

227

228 @arg a12: Spherical arc length from the first point to the

229 second point (C{degrees}).

230

231 @return: A C{GDict} with 12 items C{lat1, lon1, azi1, lat2, lon2,

232 azi2, m12, a12, s12, M12, M21, S12}.

233 '''

234 return self._GDictInvoke(self._cmdArc, True, self._Names_Direct, a12)

235

236 @Property_RO

237 def azi1(self):

238 '''Get the azimuth at the first point (compass C{degrees}).

239 '''

240 return self._lla1.azi1

241

242 azi12 = azi1 # like RhumbLineSolve

243

244 @Property_RO

245 def azi1_sincos2(self):

246 '''Get the sine and cosine of the first point's azimuth (2-tuple C{(sin, cos)}).

247 '''

248 return _sincos2d(self.azi1)

249

250 azi12_sincos2 = azi1_sincos2

251

252 @Property_RO

253 def _cmdArc(self):

254 '''(INTERNAL) Get the C{GeodSolve} I{-a -L} cmd (C{tuple}).

255 '''

256 return self._cmdDistance + ('-a',)

257

258 def Position(self, s12, outmask=_UNUSED_): # PYCHOK unused

259 '''Find the position on the line given B{C{s12}}.

260

261 @arg s12: Distance from the first point to the second (C{meter}).

262

263 @return: A C{GDict} with 12 items C{lat1, lon1, azi1, lat2, lon2,

264 azi2, m12, a12, s12, M12, M21, S12}, possibly C{a12=NAN}.

265 '''

266 return self._GDictInvoke(self._cmdDistance, True, self._Names_Direct, s12)

267

268 def toStr(self, **prec_sep): # PYCHOK signature

269 '''Return this C{GeodesicLineSolve} as string.

270

271 @kwarg prec_sep: Keyword argumens C{B{prec}=6} and C{B{sep}=', '}

272 for the C{float} C{prec}ision, number of decimal digits

273 (0..9) and the C{sep}arator string to join. Trailing

274 zero decimals are stripped for B{C{prec}} values of

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

276

277 @return: GeodesicLineSolve items (C{str}).

278 '''

279 return _SolveLineBase._toStr(self, azi1=self.azi1, geodesic=self._solve,

280 GeodSolve=self.GeodSolve, **prec_sep)

281

282

283__all__ += _ALL_DOCS(_GeodesicSolveBase)

284

285 if __name__ == '__main__':

286

287 from sys import argv

288

289 gS = GeodesicSolve(name='Test')

290 gS.verbose = '--verbose' in argv # or '-v' in argv

291

292 if gS.GeodSolve in (_PYGEODESY_GEODSOLVE_, None): # not set

293 gS.GeodSolve = '/opt/local/bin/GeodSolve' # '/opt/local/Cellar/geographiclib/1.51/bin/GeodSolve' # HomeBrew

294 printf('version: %s', gS.version)

295

296 r = gS.Direct(40.6, -73.8, 51, 5.5e6)

297 printf('Direct: %r', r, nl=1)

298 printf('Direct3: %r', gS.Direct3(40.6, -73.8, 51, 5.5e6))

299

300 printf('Inverse: %r', gS.Inverse( 40.6, -73.8, 51.6, -0.5), nl=1)

301 printf('Inverse1: %r', gS.Inverse1(40.6, -73.8, 51.6, -0.5))

302 printf('Inverse3: %r', gS.Inverse3(40.6, -73.8, 51.6, -0.5))

303

304 glS = GeodesicLineSolve(gS, 40.6, -73.8, 51, name='LineTest')

305 p = glS.Position(5.5e6)

306 printf('Position: %s %r', p == r, p, nl=1)

307 p = glS.ArcPosition(49.475527)

308 printf('ArcPosition: %s %r', p == r, p)

309

310# % python3 -m pygeodesy.geodsolve

311

312# version: /opt/local/bin/GeodSolve: GeographicLib version 1.51

313

314# version: /opt/local/bin/GeodSolve: GeographicLib version 1.51

315

316# Direct: GDict(M12=0.650911, M21=0.651229, S12=39735075134877.09375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141173, m12=4844148.703101, s12=5500000.0)

317# Direct3: Destination3Tuple(lat=51.884565, lon=-1.141173, final=107.189397)

318

319# Inverse: GDict(M12=0.64473, M21=0.645046, S12=40041368848742.53125, a12=49.94131, azi1=51.198883, azi2=107.821777, lat1=40.6, lat2=51.6, lon1=-73.8, lon2=-0.5, m12=4877684.602706, s12=5551759.400319)

320# Inverse1: 49.94131021789904

321# Inverse3: Distance3Tuple(distance=5551759.400319, initial=51.198883, final=107.821777)

322

323# Position: True GDict(M12=0.650911, M21=0.651229, S12=39735075134877.09375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141173, m12=4844148.703101, s12=5500000.0)

324# ArcPosition: False GDict(M12=0.650911, M21=0.651229, S12=39735074737272.734375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141174, m12=4844148.669561, s12=5499999.948497)

325

326

327# % python3 -m pygeodesy.geodsolve --verbose

328

329# GeodesicSolve 'Test' 1: /opt/local/bin/GeodSolve --version (invoke)

330# GeodesicSolve 'Test' 1: /opt/local/bin/GeodSolve: GeographicLib version 1.51 (0)

331# version: /opt/local/bin/GeodSolve: GeographicLib version 1.51

332# GeodesicSolve 'Test' 2: /opt/local/bin/GeodSolve -E -p 10 -f \ 40.600000000000001 -73.799999999999997 51.0 5500000.0 (Direct)

333# GeodesicSolve 'Test' 2: lat1=40.600000000000001, lon1=-73.799999999999997, azi1=51.0, lat2=51.884564505606761, lon2=-1.141172861200829, azi2=107.189397162605886, s12=5500000.0, a12=49.475527463251467, m12=4844148.703101486, M12=0.65091056699808603, M21=0.65122865892196558, S12=39735075134877.094 (0)

334

335# Direct: GDict(M12=0.650911, M21=0.651229, S12=39735075134877.09375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141173, m12=4844148.703101, s12=5500000.0)

336# GeodesicSolve 'Test' 3: /opt/local/bin/GeodSolve -E -p 10 -f \ 40.600000000000001 -73.799999999999997 51.0 5500000.0 (Direct3)

337# GeodesicSolve 'Test' 3: lat1=40.600000000000001, lon1=-73.799999999999997, azi1=51.0, lat2=51.884564505606761, lon2=-1.141172861200829, azi2=107.189397162605886, s12=5500000.0, a12=49.475527463251467, m12=4844148.703101486, M12=0.65091056699808603, M21=0.65122865892196558, S12=39735075134877.094 (0)

338# Direct3: Destination3Tuple(lat=51.884565, lon=-1.141173, final=107.189397)

339# GeodesicSolve 'Test' 4: /opt/local/bin/GeodSolve -E -p 10 -f -i \ 40.600000000000001 -73.799999999999997 51.600000000000001 -0.5 (Inverse)

340# GeodesicSolve 'Test' 4: lat1=40.600000000000001, lon1=-73.799999999999997, azi1=51.198882845579824, lat2=51.600000000000001, lon2=-0.5, azi2=107.821776735514248, s12=5551759.4003186841, a12=49.941310217899037, m12=4877684.6027061976, M12=0.64472969205948238, M21=0.64504567852134398, S12=40041368848742.531 (0)

341

342# Inverse: GDict(M12=0.64473, M21=0.645046, S12=40041368848742.53125, a12=49.94131, azi1=51.198883, azi2=107.821777, lat1=40.6, lat2=51.6, lon1=-73.8, lon2=-0.5, m12=4877684.602706, s12=5551759.400319)

343# GeodesicSolve 'Test' 5: /opt/local/bin/GeodSolve -E -p 10 -f -i \ 40.600000000000001 -73.799999999999997 51.600000000000001 -0.5 (Inverse1)

344# GeodesicSolve 'Test' 5: lat1=40.600000000000001, lon1=-73.799999999999997, azi1=51.198882845579824, lat2=51.600000000000001, lon2=-0.5, azi2=107.821776735514248, s12=5551759.4003186841, a12=49.941310217899037, m12=4877684.6027061976, M12=0.64472969205948238, M21=0.64504567852134398, S12=40041368848742.531 (0)

345# Inverse1: 49.94131021789904

346# GeodesicSolve 'Test' 6: /opt/local/bin/GeodSolve -E -p 10 -f -i \ 40.600000000000001 -73.799999999999997 51.600000000000001 -0.5 (Inverse3)

347# GeodesicSolve 'Test' 6: lat1=40.600000000000001, lon1=-73.799999999999997, azi1=51.198882845579824, lat2=51.600000000000001, lon2=-0.5, azi2=107.821776735514248, s12=5551759.4003186841, a12=49.941310217899037, m12=4877684.6027061976, M12=0.64472969205948238, M21=0.64504567852134398, S12=40041368848742.531 (0)

348# Inverse3: Distance3Tuple(distance=5551759.400319, initial=51.198883, final=107.821777)

349

350# Position: True GDict(M12=0.650911, M21=0.651229, S12=39735075134877.09375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141173, m12=4844148.703101, s12=5500000.0)

351# ArcPosition: False GDict(M12=0.650911, M21=0.651229, S12=39735074737272.734375, a12=49.475527, azi1=51.0, azi2=107.189397, lat1=40.6, lat2=51.884565, lon1=-73.8, lon2=-1.141174, m12=4844148.669561, s12=5499999.948497)

352

353

354# **) MIT License

355#

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

357#

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

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

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

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

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

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

364#

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

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

367#

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

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

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

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

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

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

374# OTHER DEALINGS IN THE SOFTWARE.