Back ported changes from astropy and made python 2 to 3 modifications

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci.sphere/trunk@38800 fe389314-cf27-0410-b35b-8c050e845b92

Former-commit-id: c12c8fef5503caab2e2e57caac8f040597326589

7 files changed, 582 insertions, 426 deletions

diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py
index 01c582a..55d5835 100644
--- a/stsci/sphere/graph.py
+++ b/stsci/sphere/graph.py
@@ -35,7 +35,6 @@
 """
 This contains the code that does the actual unioning of regions.
 """
-# TODO: Weak references for memory management problems?
 from __future__ import absolute_import, division, unicode_literals, print_function
 
 # STDLIB
@@ -594,7 +593,7 @@ class Graph:
         #
 
         edges = list(self._edges)
-        for i in xrange(len(edges)):
+        for i in range(len(edges)):
             AB = edges[i]
             if AB not in self._edges:
                 continue
@@ -703,8 +702,7 @@ class Graph:
                 # Delete CD, and push the new edges to the
                 # front so they will be tested for intersection
                 # against all remaining edges.
-                del edges[j]  # CD
-                edges = new_edges + edges
+                edges = new_edges + edges[:j] + edges[j+1:]
                 new_starts, new_ends = self._get_edge_points(new_edges)
                 starts = np.vstack(
                     (new_starts, starts[:j], starts[j+1:]))
@@ -783,7 +781,7 @@ class Graph:
             max_ary = 0
             for node in self._nodes:
                 max_ary = max(len(node._edges), max_ary)
-            order = range(max_ary + 1, 2, -1)
+            order = list(range(max_ary + 1, 2, -1))
         else:
             order = [3]
 
@@ -814,156 +812,41 @@ class Graph:
             else:
                 break
 
-    @staticmethod
-    def _fast_area(points):
-        """
-        Calculates an approximate area of a polygon.  Unlike the area
-        calculation in `SphericalPolygon`, this does not interpolate
-        along the arcs, so is not useful for getting an actual area,
-        but only whether the area is positive or negative, telling us
-        whether the points proceed clockwise or counter-clockwise
-        respectively.
-        """
-        # Rotate polygon so that center of polygon is at north pole
-        centroid = np.mean(points[:-1], axis=0)
-        centroid = vector.normalize_vector(centroid)
-        points = points - (centroid + np.array([0, 0, 1]))
-        vector.normalize_vector(points, output=points)
-
-        XY = vector.equal_area_proj(points)
-        X = XY[..., 0]
-        Y = XY[..., 1]
-
-        return np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1])
-
     def _trace(self):
         """
         Given a graph that has had cutlines removed and all
         intersections found, traces it to find a resulting single
         polygon.
         """
+        from . import polygon as mpolygon
+
         polygons = []
         edges = set(self._edges)  # copy
         while len(edges):
-            polygon = []
+            points = []
             edge = edges.pop()
             start_node = node = edge._nodes[0]
             while True:
                 # TODO: Do we need this if clause any more?
-                if len(polygon):
-                    if not np.array_equal(polygon[-1], node._point):
-                        polygon.append(node._point)
+                if len(points):
+                    if not np.array_equal(points[-1], node._point):
+                        points.append(node._point)
                 else:
-                    polygon.append(node._point)
+                    points.append(node._point)
                 edge = node.follow(edge)
                 edges.discard(edge)
                 node = edge.follow(node)
                 if node is start_node:
-                    polygon.append(node._point)
+                    points.append(node._point)
                     break
 
-            polygon = np.asarray(polygon)
-            area = self._fast_area(polygon)
-            if area < 0:
-                # Reverse the points
-                polygon = polygon[::-1, ...]
-
+            polygon = mpolygon._SingleSphericalPolygon(points)
+            area = polygon.area()
+            reverse_polygon = mpolygon._SingleSphericalPolygon(points[::-1])
+            reverse_area = reverse_polygon.area()
+            if reverse_area < area:
+                polygon = reverse_polygon
             polygons.append(polygon)
 
-        if len(polygons) == 1:
-            return polygons[0]
-        elif len(polygons) == 0:
-            return []
-        else:
-            return self._join(polygons)
-
-    def _join(self, polygons):
-        """
-        If the graph is disjunct, joins the parts with cutlines.
-
-        The closest nodes between each pair that don't intersect
-        any other edges are used as cutlines.
-
-        TODO: This is not optimal, because the closest distance
-        between two polygons may not in fact be between two vertices,
-        but may be somewhere along an edge.
-        """
-        def do_join(polygons):
-            all_polygons = polygons[:]
-
-            skipped = 0
-
-            polyA = polygons.pop(0)
-            while len(polygons):
-                polyB = polygons.pop(0)
-
-                # If fewer than 3 edges, it's not a polygon,
-                # just throw it out
-                if len(polyB) < 4:
-                    continue
-
-                # Find the closest set of vertices between polyA and
-                # polyB that don't cross any of the edges in *any* of
-                # the polygons
-                closest = np.inf
-                closest_pair_idx = (None, None)
-                for a in xrange(len(polyA) - 1):
-                    A = polyA[a]
-                    distances = great_circle_arc.length(A, polyB[:-1])
-                    b = np.argmin(distances)
-                    distance = distances[b]
-                    if distance < closest:
-                        B = polyB[b]
-                        # Does this candidate line cross other edges?
-                        crosses = False
-                        for poly in all_polygons:
-                            if np.any(
-                                great_circle_arc.intersects(
-                                    A, B, poly[:-1], poly[1:])):
-                                crosses = True
-                                break
-                        if not crosses:
-                            closest = distance
-                            closest_pair_idx = (a, b)
-
-                if not np.isfinite(closest):
-                    # We didn't find a pair of points that don't cross
-                    # something else, so we want to try to join another
-                    # polygon.  Defer the current polygon until later.
-                    if len(polygons) in (0, skipped):
-                        return None
-                    polygons.append(polyB)
-                    skipped += 1
-                else:
-                    # Splice the two polygons together using a cut
-                    # line
-                    a, b = closest_pair_idx
-                    new_poly = np.vstack((
-                        # polyA up to and including the cut point
-                        polyA[:a+1],
-
-                        # polyB starting with the cut point and
-                        # wrapping around back to the cut point.
-                        # Ignore the last point in polyB, because it
-                        # is the same as the first
-                        np.roll(polyB[:-1], -b, 0),
-
-                        # The cut point on polyB
-                        polyB[b:b+1],
-
-                        # the rest of polyA, starting with the cut
-                        # point
-                        polyA[a:]
-                        ))
-
-                    skipped = 0
-                    polyA = new_poly
-
-            return polyA
-
-        for permutation in itertools.permutations(polygons):
-            poly = do_join(list(permutation))
-            if poly is not None:
-                return poly
-
-        raise RuntimeError("Could not find cut points")
+        return mpolygon.SphericalPolygon(polygons)
+
diff --git a/stsci/sphere/great_circle_arc.py b/stsci/sphere/great_circle_arc.py
index 0b47d89..e8d421b 100644
--- a/stsci/sphere/great_circle_arc.py
+++ b/stsci/sphere/great_circle_arc.py
@@ -59,7 +59,7 @@ else:
     from numpy.core.umath_tests import inner1d
 
 
-___all__ = ['angle', 'intersection', 'intersects', 'intersects_point', 
+__all__ = ['angle', 'intersection', 'intersects', 'intersects_point', 
            'length', 'midpoint', 'interpolate']
 
 
@@ -285,9 +285,6 @@ def intersects_point(A, B, C):
     """
     Returns True if point C is along the great circle arc *AB*.
     """
-    if HAS_C_UFUNCS:
-        return math_util.intersects_point(A, B, C)
-
     total_length = length(A, B)
     left_length = length(A, C)
     right_length = length(C, B)
@@ -301,8 +298,6 @@ def angle(A, B, C, degrees=True):
     """
     Returns the angle at *B* between *AB* and *BC*.
 
-    This always returns the shortest angle < π.
-
     Parameters
     ----------
     A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples
@@ -333,8 +328,12 @@ def angle(A, B, C, degrees=True):
     ABX = _cross_and_normalize(B, ABX)
     BCX = _fast_cross(C, B)
     BCX = _cross_and_normalize(B, BCX)
+    X = _cross_and_normalize(ABX, BCX)
+    diff = inner1d(B, X)
+    inner = inner1d(ABX, BCX)
     with np.errstate(invalid='ignore'):
-        angle = np.arccos(inner1d(ABX, BCX))
+        angle = np.arccos(inner)
+    angle = np.where(diff < 0.0, (2.0 * np.pi) - angle, angle)
 
     if degrees:
         angle = np.rad2deg(angle)
@@ -392,7 +391,7 @@ def interpolate(A, B, steps=50):
 
         \frac{\sin((1 - t)\Omega)}{\sin \Omega}A + \frac{\sin(t \Omega)}{\sin \Omega}B
     """
-    steps = max(steps, 2)
+    steps = int(max(steps, 2))
     t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1))
 
     omega = length(A, B, degrees=False)
diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py
index 5167171..449bd79 100644
--- a/stsci/sphere/polygon.py
+++ b/stsci/sphere/polygon.py
@@ -51,7 +51,7 @@ from . import vector
 __all__ = ['SphericalPolygon']
 
 
-class SphericalPolygon(object):
+class _SingleSphericalPolygon(object):
     r"""
     Polygons are represented by both a set of points (in Cartesian
     (*x*, *y*, *z*) normalized on the unit sphere), and an inside
@@ -91,7 +91,7 @@ class SphericalPolygon(object):
         self._points = np.asanyarray(points)
 
         if inside is None:
-            self._inside = np.mean(points[:-1], axis=0)
+            self._inside = self._find_new_inside(points)
         else:
             self._inside = np.asanyarray(inside)
 
@@ -100,13 +100,12 @@ class SphericalPolygon(object):
     def __copy__(self):
         return deepcopy(self)
 
+    copy = __copy__
+
     def __repr__(self):
         return '%s(%r, %r)' % (self.__class__.__name__,
                                self.points, self.inside)
 
-    def copy(self):
-        return self.__class__(self._points.copy(), self._inside.copy())
-
     @property
     def points(self):
         """
@@ -123,6 +122,13 @@ class SphericalPolygon(object):
         """
         return self._inside
 
+    def iter_polygons_flat(self):
+        """
+        Iterate over all base polygons that make up this multi-polygon
+        set.
+        """
+        yield self
+
     def to_radec(self):
         """
         Convert `SphericalPolygon` footprint to RA and DEC.
@@ -141,7 +147,7 @@ class SphericalPolygon(object):
     @classmethod
     def from_radec(cls, ra, dec, center=None, degrees=True):
         r"""
-        Create a new `SphericalPolygon` from a list of (*ra*, *dec*)
+        Create a new `_SingleSphericalPolygon` from a list of (*ra*, *dec*)
         points.
 
         Parameters
@@ -176,12 +182,12 @@ class SphericalPolygon(object):
         else:
             center = vector.radec_to_vector(*center, degrees=degrees)
 
-        return cls(np.dstack((x, y, z))[0], center)
+        return cls(points, center)
 
     @classmethod
     def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0):
         r"""
-        Create a new `SphericalPolygon` from a cone (otherwise known
+        Create a new `_SingleSphericalPolygon` from a cone (otherwise known
         as a "small circle") defined using (*ra*, *dec*, *radius*).
 
         The cone is not represented as an ideal circle on the sphere,
@@ -214,13 +220,14 @@ class SphericalPolygon(object):
 
         # Get an arbitrary perpendicular vector.  This be be obtained
         # by crossing (u, v, w) with any unit vector that is not itself.
-        which_min = np.argmin([u, v, w])
+        which_min = np.argmin([u*u, v*v, w*w])
         if which_min == 0:
             perp = np.cross([u, v, w], [1., 0., 0.])
         elif which_min == 1:
             perp = np.cross([u, v, w], [0., 1., 0.])
         else:
             perp = np.cross([u, v, w], [0., 0., 1.])
+        perp = vector.normalize_vector(perp)
 
         # Rotate by radius around the perpendicular vector to get the
         # "pen"
@@ -301,119 +308,18 @@ class SphericalPolygon(object):
 
         return cls(np.dstack((x, y, z))[0], (xc, yc, zc))
 
-    def _unique_points(self):
-        """
-        Return a copy of `points` with duplicates removed.
-        Order is preserved.
-
-        .. note:: Output cannot be used to build a new
-            polygon.
-        """
-        val = []
-        for p in self.points:
-            v = tuple(p)
-            if v not in val:
-                val.append(v)
-        return np.array(val)
-
-    def _sorted_points(self, preserve_order=True, unique=False):
-        """
-        Return a copy of `points` sorted such that smallest
-        (*x*, *y*, *z*) is on top.
-
-        .. note:: Output cannot be used to build a new
-            polygon.
-
-        Parameters
-        ----------
-        preserve_order : bool
-            Preserve original order? If `True`, polygon is
-            rotated around min point. If `False`, all points
-            are sorted in ascending order.
-
-        unique : bool
-            Exclude duplicates.
-        """
-        if len(self.points) == 0:
-            return []
-
-        if unique:
-            pts = self._unique_points()
-        else:
-            pts = self.points
-
-        idx = np.lexsort((pts[:,0], pts[:,1], pts[:,2]))
-
-        if preserve_order:
-            i_min = idx[0]
-            val = np.vstack([pts[i_min:], pts[:i_min]])
-        else:
-            val = pts[idx]
-
-        return val
-
-    def same_points_as(self, other, do_sort=True, thres=0.01):
-        """
-        Determines if this `SphericalPolygon` points are the same
-        as the other. Number of points and areas are also compared.
-
-        When `do_sort` is `True`, even when *self* and *other*
-        have same points, they might not be equivalent because
-        the order of the points defines the polygon.
-
-        Parameters
-        ----------
-        other : `SphericalPolygon`
-
-        do_sort : bool
-            Compare sorted unique points.
-
-        thres : float
-            Fraction of area to use in equality decision.
-
-        Returns
-        -------
-        is_eq : bool
-            `True` or `False`.
-        """
-        self_n = len(self.points)
-
-        if self_n != len(other.points):
-            return False
-
-        if self_n == 0:
-            return True
-
-        self_a = self.area()
-        is_same_limit = thres * self_a
-
-        if np.abs(self_a - other.area()) > is_same_limit:
-            return False
-
-        if do_sort:
-            self_pts  = self._sorted_points(preserve_order=False, unique=True)
-            other_pts = other._sorted_points(preserve_order=False, unique=True)
-        else:
-            self_pts  = self.points
-            other_pts = other.points
-
-        is_eq = True
-
-        for self_p, other_p in zip(self_pts, other_pts):
-            x_sum = 0.0
-
-            for a,b in zip(self_p, other_p):
-                x_sum += (a - b) ** 2
+    @classmethod
+    def _contains_point(cls, point, P, r):
+        point = np.asanyarray(point)
 
-            if np.sqrt(x_sum) > is_same_limit:
-                is_eq = False
-                break
+        intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point)
+        crossings = np.sum(intersects)
 
-        return is_eq
+        return (crossings % 2) == 0
 
     def contains_point(self, point):
         r"""
-        Determines if this `SphericalPolygon` contains a given point.
+        Determines if this `_SingleSphericalPolygon` contains a given point.
 
         Parameters
         ----------
@@ -425,14 +331,7 @@ class SphericalPolygon(object):
         contains : bool
             Returns `True` if the polygon contains the given *point*.
         """
-        P = self._points
-        r = self._inside
-        point = np.asanyarray(point)
-
-        intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point)
-        crossings = np.sum(intersects)
-
-        return (crossings % 2) == 0
+        return self._contains_point(point, self._points, self._inside)
 
     def intersects_poly(self, other):
         r"""
@@ -471,7 +370,7 @@ class SphericalPolygon(object):
                In this case, an edge from one polygon must cross an
                edge from the other polygon.
         """
-        assert isinstance(other, SphericalPolygon)
+        assert isinstance(other, _SingleSphericalPolygon)
 
         # The easy case is in which a point of one polygon is
         # contained in the other polygon.
@@ -520,59 +419,30 @@ class SphericalPolygon(object):
 
     def area(self):
         r"""
-        Returns the area of the polygon on the unit sphere.
+        Returns the area of the polygon on the unit sphere, in
+        steradians.
 
-        The algorithm is not able to compute the area of polygons
-        that are larger than half of the sphere.  Therefore, the
-        area will always be less than 2π.
+        The area is computed using a generalization of Girard's Theorem.
 
-        The area is computed by transforming the polygon to two
-        dimensions using the `Lambert azimuthal equal-area projection
-        <http://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_
+        if :math:`\theta` is the sum of the internal angles of the
+        polygon, and *n* is the number of vertices, the area is:
 
         .. math::
 
-            X = \sqrt{\frac{2}{1-z}}x
-
-        .. math::
-
-            Y = \sqrt{\frac{2}{1-z}}y
-
-        The individual great arc circle segments are interpolated
-        before doing the transformation so that the curves are not
-        straightened in the process.
-
-        It then uses a standard 2D algorithm to compute the area.
-
-        .. math::
-
-            A = \left| \sum^n_{i=0} X_i Y_{i+1} - X_{i+1}Y_i \right|
+            S = \theta - (n - 2) \pi
         """
         if len(self._points) < 3:
-            #return np.float64(0.0)
             return np.array(0.0)
 
-        points = self._points.copy()
-
-        # Rotate polygon so that center of polygon is at north pole
-        centroid = np.mean(points[:-1], axis=0)
-        centroid = vector.normalize_vector(centroid)
-        points = self._points - (centroid + np.array([0, 0, 1]))
-        vector.normalize_vector(points, output=points)
-
-        XYs = []
-        for A, B in zip(points[:-1], points[1:]):
-            length = great_circle_arc.length(A, B, degrees=True)
-            interp = great_circle_arc.interpolate(A, B, length * 4)
-            vector.normalize_vector(interp, output=interp)
-            XY = vector.equal_area_proj(interp)
-            XYs.append(XY)
-
-        XY = np.vstack(XYs)
-        X = XY[..., 0]
-        Y = XY[..., 1]
+        points = self._points
+        angles = np.hstack([
+            great_circle_arc.angle(
+                points[:-2], points[1:-1], points[2:], degrees=False),
+            great_circle_arc.angle(
+                points[-2], points[0], points[1], degrees=False)])
+        sum = np.sum(angles) - (len(angles) - 2) * np.pi
 
-        return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi)
+        return sum
 
     def union(self, other):
         """
@@ -607,46 +477,16 @@ class SphericalPolygon(object):
         """
         from . import graph
         if len(self._points) < 3:
-            return other.copy()
+            return SphericalPolygon([other.copy()])
         elif len(other._points) < 3:
-            return self.copy()
+            return SphericalPolygon([self.copy()])
 
         g = graph.Graph([self, other])
 
-        polygon = g.union()
-
-        return self.__class__(polygon, self._inside)
+        return g.union()
 
     @classmethod
-    def multi_union(cls, polygons):
-        """
-        Return a new `SphericalPolygon` that is the union of all of the
-        polygons in *polygons*.
-
-        Parameters
-        ----------
-        polygons : sequence of `SphericalPolygon`
-
-        Returns
-        -------
-        polygon : `SphericalPolygon` object
-
-        See also
-        --------
-        union
-        """
-        assert len(polygons)
-        for polygon in polygons:
-            assert isinstance(polygon, SphericalPolygon)
-
-        from . import graph
-
-        g = graph.Graph(polygons)
-        polygon = g.union()
-        return cls(polygon, polygons[0]._inside)
-
-    @staticmethod
-    def _find_new_inside(points, polygons):
+    def _find_new_inside(cls, points):
         """
         Finds an acceptable inside point inside of *points* that is
         also inside of *polygons*.  Used by the intersection
@@ -654,7 +494,7 @@ class SphericalPolygon(object):
         it requires existing polygons with known inside points.
         """
         if len(points) < 4:
-            return np.array([0, 0, 0])
+            return points[0]
 
         # Special case for a triangle
         if len(points) == 4:
@@ -662,21 +502,20 @@ class SphericalPolygon(object):
 
         for i in range(len(points) - 1):
             A = points[i]
-            # Skip the adjacent point, since it is by definition on
-            # the edge of the polygon, not potentially running through
-            # the middle.
-            for j in range(i + 2, len(points) - 1):
-                B = points[j]
-                C = great_circle_arc.midpoint(A, B)
-                in_all = True
-                for polygon in polygons:
-                    if not polygon.contains_point(C):
-                        in_all = False
+            B = points[i+1]
+            C = points[(i+2) % (len(points) - 1)]
+            angle = great_circle_arc.angle(A, B, C, degrees=False)
+            if angle <= np.pi * 2.0:
+                inside = great_circle_arc.midpoint(A, C)
+
+                for point in points[:-1]:
+                    if not cls._contains_point(point, points, inside):
                         break
-                if in_all:
-                    return C
+                else:
+                    return inside
 
