Update geographiclib library

ext-libs/geographiclib/__init__.py

@@ -1,7 +1,7 @@
 """geographiclib: geodesic routines from GeographicLib"""
 
-__version_info__ = (1, 49, 0)
+__version_info__ = (1, 52, 0)
 """GeographicLib version as a tuple"""
 
-__version__ = "1.49"
+__version__ = "1.52"
 """GeographicLib version as a string"""

ext-libs/geographiclib/accumulator.py

@@ -7,8 +7,8 @@
 #
 #    https://geographiclib.sourceforge.io/html/annotated.html
 #
-# Copyright (c) Charles Karney (2011) <[email protected]> and licensed under
-# the MIT/X11 License.  For more information, see
+# Copyright (c) Charles Karney (2011-2019) <[email protected]> and
+# licensed under the MIT/X11 License.  For more information, see
 # https://geographiclib.sourceforge.io/
 ######################################################################
 
@@ -80,3 +80,8 @@ def Negate(self):
     """Negate sum"""
     self._s *= -1
     self._t *= -1
+
+  def Remainder(self, y):
+    """Remainder on division by y"""
+    self._s = Math.remainder(self._s, y)
+    self.Add(0.0)

ext-libs/geographiclib/geodesic.py

@@ -72,7 +72,7 @@
 #    https://doi.org/10.1007/s00190-012-0578-z
 #    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
 #
-# Copyright (c) Charles Karney (2011-2017) <[email protected]> and licensed
+# Copyright (c) Charles Karney (2011-2021) <[email protected]> and licensed
 # under the MIT/X11 License.  For more information, see
 # https://geographiclib.sourceforge.io/
 ######################################################################
@@ -810,7 +810,9 @@ def _GenInverse(self, lat1, lon1, lat2, lon2, outmask):
       # In fact, we will have sig12 > pi/2 for meridional geodesic which is
       # not a shortest path.
       if sig12 < 1 or m12x >= 0:
-        if sig12 < 3 * Geodesic.tiny_:
+        if (sig12 < 3 * Geodesic.tiny_ or
+            # Prevent negative s12 or m12 for short lines
+            (sig12 < Geodesic.tol0_ and (s12x < 0 or m12x < 0))):
           sig12 = m12x = s12x = 0.0
         m12x *= self._b
         s12x *= self._b
@@ -881,7 +883,6 @@ def _GenInverse(self, lat1, lon1, lat2, lon2, outmask):
              sbet1, cbet1, dn1, sbet2, cbet2, dn2,
              salp1, calp1, slam12, clam12, numit < Geodesic.maxit1_,
              C1a, C2a, C3a)
-          # 2 * tol0 is approximately 1 ulp for a number in [0, pi].
           # Reversed test to allow escape with NaNs
           if tripb or not (abs(v) >= (8 if tripn else 1) * Geodesic.tol0_):
             break

ext-libs/geographiclib/geodesicline.py

@@ -51,7 +51,7 @@
 #    https://doi.org/10.1007/s00190-012-0578-z
 #    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
 #
-# Copyright (c) Charles Karney (2011-2016) <[email protected]> and licensed
+# Copyright (c) Charles Karney (2011-2019) <[email protected]> and licensed
 # under the MIT/X11 License.  For more information, see
 # https://geographiclib.sourceforge.io/
 ######################################################################
