mod for d2to1 install

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

Former-commit-id: 91667828b7e01b5aae55679093564473a976e8a9

18 files changed, 3336 insertions, 0 deletions

diff --git a/stsci/__init__.py b/stsci/__init__.py
new file mode 100644
index 0000000..a5133cc
--- /dev/null
+++ b/stsci/__init__.py
@@ -0,0 +1,8 @@
+# This is a special __init__.py required to place sample_package under the
+# stsci namespace package.  There should be no other code in this module.
+try:
+    from pkg_resources import declare_namespace
+    declare_namespace(__name__)
+except ImportError:
+    from pkgutil import extend_path
+    __path__ = extend_path(__path__, __name__)
diff --git a/stsci/sphere/__init__.py b/stsci/sphere/__init__.py
new file mode 100644
index 0000000..6a03273
--- /dev/null
+++ b/stsci/sphere/__init__.py
@@ -0,0 +1,5 @@
+import sys
+
+if sys.version_info[0] >= 3:
+    # Python 3 compatibility
+    __builtins__['xrange'] = range
diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py
new file mode 100644
index 0000000..d8ccf75
--- /dev/null
+++ b/stsci/sphere/graph.py
@@ -0,0 +1,909 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+#     1. Redistributions of source code must retain the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer.
+#
+#     2. Redistributions in binary form must reproduce the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer in the documentation and/or other materials
+#       provided with the distribution.
+#
+#     3. The name of AURA and its representatives may not be used to
+#       endorse or promote products derived from this software without
+#       specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR
+# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+# OF THE POSSIBILITY OF SUCH DAMAGE.
+
+"""
+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
+import itertools
+import weakref
+
+# THIRD-PARTY
+import numpy as np
+
+# LOCAL
+from . import great_circle_arc
+from . import vector
+
+# Set to True to enable some sanity checks
+DEBUG = True
+
+
+class Graph:
+    """
+    A graph of nodes connected by edges.  The graph is used to build
+    unions between polygons.
+
+    .. note::
+       This class is not meant to be used directly.  Instead, use
+       `~sphere.polygon.SphericalPolygon.union` and
+       `~sphere.polygon.SphericalPolygon.intersection`.
+    """
+
+    class Node:
+        """
+        A `~Graph.Node` represents a single point, connected by an arbitrary
+        number of `~Graph.Edge` objects to other `~Graph.Node` objects.
+        """
+        def __init__(self, point, source_polygons=[]):
+            """
+            Parameters
+            ----------
+            point : 3-sequence (*x*, *y*, *z*) coordinate
+
+            source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+                The polygon(s) this node came from.  Used for bookkeeping.
+            """
+            point = vector.normalize_vector(*point)
+
+            self._point = np.asanyarray(point)
+            self._source_polygons = set(source_polygons)
+            self._edges = weakref.WeakSet()
+
+        def __repr__(self):
+            return "Node(%s %d)" % (str(self._point), len(self._edges))
+
+        def follow(self, edge):
+            """
+            Follows from one edge to another across this node.
+
+            Parameters
+            ----------
+            edge : `~Graph.Edge` instance
+                The edge to follow away from.
+
+            Returns
+            -------
+            other : `~Graph.Edge` instance
+                The other edge.
+            """
+            edges = list(self._edges)
+            assert len(edges) == 2
+            try:
+                return edges[not edges.index(edge)]
+            except IndexError:
+                raise ValueError("Following from disconnected edge")
+
+        def equals(self, other, thres=2e-8):
+            """
+            Returns `True` if the location of this and the *other*
+            `~Graph.Node` are the same.
+
+            Parameters
+            ----------
+            other : `~Graph.Node` instance
+                The other node.
+
+            thres : float
+                If difference is smaller than this, points are equal.
+                The default value of 2e-8 radians is set based on
+                emphirical test cases. Relative threshold based on
+                the actual sizes of polygons is not implemented.
+            """
+            # return np.array_equal(self._point, other._point)
+            return great_circle_arc.length(self._point, other._point,
+                                           degrees=False) < thres
+
+
+    class Edge:
+        """
+        An `~Graph.Edge` represents a connection between exactly two
+        `~Graph.Node` objects.  This `~Graph.Edge` class has no direction.
+        """
+        def __init__(self, A, B, source_polygons=[]):
+            """
+            Parameters
+            ----------
+            A, B : `~Graph.Node` instances
+
+            source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+                The polygon this edge came from.  Used for bookkeeping.
+            """
+            self._nodes = [A, B]
+            for node in self._nodes:
+                node._edges.add(self)
+            self._source_polygons = set(source_polygons)
+
+        def __repr__(self):
+            nodes = self._nodes
+            return "Edge(%s -> %s)" % (nodes[0]._point, nodes[1]._point)
+
+        def follow(self, node):
+            """
+            Follow along the edge from the given *node* to the other
+            node.
+
+            Parameters
+            ----------
+            node : `~Graph.Node` instance
+
+            Returns
+            -------
+            other : `~Graph.Node` instance
+            """
+            nodes = self._nodes
+            try:
+                return nodes[not nodes.index(node)]
+            except IndexError:
+                raise ValueError("Following from disconnected node")
+
+        def equals(self, other):
+            """
+            Returns `True` if the other edge is between the same two nodes.
+
+            Parameters
+            ----------
+            other : `~Graph.Edge` instance
+
+            Returns
+            -------
+            equals : bool
+            """
+            if (self._nodes[0].equals(other._nodes[0]) and
+                self._nodes[1].equals(other._nodes[1])):
+                return True
+            if (self._nodes[1].equals(other._nodes[0]) and
+                self._nodes[0].equals(other._nodes[1])):
+                return True
+            return False
+
+
+    def __init__(self, polygons):
+        """
+        Parameters
+        ----------
+        polygons : sequence of `~sphere.polygon.SphericalPolygon` instances
+            Build a graph from this initial set of polygons.
+        """
+        self._nodes = set()
+        self._edges = set()
+        self._source_polygons = set()
+        self._start_node = None
+
+        self.add_polygons(polygons)
+
+    def add_polygons(self, polygons):
+        """
+        Add more polygons to the graph.
+
+        .. note::
+            Must be called before `union` or `intersection`.
+
+        Parameters
+        ----------
+        polygons : sequence of `~sphere.polygon.SphericalPolygon` instances
+            Set of polygons to add to the graph
+        """
+        for polygon in polygons:
+            self.add_polygon(polygon)
+
+    def add_polygon(self, polygon):
+        """
+        Add a single polygon to the graph.
+
+        .. note::
+            Must be called before `union` or `intersection`.
+
+        Parameters
+        ----------
+        polygon : `~sphere.polygon.SphericalPolygon` instance
+            Polygon to add to the graph
+        """
+        points = polygon._points
+
+        if len(points) < 3:
+            raise ValueError("Too few points in polygon")
+
+        self._source_polygons.add(polygon)
+
+        start_node = nodeA = self.add_node(points[0], [polygon])
+        if self._start_node is None:
+            self._start_node = start_node
+        for i in range(1, len(points) - 1):
+            nodeB = self.add_node(points[i], [polygon])
+            # Don't create self-pointing edges
+            if nodeB is not nodeA:
+                self.add_edge(nodeA, nodeB, [polygon])
+                nodeA = nodeB
+        # Close the polygon
+        self.add_edge(nodeA, start_node, [polygon])
+
+    def add_node(self, point, source_polygons=[]):
+        """
+        Add a node to the graph.  It will be disconnected until used
+        in a call to `add_edge`.
+
+        Parameters
+        ----------
+        point : 3-sequence (*x*, *y*, *z*) coordinate
+
+        source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional
+            The polygon this node came from.  Used for bookkeeping.
+
+        Returns
+        -------
+        node : `~Graph.Node` instance
+            The new node
+        """
+        new_node = self.Node(point, source_polygons)
+
+        # Don't add nodes that already exist.  Update the existing
+        # node's source_polygons list to include the new polygon.
+        for node in self._nodes:
+            if node.equals(new_node):
+                node._source_polygons.update(source_polygons)
+                return node
+
+        self._nodes.add(new_node)
+        return new_node
+
+    def remove_node(self, node):
+        """
+        Removes a node and all of the edges that touch it.
+
+        .. note::
+            It is assumed that *Node* is already a part of the graph.
+
+        Parameters
+        ----------
+        node : `~Graph.Node` instance
+        """
+        assert node in self._nodes
+
+        for edge in list(node._edges):
+            nodeB = edge.follow(node)
+            nodeB._edges.remove(edge)
+            self._edges.remove(edge)
+        self._nodes.remove(node)
+
+    def add_edge(self, A, B, source_polygons=[]):
+        """
+        Add an edge between two nodes.
+
+        .. note::
+            It is assumed both nodes already belong to the graph.
+
+        Parameters
+        ----------
+        A, B : `~Graph.Node` instances
+
+        source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional
+            The polygon(s) this edge came from.  Used for bookkeeping.
+
+        Returns
+        -------
+        edge : `~Graph.Edge` instance
+            The new edge
+        """
+        assert A in self._nodes
+        assert B in self._nodes
+
+        # Don't add any edges that already exist.  Update the edge's
+        # source polygons list to include the new polygon.  Care needs
+        # to be taken here to not create an Edge until we know we need
+        # one, otherwise the Edge will get hooked up to the nodes but
+        # be orphaned.
+        for edge in self._edges:
+            if ((A.equals(edge._nodes[0]) and
+                 B.equals(edge._nodes[1])) or
+                (B.equals(edge._nodes[0]) and
+                 A.equals(edge._nodes[1]))):
+                edge._source_polygons.update(source_polygons)
+                return edge
+
+        new_edge = self.Edge(A, B, source_polygons)
+        self._edges.add(new_edge)
+        return new_edge
+
+    def remove_edge(self, edge):
+        """
+        Remove an edge from the graph.  The nodes it points to remain intact.
+
+        .. note::
+            It is assumed that *edge* is already a part of the graph.
+
+        Parameters
+        ----------
+        edge : `~Graph.Edge` instance
+        """
+        assert edge in self._edges
+
+        A, B = edge._nodes
+        A._edges.remove(edge)
+        if len(A._edges) == 0:
+            self.remove_node(A)
+        if A is not B:
+            B._edges.remove(edge)
+            if len(B._edges) == 0:
+                self.remove_node(B)
+        self._edges.remove(edge)
+
+    def split_edge(self, edge, node):
+        """
+        Splits an `~Graph.Edge` *edge* at `~Graph.Node` *node*, removing
+        *edge* and replacing it with two new `~Graph.Edge` instances.
+        It is intended that *E* is along the original edge, but that is
+        not enforced.
+
+        Parameters
+        ----------
+        edge : `~Graph.Edge` instance
+            The edge to split
+
+        node : `~Graph.Node` instance
+            The node to insert
+
+        Returns
+        -------
+        edgeA, edgeB : `~Graph.Edge` instances
+            The two new edges on either side of *node*.
+        """
+        assert edge in self._edges
+        assert node in self._nodes
+
+        A, B = edge._nodes
+        edgeA = self.add_edge(A, node, edge._source_polygons)
+        edgeB = self.add_edge(node, B, edge._source_polygons)
+        if edge not in (edgeA, edgeB):
+            self.remove_edge(edge)
+        return [edgeA, edgeB]
+
+    def _sanity_check(self, title, node_is_2=False):
+        """
+        For debugging purposes: assert that edges and nodes are
+        connected to each other correctly and there are no orphaned
+        edges or nodes.
+        """
+        if not DEBUG:
+            return
+
+        try:
+            unique_edges = set()
+            for edge in self._edges:
+                for node in edge._nodes:
+                    assert edge in node._edges
+                    assert node in self._nodes
+                edge_repr = [tuple(x._point) for x in edge._nodes]
+                edge_repr.sort()
+                edge_repr = tuple(edge_repr)
+                # assert edge_repr not in unique_edges
+                unique_edges.add(edge_repr)
+
+            for node in self._nodes:
+                if node_is_2:
+                    assert len(node._edges) == 2
+                else:
+                    assert len(node._edges) >= 2
+                for edge in node._edges:
+                    assert node in edge._nodes
+                    assert edge in self._edges
+        except AssertionError as e:
+            import traceback
+            traceback.print_exc()
+            self._dump_graph(title=title)
+            raise
+
+    def _dump_graph(self, title=None, lon_0=0, lat_0=90,
+                    projection='vandg', func=lambda x: len(x._edges)):
+        from mpl_toolkits.basemap import Basemap
+        from matplotlib import pyplot as plt
+        fig = plt.figure()
+        m = Basemap()
+
+        counts = {}
+        for node in self._nodes:
+            count = func(node)
+            counts.setdefault(count, [])
+            counts[count].append(list(node._point))
+
+        minx = np.inf
+        miny = np.inf
+        maxx = -np.inf
+        maxy = -np.inf
+        for k, v in counts.items():
+            v = np.array(v)
+            ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2])
+            x, y = m(ra, dec)
+            m.plot(x, y, 'o', label=str(k))
+            for x0 in x:
+                minx = min(x0, minx)
+                maxx = max(x0, maxx)
+            for y0 in y:
+                miny = min(y0, miny)
+                maxy = max(y0, maxy)
+
+        for edge in list(self._edges):
+            A, B = [x._point for x in edge._nodes]
+            r0, d0 = vector.vector_to_radec(A[0], A[1], A[2])
+            r1, d1 = vector.vector_to_radec(B[0], B[1], B[2])
+            m.drawgreatcircle(r0, d0, r1, d1, color='blue')
+
+        plt.xlim(minx, maxx)
+        plt.ylim(miny, maxy)
+        if title:
+            plt.title("%s, %d v, %d e" % (
+                title, len(self._nodes), len(self._edges)))
+        plt.legend()
+        plt.show()
+
+    def union(self):
+        """
+        Once all of the polygons have been added to the graph,
+        join the polygons together.
+
+        Returns
+        -------
+        points : Nx3 array of (*x*, *y*, *z*) points
+            This is a list of points outlining the union of the
+            polygons that were given to the constructor.  If the
+            original polygons are disjunct or contain holes, cut lines
+            will be included in the output.
+        """
+        self._remove_cut_lines()
+        self._sanity_check("union - remove cut lines")
+        self._find_all_intersections()
+        self._sanity_check("union - find all intersections")
+        self._remove_interior_edges()
+        self._sanity_check("union - remove interior edges")
+        self._remove_3ary_edges()
+        self._sanity_check("union - remove 3ary edges")
+        self._remove_orphaned_nodes()
+        self._sanity_check("union - remove orphan nodes", True)
+        return self._trace()
+
+    def intersection(self):
+        """
+        Once all of the polygons have been added to the graph,
+        calculate the intersection.
+
+        Returns
+        -------
+        points : Nx3 array of (*x*, *y*, *z*) points
+            This is a list of points outlining the intersection of the
+            polygons that were given to the constructor.  If the
+            resulting polygons are disjunct or contain holes, cut lines
+            will be included in the output.
+        """
+        self._remove_cut_lines()
+        self._sanity_check("intersection - remove cut lines")
+        self._find_all_intersections()
+        self._sanity_check("intersection - find all intersections")
+        self._remove_exterior_edges()
+        self._sanity_check("intersection - remove exterior edges")
+        self._remove_3ary_edges(large_first=True)
+        self._sanity_check("intersection - remove 3ary edges")
+        self._remove_orphaned_nodes()
+        self._sanity_check("intersection - remove orphan nodes", True)
+        return self._trace()
+
+    def _remove_cut_lines(self):
+        """
+        Removes any cutlines that may already have existed in the
+        input polygons.  This is so any cutlines in the final result
+        will be optimized to be as short as possible and won't
+        intersect each other.
+
+        This works by finding coincident edges that are reverse to
+        each other, and then splicing around them.
+        """
+        # As this proceeds, edges are removed from the graph.  It
+        # iterates over a static list of all edges that exist at the
+        # start, so each time one is selected, we need to ensure it
+        # still exists as part of the graph.
+
+        # This transforms the following (where = is the cut line)
+        #
+        #     \                    /
+        #  A' +                    + B'
+        #     |                    |
+        #  A  +====================+ B
+        #
+        #  D  +====================+ C
+        #     |                    |
+        #  D' +                    + C'
+        #     /                    \
+        #
+        # to this:
+        #
+        #     \                    /
+        #  A' +                    + B'
+        #     |                    |
+        #  A  +                    + C
+        #     |                    |
+        #  D' +                    + C'
+        #     /                    \
+        #
+
+        edges = list(self._edges)
+        for i in xrange(len(edges)):
+            AB = edges[i]
+            if AB not in self._edges:
+                continue
+            A, B = AB._nodes
+            for j in xrange(i + 1, len(edges)):
+                CD = edges[j]
+                if CD not in self._edges:
+                    continue
+                C, D = CD._nodes
+                # To be a cutline, the candidate edges need to run in
+                # the opposite direction, hence A == D and B == C, not
+                # A == C and B == D.
+                if (A.equals(D) and B.equals(C)):
+                    # Create new edges A -> D' and C -> B'
+                    self.add_edge(
+                        A, D.follow(CD).follow(D),
+                        AB._source_polygons | CD._source_polygons)
+                    self.add_edge(
+                        C, B.follow(AB).follow(B),
+                        AB._source_polygons | CD._source_polygons)
+
+                    # Remove B and D which are identical to C and A
+                    # respectively.  We do not need to remove AB and
+                    # CD because this will remove it for us.
+                    self.remove_node(D)
+                    self.remove_node(B)
+
+                    break
+
+    def _find_all_intersections(self):
+        """
+        Find all the intersecting edges in the graph.  For each
+        intersecting pair, four new edges are created around the
+        intersection point.
+        """
+        def get_edge_points(edges):
+            return (np.array([x._nodes[0]._point for x in edges]),
+                    np.array([x._nodes[1]._point for x in edges]))
+
+        # For speed, we want to vectorize all of the intersection
+        # calculations.  Therefore, there is a list of edges, and two
+        # arrays containing the end points of those edges.  They all
+        # need to have things added and removed from them at the same
+        # time to keep them in sync, but of course the interface for
+        # doing so is different between Python lists and numpy arrays.
+
+        edges = list(self._edges)
+        starts, ends = get_edge_points(edges)
+
+        # First, split all edges by any nodes that intersect them
+        while len(edges) > 1:
+            AB = edges.pop(0)
+            A = starts[0]; starts = starts[1:]  # numpy equiv of "pop(0)"
+            B = ends[0];   ends = ends[1:]      # numpy equiv of "pop(0)"
+
+            distance = great_circle_arc.length(A, B)
+            for node in self._nodes:
+                if node not in AB._nodes:
+                    distanceA = great_circle_arc.length(node._point, A)
+                    distanceB = great_circle_arc.length(node._point, B)
+                    if np.abs((distanceA + distanceB) - distance) < 1e-8:
+                        newA, newB = self.split_edge(AB, node)
+
+                        new_edges = [
+                            edge for edge in (newA, newB)
+                            if edge not in edges]
+
+                        edges = new_edges + edges
+                        new_starts, new_ends = get_edge_points(new_edges)
+                        starts = np.vstack(
+                            (new_starts, starts))
+                        ends = np.vstack(
+                            (new_ends, ends))
+                        break
+
+        edges = list(self._edges)
+        starts, ends = get_edge_points(edges)
+
+        # Next, calculate edge-to-edge intersections and break
+        # edges on the intersection point.
+        while len(edges) > 1:
+            AB = edges.pop(0)
+            A = starts[0]; starts = starts[1:]  # numpy equiv of "pop(0)"
+            B = ends[0];   ends = ends[1:]      # numpy equiv of "pop(0)"
+
+            # Calculate the intersection points between AB and all
+            # other remaining edges
+            with np.errstate(invalid='ignore'):
+                intersections = great_circle_arc.intersection(
+                    A, B, starts, ends)
+            # intersects is `True` everywhere intersections has an
+            # actual intersection
+            intersects = np.isfinite(intersections[..., 0])
+
+            intersection_indices = np.nonzero(intersects)[0]
+
+            # Iterate through the candidate intersections, if any --
+            # we want to eliminate intersections that only intersect
+            # at the end points
+            for j in intersection_indices:
+                C = starts[j]
+                D = ends[j]
+                CD = edges[j]
+                E = intersections[j]
+
+                # This is a bona-fide intersection, and E is the
+                # point at which the two lines intersect.  Make a
+                # new node for it -- this must belong to the all
+                # of the source polygons of both of the edges that
+                # crossed.
+
+                #                A
+                #                |
+                #             C--E--D
+                #                |
+                #                B
+
+                E = self.add_node(
+                    E, AB._source_polygons | CD._source_polygons)
+                newA, newB = self.split_edge(AB, E)
+                newC, newD = self.split_edge(CD, E)
+
+                new_edges = [
+                    edge for edge in (newA, newB, newC, newD)
+                    if edge not in edges]
+
+                # 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
+                new_starts, new_ends = get_edge_points(new_edges)
+                starts = np.vstack(
+                    (new_starts, starts[:j], starts[j+1:]))
+                ends = np.vstack(
+                    (new_ends, ends[:j], ends[j+1:]))
+                break
+
+    def _remove_interior_edges(self):
+        """
+        Removes any nodes that are contained inside other polygons.
+        What's left is the (possibly disjunct) outline.
+        """
+        polygons = self._source_polygons
+
+        for edge in self._edges:
+            edge._count = 0
+            for polygon in polygons:
+                if (not polygon in edge._source_polygons and
+                    polygon.intersects_arc(
+                        edge._nodes[0]._point, edge._nodes[1]._point)):
+                    edge._count += 1
+
+        for edge in list(self._edges):
+            if edge._count >= 1:
+                self.remove_edge(edge)
+
+    def _remove_exterior_edges(self):
+        """
+        Removes any edges that are not contained in all of the source
+        polygons.  What's left is the (possibly disjunct) outline.
+        """
+        polygons = self._source_polygons
+
+        for edge in self._edges:
+            edge._count = 0
+            for polygon in polygons:
+                if (polygon in edge._source_polygons or
+                    polygon.intersects_arc(
+                        edge._nodes[0]._point, edge._nodes[1]._point)):
+                    edge._count += 1
+
+        for edge in list(self._edges):
+            if edge._count < len(polygons):
+                self.remove_edge(edge)
+
+    def _remove_3ary_edges(self, large_first=False):
+        """
+        Remove edges between pairs of nodes that have 3 or more edges.
+        This removes triangles that can't be traced.
+        """
+        if large_first:
+            max_ary = 0
+            for node in self._nodes:
+                max_ary = max(len(node._edges), max_ary)
+            order = range(max_ary + 1, 2, -1)
+        else:
+            order = [3]
+
+        for i in order:
+            removals = []
+            for edge in list(self._edges):
+                if (len(edge._nodes[0]._edges) >= i and
+                    len(edge._nodes[1]._edges) >= i):
+                    removals.append(edge)
+
+            for edge in removals:
+                if edge in self._edges:
+                    self.remove_edge(edge)
+
+    def _remove_orphaned_nodes(self):
+        """
+        Remove nodes with fewer than 2 edges.
+        """
+        while True:
+            removes = []
+            for node in list(self._nodes):
+                if len(node._edges) < 2:
+                    removes.append(node)
+            if len(removes):
+                for node in removes:
+                    self.remove_node(node)
+            else:
+                break
+
+    def _trace(self):
+        """
+        Given a graph that has had cutlines removed and all
+        intersections found, traces it to find a resulting single
+        polygon.
+        """
+        polygons = []
+        edges = set(self._edges)  # copy
+        seen_nodes = set()
+        while len(edges):
+            polygon = []
+            # Carefully pick out an "original" edge first.  Synthetic
+            # edges may not be pointing in the right direction to
+            # properly calculate the area.
+            for edge in edges:
+                if len(edge._source_polygons) == 1:
+                    break
+            edges.remove(edge)
+            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)
+                else:
+                    polygon.append(node._point)
+                edge = node.follow(edge)
+                edges.discard(edge)
+                node = edge.follow(node)
+                if node is start_node:
+                    polygon.append(node._point)
+                    break
+
+            polygons.append(np.asarray(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")
diff --git a/stsci/sphere/great_circle_arc.py b/stsci/sphere/great_circle_arc.py
new file mode 100644
index 0000000..d5e78f2
--- /dev/null
+++ b/stsci/sphere/great_circle_arc.py
@@ -0,0 +1,368 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+#     1. Redistributions of source code must retain the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer.
+#
+#     2. Redistributions in binary form must reproduce the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer in the documentation and/or other materials
+#       provided with the distribution.
+#
+#     3. The name of AURA and its representatives may not be used to
+#       endorse or promote products derived from this software without
+#       specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR
+# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+# OF THE POSSIBILITY OF SUCH DAMAGE.
+
+"""
+The `sphere.great_circle_arc` module contains functions for computing
+the length, intersection, angle and midpoint of great circle arcs.
+
+Great circles are circles on the unit sphere whose center is
+coincident with the center of the sphere.  Great circle arcs are the
+section of those circles between two points on the unit sphere.
+"""
+
+from __future__ import with_statement, division, absolute_import, unicode_literals
+
+# THIRD-PARTY
+import numpy as np
+
+
+try:
+    from . import math_util
+    HAS_C_UFUNCS = True
+except ImportError:
+    HAS_C_UFUNCS = False
+
+if HAS_C_UFUNCS:
+    inner1d = math_util.inner1d
+else:
+    from numpy.core.umath_tests import inner1d
+
+
+__all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint',
+           'interpolate']
+
+
+def _fast_cross(a, b):
+    """
+    This is a reimplementation of `numpy.cross` that only does 3D x
+    3D, and is therefore faster since it doesn't need any
+    conditionals.
+    """
+    if HAS_C_UFUNCS:
+        return math_util.cross(a, b)
+    
+    cp = np.empty(np.broadcast(a, b).shape)
+    aT = a.T
+    bT = b.T
+    cpT = cp.T
+
+    cpT[0] = aT[1]*bT[2] - aT[2]*bT[1]
+    cpT[1] = aT[2]*bT[0] - aT[0]*bT[2]
+    cpT[2] = aT[0]*bT[1] - aT[1]*bT[0]
+
+    return cp
+
+def _cross_and_normalize(A, B):
+    if HAS_C_UFUNCS:
+        return math_util.cross_and_norm(A, B)
+
+    T = _fast_cross(A, B)
+    # Normalization
+    l = np.sqrt(np.sum(T ** 2, axis=-1))
+    l = np.expand_dims(l, 2)
+    # Might get some divide-by-zeros, but we don't care
+    with np.errstate(invalid='ignore'):
+        TN = T / l
+    return TN
+
+
+def intersection(A, B, C, D):
+    r"""
+    Returns the point of intersection between two great circle arcs.
+    The arcs are defined between the points *AB* and *CD*.  Either *A*
+    and *B* or *C* and *D* may be arrays of points, but not both.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        Endpoints of the first great circle arc.
+
+    C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        Endpoints of the second great circle arc.
+
+    Returns
+    -------
+    T : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        If the given arcs intersect, the intersection is returned.  If
+        the arcs do not intersect, the triple is set to all NaNs.
+
+    Notes
+    -----
+    The basic intersection is computed using linear algebra as follows
+    [1]_:
+
+    .. math::
+
+        T = \lVert(A × B) × (C × D)\rVert
+
+    To determine the correct sign (i.e. hemisphere) of the
+    intersection, the following four values are computed:
+
+    .. math::
+
+        s_1 = ((A × B) × A) · T
+
+        s_2 = (B × (A × B)) · T
+
+        s_3 = ((C × D) × C) · T
+
+        s_4 = (D × (C × D)) · T
+
+    For :math:`s_n`, if all positive :math:`T` is returned as-is.  If
+    all negative, :math:`T` is multiplied by :math:`-1`.  Otherwise
+    the intersection does not exist and is undefined.
+
+    References
+    ----------
+
+    .. [1] Method explained in an `e-mail
+        <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_
+        by Roger Stafford.
+
+    http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271
+    """
+    if HAS_C_UFUNCS:
+        return math_util.intersection(A, B, C, D)
+    
+    A = np.asanyarray(A)
+    B = np.asanyarray(B)
+    C = np.asanyarray(C)
+    D = np.asanyarray(D)
+
+    A, B = np.broadcast_arrays(A, B)
+    C, D = np.broadcast_arrays(C, D)
+
+    ABX = _fast_cross(A, B)
+    CDX = _fast_cross(C, D)
+    T = _cross_and_normalize(ABX, CDX)
+    T_ndim = len(T.shape)
+
+    if T_ndim > 1:
+        s = np.zeros(T.shape[0])
+    else:
+        s = np.zeros(1)
+    s += np.sign(inner1d(_fast_cross(ABX, A), T))
+    s += np.sign(inner1d(_fast_cross(B, ABX), T))
+    s += np.sign(inner1d(_fast_cross(CDX, C), T))
+    s += np.sign(inner1d(_fast_cross(D, CDX), T))
+    if T_ndim > 1:
+        s = np.expand_dims(s, 2)
+
+    cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan))
+
+    # If they share a common point, it's not an intersection.  This
+    # gets around some rounding-error/numerical problems with the
+    # above.
+    equals = (np.all(A == C, axis = -1) |
+              np.all(A == D, axis = -1) |
+              np.all(B == C, axis = -1) |
+              np.all(B == D, axis = -1))
+
+    equals = np.expand_dims(equals, 2)
+
+    return np.where(equals, np.nan, cross)
+
+
+def length(A, B, degrees=True):
+    r"""
+    Returns the angular distance between two points (in vector space)
+    on the unit sphere.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+       The endpoints of the great circle arc, in vector space.
+
+    degrees : bool, optional
+        If `True` (default) the result is returned in decimal degrees,
+        otherwise radians.
+
+    Returns
+    -------
+    length : scalar or array of scalars
+        The angular length of the great circle arc.
+
+    Notes
+    -----
+    The length is computed using the following:
+
+    .. math::
+
+       \Delta = \arccos(A \dot B)
+    """
+    A = np.asanyarray(A)
+    B = np.asanyarray(B)
+
+    dot = inner1d(A, B)
+    dot = np.clip(dot, -1.0, 1.0)
+    with np.errstate(invalid='ignore'):
+        result = np.arccos(dot)
+
+    if degrees:
+        return np.rad2deg(result)
+    else:
+        return result
+
+
+def intersects(A, B, C, D):
+    """
+    Returns `True` if the great circle arcs between *AB* and *CD*
+    intersect.  Either *A* and *B* or *C* and *D* may be arrays of
+    points, but not both.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        Endpoints of the first great circle arc.
+
+    C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        Endpoints of the second great circle arc.
+
+    Returns
+    -------
+    intersects : bool or array of bool
+        If the given arcs intersect, the intersection is returned as
+        `True`.
+    """
+    with np.errstate(invalid='ignore'):
+        intersections = intersection(A, B, C, D)
+    return np.isfinite(intersections[..., 0])
+
+
+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
+        Points on sphere.
+
+    degrees : bool, optional
+        If `True` (default) the result is returned in decimal degrees,
+        otherwise radians.
+
+    Returns
+    -------
+    angle : float or array of floats
+        The angle at *B* between *AB* and *BC*.
+
+    References
+    ----------
+
+    .. [1] Miller, Robert D.  Computing the area of a spherical
+       polygon.  Graphics Gems IV.  1994.  Academic Press.
+    """
+    A = np.asanyarray(A)
+    B = np.asanyarray(B)
+    C = np.asanyarray(C)
+
+    A, B, C = np.broadcast_arrays(A, B, C)
+
+    ABX = _fast_cross(A, B)
+    ABX = _cross_and_normalize(B, ABX)
+    BCX = _fast_cross(C, B)
+    BCX = _cross_and_normalize(B, BCX)
+    with np.errstate(invalid='ignore'):
+        angle = np.arccos(inner1d(ABX, BCX))
+
+    if degrees:
+        angle = np.rad2deg(angle)
+    return angle
+
+
+def midpoint(A, B):
+    """
+    Returns the midpoint on the great circle arc between *A* and *B*.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        The endpoints of the great circle arc.  It is assumed that
+        these points are already normalized.
+
+    Returns
+    -------
+    midpoint : (*x*, *y*, *z*) triple or Nx3 arrays of triples
+        The midpoint between *A* and *B*, normalized on the unit
+        sphere.
+    """
+    P = (A + B) / 2.0
+    # Now normalize...
+    l = np.sqrt(np.sum(P * P, axis=-1))
+    l = np.expand_dims(l, 2)
+    return P / l
+
+
+def interpolate(A, B, steps=50):
+    r"""
+    Interpolate along the great circle arc.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        The endpoints of the great circle arc.  It is assumed thats
+        these points are already normalized.
+
+    steps : int
+        The number of interpolation steps
+
+    Returns
+    -------
+    array : (*x*, *y*, *z*) triples
+        The points interpolated along the great circle arc
+
+    Notes
+    -----
+
+    This uses Slerp interpolation where *Ω* is the angle subtended by
+    the arc, and *t* is the parameter 0 <= *t* <= 1.
+
+    .. math::
+
+        \frac{\sin((1 - t)\Omega)}{\sin \Omega}A + \frac{\sin(t \Omega)}{\sin \Omega}B
+    """
+    steps = max(steps, 2)
+    t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1))
+
+    omega = length(A, B, degrees=False)
+    if omega == 0.0:
+        offsets = t
+    else:
+        sin_omega = np.sin(omega)
+        offsets = np.sin(t * omega) / sin_omega
+
+    return offsets[::-1] * A + offsets * B
diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py
new file mode 100644
index 0000000..e9ad379
--- /dev/null
+++ b/stsci/sphere/polygon.py
@@ -0,0 +1,831 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+#     1. Redistributions of source code must retain the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer.
+#
+#     2. Redistributions in binary form must reproduce the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer in the documentation and/or other materials
+#       provided with the distribution.
+#
+#     3. The name of AURA and its representatives may not be used to
+#       endorse or promote products derived from this software without
+#       specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR
+# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+# OF THE POSSIBILITY OF SUCH DAMAGE.
+
+"""
+The `polygon` module defines the `SphericalPolygon` class for managing
+polygons on the unit sphere.
+"""
+from __future__ import division, print_function, unicode_literals, absolute_import
+
+# STDLIB
+from copy import copy, deepcopy
+
+# THIRD-PARTY
+import numpy as np
+
+# LOCAL
+from . import great_circle_arc
+from . import vector
+
+__all__ = ['SphericalPolygon']
+
+
+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.
+    """
+
+    def __init__(self, points, inside=None):
+        r"""
+        Parameters
+        ----------
+        points : An Nx3 array of (*x*, *y*, *z*) triples in vector 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
+            This point must be inside the polygon.  If not provided, the
+            mean of the points will be used.
+        """
+        if len(points) == 0:
+            # handle special case of initializing with an empty list of
+            # vertices (ticket #1079).
+            self._inside = np.zeros(3)
+            self._points = np.asanyarray(points)
+            return
+        elif len(points) < 3:
+            raise ValueError("Polygon made of too few points")
+        else:
+            assert np.array_equal(points[0], points[-1]), 'Polygon is not closed'
+
+        self._points = np.asanyarray(points)
+
+        if inside is None:
+            self._inside = np.mean(points[:-1], axis=0)
+        else:
+            self._inside = np.asanyarray(inside)
+
+        # TODO: Detect self-intersection and fix
+
+    def __copy__(self):
+        return deepcopy(self)
+
+    def __repr__(self):
+        return '%s(%r, %r)' % (self.__class__.__name__,
+                               self.points, self.inside)
+
+    @property
+    def points(self):
+        """
+        The points defining the polygon.  It is an Nx3 array of
+        (*x*, *y*, *z*) vectors.  The polygon will be explicitly
+        closed, i.e., the first and last points are the same.
+        """
+        return self._points
+
+    @property
+    def inside(self):
+        """
+        Get the inside point of the polygon.
+        """
+        return self._inside
+
+    def to_radec(self):
+        """
+        Convert `SphericalPolygon` footprint to RA and DEC.
+
+        Returns
+        -------
+        ra, dec : list of float
+            List of *ra* and *dec* in degrees corresponding
+            to `points`.
+        """
+        return vector.vector_to_radec(self.points[:,0], self.points[:,1],
+                                      self.points[:,2], degrees=True)
+
+    @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
+        """
+        # Convert to Cartesian
+        x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees)
+
+        if center is None:
+            xc = x.mean()
+            yc = y.mean()
+            zc = z.mean()
+            center = vector.normalize_vector(xc, yc, zc)
+        else:
+            center = vector.radec_to_vector(*center, degrees=degrees)
+
+        return cls(np.dstack((x, y, z))[0], center)
+
+    @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
+        """
+        u, v, w = vector.radec_to_vector(ra, dec, degrees=degrees)
+        if degrees:
+            radius = np.deg2rad(radius)
+
+        # 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])
+        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.])
+
+        # Rotate by radius around the perpendicular vector to get the
+        # "pen"
+        x, y, z = vector.rotate_around(
+            u, v, w, perp[0], perp[1], perp[2], radius, degrees=False)
+
+        # Then rotate the pen around the center point all 360 degrees
+        C = np.linspace(0, np.pi * 2.0, steps)
+        # Ensure that the first and last elements are exactly the
+        # same.  2π should equal 0, but with rounding error that isn't
+        # always the case.
+        C[-1] = 0
+        C = C[::-1]
+        X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False)
+
+        return cls(np.dstack((X, Y, Z))[0], (u, v, w))
+
+    @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 `pywcs` and `pyfits` installed.
+
+        Parameters
+        ----------
+        fitspath : path to a FITS file, `pyfits.Header`, or `pywcs.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
+        """
+        import pywcs
+        import pyfits
+
+        if isinstance(fitspath, pyfits.Header):
+            header = fitspath
+            wcs = pywcs.WCS(header)
+        elif isinstance(fitspath, pywcs.WCS):
+            wcs = fitspath
+        else:
+            wcs = pywcs.WCS(fitspath)
+        if crval is not None:
+            wcs.wcs.crval = crval
+        xa, ya = [wcs.naxis1, wcs.naxis2]
+
+        length = steps * 4 + 1
+        X = np.empty(length)
+        Y = np.empty(length)
+
+        # Now define each of the 4 edges of the quadrilateral
+        X[0      :steps  ] = np.linspace(1, xa, steps, False)
+        Y[0      :steps  ] = 1
+        X[steps  :steps*2] = xa
+        Y[steps  :steps*2] = np.linspace(1, ya, steps, False)
+        X[steps*2:steps*3] = np.linspace(xa, 1, steps, False)
+        Y[steps*2:steps*3] = ya
+        X[steps*3:steps*4] = 1
+        Y[steps*3:steps*4] = np.linspace(ya, 1, steps, False)
+        X[-1]              = 1
+        Y[-1]              = 1
+
+        # Use wcslib to convert to (ra, dec)
+        ra, dec = wcs.all_pix2sky(X, Y, 1)
+
+        # Convert to Cartesian
+        x, y, z = vector.radec_to_vector(ra, dec)
+
+        # Calculate an inside point
+        ra, dec = wcs.all_pix2sky(xa / 2.0, ya / 2.0, 1)
+        xc, yc, zc = vector.radec_to_vector(ra, dec)
+
+        return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0]))
+
+    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
+
+            if np.sqrt(x_sum) > is_same_limit:
+                is_eq = False
+                break
+
+        return is_eq
+
+    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*.
+        """
+        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
+
+    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.
+
+        Notes
+        -----
+
+        The algorithm proceeds as follows:
+
+            1. Determine if any single point of one polygon is contained
+               within the other.
+
+            2. Deal with the case where only the edges overlap as in::
+
+               :       o---------o
+               :  o----+---------+----o
+               :  |    |         |    |
+               :  o----+---------+----o
+               :       o---------o
+
+               In this case, an edge from one polygon must cross an
+               edge from the other polygon.
+        """
+        assert isinstance(other, SphericalPolygon)
+
+        # The easy case is in which a point of one polygon is
+        # contained in the other polygon.
+        for point in other._points:
+            if self.contains_point(point):
+                return True
+        for point in self._points:
+            if other.contains_point(point):
+                return True
+
+        # The hard case is when only the edges overlap, as in:
+        #
+        #         o---------o
+        #    o----+---------+----o
+        #    |    |         |    |
+        #    o----+---------+----o
+        #         o---------o
+        #
+        for i in range(len(self._points) - 1):
+            A = self._points[i]
+            B = self._points[i+1]
+            if np.any(great_circle_arc.intersects(
+                A, B, other._points[:-1], other._points[1:])):
+                return True
+        return False
+
+    def intersects_arc(self, a, b):
+        """
+        Determines if this `SphericalPolygon` intersects or contains
+        the given arc.
+        """
+        return self.contains_point(great_circle_arc.midpoint(a, b))
+
+        if self.contains_point(a) and self.contains_point(b):
+            return True
+
+        P = self._points
+        intersects = great_circle_arc.intersects(P[:-1], P[1:], a, b)
+        if np.any(intersects):
+            return True
+        return False
+
+    def area(self):
+        r"""
+        Returns the area of the polygon on the unit sphere.
+
+        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 by transforming the polygon to two
+        dimensions using the `Lambert azimuthal equal-area projection
+        <http://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_
+
+        .. 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|
+        """
+        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[:, 0], points[:, 1], points[:, 2], inplace=True)
+
+        X = []
+        Y = []
+        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)
+            x, y, z = vector.normalize_vector(
+                interp[:, 0], interp[:, 1], interp[:, 2], inplace=True)
+            x, y = vector.equal_area_proj(x, y, z)
+            X.extend(x)
+            Y.extend(y)
+
+        X = np.array(X)
+        Y = np.array(Y)
+
+        return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi)
+
+    def union(self, other):
+        """
+        Return a new `SphericalPolygon` that is the union of *self*
+        and *other*.
+
+        If the polygons are disjoint, they result will be connected
+        using cut lines.  For example::
+
+            : o---------o
+            : |         |
+            : o---------o=====o----------o
+            :                 |          |
+            :                 o----------o
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+
+        See also
+        --------
+        multi_union
+
+        Notes
+        -----
+        For implementation details, see the :mod:`~sphere.graph`
+        module.
+        """
+        from . import graph
+        g = graph.Graph([self, other])
+
+        polygon = g.union()
+
+        return self.__class__(polygon, self._inside)
+
+    @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):
+        """
+        Finds an acceptable inside point inside of *points* that is
+        also inside of *polygons*.  Used by the intersection
+        algorithm, and is really only useful in that context because
+        it requires existing polygons with known inside points.
+        """
+        if len(points) < 4:
+            return np.array([0, 0, 0])
+
+        # Special case for a triangle
+        if len(points) == 4:
+            return np.sum(points[:3]) / 3.0
+
+        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
+                        break
+                if in_all:
+                    return C
+
+        raise RuntimeError("Suitable inside point could not be found")
+
+    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
+        points will be returned.
+
+        If the result is disjoint, the pieces will be connected using
+        cut lines.  For example::
+
+            : o---------o
+            : |         |
+            : o---------o=====o----------o
+            :                 |          |
+            :                 o----------o
+
+        Parameters
+        ----------
+        other : `SphericalPolygon`
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+
+        Notes
+        -----
+        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])
+
+        g = graph.Graph([self, other])
+
+        polygon = g.intersection()
+
+        inside = self._find_new_inside(polygon, [self, other])
+
+        return self.__class__(polygon, inside)
+
+    @classmethod
+    def multi_intersection(cls, polygons, method='parallel'):
+        """
+        Return a new `SphericalPolygon` that is the intersection of
+        all of the polygons in *polygons*.
+
+        Parameters
+        ----------
+        polygons : sequence of `SphericalPolygon`
+
+        method : 'parallel' or 'serial', optional
+            Specifies the method that is used to perform the
+            intersections:
+
+               - 'parallel' (default): A graph is built using all of
+                 the polygons, and the intersection operation is computed on
+                 the entire thing globally.
+
+               - 'serial': The polygon is built in steps by adding one
+                 polygon at a time and computing the intersection at
+                 each step.
+
+            This option is provided because one may be faster than the
+            other depending on context, but it primarily exposed for
+            testing reasons.  Both modes should theoretically provide
+            equivalent results.
+
+        Returns
+        -------
+        polygon : `SphericalPolygon` object
+        """
+        assert len(polygons)
+        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
+
+        if method.lower() == 'parallel':
+            g = graph.Graph(polygons)
+            polygon = g.intersection()
+            inside = cls._find_new_inside(polygon, polygons)
+            return cls(polygon, inside)
+        elif method.lower() == 'serial':
+            result = copy(polygons[0])
+            for polygon in polygons[1:]:
+                result = result.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
+        else:
+            raise ValueError("method must be 'parallel' or 'serial'")
+
+    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
+
+        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)
diff --git a/stsci/sphere/test/__init__.py b/stsci/sphere/test/__init__.py
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/stsci/sphere/test/__init__.py
@@ -0,0 +1 @@
+
diff --git a/stsci/sphere/test/benchmarks.py b/stsci/sphere/test/benchmarks.py
new file mode 100644
index 0000000..91ef8db
--- /dev/null
+++ b/stsci/sphere/test/benchmarks.py
@@ -0,0 +1,39 @@
+import os
+import sys
+import time
+
+import numpy as np
+from sphere import *
+from test_util import *
+from test_shared import resolve_imagename
+
+def point_in_poly_lots():
+    image_name = resolve_imagename(ROOT_DIR,'1904-66_TAN.fits')
+    
+    poly1 = SphericalPolygon.from_wcs(image_name, 64, crval=[0, 87])
+    poly2 = SphericalPolygon.from_wcs(image_name, 64, crval=[20, 89])
+    poly3 = SphericalPolygon.from_wcs(image_name, 64, crval=[180, 89])
+
+    points = get_point_set(density=25)
+
+    count = 0
+    for point in points:
+        if poly1.contains_point(point) or poly2.contains_point(point) or \
+               poly3.contains_point(point):
+            count += 1
+
+    assert count == 5
+    assert poly1.intersects_poly(poly2)
+    assert not poly1.intersects_poly(poly3)
+    assert not poly2.intersects_poly(poly3)
+
+if __name__ == '__main__':
+    for benchmark in [point_in_poly_lots]:
+        t = time.time()
+        sys.stdout.write(benchmark.__name__)
+        sys.stdout.write('...')
+        sys.stdout.flush()
+        benchmark()
+        sys.stdout.write(' %.03fs\n' % (time.time() - t))
+        sys.stdout.flush()
+
diff --git a/stsci/sphere/test/data/1904-66_TAN.fits.gz b/stsci/sphere/test/data/1904-66_TAN.fits.gz
new file mode 100644
index 0000000..34ff168
--- /dev/null
+++ b/stsci/sphere/test/data/1904-66_TAN.fits.gz
diff --git a/stsci/sphere/test/data/2chipA.fits.gz b/stsci/sphere/test/data/2chipA.fits.gz
new file mode 100644
index 0000000..27a68ab
--- /dev/null
+++ b/stsci/sphere/test/data/2chipA.fits.gz
diff --git a/stsci/sphere/test/data/2chipB.fits.gz b/stsci/sphere/test/data/2chipB.fits.gz
new file mode 100644
index 0000000..209d5f9
--- /dev/null
+++ b/stsci/sphere/test/data/2chipB.fits.gz
diff --git a/stsci/sphere/test/data/2chipC.fits.gz b/stsci/sphere/test/data/2chipC.fits.gz
new file mode 100644
index 0000000..10af627
--- /dev/null
+++ b/stsci/sphere/test/data/2chipC.fits.gz
diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py
new file mode 100644
index 0000000..b551090
--- /dev/null
+++ b/stsci/sphere/test/test.py
@@ -0,0 +1,277 @@
+from __future__ import absolute_import
+
+import os
+import random
+
+import numpy as np
+from numpy.testing import assert_almost_equal, assert_array_less
+
+from .. import graph
+from .. import great_circle_arc
+from .. import polygon
+from .. import vector
+
+from .test_util import *
+from .test_shared import resolve_imagename
+
+graph.DEBUG = True
+
+
+def test_normalize_vector():
+    x, y, z = np.ogrid[-100:100:11,-100:100:11,-100:100:11]
+    xn, yn, zn = vector.normalize_vector(x.flatten(), y.flatten(), z.flatten())
+    l = np.sqrt(xn ** 2 + yn ** 2 + zn ** 2)
+    assert_almost_equal(l, 1.0)
+
+
+def test_radec_to_vector():
+    npx, npy, npz = vector.radec_to_vector(np.arange(-360, 360, 1), 90)
+    assert_almost_equal(npx, 0.0)
+    assert_almost_equal(npy, 0.0)
+    assert_almost_equal(npz, 1.0)
+
+    spx, spy, spz = vector.radec_to_vector(np.arange(-360, 360, 1), -90)
+    assert_almost_equal(spx, 0.0)
+    assert_almost_equal(spy, 0.0)
+    assert_almost_equal(spz, -1.0)
+
+    eqx, eqy, eqz = vector.radec_to_vector(np.arange(-360, 360, 1), 0)
+    assert_almost_equal(eqz, 0.0)
+
+
+def test_vector_to_radec():
+    ra, dec = vector.vector_to_radec(0, 0, 1)
+    assert_almost_equal(dec, 90)
+
+    ra, dec = vector.vector_to_radec(0, 0, -1)
+    assert_almost_equal(dec, -90)
+
+    ra, dec = vector.vector_to_radec(1, 1, 0)
+    assert_almost_equal(ra, 45.0)
+    assert_almost_equal(dec, 0.0)
+
+
+def test_intersects_poly_simple():
+    ra1 = [-10, 10, 10, -10, -10]
+    dec1 = [30, 30, 0, 0, 30]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = [-5, 15, 15, -5, -5]
+    dec2 = [20, 20, -10, -10, 20]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+    # Make sure it isn't order-dependent
+    ra1 = ra1[::-1]
+    dec1 = dec1[::-1]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = ra2[::-1]
+    dec2 = dec2[::-1]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+
+def test_intersects_poly_fully_contained():
+    ra1 = [-10, 10, 10, -10, -10]
+    dec1 = [30, 30, 0, 0, 30]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = [-5, 5, 5, -5, -5]
+    dec2 = [20, 20, 10, 10, 20]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+    # Make sure it isn't order-dependent
+    ra1 = ra1[::-1]
+    dec1 = dec1[::-1]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = ra2[::-1]
+    dec2 = dec2[::-1]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+
+def test_hard_intersects_poly():
+    ra1 = [-10, 10, 10, -10, -10]
+    dec1 = [30, 30, 0, 0, 30]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = [-20, 20, 20, -20, -20]
+    dec2 = [20, 20, 10, 10, 20]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+    # Make sure it isn't order-dependent
+    ra1 = ra1[::-1]
+    dec1 = dec1[::-1]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = ra2[::-1]
+    dec2 = dec2[::-1]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert poly1.intersects_poly(poly2)
+
+
+def test_not_intersects_poly():
+    ra1 = [-10, 10, 10, -10, -10]
+    dec1 = [30, 30, 5, 5, 30]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = [-20, 20, 20, -20, -20]
+    dec2 = [-20, -20, -10, -10, -20]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert not poly1.intersects_poly(poly2)
+
+    # Make sure it isn't order-dependent
+    ra1 = ra1[::-1]
+    dec1 = dec1[::-1]
+    poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1)
+
+    ra2 = ra2[::-1]
+    dec2 = dec2[::-1]
+    poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2)
+
+    assert not poly1.intersects_poly(poly2)
+
+
+def test_point_in_poly():
+    point = np.asarray([-0.27475449, 0.47588873, -0.83548781])
+    points = np.asarray([[ 0.04821217, -0.29877206, 0.95310589],
+                         [ 0.04451801, -0.47274119, 0.88007608],
+                         [-0.14916503, -0.46369786, 0.87334649],
+                         [-0.16101648, -0.29210164, 0.94273555],
+                         [ 0.04821217, -0.29877206, 0.95310589]])
+    inside = np.asarray([-0.03416009, -0.36858623, 0.9289657])
+    poly = polygon.SphericalPolygon(points, inside)
+    assert not poly.contains_point(point)
+
+
+def test_point_in_poly_lots():
+    import pyfits
+    fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[0, 87])
+    poly2 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[20, 89])
+    poly3 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[180, 89])
+
+    points = get_point_set()
+    count = 0
+    for point in points:
+        if poly1.contains_point(point) or poly2.contains_point(point) or \
+               poly3.contains_point(point):
+            count += 1
+
+    assert count == 5
+    assert poly1.intersects_poly(poly2)
+    assert not poly1.intersects_poly(poly3)
+    assert not poly2.intersects_poly(poly3)
+
+
+def test_great_circle_arc_intersection():
+    A = [-10, -10]
+    B = [10, 10]
+
+    C = [-25, 10]
+    D = [15, -10]
+
+    E = [-20, 40]
+    F = [20, 40]
+
+    correct = [0.99912414, -0.02936109, -0.02981403]
+
+    A = vector.radec_to_vector(*A)
+    B = vector.radec_to_vector(*B)
+    C = vector.radec_to_vector(*C)
+    D = vector.radec_to_vector(*D)
+    E = vector.radec_to_vector(*E)
+    F = vector.radec_to_vector(*F)
+
+    assert great_circle_arc.intersects(A, B, C, D)
+    r = great_circle_arc.intersection(A, B, C, D)
+    assert r.shape == (3,)
+    assert_almost_equal(r, correct)
+
+    assert np.all(great_circle_arc.intersects([A], [B], [C], [D]))
+    r = great_circle_arc.intersection([A], [B], [C], [D])
+    assert r.shape == (1, 3)
+    assert_almost_equal(r, [correct])
+
+    assert np.all(great_circle_arc.intersects([A], [B], C, D))
+    r = great_circle_arc.intersection([A], [B], C, D)
+    assert r.shape == (1, 3)
+    assert_almost_equal(r, [correct])
+
+    assert not np.all(great_circle_arc.intersects([A, E], [B, F], [C], [D]))
+    r = great_circle_arc.intersection([A, E], [B, F], [C], [D])
+    assert r.shape == (2, 3)
+    assert_almost_equal(r[0], correct)
+    assert np.all(np.isnan(r[1]))
+
+    # Test parallel arcs
+    r = great_circle_arc.intersection(A, B, A, B)
+    assert np.all(np.isnan(r))
+
+
+def test_great_circle_arc_length():
+    A = [90, 0]
+    B = [-90, 0]
+    A = vector.radec_to_vector(*A)
+    B = vector.radec_to_vector(*B)
+    assert great_circle_arc.length(A, B) == 180.0
+
+    A = [135, 0]
+    B = [-90, 0]
+    A = vector.radec_to_vector(*A)
+    B = vector.radec_to_vector(*B)
+    assert great_circle_arc.length(A, B) == 135.0
+
+    A = [0, 0]
+    B = [0, 90]
+    A = vector.radec_to_vector(*A)
+    B = vector.radec_to_vector(*B)
+    assert great_circle_arc.length(A, B) == 90.0
+
+
+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
+
+    # TODO: More angle tests
+
+
+def test_cone():
+    random.seed(0)
+    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)
+
+
+def test_area():
+    triangles = [
+        ([[90, 0], [0, -45], [0, 45], [90, 0]], np.pi * 0.5),
+        ([[90, 0], [0, -22.5], [0, 22.5], [90, 0]], np.pi * 0.25),
+        ([[90, 0], [0, -11.25], [0, 11.25], [90, 0]], np.pi * 0.125)
+        ]
+
+    for tri, area in triangles:
+        tri = np.array(tri)
+        x, y, z = vector.radec_to_vector(tri[:,1], tri[:,0])
+        points = np.dstack((x, y, z))[0]
+        poly = polygon.SphericalPolygon(points)
+        calc_area = poly.area()
diff --git a/stsci/sphere/test/test_intersection.py b/stsci/sphere/test/test_intersection.py
new file mode 100644
index 0000000..6dacb3d
--- /dev/null
+++ b/stsci/sphere/test/test_intersection.py
@@ -0,0 +1,189 @@
+from __future__ import print_function, absolute_import
+
+# STDLIB
+import functools
+import itertools
+import math
+import os
+import random
+import sys
+
+# THIRD-PARTY
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+
+# LOCAL
+from .. import polygon
+from .test_shared import resolve_imagename
+
+GRAPH_MODE = False
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+
+class intersection_test:
+    def __init__(self, lon_0, lat_0, proj='ortho'):
+        self._lon_0 = lon_0
+        self._lat_0 = lat_0
+        self._proj = proj
+
+    def __call__(self, func):
+        @functools.wraps(func)
+        def run(*args, **kwargs):
+            if GRAPH_MODE:
+                from mpl_toolkits.basemap import Basemap
+                from matplotlib import pyplot as plt
+
+            polys = func(*args, **kwargs)
+
+            intersections = []
+            num_permutations = math.factorial(len(polys))
+            step_size = int(max(float(num_permutations) / 20.0, 1.0))
+
+            areas = [x.area() for x in polys]
+
+            if GRAPH_MODE:
+                print("%d permutations" % num_permutations)
+            for method in ('parallel', 'serial'):
+                for i, permutation in enumerate(
+                    itertools.islice(
+                        itertools.permutations(polys),
+                        None, None, step_size)):
+                    filename = '%s_%s_intersection_%04d.svg' % (
+                        func.__name__, method, i)
+                    print(filename)
+
+                    intersection = polygon.SphericalPolygon.multi_intersection(
+                        permutation, method=method)
+                    intersections.append(intersection)
+                    intersection_area = intersection.area()
+                    if GRAPH_MODE:
+                        fig = plt.figure()
+                        m = Basemap(projection=self._proj,
+                                    lon_0=self._lon_0,
+                                    lat_0=self._lat_0)
+                        m.drawparallels(np.arange(-90., 90., 20.))
+                        m.drawmeridians(np.arange(0., 420., 20.))
+                        m.drawmapboundary(fill_color='white')
+
+                        intersection.draw(m, color='red', linewidth=3)
+                        for poly in permutation:
+                            poly.draw(m, color='blue', alpha=0.5)
+                        plt.savefig(filename)
+                        fig.clear()
+
+                    assert np.all(intersection_area * 0.9 <= areas)
+
+            lengths = np.array([len(x._points) for x 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)
+
+        return run
+
+
+@intersection_test(0, 90)
+def test1():
+    import pyfits
+    import os
+
+    fits = pyfits.open(resolve_imagename(ROOT_DIR,'1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[0, 87])
+    poly2 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[20, 89])
+
+    return [poly1, poly2]
+
+
+@intersection_test(0, 90)
+def test2():
+    poly1 = polygon.SphericalPolygon.from_cone(0, 60, 8, steps=16)
+    poly2 = polygon.SphericalPolygon.from_cone(0, 68, 8, steps=16)
+    poly3 = polygon.SphericalPolygon.from_cone(12, 66, 8, steps=16)
+    return [poly1, poly2, poly3]
+
+
+@intersection_test(0, 90)
+def test3():
+    import pyfits
+    fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[0, 87])
+    poly3 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[175, 89])
+
+    return [poly1, poly3]
+
+
+def test4():
+    import pyfits
+    import pywcs
+
+    A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+    B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz'))
+
+    wcs = pywcs.WCS(A[1].header, fobj=A)
+    chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(A[4].header, fobj=A)
+    chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(B[1].header, fobj=B)
+    chipB1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(B[4].header, fobj=B)
+    chipB2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    Apoly = chipA1.union(chipA2)
+    Bpoly = chipB1.union(chipB2)
+
+    X = Apoly.intersection(Bpoly)
+
+
+def test5():
+    import pyfits
+    import pywcs
+
+    A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+    B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz'))
+
+    wcs = pywcs.WCS(A[1].header, fobj=A)
+    chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(A[4].header, fobj=A)
+    chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(B[1].header, fobj=B)
+    chipB1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(B[4].header, fobj=B)
+    chipB2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    Apoly = chipA1.union(chipA2)
+    Bpoly = chipB1.union(chipB2)
+
+    Apoly.overlap(chipB1)
+
+
+@intersection_test(0, 90)
+def test6():
+    import pyfits
+    fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1)
+    poly2 = polygon.SphericalPolygon.from_wcs(
+        header, 1)
+
+    return [poly1, poly2]
+
+
+if __name__ == '__main__':
+    if '--profile' not in sys.argv:
+        GRAPH_MODE = True
+        from mpl_toolkits.basemap import Basemap
+        from matplotlib import pyplot as plt
+
+    functions = [(k, v) for k, v in globals().items() if k.startswith('test')]
+    functions.sort()
+    for k, v in functions:
+        v()
diff --git a/stsci/sphere/test/test_shared.py b/stsci/sphere/test/test_shared.py
new file mode 100644
index 0000000..0853d3a
--- /dev/null
+++ b/stsci/sphere/test/test_shared.py
@@ -0,0 +1,12 @@
+import os
+
+def resolve_imagename(ROOT_DIR, base_name):
+    """Resolve image name for tests."""
+
+    image_name = os.path.join(ROOT_DIR, base_name)
+
+    # Is it zipped?
+    if not os.path.exists(image_name):
+        image_name = image_name.replace('.fits', '.fits.gz')
+
+    return image_name
diff --git a/stsci/sphere/test/test_skyline.py b/stsci/sphere/test/test_skyline.py
new file mode 100644
index 0000000..929f0ba
--- /dev/null
+++ b/stsci/sphere/test/test_skyline.py
@@ -0,0 +1,225 @@
+"""
+SkyLine tests.
+
+:Author: Pey Lian Lim
+
+:Organization: Space Telescope Science Institute
+
+Examples
+--------
+>>> cd path/to/project
+>>> nosetests
+
+"""
+from __future__ import absolute_import
+
+from copy import copy
+from numpy.testing import assert_almost_equal
+
+from ..vector import vector_to_radec
+from ..polygon import SphericalPolygon
+from ..skyline import SkyLine
+
+from .test_util import ROOT_DIR
+from .test_shared import resolve_imagename
+
+
+#---------------------------------------#
+# Load footprints used by all the tests #
+#---------------------------------------#
+f_2chipA  = resolve_imagename(ROOT_DIR, '2chipA.fits') # ACS/WFC #1
+im_2chipA = SkyLine(f_2chipA)
+f_2chipB  = resolve_imagename(ROOT_DIR, '2chipB.fits') # ACS/WFC #2
+im_2chipB = SkyLine(f_2chipB)
+f_2chipC  = resolve_imagename(ROOT_DIR, '2chipC.fits') # WFC3/UVIS
+im_2chipC = SkyLine(f_2chipC)
+f_66_tan  = resolve_imagename(ROOT_DIR, '1904-66_TAN.fits')
+im_66_tan = SkyLine(f_66_tan, extname='primary')
+
+
+#----- SHARED FUNCTIONS -----
+
+def same_members(mem1, mem2):
+    assert len(mem1) == len(mem2)
+
+    for m in mem1:
+        assert m in mem2
+
+    for m in mem2:
+        assert m in mem1
+
+def subset_members(mem_child, mem_parent):
+    assert len(mem_parent) > len(mem_child)
+
+    for m in mem_child:
+        assert m in mem_parent
+
+def subset_polygon(p_child, p_parent):
+    """Overlap not working. Do this instead until fixed."""
+    assert p_parent.area() >= p_child.area()
+    assert p_parent.contains_point(p_child.inside)
+
+def no_polygon(p_child, p_parent):
+    """Overlap not working. Do this instead until fixed."""
+    assert not p_parent.contains_point(p_child.inside)
+
+
+#----- MEMBERSHIP -----
+
+def do_member_overlap(im):
+    for m in im.members:
+        assert_almost_equal(m.polygon.overlap(im), 1.0)
+
+def test_membership():
+    do_member_overlap(im_2chipA)
+    do_member_overlap(im_2chipB)
+    do_member_overlap(im_66_tan)
+
+    assert len(im_2chipA.members) == 1
+    assert im_2chipA.members[0].fname == f_2chipA
+    assert im_2chipA.members[0].ext == (1,4)
+
+
+#----- COPY -----
+
+def test_copy():
+    a_copy = copy(im_2chipA)
+    assert a_copy is not im_2chipA
+
+
+#----- SPHERICAL POLYGON RELATED -----
+
+def test_sphericalpolygon():
+    assert im_2chipA.contains_point(im_2chipA.inside)
+    
+    assert im_2chipA.intersects_poly(im_2chipB.polygon)
+    
+    assert im_2chipA.intersects_arc(im_2chipA.inside, im_2chipB.inside)
+    
+    assert im_2chipA.overlap(im_2chipB) < im_2chipA.overlap(im_2chipA)
+    
+    assert_almost_equal(im_2chipA.area(), im_2chipB.area())
+
+    ra_A, dec_A = im_2chipA.to_radec()
+    for i in xrange(len(im_2chipA.points)):
+        p = im_2chipA.points[i]
+        ra, dec = vector_to_radec(p[0], p[1], p[2], degrees=True)
+        assert_almost_equal(ra_A[i], ra)
+        assert_almost_equal(dec_A[i], dec)
+
+
+#----- WCS -----
+
+def test_wcs():
+    wcs = im_2chipA.to_wcs()
+    new_p = SphericalPolygon.from_wcs(wcs)
+    subset_polygon(im_2chipA, new_p)
+
+
+#----- UNION -----
+
+def do_add_image(im1, im2):
+    u1 = im1.add_image(im2)
+    u2 = im2.add_image(im1)
+
+    assert u1.same_points_as(u2)
+    same_members(u1.members, u2.members)
+
+    all_mems = im1.members + im2.members
+    same_members(u1.members, all_mems)
+
+    subset_polygon(im1, u1)
+    subset_polygon(im2, u1)
+
+def test_add_image():
+    # Dithered
+    do_add_image(im_2chipA, im_2chipB)
+
+    # Not related
+    do_add_image(im_2chipA, im_66_tan)
+
+
+#----- INTERSECTION -----
+
+def do_intersect_image(im1, im2):
+    i1 = im1.find_intersection(im2)
+    i2 = im2.find_intersection(im1)
+
+    assert i1.same_points_as(i2)
+    same_members(i1.members, i2.members)
+
+    if len(i1.points) > 0:
+        subset_members(im1.members, i1.members)
+        subset_members(im2.members, i1.members)
+    
+        subset_polygon(i1, im1)
+        subset_polygon(i1, im2)
+
+def test_find_intersection():   
+    # Dithered
+    do_intersect_image(im_2chipA, im_2chipB)
+
+    # Not related
+    do_intersect_image(im_2chipA, im_66_tan)
+
+
+#----- SKYLINE OVERLAP -----
+
+def test_max_overlap():
+    max_s, max_a = im_2chipA.find_max_overlap([im_2chipB, im_2chipC, im_66_tan])
+    assert max_s is im_2chipB
+    assert_almost_equal(max_a, im_2chipA.intersection(im_2chipB).area())
+
+    max_s, max_a = im_2chipA.find_max_overlap([im_2chipB, im_2chipA])
+    assert max_s is im_2chipA
+    assert_almost_equal(max_a, im_2chipA.area())
+
+def test_max_overlap_pair():
+    assert SkyLine.max_overlap_pair(
+        [im_2chipB, im_2chipC, im_2chipA, im_66_tan]) == (im_2chipB, im_2chipA)
+
+    assert SkyLine.max_overlap_pair([im_2chipC, im_2chipA, im_66_tan]) is None
+
+
+#----- INTENDED USE CASE -----
+
+def test_science_1():
+    mos, inc, exc = SkyLine.mosaic([im_2chipA, im_2chipB, im_2chipC, im_66_tan])
+
+    assert inc == [f_2chipA, f_2chipB]
+    assert exc == [f_2chipC, f_66_tan]
+
+    subset_polygon(im_2chipA, mos)
+    subset_polygon(im_2chipB, mos)
+
+    no_polygon(im_2chipC, mos)
+    no_polygon(im_66_tan, mos)
+
+def test_science_2():
+    """Like `test_science_1` but different input order."""
+    mos, inc, exc = SkyLine.mosaic([im_2chipB, im_66_tan, im_2chipC, im_2chipA])
+
+    assert inc == [f_2chipB, f_2chipA]
+    assert exc == [f_66_tan, f_2chipC]
+
+    subset_polygon(im_2chipA, mos)
+    subset_polygon(im_2chipB, mos)
+
+    no_polygon(im_2chipC, mos)
+    no_polygon(im_66_tan, mos)
+
+
+#----- UNSTABLE -----
+
+def DISABLED_unstable_overlap():
+    i1 = im_2chipA.find_intersection(im_2chipB)
+    i2 = im_2chipB.find_intersection(im_2chipA)
+    
+    u1 = im_2chipA.add_image(im_2chipB)
+    u2 = im_2chipB.add_image(im_2chipA)
+
+    # failed here before - known bug
+    # failure not always the same due to hash mapping
+    assert_almost_equal(i1.overlap(u1), 1.0)
+    assert_almost_equal(i1.overlap(i2), 1.0)
+    assert_almost_equal(u1.overlap(u2), 1.0)
diff --git a/stsci/sphere/test/test_union.py b/stsci/sphere/test/test_union.py
new file mode 100644
index 0000000..4b649f6
--- /dev/null
+++ b/stsci/sphere/test/test_union.py
@@ -0,0 +1,215 @@
+from __future__ import print_function, absolute_import
+
+# STDLIB
+import functools
+import itertools
+import math
+import os
+import random
+import sys
+
+# THIRD-PARTY
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+
+# LOCAL
+from .. import polygon
+from .test_shared import resolve_imagename
+
+GRAPH_MODE = False
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+
+class union_test:
+    def __init__(self, lon_0, lat_0, proj='ortho'):
+        self._lon_0 = lon_0
+        self._lat_0 = lat_0
+        self._proj = proj
+
+    def __call__(self, func):
+        @functools.wraps(func)
+        def run(*args, **kwargs):
+            if GRAPH_MODE:
+                from mpl_toolkits.basemap import Basemap
+                from matplotlib import pyplot as plt
+
+            polys = func(*args, **kwargs)
+
+            unions = []
+            num_permutations = math.factorial(len(polys))
+            step_size = int(max(float(num_permutations) / 20.0, 1.0))
+            if GRAPH_MODE:
+                print("%d permutations" % num_permutations)
+
+            areas = [x.area() for x in polys]
+
+            for i, permutation in enumerate(
+                    itertools.islice(
+                    itertools.permutations(polys),
+                    None, None, step_size)):
+                filename = '%s_union_%04d.svg' % (
+                    func.__name__, i)
+                print(filename)
+
+                union = polygon.SphericalPolygon.multi_union(
+                    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,
+                                lon_0=self._lon_0,
+                                lat_0=self._lat_0)
+                    m.drawmapboundary(fill_color='white')
+                    m.drawparallels(np.arange(-90., 90., 20.))
+                    m.drawmeridians(np.arange(0., 420., 20.))
+
+                    union.draw(m, color='red', linewidth=3)
+                    for poly in permutation:
+                        poly.draw(m, color='blue', alpha=0.5)
+                    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])
+            assert np.all(lengths == [lengths[0]])
+            areas = np.array([x.area() for x in unions])
+            assert_array_almost_equal(areas, areas[0], 1)
+
+        return run
+
+
+@union_test(0, 90)
+def test1():
+    import pyfits
+    fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[0, 87])
+    poly2 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[20, 89])
+    poly3 = polygon.SphericalPolygon.from_wcs(
+        header, 1, crval=[175, 89])
+    poly4 = polygon.SphericalPolygon.from_cone(
+        90, 70, 10, steps=50)
+
+    return [poly1, poly2, poly3, poly4]
+
+
+@union_test(0, 90)
+def test2():
+    poly1 = polygon.SphericalPolygon.from_cone(0, 60, 7, steps=16)
+    poly2 = polygon.SphericalPolygon.from_cone(0, 72, 7, steps=16)
+    poly3 = polygon.SphericalPolygon.from_cone(20, 60, 7, steps=16)
+    poly4 = polygon.SphericalPolygon.from_cone(20, 72, 7, steps=16)
+    poly5 = polygon.SphericalPolygon.from_cone(35, 55, 7, steps=16)
+    poly6 = polygon.SphericalPolygon.from_cone(60, 60, 3, steps=16)
+    return [poly1, poly2, poly3, poly4, poly5, poly6]
+
+
+@union_test(0, 90)
+def test3():
+    random.seed(0)
+    polys = []
+    for i in range(10):
+        polys.append(polygon.SphericalPolygon.from_cone(
+            random.randrange(-180, 180),
+            random.randrange(20, 90),
+            random.randrange(5, 16),
+            steps=16))
+    return polys
+
+
+@union_test(0, 15)
+def test4():
+    random.seed(64)
+    polys = []
+    for i in range(10):
+        polys.append(polygon.SphericalPolygon.from_cone(
+            random.randrange(-30, 30),
+            random.randrange(-15, 60),
+            random.randrange(5, 16),
+            steps=16))
+    return polys
+
+
+def test5():
+    import pyfits
+    import pywcs
+
+    A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+
+    wcs = pywcs.WCS(A[1].header, fobj=A)
+    chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(A[4].header, fobj=A)
+    chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    null_union = chipA1.union(chipA2)
+
+
+def test6():
+    import pyfits
+    import pywcs
+
+    A = pyfits.open(os.path.join(ROOT_DIR, '2chipC.fits.gz'))
+
+    wcs = pywcs.WCS(A[1].header, fobj=A)
+    chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(A[4].header, fobj=A)
+    chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    null_union = chipA1.union(chipA2)
+
+
+@union_test(0, 90)
+def test7():
+    import pyfits
+    import pywcs
+
+    A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz'))
+
+    wcs = pywcs.WCS(A[1].header, fobj=A)
+    chipA1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(A[4].header, fobj=A)
+    chipA2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz'))
+
+    wcs = pywcs.WCS(B[1].header, fobj=B)
+    chipB1 = polygon.SphericalPolygon.from_wcs(wcs)
+    wcs = pywcs.WCS(B[4].header, fobj=B)
+    chipB2 = polygon.SphericalPolygon.from_wcs(wcs)
+
+    return [chipA1, chipA2, chipB1, chipB2]
+
+
+@union_test(0, 90)
+def test8():
+    import pyfits
+    fits = pyfits.open(resolve_imagename(ROOT_DIR, '1904-66_TAN.fits'))
+    header = fits[0].header
+
+    poly1 = polygon.SphericalPolygon.from_wcs(
+        header, 1)
+    poly2 = polygon.SphericalPolygon.from_wcs(
+        header, 1)
+
+    return [poly1, poly2]
+
+
+if __name__ == '__main__':
+    if '--profile' not in sys.argv:
+        GRAPH_MODE = True
+        from mpl_toolkits.basemap import Basemap
+        from matplotlib import pyplot as plt
+
+    functions = [(k, v) for k, v in globals().items() if k.startswith('test')]
+    functions.sort()
+    for k, v in functions:
+        v()
diff --git a/stsci/sphere/test/test_util.py b/stsci/sphere/test/test_util.py
new file mode 100644
index 0000000..6dc0393
--- /dev/null
+++ b/stsci/sphere/test/test_util.py
@@ -0,0 +1,14 @@
+import os
+
+import numpy as np
+from .. import vector
+
+ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data')
+
+def get_point_set(density=25):
+    points = []
+    for i in np.linspace(-85, 85, density, True):
+        for j in np.linspace(-180, 180, np.cos(np.deg2rad(i)) * density):
+            points.append([j, i])
+    points = np.asarray(points)
+    return np.dstack(vector.radec_to_vector(points[:,0], points[:,1]))[0]
diff --git a/stsci/sphere/vector.py b/stsci/sphere/vector.py
new file mode 100644
index 0000000..c6887c4
--- /dev/null
+++ b/stsci/sphere/vector.py
@@ -0,0 +1,243 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2011 Association of Universities for Research in
+# Astronomy (AURA)
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+#     1. Redistributions of source code must retain the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer.
+#
+#     2. Redistributions in binary form must reproduce the above
+#       copyright notice, this list of conditions and the following
+#       disclaimer in the documentation and/or other materials
+#       provided with the distribution.
+#
+#     3. The name of AURA and its representatives may not be used to
+#       endorse or promote products derived from this software without
+#       specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR
+# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+# OF THE POSSIBILITY OF SUCH DAMAGE.
+
+"""
+The `sphere.vector` module contains the basic operations for handling
+vectors and converting them to and from other representations.
+"""
+
+from __future__ import unicode_literals
+
+# THIRD-PARTY
+import numpy as np
+
+
+__all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector',
+           'rotate_around']
+
+
+def radec_to_vector(ra, dec, degrees=True):
+    r"""
+    Converts a location on the unit sphere from right-ascension and
+    declination to an *x*, *y*, *z* vector.
+
+    Parameters
+    ----------
+    ra, dec : scalars or 1-D arrays
+
+    degrees : bool, optional
+
+       If `True`, (default) *ra* and *dec* are in decimal degrees,
+       otherwise in radians.
+
+    Returns
+    -------
+    x, y, z : tuple of scalars or 1-D arrays of the same length
+
+    Notes
+    -----
+    Where right-ascension is *α* and declination is *δ*:
+
+    .. math::
+        x = \cos\alpha \cos\delta
+
+        y = \sin\alpha \cos\delta
+
+        z = \sin\delta
+    """
+    ra = np.asanyarray(ra)
+    dec = np.asanyarray(dec)
+
+    if degrees:
+        ra_rad = np.deg2rad(ra)
+        dec_rad = np.deg2rad(dec)
+    else:
+        ra_rad = ra
+        dec_rad = dec
+
+    cos_dec = np.cos(dec_rad)
+
+    return (
+        np.cos(ra_rad) * cos_dec,
+        np.sin(ra_rad) * cos_dec,
+        np.sin(dec_rad))
+
+
+def vector_to_radec(x, y, z, degrees=True):
+    r"""
+    Converts a vector to right-ascension and declination.
+
+    Parameters
+    ----------
+    x, y, z : scalars or 1-D arrays
+        The input vectors
+
+    degrees : bool, optional
+        If `True` (default) the result is returned in decimal degrees,
+        otherwise radians.
+
+    Returns
+    -------
+    ra, dec : tuple of scalars or arrays of the same length
+
+    Notes
+    -----
+    Where right-ascension is *α* and declination is
+    *δ*:
+
+    .. math::
+        \alpha = \arctan2(y, x)
+
+        \delta = \arctan2(z, \sqrt{x^2 + y^2})
+    """
+    x = np.asanyarray(x)
+    y = np.asanyarray(y)
+    z = np.asanyarray(z)
+
+    result = (
+        np.arctan2(y, x),
+        np.arctan2(z, np.sqrt(x ** 2 + y ** 2)))
+
+    if degrees:
+        return np.rad2deg(result[0]), np.rad2deg(result[1])
+    else:
+        return result
+
+
+def normalize_vector(x, y, z, inplace=False):
+    r"""
+    Normalizes a vector so it falls on the unit sphere.
+
+    *x*, *y*, *z* may be scalars or 1-D arrays
+
+    Parameters
+    ----------
+    x, y, z : scalars or 1-D arrays of the same length
+        The input vectors
+
+    inplace : bool, optional
+        When `True`, the original arrays will be normalized in place,
+        otherwise a normalized copy is returned.
+
+    Returns
+    -------
+    X, Y, Z : scalars of 1-D arrays of the same length
+        The normalized output vectors
+    """
+    x = np.asanyarray(x)
+    y = np.asanyarray(y)
+    z = np.asanyarray(z)
+
+    l = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+
+    if inplace:
+        x /= l
+        y /= l
+        z /= l
+        return (x, y, z)
+    else:
+        return (x / l, y / l, z / l)
+
+
+def rotate_around(x, y, z, u, v, w, theta, degrees=True):
+    r"""
+    Rotates the vector (*x*, *y*, *z*) around the arbitrary axis defined by
+    vector (*u*, *v*, *w*) by *theta*.
+
+    It is assumed that both (*x*, *y*, *z*) and (*u*, *v*, *w*) are
+    already normalized.
+
+    Parameters
+    ----------
+    x, y, z : doubles
+        The normalized vector to rotate
+
+    u, v, w : doubles
+        The normalized vector to rotate around
+
+    theta : double, or array of doubles
+        The amount to rotate
+
+    degrees : bool, optional
+        When `True`, *theta* is given in degrees, otherwise radians.
+
+    Returns
+    -------
+    X, Y, Z : doubles
+        The rotated vector
+    """
+    if degrees:
+        theta = np.deg2rad(theta)
+
+    costheta = np.cos(theta)
+    sintheta = np.sin(theta)
+    icostheta = 1.0 - costheta
+
+    det = (-u*x - v*y - w*z)
+    X = (-u*det)*icostheta + x*costheta + (-w*y + v*z)*sintheta
+    Y = (-v*det)*icostheta + y*costheta + ( w*x - u*z)*sintheta
+    Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta
+
+    return X, Y, Z
+
+
+def equal_area_proj(x, y, z):
+    """
+    Transform the coordinates to a 2-dimensional plane using the
+    Lambert azimuthal equal-area projection.
+
+    Parameters
+    ----------
+    x, y, z : scalars or 1-D arrays
+        The input vectors
+
+    Returns
+    -------
+    X, Y : tuple of scalars or arrays of the same length
+
+    Notes
+    -----
+
+    .. math::
+
+        X = \sqrt{\frac{2}{1-z}}x
+
+    .. math::
+
+        Y = \sqrt{\frac{2}{1-z}}y
+    """
+    scale = np.sqrt(2.0 / (1.0 - z))
+    X = scale * x
+    Y = scale * y
+    return X, Y