Coverage for pygeodesy/rhumb/aux_.py: 96%
142 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
2# -*- coding: utf-8 -*-
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}.
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.
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>}.
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.
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.
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
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_get1
36from pygeodesy.constants import EPS_2, MANT_DIG, PI4, isinf, \
37 _0_0, _4_0, _720_0, _log2, _over
38# from pygeodesy.datums import _WGS84 # from .rhumb.bases
39# from pygeodesy.errors import _xkwds_get1 # from .basics
40from pygeodesy.karney import Caps, _polynomial
41# from pygeodesy.fmath import hypot # from .auxilats.auxAngle
42from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
43# from pygeodesy.props import Property_RO # from .rhumb.bases
44from pygeodesy.rhumb.bases import RhumbBase, RhumbLineBase, \
45 Property_RO, _WGS84
47from math import ceil as _ceil, fabs, radians
49__all__ = _ALL_LAZY.rhumb_aux_
50__version__ = '24.05.29'
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)
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.
67 @note: Package U{numpy<https://PyPI.org/project/numpy>} must be
68 installed, version 1.16 or later.
69 '''
71 def __init__(self, a_earth=_WGS84, f=None, exact=True, **TMorder_name): # PYCHOK signature
72 '''New C{RhumbAux}.
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 TMorder_name: Optional C{B{name}=NN} (C{str}) and optional
83 keyword argument C{B{TMorder}=6} for the order of
84 the L{KTransverseMercator}, see property C{TMorder}.
86 @raise ImportError: Package C{numpy} not found or not installed, only
87 required for area C{S12} when C{B{exact} is True}.
89 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{TMorder}}.
90 '''
91 RhumbBase.__init__(self, a_earth, f, exact, TMorder_name)
93 def areaux(self, **exact):
94 '''Get this ellipsoid's B{C{exact}} surface area (C{meter} I{squared}).
96 @kwarg exact: Optional C{exact} (C{bool}), overriding this rhumb's
97 C{exact} setting, if C{True}, use the exact expression
98 for the authalic radius otherwise the I{Taylor} series.
100 @return: The (signed?) surface area (C{meter} I{squared}).
102 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)} exceeds
103 property C{f_max}.
105 @note: The area of a polygon encircling a pole can be found by adding
106 C{areaux / 2} to the sum of C{S12} for each side of the polygon.
108 @see: U{The area of rhumb polygons<https://ArXiv.org/pdf/2303.03219.pdf>}
109 and method L{auxilats.AuxLat.AuthalicRadius2}.
110 '''
111 x = _xkwds_get1(exact, exact=self.exact)
112 a = (self._c2 * _720_0) if bool(x) is self.exact else (
113 self._auxD.AuthalicRadius2(exact=x, f_max=self.f_max) * PI4)
114 return a
116 @Property_RO
117 def _auxD(self):
118 return AuxDLat(self.ellipsoid)
120 @Property_RO
121 def _c2(self): # radians makes _c2 a factor per degree
122 return radians(self._auxD.AuthalicRadius2(exact=self.exact, f_max=self.f_max))
124 def _DMu_DPsi(self, Phi1, Phi2, Chi1, Chi2):
125 xD = self._auxD
126 r = xD.DRectifying(Phi1, Phi2) if self.exact else \
127 xD.CRectifying(Chi1, Chi2)
128 if r:
129 r = _over(r, xD.DIsometric(Phi1, Phi2) if self.exact else
130 _Dlam(Chi1.tan, Chi2.tan)) # not Lambertian!
131 return r
133 def _Inverse4(self, lon12, r, outmask):
134 '''(INTERNAL) See method C{RhumbBase.Inverse}.
135 '''
136 psi1, Chi1, Phi1 = self._psiChiPhi3(r.lat1)
137 psi2, Chi2, Phi2 = self._psiChiPhi3(r.lat2)
138 psi12 = psi2 - psi1 # radians
139 lam12 = radians(lon12)
140 if (outmask & Caps.DISTANCE):
141 if isinf(psi1) or isinf(psi2): # PYCHOK no cover
142 s = fabs(Phi2.toMu(self).toRadians -
143 Phi1.toMu(self).toRadians)
144 else: # dmu/dpsi = dmu/dchi/dpsi/dchi
145 s = hypot(lam12, psi12)
146 if s:
147 s *= self._DMu_DPsi(Phi1, Phi2, Chi1, Chi2)
148 s *= self._rrm
149 a = _over(s, self._mpd)
150 r.set_(a12=copysign0(a, s), s12=s)
151 return lam12, psi12, Chi1, Chi2
153 def _latPhi2(self, mu):
154 Mu = AuxMu.fromDegrees(mu)
155 Phi = Mu.toPhi(self)
156 return Phi.toDegrees, Phi
158 @Property_RO
159 def _mpd(self): # meter per degree
160 return radians(self._rrm) # == self.ellipsoid._Lpd
162 def _psiChiPhi3(self, lat):
163 Phi = AuxPhi.fromDegrees(lat)
164 Chi = Phi.toChi(self)
165 psi = Chi.toLambertianRadians
166 return psi, Chi, Phi
168 @Property_RO
169 def _RA(self): # get the coefficients for area calculation
170 return tuple(_RAintegrate(self._auxD) if self.exact else
171 _RAseries(self._auxD))
173# _RhumbLine = RhumbLineAux # see further below
175 @Property_RO
176 def _rrm(self):
177 return self._auxD.RectifyingRadius(exact=self.exact)
179 _mpr = _rrm # meter per radian, see _mpd
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
186 tx, Phix = Chix.tan, Chix.toPhi(self)
187 ty, Phiy = Chiy.tan, Chiy.toPhi(self)
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)
202class RhumbLineAux(RhumbLineBase):
203 '''Compute one or several points on a single rhumb line.
205 Class C{RhumbLineAux} facilitates the determination of points
206 on a single rhumb line. The starting point (C{lat1}, C{lon1})
207 and the azimuth C{azi12} are specified once.
208 '''
209 _Rhumb = RhumbAux # rhumb.aux_.RhumbAux
211 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, **caps_name): # PYCHOK signature
212 '''New C{RhumbLineAux}.
214 @arg rhumb: The rhumb reference (L{RhumbAux}).
215 @kwarg lat1: Latitude of the start point (C{degrees90}).
216 @kwarg lon1: Longitude of the start point (C{degrees180}).
217 @kwarg azi12: Azimuth of this rhumb line (compass C{degrees}).
218 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
219 C{B{caps}=0}, a bit-or'ed combination of L{Caps}
220 values specifying the required capabilities. Include
221 C{Caps.LINE_OFF} if updates to the B{C{rhumb}} should
222 I{not} be reflected in this rhumb line.
223 '''
224 RhumbLineBase.__init__(self, rhumb, lat1, lon1, azi12, **caps_name)
226 @Property_RO
227 def _Chi1(self):
228 return self._Phi1.toChi(self.rhumb)
230 @Property_RO
231 def _mu1(self):
232 '''(INTERNAL) Get the I{rectifying auxiliary} latitude (C{degrees}).
233 '''
234 return self._Phi1.toMu(self.rhumb).toDegrees
236 def _mu2lat(self, mu):
237 '''(INTERNAL) Get the inverse I{rectifying auxiliary} latitude (C{degrees}).
238 '''
239 lat, _ = self.rhumb._latPhi2(mu)
240 return lat
242 @Property_RO
243 def _Phi1(self):
244 return AuxPhi.fromDegrees(self.lat1)
246 def _Position4(self, a12, mu2, *unused): # PYCHOK s12, mu2
247 '''(INTERNAL) See method C{RhumbLineBase._Position}.
248 '''
249 R = self.rhumb
250 lat2, Phi2 = R._latPhi2(mu2)
251 Chi2 = Phi2.toChi(R)
252 Chi1 = self._Chi1
253 lon2 = self._salp * a12
254 if lon2:
255 m = R._DMu_DPsi(self._Phi1, Phi2, Chi1, Chi2)
256 lon2 = _over(lon2, m)
257 return lat2, lon2, Chi1, Chi2
259# @Property_RO
260# def _psi1(self):
261# return self._Chi1.toLambertianRadians
263RhumbAux._RhumbLine = RhumbLineAux # PYCHOK see RhumbBase._RhumbLine
266def _RAintegrate(auxD):
267 # Compute coefficients by Fourier transform of integrand
268 L = 2
269 fft = _MODS.auxilats.auxDST.AuxDST(L)
270 f = auxD._qIntegrand
271 c_ = fft.transform(f)
272 pP = []
273 _P = pP.append
274 # assert L < _Lbits
275 while L < _Lbits:
276 L = fft.reset(L) * 2
277 c_ = fft.refine(f, c_, _0_0) # sentine[L]
278 # assert len(c_) == L + 1
279 pP[:], k = [], -1
280 for j in range(1, L + 1):
281 # Compute Fourier coefficients of integral
282 p = (c_[j - 1] + c_[j]) / (_4_0 * j)
283 if fabs(p) > EPS_2:
284 k = -1 # run interrupted
285 else:
286 if k < 0:
287 k = 1 # mark as first small value
288 if (j - k) >= ((j + 7) // 8):
289 # run of at least (j - 1) // 8 small values
290 return pP[:j] # break while L loop
291 _P(-p)
292 return pP # no convergence, use pP as-is
295def _RAseries(auxD):
296 # Series expansions in n for Fourier coeffients of the integral
297 # @see: U{"Series expansions for computing rhumb areas"
298 # <https:#DOI.org/10.5281/zenodo.7685484>}.
299 d = n = auxD._n
300 i = 0
301 aL = auxD.ALorder
302 Cs = _RACoeffs[aL]
303 # assert len(Cs) == (aL * (aL + 1)) // 2
304 pP = []
305 _P = pP.append
306 _p = _polynomial
307 for m in _reverange(aL): # order
308 j = i + m + 1
309 _P(_p(n, Cs, i, j) * d)
310 d *= n
311 i = j
312 # assert i == len(pP)
313 return pP
316_f, _u = float, _Ufloats()
317_RACoeffs = { # Rhumb Area Coefficients in matrix Q
318 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4
319 596 / _f(2025), -398 / _f(945), 22 / _f(45), -1 / _f(3),
320 1543 / _f(4725), -118 / _f(315), 1 / _f(5),
321 152 / _f(945), -17 / _f(315),
322 5 / _f(252)),
323 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5
324 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),
325 -1 / _f(3),
326 -24562 / _f(155925), 1543 / _f(4725), -118 / _f(315), 1 / _f(5),
327 -38068 / _f(155925), 152 / _f(945), -17 / _f(315),
328 -752 / _f(10395), 5 / _f(252),
329 -101 / _f(17325)),
330 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6
331 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),
332 -398 / _f(945), 22 / _f(45), -1 / _f(3),
333 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),
334 -118 / _f(315), 1 / _f(5),
335 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),
336 -17 / _f(315),
337 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),
338 62464 / _f(2027025), -101 / _f(17325),
339 11537 / _f(4054050)),
340 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7
341 -565017322 / _f(1915538625), 138734126 / _f(638512875),
342 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),
343 -1 / _f(3),
344 -1969276 / _f(58046625), 17749373 / _f(425675250), -24562 / _f(155925),
345 1543 / _f(4725), -118 / _f(315), 1 / _f(5),
346 -58573784 / _f(638512875), 1882432 / _f(8513505), -38068 / _f(155925),
347 152 / _f(945), -17 / _f(315),
348 -6975184 / _f(42567525), 268864 / _f(2027025), -752 / _f(10395),
349 5 / _f(252),
350 -112832 / _f(1447875), 62464 / _f(2027025), -101 / _f(17325),
351 -4096 / _f(289575), 11537 / _f(4054050),
352 -311 / _f(525525)),
353 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8
354 188270561816 / _f(488462349375), -565017322 / _f(1915538625),
355 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),
356 -398 / _f(945), 22 / _f(45), -1 / _f(3),
357 2332829602 / _f(23260111875), -1969276 / _f(58046625),
358 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),
359 -118 / _f(315), 1 / _f(5),
360 -41570288 / _f(930404475), -58573784 / _f(638512875),
361 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),
362 -17 / _f(315),
363 1538774036 / _f(10854718875), -6975184 / _f(42567525),
364 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),
365 436821248 / _f(3618239625), -112832 / _f(1447875),
366 62464 / _f(2027025), -101 / _f(17325),
367 3059776 / _f(80405325), -4096 / _f(289575), 11537 / _f(4054050),
368 4193792 / _f(723647925), -311 / _f(525525),
369 1097653 / _f(1929727800))
370}
371del _f, _u, _Ufloats
373__all__ += _ALL_DOCS(Caps)
375# **) MIT License
376#
377# Copyright (C) 2023-2024 -- mrJean1 at Gmail -- All Rights Reserved.
378#
379# Permission is hereby granted, free of charge, to any person obtaining a
380# copy of this software and associated documentation files (the "Software"),
381# to deal in the Software without restriction, including without limitation
382# the rights to use, copy, modify, merge, publish, distribute, sublicense,
383# and/or sell copies of the Software, and to permit persons to whom the
384# Software is furnished to do so, subject to the following conditions:
385#
386# The above copyright notice and this permission notice shall be included
387# in all copies or substantial portions of the Software.
388#
389# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
390# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
391# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
392# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
393# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
394# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
395# OTHER DEALINGS IN THE SOFTWARE.