@@ -112,7 +112,7 @@ def __init__(self, geod, lat1, lon1, azi1,
       """the cosine of the azimuth at the first point (readonly)"""
 
     # real cbet1, sbet1
-    sbet1, cbet1 = Math.sincosd(Math.AngRound(lat1)); sbet1 *= self._f1
+    sbet1, cbet1 = Math.sincosd(Math.AngRound(self.lat1)); sbet1 *= self._f1
     # Ensure cbet1 = +epsilon at poles
     sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
     self._dn1 = math.sqrt(1 + geod._ep2 * Math.sq(sbet1))
@@ -205,6 +205,7 @@ def _GenPosition(self, arcmode, s12_a12, outmask):
     else:
       # Interpret s12_a12 as distance
       tau12 = s12_a12 / (self._b * (1 + self._A1m1))
+      tau12 = tau12 if Math.isfinite(tau12) else Math.nan
       s = math.sin(tau12); c = math.cos(tau12)
       # tau2 = tau1 + tau12
       B12 = - Geodesic._SinCosSeries(True,

ext-libs/geographiclib/geomath.py

@@ -6,7 +6,7 @@
 #
 #    https://geographiclib.sourceforge.io/html/annotated.html
 #
-# Copyright (c) Charles Karney (2011-2017) <[email protected]> and
+# Copyright (c) Charles Karney (2011-2021) <[email protected]> and
 # licensed under the MIT/X11 License.  For more information, see
 # https://geographiclib.sourceforge.io/
 ######################################################################
@@ -44,7 +44,7 @@ def cbrt(x):
     """Real cube root of a number"""
 
     y = math.pow(abs(x), 1/3.0)
-    return y if x >= 0 else -y
+    return y if x > 0 else (-y if x < 0 else x)
   cbrt = staticmethod(cbrt)
 
   def log1p(x):
@@ -70,7 +70,7 @@ def atanh(x):
 
     y = abs(x)                  # Enforce odd parity
     y = Math.log1p(2 * y/(1 - y))/2
-    return -y if x < 0 else y
+    return y if x > 0 else (-y if x < 0 else x)
   atanh = staticmethod(atanh)
 
   def copysign(x, y):
@@ -84,7 +84,13 @@ def copysign(x, y):
 
   def norm(x, y):
     """Private: Normalize a two-vector."""
-    r = math.hypot(x, y)
+    r = (math.sqrt(Math.sq(x) + Math.sq(y))
+         # hypot is inaccurate for 3.[89].  Problem reported by agdhruv
+         # https://github.com/geopy/geopy/issues/466 ; see
+         # https://bugs.python.org/issue43088
+         # Visual Studio 2015 32-bit has a similar problem.
+         if (3, 8) <= sys.version_info < (3, 10)
+         else math.hypot(x, y))
     return x/r, y/r
   norm = staticmethod(norm)
 
@@ -126,16 +132,22 @@ def AngRound(x):
     return 0.0 if x == 0 else (-y if x < 0 else y)
   AngRound = staticmethod(AngRound)
 
-  def AngNormalize(x):
-    """reduce angle to (-180,180]"""
-
-    y = math.fmod(x, 360)
+  def remainder(x, y):
+    """remainder of x/y in the range [-y/2, y/2]."""
+    z = math.fmod(x, y) if Math.isfinite(x) else Math.nan
     # On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0
     # This fixes this bug.  See also Math::AngNormalize in the C++ library.
     # sincosd has a similar fix.
-    y = x if x == 0 else y
-    return (y + 360 if y <= -180 else
-            (y if y <= 180 else y - 360))
+    z = x if x == 0 else z
+    return (z + y if z < -y/2 else
+            (z if z < y/2 else z -y))
+  remainder = staticmethod(remainder)
+
+  def AngNormalize(x):
+    """reduce angle to (-180,180]"""
+
+    y = Math.remainder(x, 360)
+    return 180 if y == -180 else y
   AngNormalize = staticmethod(AngNormalize)
 
   def LatFix(x):
@@ -155,8 +167,8 @@ def AngDiff(x, y):
   def sincosd(x):
     """Compute sine and cosine of x in degrees."""
 
-    r = math.fmod(x, 360)
-    q = Math.nan if Math.isnan(r) else int(math.floor(r / 90 + 0.5))
+    r = math.fmod(x, 360) if Math.isfinite(x) else Math.nan
+    q = 0 if Math.isnan(r) else int(round(r / 90))
     r -= 90 * q; r = math.radians(r)
     s = math.sin(r); c = math.cos(r)
     q = q % 4

ext-libs/geographiclib/polygonarea.py

@@ -39,7 +39,7 @@
 #    https://doi.org/10.1007/s00190-012-0578-z
 #    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
 #
-# Copyright (c) Charles Karney (2011-2017) <[email protected]> and licensed
+# Copyright (c) Charles Karney (2011-2019) <[email protected]> and licensed
 # under the MIT/X11 License.  For more information, see
 # https://geographiclib.sourceforge.io/
 ######################################################################
@@ -67,14 +67,60 @@ def _transit(lon1, lon2):
   def _transitdirect(lon1, lon2):
     """Count crossings of prime meridian for AddEdge."""
     # We want to compute exactly
-    #   int(floor(lon2 / 360)) - int(floor(lon1 / 360))
+    #   int(ceil(lon2 / 360)) - int(ceil(lon1 / 360))
     # Since we only need the parity of the result we can use std::remquo but
     # this is buggy with g++ 4.8.3 and requires C++11.  So instead we do
     lon1 = math.fmod(lon1, 720.0); lon2 = math.fmod(lon2, 720.0)
-    return ( (0 if ((lon2 >= 0 and lon2 < 360) or lon2 < -360) else 1) -
-             (0 if ((lon1 >= 0 and lon1 < 360) or lon1 < -360) else 1) )
+    return ( (1 if ((lon2 <= 0 and lon2 > -360) or lon2 > 360) else 0) -
+             (1 if ((lon1 <= 0 and lon1 > -360) or lon1 > 360) else 0) )
   _transitdirect = staticmethod(_transitdirect)
 
+  def _areareduceA(area, area0, crossings, reverse, sign):
+    """Reduce accumulator area to allowed range."""
+    area.Remainder(area0)
+    if crossings & 1:
+      area.Add( (1 if area.Sum() < 0 else -1) * area0/2 )
+    # area is with the clockwise sense.  If !reverse convert to
+    # counter-clockwise convention.
+    if not reverse: area.Negate()
+    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+    if sign:
+      if area.Sum() > area0/2:
+        area.Add( -area0 )
+      elif area.Sum() <= -area0/2:
+        area.Add(  area0 )
+    else:
+      if area.Sum() >= area0:
+        area.Add( -area0 )
+      elif area.Sum() < 0:
+        area.Add(  area0 )
+
+    return 0.0 + area.Sum()
+  _areareduceA = staticmethod(_areareduceA)
+
+  def _areareduceB(area, area0, crossings, reverse, sign):
+    """Reduce double area to allowed range."""
+    area = Math.remainder(area, area0)
+    if crossings & 1:
+      area += (1 if area < 0 else -1) * area0/2
+    # area is with the clockwise sense.  If !reverse convert to
+    # counter-clockwise convention.
+    if not reverse: area *= -1
+    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
+    if sign:
+      if area > area0/2:
+        area -= area0
+      elif area <= -area0/2:
+        area += area0
+    else:
+      if area >= area0:
+        area -= area0
+      elif area < 0:
+        area += area0
+
+    return 0.0 + area
+  _areareduceB = staticmethod(_areareduceB)
+
   def __init__(self, earth, polyline = False):
     """Construct a PolygonArea object
 
@@ -169,7 +215,12 @@ def Compute(self, reverse = False, sign = True):
       the area for the rest of the earth
     :return: a tuple of number, perimeter (meters), area (meters^2)
 
-    If the object is a polygon (and not a polygon), the perimeter
+    Arbitrarily complex polygons are allowed.  In the case of
+    self-intersecting polygons the area is accumulated "algebraically",
+    e.g., the areas of the 2 loops in a figure-8 polygon will partially
+    cancel.
+
+    If the object is a polygon (and not a polyline), the perimeter
     includes the length of a final edge connecting the current point to
     the initial point.  If the object is a polyline, then area is nan.
 
@@ -192,24 +243,8 @@ def Compute(self, reverse = False, sign = True):
     tempsum = Accumulator(self._areasum)
     tempsum.Add(S12)
     crossings = self._crossings + PolygonArea._transit(self.lon1, self._lon0)
-    if crossings & 1:
-      tempsum.Add( (1 if tempsum.Sum() < 0 else -1) * self.area0/2 )
-    # area is with the clockwise sense.  If !reverse convert to
-    # counter-clockwise convention.
-    if not reverse: tempsum.Negate()
-    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
-    if sign:
-      if tempsum.Sum() > self.area0/2:
-        tempsum.Add( -self.area0 )
-      elif tempsum.Sum() <= -self.area0/2:
-        tempsum.Add(  self.area0 )
-    else:
-      if tempsum.Sum() >= self.area0:
-        tempsum.Add( -self.area0 )
-      elif tempsum.Sum() < 0:
-        tempsum.Add(  self.area0 )
-
-    area = 0.0 + tempsum.Sum()
+    area = PolygonArea._areareduceA(tempsum, self.area0, crossings,
+                                    reverse, sign)
     return self.num, perimeter, area
 
   # return number, perimeter, area
@@ -249,24 +284,8 @@ def TestPoint(self, lat, lon, reverse = False, sign = True):
     if self.polyline:
       return num, perimeter, area
 
-    if crossings & 1:
-      tempsum += (1 if tempsum < 0 else -1) * self.area0/2
-    # area is with the clockwise sense.  If !reverse convert to
-    # counter-clockwise convention.
-    if not reverse: tempsum *= -1
-    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
-    if sign:
-      if tempsum > self.area0/2:
-        tempsum -= self.area0
-      elif tempsum <= -self.area0/2:
-        tempsum += self.area0
-    else:
-      if tempsum >= self.area0:
-        tempsum -= self.area0
-      elif tempsum < 0:
-        tempsum += self.area0
-
-    area = 0.0 + tempsum
+    area = PolygonArea._areareduceB(tempsum, self.area0, crossings,
+                                    reverse, sign)
     return num, perimeter, area
 
   # return num, perimeter, area
@@ -303,22 +322,6 @@ def TestEdge(self, azi, s, reverse = False, sign = True):
     tempsum += S12
     crossings += PolygonArea._transit(lon, self._lon0)
 
-    if crossings & 1:
-      tempsum += (1 if tempsum < 0 else -1) * self.area0/2
-    # area is with the clockwise sense.  If !reverse convert to
-    # counter-clockwise convention.
-    if not reverse: tempsum *= -1
-    # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
-    if sign:
-      if tempsum > self.area0/2:
-        tempsum -= self.area0
-      elif tempsum <= -self.area0/2:
-        tempsum += self.area0
-    else:
-      if tempsum >= self.area0:
-        tempsum -= self.area0
-      elif tempsum < 0:
-        tempsum += self.area0
-
-    area = 0.0 + tempsum
+    area = PolygonArea._areareduceB(tempsum, self.area0, crossings,
+                                    reverse, sign)
     return num, perimeter, area