Coverage for pygeodesy/rhumbBase.py: 96%

317 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-09-20 13:43 -0400

1 

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

3 

4u'''Base classes C{RhumbBase} and C{RhumbLineBase}, pure Python version of I{Karney}'s C++ 

5classes U{Rhumb<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Rhumb.html>} 

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

7from I{GeographicLib versions 2.0} and I{2.2} and I{Karney}'s C++ example U{Rhumb intersect 

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

9 

10Class L{RhumbLine} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively 

11find the intersection of two rhumb lines, respectively the nearest point on a rumb line along a 

12geodesic or perpendicular rhumb line from an other point. 

13 

14For more details, see the C++ U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>} 

15documentation, especially the U{Class List<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>}, 

16the background information on U{Rhumb lines<https://GeographicLib.SourceForge.io/C++/doc/rhumb.html>}, 

17the utily U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online 

18rhumb line calculations<https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve>}. 

19 

20Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2014-2023) and licensed under the MIT/X11 

21License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

22''' 

23# make sure int/int division yields float quotient 

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

25 

26# from pygeodesy.basics import unsigned0, _xinstanceof # from .karney 

27from pygeodesy.constants import EPS, EPS1, INT0, _EPSqrt as _TOL, NAN, \ 

28 _0_0, _0_01, _1_0, _90_0, _over 

29# from pygeodesy.datums import _spherical_datum # from .formy 

30# from pygeodesy.ellipsoids import _EWGS84 # from .karney 

31from pygeodesy.errors import IntersectionError, itemsorted, RhumbError, \ 

32 _xdatum, _xkwds, _Xorder 

33# from pygeodesy.etm import ExactTransverseMercator # _MODS 

34from pygeodesy.fmath import euclid, favg, fabs, Fsum 

35from pygeodesy.formy import opposing, _spherical_datum 

36# from pygeodesy.fsums import Fsum # from .fmath 

37from pygeodesy.interns import NN, _coincident_, _COMMASPACE_, _Dash, \ 

38 _intersection_, _no_, _parallel_, _under 

39from pygeodesy.karney import _atan2d, Caps, _CapsBase, _diff182, _fix90, \ 

40 GDict, _norm180, _EWGS84, unsigned0, _xinstanceof 

41# from pygeodesy.ktm import KTransverseMercator, _AlpCoeffs # _MODS 

42from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS 

43# from pygeodesy.named import notOverloaded # _MODS 

44from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple 

45from pygeodesy.props import deprecated_method, Property, Property_RO, property_RO, \ 

46 _update_all 

47from pygeodesy.streprs import Fmt, pairs, unstr 

48from pygeodesy.units import Float_, Lat, Lon, Meter, Int # PYCHOK shared 

49from pygeodesy.utily import atan2d, _azireversed, _loneg, sincos2d, sincos2d_, \ 

50 _unrollon, _Wrap, sincos2_ # PYCHOK shared 

51from pygeodesy.vector3d import _intersect3d3, Vector3d # in .intersection2 below 

52 

53# from math import fabs # from .fmath 

54 

55__all__ = () 

56__version__ = '23.09.20' 

57 

58_anti_ = _Dash('anti') 

59_rls = [] # instances of C{RbumbLine...} to be updated 

60_TRIPS = 65 # .intersection2, .nearestOn4, 19+ 

61 

62 

63class _Lat(Lat): 

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

65 ''' 

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

67 kwds = _xkwds(Error_name, clip=0, Error=RhumbError) 

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

69 

70 

71class _Lon(Lon): 

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

73 ''' 

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

75 kwds = _xkwds(Error_name, clip=0, Error=RhumbError) 

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

77 

78 

79def _update_all_rls(r): 

80 '''(INTERNAL) Zap cached/memoized C{Property[_RO]}s 

81 of any C{RhumbLine} instances tied to the given 

82 C{Rhumb} instance B{C{r}}. 

83 ''' 

84 # _xinstanceof(_MODS.rhumbaux.RhumbAux, _MODS.rhumbx.Rhumb, r=r) 

85 _update_all(r) 

86 for rl in _rls: # PYCHOK use weakref? 

87 if rl._rhumb is r: 

88 _update_all(rl) 

89 

90 

91class RhumbBase(_CapsBase): 

92 '''(INTERNAL) Base class for C{rhumbaux.RhumbAux} and C{rhumbx.Rhumb}. 

93 ''' 

94 _E = _EWGS84 

95 _exact = True 

96 _f_max = _0_01 

97 _mTM = 6 # see .TMorder 

98 

99 def __init__(self, a_earth, f, exact, name): 

100 '''New C{rhumbaux.RhumbAux} or C{rhumbx.Rhum}. 

101 ''' 

102 if f is not None: 

103 self.ellipsoid = a_earth, f 

104 elif a_earth not in (_EWGS84, None): 

105 self.ellipsoid = a_earth 

106 if not exact: 

107 self.exact = False 

108 if name: 

109 self.name = name 

110 

111 @Property_RO 

112 def a(self): 

113 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}). 

114 ''' 

115 return self.ellipsoid.a 

116 

117 equatoradius = a 

118 

119 def ArcDirect(self, lat1, lon1, azi12, a12, outmask=Caps.LATITUDE_LONGITUDE): 

120 '''Solve the I{direct rhumb} problem, optionally with area. 

121 

122 @arg lat1: Latitude of the first point (C{degrees90}). 

123 @arg lon1: Longitude of the first point (C{degrees180}). 

124 @arg azi12: Azimuth of the rhumb line (compass C{degrees}). 

125 @arg a12: Angle along the rhumb line from the given to the 

126 destination point (C{degrees}), can be negative. 

127 

128 @return: L{GDict} with 2 up to 8 items C{lat2, lon2, a12, S12, 

129 lat1, lon1, azi12, s12} with the destination point's 

130 latitude C{lat2} and longitude C{lon2} in C{degrees}, 

131 the rhumb angle C{a12} in C{degrees} and area C{S12} 

132 under the rhumb line in C{meter} I{squared}. 

133 

134 @raise ImportError: Package C{numpy} not found or not installed, 

135 only required for area C{S12} when C{B{exact} 

136 is True} and L{RhumbAux}. 

137 

138 @note: If B{C{a12}} is large enough that the rhumb line crosses 

139 a pole, the longitude of the second point is indeterminate 

140 and C{NAN} is returned for C{lon2} and area C{S12}. 

141 

142 @note: If the given point is a pole, the cosine of its latitude is 

143 taken to be C{sqrt(L{EPS})}. This position is extremely 

144 close to the actual pole and allows the calculation to be 

145 carried out in finite terms. 

146 ''' 

147 s12 = a12 * self._mpd 

148 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask) 

149 

150 @Property_RO 

151 def b(self): 

152 '''Get the C{ellipsoid}'s polar radius, semi-axis (C{meter}). 

153 ''' 

154 return self.ellipsoid.b 

155 

156 polaradius = b 

157 

158 def _Direct(self, ll1, azi12, s12, **outmask): 

159 '''(INTERNAL) Short-cut version, see .latlonBase. 

160 ''' 

161 return self.Direct(ll1.lat, ll1.lon, azi12, s12, **outmask) 

162 

163 def Direct(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE): 

164 '''Solve the I{direct rhumb} problem, optionally with area. 

165 

166 @arg lat1: Latitude of the first point (C{degrees90}). 

167 @arg lon1: Longitude of the first point (C{degrees180}). 

168 @arg azi12: Azimuth of the rhumb line (compass C{degrees}). 

169 @arg s12: Distance along the rhumb line from the given to 

170 the destination point (C{meter}), can be negative. 

171 

172 @return: L{GDict} with 2 up to 8 items C{lat2, lon2, a12, S12, 

173 lat1, lon1, azi12, s12} with the destination point's 

174 latitude C{lat2} and longitude C{lon2} in C{degrees}, 

175 the rhumb angle C{a12} in C{degrees} and area C{S12} 

176 under the rhumb line in C{meter} I{squared}. 

177 

178 @raise ImportError: Package C{numpy} not found or not installed, 

179 only required for area C{S12} when C{B{exact} 

180 is True} and L{RhumbAux}. 

181 

182 @note: If B{C{s12}} is large enough that the rhumb line crosses 

183 a pole, the longitude of the second point is indeterminate 

184 and C{NAN} is returned for C{lon2} and area C{S12}. 

185 

186 @note: If the given point is a pole, the cosine of its latitude is 

187 taken to be C{sqrt(L{EPS})}. This position is extremely 

188 close to the actual pole and allows the calculation to be 

189 carried out in finite terms. 

190 ''' 

191 a12 = _over(s12, self._mpd) 

192 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask) 

193 

194 def Direct8(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA): 

195 '''Like method L{Rhumb.Direct} but returning a L{Rhumb8Tuple} with area C{S12}. 

196 ''' 

197 return self.Direct(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple() 

198 

199 def _DirectLine(self, ll1, azi12, **caps_name): 

200 '''(INTERNAL) Short-cut version, see .latlonBase. 

201 ''' 

202 return self.DirectLine(ll1.lat, ll1.lon, azi12, **caps_name) 

203 

204 def DirectLine(self, lat1, lon1, azi12, **caps_name): 

205 '''Define a C{RhumbLine} in terms of the I{direct} rhumb 

206 problem to compute several points on a single rhumb line. 

207 

208 @arg lat1: Latitude of the first point (C{degrees90}). 

209 @arg lon1: Longitude of the first point (C{degrees180}). 

210 @arg azi12: Azimuth of the rhumb line (compass C{degrees}). 

211 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and 

212 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of 

213 L{Caps} values specifying the required capabilities. 

214 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}} 

215 should I{not} be reflected in this rhumb line. 

216 

217 @return: A C{RhumbLine...} instance and invoke its method 

218 C{.Position} to compute each point. 

219 

220 @note: Updates to this rhumb are reflected in the returned 

221 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}. 

222 ''' 

223 return self._RhumbLine(self, lat1, lon1, azi12, **caps_name) 

224 

225 Line = DirectLine # synonyms 

226 

227 def _DirectRhumb(self, lat1, lon1, azi12, a12, s12, outmask): 

228 '''(INTERNAL) See methods C{.ArcDirect} and C{.Direct}. 

229 ''' 

230 rl = self._RhumbLine(self, lat1, lon1, azi12, caps=Caps.LINE_OFF, 

231 name=self.name) 

232 return rl._Position(a12, s12, outmask | self._debug) # lat2, lon2, S12 

233 

234 @Property 

235 def ellipsoid(self): 

236 '''Get this rhumb's ellipsoid (L{Ellipsoid}). 

237 ''' 

238 return self._E 

239 

240 @ellipsoid.setter # PYCHOK setter! 

241 def ellipsoid(self, a_earth_f): 

242 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or 

243 L{a_f2Tuple}) or (equatorial) radius and flattening (2-tuple C{(a, f)}). 

244 

245 @raise RhumbError: If C{abs(B{f}} exceeds non-zero C{f_max} and C{exact=False}. 

246 ''' 

247 E = _spherical_datum(a_earth_f, Error=RhumbError).ellipsoid 

248 if self._E != E: 

249 self._exactest(self.exact, E, self.f_max) 

250 _update_all_rls(self) 

251 self._E = E 

252 

253 @Property 

254 def exact(self): 

255 '''Get the I{exact} option (C{bool}). 

256 ''' 

257 return self._exact 

258 

259 @exact.setter # PYCHOK setter! 

260 def exact(self, exact): 

261 '''Set the I{exact} option (C{bool}). If C{True}, use I{exact} rhumb 

262 expressions, otherwise a series expansion (accurate for oblate or 

263 prolate ellipsoids with C{abs(flattening)} below C{f_max}. 

264 

265 @raise RhumbError: If C{B{exact}=False} and C{abs(flattening}) 

266 exceeds non-zero C{f_max}. 

267 

268 @see: Option U{B{-s}<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} 

269 and U{ACCURACY<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html#ACCURACY>}. 

270 ''' 

271 x = bool(exact) 

272 if self._exact != x: 

273 self._exactest(x, self.ellipsoid, self.f_max) 

274 _update_all_rls(self) 

275 self._exact = x 

276 

277 def _exactest(self, exact, ellipsoid, f_max): 

278 # Helper for property setters C{ellipsoid}, C{exact} and C{f_max} 

279 if fabs(ellipsoid.f) > f_max > 0 and not exact: 

280 raise RhumbError(exact=exact, f=ellipsoid.f, f_max=f_max) 

281 

282 @Property_RO 

283 def f(self): 

284 '''Get the C{ellipsoid}'s flattening (C{float}). 

285 ''' 

286 return self.ellipsoid.f 

287 

288 flattening = f 

289 

290 @property 

291 def f_max(self): 

292 '''Get the I{max.} flattening (C{float}). 

293 ''' 

294 return self._f_max 

295 

296 @f_max.setter # PYCHOK setter! 

297 def f_max(self, f_max): # PYCHOK no cover 

298 '''Set the I{max.} flattening, not to exceed (C{float}). 

299 

300 @raise RhumbError: If C{exact=False} and C{abs(flattening}) 

301 exceeds non-zero C{f_max}. 

302 ''' 

303 f = Float_(f_max=f_max, low=_0_0, high=EPS1) 

304 if self._f_max != f: 

305 self._exactest(self.exact, self.ellipsoid, f) 

306 self._f_max = f 

307 

308 def _Inverse(self, ll1, ll2, wrap, **outmask): 

309 '''(INTERNAL) Short-cut version, see .latlonBase. 

310 ''' 

311 if wrap: 

312 ll2 = _unrollon(ll1, _Wrap.point(ll2)) 

313 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask) 

314 

315 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH_DISTANCE): 

316 '''Solve the I{inverse rhumb} problem. 

317 

318 @arg lat1: Latitude of the first point (C{degrees90}). 

319 @arg lon1: Longitude of the first point (C{degrees180}). 

320 @arg lat2: Latitude of the second point (C{degrees90}). 

321 @arg lon2: Longitude of the second point (C{degrees180}). 

322 

323 @return: L{GDict} with 4 to 8 items C{lat1, lon1, lat2, lon2, 

324 azi12, s12, a12, S12, }, the rhumb line's azimuth 

325 C{azi12} in compass C{degrees} between C{-180} and 

326 C{+180}, the rhumb distance C{s12} and rhumb angle 

327 C{a12} between both points in C{meter} respectively 

328 C{degrees} and the area C{S12} under the rhumb line 

329 in C{meter} I{squared}. 

330 

331 @raise ImportError: Package C{numpy} not found or not installed, 

332 only required for L{RhumbAux} area C{S12} 

333 when C{B{exact} is True}. 

334 

335 @note: The shortest rhumb line is found. If the end points are 

336 on opposite meridians, there are two shortest rhumb lines 

337 and the East-going one is chosen. 

338 

339 @note: If either point is a pole, the cosine of its latitude is 

340 taken to be C{sqrt(L{EPS})}. This position is extremely 

341 close to the actual pole and allows the calculation to be 

342 carried out in finite terms. 

343 ''' 

344 r = GDict(lat1=lat1, lon1=lon1, lat2=lat2, lon2=lon2, name=self.name) 

345 Cs = Caps 

346 if (outmask & Cs.AZIMUTH_DISTANCE_AREA): 

347 lon12, _ = _diff182(lon1, lon2, K_2_0=True) 

348 y, x, s1, s2 = self._Inverse4(lon12, r, outmask) 

349 if (outmask & Cs.AZIMUTH): 

350 r.set_(azi12=_atan2d(y, x)) 

351 if (outmask & Cs.AREA): 

352 S12 = self._S12d(s1, s2, lon12) 

353 r.set_(S12=unsigned0(S12)) # like .gx 

354 return r 

355 

356 def _Inverse4(self, lon12, r, outmask): # PYCHOK no cover 

357 '''I{Must be overloaded}, see function C{notOverloaded}. 

358 ''' 

359 _MODS.named.notOverloaded(self, lon12, r, outmask) 

360 

361 def Inverse8(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA): 

362 '''Like method L{Rhumb.Inverse} but returning a L{Rhumb8Tuple} with area C{S12}. 

363 ''' 

364 return self.Inverse(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple() 

365 

366 def _InverseLine(self, ll1, ll2, wrap, **caps_name): 

367 '''(INTERNAL) Short-cut version, see .latlonBase. 

368 ''' 

369 if wrap: 

370 ll2 = _unrollon(ll1, _Wrap.point(ll2)) 

371 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name) 

372 

373 def InverseLine(self, lat1, lon1, lat2, lon2, **caps_name): 

374 '''Define a C{RhumbLine} in terms of the I{inverse} rhumb problem. 

375 

376 @arg lat1: Latitude of the first point (C{degrees90}). 

377 @arg lon1: Longitude of the first point (C{degrees180}). 

378 @arg lat2: Latitude of the second point (C{degrees90}). 

379 @arg lon2: Longitude of the second point (C{degrees180}). 

380 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and 

381 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of 

382 L{Caps} values specifying the required capabilities. 

383 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}} 

384 should I{not} be reflected in this rhumb line. 

385 

386 @return: A C{RhumbLine...} instance and invoke its method 

387 C{ArcPosition} or C{Position} to compute points. 

388 

389 @note: Updates to this rhumb are reflected in the returned 

390 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}. 

391 ''' 

392 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH) 

393 return self._RhumbLine(self, lat1, lon1, r.azi12, **caps_name) 

394 

395 @Property_RO 

396 def _mpd(self): # PYCHOK no cover 

397 '''I{Must be overloaded}, see function C{notOverloaded}. 

398 ''' 

399 _MODS.named.notOverloaded(self) 

400 

401 @property_RO 

402 def _RhumbLine(self): # PYCHOK no cover 

403 '''I{Must be overloaded}, see function C{notOverloaded}. 

404 ''' 

405 _MODS.named.notOverloaded(self) 

406 

407 def _S12d(self, s1, s2, lon): # PYCHOK no cover 

408 '''I{Must be overloaded}, see function C{notOverloaded}. 

409 ''' 

410 _MODS.named.notOverloaded(self, s1, s2, lon) 

411 

412 @Property 

413 def TMorder(self): 

414 '''Get the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

415 ''' 

416 return self._mTM 

417 

418 @TMorder.setter # PYCHOK setter! 

419 def TMorder(self, order): 

420 '''Set the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

421 

422 @note: Setting C{TMorder} turns property C{exact} off, but only 

423 for L{Rhumb} instances. 

424 ''' 

425 m = _Xorder(_MODS.ktm._AlpCoeffs, RhumbError, TMorder=order) 

426 if self._mTM != m: 

427 _update_all_rls(self) 

428 self._mTM = m 

429 if self.exact and isinstance(self, _MODS.rhumbx.Rhumb): 

430 self.exact = False 

431 

432 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK no cover 

433 '''I{Must be overloaded}, see function C{notOverloaded}. 

434 ''' 

435 _MODS.named.notOverloaded(self, prec=prec, sep=sep) 

436 

437 

438class RhumbLineBase(_CapsBase): 

439 '''(INTERNAL) Base class for C{rhumbaux.RhumbLineAux} and C{rhumbx.RhumbLine}. 

440 ''' 

441 _azi12 = _0_0 

442 _calp = _1_0 

443# _caps = 0 

444# _debug = 0 

445# _lat1 = _0_0 

446# _lon1 = _0_0 

447# _lon12 = _0_0 

448 _Rhumb = RhumbBase # compatible C{Rhumb} class 

449 _rhumb = None # C{Rhumb} instance 

450 _salp = _0_0 

451 _talp = _0_0 

452 

453 def __init__(self, rhumb, lat1, lon1, azi12, caps=Caps.STANDARD, name=NN): 

454 '''New C{RhumbLine}. 

455 ''' 

456 _xinstanceof(self._Rhumb, rhumb=rhumb) 

457 

458 self._lat1 = _Lat(lat1=_fix90(lat1)) 

459 self._lon1 = _Lon(lon1= lon1) 

460 self._lon12 = _norm180(self._lon1) 

461 if azi12: # non-zero, non-None 

462 self.azi12 = _norm180(azi12) 

463 

464 n = name or rhumb.name 

465 if n: 

466 self.name=n 

467 

468 self._caps = caps 

469 self._debug |= (caps | rhumb._debug) & Caps._DEBUG_DIRECT_LINE 

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

471 self._rhumb = rhumb.copy(deep=False, name=_under(rhumb.name)) 

472 else: 

473 self._rhumb = rhumb 

474 _rls.append(self) 

475 

476 def __del__(self): # XXX use weakref? 

477 if _rls: # may be empty or None 

478 try: # PYCHOK no cover 

479 _rls.remove(self) 

480 except (TypeError, ValueError): 

481 pass 

482 self._rhumb = None 

483 # _update_all(self) # throws TypeError during Python 2 cleanup 

484 

485 def ArcPosition(self, a12, outmask=Caps.LATITUDE_LONGITUDE): 

486 '''Compute a point at a given angular distance on this rhumb line. 

487 

488 @arg a12: The angle along this rhumb line from its origin to the 

489 point (C{degrees}), can be negative. 

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

491 the quantities to be returned. 

492 

493 @return: L{GDict} with 4 to 8 items C{azi12, a12, s12, S12, lat2, 

494 lon2, lat1, lon1} with latitude C{lat2} and longitude 

495 C{lon2} of the point in C{degrees}, the rhumb distance 

496 C{s12} in C{meter} from the start point of and the area 

497 C{S12} under this rhumb line in C{meter} I{squared}. 

498 

499 @raise ImportError: Package C{numpy} not found or not installed, 

500 only required for L{RhumbLineAux} area C{S12} 

501 when C{B{exact} is True}. 

502 

503 @note: If B{C{a12}} is large enough that the rhumb line crosses a 

504 pole, the longitude of the second point is indeterminate and 

505 C{NAN} is returned for C{lon2} and area C{S12}. 

506 

507 If the first point is a pole, the cosine of its latitude is 

508 taken to be C{sqrt(L{EPS})}. This position is extremely 

509 close to the actual pole and allows the calculation to be 

510 carried out in finite terms. 

511 ''' 

512 return self._Position(a12, self.degrees2m(a12), outmask) 

513 

514 @Property 

515 def azi12(self): 

516 '''Get this rhumb line's I{azimuth} (compass C{degrees}). 

517 ''' 

518 return self._azi12 

519 

520 @azi12.setter # PYCHOK setter! 

521 def azi12(self, azi12): 

522 '''Set this rhumb line's I{azimuth} (compass C{degrees}). 

523 ''' 

524 z = _norm180(azi12) 

525 if self._azi12 != z: 

526 if self._rhumb: 

527 _update_all(self) 

528 self._azi12 = z 

529 self._salp, self._calp = t = sincos2d(z) # no NEG0 

530 self._talp = _over(*t) 

531 

532 @property_RO 

533 def azi12_sincos2(self): # PYCHOK no cover 

534 '''Get the sine and cosine of this rhumb line's I{azimuth} (2-tuple C{(sin, cos)}). 

535 ''' 

536 return self._scalp, self._calp 

537 

538 def degrees2m(self, angle): 

539 '''Convert an angular distance along this rhumb line to C{meter}. 

540 

541 @arg angle: Angular distance (C{degrees}). 

542 

543 @return: Distance (C{meter}). 

544 ''' 

545 return float(angle) * self.rhumb._mpd 

546 

547 @deprecated_method 

548 def distance2(self, lat, lon): # PYCHOK no cover 

549 '''DEPRECATED, use method L{RhumbLineAux.Inverse} or L{RhumbLine.Inverse}. 

550 

551 @return: A L{Distance2Tuple}C{(distance, initial)} with the C{distance} 

552 in C{meter} and C{initial} bearing (azimuth) in C{degrees}. 

553 ''' 

554 r = self.Inverse(lat, lon) 

555 return Distance2Tuple(r.s12, r.azi12) 

556 

557 @property_RO 

558 def ellipsoid(self): 

559 '''Get this rhumb line's ellipsoid (L{Ellipsoid}). 

560 ''' 

561 return self.rhumb.ellipsoid 

562 

563 @property_RO 

564 def exact(self): 

565 '''Get this rhumb line's I{exact} option (C{bool}). 

566 ''' 

567 return self.rhumb.exact 

568 

569 def intersection2(self, other, tol=_TOL, **eps): 

570 '''I{Iteratively} find the intersection of this and an other rhumb line. 

571 

572 @arg other: The other rhumb line (C{RhumbLine}). 

573 @kwarg tol: Tolerance for longitudinal convergence and parallel 

574 error (C{degrees}). 

575 @kwarg eps: Tolerance for L{pygeodesy.intersection3d3} (C{EPS}). 

576 

577 @return: A L{LatLon2Tuple}{(lat, lon)} with the C{lat}- and 

578 C{lon}gitude of the intersection point. 

579 

580 @raise IntersectionError: No convergence for this B{C{tol}} or 

581 no intersection for an other reason. 

582 

583 @see: Methods C{distance2} and C{nearestOn4} and function 

584 L{pygeodesy.intersection3d3}. 

585 

586 @note: Each iteration involves a round trip to this rhumb line's 

587 L{ExactTransverseMercator} or L{KTransverseMercator} 

588 projection and invoking function L{pygeodesy.intersection3d3} 

589 in that domain. 

590 ''' 

591 _xinstanceof(RhumbLineBase, other=other) 

592 _xdatum(self.rhumb, other.rhumb, Error=RhumbError) 

593 try: 

594 if other is self: 

595 raise ValueError(_coincident_) 

596 # make invariants and globals locals 

597 _s_3d, s_az = self._xTM3d, self.azi12 

598 _o_3d, o_az = other._xTM3d, other.azi12 

599 p = opposing(s_az, o_az, margin=tol) 

600 if p is not None: # == t in (True, False) 

601 raise ValueError(_anti_(_parallel_) if p else _parallel_) 

602 _diff = euclid # approximate length 

603 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3 

604 _LL2T = LatLon2Tuple 

605 _xTMr = self.xTM.reverse # ellipsoidal or spherical 

606 # use halfway point as initial estimate 

607 p = _LL2T(favg(self.lat1, other.lat1), 

608 favg(self.lon1, other.lon1)) 

609 for i in range(1, _TRIPS): 

610 v = _i3d3(_s_3d(p), s_az, # point + bearing 

611 _o_3d(p), o_az, useZ=False, **eps)[0] 

612 t = _xTMr(v.x, v.y, lon0=p.lon) # PYCHOK Reverse4Tuple 

