SEP uses the "core algorithms" of Source Extractor as a python library of functions and classes (see Barbary 2016).

Install

with conda:

conda install -c openastronomy sep

with pip

pip install --no-deps set

Read image


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:

ValueError: Input array with dtype '>f4' has non-native byte order. Only native byte order arrays are supported. To change the byte order of the array 'data', do 'data = data.byteswap().newbyteorder()'

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

Plot original image


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]:
<matplotlib.colorbar.Colorbar at 0x10e72ffd0>

Subtract Background

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


0.137105196714
0.0122048910707
(2036, 2036)

Plot the background


In [6]:
plt.figure(figsize=[15,15])
plt.imshow(bkg_img, interpolation='nearest', cmap='gray', origin = 'upper')
plt.colorbar()


Out[6]:
<matplotlib.colorbar.Colorbar at 0x10f344350>

Extraction

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)


Number of objects detected:  155

objects will have the following object parameters:

  • thresh (float) Threshold at object location.
  • npix (int) Number of pixels belonging to the object.
  • tnpix (int) Number of pixels above threshold (unconvolved data).
  • xmin, xmax (int) Minimum, maximum x coordinates of pixels.
  • ymin, ymax (int) Minimum, maximum y coordinates of pixels.
  • x, y (float) object barycenter (first moments).
  • x2, y2, xy (float) Second moments.
  • errx2, erry2, errxy (float) Second moment errors. Note that these will be zero if error is not given.
  • a, b, theta (float) Ellipse parameters, scaled as described by Section 8.4.2 in “The Source Extractor Guide” or Section 10.1.5-6 of v2.13 of SExtractor’s User Manual.
  • cxx, cyy, cxy (float) Alternative ellipse parameters.
  • cflux (float) Sum of member pixels in convolved data.
  • flux (float) Sum of member pixels in unconvolved data.
  • cpeak (float) Peak value in convolved data.
  • peak (float) Peak value in unconvolved data.
  • xcpeak, ycpeak (int) Coordinate of convolved peak pixel.
  • xpeak, ypeak (int) Coordinate of unconvolved peak pixel.
  • flag (int) Extraction flags.

(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)


Let's do some photometry

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]


