Merge pull request #421 from ghiggi/bugfix-arc-inplace-modification

Fix inplace modification occuring in Arc.intersections

pyresample/spherical.py

@@ -26,6 +26,13 @@
 
 logger = logging.getLogger(__name__)
 
+EPSILON = 0.0000001
+
+
+def _unwrap_radians(val, mod=np.pi):
+    """Put *val* between -*mod* and *mod*."""
+    return (val + mod) % (2 * mod) - mod
+
 
 class SCoordinate(object):
     """Spherical coordinates.
@@ -35,7 +42,7 @@ class SCoordinate(object):
     """
 
     def __init__(self, lon, lat):
-        self.lon = lon
+        self.lon = float(_unwrap_radians(lon))
         self.lat = lat
 
     def cross2cart(self, point):
@@ -171,14 +178,6 @@ def to_spherical(self):
                            np.arcsin(self.cart[2]))
 
 
-EPSILON = 0.0000001
-
-
-def modpi(val, mod=np.pi):
-    """Put *val* between -*mod* and *mod*."""
-    return (val + mod) % (2 * mod) - mod
-
-
 class Arc(object):
     """An arc of the great circle between two points."""
 
@@ -252,25 +251,31 @@ def intersections(self, other_arc):
 
         From http://williams.best.vwh.net/intersect.htm
         """
+        end_lon = self.end.lon
+        other_end_lon = other_arc.end.lon
+
         if self.end.lon - self.start.lon > np.pi:
-            self.end.lon -= 2 * np.pi
+            end_lon -= 2 * np.pi
         if other_arc.end.lon - other_arc.start.lon > np.pi:
-            other_arc.end.lon -= 2 * np.pi
+            other_end_lon -= 2 * np.pi
         if self.end.lon - self.start.lon < -np.pi:
-            self.end.lon += 2 * np.pi
+            end_lon += 2 * np.pi
         if other_arc.end.lon - other_arc.start.lon < -np.pi:
-            other_arc.end.lon += 2 * np.pi
+            other_end_lon += 2 * np.pi
 
-        ea_ = self.start.cross2cart(self.end).normalize()
-        eb_ = other_arc.start.cross2cart(other_arc.end).normalize()
+        end_point = SCoordinate(end_lon, self.end.lat)
+        other_end_point = SCoordinate(other_end_lon, other_arc.end.lat)
+
+        ea_ = self.start.cross2cart(end_point).normalize()
+        eb_ = other_arc.start.cross2cart(other_end_point).normalize()
 
         cross = ea_.cross(eb_)
         lat = np.arctan2(cross.cart[2],
                          np.sqrt(cross.cart[0] ** 2 + cross.cart[1] ** 2))
         lon = np.arctan2(cross.cart[1], cross.cart[0])
 
         return (SCoordinate(lon, lat),
-                SCoordinate(modpi(lon + np.pi), -lat))
+                SCoordinate(_unwrap_radians(lon + np.pi), -lat))
 
     def intersects(self, other_arc):
         """Check if the current arc and the *other_arc* intersect.
@@ -356,7 +361,7 @@ def __init__(self, vertices, radius=1):
             radius (optional, number): Radius of spherical planet.
         """
         self.vertices = vertices.astype(np.float64, copy=False)
-        self.lon = self.vertices[:, 0]
+        self.lon = _unwrap_radians(self.vertices[:, 0])
         self.lat = self.vertices[:, 1]
         self.radius = radius
         self.cvertices = np.array([np.cos(self.lat) * np.cos(self.lon),

pyresample/spherical_geometry.py

@@ -34,7 +34,11 @@
 
 
 class Coordinate(object):
-    """Point on earth in terms of lat and lon."""
+    """Point on earth in terms of lat and lon.
+
+    It expects lon,lat in degrees
+    But self.lat and self.lon are returned in radians !
+    """
 
     lat = None
     lon = None
@@ -236,25 +240,30 @@ def angle(self, other_arc, snap=True):
 
     def intersections(self, other_arc):
         """Get the two intersections of the greats circles defined by the current arc and *other_arc*."""
+        end_lon = self.end.lon
+        other_end_lon = other_arc.end.lon
+
         if self.end.lon - self.start.lon > math.pi:
-            self.end.lon -= 2 * math.pi
+            end_lon -= 2 * math.pi
         if other_arc.end.lon - other_arc.start.lon > math.pi:
-            other_arc.end.lon -= 2 * math.pi
+            other_end_lon -= 2 * math.pi
         if self.end.lon - self.start.lon < -math.pi:
-            self.end.lon += 2 * math.pi
+            end_lon += 2 * math.pi
         if other_arc.end.lon - other_arc.start.lon < -math.pi:
-            other_arc.end.lon += 2 * math.pi
+            other_end_lon += 2 * math.pi
+
+        end_point = Coordinate(math.degrees(modpi(end_lon)), math.degrees(self.end.lat))
+        other_end_point = Coordinate(math.degrees(modpi(other_end_lon)), math.degrees(other_arc.end.lat))
 
-        ea_ = self.start.cross2cart(self.end).normalize()
-        eb_ = other_arc.start.cross2cart(other_arc.end).normalize()
+        ea_ = self.start.cross2cart(end_point).normalize()
+        eb_ = other_arc.start.cross2cart(other_end_point).normalize()
 
         cross = ea_.cross(eb_)
         lat = math.atan2(cross.z__, math.sqrt(cross.x__ ** 2 + cross.y__ ** 2))
         lon = math.atan2(-cross.y__, cross.x__)
 
         return (Coordinate(math.degrees(lon), math.degrees(lat)),
-                Coordinate(math.degrees(modpi(lon + math.pi)),
-                           math.degrees(-lat)))
+                Coordinate(math.degrees(modpi(lon + math.pi)), math.degrees(-lat)))
 
     def intersects(self, other_arc):
         """Determine if this arc and another arc intersect.