A convenience function to detect astronomical objects with SEP

sep_demo.py

#! /usr/bin/env python

# Author: Victor Terron (c) 2015
# Email: `echo vt2rron1iaa32s | tr 132 @.e`
# License: GNU GPLv3

import sep as sextractor
import sys
import astropy.io.fits
import astropy.wcs
from astropy import log
import warnings

def detect_sources(image_path):
    """ A generator of (ra, dec) tuples. """

    log.debug("Reading FITS file")
    with astropy.io.fits.open(path) as hdulist:
        data = hdulist[0].data
        header = hdulist[0].header

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        log.info("Loading WCS data")
        wcs = astropy.wcs.WCS(header)

    try:
        log.info("Estimating background")
        background = sextractor.Background(data)
    except ValueError:
        # Fix for error "Input array with dtype '>f4' has non-native
        # byte order. Only native byte order arrays are supported".
        log.debug("Converting data to native byte order")
        data = data.byteswap(True).newbyteorder()
        background = sextractor.Background(data)

    log.info("Subtracting background")
    background.subfrom(data) # in-place

    # Global "average" RMS of background
    log.info("Detecting sources")
    rms_ntimes = 1.5
    while True:
        try:
            threshold = rms_ntimes * background.globalrms
            objects = sextractor.extract(data, threshold)

        except Exception as e:

            # Fix for error "internal pixel buffer full: The limit of 300000
            # active object pixels over the detection threshold was reached.
            # Check that the image is background subtracted and the detection
            # threshold is not too low. detection threshold". If a different
            # exception is raised, just re-raise it.

            if "threshold is not too low" in str(e):
                rms_ntimes *= 1.25
                log.debug("Internal pixel buffer full")
                log.debug("Retrying with threshold {0:.2} times RMS of background".format(rms_ntimes))
            else:
                raise
        else:
            break

    pixel_coords = []

    # https://github.com/kbarbary/sep/blob/master/sep.pyx#L555
    log.info("Transforming pixel to celestial coordinates")
    for index in range(len(objects)):
        x, y = objects[index]['x'], objects[index]['y']
        pixel_coords.append((x, y))

    return iter(wcs.all_pix2world(pixel_coords, 1))

if __name__ == "__main__":

    if len(sys.argv) != 2:
        sys.exit("usage: {0} FITS_FILE".format(sys.argv[0]))

    path = sys.argv[1]
    for ra, dec in detect_sources(path):
        print ra, dec