1

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

3

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

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

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

7 from 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

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

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

12 geodesic or perpendicular rhumb line.

13

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

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

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

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

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

19

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

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

22'''

23# make sure int/int division yields float quotient

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

25

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

27 from 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

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

31 _xdatum, _xkwds, _Xorder

32# from pygeodesy.etm import ExactTransverseMercator # _MODS

33 from pygeodesy.fmath import euclid, favg, fabs

34# from pygeodesy.fsums import Fsum # _MODS

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

36 _no_, _parallel_, _under

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

38 _norm180, _xinstanceof

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

40 from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS

41# from pygeodesy.named import notOverloaded # _MODS

42 from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple

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

44 from pygeodesy.streprs import Fmt, pairs, unstr

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

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

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

48

49# from math import fabs # from .fmath

50

51__all__ = ()

52__version__ = '23.08.09'

53

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

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

56

57

58 class _Lat(Lat):

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

60 '''

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

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

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

64

65

66 class _Lon(Lon):

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

68 '''

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

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

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

72

73

74 def _update_all_rls(r):

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

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

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

78 '''

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

80 _update_all(r)

81 for rl in _rls: # PYCHOK use weakref?

82 if rl._rhumb is r:

83 _update_all(rl)

84

85

86 class RhumbBase(_CapsBase):

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

88 '''

89 _E = _EWGS84

90 _exact = True

91 _f_max = _0_01

92 _mTM = 6 # see .TMorder

93

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

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

96 '''

97 if f is not None:

98 self.ellipsoid = a_earth, f

99 elif a_earth not in (_EWGS84, None):

100 self.ellipsoid = a_earth

101 if not exact:

102 self.exact = False

103 if name:

104 self.name = name

105

106 @Property_RO

107 def a(self):

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

109 '''

110 return self.ellipsoid.a

111

112 equatoradius = a

113

114 @Property_RO

115 def b(self):

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

117 '''

118 return self.ellipsoid.b

119

120 polaradius = b

121

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

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

124 '''

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

126

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

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

129 '''

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

131

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

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

134 '''

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

136

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

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

139 '''

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

141

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

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

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

145

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

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

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

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

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

151 L{Caps} values specifying the required capabilities.

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

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

154

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

156 C{.Position} to compute each point.

157

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

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

160 '''

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

162 **caps_name)

163

164 Line = DirectLine # synonyms

165

166 @Property

167 def ellipsoid(self):

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

169 '''

170 return self._E

171

172 @ellipsoid.setter # PYCHOK setter!

173 def ellipsoid(self, a_earth_f):

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

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

176

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

178 '''

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

180 if self._E != E:

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

182 _update_all_rls(self)

183 self._E = E

184

185 @Property

186 def exact(self):

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

188 '''

189 return self._exact

190

191 @exact.setter # PYCHOK setter!

192 def exact(self, exact):

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

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

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

196

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

198 exceeds non-zero C{f_max}.

199

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

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

202 '''

203 x = bool(exact)

204 if self._exact != x:

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

206 _update_all_rls(self)

207 self._exact = x

208

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

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

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

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

213

214 @Property_RO

215 def f(self):

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

217 '''

218 return self.ellipsoid.f

219

220 flattening = f

221

222 @property

223 def f_max(self):

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

225 '''

226 return self._f_max

227

228 @f_max.setter # PYCHOK setter!

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

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

231

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

233 exceeds non-zero C{f_max}.

234 '''

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

236 if self._f_max != f:

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

238 self._f_max = f

239

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

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

242 '''

243 if wrap:

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

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

246

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

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

249 '''

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

251

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

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

254 '''

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

256

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

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

259 '''

260 if wrap:

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

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

263

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

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

266

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

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

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

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

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

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

273 L{Caps} values specifying the required capabilities.

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

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

276

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

278 C{.Position} to compute each point.

279

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

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

282 '''

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

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

285 **caps_name)

286

287 @Property_RO

288 def _RhumbLine(self): # PYCHOK no cover

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

290 '''

291 _MODS.named.notOverloaded(self)

292

293 def _TMorder(self, order):

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

295

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

297 '''

298 x = self.exact

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

300 if self._mTM != n:

301 _update_all_rls(self)

302 self._mTM = n

303 x = False

304 return x

305

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

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

308 '''

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

310

311

312 class RhumbLineBase(RhumbBase):

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

314 '''

315 _azi12 = _0_0

316# _caps = 0

317# _debug = 0

318# _lat1 = _0_0

319# _lon1 = _0_0

320# _lon12 = _0_0

321 _salp = _0_0

322 _calp = _1_0

323 _Rhumb = RhumbBase # compatible C{Rhumb} class

324 _rhumb = None # C{Rhumb} instance

325

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

327 '''New C{RhumbLine}.

328 '''

329 _xinstanceof(self._Rhumb, rhumb=rhumb)

330

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

332 self._lon1 = _Lon(lon1= lon1)

333 self._lon12 = _norm180(self._lon1)

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

335 self.azi12 = _norm180(azi12)

336

337 n = name or rhumb.name

338 if n:

339 self.name=n

340

341 self._caps = caps

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

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

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

345 else:

346 self._rhumb = rhumb

347 _rls.append(self)

348

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

350 if _rls: # may be empty or None

351 try: # PYCHOK no cover

352 _rls.remove(self)

353 except (TypeError, ValueError):

354 pass

355 self._rhumb = None

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

357

358 @Property

359 def azi12(self):

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

361 '''

362 return self._azi12

363

364 @azi12.setter # PYCHOK setter!

365 def azi12(self, azi12):

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