-        raise RuntimeError("Suitable inside point could not be found")
+        # Fallback to the mean
+        return np.sum(points[:-1]) / (len(points) - 1)
 
     def intersection(self, other):
         """
@@ -708,20 +547,449 @@ class SphericalPolygon(object):
         For implementation details, see the :mod:`~sphere.graph`
         module.
         """
-        # if not self.intersects_poly(other):
-        #     return self.__class__([], [0, 0, 0])
-
         from . import graph
         if len(self._points) < 3 or len(other._points) < 3:
-            return self.__class__([], [0, 0, 0])
+            return SphericalPolygon([])
 
         g = graph.Graph([self, other])
 
-        polygon = g.intersection()
+        return g.intersection()
+
+    def overlap(self, other):
+        r"""
+        Returns the fraction of *self* that is overlapped by *other*.
+
+        Let *self* be *a* and *other* be *b*, then the overlap is
+        defined as:
+
+        .. math::
+
+            \frac{S_a}{S_{a \cap b}}
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        frac : float
+            The fraction of *self* that is overlapped by *other*.
+        """
+        s1 = self.area()
+        intersection = self.intersection(other)
+        s2 = intersection.area()
+        return s2 / s1
+
+    def draw(self, m, **plot_args):
+        """
+        Draws the polygon in a matplotlib.Basemap axes.
+
+        Parameters
+        ----------
+        m : Basemap axes object
+
+        **plot_args : Any plot arguments to pass to basemap
+        """
+        if not len(self._points):
+            return
+        if not len(plot_args):
+            plot_args = {'color': 'blue'}
+        points = self._points
+        if 'alpha' in plot_args:
+            del plot_args['alpha']
 
-        inside = self._find_new_inside(polygon, [self, other])
+        alpha = 1.0
+        for A, B in zip(points[0:-1], points[1:]):
+            length = great_circle_arc.length(A, B, degrees=True)
+            if not np.isfinite(length):
+                length = 2
+            interpolated = great_circle_arc.interpolate(A, B, length * 4)
+            ra, dec = vector.vector_to_radec(
+                interpolated[:, 0], interpolated[:, 1], interpolated[:, 2],
+                degrees=True)
+            for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]):
+                m.drawgreatcircle(r0, d0, r1, d1, alpha=alpha, **plot_args)
+            alpha -= 1.0 / len(points)
 
-        return self.__class__(polygon, inside)
+        ra, dec = vector.vector_to_radec(
+            *self._inside, degrees=True)
+        x, y = m(ra, dec)
+        m.scatter(x, y, 1, **plot_args)
+
+
+class SphericalPolygon(object):
+    r"""
+    Polygons are represented by both a set of points (in Cartesian
+    (*x*, *y*, *z*) normalized on the unit sphere), and an inside
+    point.  The inside point is necessary, because both the inside and
+    outside of the polygon are finite areas on the great sphere, and
+    therefore we need a way of specifying which is which.
+
+    This class contains a list of disjoint closed polygons.
+    """
+
+    def __init__(self, init, inside=None):
+        r"""
+        Parameters
+        ----------
+        init : object
+            May be either:
+               - A list of disjoint `SphericalPolygon` objects.
+
+               - An Nx3 array of (*x*, *y*, *z*) triples in Cartesian
+                 space.  These points define the boundary of the
+                 polygon.  It must be "closed", i.e., the last point
+                 is the same as the first.
+
+                 It may contain zero points, in which it defines the
+                 null polygon.  It may not contain one, two or three
+                 points.  Four points are needed to define a triangle,
+                 since the polygon must be closed.
+
+        inside : An (*x*, *y*, *z*) triple, optional
+            If *init* is an array of points, this point must be inside
+            the polygon.  If not provided, the mean of the points will
+            be used.
+        """
+        for polygon in init:
+            if not isinstance(polygon, (SphericalPolygon, _SingleSphericalPolygon)):
+                break
+        else:
+            self._polygons = tuple(init)
+            return
+
+        self._polygons = (_SingleSphericalPolygon(init, inside),)
+
+    def __copy__(self):
+        return deepcopy(self)
+
+    copy = __copy__
+
+    def iter_polygons_flat(self):
+        """
+        Iterate over all base polygons that make up this multi-polygon
+        set.
+        """
+        for polygon in self._polygons:
+            for subpolygon in polygon.iter_polygons_flat():
+                yield subpolygon
+
+    @property
+    def points(self):
+        """
+        The points defining the polygons.  It is an iterator over
+        disjoint closed polygons, where each element is an Nx3 array
+        of (*x*, *y*, *z*) vectors.  Each polygon is explicitly
+        closed, i.e., the first and last points are the same.
+        """
+        for polygon in self.iter_polygons_flat():
+            yield polygon.points
+
+    @property
+    def inside(self):
+        """
+        Iterate over the inside point of each of the polygons.
+        """
+        for polygon in self.iter_polygons_flat():
+            yield polygon.points
+
+    @property
+    def polygons(self):
+        """
+        Get a sequence of all of the subpolygons.  Each subpolygon may
+        itself have subpolygons.  To get a flattened sequence of all
+        base polygons, use `iter_polygons_flat`.
+        """
+        return self._polygons
+
+    def to_radec(self):
+        """
+        Convert the `SphericalPolygon` footprint to RA and DEC
+        coordinates.
+
+        Returns
+        -------
+        polyons : iterator
+            Each element in the iterator is a tuple of the form (*ra*,
+            *dec*), where each is an array of points.
+        """
+        for polygon in self.iter_polygons_flat():
+            yield polygon.to_radec()
+
+    @classmethod
+    def from_radec(cls, ra, dec, center=None, degrees=True):
+        r"""
+        Create a new `SphericalPolygon` from a list of (*ra*, *dec*)
+        points.
+
+        Parameters
+        ----------
+        ra, dec : 1-D arrays of the same length
+            The vertices of the polygon in right-ascension and
+            declination.  It must be \"closed\", i.e., that is, the
+            last point is the same as the first.
+
+        center : (*ra*, *dec*) pair, optional
+            A point inside of the polygon to define its inside.  If no
+            *center* point is provided, the mean of the polygon's
+            points in vector space will be used.  That approach may
+            not work for concave polygons.
+
+        degrees : bool, optional
+            If `True`, (default) *ra* and *dec* are in decimal degrees,
+            otherwise in radians.
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+        """
+        return cls([
+            _SingleSphericalPolygon.from_radec(
+                ra, dec, center=center, degrees=degrees)])
+
+    @classmethod
+    def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0):
+        r"""
+        Create a new `SphericalPolygon` from a cone (otherwise known
+        as a "small circle") defined using (*ra*, *dec*, *radius*).
+
+        The cone is not represented as an ideal circle on the sphere,
+        but as a series of great circle arcs.  The resolution of this
+        conversion can be controlled using the *steps* parameter.
+
+        Parameters
+        ----------
+        ra, dec : float scalars
+            This defines the center of the cone
+
+        radius : float scalar
+            The radius of the cone
+
+        degrees : bool, optional
+            If `True`, (default) *ra*, *dec* and *radius* are in
+            decimal degrees, otherwise in radians.
+
+        steps : int, optional
+            The number of steps to use when converting the small
+            circle to a polygon.
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+        """
+        return cls([
+            _SingleSphericalPolygon.from_cone(
+                ra, dec, radius, degrees=degrees, steps=steps)])
+
+    @classmethod
+    def from_wcs(cls, fitspath, steps=1, crval=None):
+        r"""
+        Create a new `SphericalPolygon` from the footprint of a FITS
+        WCS specification.
+
+        This method requires having `astropy <http://astropy.org>`__
+        installed.
+
+        Parameters
+        ----------
+        fitspath : path to a FITS file, `astropy.io.fits.Header`, or `astropy.wcs.WCS`
+            Refers to a FITS header containing a WCS specification.
+
+        steps : int, optional
+            The number of steps along each edge to convert into
+            polygon edges.
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+        """
+        return cls([
+            _SingleSphericalPolygon.from_wcs(
+                fitspath, steps=steps, crval=crval)])
+
+    def contains_point(self, point):
+        r"""
+        Determines if this `SphericalPolygon` contains a given point.
+
+        Parameters
+        ----------
+        point : an (*x*, *y*, *z*) triple
+            The point to test.
+
+        Returns
+        -------
+        contains : bool
+            Returns `True` if the polygon contains the given *point*.
+        """
+        for polygon in self.iter_polygons_flat():
+            if polygon.contains_point(point):
+                return True
+        return False
+
+    def intersects_poly(self, other):
+        r"""
+        Determines if this `SphericalPolygon` intersects another
+        `SphericalPolygon`.
+
+        This method is much faster than actually computing the
+        intersection region between two polygons.
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        intersects : bool
+            Returns `True` if this polygon intersects the *other*
+            polygon.
+        """
+        assert isinstance(other, SphericalPolygon)
+
+        for polya in self.iter_polygons_flat():
+            for polyb in other.iter_polygons_flat():
+                if polya.intersects_poly(polyb):
+                    return True
+        return False
+
+    def intersects_arc(self, a, b):
+        """
+        Determines if this `SphericalPolygon` intersects or contains
+        the given arc.
+        """
+        for subpolygon in self.iter_polygons_flat():
+            if subpolygon.intersects_arc(a, b):
+                return True
+        return False
+
+    def contains_arc(self, a, b):
+        """
+        Returns `True` if the polygon fully encloses the arc given by a
+        and b.
+        """
+        for subpolygon in self.iter_polygons_flat():
+            if subpolygon.contains_arc(a, b):
+                return True
+        return False
+
+    def area(self):
+        r"""
+        Returns the area of the polygon on the unit sphere in
+        steradians.
+
+        The area is computed using a generalization of Girard's Theorem.
+
+        if :math:`\theta` is the sum of the internal angles of the
+        polygon, and *n* is the number of vertices, the area is:
+
+        .. math::
+
+            S = \theta - (n - 2) \pi
+        """
+        area = 0.0
+        for subpoly in self.iter_polygons_flat():
+            area += subpoly.area()
+        return area
+
+    def union(self, other):
+        """
+        Return a new `SphericalPolygon` that is the union of *self*
+        and *other*.
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+
+        See also
+        --------
+        multi_union
+
+        Notes
+        -----
+        For implementation details, see the :mod:`~spherical_geometry.graph`
+        module.
+        """
+        from . import graph
+
+        if self.area() == 0.0:
+            return other.copy()
+        elif other.area() == 0.0:
+            return self.copy()
+
+        g = graph.Graph(
+            list(self.iter_polygons_flat()) +
+            list(other.iter_polygons_flat()))
+
+        return g.union()
+
+    @classmethod
+    def multi_union(cls, polygons):
+        """
+        Return a new `SphericalPolygon` that is the union of all of the
+        polygons in *polygons*.
+
+        Parameters
+        ----------
+        polygons : sequence of `SphericalPolygon`
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+
+        See also
+        --------
+        union
+        """
+        assert len(polygons)
+        for polygon in polygons:
+            assert isinstance(polygon, SphericalPolygon)
+
+        from . import graph
+
+        all_polygons = []
+        for polygon in polygons:
+            all_polygons.extend(polygon.iter_polygons_flat())
+
+        g = graph.Graph(all_polygons)
+        return g.union()
+
+    def intersection(self, other):
+        """
+        Return a new `SphericalPolygon` that is the intersection of
+        *self* and *other*.
+
+        If the intersection is empty, a `SphericalPolygon` with zero
+        subpolygons will be returned.
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+
+        Notes
+        -----
+        For implementation details, see the
+        :mod:`~spherical_geometry.graph` module.
+        """
+        from . import graph
+
+        if self.area() == 0.0:
+            return SphericalPolygon([])
+        elif other.area() == 0.0:
+            return SphericalPolygon([])
+
+        g = graph.Graph(
+            list(self.iter_polygons_flat()) +
+            list(other.iter_polygons_flat()))
+
+        return g.intersection()
 
     @classmethod
     def multi_intersection(cls, polygons, method='parallel'):
@@ -758,29 +1026,26 @@ class SphericalPolygon(object):
         for polygon in polygons:
             assert isinstance(polygon, SphericalPolygon)
 
-        # for i in range(len(polygons)):
-        #     polyA = polygons[i]
-        #     for j in range(i + 1, len(polygons)):
-        #         polyB = polygons[j]
-        #         if not polyA.intersects_poly(polyB):
-        #             return cls([], [0, 0, 0])
-
         from . import graph
 
+        all_polygons = []
+        for polygon in polygons:
+            if polygon.area() == 0.0:
+                return SphericalPolygon([])
+            all_polygons.extend(polygon.iter_polygons_flat())
+
         if method.lower() == 'parallel':
-            g = graph.Graph(polygons)
-            polygon = g.intersection()
-            inside = cls._find_new_inside(polygon, polygons)
-            return cls(polygon, inside)
+            g = graph.Graph(all_polygons)
+            return g.intersection()
         elif method.lower() == 'serial':
-            result = copy(polygons[0])
-            for polygon in polygons[1:]:
-                result = result.intersection(polygon)
+            results = copy(all_polygons[0])
+            for polygon in all_polygons[1:]:
+                results = results.intersection(polygon)
                 # If we have a null intersection already, we don't
                 # need to go any further.
-                if len(result._points) < 3:
-                    return result
-            return result
+                if len(results.polygons) == 0:
+                    return results
+            return results
         else:
             raise ValueError("method must be 'parallel' or 'serial'")
 
@@ -819,24 +1084,5 @@ class SphericalPolygon(object):
 
         **plot_args : Any plot arguments to pass to basemap
         """
-        if not len(self._points):
-            return
-        if not len(plot_args):
-            plot_args = {'color': 'blue'}
-        points = self._points
-
-        for A, B in zip(points[0:-1], points[1:]):
-            length = great_circle_arc.length(A, B, degrees=True)
-            if not np.isfinite(length):
-                length = 2
-            interpolated = great_circle_arc.interpolate(A, B, length * 4)
-            ra, dec = vector.vector_to_radec(
-                interpolated[:, 0], interpolated[:, 1], interpolated[:, 2],
-                degrees=True)
-            for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]):
-                m.drawgreatcircle(r0, d0, r1, d1, **plot_args)
-
-        ra, dec = vector.vector_to_radec(
-            *self._inside, degrees=True)
-        x, y = m(ra, dec)
-        m.scatter(x, y, 1, **plot_args)
+        for polygon in self._polygons:
+            polygon.draw(m, **plot_args)
diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py
index 567ec1d..5b2d47e 100644
--- a/stsci/sphere/test/test.py
+++ b/stsci/sphere/test/test.py
@@ -24,6 +24,13 @@ def test_normalize_vector():
     l = np.sqrt(np.sum(xyzn * xyzn, axis=-1))
     assert_almost_equal(l, 1.0)
 
+def test_normalize_unit_vector():
+    for i in range(3):
+        xyz = [0.0, 0.0, 0.0]
+        xyz[i] = 1.0
+        xyzn = vector.normalize_vector(xyz)
+        l = np.sqrt(np.sum(xyzn * xyzn, axis=-1))
+        assert_almost_equal(l, 1.0)
 
 def test_radec_to_vector():
     npx, npy, npz = vector.radec_to_vector(np.arange(-360, 360, 1), 90)
@@ -250,7 +257,7 @@ def test_great_circle_arc_angle():
     A = [1, 0, 0]
     B = [0, 1, 0]
     C = [0, 0, 1]
-    assert great_circle_arc.angle(A, B, C) == 90.0
+    assert great_circle_arc.angle(A, B, C) == 270.0
 
     # TODO: More angle tests
 
@@ -260,8 +267,7 @@ def test_cone():
     for i in range(50):
         ra = random.randrange(-180, 180)
         dec = random.randrange(20, 90)
-        cone = polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64)
-
+        polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64)
 
 def test_area():
     triangles = [
@@ -277,6 +283,13 @@ def test_area():
         poly = polygon.SphericalPolygon(points)
         calc_area = poly.area()
 
+def test_cone_area():
+    saved_area = None
+    for ra in  (0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330):
+        for dec in (0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330):
+            area = polygon.SphericalPolygon.from_cone(ra, dec, 30, steps=64).area()
+            if saved_area is None: saved_area = area 
+            assert_almost_equal(area, saved_area)
 
 def test_fast_area():
     a = np.array(  # Clockwise
@@ -304,6 +317,14 @@ def test_fast_area():
         [ 0.3536442 ,  0.63515101, -0.68667239],
         [ 0.35331737,  0.6351013 , -0.68688658]])
 
-    assert graph.Graph._fast_area(a) > 0
-    assert graph.Graph._fast_area(b) > 0
-    assert graph.Graph._fast_area(c) < 0
+    apoly = polygon._SingleSphericalPolygon(a)
+    bpoly = polygon._SingleSphericalPolygon(b)
+    cpoly = polygon._SingleSphericalPolygon(c)
+
+    aarea = apoly.area()
+    barea = bpoly.area()
+    carea = cpoly.area()
+    
+    assert aarea > 0 and aarea < np.pi * 2.0
+    assert barea > 0 and barea < np.pi * 2.0
+    assert carea > np.pi * 2.0 and carea < np.pi * 4.0
diff --git a/stsci/sphere/test/test_intersection.py b/stsci/sphere/test/test_intersection.py
index 9b01d8b..dc413a7 100644
--- a/stsci/sphere/test/test_intersection.py
+++ b/stsci/sphere/test/test_intersection.py
@@ -40,7 +40,7 @@ class intersection_test:
             num_permutations = math.factorial(len(polys))
             step_size = int(max(float(num_permutations) / 20.0, 1.0))
 
-            areas = [x.area() for x in polys]
+            areas = np.array([x.area() for x in polys])
 
             if GRAPH_MODE:
                 print("%d permutations" % num_permutations)
@@ -72,9 +72,11 @@ class intersection_test:
                         plt.savefig(filename)
                         fig.clear()
 
-                    assert np.all(intersection_area * 0.9 <= areas)
+                    assert np.all(areas >= intersection_area * 0.9)
 
-            lengths = np.array([len(x._points) for x in intersections])
+            lengths = np.array([
+                np.sum(len(x._points) for x in y.iter_polygons_flat())
+                for y in intersections])
             assert np.all(lengths == [lengths[0]])
             areas = np.array([x.area() for x in intersections])
             assert_array_almost_equal(areas, areas[0], decimal=1)
@@ -164,7 +166,7 @@ def test_intersection_empty():
 
     p2 = p.intersection(polygon.SphericalPolygon([]))
 
-    assert_array_almost_equal(p2._points, [])
+    assert len(p2.polygons) == 0
 
 
 def test_difficult_intersections():
@@ -236,10 +238,14 @@ def test_ordering():
     assert_array_almost_equal(areas[:-1], areas[1:])
 
     def roll_polygon(P, i):
-        points = P.points
-        points = np.roll(points[:-1], i, 0)
-        points = np.append(points, [points[0]], 0)
-        return polygon.SphericalPolygon(points, P.inside)
+        polygons = []
+        for p in P.polygons:
+            points = p.points
+            points = np.roll(points[:-1], i, 0)
+            points = np.append(points, [points[0]], 0)
+            p = polygon._SingleSphericalPolygon(points, p.inside)
+            polygons.append(p)
+        return polygon.SphericalPolygon(polygons)
 
     Aareas = []
     Bareas = []
diff --git a/stsci/sphere/test/test_union.py b/stsci/sphere/test/test_union.py
index 5559c87..b172e0f 100644
--- a/stsci/sphere/test/test_union.py
+++ b/stsci/sphere/test/test_union.py
@@ -56,9 +56,7 @@ class union_test:
                     permutation)
                 unions.append(union)
                 union_area = union.area()
-                print(union._points)
-                print(permutation[0]._points)
-
+                
                 if GRAPH_MODE:
                     fig = plt.figure()
                     m = Basemap(projection=self._proj,
@@ -74,10 +72,11 @@ class union_test:
                     plt.savefig(filename)
                     fig.clear()
 
-                print(union_area, areas)
                 assert np.all(union_area * 1.1 >= areas)
 
-            lengths = np.array([len(x._points) for x in unions])
+            lengths = np.array([
+                np.sum(len(x._points) for x in y.iter_polygons_flat())
+                for y in unions])
             assert np.all(lengths == [lengths[0]])
             areas = np.array([x.area() for x in unions])
             assert_array_almost_equal(areas, areas[0], 1)
@@ -188,7 +187,8 @@ def test_union_empty():
 
     p2 = p.union(polygon.SphericalPolygon([]))
 
-    assert_array_almost_equal(p2._points, p._points)
+    assert len(p2.polygons) == 1
+    assert_array_almost_equal(p2.polygons[0].points, p.polygons[0].points)
 
 
 def test_difficult_unions():
diff --git a/stsci/sphere/vector.py b/stsci/sphere/vector.py
index bb444e6..f13b9e7 100644
--- a/stsci/sphere/vector.py
+++ b/stsci/sphere/vector.py
@@ -44,6 +44,7 @@ import numpy as np
 
 try:
     from . import math_util
+    import math_util
     HAS_C_UFUNCS = True
 except ImportError:
     HAS_C_UFUNCS = False