Coverage for pygeodesy/rhumb/aux_.py: 95%

145 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-15 16:36 -0400

1 

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

3 

4u'''A pure Python version of I{Karney}'s I{Auxiliary Latitudes}, C++ classes U{Rhumb 

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

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

7I{GeographicLib version 2.2+} renamed to L{RhumbAux} respectively L{RhumbLineAux}. 

8 

9Class L{RhumbLineAux} has been enhanced with methods C{Intersecant2}, C{Intersection} and C{PlumbTo} 

10to iteratively find the intersection of a rhumb line and a circle or an other rhumb line, respectively 

11a perpendicular geodesic or other rhumb line. 

12 

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

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

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

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

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

18 

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

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

21 

22@note: C{S12} area calculations in classes L{RhumbAux} and L{RhumbLineAux} depend on class L{AuxDST} which 

23 requires U{numpy<https://PyPI.org/project/numpy>} to be installed, version 1.16 or newer. 

24 

25@note: Windows reserves file names U{AUX, COM[#], CON, LPT[#], NUL, PRN<https://learn.Microsoft.com/en-us/ 

26 windows/win32/fileio/naming-a-file#naming-conventions>} with and without extension. 

27''' 

28# make sure int/int division yields float quotient 

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

30 

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

32from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw 

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

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

35from pygeodesy.basics import copysign0, _reverange, _xkwds_get 

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

37 _0_0, _4_0, _720_0, _log2, _over 

38from pygeodesy.datums import _WGS84, NN 

39# from pygeodesy.errors import _xkwds_get # from .basics 

40from pygeodesy.karney import Caps, _polynomial 

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

42# from pygeodesy.interns import NN # from .datums 

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

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

45from pygeodesy.rhumb.bases import RhumbBase, RhumbLineBase, Property_RO 

46 

47from math import ceil as _ceil, fabs, radians 

48 

49__all__ = _ALL_LAZY.rhumb_aux_ 

50__version__ = '23.12.09' 

51 

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

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

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

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

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

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

58# 

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

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

61 

62 

63class RhumbAux(RhumbBase): 

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

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

66 

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

68 installed, version 1.16 or later. 

69 ''' 

70 

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

72 '''New C{RhumbAux}. 

73 

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

75 L{Ellipsoid2}, L{a_f2Tuple}, 2-tuple C{(a, f)}) or 

76 the (equatorial) radius (C{meter}, conventionally). 

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

78 C{scalar}, ignored otherwise. 

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

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

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

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

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

84 C{TMorder}. 

85 

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

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

88 

89 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{TMorder}}. 

90 ''' 

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

92 if TMorder: 

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

94 

95 def areaux(self, **exact): 

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

97 

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

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

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

101 

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

103 

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

105 property C{f_max}. 

106 

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

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

109 

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

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

112 ''' 

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

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

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

116 return a 

117 

118 @Property_RO 

119 def _auxD(self): 

120 return AuxDLat(self.ellipsoid) 

121 

122 @Property_RO 

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

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

125 

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

127 xD = self._auxD 

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

129 xD.CRectifying(Chi1, Chi2) 

130 if r: 

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

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

133 return r 

134 

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

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

137 ''' 

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

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

140 psi12 = psi2 - psi1 # radians 

141 lam12 = radians(lon12) 

142 if (outmask & Caps.DISTANCE): 

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

144 s = fabs(Phi2.toMu(self).toRadians - 

145 Phi1.toMu(self).toRadians) 

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 _mpr = _rrm # meter per radian, see _mpd 

182 

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

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

185 ''' 

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

187 

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

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

190 

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

192 xD.CParametric(Chix, Chiy) 

193 if dD: 

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

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

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

197 Phiy.toBeta(self).normalized, 

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

199 dD += _Dp0Dpsi(tx, ty) 

200 dD *= self._c2 * lon12 

201 return float(dD) 

202 

203 

204class RhumbLineAux(RhumbLineBase): 

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

206 

207 Class C{RhumbLineAux} facilitates the determination of points 

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

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

210 ''' 

211 _Rhumb = RhumbAux # rhumb.aux_.RhumbAux 

212 

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

214 '''New C{RhumbLineAux}. 

215 

216 @arg rhumb: The rhumb reference (L{RhumbAux}). 

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

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

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

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

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

222 values specifying the required capabilities. Include 

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

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

225 ''' 

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

227 

228 @Property_RO 

229 def _Chi1(self): 

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

231 

232 @Property_RO 

233 def _mu1(self): 

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

235 ''' 

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

237 

238 def _mu2lat(self, mu): 

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

240 ''' 

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

242 return lat 

243 

244 @Property_RO 

245 def _Phi1(self): 

246 return AuxPhi.fromDegrees(self.lat1) 

247 

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

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

250 ''' 

251 R = self.rhumb 

252 lat2, Phi2 = R._latPhi2(mu2) 

253 Chi2 = Phi2.toChi(R) 

254 Chi1 = self._Chi1 

255 lon2 = self._salp * a12 

256 if lon2: 

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

258 lon2 = _over(lon2, m) 

259 return lat2, lon2, Chi1, Chi2 

260 

261# @Property_RO 

262# def _psi1(self): 

263# return self._Chi1.toLambertianRadians 

264 

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

266 

267 

268def _RAintegrate(auxD): 

269 # Compute coefficients by Fourier transform of integrand 

270 L = 2 

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

272 f = auxD._qIntegrand 

273 c_ = fft.transform(f) 

274 pP = [] 

275 _P = pP.append 

276 # assert L < _Lbits 

277 while L < _Lbits: 

278 L = fft.reset(L) * 2 

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

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

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

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

283 # Compute Fourier coefficients of integral 

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

285 if fabs(p) > EPS_2: 

286 k = -1 # run interrupted 

287 else: 

288 if k < 0: 

289 k = 1 # mark as first small value 

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

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

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

293 _P(-p) 

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

295 

296 

297def _RAseries(auxD): 

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

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

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

301 d = n = auxD._n 

302 i = 0 

303 aL = auxD.ALorder 

304 Cs = _RACoeffs[aL] 

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

306 pP = [] 

307 _P = pP.append 

308 _p = _polynomial 

309 for m in _reverange(aL): # order 

310 j = i + m + 1 

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

312 d *= n 

313 i = j 

314 # assert i == len(pP) 

315 return pP 

316 

317 

318_f, _u = float, _Ufloats() 

319_RACoeffs = { # Rhumb Area Coefficients in matrix Q 

320 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

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

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

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

324 5 / _f(252)), 

325 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

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

327 -1 / _f(3), 

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

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

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

331 -101 / _f(17325)), 

332 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

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

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

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

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

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

338 -17 / _f(315), 

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

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

341 11537 / _f(4054050)), 

342 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

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

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

345 -1 / _f(3), 

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

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

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

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

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

351 5 / _f(252), 

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

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

354 -311 / _f(525525)), 

355 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

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

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

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

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

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

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

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

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

364 -17 / _f(315), 

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

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

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

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

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

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

371 1097653 / _f(1929727800)) 

372} 

373del _f, _u, _Ufloats 

374 

375__all__ += _ALL_DOCS(Caps) 

376 

377# **) MIT License 

378# 

379# Copyright (C) 2023-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

380# 

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

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

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

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

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

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

387# 

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

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

390# 

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

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

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

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

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

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

397# OTHER DEALINGS IN THE SOFTWARE.