367 '''

368 z = _norm180(azi12)

369 if self._azi12 != z:

370 if self._rhumb:

371 _update_all(self)

372 self._azi12 = z

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

374

375 @Property_RO

376 def azi12_sincos2(self): # PYCHOK no cover

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

378 '''

379 return self._scalp, self._calp

380

381 def distance2(self, lat, lon):

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

383 point to this rhumb line's start point.

384

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

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

387

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

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

390

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

392 L{RhumbLineAux}.

393 '''

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

395# outmask=Caps.AZIMUTH_DISTANCE)

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

397

398 @property_RO

399 def ellipsoid(self):

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

401 '''

402 return self.rhumb.ellipsoid

403

404 @property_RO

405 def exact(self):

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

407 '''

408 return self.rhumb.exact

409

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

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

412

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

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

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

416

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

418 C{lon}gitude of the intersection point.

419

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

421 no intersection for an other reason.

422

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

424 L{pygeodesy.intersection3d3}.

425

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

427 L{ExactTransverseMercator} or L{KTransverseMercator}

428 projection and invoking function L{pygeodesy.intersection3d3}

429 in that domain.

430 '''

431 _xinstanceof(RhumbLineBase, other=other)

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

433 try:

434 if other is self:

435 raise ValueError(_coincident_)

436 # make globals and invariants locals

437 _diff = euclid # approximate length

438 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3

439 _LL2T = LatLon2Tuple

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

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

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

443 t = fabs(s_az - o_az)

444 if t < EPS or fabs(t - _180_0) < EPS:

445 raise ValueError(_parallel_)

446 # use halfway point as initial estimate

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

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

449 for i in range(1, _TRIPS):

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

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

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

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

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

455 if d < tol:

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

457 name=self.intersection2.__name__)

458 except Exception as x:

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

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

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

462

463 @Property_RO

464 def lat1(self):

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

466 '''

467 return self._lat1

468

469 @Property_RO

470 def lon1(self):

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

472 '''

473 return self._lon1

474

475 @Property_RO

476 def latlon1(self):

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

478 '''

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

480

481 def _mu22(self, mu12, mu1):

482 mu2 = mu12 + mu1

483 x90 = fabs(mu2) > 90

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

485 mu2 = _norm180(mu2)

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

487 mu2 = _norm180(_180_0 - mu2)

488 return mu2, x90

489

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

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

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

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

494

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

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

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

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

499 respectively C{B{exact} is not None}.

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

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

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

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

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

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

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

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

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

509

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

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

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

513 C{normal} azimuth at the intersection.

514

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

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

517 package not found or not installed.

518

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

520 no intersection for an other reason.

521

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

523 L{pygeodesy.intersection3d3}.

524 '''

525 if exact is None:

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

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

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

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

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

531 iteration=p.iteration)

532 else: # C{rhumb-intercept}

533 azi = self.azi12

534 szi = self._salp

535 E = self.ellipsoid

536 gX = E.geodesic_(exact=exact)

537 Cs = Caps

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

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

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

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

542 _, c = sincos2d(d)

543 s12 = c * r.s12 # signed

544 else:

545 s12 = Meter(est=est)

546 try:

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

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

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

550

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

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

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

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

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

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

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

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

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

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

561 break

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

563 iteration=i)

564 except Exception as x: # Fsum Value-, ZeroDivisionError

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

566 return r

567

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

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

570 '''

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

572

573 @Property_RO

574 def rhumb(self):

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

576 '''

577 return self._rhumb

578

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

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

581

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

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

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

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

586

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

588 '''

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

590 azi12=self.azi12, exact=self.exact,

591 TMorder=self.TMorder, xTM=self.xTM)

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

593

594 @property_RO

595 def TMorder(self):

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

597 '''

598 return self.rhumb.TMorder

599

600 @Property_RO

601 def xTM(self):

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

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

604 '''

605 E = self.ellipsoid

606 # ExactTransverseMercator doesn't handle spherical earth models

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

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

609

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

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

612 as current intersection estimate and central meridian.

613 '''

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

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

616

617

618__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase)

619

620

621 if __name__ == '__main__':

622

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

624

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

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

627

628 n = 1

629 for i in range(1, 10):

630 s = .5e6 + 1e6 / i

631 a = A.Position(s).lon2

632 r = R.Position(s).lon2

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

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

635 n = 0

636

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

638 for est in (None, 1e6):

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

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

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

642 a.toRepr(), a.iteration, exact, est, r.toRepr(), r.iteration, nl=1)

643

644# % python3 -m pygeodesy.rhumbBase

645

646# Positions 11.614558469016366 vs 11.61455846901637, diff 0.00%

647# Positions 7.589823028268422 vs 7.589823028268423, diff 0.00%

648# Positions 6.285260674163688 vs 6.285260674163687, diff 0.00%

649# Positions 5.639389953251457 vs 5.6393899532514595, diff 0.00%

650# Positions 5.253855274357075 vs 5.253855274357073, diff 0.00%

651# Positions 4.9976460429038 vs 4.9976460429038045, diff 0.00%

652# Positions 4.815033637404729 vs 4.81503363740473, diff 0.00%

653# Positions 4.678288217488354 vs 4.678288217488353, diff 0.00%

654# Positions 4.572056679062825 vs 4.572056679062826, diff 0.00%

655

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

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

658

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

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

661

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

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

664

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

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

667

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

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

670

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

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

673

674# **) MIT License

675#

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

677#

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

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

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

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

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

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

684#

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

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

687#

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

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

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

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

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

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

694# OTHER DEALINGS IN THE SOFTWARE.