Merge pull request #8 from keflavich/mybranch

More changes I've made while applying dendrocat

dendrocat/aperture.py

@@ -1,6 +1,7 @@
 import regions
 import astropy.units as u
 from astropy.coordinates import SkyCoord
+import astropy.wcs
 import numpy as np
 import warnings
 
@@ -66,6 +67,9 @@ def __init__(self, center, major, minor, pa, unit=None, frame='icrs',
         self.pa = ucheck(pa, u.deg)
         self.frame = frame
 
+        if hasattr(self.center, 'frame'):
+            assert self.center.frame.name == self.frame
+
     def _refresh_xycen(self):
         if type(self.center) == SkyCoord:
             self.x_cen = ucheck(self.center.spherical.lon, self.unit)
@@ -97,6 +101,8 @@ def place(self, image, wcs=None):
         numpy.ndarray
             A boolean mask for the aperture with the same dimensions as `image`
         """
+        if wcs is not None and self.frame != astropy.wcs.utils.wcs_to_celestial_frame(wcs).name:
+            raise ValueError("Frame mismatch in aperture placement")
         self._refresh_xycen()
         if self.unit.is_equivalent(u.deg) and wcs is not None:
             pixel_scale = (np.abs(wcs.pixel_scale_matrix.diagonal()
@@ -163,7 +169,8 @@ def __init__(self, center, major, minor, pa, unit=None, frame='icrs', name=None)
             The name used in the catalog column names when photometry is
             performed with this aperture.
         """
-        Aperture.__init__(self, center, major, minor, pa, unit=unit, name=name)
+        Aperture.__init__(self, center, major, minor, pa, unit=unit, name=name,
+                          frame=frame)
 
     def place(self, image, wcs=None):
         """
@@ -224,11 +231,15 @@ def __init__(self, center, inner, outer, unit=None, frame='icrs', name=None):
         if name is not None:
             self.__name__ = name
 
-        self.aperture_inner = Aperture(center, inner, inner, 0, unit=unit)
-        self.aperture_outer = Aperture(center, outer, outer, 0, unit=unit)
+        self.aperture_inner = Aperture(center, inner, inner, 0, unit=unit, frame=frame)
+        self.aperture_outer = Aperture(center, outer, outer, 0, unit=unit, frame=frame)
         self.center = ucheck(center, self.unit)
         self.x_cen = ucheck(center[0], self.unit)
         self.y_cen = ucheck(center[1], self.unit)
+        self.frame = frame
+
+        if hasattr(self.center, 'frame'):
+            assert self.center.frame.name == self.frame
 
     @property
     def inner(self):
@@ -303,7 +314,8 @@ def __init__(self, center, radius, unit=None, frame='icrs', name=None):
             The name used in the catalog column names when photometry is
             performed with this aperture.
         """
-        Aperture.__init__(self, center, radius, radius, 0, unit=unit, name=name)
+        Aperture.__init__(self, center, radius, radius, 0, unit=unit,
+                          name=name, frame=frame)
         self.radius = ucheck(radius, self.unit)
 
     def place(self, image, wcs=None):

dendrocat/radiosource.py

@@ -284,10 +284,11 @@ def _make_cutouts(self, catalog=None, data=None, save=True):
 
             position = coordinates.SkyCoord(x_cen,
                                             y_cen,
-                                            frame='icrs',
+                                            frame=wcs.utils.wcs_to_celestial_frame(self.wcs).name,
                                             unit=(u.deg, u.deg))
 
-            pixel_position = np.array(position.to_pixel(self.wcs))
+            # commented out b/c not used
+            # pixel_position = np.array(position.to_pixel(self.wcs))
 
             try:
                 cutout = Cutout2D(data,
@@ -362,6 +363,8 @@ def get_pixels(self, aperture, catalog=None, data=None, cutouts=None,
                 masks.append(float('nan'))
                 continue
 
+            frame = wcs.utils.wcs_to_celestial_frame(cutouts[i].wcs).name
+
             x_cen = catalog['x_cen'][i]
             y_cen = catalog['y_cen'][i]
             major = catalog['major_fwhm'][i]
@@ -374,9 +377,12 @@ def get_pixels(self, aperture, catalog=None, data=None, cutouts=None,
                 # replace the center value with the centers from the sources.
 
                 if aperture.unit.is_equivalent(u.deg):
-                    aperture.center = coordinates.SkyCoord(x_cen*u.deg, y_cen*u.deg)
+                    aperture.center = coordinates.SkyCoord(x_cen*u.deg,
+                                                           y_cen*u.deg,
+                                                           frame=frame)
                 elif aperture.unit.is_equivalent(u.pix):
-                    sky = coordinates.SkyCoord(x_cen*u.deg, y_cen*u.deg)
+                    sky = coordinates.SkyCoord(x_cen*u.deg, y_cen*u.deg,
+                                               frame=frame)
                     pixel = ucheck(sky.to_pixel(cutouts[i].wcs), u.pix)
                     aperture.center = pixel
                     aperture.x_cen, aperture.y_cen = pixel[0], pixel[1]
@@ -388,23 +394,26 @@ def get_pixels(self, aperture, catalog=None, data=None, cutouts=None,
                 # DEFAULTS FOR VARIABLE APERTURES STORED HERE
                 cen = [x_cen, y_cen]
                 if aperture == Ellipse:
-                    aperture = Ellipse(cen, major, minor, pa, unit=u.deg)
+                    aperture = Ellipse(cen, major, minor, pa, unit=u.deg,
+                                       frame=frame)
 
                 elif aperture == Annulus:
                     inner_r = major*u.deg+self.annulus_padding
                     outer_r = major*u.deg+self.annulus_padding+self.annulus_width
-                    aperture = Annulus(cen, inner_r, outer_r, unit=u.deg)
+                    aperture = Annulus(cen, inner_r, outer_r, unit=u.deg, frame=frame)
 
                 elif aperture == Circle:
                     radius = major
-                    aperture = Circle(cen, radius, unit=u.deg)
+                    aperture = Circle(cen, radius, unit=u.deg, frame=frame)
 
                 else:
                     raise UnknownApertureError('Aperture not recognized. Pass'
                                                ' an instance of a custom aper'
                                                'ture instead.')
 
             this_mask = aperture.place(cutouts[i].data, wcs=cutouts[i].wcs)
+            if this_mask.sum() == 0:
+                raise ValueError("No pixels within aperture")
             pix_arrays.append(cutouts[i].data[this_mask])
             masks.append(this_mask)
             aperture = aperture_original # reset the aperture for the next source
@@ -551,10 +560,8 @@ def plot_grid(self, catalog=None, data=None, cutouts=None,
         masks = []
 
         for aperture in [source_aperture, bkg_aperture]:
-            some_pixels, a_mask = self.get_pixels(aperture,
-                                                  catalog=catalog,
-                                                  data=data,
-                                                  cutouts=cutouts)
+            some_pixels, a_mask = self.get_pixels(aperture, catalog=catalog,
+                                                  data=data, cutouts=cutouts)
             ap_names.append(aperture.__name__)
             pixels.append(some_pixels)
             masks.append(a_mask)
@@ -582,8 +589,9 @@ def plot_grid(self, catalog=None, data=None, cutouts=None,
         an = np.ones(len(cutouts), dtype='bool')
         for i in range(len(cutouts)):
             try:
+                # check whether cutouts[i] is a cutout or is NaN
                 np.isnan(cutouts[i])
-                an[i] = 0
+                an[i] = False
             except TypeError:
                 pass