12 files changed, 972 insertions, 0 deletions
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]