Coverage for pygeodesy/rhumbBase.py: 96%

254 statements  

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

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

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

28 _0_0, _0_01, _1_0, _90_0 

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_, _intersection_, \ 

38 _no_, _parallel_, _under 

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

40 _EWGS84, _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 Property, Property_RO, property_RO, _update_all 

46from pygeodesy.streprs import Fmt, pairs, unstr 

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

48from pygeodesy.utily import _loneg, sincos2d, sincos2d_, _unrollon, _Wrap, \ 

49 sincos2_ # PYCHOK shared 

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

51 

52# from math import fabs # from .fmath 

53 

54__all__ = () 

55__version__ = '23.09.15' 

56 

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

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

59 

60 

61class _Lat(Lat): 

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

63 ''' 

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

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

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

67 

68 

69class _Lon(Lon): 

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

71 ''' 

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

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

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

75 

76 

77def _update_all_rls(r): 

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

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

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

81 ''' 

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

83 _update_all(r) 

84 for rl in _rls: # PYCHOK use weakref? 

85 if rl._rhumb is r: 

86 _update_all(rl) 

87 

88 

89class RhumbBase(_CapsBase): 

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

91 ''' 

92 _E = _EWGS84 

93 _exact = True 

94 _f_max = _0_01 

95 _mTM = 6 # see .TMorder 

96 

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

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

99 ''' 

100 if f is not None: 

101 self.ellipsoid = a_earth, f 

102 elif a_earth not in (_EWGS84, None): 

103 self.ellipsoid = a_earth 

104 if not exact: 

105 self.exact = False 

106 if name: 

107 self.name = name 

108 

109 @Property_RO 

110 def a(self): 

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

112 ''' 

113 return self.ellipsoid.a 

114 

115 equatoradius = a 

116 

117 @Property_RO 

118 def b(self): 

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

120 ''' 

121 return self.ellipsoid.b 

122 

123 polaradius = b 

124 

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

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

127 ''' 

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

129 

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

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

132 ''' 

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

134 

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

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

137 ''' 

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

139 

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

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

142 ''' 

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

144 

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

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

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

148 

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

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

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

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

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

154 L{Caps} values specifying the required capabilities. 

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

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

157 

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

159 C{.Position} to compute each point. 

160 

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

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

163 ''' 

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

165 **caps_name) 

166 

167 Line = DirectLine # synonyms 

168 

169 @Property 

170 def ellipsoid(self): 

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

172 ''' 

173 return self._E 

174 

175 @ellipsoid.setter # PYCHOK setter! 

176 def ellipsoid(self, a_earth_f): 

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

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

179 

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

181 ''' 

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

183 if self._E != E: 

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

185 _update_all_rls(self) 

186 self._E = E 

187 

188 @Property 

189 def exact(self): 

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

191 ''' 

192 return self._exact 

193 

194 @exact.setter # PYCHOK setter! 

195 def exact(self, exact): 

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

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

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

199 

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

201 exceeds non-zero C{f_max}. 

202 

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

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

205 ''' 

206 x = bool(exact) 

207 if self._exact != x: 

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

209 _update_all_rls(self) 

210 self._exact = x 

211 

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

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

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

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

216 

217 @Property_RO 

218 def f(self): 

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

220 ''' 

221 return self.ellipsoid.f 

222 

223 flattening = f 

224 

225 @property 

226 def f_max(self): 

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

228 ''' 

229 return self._f_max 

230 

231 @f_max.setter # PYCHOK setter! 

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

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

234 

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

236 exceeds non-zero C{f_max}. 

237 ''' 

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

239 if self._f_max != f: 

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

241 self._f_max = f 

242 

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

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

245 ''' 

246 if wrap: 

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

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

249 

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

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

252 ''' 

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

254 

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

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

257 ''' 

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

259 

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

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

262 ''' 

263 if wrap: 

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

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

266 

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

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

269 

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

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

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

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

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

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

276 L{Caps} values specifying the required capabilities. 

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

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

279 

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

281 C{.Position} to compute each point. 

282 

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

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

285 ''' 

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

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

288 **caps_name) 

289 

290 @Property_RO 

