diff --git a/dendrocat/aperture.py b/dendrocat/aperture.py index 1e6351e..d5e21b0 100644 --- a/dendrocat/aperture.py +++ b/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): diff --git a/dendrocat/radiosource.py b/dendrocat/radiosource.py index f898e01..20aa19e 100644 --- a/dendrocat/radiosource.py +++ b/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,16 +394,17 @@ 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' @@ -405,6 +412,8 @@ def get_pixels(self, aperture, catalog=None, data=None, cutouts=None, '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