working test of is_in_window

docs/pydl/changes.rst

@@ -8,8 +8,9 @@ PyDL Changelog
 * Changes in how Mangle-polygon containing FITS files are handled, related to
   `Issue #11`_.
 * Added :func:`~pydl.pydlutils.mangle.is_in_window`.
+* Allow polygon area functions to deal with negative caps and ``use_caps``.
 
-.. _`Issue #8`: https://github.com/weaverba137/pydl/issues/11
+.. _`Issue #11`: https://github.com/weaverba137/pydl/issues/11
 
 0.5.1 (2016-06-22)
 ------------------

pydl/pydlutils/mangle.py

@@ -11,8 +11,8 @@
 idlutils, the caps information is accessed through ``polygon.x`` and
 ``polygon.cm``, not ``polygon.caps.x`` or ``polygon.caps.cm``.
 
-Window function operations are already supported by :func:`~is_in_cap` and
-:func:`~is_in_polygon`.  However, calculation of polygon area (solid angle)
+Window function operations are already supported by :func:`~is_in_polygon` and
+:func:`~is_in_window`.  However, calculation of polygon area (solid angle)
 from first principles is not yet implemented.
 
 Note that in traditional geometry "spherical polygon" means a figure
@@ -188,7 +188,7 @@ def str(self):
 
     def cmminf(self):
         """The index of the smallest cap in the polygon, accounting for
-        negative caps.
+        negative caps and ``use_caps``.
 
         Returns
         -------
@@ -200,13 +200,14 @@ def cmminf(self):
         cmmin = 2.0
         kmin = -1
         for k in range(self.ncaps):
-            if self.cm[k] >= 0:
-                cmk = self.cm[k]
-            else:
-                cmk = 2.0 + self.cm[k]
-            if cmk < cmmin:
-                cmmin = cmk
-                kmin = k
+            if is_cap_used(self.use_caps, k):
+                if self.cm[k] >= 0:
+                    cmk = self.cm[k]
+                else:
+                    cmk = 2.0 + self.cm[k]
+                if cmk < cmmin:
+                    cmmin = cmk
+                    kmin = k
         return kmin
 
     def garea(self):
@@ -236,16 +237,26 @@ def garea(self):
         if smallest_cap is None:
             return 4.0 * np.pi
         cmmin = self.cm[smallest_cap]
-        if self.ncaps >= 2 and cmmin > 1.0:
-            npl = self.ncaps + 1
+        used_caps = list()
+        for k in range(self.ncaps):
+            if is_cap_used(self.use_caps, k):
+                used_caps.append(k)
+        nused_caps = len(used_caps)
+        if nused_caps == 0:
+            return 0.0
+        if nused_caps >= 2 and cmmin > 1.0:
+            npl = nused_caps + 1
         else:
-            npl = self.ncaps
-        if npl == self.ncaps:
-            if self.ncaps == 1:
+            npl = nused_caps
+        if npl == nused_caps:
+            if nused_caps == 1:
                 #
                 # Short-circuit this case.
                 #
-                return 2.0*np.pi*cmmin
+                if cmmin < 0:
+                    return 2.0*np.pi*(2.0+cmmin)
+                else:
+                    return 2.0*np.pi*cmmin
             else:
                 #
                 # Two or more caps, at least one has area < 2.0*pi.
@@ -254,11 +265,12 @@ def garea(self):
         else:
             #
             # More than two caps, and all have area > 2.0*pi.
+            # This section still needs to account for use_caps.
             #
             dpoly = self.polyn(self, smallest_cap)
-            dpoly.cm[self.ncaps] = cmmin / 2.0
+            dpoly.cm[nused_caps] = cmmin / 2.0
             area1 = dpoly._garea_helper()
-            dpoly.cm[self.ncaps] = -1.0 * dpoly.cm[self.ncaps]
+            dpoly.cm[nused_caps] = -1.0 * dpoly.cm[nused_caps]
             area2 = dpoly._garea_helper()
             return area1 + area2
 

pydl/pydlutils/tests/test_mangle.py

@@ -152,10 +152,24 @@ def test_is_in_polygon(self):
         d = mng.is_in_polygon(p, np.array([[45.0, 45.0], [-45.0, 45.0]]),
                               ncaps=2)
         assert (d == np.array([True, False])).all()
+        poly = mng.read_fits_polygons(self.no_id_fits)
+        d = mng.is_in_polygon(poly[0],
+                              np.array([[0.0, 10.0],
+                                        [4.0, 4.5],
+                                        [90, 4.0],
+                                        [135, 3.0],
+                                        [180, 0.5],
+                                        [270, 0.25]]))
+        assert (d == np.array([False, True, False, False, False, False])).all()
 
     def test_is_in_window(self):
-        # poly = mng.read_fits_polygons(self.poly_fits)
-        pass
+        np.random.seed(271828)
+        RA = 7.0*np.random.random(1000) + 268.0
+        Dec = 90.0 - np.degrees(np.arccos(0.08*np.random.random(1000)))
+        points = np.vstack((RA, Dec)).T
+        poly = mng.read_fits_polygons(self.poly_fits)
+        i = mng.is_in_window(poly, points)
+        assert i[0].sum() == 3
 
     def test_read_fits_polygons(self):
         poly = mng.read_fits_polygons(self.poly_fits)