291 def _RhumbLine(self): # PYCHOK no cover 

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

293 ''' 

294 _MODS.named.notOverloaded(self) 

295 

296 def _TMorder(self, order): 

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

298 

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

300 ''' 

301 x = self.exact 

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

303 if self._mTM != n: 

304 _update_all_rls(self) 

305 self._mTM = n 

306 x = False 

307 return x 

308 

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

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

311 ''' 

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

313 

314 

315class RhumbLineBase(RhumbBase): 

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

317 ''' 

318 _azi12 = _0_0 

319# _caps = 0 

320# _debug = 0 

321# _lat1 = _0_0 

322# _lon1 = _0_0 

323# _lon12 = _0_0 

324 _salp = _0_0 

325 _calp = _1_0 

326 _Rhumb = RhumbBase # compatible C{Rhumb} class 

327 _rhumb = None # C{Rhumb} instance 

328 

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

330 '''New C{RhumbLine}. 

331 ''' 

332 _xinstanceof(self._Rhumb, rhumb=rhumb) 

333 

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

335 self._lon1 = _Lon(lon1= lon1) 

336 self._lon12 = _norm180(self._lon1) 

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

338 self.azi12 = _norm180(azi12) 

339 

340 n = name or rhumb.name 

341 if n: 

342 self.name=n 

343 

344 self._caps = caps 

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

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

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

348 else: 

349 self._rhumb = rhumb 

350 _rls.append(self) 

351 

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

353 if _rls: # may be empty or None 

354 try: # PYCHOK no cover 

355 _rls.remove(self) 

356 except (TypeError, ValueError): 

357 pass 

358 self._rhumb = None 

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

360 

361 @Property 

362 def azi12(self): 

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

364 ''' 

365 return self._azi12 

366 

367 @azi12.setter # PYCHOK setter! 

368 def azi12(self, azi12): 

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

370 ''' 

371 z = _norm180(azi12) 

372 if self._azi12 != z: 

373 if self._rhumb: 

374 _update_all(self) 

375 self._azi12 = z 

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

377 

378 @property_RO 

379 def azi12_sincos2(self): # PYCHOK no cover 

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

381 ''' 

382 return self._scalp, self._calp 

383 

384 def distance2(self, lat, lon): 

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

386 point to this rhumb line's start point. 

387 

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

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

390 

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

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

393 

394 @see: Methods C{intersection2} and C{nearestOn4} of L{RhumbLineAux} and 

395 L{RhumbLine}. 

396 ''' 

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

398# outmask=Caps.AZIMUTH_DISTANCE) 

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

400 

401 @property_RO 

402 def ellipsoid(self): 

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

404 ''' 

405 return self.rhumb.ellipsoid 

406 

407 @property_RO 

408 def exact(self): 

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

410 ''' 

411 return self.rhumb.exact 

412 

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

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

415 

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

417 @kwarg tol: Tolerance for longitudinal convergence and parallel 

418 error (C{degrees}). 

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

420 

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

422 C{lon}gitude of the intersection point. 

423 

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

425 no intersection for an other reason. 

426 

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

428 L{pygeodesy.intersection3d3}. 

429 

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

431 L{ExactTransverseMercator} or L{KTransverseMercator} 

432 projection and invoking function L{pygeodesy.intersection3d3} 

433 in that domain. 

434 ''' 

435 _xinstanceof(RhumbLineBase, other=other) 

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

437 try: 

438 if other is self: 

439 raise ValueError(_coincident_) 

440 # make invariants and globals locals 

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

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

443 t = opposing(s_az, o_az, margin=tol) 

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

445 raise ValueError(_parallel_) 

446 _diff = euclid # approximate length 

447 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3 

448 _LL2T = LatLon2Tuple 

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

450 # use halfway point as initial estimate 

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

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

453 for i in range(1, _TRIPS): 

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

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

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

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

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

459 if d < tol: # 19 trips 

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

461 name=self.intersection2.__name__) 

462 except Exception as x: 

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

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

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

466 

467 @Property_RO 

468 def isLoxodrome(self): 

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

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

471 

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

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

474 ''' 

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

476 

