diff options
| -rw-r--r-- | lib/polygon.py | 97 | ||||
| -rw-r--r-- | lib/skyline.py | 11 | ||||
| -rw-r--r-- | lib/test/test_skyline.py | 140 |
3 files changed, 222 insertions, 26 deletions
diff --git a/lib/polygon.py b/lib/polygon.py index fc4dad7..e0bf841 100644 --- a/lib/polygon.py +++ b/lib/polygon.py @@ -113,17 +113,17 @@ class SphericalPolygon(object): """ return self._inside - @property - def radec(self): + def to_radec(self): """Convert `SphericalPolygon` footprint to RA and DEC. Returns ------- - ra, dec: array_like - RA and DEC of `self.points`. + List of RA and DEC pairs (in degrees) corresponding to + `self.points`. """ - return vector.vector_to_radec(self.points[:,0], self.points[:,1], - self.points[:,2], degrees=True) + ra, dec = vector.vector_to_radec(self.points[:,0], self.points[:,1], + self.points[:,2], degrees=True) + return zip(ra,dec) @classmethod def from_radec(cls, ra, dec, center=None, degrees=True): @@ -288,6 +288,89 @@ class SphericalPolygon(object): return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0])) + def _unique_points(self): + """Ignore duplicate `points`. Order not preserved.""" + # http://stackoverflow.com/questions/7989722/finding-unique-points-in-numpy-array + if len(self.points) == 0: + val = [] + else: + val = np.vstack([np.array(u) for u in + set([tuple(p) for p in self.points])]) + return val + + def _sorted_points(self, unique=True): + """ + Sort `points` in the order of *x*, *y*, and *z*. + + Parameters + ---------- + unique : bool + Ignore duplicates. + """ + if len(self.points) == 0: + return [] + + if unique: + pts = self._unique_points() + else: + pts = self.points.copy() + + return pts[np.lexsort((pts[:,0], pts[:,1], pts[:,2]))] + + def same_points_as(self, other, do_sort=True, thres=0.01): + """ + Determines if this `SphericalPolygon` points are the same + as the other. + + Parameters + ---------- + other : `SphericalPolygon` + + do_sort : bool + Compare unique sorted points. + + thres : float + Fraction of area to use in equality decision. + + Returns + ------- + `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() + other_pts = other._sorted_points() + 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. @@ -422,7 +505,7 @@ class SphericalPolygon(object): A = \left| \sum^n_{i=0} X_i Y_{i+1} - X_{i+1}Y_i \right| """ if len(self._points) < 3: - return np.array(0.0) + return np.float64(0.0) points = self._points.copy() diff --git a/lib/skyline.py b/lib/skyline.py index bcc3bae..2cc70c1 100644 --- a/lib/skyline.py +++ b/lib/skyline.py @@ -161,6 +161,8 @@ class SkyLine(object): HST images. PRIMARY if image is single ext. """ + extname = extname.upper() + # Convert SCI data to SkyLineMember if fname is not None: with pyfits.open(fname) as pf: @@ -177,7 +179,14 @@ class SkyLine(object): self.polygon = SphericalPolygon([]) def __getattr__(self, what): - return getattr(self.polygon, what) + """Control attribute access to `SphericalPolygon`.""" + if what in ('from_radec', 'from_cone', 'from_wcs', + 'multi_union', 'multi_intersection', + '_find_new_inside',): + raise AttributeError('\'SkyLine\' object has no attribute \'%s\'' % + what) + else: + return getattr(self.polygon, what) def __copy__(self): return deepcopy(self) diff --git a/lib/test/test_skyline.py b/lib/test/test_skyline.py index d450e49..ef1fce5 100644 --- a/lib/test/test_skyline.py +++ b/lib/test/test_skyline.py @@ -1,34 +1,138 @@ -"""SkyLine tests. +""" +SkyLine tests. :Author: Pey Lian Lim :Organization: Space Telescope Science Institute +Examples +-------- +>>> cd path/to/project +>>> nosetests + """ from __future__ import absolute_import -import pyfits -from numpy.testing import assert_almost_equal, assert_array_less +from copy import copy +from numpy.testing import assert_almost_equal -from .. import skyline +from ..skyline import SkyLine from .test_util import ROOT_DIR from .test_shared import resolve_imagename -def test_union_simple(): - # Two similar exposures with slight offset and some rotation. - im1 = resolve_imagename(ROOT_DIR, '2chipA.fits') - im2 = resolve_imagename(ROOT_DIR, '2chipB.fits') - - skyline1 = skyline.SkyLine(im1) - skyline2 = skyline.SkyLine(im2) - - union_1_2 = skyline1.union(skyline2) - union_2_1 = skyline2.union(skyline1) - assert_almost_equal(union_1_2.area(), union_2_1.area()) +#---------------------------------------# +# 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 compare_members(im1, im2): + assert len(im1.members) == len(im2.members) + + for m in im1.members: + assert m in im2.members + + for m in im2.members: + assert m in im1.members + + +#----- 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) == 2 + assert im_2chipA.members[0].fname == f_2chipA + assert im_2chipA.members[0].ext == 1 + assert im_2chipA.members[1].fname == f_2chipA + assert im_2chipA.members[1].ext == 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()) - for m in union_1_2.members: - assert m in union_2_1.members - assert len(union_1_2.members) == 4 +#----- UNION ----- + +def do_add_image(im1, im2): + u1 = im1.add_image(im2) + u2 = im2.add_image(im1) + + assert u1.same_points_as(u2) + compare_members(u1, u2) + +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) + compare_members(i1, i2) + +def test_find_intersection(): + # Dithered + do_intersect_image(im_2chipA, im_2chipB) + + # Not related + do_intersect_image(im_2chipA, im_66_tan) + + +# ----- INTENDED USE CASE ----- + +def test_science(): + skylines = [im_2chipA, im_2chipB, im_2chipC, im_66_tan] + + # TODO: Add Warren's example use case + + +# ----- 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) + + assert_almost_equal(i1.overlap(u1), 1.0) + assert_almost_equal(i1.overlap(i2), 1.0) # failed here - known bug + assert_almost_equal(u1.overlap(u2), 1.0) |