Merge pull request #11 from madelyn914/unit-fix

Debug aij photometry

stellarphot/photometry.py

@@ -517,13 +517,24 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
                     exposure=0,
                     include_digitization=False):
 
-    if annulus_area == 0:
+    try:
+        no_annulus = (annulus_area == 0).all()
+    except AttributeError:
+        no_annulus = annulus_area == 0
+
+    if no_annulus:
         area_ratio = aperture_area
     else:
         area_ratio = aperture_area * (1 + aperture_area / annulus_area)
 
+    # Units in here are an absolute mess
     poisson_source = gain * flux
 
+    try:
+        poisson_source = poisson_source.value
+    except AttributeError:
+        pass
+
     sky = area_ratio * gain * sky_per_pix
     dark = area_ratio * dark_current_per_sec * exposure
     rn_error = area_ratio * read_noise ** 2

stellarphot/source_detection.py

@@ -5,17 +5,27 @@
 from photutils import fit_2dgaussian
 from astropy.nddata.utils import Cutout2D
 from astropy.stats import gaussian_sigma_to_fwhm
+from astropy import units as u
 
 __all__ = ['source_detection', 'compute_fwhm']
 
 
 def compute_fwhm(ccd, sources, fwhm_estimate=5,
-                 x_column='xcentroid', y_column='ycentroid'):
+                 x_column='xcenter', y_column='ycenter'):
     fwhm_x = []
     fwhm_y = []
     for source in sources:
         x = source[x_column]
         y = source[y_column]
+
+        # Cutout2D needs no units on the center position, so remove unit
+        # if it is present.
+        try:
+            x = x.value
+            y = y.value
+        except AttributeError:
+            pass
+
         cutout = Cutout2D(ccd, (x, y), 5 * fwhm_estimate)
         fit = fit_2dgaussian(cutout.data)
         fwhm_x.append(gaussian_sigma_to_fwhm * fit.x_stddev)
@@ -62,7 +72,9 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
     mean, median, std = sigma_clipped_stats(ccd, sigma=sigma, maxiters=iters)
     daofind = DAOStarFinder(fwhm=fwhm, threshold=threshold * std)
     sources = daofind(ccd - median)
+    print(sources)
     if find_fwhm:
-        x, y = compute_fwhm(ccd, sources, fwhm_estimate=fwhm)
+        x, y = compute_fwhm(ccd, sources, fwhm_estimate=fwhm,
+                            x_column='xcentroid', y_column='ycentroid')
         sources['FWHM'] = (x + y) / 2
     return sources

stellarphot/tests/test_detection.py

@@ -4,8 +4,11 @@
 from astropy.utils.data import get_pkg_data_filename
 from astropy.table import Table
 from astropy.stats import gaussian_sigma_to_fwhm
+from astropy import units as u
+
 from photutils.datasets import make_gaussian_sources_image, make_noise_image
-from stellarphot.source_detection import source_detection
+
+from stellarphot.source_detection import source_detection, compute_fwhm
 from stellarphot.photometry import photutils_stellar_photometry
 
 
@@ -31,6 +34,23 @@ def image(self):
         return self._stars + self._noise
 
 
+@pytest.mark.parametrize('units', [u.pixel, None])
+def test_compute_fwhm(units):
+    fake_image = FakeImage()
+    sources = fake_image.sources
+    if units is not None:
+        # It turns out having a unit on a column is not the same as
+        # things in the column having units. The construct below ensures
+        # that the source table values have units.
+        # Do not try: sources['x_mean'] = sources['x_mean'] * units
+        # Turns out individual values do NOT have units in that case.
+        sources['x_mean'] = [v * units for v in sources['x_mean']]
+        sources['y_mean'] = [v * units for v in sources['y_mean']]
+
+    fwhm_x, fwhm_y = compute_fwhm(fake_image.image, sources,
+                                  x_column='x_mean', y_column='y_mean')
+
+
 def test_detect_source_number_location():
     """
     Make sure we detect the sources in the input table....
@@ -40,12 +60,10 @@ def test_detect_source_number_location():
     found_sources = source_detection(fake_image.image,
                                      fwhm=2 * sources['x_stddev'].mean(),
                                      threshold=10)
-    print(found_sources.colnames)
     # Sort by flux so we can reliably match them
     sources.sort('amplitude')
     found_sources.sort('flux')
 
-    print(found_sources)
     # Do we have the right number of sources?
     assert len(sources) == len(found_sources)
 
@@ -74,8 +92,7 @@ def test_aperture_photometry_no_outlier_rejection():
     phot.sort('aperture_sum')
     sources.sort('amplitude')
     found_sources.sort('flux')
-    print(found_sources['xcentroid'])
-    print(phot['background_per_pixel', 'aperture_sum_err', 'net_flux', 'aperture_sum'])
+
     for inp, out in zip(sources, phot):
         stdev = inp['x_stddev']
         expected_flux = (inp['amplitude'] * 2 * np.pi *
@@ -131,8 +148,7 @@ def test_aperture_photometry_with_outlier_rejection(reject):
     phot.sort('aperture_sum')
     sources.sort('amplitude')
     found_sources.sort('flux')
-    print(found_sources['xcentroid'])
-    print(phot['background_per_pixel', 'aperture_sum_err', 'net_flux', 'aperture_sum'])
+
     for inp, out in zip(sources, phot):
         stdev = inp['x_stddev']
         expected_flux = (inp['amplitude'] * 2 * np.pi *