477 @Property_RO 

478 def lat1(self): 

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

480 ''' 

481 return self._lat1 

482 

483 @Property_RO 

484 def lon1(self): 

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

486 ''' 

487 return self._lon1 

488 

489 @Property_RO 

490 def latlon1(self): 

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

492 ''' 

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

494 

495 def _mu22(self, mu12, mu1): 

496 mu2 = mu12 + mu1 

497 x90 = fabs(mu2) > 90 

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

499 mu2 = _norm180(mu2) 

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

501 mu2 = _norm180(_loneg(mu2)) 

502 return mu2, x90 

503 

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

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

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

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

508 

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

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

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

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

513 respectively C{not None}. 

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

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

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

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

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

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

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

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

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

523 

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

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

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

527 C{normal} azimuth at the intersection. 

528 

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

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

531 package not found or not installed. 

532 

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

534 no intersection for an other reason. 

535 

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

537 L{pygeodesy.intersection3d3}. 

538 ''' 

539 if exact is None: 

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

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

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

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

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

545 iteration=p.iteration) 

546 else: # C{rhumb-intercept} 

547 azi = self.azi12 

548 szi = self._salp 

549 E = self.ellipsoid 

550 gX = E.geodesic_(exact=exact) 

551 Cs = Caps 

552 gm = Cs.AZIMUTH | Cs.LATITUDE_LONGITUDE | Cs.REDUCEDLENGTH | Cs.GEODESICSCALE 

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

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

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

556 _, c = sincos2d(d) 

557 s12 = c * r.s12 # signed 

558 else: 

559 s12 = Meter(est=est) 

560 try: 

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

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

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

564 

565 _s12 = Fsum(s12).fsum2_ 

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

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

568 r = gX.Inverse(lat, lon, p.lat2, p.lon2, outmask=gm) 

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

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

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

572 s *= (s2 * szi / c2) - (s * r.M21 / r.m12) # XXX _over? 

573 s12, t = _s12(c / s) # XXX _over? 

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

575 break 

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

577 iteration=i) 

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

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

580 return r 

581 

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

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

584 ''' 

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

586 

587 @Property_RO 

588 def rhumb(self): 

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

590 ''' 

591 return self._rhumb 

592 

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

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

595 

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

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

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

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

600 

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

602 ''' 

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

604 azi12=self.azi12, exact=self.exact, 

605 TMorder=self.TMorder, xTM=self.xTM) 

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

607 

608 @property_RO 

609 def TMorder(self): 

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

611 ''' 

612 return self.rhumb.TMorder 

613 

614 @Property_RO 

615 def xTM(self): 

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

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

618 ''' 

619 E = self.ellipsoid 

620 # ExactTransverseMercator doesn't handle spherical earth models 

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

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

623 

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

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

626 as current intersection estimate and central meridian. 

627 ''' 

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

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

630 

631 

632__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase) 

633 

634 

635if __name__ == '__main__': 

636 

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

638 

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

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

641 

642 for i in range(1, 10): 

643 s = .5e6 + 1e6 / i 

644 a = A.Position(s).lon2 

645 r = R.Position(s).lon2 

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

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

648 

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

650 for est in (None, 1e6): 

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

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

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

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

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

656 

657# % python3 -m pygeodesy.rhumbBase 

658 

659# Position.lon2 11.61455846901637 vs 11.61455846901637, diff 3.05885e-16 

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

661# Position.lon2 6.28526067416369 vs 6.28526067416369, diff 1.41311e-16 

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

663# Position.lon2 5.25385527435707 vs 5.25385527435707, diff 3.38105e-16 

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

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

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

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

668 

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

670# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9 

671 

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

673# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9 

674 

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

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

677 

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

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

680 

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

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

683 

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

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

686 

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

688# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9 

689 

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

691# NearestOn4Tuple(lat=45.0, lon=15.830286, distance=1977981.142985, normal=135.0), iteration=9 

692 

693# **) MIT License 

694# 

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

696# 

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

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

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

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

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

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

703# 

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

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

706# 

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

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

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

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

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

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

713# OTHER DEALINGS IN THE SOFTWARE.