Coverage for pygeodesy/rhumbBase.py: 96%

248 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-07 07:28 -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. 

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:Charles@Karney.com>} (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 _xinstanceof # from .karney 

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

28 _0_0, _0_01, _1_0, _90_0, _180_0 

29# from pygeodesy.datums import _spherical_datum # _MODS 

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

31 _xdatum, _xkwds, _Xorder 

32# from pygeodesy.etm import ExactTransverseMercator # _MODS 

33from pygeodesy.fmath import euclid, favg, fabs 

34# from pygeodesy.fsums import Fsum # _MODS 

35from pygeodesy.interns import NN, _coincident_, _COMMASPACE_, _intersection_, \ 

36 _no_, _under 

37from pygeodesy.karney import Caps, _CapsBase, _diff182, _EWGS84, _fix90, \ 

38 _norm180, _xinstanceof 

39from pygeodesy.ktm import KTransverseMercator, _AlpCoeffs # PYCHOK used! 

40from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS 

41# from pygeodesy.named import notOverloaded # _MODS 

42from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple 

43from pygeodesy.props import Property, Property_RO, property_RO, _update_all 

44from pygeodesy.auxilats.auxily import _over 

45from pygeodesy.streprs import Fmt, pairs, unstr 

46from pygeodesy.units import Float_, Lat, Lon, Meter 

47from pygeodesy.utily import sincos2d, sincos2d_, _unrollon, _Wrap 

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

49 

50# from math import fabs # from .fmath 

51 

52__all__ = () 

53__version__ = '23.08.05' 

54 

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

56_TRIPS = 65 # .intersection2, .nearestOn4, 18+ 

57 

58 

59class _Lat(Lat): 

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

61 ''' 

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

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

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

65 

66 

67class _Lon(Lon): 

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

69 ''' 

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

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

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

73 

74 

75def _update_all_rls(r): 

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

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

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

79 ''' 

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

81 _update_all(r) 

82 for rl in _rls: # PYCHOK use weakref? 

83 if rl._rhumb is r: 

84 _update_all(rl) 

85 

86 

87class RhumbBase(_CapsBase): 

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

89 ''' 

90 _E = _EWGS84 

91 _exact = True 

92 _f_max = _0_01 

93 _mTM = 6 # see .TMorder 

94 

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

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

97 ''' 

98 if f is not None: 

99 self.ellipsoid = a_earth, f 

100 elif a_earth not in (_EWGS84, None): 

101 self.ellipsoid = a_earth 

102 if not exact: 

103 self.exact = False 

104 if name: 

105 self.name = name 

106 

107 @Property_RO 

108 def a(self): 

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

110 ''' 

111 return self.ellipsoid.a 

112 

113 equatoradius = a 

114 

115 @Property_RO 

116 def b(self): 

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

118 ''' 

119 return self.ellipsoid.b 

120 

121 polaradius = b 

122 

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

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

125 ''' 

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

127 

128 def Direct(self, lat1, lon1, azi12, s12, outmask=0): # PYCHOK no cover 

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

130 ''' 

131 _MODS.named.notOverloaded(self, lat1, lon1, azi12, s12, outmask=outmask) 

132 

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

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

135 ''' 

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

137 

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

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

140 ''' 

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

142 

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

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

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

146 

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

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

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

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

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

152 L{Caps} values specifying the required capabilities. 

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

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

155 

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

157 C{.Position} to compute each point. 

158 

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

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

161 ''' 

162 return self._RhumbLine(self, lat1=lat1, lon1=lon1, azi12=azi12, 

163 **caps_name) 

164 

165 Line = DirectLine # synonyms 

166 

167 @Property 

168 def ellipsoid(self): 

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

170 ''' 

171 return self._E 

172 

173 @ellipsoid.setter # PYCHOK setter! 

174 def ellipsoid(self, a_earth_f): 

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

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

177 

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

179 ''' 

180 E = _MODS.datums._spherical_datum(a_earth_f, Error=RhumbError).ellipsoid 

181 if self._E != E: 

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

183 _update_all_rls(self) 

184 self._E = E 

185 

186 @Property 

187 def exact(self): 

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

189 ''' 

190 return self._exact 

191 

192 @exact.setter # PYCHOK setter! 

193 def exact(self, exact): 

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

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

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

197 

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

199 exceeds non-zero C{f_max}. 

200 

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

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

203 ''' 

204 x = bool(exact) 

205 if self._exact != x: 

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

207 _update_all_rls(self) 

208 self._exact = x 

209 

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

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

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

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

214 

215 @Property_RO 

216 def f(self): 

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

218 ''' 

219 return self.ellipsoid.f 

220 

221 flattening = f 

222 

223 @property 

224 def f_max(self): 

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

226 ''' 

227 return self._f_max 

228 

229 @f_max.setter # PYCHOK setter! 

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

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

232 

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

234 exceeds non-zero C{f_max}. 

235 ''' 

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

237 if self._f_max != f: 

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

239 self._f_max = f 

240 

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

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

243 ''' 

244 if wrap: 

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

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

247 

248 def Inverse(self, lat1, lon1, lat2, lon2, outmask=0): # PYCHOK no cover 

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

250 ''' 

251 _MODS.named.notOverloaded(self, lat1, lon1, lat2, lon2, outmask=outmask) 

252 

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

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

255 ''' 

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

257 

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

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

260 ''' 

261 if wrap: 

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

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

264 

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

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

267 

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

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

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

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

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

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

274 L{Caps} values specifying the required capabilities. 

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

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

277 

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

279 C{.Position} to compute each point. 

280 

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

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

283 ''' 

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

285 return self._RhumbLine(self, lat1=lat1, lon1=lon1, azi12=r.azi12, 

286 **caps_name) 

287 

288 @Property_RO 

289 def _RhumbLine(self): # PYCHOK no cover 

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

291 ''' 

292 _MODS.named.notOverloaded(self) 

293 

294 def _TMorder(self, order): 

295 '''(INTERNAL) Set the I{Transverse Mercator} order (C{int}). 

296 

297 @note: Setting C{TMorder} turns property C{exact} off. 

298 ''' 

299 x = self.exact 

300 n = _Xorder(_MODS.ktm._AlpCoeffs, RhumbError, TMorder=order) 

301 if self._mTM != n: 

302 _update_all_rls(self) 

303 self._mTM = n 

304 x = False 

305 return x 

306 

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

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

309 ''' 

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

311 

312 

313class RhumbLineBase(RhumbBase): 

314 '''(INTERNAL) Class C{RhumbLine} 

315 ''' 

316 _azi12 = _0_0 

317# _caps = 0 

318# _debug = 0 

319# _lat1 = _0_0 

320# _lon1 = _0_0 

321# _lon12 = _0_0 

322 _salp = _0_0 

323 _calp = _1_0 

324 _Rhumb = RhumbBase # compatible C{Rhumb} class 

325 _rhumb = None # C{Rhumb} instance 

326 

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

328 '''New C{RhumbLine}. 

329 ''' 

330 _xinstanceof(self._Rhumb, rhumb=rhumb) 

331 

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

333 self._lon1 = _Lon(lon1= lon1) 

334 self._lon12 = _norm180(self._lon1) 

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

336 self.azi12 = _norm180(azi12) 

337 

338 n = name or rhumb.name 

339 if n: 

340 self.name=n 

341 

342 self._caps = caps 

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

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

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

346 else: 

347 self._rhumb = rhumb 

348 _rls.append(self) 

349 

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

351 if _rls: # may be empty or None 

352 try: # PYCHOK no cover 

353 _rls.remove(self) 

354 except (TypeError, ValueError): 

355 pass 

356 self._rhumb = None 

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

358 

359 @Property 

360 def azi12(self): 

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

362 ''' 

363 return self._azi12 

364 

365 @azi12.setter # PYCHOK setter! 

366 def azi12(self, azi12): 

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

368 ''' 

369 z = _norm180(azi12) 

370 if self._azi12 != z: 

371 if self._rhumb: 

372 _update_all(self) 

373 self._azi12 = z 

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

375 

376 @Property_RO 

377 def azi12_sincos2(self): # PYCHOK no cover 

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

379 ''' 

380 return self._scalp, self._calp 

381 

382 def distance2(self, lat, lon): 

383 '''Return the distance from and (initial) bearing at the given 

384 point to this rhumb line's start point. 

385 

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

387 @arg lon: Longitude of the points (C{degrees}). 

388 

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

390 in C{meter} and C{initial} bearing in C{degrees}. 

391 

392 @see: Methods C{intersection2} and C{nearestOn4} of L{RhumbLine} and 

393 L{RhumbLineAux}. 

394 ''' 

395 r = self.rhumb.Inverse(self.lat1, self.lon1, lat, lon) 

396# outmask=Caps.AZIMUTH_DISTANCE) 

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

398 

399 @property_RO 

400 def ellipsoid(self): 

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

402 ''' 

403 return self.rhumb.ellipsoid 

404 

405 @property_RO 

406 def exact(self): 

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

408 ''' 

409 return self.rhumb.exact 

410 

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

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

413 

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

415 @kwarg tol: Tolerance for longitudinal convergence (C{degrees}). 

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

417 

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

419 C{lon}gitude of the intersection point. 

420 

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

422 no intersection for an other reason. 

423 

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

425 L{pygeodesy.intersection3d3}. 

426 

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

428 L{ExactTransverseMercator} or L{KTransverseMercator} 

429 projection and invoking function L{pygeodesy.intersection3d3} 

430 in that domain. 

431 ''' 

432 _xinstanceof(RhumbLineBase, other=other) 

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

434 try: 

435 if other is self: 

436 raise ValueError(_coincident_) 

437 # make globals and invariants locals 

438 _diff = euclid # approximate length 

439 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3 

440 _LL2T = LatLon2Tuple 

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

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

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

444 # use halfway point as initial estimate 

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

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

447 for i in range(1, _TRIPS): 

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

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

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

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

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

453 if d < tol: 

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

455 name=self.intersection2.__name__) 

456 except Exception as x: 

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

458 t = unstr(self.intersection2, tol=tol, **eps) 

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

460 

461 @Property_RO 

462 def lat1(self): 

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

464 ''' 

465 return self._lat1 

466 

467 @Property_RO 

468 def lon1(self): 

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

470 ''' 

471 return self._lon1 

472 

473 @Property_RO 

474 def latlon1(self): 

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

476 ''' 

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

478 

479 def _mu22(self, mu12, mu1): 

480 mu2 = mu12 + mu1 

481 x90 = fabs(mu2) > 90 

482 if x90: # reduce to [-180, 180) 

483 mu2 = _norm180(mu2) 

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

485 mu2 = _norm180(_180_0 - mu2) 

486 return mu2, x90 

487 

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

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

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

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

492 

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

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

495 @kwarg tol: Longitudinal convergence tolerance (C{degrees}) only if 

496 C{B{exact} is None}, otherwise the distance tolerance 

497 (C(meter)) and when C{B{exact} is not None}. 

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

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

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

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

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

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

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

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

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

507 

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

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

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

511 azimuth of the C{normal}, perpendicular to this rhumb line. 

512 

513 @raise ImportError: I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>} 

514 package not found or not installed iff C{B{exact} is False}. 

515 

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

517 no intersection for an other reason. 

518 

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

520 L{pygeodesy.intersection3d3}. 

521 ''' 