613 d = _diff(t.lon - p.lon, t.lat) # PYCHOK t.lat + p.lat - p.lat 

614 p = _LL2T(t.lat + p.lat, t.lon) # PYCHOK t.lon + p.lon = lon0 

615 if d < tol: # 19 trips 

616 return _LL2T(p.lat, p.lon, iteration=i, # PYCHOK p... 

617 name=self.intersection2.__name__) 

618 except Exception as x: 

619 raise IntersectionError(_no_(_intersection_), cause=x) 

620 t = unstr(self.intersection2, other, **eps) 

621 raise IntersectionError(Fmt.no_convergence(d, tol), txt=t) 

622 

623 def Inverse(self, lat2, lon2, wrap=False): 

624 '''Return the rhumb angle, distance, azimuth, I{reverse} azimuth, etc. of 

625 a rhumb line between the given point and this rhumb line's start point. 

626 

627 @arg lat2: Latitude of the point (C{degrees}). 

628 @arg lon2: Longitude of the points (C{degrees}). 

629 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}} 

630 and B{C{lon2}} (C{bool}). 

631 

632 @return: L{GDict} with 8 items C{a12, s12, azi12, azi21, lat1, lon1, 

633 lat2, lon2}, the rhumb angle C{a12} and rhumb distance C{s12} 

634 between both points in C{degrees} respectively C{meter}, the 

635 rhumb line's azimuth C{azi12} and I{reverse} azimuth C{azi21} 

636 both in compass C{degrees} between C{-180} and C{+180}. 

637 ''' 

638 if wrap: 

639 _, lat2, lon2 = _Wrap.latlon3(self.lon1, _fix90(lat2), lon2) 

640 r = self.rhumb.Inverse(self.lat1, self.lon1, lat2, lon2) 

641 r.set_(azi21=_azireversed(r.azi12)) 

642 return r 

643 

644 @Property_RO 

645 def isLoxodrome(self): 

646 '''Is this rhumb line a loxodrome (C{True} if non-meridional and 

647 non-parallel, C{False} if parallel or C{None} if meridional)? 

648 

649 @see: I{Osborne's} U{2.5 Rumb lines and loxodromes 

650 <https://Zenodo.org/record/35392>}, page 37. 

651 ''' 

652 return bool(self._salp) if self._calp else None 

653 

654 @Property_RO 

655 def lat1(self): 

656 '''Get this rhumb line's latitude (C{degrees90}). 

657 ''' 

658 return self._lat1 

659 

660 @Property_RO 

661 def lon1(self): 

662 '''Get this rhumb line's longitude (C{degrees180}). 

663 ''' 

664 return self._lon1 

665 

666 @Property_RO 

667 def latlon1(self): 

668 '''Get this rhumb line's lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}). 

669 ''' 

670 return LatLon2Tuple(self.lat1, self.lon1) 

671 

672 def m2degrees(self, distance): 

673 '''Convert a distance along this rhumb line to an angular distance. 

674 

675 @arg distance: Distance (C{meter}). 

676 

677 @return: Angular distance (C{degrees}). 

678 ''' 

679 return _over(float(distance), self.rhumb._mpd) 

680 

681 @property_RO 

682 def _mu1(self): # PYCHOK no cover 

683 '''I{Must be overloaded}, see function C{notOverloaded}. 

684 ''' 

685 _MODS.named.notOverloaded(self) 

686 

687 def _mu2lat(self, mu2): # PYCHOK no cover 

688 '''I{Must be overloaded}, see function C{notOverloaded}. 

689 ''' 

690 _MODS.named.notOverloaded(self, mu2) 

691 

692 def nearestOn4(self, lat, lon, tol=_TOL, exact=None, eps=EPS, est=None): 

693 '''I{Iteratively} locate the point on this rhumb line nearest to the 

694 given point, in part transcoded from I{Karney}'s C++ U{rhumb-intercept 

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

696 

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

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

699 @kwarg tol: Longitudinal convergence tolerance (C{degrees}) or the 

700 distance tolerance (C(meter)) when C{B{exact} is None}, 

701 respectively C{not None}. 

702 @kwarg exact: If C{None}, use a rhumb line perpendicular to this rhumb 

703 line, otherwise use an I{exact} C{Geodesic...} from the 

704 given point perpendicular to this rhumb line (C{bool} or 

705 C{Geodesic...}), see method L{Ellipsoid.geodesic_}. 

706 @kwarg eps: Optional tolerance for L{pygeodesy.intersection3d3} 

707 (C{EPS}), used only if C{B{exact} is None}. 

708 @kwarg est: Optional, initial estimate for the distance C{s13} of 

709 the intersection I{along} this rhumb line (C{meter}), 

710 used only if C{B{exact} is not None}. 

711 

712 @return: A L{NearestOn4Tuple}C{(lat, lon, distance, normal)} with 

713 the C{lat}- and C{lon}gitude of the nearest point on and 

714 the C{distance} in C{meter} to this rhumb line and the 

715 C{normal} azimuth at the intersection. 

716 

717 @raise ImportError: I{Karney}'s U{geographiclib 

718 <https://PyPI.org/project/geographiclib>} 

719 package not found or not installed. 

720 

721 @raise IntersectionError: No convergence for this B{C{eps}} or 

722 no intersection for an other reason. 

723 

724 @see: Methods C{distance2} and C{intersection2} and function 

725 L{pygeodesy.intersection3d3}. 

726 ''' 

727 if exact is None: 

728 z = _norm180(self.azi12 + _90_0) # perpendicular azimuth 

729 rl = RhumbLineBase(self.rhumb, lat, lon, z, caps=Caps.LINE_OFF) 

730 p = self.intersection2(rl, tol=tol, eps=eps) 

731 r = rl.Inverse(p.lat, p.lon) 

732 r = NearestOn4Tuple(p.lat, p.lon, r.s12, z, 

733 iteration=p.iteration) 

734 else: # C{rhumb-intercept} 

735 azi = self.azi12 

736 szi = self._salp 

737 E = self.ellipsoid 

738 _gI = E.geodesic_(exact=exact).Inverse 

739 Cs = Caps 

740 gm = Cs._AZIMUTH_LATITUDE_LONGITUDE | Cs._REDUCEDLENGTH_GEODESICSCALE 

741 if est is None: # get an estimate from the perpendicular 

742 r = _gI(self.lat1, self.lon1, lat, lon, outmask=Cs.AZIMUTH_DISTANCE) 

743 d, _ = _diff182(r.azi2, azi, K_2_0=True) 

744 _, c = sincos2d(d) 

745 s12 = c * r.s12 # signed 

746 else: 

747 s12 = Meter(est=est) 

748 try: 

749 tol = Float_(tol=tol, low=EPS, high=None) 

750 # def _over(p, q): # see @note at RhumbLine[Aux].Position 

751 # return (p / (q or _copysign(tol, q))) if isfinite(q) else NAN 

752 

753 _S12 = Fsum(s12).fsum2_ 

754 for i in range(1, _TRIPS): # suffix 1 == C++ 2, 2 == C++ 3 

755 p = self.Position(s12) # outmask=Cs.LATITUDE_LONGITUDE 

756 r = _gI(lat, lon, p.lat2, p.lon2, outmask=gm) 

757 d, _ = _diff182(azi, r.azi2, K_2_0=True) 

758 s, c, s2, c2 = sincos2d_(d, r.lat2) 

759 c2 *= E.rocPrimeVertical(r.lat2) # aka rocTransverse 

760 s *= _over(s2 * szi, c2) - _over(s * r.M21, r.m12) 

761 s12, t = _S12(c / s) # XXX _over? 

762 if fabs(t) < tol: # or fabs(c) < EPS 

763 break 

764 r = NearestOn4Tuple(r.lat2, r.lon2, s12, r.azi2, 

765 iteration=i) 

766 except Exception as x: # Fsum(NAN) Value-, ZeroDivisionError 

767 t = unstr(self.nearestOn4, lat, lon, tol=tol, exact=exact, 

768 iteration=i, eps=eps, est=est) 

769 raise IntersectionError(t, cause=x) 

770 return r 

771 

772 def Position(self, s12, outmask=Caps.LATITUDE_LONGITUDE): 

773 '''Compute a point at a given distance on this rhumb line. 

774 

775 @arg s12: The distance along this rhumb line from its origin to 

776 the point (C{meters}), can be negative. 

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

778 the quantities to be returned. 

779 

780 @return: L{GDict} with 4 to 8 items C{azi12, a12, s12, S12, lat2, 

781 lon2, lat1, lon1} with latitude C{lat2} and longitude 

782 C{lon2} of the point in C{degrees}, the rhumb angle C{a12} 

783 in C{degrees} from the start point of and the area C{S12} 

784 under this rhumb line in C{meter} I{squared}. 

785 

786 @raise ImportError: Package C{numpy} not found or not installed, 

787 only required for L{RhumbLineAux} area C{S12} 

788 when C{B{exact} is True}. 

789 

790 @note: If B{C{s12}} is large enough that the rhumb line crosses a 

791 pole, the longitude of the second point is indeterminate and 

792 C{NAN} is returned for C{lon2} and area C{S12}. 

793 

794 If the first point is a pole, the cosine of its latitude is 

795 taken to be C{sqrt(L{EPS})}. This position is extremely 

796 close to the actual pole and allows the calculation to be 

797 carried out in finite terms. 

798 ''' 

799 return self._Position(self.m2degrees(s12), s12, outmask) 

800 

801 def _Position(self, a12, s12, outmask): 

802 '''(INTERNAL) C{Arc-/Position} helper. 

803 ''' 

804 r = GDict(azi12=self.azi12, a12=a12, s12=s12, name=self.name) 

805 Cs = Caps 

806 if (outmask & Cs.LATITUDE_LONGITUDE_AREA): 

807 if a12 or s12: 

808 mu12 = self._calp * a12 

809 mu2 = self._mu1 + mu12 

810 if fabs(mu2) > 90: # past pole 

811 mu2 = _norm180(mu2) # reduce to [-180, 180) 

812 if fabs(mu2) > 90: # point on anti-meridian 

813 mu2 = _norm180(_loneg(mu2)) 

814 lat2 = self._mu2lat(mu2) 

815 lon2 = S12 = NAN 

816 else: 

817 lat2, lon2, S1, S2 = self._Position4(a12, mu2, s12, mu12) 

818 if (outmask & Cs.AREA): 

819 S12 = self.rhumb._S12d(S1, S2, lon2) 

820 S12 = unsigned0(S12) # like .gx 

821# else: 

822# S12 = None # unused 

823 if (outmask & Cs.LONGITUDE): 

824 if (outmask & Cs.LONG_UNROLL): 

825 lon2 += self.lon1 

826 else: 

827 lon2 = _norm180(self._lon12 + lon2) 

828 else: # coincident or outmask 0 

829 lat2, lon2 = self.latlon1 

830 S12 = _0_0 

831 

832 if (outmask & Cs.LATITUDE): 

833 r.set_(lat2=lat2, lat1=self.lat1) 

834 if (outmask & Cs.LONGITUDE): 

835 r.set_(lon2=lon2, lon1=self.lon1) 

836 if (outmask & Cs.AREA): 

837 r.set_(S12=S12) 

838 return r 

839 

840 def _Position4(self, a12, mu2, s12, mu12): # PYCHOK no cover 

841 '''I{Must be overloaded}, see function C{notOverloaded}. 

842 ''' 

843 _MODS.named.notOverloaded(self, a12, s12, mu2, mu12) 

844 

845 @Property_RO 

846 def rhumb(self): 

847 '''Get this rhumb line's rhumb (L{RhumbAux} or L{Rhumb}). 

848 ''' 

849 return self._rhumb 

850 

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

852 '''Return this C{RhumbLine} as string. 

853 

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

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

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

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

858 

859 @return: C{RhumbLine} (C{str}). 

860 ''' 

861 d = dict(rhumb=self.rhumb, lat1=self.lat1, lon1=self.lon1, 

862 azi12=self.azi12, exact=self.exact, 

863 TMorder=self.TMorder, xTM=self.xTM) 

864 return sep.join(pairs(itemsorted(d, asorted=False), prec=prec)) 

865 

866 @property_RO 

867 def TMorder(self): 

868 '''Get this rhumb line's I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

869 ''' 

870 return self.rhumb.TMorder 

871 

872 @Property_RO 

873 def xTM(self): 

874 '''Get this rhumb line's I{Transverse Mercator} projection (L{ExactTransverseMercator} 

875 if I{exact} and I{ellipsoidal}, otherwise L{KTransverseMercator} for C{TMorder}). 

876 ''' 

877 E = self.ellipsoid 

878 # ExactTransverseMercator doesn't handle spherical earth models 

879 return _MODS.etm.ExactTransverseMercator(E) if self.exact and E.isEllipsoidal else \ 

880 _MODS.ktm.KTransverseMercator(E, TMorder=self.TMorder) 

881 

882 def _xTM3d(self, latlon0, z=INT0, V3d=Vector3d): 

883 '''(INTERNAL) C{xTM.forward} this C{latlon1} to C{V3d} with B{C{latlon0}} 

884 as current intersection estimate and central meridian. 

885 ''' 

886 t = self.xTM.forward(self.lat1 - latlon0.lat, self.lon1, lon0=latlon0.lon) 

887 return V3d(t.easting, t.northing, z) 

888 

889 

890__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase) 

891 

892if __name__ == '__main__': 

893 

894 from pygeodesy import printf, Rhumb as R, RhumbAux as A 

895 

896 A = A(_EWGS84).Line(30, 0, 45) 

897 R = R(_EWGS84).Line(30, 0, 45) 

898 

899 for i in range(1, 10): 

900 s = .5e6 + 1e6 / i 

901 a = A.Position(s).lon2 

902 r = R.Position(s).lon2 

903 e = (fabs(a - r) / a) if a else 0 

904 printf('Position.lon2 %.14f vs %.14f, diff %g', r, a, e) 

905 

906 for exact in (None, False, True): 

907 for est in (None, 1e6): 

908 a = A.nearestOn4(60, 0, exact=exact, est=est) 

909 r = R.nearestOn4(60, 0, exact=exact, est=est) 

910 printf('%s, iteration=%s, exact=%s, est=%s\n%s, iteration=%s', 

911 a.toRepr(), a.iteration, exact, est, 

912 r.toRepr(), r.iteration, nl=1) 

913 

914# % python3 -m pygeodesy.rhumbBase 

915 

916# Position.lon2 11.61455846901637 vs 11.61455846901637, diff 1.52942e-16 

917# Position.lon2 7.58982302826842 vs 7.58982302826842, diff 2.34045e-16 

918# Position.lon2 6.28526067416369 vs 6.28526067416369, diff 2.82623e-16 

919# Position.lon2 5.63938995325146 vs 5.63938995325146, diff 1.57495e-16 

920# Position.lon2 5.25385527435707 vs 5.25385527435707, diff 0 

921# Position.lon2 4.99764604290380 vs 4.99764604290380, diff 8.88597e-16 

922# Position.lon2 4.81503363740473 vs 4.81503363740473, diff 1.84459e-16 

923# Position.lon2 4.67828821748836 vs 4.67828821748835, diff 5.69553e-16 

924# Position.lon2 4.57205667906283 vs 4.57205667906283, diff 1.94262e-16 

925 

926# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=None 

927# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=3278714.911694, normal=135.0), iteration=9 

928 

929# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9, exact=None, est=1000000.0 

930# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=3278714.911694, normal=135.0), iteration=9 

931 

932# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=False, est=None 

933# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5 

934 

935# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=False, est=1000000.0 

936# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7 

937 

938# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5, exact=True, est=None 

939# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=5 

940 

941# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7, exact=True, est=1000000.0 

942# NearestOn4Tuple(lat=49.634582, lon=25.767876, distance=3083112.636236, normal=135.0), iteration=7 

943 

944# **) MIT License 

945# 

946# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

947# 

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

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

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

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

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

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

954# 

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

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

957# 

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

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

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

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

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

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

964# OTHER DEALINGS IN THE SOFTWARE.