Coverage for pygeodesy/rhumbaux.py: 95%

147 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-10-11 16:04 -0400

1# -*- coding: utf-8 -*- 

2 

3u'''A pure Python version of I{Karney}'s C++ classes U{Rhumb 

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

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1RhumbLine.html>} from 

6I{GeographicLib version 2.2+} renamed as L{RhumbAux} respectively L{RhumbLineAux}. 

7 

8Class L{RhumbLineAux} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively 

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

10geodesic or perpendicular rhumb line from an other point. 

11 

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

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

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

15utility U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online rhumb 

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

17 

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

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

20 

21@note: Class L{AuxDST} requires package U{numpy<https://PyPI.org/project/numpy>} to be installed, 

22 version 1.16 or newer, needed only for I{exact} area calculations C{S12} in L{RhumbAux} 

23 and L{RhumbLineAux}. 

24''' 

25# make sure int/int division yields float quotient 

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

27 

28from pygeodesy.auxilats.auxAngle import AuxMu, AuxPhi, hypot 

29from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw 

30# from pygeodesy.auxilats.auxDST import AuxDST # _MODS 

31from pygeodesy.auxilats.auxily import _Dlam, _Dp0Dpsi, _Ufloats 

32from pygeodesy.basics import copysign0, _reverange, _zip, _xkwds_get 

33from pygeodesy.constants import EPS_2, MANT_DIG, PI4, isinf, \ 

34 _0_0, _4_0, _720_0, _log2, _over 

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

36# from pygeodesy.errors import itemsorted, _xkwds_get # from .basics, ... 

37from pygeodesy.karney import Caps, _polynomial, _EWGS84 

38# from pygeodesy.fmath import hypot # from .auxilats.auxAngle 

39from pygeodesy.interns import NN, _COMMASPACE_ 

40from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS 

41# from pygeodesy.props import Property_RO # from .rhumbBase 

42from pygeodesy.rhumbBase import RhumbBase, RhumbLineBase, \ 

43 itemsorted, pairs, Property_RO 

44# from pygeodesy.streprs import pairs # from .rhumbBase 

45 

46from math import ceil as _ceil, fabs, radians 

47 

48__all__ = _ALL_LAZY.rhumbaux 

49__version__ = '23.09.20' 

50 

51# DIGITS = (sizeof(real) * 8) bits 

52# = (ctypes.sizeof(ctypes.c_double(1.0)) * 8) bits 

53# For |n| <= 0.99, actual max for doubles is 2163. This scales 

54# as DIGITS and for long doubles (GEOGRAPHICLIB_PRECISION = 3, 

55# DIGITS = 64), this becomes 2163 * 64 / 53 = 2612. Round this 

56# up to 2^12 = 4096 and scale this by DIGITS//64 if DIGITS > 64. 

57# 

58# 64 = DIGITS for long double, 6 = 12 - _log2(64) 

59_Lbits = 1 << (int(_ceil(_log2(max(MANT_DIG, 64)))) + 6) 

60 

61 

62class RhumbAux(RhumbBase): 

63 '''Class to solve the I{direct} and I{inverse rhumb} problems, based 

64 on I{Auxiliary Latitudes} for accuracy near the poles. 

65 

66 @note: Package U{numpy<https://PyPI.org/project/numpy>} must be 

67 installed, version 1.16 or later. 

68 ''' 

69 

70 def __init__(self, a_earth=_EWGS84, f=None, exact=True, name=NN, **TMorder): # PYCHOK signature 

71 '''New C{rhumbaux.RhumbAux}. 

72 

73 @kwarg a_earth: This rhumb's earth model (L{Ellipsoid}, L{Ellipsoid2}, 

74 L{a_f2Tuple}, L{Datum}, 2-tuple C{(a, f)}) or the 

75 (equatorial) radius (C{scalar}). 

76 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is 

77 a C{scalar}, ignored otherwise. 

78 @kwarg exact: If C{True}, use the exact expressions for the I{Auxiliary 

79 Latitudes}, otherwise use the I{Fourier} series expansion 

80 (C{bool}), see also property C{exact}. 

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

82 @kwarg TMorder: Optional keyword argument B{C{TMorder}}, see property 

83 C{TMorder}. 

84 

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

86 required for area C{S12} when C{B{exact} is True}. 

87 

88 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{RA_TMorder}}. 

89 ''' 

90 RhumbBase.__init__(self, a_earth, f, exact, name) 

91 if TMorder: 

92 self.Tmorder = _xkwds_get(TMorder, TMorder=RhumbBase._mTM) 

93 

94 def areaux(self, **exact): 

95 '''Get this ellipsoid's B{C{exact}} surface area (C{meter} I{squared}). 

96 

97 @kwarg exact: Optional C{exact} (C{bool}), overriding this rhumb's 

98 C{exact} setting, if C{True}, use the exact expression 

99 for the authalic radius otherwise the I{Taylor} series. 

100 

101 @return: The (signed?) surface area (C{meter} I{squared}). 

102 

103 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)} exceeds 

104 property C{f_max}. 

105 

106 @note: The area of a polygon encircling a pole can be found by adding 

107 C{areaux / 2} to the sum of C{S12} for each side of the polygon. 

108 

109 @see: U{The area of rhumb polygons<https://ArXiv.org/pdf/2303.03219.pdf>} 

110 and method L{auxilats.AuxLat.AuthalicRadius2}. 

111 ''' 

112 x = _xkwds_get(exact, exact=self.exact) 

113 a = (self._c2 * _720_0) if bool(x) is self.exact else \ 

114 (self._auxD.AuthalicRadius2(exact=x, f_max=self.f_max) * PI4) 

115 return a 

116 

117 @Property_RO 

118 def _auxD(self): 

119 return AuxDLat(self.ellipsoid) 

120 

121 @Property_RO 

122 def _c2(self): # radians makes _c2 a factor per degree 

123 return radians(self._auxD.AuthalicRadius2(exact=self.exact, f_max=self.f_max)) 

124 

125 def _DMu_DPsi(self, Phi1, Phi2, Chi1, Chi2): 

126 xD = self._auxD 

127 r = xD.DRectifying(Phi1, Phi2) if self.exact else \ 

128 xD.CRectifying(Chi1, Chi2) 

129 if r: 

130 r = _over(r, xD.DIsometric(Phi1, Phi2) if self.exact else 

131 _Dlam(Chi1.tan, Chi2.tan)) # not Lambertian! 

132 return r 

133 

134 def _Inverse4(self, lon12, r, outmask): 

135 '''(INTERNAL) See method C{RhumbBase.Inverse}. 

136 ''' 

137 psi1, Chi1, Phi1 = self._psiChiPhi3(r.lat1) 

138 psi2, Chi2, Phi2 = self._psiChiPhi3(r.lat2) 

139 psi12 = psi2 - psi1 # radians 

140 lam12 = radians(lon12) 

141 if (outmask & Caps.DISTANCE): 

142 if isinf(psi1) or isinf(psi2): # PYCHOK no cover 

143 d = Phi2.toMu(self).toRadians 

144 d -= Phi1.toMu(self).toRadians 

145 s = fabs(d) 

146 else: # dmu/dpsi = dmu/dchi/dpsi/dchi 

147 s = hypot(lam12, psi12) 

148 if s: 

149 s *= self._DMu_DPsi(Phi1, Phi2, Chi1, Chi2) 

150 s *= self._rrm 

151 a = _over(s, self._mpd) 

152 r.set_(a12=copysign0(a, s), s12=s) 

153 return lam12, psi12, Chi1, Chi2 

154 

155 def _latPhi2(self, mu): 

156 Mu = AuxMu.fromDegrees(mu) 

157 Phi = Mu.toPhi(self) 

158 return Phi.toDegrees, Phi 

159 

160 @Property_RO 

161 def _mpd(self): # meter per degree 

162 return radians(self._rrm) # == self.ellipsoid._Lpd 

163 

164 def _psiChiPhi3(self, lat): 

165 Phi = AuxPhi.fromDegrees(lat) 

166 Chi = Phi.toChi(self) 

167 psi = Chi.toLambertianRadians 

168 return psi, Chi, Phi 

169 

170 @Property_RO 

171 def _RA(self): # get the coefficients for area calculation 

172 return tuple(_RAintegrate(self._auxD) if self.exact else 

173 _RAseries(self._auxD)) 

174 

175# _RhumbLine = RhumbLineAux # see further below 

176 

177 @Property_RO 

178 def _rrm(self): 

179 return self._auxD.RectifyingRadius(exact=self.exact) 

180 

181 def _S12d(self, Chix, Chiy, lon12): # degrees 

182 '''(INTERNAL) Compute the area C{S12} from C{._meanSinXi(Chix, Chiy) * .c2 * lon12}. 

183 ''' 

184 pP, xD = self._RA, self._auxD 

185 

186 tx, Phix = Chix.tan, Chix.toPhi(self) 

187 ty, Phiy = Chiy.tan, Chiy.toPhi(self) 

188 

189 dD = xD.DParametric(Phix, Phiy) if self.exact else \ 

190 xD.CParametric(Chix, Chiy) 

191 if dD: 

192 dD = _over(dD, xD.DIsometric(Phix, Phiy) if self.exact else 

193 _Dlam(tx, ty)) # not Lambertian! 

194 dD *= _DClenshaw(False, Phix.toBeta(self).normalized, 

195 Phiy.toBeta(self).normalized, 

196 pP, min(len(pP), xD.ALorder)) # Fsum 

197 dD += _Dp0Dpsi(tx, ty) 

198 dD *= self._c2 * lon12 

199 return float(dD) 

200 

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

202 '''Return this C{Rhumb} as string. 

203 

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

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

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

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

208 

209 @return: Tuple items (C{str}). 

210 ''' 

211 d = dict(ellipsoid=self.ellipsoid, exact=self.exact, 

212 TMorder=self.TMorder) 

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

214 

215 

216class RhumbLineAux(RhumbLineBase): 

217 '''Compute one or several points on a single rhumb line. 

218 

219 Class C{RhumbLineAux} facilitates the determination of points 

220 on a single rhumb line. The starting point (C{lat1}, C{lon1}) 

221 and the azimuth C{azi12} are specified once. 

222 

223 Method C{RhumbLineAux.Position} returns the location of an 

224 other point and optionally the distance C{s12} along and the 

225 area C{S12} under the rhumb line. 

226 

227 Method C{RhumbLineAux.intersection2} finds the intersection 

228 between two rhumb lines. 

229 

230 Method C{RhumbLineAux.nearestOn4} computes the nearest point 

231 on and the distance to a rhumb line in different ways. 

232 ''' 

233 _Rhumb = RhumbAux # rhumbaux.RhumbAux 

234 

235 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, **caps_name): # PYCHOK signature 

236 '''New C{rhumbaux.RhumbLineAux}. 

237 

238 @arg rhumb: The rhumb reference (C{rhumbaux.RhumbAux}). 

239 @kwarg lat1: Latitude of the start point (C{degrees90}). 

240 @kwarg lon1: Longitude of the start point (C{degrees180}). 

241 @kwarg azi12: Azimuth of this rhumb line (compass C{degrees}). 

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

243 C{B{caps}=0}, a bit-or'ed combination of L{Caps} 

244 values specifying the required capabilities. Include 

245 C{Caps.LINE_OFF} if updates to the B{C{rhumb}} should 

246 I{not} be reflected in this rhumb line. 

247 ''' 

248 RhumbLineBase.__init__(self, rhumb, lat1, lon1, azi12, **caps_name) 

249 

250 @Property_RO 

251 def _Chi1(self): 

252 return self._Phi1.toChi(self.rhumb) 

253 

254 @Property_RO 

255 def _mu1(self): 

256 '''(INTERNAL) Get the I{rectifying auxiliary} latitude (C{degrees}). 

257 ''' 

258 return self._Phi1.toMu(self.rhumb).toDegrees 

259 

260 def _mu2lat(self, mu): 

261 '''(INTERNAL) Get the inverse I{rectifying auxiliary} latitude (C{degrees}). 

262 ''' 

263 lat, _ = self.rhumb._latPhi2(mu) 

264 return lat 

265 

266 @Property_RO 

267 def _Phi1(self): 

268 return AuxPhi.fromDegrees(self.lat1) 

269 

270 def _Position4(self, a12, mu2, *unused): # PYCHOK s12, mu2 

271 '''(INTERNAL) See method C{RhumbLineBase._Position}. 

272 ''' 

273 R = self.rhumb 

274 lat2, Phi2 = R._latPhi2(mu2) 

275 Chi2 = Phi2.toChi(R) 

276 Chi1 = self._Chi1 

277 lon2 = self._salp * a12 

278 if lon2: 

279 m = R._DMu_DPsi(self._Phi1, Phi2, Chi1, Chi2) 

280 lon2 = _over(lon2, m) 

281 return lat2, lon2, Chi1, Chi2 

282 

283# @Property_RO 

284# def _psi1(self): 

285# return self._Chi1.toLambertianRadians 

286 

287RhumbAux._RhumbLine = RhumbLineAux # PYCHOK see RhumbBase._RhumbLine 

288 

289 

290def _RAintegrate(auxD): 

291 # Compute coefficients by Fourier transform of integrand 

292 L = 2 

293 fft = _MODS.auxilats.auxDST.AuxDST(L) 

294 f = auxD._qIntegrand 

295 c_ = fft.transform(f) 

296 pP = [] 

297 _P = pP.append 

298 # assert L < _Lbits 

299 while L < _Lbits: 

300 L = fft.reset(L) * 2 

301 c_ = fft.refine(f, c_, _0_0) # sentine[L] 

302 # assert len(c_) == L + 1 

303 pP[:], k = [], -1 

304 for j in range(1, L + 1): 

305 # Compute Fourier coefficients of integral 

306 p = (c_[j - 1] + c_[j]) / (_4_0 * j) 

307 if fabs(p) > EPS_2: 

308 k = -1 # run interrupted 

309 else: 

310 if k < 0: 

311 k = 1 # mark as first small value 

312 if (j - k) >= ((j + 7) // 8): 

313 # run of at least (j - 1) // 8 small values 

314 return pP[:j] # break while L loop 

315 _P(-p) 

316 return pP # no convergence, use pP as-is 

317 

318 

319def _RAseries(auxD): 

320 # Series expansions in n for Fourier coeffients of the integral 

321 # @see: U{"Series expansions for computing rhumb areas" 

322 # <https:#DOI.org/10.5281/zenodo.7685484>}. 

323 d = n = auxD._n 

324 i = 0 

325 aL = auxD.ALorder 

326 Cs = _RACoeffs[aL] 

327 # assert len(Cs) == (aL * (aL + 1)) // 2 

328 pP = [] 

329 _P = pP.append 

330 _p = _polynomial 

331 for m in _reverange(aL): # order 

332 j = i + m + 1 

333 _P(_p(n, Cs, i, j) * d) 

334 d *= n 

335 i = j 

336 # assert i == len(pP) 

337 return pP 

338 

339 

340_f, _u = float, _Ufloats() 

341_RACoeffs = { # Rhumb Area Coefficients in matrix Q 

342 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

343 596 / _f(2025), -398 / _f(945), 22 / _f(45), -1 / _f(3), 

344 1543 / _f(4725), -118 / _f(315), 1 / _f(5), 

345 152 / _f(945), -17 / _f(315), 

346 5 / _f(252)), 

347 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

348 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45), 

349 -1 / _f(3), 

350 -24562 / _f(155925), 1543 / _f(4725), -118 / _f(315), 1 / _f(5), 

351 -38068 / _f(155925), 152 / _f(945), -17 / _f(315), 

352 -752 / _f(10395), 5 / _f(252), 

353 -101 / _f(17325)), 

354 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

355 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025), 

356 -398 / _f(945), 22 / _f(45), -1 / _f(3), 

357 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725), 

358 -118 / _f(315), 1 / _f(5), 

359 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945), 

360 -17 / _f(315), 

361 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252), 

362 62464 / _f(2027025), -101 / _f(17325), 

363 11537 / _f(4054050)), 

364 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

365 -565017322 / _f(1915538625), 138734126 / _f(638512875), 

366 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45), 

367 -1 / _f(3), 

368 -1969276 / _f(58046625), 17749373 / _f(425675250), -24562 / _f(155925), 

369 1543 / _f(4725), -118 / _f(315), 1 / _f(5), 

370 -58573784 / _f(638512875), 1882432 / _f(8513505), -38068 / _f(155925), 

371 152 / _f(945), -17 / _f(315), 

372 -6975184 / _f(42567525), 268864 / _f(2027025), -752 / _f(10395), 

373 5 / _f(252), 

374 -112832 / _f(1447875), 62464 / _f(2027025), -101 / _f(17325), 

375 -4096 / _f(289575), 11537 / _f(4054050), 

376 -311 / _f(525525)), 

377 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

378 188270561816 / _f(488462349375), -565017322 / _f(1915538625), 

379 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025), 

380 -398 / _f(945), 22 / _f(45), -1 / _f(3), 

381 2332829602 / _f(23260111875), -1969276 / _f(58046625), 

382 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725), 

383 -118 / _f(315), 1 / _f(5), 

384 -41570288 / _f(930404475), -58573784 / _f(638512875), 

385 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945), 

386 -17 / _f(315), 

387 1538774036 / _f(10854718875), -6975184 / _f(42567525), 

388 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252), 

389 436821248 / _f(3618239625), -112832 / _f(1447875), 

390 62464 / _f(2027025), -101 / _f(17325), 

391 3059776 / _f(80405325), -4096 / _f(289575), 11537 / _f(4054050), 

392 4193792 / _f(723647925), -311 / _f(525525), 

393 1097653 / _f(1929727800)) 

394} 

395del _f, _u, _Ufloats 

396 

397__all__ += _ALL_DOCS(Caps, RhumbAux, RhumbLineAux) 

398 

399if __name__ == '__main__': 

400 

401 from pygeodesy import printf, RhumbAux # PYCHOK RhumbAux is __main__.RhumbAux 

402 

403 def _re(fmt, r3, x3): 

404 e3 = [] 

405 for r, x in _zip(r3, x3): # strict=True 

406 e = fabs(r - x) / fabs(x) 

407 e3.append('%.g' % (e,)) 

408 printf((fmt % r3) + ' rel errors: ' + ', '.join(e3)) 

409 

410 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve -p 9> version 2.2 

411 rhumb = RhumbAux(exact=True) # WGS84 default 

412 printf('# %r\n', rhumb) 

413 r = rhumb.Direct8(40.6, -73.8, 51, 5.5e6) # from JFK about NE 

414 _re('# JFK NE lat2=%.8f, lon2=%.8f, S12=%.1f', (r.lat2, r.lon2, r.S12), (71.688899882813, 0.2555198244234, 44095641862956.11)) 

415 r = rhumb.Inverse8(40.6, -73.8, 51.6, -0.5) # JFK to LHR 

416 _re('# JFK-LHR azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (77.7683897102557, 5771083.38332803, 37395209100030.39)) 

417 r = rhumb.Inverse8(40.6, -73.8, 35.8, 140.3) # JFK to Tokyo Narita 

418 _re('# JFK-NRT azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (-92.38888798169965, 12782581.067684170, -63760642939072.50)) 

419 

420# % python3 -m pygeodesy.rhumbaux 

421 

422# RhumbAux(TMorder=6, ellipsoid=Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594), exact=True) 

423 

424# JFK NE lat2=71.68889988, lon2=0.25551982, S12=44095641862956.1 rel errors: 4e-11, 2e-08, 5e-16 

425# JFK-LHR azi12=77.76838971, s12=5771083.383 S12=37395209100030.3 rel errors: 3e-12, 5e-15, 6e-16 

426# JFK-NRT azi12=-92.38888798, s12=12782581.068 S12=-63760642939072.5 rel errors: 2e-16, 3e-16, 6e-16