522 if exact is None: 

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

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

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

526 t = rl.distance2(p.lat, p.lon) 

527 r = NearestOn4Tuple(p.lat, p.lon, t.distance, z, 

528 iteration=p.iteration) 

529 else: # C{rhumb-intercept} 

530 azi = self.azi12 

531 szi = self._salp 

532 E = self.ellipsoid 

533 gX = E.geodesic_(exact=exact) 

534 Cs = Caps 

535 m = Cs.STANDARD | Cs.REDUCEDLENGTH | Cs.GEODESICSCALE 

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

537 r = gX.Inverse(self.lat1, self.lon1, lat, lon, outmask=Cs.AZIMUTH_DISTANCE) 

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

539 _, c = sincos2d(d) 

540 s12 = c * r.s12 

541 else: 

542 s12 = Meter(est=est) 

543 _s12 = _MODS.fsums.Fsum(s12).fsum2_ 

544 # suffix 1 == C++2, 2 == C++3, 3 == C++2 

545 for i in range(1, _TRIPS): 

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

547 r = gX.Inverse(lat, lon, p.lat2, p.lon2, outmask=m) 

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

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

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

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

552 s12, t = _s12(_over(c, s)) 

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

554 break 

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

556 iteration=i) 

557 return r 

558 

559 def Position(self, s12, outmask=0): # PYCHOK no cover 

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

561 ''' 

562 _MODS.named.notOverloaded(self, s12, outmask=outmask) 

563 

564 @Property_RO 

565 def rhumb(self): 

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

567 ''' 

568 return self._rhumb 

569 

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

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

572 

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

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

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

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

577 

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

579 ''' 

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

581 azi12=self.azi12, exact=self.exact, 

582 TMorder=self.TMorder, xTM=self.xTM) 

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

584 

585 @property_RO 

586 def TMorder(self): 

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

588 ''' 

589 return self.rhumb.TMorder 

590 

591 @Property_RO 

592 def xTM(self): 

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

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

595 ''' 

596 E = self.ellipsoid 

597 # ExactTransverseMercator doesn't handle spherical earth models 

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

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

600 

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

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

603 as current intersection estimate and central meridian. 

604 ''' 

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

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

607 

608 

609__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase) 

610 

611 

612if __name__ == '__main__': 

613 import sys 

614 

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

616 

617 if 'aux' in sys.argv: 

618 r = A(_EWGS84).Line(30, 0, 45) 

619 elif 'x' in sys.argv: 

620 r = R(_EWGS84).Line(30, 0, 45) 

621 else: 

622 sys.exit('aux or x') 

623 

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

625 for est in (None, 1e6): 

626 p = r.nearestOn4(60, 0, exact=exact, est=est) 

627 printf('%s, exact=%s, est=%s, iteration=%s', 

628 p.toRepr(), exact, est, p.iteration) 

629 

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

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

632 n = 1 

633 for i in range(1, 10): 

634 s = .5e6 + 1e6 / i 

635 a = A.Position(s).lon2 

636 r = R.Position(s).lon2 

637 e = (fabs(a - r) * 100.0 / a) if a else 0 

638 printf('Positions %s vs %s, diff %.2f%%', r, a, e, nl=n) 

639 n = 0 

640 

641# % python3 -m pygeodesy.rhumbBase x 

642 

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

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

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

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

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

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

649 

650# Positions 11.614558469016366 vs 11.61455846901637, diff 0.00% 

651# Positions 7.589823028268422 vs 7.589823028268423, diff 0.00% 

652# Positions 6.285260674163688 vs 6.285260674163687, diff 0.00% 

653# Positions 5.639389953251457 vs 5.6393899532514595, diff 0.00% 

654# Positions 5.253855274357075 vs 5.253855274357073, diff 0.00% 

655# Positions 4.9976460429038 vs 4.9976460429038045, diff 0.00% 

656# Positions 4.815033637404729 vs 4.81503363740473, diff 0.00% 

657# Positions 4.678288217488354 vs 4.678288217488353, diff 0.00% 

658# Positions 4.572056679062825 vs 4.572056679062826, diff 0.00% 

659 

660 

661# % python3 -m pygeodesy.rhumbBase aux 

662 

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

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

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

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

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

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

669 

670# Positions 11.614558469016366 vs 11.61455846901637, diff 0.00% 

671# Positions 7.589823028268422 vs 7.589823028268423, diff 0.00% 

672# Positions 6.285260674163688 vs 6.285260674163687, diff 0.00% 

673# Positions 5.639389953251457 vs 5.6393899532514595, diff 0.00% 

674# Positions 5.253855274357075 vs 5.253855274357073, diff 0.00% 

675# Positions 4.9976460429038 vs 4.9976460429038045, diff 0.00% 

676# Positions 4.815033637404729 vs 4.81503363740473, diff 0.00% 

677# Positions 4.678288217488354 vs 4.678288217488353, diff 0.00% 

678# Positions 4.572056679062825 vs 4.572056679062826, diff 0.00% 

679 

680# **) MIT License 

681# 

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

683# 

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

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

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

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

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

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

690# 

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

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

693# 

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

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

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

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

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

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

700# OTHER DEALINGS IN THE SOFTWARE.