Repository: feder-observatory/stellarphot

Issue:
Function `find_too_close()` never used
Function `find_too_close(star_locs, aperture_rad, pixel_scale=0.563)` in `photometry.py` doesn't appear to be called from anywhere in the stellarphot package.  Purge or do we need to use it!

Discussion:

Patch:
diff --git a/stellarphot/core.py b/stellarphot/core.py
index 525b0cb..4515fd7 100644
--- a/stellarphot/core.py
+++ b/stellarphot/core.py
@@ -1,11 +1,9 @@
+import numpy as np
 from astropy import units as u
-from astropy.table import QTable, Table, Column
 from astropy.coordinates import EarthLocation, SkyCoord
+from astropy.table import Column, QTable, Table
 from astropy.time import Time
 
-
-import numpy as np
-
 __all__ = ['Camera', 'BaseEnhancedTable', 'PhotometryData', 'CatalogData',
            'SourceListData']
 
@@ -106,6 +104,10 @@ class Camera:
     def __copy__(self):
         return self.copy()
 
+    def __repr__(self):
+        return f"Camera(gain={self.gain}, read_noise={self.read_noise}, " \
+                f"dark_current={self.dark_current}, pixel_scale={self.pixel_scale})"
+
 
 class BaseEnhancedTable(QTable):
     """
@@ -256,7 +258,9 @@ class PhotometryData(BaseEnhancedTable):
 
     input_data: `astropy.table.QTable`, optional (Default: None)
         A table containing all the instrumental aperture photometry results
-        to be validated.
+        to be validated.  Note: It is allowed for the 'ra' and 'dec' columns
+        to have np.nan values, but if they do, the 'bjd' column will not be
+        computed and will also be left with 'np.nan'.
 
     observatory: `astropy.coordinates.EarthLocation`, optional (Default: None)
         The location of the observatory.
@@ -293,8 +297,10 @@ class PhotometryData(BaseEnhancedTable):
     if you do not, an empty table will be returned.
 
     To be accepted as valid, the  `input_data` must MUST contain the following columns
-    with the following units.  The 'consistent count units' simply means it can be any
-    unit, but it must be the same for all the 'consistent count units' columns.
+    with the following units.  The data in those columns is NOT validated, the values in
+    those columns could be invalid.  Furthermore, the 'consistent count units' below
+    simply means it can be any unit, but it must be the same for all the columns with
+    'consistent count units'.
 
     name                  unit
     -----------------     -------
@@ -307,15 +313,15 @@ class PhotometryData(BaseEnhancedTable):
     fwhm_y                u.pix
     width                 u.pix
     aperture              u.pix
-    aperture_area         u.pix * u.pix
+    aperture_area         u.pix
     annulus_inner         u.pix
     annulus_outer         u.pix
-    annulus_area          u.pix * u.pix
+    annulus_area          u.pix
     aperture_sum          consistent count units
     annulus_sum           consistent count units
-    sky_per_pix_avg       consistent count units
-    sky_per_pix_med       consistent count units
-    sky_per_pix_std       consistent count units
+    sky_per_pix_avg       consistent count units (per pixel)
+    sky_per_pix_med       consistent count units (per pixel)
+    sky_per_pix_std       consistent count units (per pixel)
     aperture_net_cnts     consistent count units
     noise_cnts            consistent count units
     noise_electrons       u.electron
@@ -328,10 +334,10 @@ class PhotometryData(BaseEnhancedTable):
     passband              None
     file                  None
 
-    In addition to these required columns, the following columns are computed based
+    In addition to these required columns, the following columns are created based
     on the input data during creation.
 
-    bjd
+    bjd    (only if ra and dec are all real numbers, otherwise set to np.nan)
     night
 
     If these computed columns already exist in `data` class the class
@@ -352,10 +358,10 @@ class PhotometryData(BaseEnhancedTable):
         'fwhm_y' : u.pix,
         'width' : u.pix,
         'aperture' : u.pix,
-        'aperture_area' : u.pix * u.pix,
+        'aperture_area' : u.pix,
         'annulus_inner' : u.pix,
         'annulus_outer' : u.pix,
-        'annulus_area' : u.pix * u.pix,
+        'annulus_area' : u.pix,
         'aperture_sum' : None,
         'annulus_sum' : None,
         'sky_per_pix_avg' : None,
@@ -403,7 +409,8 @@ class PhotometryData(BaseEnhancedTable):
                                  "astropy.time.Time entries.")
 
             # Convert input data to QTable (while also checking for required columns)
-            super().__init__(input_data=input_data, table_description=self.phot_descript,
+            super().__init__(input_data=input_data,
+                             table_description=self.phot_descript,
                              colname_map=colname_map, **kwargs)
 
             # Add the TableAttributes directly to meta (and adding attribute
@@ -419,9 +426,10 @@ class PhotometryData(BaseEnhancedTable):
             self.meta['pixel_scale'] = camera.pixel_scale
 
             # Check for consistency of counts-related columns
-            counts_columns = ['aperture_sum', 'annulus_sum', 'sky_per_pix_avg',
-                            'sky_per_pix_med', 'sky_per_pix_std', 'aperture_net_cnts',
-                            'noise_cnts']
+            counts_columns = ['aperture_sum', 'annulus_sum', 'aperture_net_cnts',
+                              'noise_cnts']
+            counts_per_pixel_columns = ['sky_per_pix_avg', 'sky_per_pix_med',
+                                            'sky_per_pix_std']
             cnts_unit = self[counts_columns[0]].unit
             for this_col in counts_columns[1:]:
                 if input_data[this_col].unit != cnts_unit:
@@ -429,6 +437,16 @@ class PhotometryData(BaseEnhancedTable):
                                     f"with input_data['{counts_columns[0]}'] (should "
                                     f"be {cnts_unit} but it's "
                                     f"{input_data[this_col].unit}).")
+            for this_col in counts_per_pixel_columns:
+                if cnts_unit is None:
+                    perpixel = u.pixel**-1
+                else:
+                    perpixel = cnts_unit * u.pixel**-1
+                if input_data[this_col].unit != perpixel:
+                    raise ValueError(f"input_data['{this_col}'] has inconsistent units "
+                                    f"with input_data['{counts_columns[0]}'] (should "
+                                    f"be {perpixel} but it's "
+                                    f"{input_data[this_col].unit}).")
 
             # Compute additional columns (not done yet)
             computed_columns = ['bjd', 'night']
@@ -488,18 +506,24 @@ class PhotometryData(BaseEnhancedTable):
         the input observations.  It modifies that table in place.
         """
 
-        # Convert times at start of each observation to TDB (Barycentric Dynamical Time)
-        times = Time(self['date-obs'])
-        times_tdb = times.tdb
-        times_tdb.format='jd' # Switch to JD format
+        if ( np.isnan(self['ra']).any() or np.isnan(self['dec']).any() ):
+            print("WARNING: BJD could not be computed in output PhotometryData object "
+                  "because some RA or Dec values are missing.")
+            return np.full(len(self), np.nan)
+        else:
+            # Convert times at start of each observation to TDB (Barycentric Dynamical
+            # Time)
+            times = Time(self['date-obs'])
+            times_tdb = times.tdb
+            times_tdb.format='jd' # Switch to JD format
 
-        # Compute light travel time corrections
-        ip_peg = SkyCoord(ra=self['ra'], dec=self['dec'], unit='degree')
-        ltt_bary = times.light_travel_time(ip_peg, location=observatory)
-        time_barycenter = times_tdb + ltt_bary
+            # Compute light travel time corrections
+            ip_peg = SkyCoord(ra=self['ra'], dec=self['dec'], unit='degree')
+            ltt_bary = times.light_travel_time(ip_peg, location=observatory)
+            time_barycenter = times_tdb + ltt_bary
 
-        # Return BJD at midpoint of exposure at each location
-        return Time(time_barycenter + self['exposure'] / 2, scale='tdb')
+            # Return BJD at midpoint of exposure at each location
+            return Time(time_barycenter + self['exposure'] / 2, scale='tdb')
 
     @property
     def camera(self):
@@ -516,8 +540,8 @@ class PhotometryData(BaseEnhancedTable):
 class CatalogData(BaseEnhancedTable):
     """
     A class to hold astronomical catalog data while performing validation
-    to confirm the minumum required columns ('id', 'ra', and 'dec') are present
-    and have the correct units.
+    to confirm the minumum required columns ('id', 'ra', 'dec', 'mag', and
+    'passband') are present and have the correct units.
 
     As a convience function, when the user passes in an astropy table to validate,
     the user can also pass in a col_rename dict which can be used to rename columns
@@ -731,9 +755,14 @@ class SourceListData(BaseEnhancedTable):
             ra_dec_present = True
             x_y_present = True
             nosky_pos = ('ra' not in data.colnames or
-                        'dec' not in data.colnames)
+                        'dec' not in data.colnames or
+                        np.isnan(data['ra'].value).all() or
+                        np.isnan(data['dec'].value).all())
             noimg_pos = ('xcenter' not in data.colnames or
-                        'ycenter' not in data.colnames)
+                        'ycenter' not in data.colnames or
+                        np.isnan(data['xcenter'].value).all() or
+                        np.isnan(data['ycenter'].value).all())
+
             if nosky_pos:
                 ra_dec_present = False
             if noimg_pos:
@@ -764,3 +793,17 @@ class SourceListData(BaseEnhancedTable):
     @property
     def has_x_y(self):
         return self.meta['has_x_y']
+
+    def drop_ra_dec(self):
+        # drop sky-based positions from existing SourceListData structure
+        self.meta['has_ra_dec'] = False
+        self['ra'] = Column(data=np.full(len(self), np.nan), name='ra', unit=u.deg)
+        self['dec'] = Column(data=np.full(len(self), np.nan), name='dec', unit=u.deg)
+
+    def drop_x_y(self):
+        # drop image-based positionsfrom existing SourceListData structure
+        self.meta['has_x_y'] = False
+        self['xcenter'] = Column(data=np.full(len(self), np.nan), name='ra',
+                                 unit=u.deg)
+        self['ycenter'] = Column(data=np.full(len(self), np.nan), name='dec', 
+                                 unit=u.deg)
diff --git a/stellarphot/photometry/photometry.py b/stellarphot/photometry/photometry.py
index cc32f7f..0bffcd5 100644
--- a/stellarphot/photometry/photometry.py
+++ b/stellarphot/photometry/photometry.py
@@ -1,138 +1,447 @@
 import bottleneck as bn
 import numpy as np
-from photutils import (CircularAnnulus, aperture_photometry, CircularAperture)
-from photutils.centroids import centroid_sources
-
-from astropy.coordinates import SkyCoord, EarthLocation
-from astropy.table import vstack
 from astropy import units as u
+from astropy.coordinates import EarthLocation, SkyCoord
+from astropy.nddata import CCDData, NoOverlapError
+from astropy.table import Column, vstack
 from astropy.time import Time
-
 from ccdproc import ImageFileCollection
+from photutils.aperture import (CircularAnnulus, CircularAperture,
+                                aperture_photometry)
+from photutils.centroids import centroid_sources
+from scipy.spatial.distance import cdist
 
-from astropy.stats import sigma_clipped_stats
-from astropy.nddata import NoOverlapError
+from stellarphot import Camera, PhotometryData, SourceListData
 
 from .source_detection import compute_fwhm
 
-__all__ = ['photutils_stellar_photometry',
+__all__ = ['single_image_photometry',
            'faster_sigma_clip_stats',
            'find_too_close', 'clipped_sky_per_pix_stats',
            'add_to_photometry_table', 'photometry_on_directory',
            'calculate_noise', 'find_bjd']
 
-
-def photutils_stellar_photometry(ccd_image, sources,
-                                 aperture_radius, inner_annulus,
-                                 outer_annulus, gain=1.0, N_R=0, N_dark_pp=0,
-                                 reject_background_outliers=True):
+# Allowed FITS header keywords for exposure values
+EXPOSURE_KEYWORDS = ["EXPOSURE", "EXPTIME", "TELAPSE", "ELAPTIME", "ONTIME",
+                     "LIVETIME"]
+
+def single_image_photometry(ccd_image, sourcelist, camera, observatory_location,
+                            aperture_radius, inner_annulus, outer_annulus,
+                            shift_tolerance, max_adu, fwhm_estimate,
+                            use_coordinates='pixel',
+                            include_dig_noise=True,
+                            reject_too_close=True,
+                            reject_background_outliers=True,
+                            passband_map=None,
+                            fwhm_by_fit=True, fname=None,
+                            logline="single_image_photometry:"):
     """
-    Perform aperture photometry on an image, with a few options for estimating
-    the local background from an annulus around the aperture.
+    Perform aperture photometry on a single image, with an options for estimating
+    the local background from sigma clipped stats of the counts in an annulus around
+    the aperture.  If the sky positions of the sources are in the sourcelist and not
+    the image positions (x,y), then the function will compute the image positions from
+    the sky positions using the WCS of the image.  It assumes the input sky positions
+    are in decimal degree units and in the ICRS frame.
 
     Parameters
     ----------
-    ccd_image : `ccdproc.CCDData`
-        Image on which to perform aperture photometry.
+    ccd_image : `astropy.nddata.CCDData`
+        Image on which to perform aperture photometry.  It's headers must contain
+        DATE-OBS, exposure time (which can be any of the following: EXPOSURE,
+        EXPTIME, TELAPSE, ELAPTIME, ONTIME, or LIVETIME), and FILTER.  If AIRMASS is
+        available it will be added to `phot_table`.  This image  must also have
+        a WCS header associated with it if you want to use sky positions as inputs.
+
+    sourcelist : `stellarphot.SourceList`
+        Table of extracted sources with positions in terms of pixel coordinates OR
+        RA/Dec coordinates. If both positions provided, pixel coordinates will be used.
+        For RA/Dec coordinates to be used, `ccd_image` must have a valid WCS.
+
+    camera : `stellarphot.Camera`
+        Camera object which has gain, read noise and dark current set.
 
-    sources : `astropy.table.Table`
-        Table of extracted sources. Assumed to be the output of
-        `photutils.daofind()` source extraction function.
+    observatory_location : `astropy.coordinates.EarthLocation`
+        Location of the observatory where the images were taken.  Used for calculating
+        the BJD.
 
-    aperture_radius : float
+    aperture_radius : int
         Radius of aperture(s) in pixels.
 
-    inner_annulus : float
+    inner_annulus : int
         Inner radius of the annulus in pixels.
 
-    outer_annulus : float
+    outer_annulus : int
         Outer radius of the annulus in pixels.
 
-    gain : float
-        Gain of the CCD. In units of electrons per DN.
+    shift_tolerance : float
+        If the x/y position needs to be computed using WCS, then computed
+        x/y positions are refined by afterward centroiding the sources.
+        This reflects the tolerance in pixels for the shift between the
+        the computed positions and the refined positions, in pixels.
+        The expected shift shift should not be more than the FWHM, so a
+        measured FWHM might be a good value to provide here.
 
-    N_R : float
-        Read noise of the CCD in electrons per pixel.
+    max_adu : float
+        Maximum allowed pixel value before a source is considered saturated.
 
-    N_dark_pp : float
-        Number of dark counts per pixel.
+    fwhm_estimate : float
+        Initial estimate of the FWHM in pixels for sources in the image.
+        This is used to determine the size of the box used to fit the FWHM
+        (which is 5 times the FWHM estimate in width).
 
-    reject_background_outliers : bool, optional
+    use_coordinates : str, optional (Default: 'pixel')
+        If ``'pixel'``, use the x/y positions in the sourcelist for
+        performing aperture photometry.  If ``'sky'``, use the ra/dec
+        positions in the sourcelist and the WCS of the `ccd_image` to
+        compute the x/y positions on the image.
+
+    reject_background_outliers : bool, optional (Default: True)
         If ``True``, sigma clip the pixels in the annulus to reject outlying
         pixels (e.g. like stars in the annulus)
 
+    reject_too_close : bool, optional (Default: True)
+        If ``True``, any sources that are closer than twice the aperture radius
+        are rejected.  If ``False``, all sources in field are used.
+
+    include_dig_noise : bool, optional (Default: True)
+        If ``True``, include the digitization noise in the calculation of the
+        noise for each observation.  If ``False``, only the Poisson noise from
+        the source and the sky will be included.
+
+    passband_map: dict, optional (Default: None)
+        A dictionary containing instrumental passband names as keys and
+        AAVSO passband names as values. This is used to rename the passband
+        entries in the output photometry table to be AAVSO standard versus
+        whatever is in the source list.
+
+    fwhm_by_fit : bool, optional (default: True)
+        If ``True``, the FWHM will be calculated by fitting a Gaussian to
+        the star. If ``False``, the FWHM will be calculated by finding the
+        second moments of the light distribution. Default is ``True``.
+
+    fname : str, optional ()
+        Name of the image file on which photometry is being performed.
+
+    logline : str, optional (Default: "single_image_photometry:")
+        String to prepend to all log messages.
+
     Returns
     -------
-    phot_table : `astropy.table.Table`
-        Astropy table with columns for flux, x/y coordinates of center,
-        RA/dec coordinates of center, sky background per pixel,
-        net flux, aperture and annulus radii used, and flux error.
+
+    phot_table : `stellarphot.PhotometryData`
+        Photometry data for all the locations in which aperture photometry was
+        performed.  This may be a subset of the sources in the sourcelist if
+        locations were too close to the edge of the image or to each other for
+        successful apeture photometry.  If pixel (x/y) positions were used for
+        the photometry, but a valid WCS header was not available for `ccd_image`,
+        the output 'ra', 'dec', and 'bjd' columns will have np.nan values
+
+    dropped_sources : list
+        This of the star_ids of the sources that fell outside the image or were
+        too close together and did not have photometry performed.
+
+    Notes
+    -----
+    The default behavior (set by `use_coordinates`="pixel") for determining
+    the placement of apertures is to use the x/y positions in the `sourcelist`.
+    This is because in most cases, ra/dec positions are derived from the x/y
+    positions and the WCS of the `ccd_image`.  This default reduces unnecessary
+    computation and works if the WCS is not very accurate.
+
+    When attempting to process multiple images of the same field, it makes sense
+    to use the ra/dec positions in the sourcelist and the WCS of the `ccd_image`
+    to compute the x/y positions on each image individually. In this scenario,
+    the `use_coordinates` parameter should be set to "sky".
     """
 
-    # check that the outer radius is greater or equal the inner radius
-    # for annulus
+    # Check that the input parameters are valid
+    if not isinstance(ccd_image, CCDData):
+        raise TypeError("ccd_image must be a CCDData object, but it is "
+                        f"'{type(ccd_image)}'.")
+    if not isinstance(sourcelist, SourceListData):
+        raise TypeError("sourcelist must be a SourceListData object, but it is "
+                        f"'{type(sourcelist)}'.")
+    if not isinstance(camera, Camera):
+        raise TypeError(f"camera must be a Camera object, but it is '{type(camera)}'.")
+    if not isinstance(observatory_location, EarthLocation):
+        raise TypeError("observatory_location must be a EarthLocation object, but it "
+                        f"is '{type(observatory_location)}'.")
     if inner_annulus >= outer_annulus:
-        raise ValueError("outer_annulus must be greater than inner_annulus")
-
-    # check that the annulus inner radius is greater or equal
-    # the aperture radius
+        raise ValueError(f"outer_annulus ({outer_annulus}) must be greater than "
+                         f"inner_annulus ({inner_annulus}).")
     if aperture_radius >= inner_annulus:
-        raise ValueError("inner_radius must be greater than aperture_radius")
+        raise ValueError(f"aperture_radius ({aperture_radius}) must be greater than "
+                         f"inner_annulus ({inner_annulus}).")
+    if (shift_tolerance<=0):
+        raise ValueError(f"shift_tolerance ({shift_tolerance}) must be greater than 0 "
+                         "(should be on order of FWHM).")
+    if (max_adu<=0):
+        raise ValueError(f"max_adu ({max_adu}) must be greater than 0.")
+    if (use_coordinates not in ['pixel', 'sky']):
+        raise ValueError(f"input_coordinates ({use_coordinates}) must be either "
+                         "'pixel' or 'sky'.")
+
+    # # Set high pixels to NaN
+    # ccd_image.data[ccd_image.data > max_adu] = np.nan
+
+    # Extract necessary values from sourcelist structure
+    star_ids = sourcelist['star_id'].value
+    xs = sourcelist['xcenter'].value
+    ys = sourcelist['ycenter'].value
+    ra = sourcelist['ra'].value
+    dec = sourcelist['dec'].value
+    src_cnt = len(sourcelist)
+
+    # If RA/Dec are available attempt to use them to determine the source positions
+    if use_coordinates == 'sky' and sourcelist.has_ra_dec:
+        try:
+            imgpos = ccd_image.wcs.world_to_pixel(SkyCoord(ra, dec, unit=u.deg,
+                                                           frame='icrs'))
+            xs, ys = imgpos[0], imgpos[1]
+        except AttributeError:
+            # No WCS
+            raise ValueError("ccd_image must have a valid WCS to use RA/Dec!")
+    elif use_coordinates == 'sky' and not sourcelist.has_ra_dec:
+        raise ValueError("use_coordinates='sky' but sourcelist does not have"
+                         "RA/Dec coordinates!")
+
+    # Reject sources that are within an aperture diameter of each other.
+    dropped_sources = []
+    try:
+        too_close = find_too_close(sourcelist, aperture_radius,
+                                pixel_scale=camera.pixel_scale.value)
+    except Exception as e:
+        # Any failure here is BAD, so raise an error
+        raise RuntimeError(
+            f"Call to find_too_close() returned {type(e).__name__}: {str(e)}")
+    too_close_cnt = np.sum(too_close)
+    non_overlap = ~too_close
+    print(f"{logline} {too_close_cnt} of {src_cnt} sources within 2 aperture radii of "
+          "nearest neighbor", end="")
+    if reject_too_close:
+        # Track dropped sources due to being too close together
+        dropped_sources.extend(star_ids[too_close].tolist())
+        # Remove sources too close together
+        star_ids = star_ids[non_overlap]
+        xs = xs[non_overlap]
+        ys = ys[non_overlap]
+        ra = ra[non_overlap]
+        dec = dec[non_overlap]
+        print(" ... removed them.")
+    else:
+        print(" ... keeping them.")
+
+    # Remove all source positions too close to edges of image (where the annulus would
+    # extend beyond the image boundaries).
+    padding = outer_annulus
+    out_of_bounds = ( (xs < padding) | (xs > (ccd_image.shape[1] - padding)) |
+                    (ys < padding) | (ys  > (ccd_image.shape[0] - padding)) )
+    in_bounds = ~out_of_bounds
+    # Track dropped sources due to out of bounds positions
+    dropped_sources.extend(star_ids[out_of_bounds].tolist())
+    # Remove sources too close to the edges
+    star_ids = star_ids[in_bounds]
+    xs = xs[in_bounds]
+    ys = ys[in_bounds]
+    ra = ra[in_bounds]
+    dec = dec[in_bounds]
+    in_cnt = np.sum(in_bounds)
+    out_cnt = np.sum(out_of_bounds)
+    print(f"{logline} {in_cnt} sources kept of {src_cnt} original sources ({out_cnt} "
+          "sources too close to image edge).")
+
+    # If we are using x/y positions previously obtained from the ra/dec positions and
+    # WCS, then recentroid the sources to refine the positions. This is
+    # particularly useful is processing multiple images of the same field
+    # and just passing the same sourcelist when calling single_image_photometry
+    # on each image.
+    if use_coordinates == 'sky':
+        try:
+            xcen, ycen = centroid_sources(ccd_image.data, xs, ys,
+                                        box_size=2 * aperture_radius + 1)
+        except NoOverlapError:
+            raise NoOverlapError("Determining new centroid failed!")
+        except Exception as e:
+            raise RuntimeError(
+                f"Call to centroid_sources() returned {type(e).__name__}: {str(e)}")
+        else:
+            # Calculate offset between centroid in this image and the positions
+            # based on input RA/Dec.
+            center_diff = np.sqrt((xs - xcen)**2 + (ys - ycen)**2)
+
+            # The center really shouldn't move more than about the fwhm, could
+            # rework this in the future to use that instead.
+            too_much_shift = center_diff > shift_tolerance
+
+            # If the shift is too large, use the WCS-derived positions instead
+            # (these sources are probably too faint for centroiding to work well)
+            xcen[too_much_shift] = xs[too_much_shift]
+            ycen[too_much_shift] = ys[too_much_shift]
+            xs, ys = xcen, ycen
+
+    # Compute RA/Dec if not already provided
+    if not sourcelist.has_ra_dec:
+        try:
+            skypos = ccd_image.wcs.pixel_to_world(xs, ys)
+            ra = skypos.ra.value
+            dec = skypos.dec.value
+        except AttributeError:
+            ra = [np.nan] * len(xs)
+            dec = [np.nan] * len(ys)
 
-    # Extract x,y coordinates from sources table, construct aperture and
-    # annulus objects from coordinates, and perform aperture photometry
-    coords = np.array([sources['xcentroid'], sources['ycentroid']]).T
-    apertures = CircularAperture(coords, aperture_radius)
-    annulus = CircularAnnulus(coords, inner_annulus, outer_annulus)
-    phot_table = aperture_photometry(ccd_image, apertures)
-    phot_table_1 = aperture_photometry(ccd_image, annulus)
+    # Define apertures and annuli for the aperture photometry
+    aper_locs = np.array([xs, ys]).T
+    apers = CircularAperture(aper_locs, r=aperture_radius)
+    anuls = CircularAnnulus(aper_locs, r_in=inner_annulus, r_out=outer_annulus)
 
-    # Obtain the local background/pixel and net flux between the aperture and
-    # annulus objects
-    n_pix_ap = apertures.area
-    n_pix_ann = annulus.area
+    # Perform the aperture photometry
+    photom = aperture_photometry(ccd_image.data, (apers, anuls),
+                                mask=ccd_image.mask, method='center')
 
-    if reject_background_outliers:
-        annulus_masks = annulus.to_mask()
-        bkgd_pp = []
-        for mask in annulus_masks:
-            annulus_data = mask.cutout(ccd_image)
-            bool_mask = mask.data < 1
-            # Only include whole pixels in the estimate of the background
-            bkgd, _, _ = sigma_clipped_stats(annulus_data * mask,
-                                             mask=bool_mask)
-            bkgd_pp.append(bkgd)
-        bkgd_pp = np.array(bkgd_pp)
+    # Add source ids to the photometry table
+    photom['star_id'] = star_ids
+    photom['ra'] = ra * u.deg
+    photom['dec'] = dec * u.deg
+
+    # Drop ID column from aperture_photometry()
+    del photom['id']
+
+    # Add various CCD image parameters to the photometry table
+    if fname is not None:
+        photom['file'] = fname
     else:
-        bkgd_pp = phot_table_1['aperture_sum'] / n_pix_ann
-
-    net_flux = phot_table['aperture_sum'] - (n_pix_ap * bkgd_pp)
-    phot_table['background_per_pixel'] = bkgd_pp
-    phot_table['net_flux'] = net_flux
-
-    # Return a columns with the aperture radius and
-    # the inner/outer annulus radii
-    phot_table['aperture_radius'] = \
-        np.ones(len(phot_table['aperture_sum'])) * aperture_radius
-    phot_table['inner_radius'] = \
-        np.ones(len(phot_table['aperture_sum'])) * inner_annulus
-    phot_table['outer_radius'] = \
-        np.ones(len(phot_table['aperture_sum'])) * outer_annulus
-
-    # Obtain RA/Dec coordinates and add them to table
+        photom['file'] = [''] * len(photom)
+
+    # Search for exposure keyword in header and set exposure if found
+    matched_kw = None
+    for kw in EXPOSURE_KEYWORDS:
+        if kw in ccd_image.header:
+            matched_kw = kw
+            break
+
+    if matched_kw is None:
+        raise ValueError(f"None of the allowed exposure keywords "
+                         f"({format(', '.join(EXPOSURE_KEYWORDS))}) found in the header!")
+    else:
+        photom['exposure'] = [ccd_image.header[matched_kw]] * len(photom) * u.second
+
     try:
-        ra, dec = ccd_image.wcs.all_pix2world(coords[0], coords[1], 0)
-        phot_table['RA_center'] = ra
-        phot_table['Dec_center'] = dec
-    except AttributeError:
-        pass
+        photom['date-obs'] = Column(data=Time([ccd_image.header['DATE-OBS']]
+                                              * len(photom),
+                                    format='isot', scale='utc'), name='date-obs')
+    except KeyError:
+        raise ValueError("'DATE-OBS' not found in CCD image header!")
+
+    # Check for airmass keyword in header and set 'airmass' if found,
+    # but accept it may not be available
+    try:
+        photom['airmass'] = [ccd_image.header['AIRMASS']] * len(photom)
+    except KeyError:
+        print(f"{logline} 'AIRMASS' not found in CCD "
+                            "image header ... setting to NaN!")
+        photom['airmass'] = [np.nan] * len(photom)
+
+    # Check for filter keyword in header and set 'filter' if found
+    try:
+        photom['filter'] = [ccd_image.header['FILTER']] * len(photom)
+    except KeyError:
+        raise ValueError("'FILTER' not found in CCD image header!")
+    photom.rename_column('filter', 'passband')
+
+    # Save aperture and annulus information
+    photom.rename_column('aperture_sum_0', 'aperture_sum')
+    photom.rename_column('aperture_sum_1', 'annulus_sum')
+    photom['aperture_sum'].unit = ccd_image.unit
+    photom['annulus_sum'].unit = ccd_image.unit
+    photom['aperture'] = apers.r * u.pixel
+    photom['annulus_inner'] = anuls.r_in * u.pixel
+    photom['annulus_outer'] = anuls.r_out * u.pixel
+    # By convention, area is in units of pixels (not pixels squared) in a digital image
+    photom['aperture_area'] = apers.area * u.pixel
+    photom['annulus_area'] = anuls.area  * u.pixel
+
+    if reject_background_outliers:
+        print(f"{logline} computing clipped sky stats ...")
+        try:
+            avg_sky_per_pix, med_sky_per_pix, std_sky_per_pix = \
+                    clipped_sky_per_pix_stats(ccd_image, anuls)
+        except AttributeError:
+            print(f"{logline} BAD ANNULUS, stats set to np.nan!")
+            avg_sky_per_pix, med_sky_per_pix, std_sky_per_pix = \
+                np.nan, np.nan, np.nan
+        photom['sky_per_pix_avg'] = avg_sky_per_pix / u.pixel
+        photom['sky_per_pix_med'] = med_sky_per_pix / u.pixel
+        photom['sky_per_pix_std'] = std_sky_per_pix / u.pixel
+        print("                          ...DONE computing clipped sky stats")
+    else: # Don't reject outliers (but why would you do this?)
+        print(f"{logline} WARNING: Computing sky per pixel without clipping (set "
+              "reject_background_outliers=True to clip).")
+        med_pp = []
+        std_pp = []
+        for mask in anuls.to_mask():
+            annulus_data = mask.cutout(ccd_image)
+            med_pp.append(np.median(annulus_data))
+            std_pp.append(np.std(annulus_data))
+        photom['sky_per_pix_avg'] = photom['annulus_sum'] / photom['annulus_area']
+        photom['sky_per_pix_med'] = np.array(med_pp) * ccd_image.unit / u.pixel
+        photom['sky_per_pix_std'] = np.array(std_pp) * ccd_image.unit / u.pixel
+
+    # Compute counts using clipped stats on sky per pixel
+    photom['aperture_net_cnts'] = (photom['aperture_sum'].value -
+                                    (photom['aperture_area'].value *
+                                    photom['sky_per_pix_avg'].value))
+    photom['aperture_net_cnts'].unit = ccd_image.unit
+
+    # Fit the FWHM of the sources
+    print(f"{logline} fitting fwhm of all sources ...")
+    fwhm_x, fwhm_y = compute_fwhm(ccd_image, photom,
+                                  fwhm_estimate=fwhm_estimate, fit=fwhm_by_fit)
+    print("                          ...DONE fitting fwhm of all sources")
+
+    # Deal with bad FWHM values
+    bad_fwhm = (fwhm_x < 1) | (fwhm_y < 1) # Set bad values to NaN now
+    fwhm_x[bad_fwhm] = np.nan
+    fwhm_y[bad_fwhm] = np.nan
+    photom['fwhm_x'] = fwhm_x * u.pixel
+    photom['fwhm_y'] = fwhm_y * u.pixel
+    photom['width'] = ((fwhm_x + fwhm_y) / 2) * u.pixel
 
-    # Obtain flux error and add column to return table
-    noise = np.sqrt(gain * net_flux + n_pix_ap * (1 + (n_pix_ap / n_pix_ann)) *
-                    (gain * (bkgd_pp + N_dark_pp) + (N_R**2)))
-    phot_table['aperture_sum_err'] = noise
+    # Flag sources with bad counts before computing noise.
+    # This can happen, for example, when the object is faint and centroiding is
+    # bad.  It can also happen when the sky background is low.
+    bad_cnts = photom['aperture_net_cnts'].value < 0
+    all_bads = bad_cnts | bad_fwhm
+    photom['aperture_net_cnts'][all_bads] = np.nan
 
-    return phot_table
+    # Compute and save noise
+    noise = calculate_noise(camera=camera,
+                            counts=photom['aperture_net_cnts'].value,
+                            sky_per_pix=photom['sky_per_pix_avg'].value,
+                            aperture_area=photom['aperture_area'].value,
+                            annulus_area=photom['annulus_area'].value,
+                            exposure=photom['exposure'].value,
+                            include_digitization=include_dig_noise)
+    photom['noise_electrons'] = noise  # Noise in electrons
+    photom['noise_electrons'].unit = u.electron
+    photom['noise_cnts'] = noise / camera.gain.value  # Noise in counts
+    photom['noise_cnts'].unit = ccd_image.unit
+
+    # Compute and save SNR
+    snr = camera.gain.value * photom['aperture_net_cnts'] / noise
+    photom['snr'] = snr
+    photom['mag_error'] = 1.085736205 / snr
+
+    # Compute instrumental magnitudes
+    photom['mag_inst'] = (
+        -2.5 * np.log10(camera.gain.value * photom['aperture_net_cnts'].value /
+                                photom['exposure'].value)
+    )
+
+    # Create PhotometryData object to return
+    photom_data = PhotometryData(observatory=observatory_location, camera=camera,
+                                input_data=photom, passband_map=passband_map)
+
+    return photom_data, dropped_sources
 
 
 def faster_sigma_clip_stats(data, sigma=5, iters=5, axis=None):
@@ -183,23 +492,27 @@ def faster_sigma_clip_stats(data, sigma=5, iters=5, axis=None):
             bn.nanstd(data, axis=axis))
 
 
-def find_too_close(star_locs, aperture_rad, pixel_scale=0.563):
+def find_too_close(sourcelist, aperture_rad, pixel_scale=None):
     """
     Identify sources that are closer together than twice the aperture radius.
 
+    If 'xcenter' and 'ycenter' are available in the sourcelist (as determined by
+    the value of sourcelist.has_x_y), they are used to determine separation of sources
+    in units of pixels, otherwise if only ra/dec are available, the pixel_scale
+    is necessary to determine the separation.
+
     Parameters
     ----------
 
-    star_locs : tuple of numpy array
-        The first entry in the tuple should be the right ascension of the
-        sources, in degrees. The second should be the declination of
-        the sources, in degrees.
+    sourcelist : `stellarphot.SourceListData`
+        A list of sources with x/y and/or RA/Dec coordinates.
 
     aperture_rad : int
         Radius of the aperture, in pixels.
 
-    pixel_scale : float, optional
-        Pixel scale of the image in arcsec/pixel
+    pixel_scale : float, optional (Default: None)
+        Pixel scale of the image in arcsec/pixel. Only required
+        if x/y coordinates are NOT provided.
 
     Returns
     -------
@@ -208,11 +521,26 @@ def find_too_close(star_locs, aperture_rad, pixel_scale=0.563):
         Array the same length as the RA/Dec that is ``True`` where the sources
         are closer than two aperture radii, ``False`` otherwise.
     """
-    star_coords = SkyCoord(ra=star_locs[0], dec=star_locs[1],
+    if not isinstance(sourcelist, SourceListData):
+        raise TypeError("sourcelist must be of type SourceListData not "
+                        f"{type(sourcelist)}")
+
+    if sourcelist.has_x_y:
+        x, y = sourcelist['xcenter'], sourcelist['ycenter']
+        # Find the pixel distance to the nearest neighbor for each source
+        dist_mat = cdist(np.array([x, y]).T, np.array([x, y]).T, metric='euclidean')
+        np.fill_diagonal(dist_mat, np.inf)
+        # Return array with True where the distance is less than twice the aperture
+        # radius
+        return (dist_mat.min(0) < 2 * aperture_rad)
+    elif sourcelist.has_ra_dec:
+        star_coords = SkyCoord(ra=sourcelist['ra'], dec=sourcelist['dec'],
                            frame='icrs', unit='degree')
-    idxc, d2d, d3d = star_coords.match_to_catalog_sky(star_coords,
-                                                      nthneighbor=2)
-    return d2d < (aperture_rad * 2 * 0.563 * u.arcsec)
+        idxc, d2d, d3d = star_coords.match_to_catalog_sky(star_coords,
+                                                        nthneighbor=2)
+        return (d2d < (aperture_rad * 2 * pixel_scale * u.arcsec))
+    else:
+        raise ValueError("sourcelist must have x/y or ra/dec coordinates")
 
 
 def clipped_sky_per_pix_stats(data, annulus, sigma=5, iters=5):
@@ -563,8 +891,7 @@ def photometry_on_directory(directory_with_images, object_of_interest,
     gain = camera.gain
 
     # Compute and save noise
-    noise = calculate_noise(gain=camera.gain, read_noise=camera.read_noise,
-                            dark_current_per_sec=camera.dark_current,
+    noise = calculate_noise(camera=camera,
                             counts=all_phot['aperture_net_counts'],
                             sky_per_pix=all_phot['sky_per_pix_avg'].value,
                             aperture_area=all_phot['aperture_area'],
@@ -582,17 +909,17 @@ def photometry_on_directory(directory_with_images, object_of_interest,
     return all_phot
 
 
-def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
-                    counts=0.0, sky_per_pix=0.0, aperture_area=0, annulus_area=0,
+def calculate_noise(camera=None, counts=0.0, sky_per_pix=0.0,
+                    aperture_area=0, annulus_area=0,
                     exposure=0, include_digitization=False):
     """
     Computes the noise in a photometric measurement.
 
     This function computes the noise (in units of electrons) in a photometric
     measurement using the revised CCD equation from Collins et al (2017) AJ, 153, 77
-    who based their expression on the one originally proposed for SNR by 
+    who based their expression on the one originally proposed for SNR by
     Merline, W. J., & Howell, S. B. 1995, Experimental Astronomy, 6, 163.
-    
+
     The equation is:
 
     .. math::
@@ -600,12 +927,12 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
         \\sigma = \\sqrt{G \\cdot N_C + A \\cdot \\left(1 + \\frac{A}{B}\\right)\\cdot \\left[ G\\cdot S + D \\cdot t + R^2 + (0.289 G)^2\\right]}
 
     where :math:`\sigma` is the noise, :math:`G` is the gain,
-    :math:`N_C` is the source counts (which is photon/electron counts divided by gain), 
+    :math:`N_C` is the source counts (which is photon/electron counts divided by gain),
     :math:`A` is the aperture area in pixels,
     :math:`B` is the annulus area in pixels,
     :math:`S` is the sky counts per pixel,
     :math:`D` is the dark current in electrons per second,
-    :math:`R` is the read noise in electrons per pixel per read, 
+    :math:`R` is the read noise in electrons per pixel per read,
     and :math:`t` is exposure time in seconds.
 
     Note: The :math:`(0.289 G)^2` term is "digitization noise" and is optional.
@@ -613,14 +940,9 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
     Parameters
     ----------
 
-    gain : float, optional
-        Gain of the CCD. In units of electrons per count.
-
-    read_noise : float, optional
-        Read noise of the CCD in electrons.
-
-    dark_current_per_sec : float, optional
-        Dark current of the CCD in electrons per second.
+    camera : `stellarphot.Camera`
+        The camera object that contains the gain, read noise and dark current
+        of the CCD.
 
     counts : float, optional
         Counts of the source.
@@ -646,6 +968,15 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
     noise : float
         The noise in the photometric measurement in electrons.
     """
+    if camera is None:
+        raise ValueError("camera must be provided")
+    elif not isinstance(camera, Camera):
+        raise ValueError(f"camera must be of type Camera not '{type(camera)}'")
+
+    # Extract values from camera object
+    gain = camera.gain.value
+    dark_current_per_sec = camera.dark_current.value
+    read_noise = camera.read_noise.value
 
     try:
         no_annulus = (annulus_area == 0).all()
diff --git a/stellarphot/photometry/source_detection.py b/stellarphot/photometry/source_detection.py
index dc84386..b4d096d 100644
--- a/stellarphot/photometry/source_detection.py
+++ b/stellarphot/photometry/source_detection.py
@@ -1,10 +1,11 @@
 import numpy as np
 
-from astropy.stats import sigma_clipped_stats
+from astropy import units as u
 from astropy.modeling.models import Const2D, Gaussian2D
 from astropy.modeling.fitting import LevMarLSQFitter
+from astropy.nddata import CCDData
+from astropy.stats import sigma_clipped_stats
 from astropy.table import Table
-from astropy import units as u
 
 from photutils.detection import DAOStarFinder
 from photutils.morphology import data_properties
@@ -155,8 +156,9 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
     Parameters
     ----------
 
-    ccd : `astropy.nddata.CCDData`
-        The CCD Image array.
+    ccd : `numpy.ndarray` or `astropy.nddata.CCDData`
+        This is the 2-D array of the image to be analyzed.  If
+        a CCDData object is passed, its data attribute will be used.
 
     fwhm : float, optional (default=8)
         Initial estimate of full-width half-max of stars in the image,
@@ -195,6 +197,9 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
         the input CCDData object has WCS information, the columns `ra` and
         `dec` are also populated, otherwise they are set to NaN.
     """
+    # if image is a CCDData object, extract the data array
+    if not isinstance(ccd, CCDData) and not isinstance(ccd, np.ndarray):
+        raise ValueError("ccd must be a numpy array or CCDData object")
 
     # If user uses units for fwhm or sky_per_pix_avg, convert to value
     if isinstance(sky_per_pix_avg, u.Quantity):
@@ -206,8 +211,11 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
     # if not provided).  Using clipped stats should hopefully get rid of any
     # bright stars that might be in the image, so the mean should be a good
     # estimate of the sky background.
+    print("source_detection: You may see a warning about invalid values in the "
+          "input image.  This is expected if any pixels are saturated and can be "
+          "ignored.")
     mean, median, std = sigma_clipped_stats(ccd, sigma=sigma, maxiters=iters)
-    print(f"source_detection: mean={mean:.4f}, median={median:.4f}, std={std:.4f}")
+    print(f"source_detection: sigma_clipped_stats mean={mean:.4f}, median={median:.4f}, std={std:.4f}")
     if sky_per_pix_avg is None:
         sky_per_pix_avg = mean
         print(f"source_detection: sky_per_pix_avg set to {sky_per_pix_avg:.4f}")
@@ -220,7 +228,10 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
     # daofind should be run on background subtracted image
     # (fails, or at least returns garbage, if sky_per_pix_avg is too low)
     daofind = DAOStarFinder(fwhm = fwhm, threshold = threshold * std)
-    sources = daofind(ccd - sky_per_pix_avg)
+    if isinstance(ccd, CCDData):
+        sources = daofind(ccd.data - sky_per_pix_avg)
+    else:
+        sources = daofind(ccd - sky_per_pix_avg)
 
     # Identify sources near the edge of the image and remove them
     # from the source list.


Test patch:
diff --git a/stellarphot/photometry/tests/data/test_corner.ecsv b/stellarphot/photometry/tests/data/test_corner.ecsv
new file mode 100644
index 0000000..5dad888
--- /dev/null
+++ b/stellarphot/photometry/tests/data/test_corner.ecsv
@@ -0,0 +1,65 @@
+# %ECSV 1.0
+# ---
+# datatype:
+# - {name: star_id, datatype: int64}
+# - {name: ra, unit: deg, datatype: float64}
+# - {name: dec, unit: deg, datatype: float64}
+# - {name: xcenter, unit: pix, datatype: float64}
+# - {name: ycenter, unit: pix, datatype: float64}
+# meta:
+#   __serialized_columns__:
+#     dec:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: &id001 !astropy.units.Unit {unit: deg}
+#       value: !astropy.table.SerializedColumn {name: dec}
+#     ra:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id001
+#       value: !astropy.table.SerializedColumn {name: ra}
+#     xcenter:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: &id002 !astropy.units.Unit {unit: pix}
+#       value: !astropy.table.SerializedColumn {name: xcenter}
+#     ycenter:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id002
+#       value: !astropy.table.SerializedColumn {name: ycenter}
+#   aperture_photometry_args: method='center', subpixels=5
+#   dark_current: !astropy.units.Quantity
+#     unit: !astropy.units.Unit {unit: electron / s}
+#     value: 0.01
+#   gain: !astropy.units.Quantity
+#     unit: !astropy.units.Unit {unit: electron / adu}
+#     value: 1.5
+#   height: !astropy.units.Quantity
+#     unit: !astropy.units.Unit {unit: m}
+#     value: 310.9999999998063
+#   lat: !astropy.coordinates.Latitude
+#     unit: *id001
+#     value: 46.86678
+#   lon: !astropy.coordinates.Longitude
+#     unit: *id001
+#     value: -96.45328
+#     wrap_angle: !astropy.coordinates.Angle
+#       unit: *id001
+#       value: 180.0
+#   name: Aperture photometry results
+#   pixel_scale: !astropy.units.Quantity
+#     unit: !astropy.units.Unit {unit: arcsec / pix}
+#     value: 0.563
+#   read_noise: !astropy.units.Quantity
+#     unit: !astropy.units.Unit {unit: electron}
+#     value: 10.0
+#   version: {Python: 3.11.4, astropy: 5.3.1, bottleneck: 1.3.7, gwcs: 0.18.3, matplotlib: 3.7.2, numpy: 1.25.1, photutils: 1.8.0, scipy: 1.11.1,
+#     skimage: 0.21.0, sklearn: 1.3.0}
+# schema: astropy-2.0
+star_id ra dec xcenter ycenter
+46 78.24000741682637 22.20399499365103 2901.1665129043017 277.5520148018723
+47 78.24439733525055 22.20387400138809 2926.5832753701025 283.0833495082689
+67 78.09260250545755 22.255900218799265 1974.1108195504123 389.221188601664
+69 78.09549601089884 22.256816929145405 1989.3171503776648 399.0300672006914
+91 77.96114692249174 22.30477706450149 1144.0979389304412 505.81064686500184
+92 78.16622063517367 22.258255458823285 2393.335602710066 508.77794245798015
+93 77.96096642411284 22.30587151471993 1141.3812664836696 512.3480118471604
+94 77.9588585863262 22.307161687220702 1127.3004747217678 517.3683506739007
+96 78.16803226534509 22.25985518912216 2401.271927701757 521.2898580662454
diff --git a/stellarphot/photometry/tests/fake_image.py b/stellarphot/photometry/tests/fake_image.py
new file mode 100644
index 0000000..63813c3
--- /dev/null
+++ b/stellarphot/photometry/tests/fake_image.py
@@ -0,0 +1,65 @@
+
+import astropy.io.fits as fits
+from astropy.nddata import CCDData
+from astropy.table import Table
+from astropy.utils.data import get_pkg_data_filename
+
+from photutils.datasets import make_gaussian_sources_image, make_noise_image
+
+
+class FakeImage:
+    """
+    Creates a fake image with a set of sources using the datafile stored
+    at `data/test_sources.csv`.
+
+    Parameters
+    ----------
+
+    noise_dev : float, optional
+        The standard deviation of the noise in the image.
+    """
+    def __init__(self, noise_dev=1.0):
+        self.image_shape = [400, 500]
+        data_file = get_pkg_data_filename('data/test_sources.csv')
+        self._sources = Table.read(data_file)
+        self.mean_noise = self.sources['amplitude'].max() / 100
+        self.noise_dev = noise_dev
+        self._stars = make_gaussian_sources_image(self.image_shape,
+                                                  self.sources)
+        self._noise = make_noise_image(self._stars.shape,
+                                       mean=self.mean_noise,
+                                       stddev=noise_dev)
+        # Sky background per pixel should be the mean level of the noise.
+        self._sources['sky_per_pix_avg'] = self.mean_noise
+
+    @property
+    def sources(self):
+        """
+        Return the table of sources used to create the fake image.
+        """
+        return self._sources
+
+    @property
+    def image(self):
+        """
+        Return the fake image.
+        """
+        return self._stars + self._noise
+
+
+class FakeCCDImage(CCDData):
+    # Generates a fake CCDData object for testing purposes.
+    def __init__(self):
+        base_data = FakeImage()
+        super().__init__(base_data.image.copy(), unit='adu')
+        # Add some additional features to the CCDData object, like
+        # a header and the sources used to create the image.abs
+        self.header = fits.Header()
+        self.header['EXPOSURE'] = 1.0
+        self.header['DATE-OBS'] = '2018-01-01T00:00:00.0'
+        self.header['AIRMASS'] = 1.2
+        self.header['FILTER'] = 'V'
+        self.sources = base_data.sources.copy()
+        self.noise_dev = base_data.noise_dev
+        self.mean_noise = base_data.mean_noise
+        self.image_shape = base_data.image_shape.copy()
diff --git a/stellarphot/photometry/tests/test_detection.py b/stellarphot/photometry/tests/test_detection.py
index 3472cb7..cbb81d7 100644
--- a/stellarphot/photometry/tests/test_detection.py
+++ b/stellarphot/photometry/tests/test_detection.py
@@ -1,56 +1,13 @@
 import numpy as np
 import pytest
 
-from astropy.utils.data import get_pkg_data_filename
-from astropy.table import Table, QTable
+from astropy.table import QTable
 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.photometry import (source_detection, compute_fwhm,
-                                    photutils_stellar_photometry)
-
-
-class FakeImage:
-    """
-    Creates a fake image with a set of sources using the datafile stored
-    at `data/test_sources.csv`.
-
-    Parameters
-    ----------
-
-    noise_dev : float, optional
-        The standard deviation of the noise in the image.
-    """
-    def __init__(self, noise_dev=1.0):
-        self.image_shape = [400, 500]
-        data_file = get_pkg_data_filename('data/test_sources.csv')
-        self._sources = Table.read(data_file)
-        self.mean_noise = self.sources['amplitude'].max() / 100
-        self.noise_dev = noise_dev
-        self._stars = make_gaussian_sources_image(self.image_shape,
-                                                  self.sources)
-        self._noise = make_noise_image(self._stars.shape,
-                                       mean=self.mean_noise,
-                                       stddev=noise_dev)
-        # Sky background per pixel should be the mean level of the noise.
-        self._sources['sky_per_pix_avg'] = self.mean_noise
-
-    @property
-    def sources(self):
-        """
-        Return the table of sources used to create the fake image.
-        """
-        return self._sources
-
-    @property
-    def image(self):
-        """
-        Return the fake image.
-        """
-        return self._stars + self._noise
+from stellarphot.photometry import (source_detection, compute_fwhm)
 
+from fake_image import FakeImage
 
 @pytest.mark.parametrize('units', [u.pixel, None])
 def test_compute_fwhm(units):
@@ -122,8 +79,8 @@ def test_detect_source_with_padding():
 
 
 ##
-## Disabled BOTH FINAL TESTS until photutils_stellar_photometry re-written
-## to access source_detection output properly.
+## Disabled BOTH FINAL TESTS until new single_image_photometry
+## is implemented.  Also moved tests to test_photometry.py
 ##
 # def test_aperture_photometry_no_outlier_rejection():
 #     fake_image = FakeImage()
@@ -148,7 +105,7 @@ def test_detect_source_with_padding():
 #         This expected flux is correct IF there were no noise. With noise, the
 #         standard deviation in the sum of the noise within in the aperture is
 #         n_pix_in_aperture times the single-pixel standard deviation.
-        
+
 #         We could require that the result be within some reasonable
 #         number of those expected variations or we could count up the
 #         actual number of background counts at each of the source
diff --git a/stellarphot/photometry/tests/test_photometry.py b/stellarphot/photometry/tests/test_photometry.py
index b2815da..6f7038b 100644
--- a/stellarphot/photometry/tests/test_photometry.py
+++ b/stellarphot/photometry/tests/test_photometry.py
@@ -1,14 +1,26 @@
-import pytest
 import numpy as np
+import pytest
+from astropy import units as u
+from astropy.coordinates import EarthLocation
+from astropy.io import ascii
+from astropy.utils.data import get_pkg_data_filename
+from fake_image import FakeCCDImage
 
-from stellarphot.photometry import calculate_noise
+from stellarphot.core import Camera, SourceListData
+from stellarphot.photometry import (
+    calculate_noise,
+    find_too_close,
+    single_image_photometry,
+    source_detection,
+)
 
 GAINS = [1.0, 1.5, 2.0]
 
-
 def test_calc_noise_defaults():
-    # If we put in nothing we should get zero back
-    assert calculate_noise() == 0
+    # If we put in nothing we should get an error about is missing camera
+    # instance.
+    with pytest.raises(ValueError):
+        assert calculate_noise() == 0
 
 
 @pytest.mark.parametrize('aperture_area', [5, 20])
@@ -19,7 +31,12 @@ def test_calc_noise_source_only(gain, aperture_area):
     counts = 100
     expected = np.sqrt(gain * counts)
 
-    np.testing.assert_allclose(calculate_noise(gain=gain,
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu, read_noise=0*u.electron,
+                    dark_current=0*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
+    np.testing.assert_allclose(calculate_noise(camera,
                                                counts=counts,
                                                aperture_area=aperture_area),
                                expected)
@@ -32,10 +49,16 @@ def test_calc_noise_dark_only(gain, aperture_area):
     # but this is basically Poisson error.
     dark_current = 10
     exposure = 20
+
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu,
+                    read_noise=0*u.electron,
+                    dark_current=dark_current*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
     expected = np.sqrt(dark_current * aperture_area * exposure)
 
-    np.testing.assert_allclose(calculate_noise(gain=gain,
-                                               dark_current_per_sec=dark_current,
+    np.testing.assert_allclose(calculate_noise(camera,
                                                aperture_area=aperture_area,
                                                exposure=exposure),
                                expected)
@@ -49,8 +72,13 @@ def test_calc_read_noise_only(gain, aperture_area):
     read_noise = 10
     expected = np.sqrt(aperture_area * read_noise**2)
 
-    np.testing.assert_allclose(calculate_noise(gain=gain,
-                                               read_noise=read_noise,
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu,
+                    read_noise=read_noise*u.electron,
+                    dark_current=0*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
+    np.testing.assert_allclose(calculate_noise(camera,
                                                aperture_area=aperture_area),
                                expected)
 
@@ -62,7 +90,13 @@ def test_calc_sky_only(gain, aperture_area):
     sky = 10
     expected = np.sqrt(gain * aperture_area * sky)
 
-    np.testing.assert_allclose(calculate_noise(gain=gain,
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu,
+                    read_noise=0*u.electron,
+                    dark_current=0*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
+    np.testing.assert_allclose(calculate_noise(camera,
                                                aperture_area=aperture_area,
                                                sky_per_pix=sky),
                                expected)
@@ -78,7 +112,14 @@ def test_annulus_area_term():
     sky = 10
     expected = np.sqrt(gain * aperture_area *
                        (1 + aperture_area / annulus_area) * sky)
-    np.testing.assert_allclose(calculate_noise(gain=gain,
+
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu,
+                    read_noise=0*u.electron,
+                    dark_current=0*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
+    np.testing.assert_allclose(calculate_noise(camera,
                                                aperture_area=aperture_area,
                                                annulus_area=annulus_area,
                                                sky_per_pix=sky),
@@ -101,11 +142,15 @@ def test_calc_noise_messy_case(digit, expected):
     exposure = 18
     read_noise = 12
 
+    # Create camera instance
+    camera = Camera(gain=gain*u.electron/u.adu,
+                    read_noise=read_noise*u.electron,
+                    dark_current=dark_current*u.electron/u.second,
+                    pixel_scale=1*u.arcsec/u.pixel)
+
     np.testing.assert_allclose(
-        calculate_noise(counts=counts,
-                        gain=gain,
-                        dark_current_per_sec=dark_current,
-                        read_noise=read_noise,
+        calculate_noise(camera,
+                        counts=counts,
                         sky_per_pix=sky,
                         exposure=exposure,
                         aperture_area=aperture_area,
@@ -113,3 +158,179 @@ def test_calc_noise_messy_case(digit, expected):
                         include_digitization=digit),
         expected
     )
+
+
+def test_find_too_close():
+    # Load test sourcelist into memory
+    test_sl_data = ascii.read(get_pkg_data_filename('data/test_corner.ecsv'),
+                                format='ecsv',
+                                fast_reader=False)
+
+    # Create no sky position sourcelist
+    test_sl_data_nosky = test_sl_data.copy()
+    test_sl_data_nosky.remove_column('ra')
+    test_sl_data_nosky.remove_column('dec')
+
+    # Create no image position sourcelist
+    test_sl_data_noimgpos = test_sl_data.copy()
+    test_sl_data_noimgpos.remove_column('xcenter')
+    test_sl_data_noimgpos.remove_column('ycenter')
+
+    # Create SourceListData objects
+    sl_test = SourceListData(input_data=test_sl_data, colname_map=None)
+    sl_test_nosky = SourceListData(input_data=test_sl_data_nosky, colname_map=None)
+    sl_test_noimgpos = SourceListData(input_data=test_sl_data_noimgpos,
+                                      colname_map=None)
+
+    assert sl_test.has_ra_dec is True
+    assert sl_test.has_x_y is True
+    assert sl_test_nosky.has_ra_dec is False
+    assert sl_test_nosky.has_x_y is True
+    assert sl_test_noimgpos.has_ra_dec is True
+    assert sl_test_noimgpos.has_x_y is False
+
+    # Test full positions available
+    ap_in_asec = 5
+    feder_scale = 0.563
+    aperture_rad = ap_in_asec / feder_scale
+
+    rejects = find_too_close(sl_test, aperture_rad, pixel_scale=feder_scale)
+    assert np.sum(rejects) == 5
+
+    # Test only sky positions available
+    rejects = find_too_close(sl_test_noimgpos, aperture_rad, pixel_scale=feder_scale)
+    assert np.sum(rejects) == 5
+
+    # Test only image positions available
+    rejects = find_too_close(sl_test_nosky, aperture_rad, pixel_scale=feder_scale)
+    assert np.sum(rejects) == 5
+
+# Constants for the following tests
+shift_tolerance = 6
+max_adu = 60000
+fwhm_estimate = 5
+fake_camera = Camera(gain = 1.0*u.electron/u.adu,
+                     read_noise = 0*u.electron,
+                     dark_current = 0.1*u.electron/u.second,
+                     pixel_scale = 1*u.arcsec/u.pixel)
+fake_obs = EarthLocation(lat = 0*u.deg,
+                         lon = 0*u.deg,
+                         height = 0*u.m)
+coords2use='pixel'
+
+def test_aperture_photometry_no_outlier_rejection():
+    fake_CCDimage = FakeCCDImage()
+    sources = fake_CCDimage.sources
+    aperture = sources['aperture'][0]
+    inner_annulus = 2 * aperture
+    outer_annulus = 3 * aperture
+    found_sources = source_detection(fake_CCDimage,
+                                    fwhm=sources['x_stddev'].mean(),
+                                    threshold=10)
+    phot, missing_sources = single_image_photometry(fake_CCDimage,
+                                                    found_sources,
+                                                    fake_camera,
+                                                    fake_obs,
+                                                    aperture, inner_annulus,
+                                                    outer_annulus,
+                                                    shift_tolerance,
+                                                    max_adu, fwhm_estimate,
+                                                    use_coordinates=coords2use,
+                                                    include_dig_noise=True,
+                                                    reject_too_close=False,
+                                                    reject_background_outliers=False)
+
+    phot.sort('aperture_sum')
+    sources.sort('amplitude')
+    found_sources.sort('flux')
+
+    for inp, out in zip(sources, phot):
+        stdev = inp['x_stddev']
+        expected_flux = (inp['amplitude'] * 2 * np.pi *
+                            stdev**2 *
+                            (1 - np.exp(-aperture**2 / (2 * stdev**2))))
+        # This expected flux is correct IF there were no noise. With noise, the
+        # standard deviation in the sum of the noise within in the aperture is
+        # n_pix_in_aperture times the single-pixel standard deviation.
+
+        # We could require that the result be within some reasonable
+        # number of those expected variations or we could count up the
+        # actual number of background counts at each of the source
+        # positions.
+
+        # Here we just check whether any difference is consistent with
+        # less than the expected one sigma deviation.
+        assert (np.abs(expected_flux - out['aperture_net_cnts'].value) <
+                np.pi * aperture**2 * fake_CCDimage.noise_dev)
+
+
+@pytest.mark.parametrize('reject', [True, False])
+def test_aperture_photometry_with_outlier_rejection(reject):
+    """
+    Insert some really large pixel values in the annulus and check that
+    the photometry is correct when outliers are rejected and is
+    incorrect when outliers are not rejected.
+    """
+    fake_CCDimage = FakeCCDImage()
+    sources = fake_CCDimage.sources
+    aperture = sources['aperture'][0]
+    inner_annulus = 2 * aperture
+    outer_annulus = 3 * aperture
+    image = fake_CCDimage.data
+
+    found_sources = source_detection(fake_CCDimage,
+                                    fwhm=sources['x_stddev'].mean(),
+                                    threshold=10)
+
+    # Add some large pixel values to the annulus for each source.
+    # adding these moves the average pixel value by quite a bit,
+    # so we'll only get the correct net flux if these are removed.
+    for source in fake_CCDimage.sources:
+        center_px = (int(source['x_mean']), int(source['y_mean']))
+        begin = center_px[0] + inner_annulus + 1
+        end = begin + (outer_annulus - inner_annulus - 1)
+        # Yes, x and y are deliberately reversed below.
+        image[center_px[1], begin:end] = 100 * fake_CCDimage.mean_noise
+
+    phot, missing_sources = single_image_photometry(fake_CCDimage,
+                                                    found_sources,
+                                                    fake_camera,
+                                                    fake_obs,
+                                                    aperture, inner_annulus,
+                                                    outer_annulus,
+                                                    shift_tolerance,
+                                                    max_adu, fwhm_estimate,
+                                                    use_coordinates=coords2use,
+                                                    include_dig_noise=True,
+                                                    reject_too_close=False,
+                                                    reject_background_outliers=reject)
+
+    phot.sort('aperture_sum')
+    sources.sort('amplitude')
+    found_sources.sort('flux')
+
+    for inp, out in zip(sources, phot):
+        stdev = inp['x_stddev']
+        expected_flux = (inp['amplitude'] * 2 * np.pi *
+                         stdev**2 *
+                         (1 - np.exp(-aperture**2 / (2 * stdev**2))))
+        # This expected flux is correct IF there were no noise. With noise, the
+        # standard deviation in the sum of the noise within in the aperture is
+        # n_pix_in_aperture times the single-pixel standard deviation.
+        #
+
+        expected_deviation = np.pi * aperture**2 * fake_CCDimage.noise_dev
+        # We could require that the result be within some reasonable
+        # number of those expected variations or we could count up the
+        # actual number of background counts at each of the source
+        # positions.
+
+        # Here we just check whether any difference is consistent with
+        # less than the expected one sigma deviation.
+        if reject:
+            assert (np.abs(expected_flux - out['aperture_net_cnts'].value) <
+                    expected_deviation)
+        else:
+            with pytest.raises(AssertionError):
+                assert (np.abs(expected_flux - out['aperture_net_cnts'].value) <
+                        expected_deviation)
\ No newline at end of file
diff --git a/stellarphot/tests/test_core.py b/stellarphot/tests/test_core.py
index 143784e..7637fd6 100644
--- a/stellarphot/tests/test_core.py
+++ b/stellarphot/tests/test_core.py
@@ -151,7 +151,7 @@ photcoltypes = ['int', 'float', 'float', 'float', 'float', 'float', 'float', 'fl
                 'float', 'float', 'float', 'str', 'int', 'float', 'float', 'float',
                 'float', 'str', 'str', 'int', 'float', 'float', 'float', 'float']
 photcolunits = [None, u.pix, u.pix, u.adu, None, u.deg, u.deg, u.adu, u.adu, u.adu,
-                None, None, None, u.pix, u.pix*u.pix, u.pix, u.pix, u.pix*u.pix, u.s,
+                None, None, None, u.pix, u.pix, u.pix, u.pix, u.pix, u.s,
                 None, None, u.adu, None, None, None, None, None, None, None,
                 u.electron, None, u.adu]
 
@@ -165,13 +165,6 @@ testphot_goodTime['date-obs'] = Column(data=Time(testphot_goodTime['date-obs'],
                                                  format='isot', scale='utc'),
                                        name='date-obs')
 
-# Fix the units for the counts-related columns
-counts_columns = ['aperture_sum', 'annulus_sum', 'sky_per_pix_avg', 'sky_per_pix_med',
-                          'sky_per_pix_std', 'aperture_net_cnts', 'noise_cnts']
-testphot_goodCounts = testphot_goodTime.copy()
-for this_col in counts_columns:
-    testphot_goodCounts[this_col].unit = u.adu
-
 # Fix all the units for PhotometryData
 phot_descript = {
     'star_id' : None,
@@ -203,6 +196,15 @@ testphot_goodUnits = testphot_goodTime.copy()
 for this_col, this_unit in phot_descript.items():
     testphot_goodUnits[this_col].unit = this_unit
 
+# Fix the units for the counts-related columns
+counts_columns = ['aperture_sum', 'annulus_sum', 'aperture_net_cnts', 'noise_cnts']
+counts_per_pixel_sqr_columns = ['sky_per_pix_avg', 'sky_per_pix_med',
+                                            'sky_per_pix_std']
+for this_col in counts_columns:
+    testphot_goodUnits[this_col].unit = u.adu
+for this_col in counts_per_pixel_sqr_columns:
+    testphot_goodUnits[this_col].unit = u.adu * u.pixel**-1
+
 # Remove calculated columns from the test data to produce clean data
 computed_columns = ['bjd', 'night']
 testphot_clean = testphot_goodUnits.copy()
@@ -253,17 +255,17 @@ def test_photometry_slicing():
                                passband_map=feder_passbands, input_data=testphot_clean)
 
     # Test slicing works as expected, leaving attributes intact
-    just_cols = phot_data[['ra','dec']]
-    assert just_cols.camera.gain == 1.5 *  u.electron / u.adu
-    assert just_cols.camera.read_noise == 10.0 * u.electron
-    assert just_cols.camera.dark_current == 0.01 * u.electron / u.second
-    assert just_cols.camera.pixel_scale == 0.563 * u.arcsec / u.pix
-    assert just_cols.observatory.lat.value == 46.86678
-    assert just_cols.observatory.lat.unit == u.deg
-    assert just_cols.observatory.lon.value == -96.45328
-    assert just_cols.observatory.lon.unit == u.deg
-    assert round(just_cols.observatory.height.value) == 311
-    assert just_cols.observatory.height.unit == u.m
+    two_cols = phot_data[['ra','dec']]
+    assert two_cols.camera.gain == 1.5 *  u.electron / u.adu
+    assert two_cols.camera.read_noise == 10.0 * u.electron
+    assert two_cols.camera.dark_current == 0.01 * u.electron / u.second
+    assert two_cols.camera.pixel_scale == 0.563 * u.arcsec / u.pix
+    assert two_cols.observatory.lat.value == 46.86678
+    assert two_cols.observatory.lat.unit == u.deg
+    assert two_cols.observatory.lon.value == -96.45328
+    assert two_cols.observatory.lon.unit == u.deg
+    assert round(two_cols.observatory.height.value) == 311
+    assert two_cols.observatory.height.unit == u.m
 
 
 def test_photometry_recursive():
@@ -290,14 +292,6 @@ def test_photometry_inconsistent_count_units():
                                    passband_map=feder_passbands,
                                    input_data=testphot_goodTime)
 
-
-def test_photometry_inconsistent_badunits():
-    with pytest.raises(ValueError):
-        phot_data = PhotometryData(observatory=feder_obs, camera=feder_cg_16m,
-                                   passband_map=feder_passbands,
-                                   input_data=testphot_goodCounts)
-
-
 def test_photometry_inconsistent_computed_col_exists():
     with pytest.raises(ValueError):
         phot_data = PhotometryData(observatory=feder_obs, camera=feder_cg_16m,
@@ -425,6 +419,26 @@ def test_sourcelist_recursive():
         sl_test2 = SourceListData(input_data=sl_test, colname_map=None)
 
 
+def test_sourcelist_dropping_skycoords():
+     # Create good sourcelist data instance
+    sl_test = SourceListData(input_data=test_sl_data, colname_map=None)
+
+    # Drop sky coordinates
+    sl_test.drop_ra_dec()
+    assert sl_test.has_ra_dec == False
+    assert sl_test.has_x_y == True
+
+
+def test_sourcelist_dropping_imagecoords():
+     # Create good sourcelist data instance
+    sl_test = SourceListData(input_data=test_sl_data, colname_map=None)
+
+    # Drop sky coordinates
+    sl_test.drop_x_y()
+    assert sl_test.has_ra_dec == True
+    assert sl_test.has_x_y == False
+
+
 def test_sourcelist_slicing():
     # Create good sourcelist data instance
     sl_test = SourceListData(input_data=test_sl_data, colname_map=None)
@@ -439,4 +453,4 @@ def test_sourcelist_slicing():
     assert slicing_test['xcenter'][1] == sl_test['xcenter'][2]
     # Checking attributes survive slicing
     assert slicing_test.has_ra_dec == True
-    assert slicing_test.has_x_y == True
\ No newline at end of file
+    assert slicing_test.has_x_y == True