ABMag 		RA 		DEC
18.4338175303 	53.1179572388 	-27.8021195522
18.6922874897 	53.1038187917 	-27.8041250348
19.4646746852 	53.1016530735 	-27.8023052686
21.0839726351 	53.1186003597 	-27.8074310079
21.1223685721 	53.1064193419 	-27.8010303207
21.1612215212 	53.1110567529 	-27.7997004726
21.5547758626 	53.1201448023 	-27.7988975445
21.7200067276 	53.1198153904 	-27.7987930171
21.8422955969 	53.1124821238 	-27.8080752336
22.5117728202 	53.1164486117 	-27.8067989968
22.5919543589 	53.1058802678 	-27.808455173
22.6971715528 	53.103111525 	-27.8068496128
22.742517473 	53.1071695087 	-27.8058111926
22.8568470333 	53.120038775 	-27.8083262966
22.8640084283 	53.1139730663 	-27.8117612544
22.9021623923 	53.1087836893 	-27.8117738983
22.9541450837 	53.1117540581 	-27.8052975087
23.0088984105 	53.1110830653 	-27.8097221797
23.0268294414 	53.1082792878 	-27.7977005734
23.2603633296 	53.1070428929 	-27.8048023407
23.2724052034 	53.1017358241 	-27.8123211407
23.2837868966 	53.11065184 	-27.8006818351
23.3603667548 	53.1136202931 	-27.8127815252
23.4252163301 	53.1119776909 	-27.8010718972
23.4374732142 	53.1090918972 	-27.8118836546
23.5920143503 	53.1192116835 	-27.7970428192
23.6000536164 	53.1059608463 	-27.8079992082
23.6822692885 	53.1183854017 	-27.8053682187
23.7771946838 	53.1107458447 	-27.7973409537
23.9325989729 	53.1158060288 	-27.8034331343
23.9818859632 	53.1088147144 	-27.8073283915
24.110214653 	53.1108229386 	-27.8005383214
24.1173721941 	53.113647547 	-27.8096271278
24.148308445 	53.1160803569 	-27.8072500305
24.1878577663 	53.1172662739 	-27.8126628369
24.2501089534 	53.1076295533 	-27.8038291353
24.2734933984 	53.1094447947 	-27.8135700816
24.2757888423 	53.1132000808 	-27.8070176662
24.2852717717 	53.1079803706 	-27.7975957273
24.3938032569 	53.115451008 	-27.8051949158
24.4384242318 	53.1160986062 	-27.8049043187
24.4766185585 	53.1142180145 	-27.8128942183
24.5422997269 	53.1078004362 	-27.809970732
24.6615998965 	53.1147390551 	-27.8132150639
24.6860373124 	53.1207413522 	-27.8105804772
24.7916730306 	53.103972743 	-27.8023967752
24.8305641147 	53.1189438508 	-27.8024914715
24.8541090128 	53.1101409154 	-27.8068919016
24.8676122732 	53.1107358128 	-27.7981064549
24.8708158712 	53.1207054637 	-27.8033752699
24.940692134 	53.1187366162 	-27.8108837932
24.9460136426 	53.1188718701 	-27.8000879873
24.9946890647 	53.1082169266 	-27.7987237685
25.0443029322 	53.1155873844 	-27.8035868283
25.0629242478 	53.1040204096 	-27.8127991877
25.0792763487 	53.1104526312 	-27.8074350542
25.2191296982 	53.1204239163 	-27.8006960758
25.248217996 	53.1151812308 	-27.8058613447
25.2765501087 	53.1197090734 	-27.8064070336
25.2849104567 	53.121243765 	-27.8108890408
25.3340702828 	53.1081786455 	-27.8031694952
25.3823457636 	53.1093261001 	-27.808477345
25.3841941096 	53.1199807847 	-27.8021491134
25.4191824549 	53.1041537934 	-27.8114928809
25.4378767187 	53.116677446 	-27.7996805749
25.4536842388 	53.1181322804 	-27.8074343518
25.4922891691 	53.110271965 	-27.8081031237
25.5325897358 	53.1076325442 	-27.8045466884
25.5509665875 	53.119211108 	-27.798991343
25.5626831051 	53.1055845754 	-27.8095331881
25.5649302771 	53.1205308386 	-27.8050764745
25.6336759325 	53.1197214538 	-27.8118184303
25.6348077789 	53.1022537353 	-27.8095560143
25.6481839572 	53.1148149032 	-27.8102833025
25.6486872069 	53.1146693957 	-27.8007152143
25.6815461162 	53.1188504817 	-27.8042835533
25.6941111332 	53.1038019106 	-27.8003744799
25.7111588708 	53.114529045 	-27.8052577671
25.7395234515 	53.1033780847 	-27.8133287935
25.7781529737 	53.119041847 	-27.8119889904
25.8454784876 	53.1187972982 	-27.8060493831
25.8559298699 	53.1144466798 	-27.8014901968
25.8716642602 	53.1164459402 	-27.805892914
25.9216141861 	53.1060064947 	-27.8096538705
25.9307584822 	53.1025088329 	-27.8085977311
25.956161662 	53.1043222138 	-27.8103677073
25.9614209448 	53.1015421784 	-27.813247101
26.0058106855 	53.1195677557 	-27.8027586057
26.0311328981 	53.1030883773 	-27.8083924307
26.052825069 	53.1062616484 	-27.8069204022
26.067082083 	53.1115793832 	-27.8095887872
26.1243878897 	53.1165581938 	-27.8097556003
26.1441587229 	53.1210229038 	-27.8044259368
26.1724515271 	53.1017233597 	-27.7994550456
26.1786085533 	53.1072067663 	-27.80411162
26.1955808315 	53.1036394486 	-27.8099305892
26.2123228501 	53.1185761903 	-27.8057880022
26.2169310722 	53.1138482166 	-27.8061772326
26.228162971 	53.1022848884 	-27.8066916167
26.2660803827 	53.1094335085 	-27.8009980267
26.3082505009 	53.1133022915 	-27.7992157412
26.3886279515 	53.1173931524 	-27.8052437152
26.3972020125 	53.1121566674 	-27.8016385346
26.4400184156 	53.1118186758 	-27.8068259777
26.4481714819 	53.1145577599 	-27.8077881452
26.4831580535 	53.1104103815 	-27.8049796139
26.4976377947 	53.1112291525 	-27.8120793756
26.5579733363 	53.1069461098 	-27.8083779662
26.5977621309 	53.1186612223 	-27.8067529624
26.6189825932 	53.1058387019 	-27.8012114615
26.637994402 	53.1086080336 	-27.8002536888
26.6479793508 	53.1113660998 	-27.7973658534
26.6614782415 	53.1027372217 	-27.8069660038
26.6697970963 	53.1027044253 	-27.8093216268
26.6801516502 	53.1039250882 	-27.8058289629
26.776525823 	53.1123716509 	-27.7979631377
26.8165537692 	53.1153635597 	-27.8039606434
26.8670252395 	53.1147210521 	-27.8128092571
26.9003342065 	53.11828233 	-27.8127908442
26.9267283466 	53.1088421091 	-27.8092337372
26.9321785739 	53.1203518618 	-27.8039949552
26.9565590081 	53.1044384902 	-27.8064467384
26.9585844061 	53.1154324518 	-27.8041753635
26.9674466288 	53.1136955885 	-27.8052011538
26.9708598946 	53.1147853034 	-27.8052675959

In [ ]: