SEP uses the "core algorithms" of Source Extractor as a python library of functions and classes (see Barbary 2016).
with conda:
with pip
In [2]:
%matplotlib inline
import numpy as np
import sep
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy import wcs
from matplotlib.patches import Ellipse
We'll read in the fits image using astropy.io.fits. In doing so you might get the error:
You should then perform the byte swap operation below. If you're using fitsio this error probably won't occur
In [3]:
fname = 'simple_sim_cube_F090W_487_01.slp.fits'
hdu_list = fits.open(fname, do_not_scale_image_data=True)
tbdat = hdu_list[0].data
hdu_list.close()
tbdat = tbdat.byteswap().newbyteorder() #byte swap operation
In [4]:
plt.figure(figsize=[15,15])
plt.imshow(tbdat, interpolation='nearest', cmap='gray', vmin=0.093537331, vmax=0.21167476, origin='upper')
plt.colorbar()
Out[4]:
The background subtraction, sep.Background(), does allow for masking as well as filter dimension parameters. Let's find the background first.
In [5]:
bkg = sep.Background(tbdat)
print(bkg.globalback) #gives global mean of background
print(bkg.globalrms) #gives global noise of background
bkg_img = bkg.back() #2d array of background
print(bkg_img.shape)
tbdat_sub = tbdat - bkg
Plot the background
In [6]:
plt.figure(figsize=[15,15])
plt.imshow(bkg_img, interpolation='nearest', cmap='gray', origin = 'upper')
plt.colorbar()
Out[6]:
Here we perform the object detection, where 'thresh' is the detection threshold. A detected object is distinguished at this threshold (e.g 2.0*err).
In [7]:
objects = sep.extract(tbdat_sub, thresh=2.0, err = bkg.globalrms)
print'Number of objects detected: ', len(objects)
objects will have the following object parameters:
(see Baraby 2016)
We can also plot ellipses over the objects for a nice visual.
In [8]:
#pulled from Baraby2016
fig, ax = plt.subplots(figsize=[15,15])
m, s = np.mean(tbdat_sub), np.std(tbdat_sub)
im = ax.imshow(tbdat_sub, interpolation='nearest', cmap='gray', vmin=m-s, vmax=m+s, origin='upper')
for i in range(len(objects)):
e = Ellipse(xy=(objects['x'][i], objects['y'][i]), width = 6*objects['a'][i], height = 6*objects['b'][i], angle = objects['theta'][i]*108/np.pi)
e.set_edgecolor('red')
e.set_facecolor('none')
ax.add_artist(e)
We'll perform elliptical aperture photometry for objects within the Kron radius defined by $$\sum_{i} r_{i}I(r_{i})/\sum_{i}I(r_{i})$$ Here $r_{i}$ is the distance to the pixeld from the ellipse. he Kron aperature photometry is a proposed technique that captures the majority of the flux.
sep.kron_radius() uses an object's x, y, and elliptical properties, along with the elliptical radius (r) which is integrated over.
First let's find the Kron radius.
In [9]:
kronrad, kronflag = sep.kron_radius(tbdat_sub, objects['x'], objects['y'], objects['a'], objects['b'], objects['theta'], r=6.0)
Next we'll perform the elliptical aperture photometry.
In [10]:
flux, fluxerr, flag = sep.sum_ellipse(tbdat_sub, objects['x'], objects['y'], objects['a'], objects['b'], objects['theta'], r = (2.5*kronrad), err = bkg.globalrms, subpix=1)
flag |=kronflag #combines flags
r_min = 1.75 #minimum diameter = 3.5
Now we'll perform circular aperture photometry for sources that don't meet the Kron radius criteria. This process is equivalent to the FLUX_AUTO option in SExtractor.
In [13]:
use_circle = kronrad*np.sqrt(objects['a']*objects['b']) < r_min
cflux, cfluxerr, cflag = sep.sum_circle(tbdat_sub, objects['x'][use_circle], objects['y'][use_circle], r_min, subpix=1)
flux[use_circle] = cflux
fluxerr[use_circle] = cfluxerr
flag[use_circle] = cflag
We can find the radius of a circle contianing a percentage of the total flux of the object.
In [14]:
r, rflag = sep.flux_radius(tbdat_sub, objects['x'], objects['y'], rmax = 6.0*objects['a'], frac = 0.5, normflux = flux, subpix =5)
We can convert the flux to AB mags.
In [15]:
hdu_list = fits.open(fname, do_not_scale_image_data=True)
zero_point = hdu_list[0].header['ABMAG']
mag = -2.5*np.log10(flux) + zero_point
Let's get the RA & DEC of our extracted objects
In [22]:
hdu_list = fits.open(fname)
w = wcs.WCS(hdu_list[0].header)
hdu_list.close()
wrd = w.wcs_pix2world(objects['x'], objects['y'], np.zeros(len(objects['x'])), 0)
ra = wrd[:][0]
dec = wrd[:][1]
Print objects by increasing magnitude with coordinates.
In [29]:
index = sorted(range(len(mag)), key=lambda k: mag[k])
print "ABMag", "\t\t", "RA", "\t", "\t", "DEC"
for j in range(0,125):
k = index[j]
print mag[k], "\t", ra[k], "\t", dec[k]
In [ ]: