from IPython.core.display import SVG
url = 'https://stsci.box.com/shared/static/p2pmd5pzmmhs8hscoojnvlb4bcvnstck.svg'
SVG(url=url)
# initial imports
import numpy as np
import matplotlib.pyplot as plt
# run the %matplotlib magic command to enable inline plotting
# in the current Notebook
%matplotlib inline
import os
import os.path as op
from astropy.io import fits
#path = /home/analyst/data/glue'
path = op.join(os.environ['HOME'], 'data/glue')
sci_fn = 'hlsp_xdf_hst_wfc3ir-60mas_hudf_f160w_v1_sci.fits'
sci_fn = op.join(path, sci_fn)
sci_hdulist = fits.open(sci_fn)
# define the data array
data = sci_hdulist[0].data.astype(np.float)
print(data.shape)
print(data.dtype)
(5250, 5250) float64
# display the data
from astropy.visualization import scale_image
plt.imshow(scale_image(data, scale='sqrt', percent=98.),
origin='lower', cmap='Greys_r')
<matplotlib.image.AxesImage at 0x123fe21d0>
# load the weight file
wht_fn = 'hlsp_xdf_hst_wfc3ir-60mas_hudf_f160w_v1_wht.fits'
wht_fn = op.join(path, wht_fn)
#wht_fn = 'https://archive.stsci.edu/pub/hlsp/xdf/hlsp_xdf_hst_wfc3ir-60mas_hudf_f160w_v1_wht.fits'
wht = fits.getdata(wht_fn).astype(np.float)
# create an RMS map from the inverse-variance WHT map
err = 1. / np.sqrt(wht)
-c:9: RuntimeWarning: divide by zero encountered in divide
# extract the data header and create a WCS object
from astropy.wcs import WCS
hdr = sci_hdulist[0].header
wcs = WCS(hdr)
# plot a histogram of data values
idx_good = np.isfinite(err) # exclude regions with no coverage
r = plt.hist(data[idx_good].ravel(), bins=50, range=(-0.01, 0.01))
url = 'http://photutils.readthedocs.org/en/latest/_images/inheritance-41bd5c4e22ff4e2f493f41359f854c831e8e9c86.svg'
SVG(url=url)
url = 'https://stsci.box.com/shared/static/bhe6eyer7w5pmn9spxfjdx3celj9irzt.svg'
SVG(url=url)
from photutils import CircularAperture, aperture_photometry
# define the aperture
position = (2815.73, 4209.43)
radius = 5.
aperture = CircularAperture(position, r=radius)
# perform the photometry; the default method is 'exact'
phot = aperture_photometry(data, aperture)
phot
aperture_sum | xcenter [1] | ycenter [1] |
---|---|---|
float64 | float64 | float64 |
0.120391545825 | 2815.73 | 4209.43 |
# center method
phot = aperture_photometry(data, aperture, method='center')
phot
aperture_sum | xcenter [1] | ycenter [1] |
---|---|---|
float64 | float64 | float64 |
0.120292291098 | 2815.73 | 4209.43 |
# subpixel method, subpixels=5 (same as SExtractor)
phot = aperture_photometry(data, aperture, method='subpixel', subpixels=5)
phot
aperture_sum | xcenter [1] | ycenter [1] |
---|---|---|
float64 | float64 | float64 |
0.12040869806 | 2815.73 | 4209.43 |
# include units
import astropy.units as u
unit = u.electron / u.s
phot = aperture_photometry(data, aperture, unit=unit)
phot
aperture_sum | xcenter [1] | ycenter [1] |
---|---|---|
electron / s | ||
float64 | float64 | float64 |
0.120391545825 | 2815.73 | 4209.43 |
help(aperture_photometry)
Help on function aperture_photometry in module photutils.aperture_core: aperture_photometry(data, apertures, unit=None, wcs=None, error=None, effective_gain=None, mask=None, method=u'exact', subpixels=5, pixelwise_error=True) Sum flux within an aperture at the given position(s). Parameters ---------- data : array_like, `~astropy.io.fits.ImageHDU`, `~astropy.io.fits.HDUList` The 2-d array on which to perform photometry. ``data`` should be background-subtracted. Units are used during the photometry, either provided along with the data array, or stored in the header keyword ``'BUNIT'``. apertures : `~photutils.Aperture` instance The apertures to use for the photometry. unit : `~astropy.units.UnitBase` instance, str An object that represents the unit associated with ``data``. Must be an `~astropy.units.UnitBase` object or a string parseable by the :mod:`~astropy.units` package. An error is raised if ``data`` already has a different unit. wcs : `~astropy.wcs.WCS`, optional Use this as the wcs transformation. It overrides any wcs transformation passed along with ``data`` either in the header or in an attribute. error : float or array_like, optional Error in each pixel, interpreted as Gaussian 1-sigma uncertainty. effective_gain : float or array_like, optional Ratio of counts (e.g., electrons or photons) to units of the data (e.g., ADU), for the purpose of calculating Poisson error from the object itself. If ``effective_gain`` is `None` (default), ``error`` is assumed to include all uncertainty in each pixel. If ``effective_gain`` is given, ``error`` is assumed to be the "background error" only (not accounting for Poisson error in the flux in the apertures). mask : array_like (bool), optional Mask to apply to the data. Masked pixels are excluded/ignored. method : str, optional Method to use for determining overlap between the aperture and pixels. Options include ['center', 'subpixel', 'exact'], but not all options are available for all types of apertures. More precise methods will generally be slower. * ``'center'`` A pixel is considered to be entirely in or out of the aperture depending on whether its center is in or out of the aperture. * ``'subpixel'`` A pixel is divided into subpixels and the center of each subpixel is tested (as above). With ``subpixels`` set to 1, this method is equivalent to 'center'. Note that for subpixel sampling, the input array is only resampled once for each object. * ``'exact'`` (default) The exact overlap between the aperture and each pixel is calculated. subpixels : int, optional For the ``'subpixel'`` method, resample pixels by this factor (in each dimension). That is, each pixel is divided into ``subpixels ** 2`` subpixels. pixelwise_error : bool, optional For ``error`` and/or ``effective_gain`` arrays. If `True`, assume ``error`` and/or ``effective_gain`` vary significantly within an aperture: sum contribution from each pixel. If `False`, assume ``error`` and ``effective_gain`` do not vary significantly within an aperture. Use the single value of ``error`` and/or ``effective_gain`` at the center of each aperture as the value for the entire aperture. Default is `True`. Returns ------- phot_table : `~astropy.table.Table` A table of the photometry with the following columns: * ``'aperture_sum'``: Sum of the values within the aperture. * ``'aperture_sum_err'``: Corresponding uncertainty in ``'aperture_sum'`` values. Returned only if input ``error`` is not `None`. * ``'xcenter'``, ``'ycenter'``: x and y pixel coordinates of the center of the apertures. Unit is pixel. * ``'xcenter_input'``, ``'ycenter_input'``: input x and y coordinates as they were given in the input ``positions`` parameter. The metadata of the table stores the version numbers of both astropy and photutils, as well as the calling arguments. Notes ----- This function is decorated with `~astropy.nddata.support_nddata` and thus supports `~astropy.nddata.NDData` objects as input.
positions = [(2815.73, 4209.43), (2798.63, 4289.41), (2768.62, 4211.63)]
radius = 5.
apertures = CircularAperture(positions, r=radius)
phot = aperture_photometry(data, apertures, error=err)
phot
aperture_sum | aperture_sum_err | xcenter | ycenter |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | float64 |
0.120391545825 | 0.00802117560768 | 2815.73 | 4209.43 |
0.615024363298 | 0.00805158065863 | 2798.63 | 4289.41 |
0.213526159949 | 0.00802964643742 | 2768.62 | 4211.63 |
# this time include the Poisson error of the sources
eff_gain = hdr['TEXPTIME']
#toterr = np.sqrt(err**2 + (data / eff_gain))
#phot = aperture_photometry(data, apertures, error=toterr)
phot = aperture_photometry(data, apertures, error=err, effective_gain=eff_gain)
phot
aperture_sum | aperture_sum_err | xcenter | ycenter |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | float64 |
0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 |
0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 |
0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 |
# create a bad pixel
data2 = data.copy()
y, x = 4209, 2816
data2[y, x] = 100.
aperture_photometry(data2, apertures, error=err)
aperture_sum | aperture_sum_err | xcenter | ycenter |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | float64 |
100.113190738 | 0.00802117560768 | 2815.73 | 4209.43 |
0.615024363298 | 0.00805158065863 | 2798.63 | 4289.41 |
0.213526159949 | 0.00802964643742 | 2768.62 | 4211.63 |
# mask the bad pixel
mask = np.zeros_like(data2, dtype=bool)
mask[y, x] = True
aperture_photometry(data2, apertures, error=err, mask=mask)
aperture_sum | aperture_sum_err | xcenter | ycenter |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | float64 |
0.113190738027 | 0.00797011854574 | 2815.73 | 4209.43 |
0.615024363298 | 0.00805158065863 | 2798.63 | 4289.41 |
0.213526159949 | 0.00802964643742 | 2768.62 | 4211.63 |
# Add an ID column
phot['id'] = np.arange(len(phot)) + 1
phot
aperture_sum | aperture_sum_err | xcenter | ycenter | id |
---|---|---|---|---|
pix | pix | |||
float64 | float64 | float64 | float64 | int64 |
0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 | 1 |
0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 | 2 |
0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 | 3 |
# rearrange the column order to put 'id' first
colnames = phot.colnames
colnames.insert(0, colnames.pop(colnames.index('id')))
phot = phot[colnames]
phot
id | aperture_sum | aperture_sum_err | xcenter | ycenter |
---|---|---|---|---|
pix | pix | |||
int64 | float64 | float64 | float64 | float64 |
1 | 0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 |
2 | 0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 |
3 | 0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 |
# calculate the SNR and add it to the table
snr = phot['aperture_sum'] / phot['aperture_sum_err']
phot['snr'] = snr
phot
id | aperture_sum | aperture_sum_err | xcenter | ycenter | snr |
---|---|---|---|---|---|
pix | pix | ||||
int64 | float64 | float64 | float64 | float64 | float64 |
1 | 0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 | 14.9669210038 |
2 | 0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 | 75.312664422 |
3 | 0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 | 26.4600338963 |
# calculate the F160W AB magnitude and add it to the table
f160w_zpt = 25.9463
abmag = -2.5 * np.log10(phot['aperture_sum']) + f160w_zpt
phot['abmag'] = abmag
phot
id | aperture_sum | aperture_sum_err | xcenter | ycenter | snr | abmag |
---|---|---|---|---|---|---|
pix | pix | |||||
int64 | float64 | float64 | float64 | float64 | float64 | float64 |
1 | 0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 | 14.9669210038 | 28.2448100229 |
2 | 0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 | 75.312664422 | 26.4740691998 |
3 | 0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 | 26.4600338963 | 27.6226722755 |
# calculate the ICRS RA and Dec and add them to the table
from astropy.wcs.utils import pixel_to_skycoord
x, y = np.transpose(positions)
coord = pixel_to_skycoord(x, y, wcs)
phot['ra_icrs'] = coord.icrs.ra
phot['dec_icrs'] = coord.icrs.dec
phot
id | aperture_sum | aperture_sum_err | xcenter | ycenter | snr | abmag | ra_icrs | dec_icrs |
---|---|---|---|---|---|---|---|---|
pix | pix | deg | deg | |||||
int64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
1 | 0.120391545825 | 0.00804384186934 | 2815.73 | 4209.43 | 14.9669210038 | 28.2448100229 | 53.1588879497 | -27.7649992332 |
2 | 0.615024363298 | 0.00816628077121 | 2798.63 | 4289.41 | 75.312664422 | 26.4740691998 | 53.1592095849 | -27.76366615 |
3 | 0.213526159949 | 0.0080697613913 | 2768.62 | 4211.63 | 26.4600338963 | 27.6226722755 | 53.1597752653 | -27.7649623329 |
coord
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg [(53.1588915, -27.76500214), (53.15921314, -27.76366906), (53.15977882, -27.76496524)]>
# define circular aperture using sky coordinates
from photutils import SkyCircularAperture
radius2 = 5. * u.pix # sky aperture radius must be a Quantity
sky_apers = SkyCircularAperture(coord, r=radius2)
# aperture_photometry needs either an wcs object
phot2 = aperture_photometry(data, sky_apers, wcs=wcs)
phot2
aperture_sum | xcenter | ycenter | center_input |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | object |
0.120391545826 | 2815.73 | 4209.43 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.1588915, -27.76500214)> |
0.615024363294 | 2798.63 | 4289.41 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15921314, -27.76366906)> |
0.213526159949 | 2768.62 | 4211.63 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15977882, -27.76496524)> |
# ...or an hdu (i.e. header and data)
sci_hdulist[0].header['BUNIT'] = 'electron/s'
phot2 = aperture_photometry(sci_hdulist[0], sky_apers)
phot2
aperture_sum | xcenter | ycenter | center_input |
---|---|---|---|
electron / s | pix | pix | |
float64 | float64 | float64 | object |
0.120391545826 | 2815.73 | 4209.43 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.1588915, -27.76500214)> |
0.615024363294 | 2798.63 | 4289.41 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15921314, -27.76366906)> |
0.213526159949 | 2768.62 | 4211.63 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15977882, -27.76496524)> |
# SkyApertures can take angular radii
radius3 = 5. * u.arcsec
sky_apers3 = SkyCircularAperture(coord, r=radius3)
phot3 = aperture_photometry(data, sky_apers3, wcs=wcs)
phot3
aperture_sum | xcenter | ycenter | center_input |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | object |
4.54166803441 | 2815.73 | 4209.43 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.1588915, -27.76500214)> |
5.18998702492 | 2798.63 | 4289.41 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15921314, -27.76366906)> |
2.91213857764 | 2768.62 | 4211.63 | <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg (53.15977882, -27.76496524)> |
# can convert SkyApertures to PixelApertures
apers4 = sky_apers3.to_pixel(wcs)
apers4
<photutils.aperture_core.CircularAperture at 0x10c7dc350>
# get the radius in pixels
print(sky_apers3.r)
print(apers4.r) # radius in pixels
print(sky_apers3.r) / (0.060 * u.arcsec / u.pix)
5.0 arcsec 83.3333333283 83.3333333333 pix
plt.imshow(scale_image(data, scale='sqrt', percent=98.),
origin='lower', cmap='Greys_r')
apertures.plot(color='blue')
plt.xlim(2700, 2900)
plt.ylim(4150, 4350)
(4150, 4350)
unit = u.electron / u.s
radii = [3, 4, 5]
tbl_list = []
for radius in radii:
tbl_list.append(aperture_photometry(data,
CircularAperture(positions, radius),
error=err, unit=unit))
# combine the table list using hstack
from astropy.table import hstack
phot = hstack(tbl_list)
phot
aperture_sum_1 | aperture_sum_err_1 | xcenter_1 | ycenter_1 | aperture_sum_2 | aperture_sum_err_2 | xcenter_2 | ycenter_2 | aperture_sum_3 | aperture_sum_err_3 | xcenter_3 | ycenter_3 |
---|---|---|---|---|---|---|---|---|---|---|---|
electron / s | electron / s | pix | pix | electron / s | electron / s | pix | pix | electron / s | electron / s | pix | pix |
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
0.0866436609693 | 0.00480744058841 | 2815.73 | 4209.43 | 0.107911591525 | 0.00641468014996 | 2815.73 | 4209.43 | 0.120391545825 | 0.00802117560768 | 2815.73 | 4209.43 |
0.393646538117 | 0.00483447616294 | 2798.63 | 4289.41 | 0.520937788127 | 0.00643953948754 | 2798.63 | 4289.41 | 0.615024363298 | 0.00805158065863 | 2798.63 | 4289.41 |
0.130109734904 | 0.00480517791676 | 2768.62 | 4211.63 | 0.175370895678 | 0.0064132830568 | 2768.62 | 4211.63 | 0.213526159949 | 0.00802964643742 | 2768.62 | 4211.63 |
# instead of 1, 2, 3 names, use r3, r4, r5
from astropy.table import hstack
names = ['r3', 'r4', 'r5']
phot2 = hstack(tbl_list, table_names=names)
phot2
aperture_sum_r3 | aperture_sum_err_r3 | xcenter_r3 | ycenter_r3 | aperture_sum_r4 | aperture_sum_err_r4 | xcenter_r4 | ycenter_r4 | aperture_sum_r5 | aperture_sum_err_r5 | xcenter_r5 | ycenter_r5 |
---|---|---|---|---|---|---|---|---|---|---|---|
electron / s | electron / s | pix | pix | electron / s | electron / s | pix | pix | electron / s | electron / s | pix | pix |
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
0.0866436609693 | 0.00480744058841 | 2815.73 | 4209.43 | 0.107911591525 | 0.00641468014996 | 2815.73 | 4209.43 | 0.120391545825 | 0.00802117560768 | 2815.73 | 4209.43 |
0.393646538117 | 0.00483447616294 | 2798.63 | 4289.41 | 0.520937788127 | 0.00643953948754 | 2798.63 | 4289.41 | 0.615024363298 | 0.00805158065863 | 2798.63 | 4289.41 |
0.130109734904 | 0.00480517791676 | 2768.62 | 4211.63 | 0.175370895678 | 0.0064132830568 | 2768.62 | 4211.63 | 0.213526159949 | 0.00802964643742 | 2768.62 | 4211.63 |
# do some table cleanup
for key in ['xcenter_r4', 'ycenter_r4', 'xcenter_r5', 'ycenter_r5']:
del phot2[key]
phot2.rename_column('xcenter_r3', 'xcenter')
phot2.rename_column('ycenter_r3', 'ycenter')
colnames = phot2.colnames
colnames.insert(0, colnames.pop(colnames.index('ycenter')))
colnames.insert(0, colnames.pop(colnames.index('xcenter')))
phot2 = phot2[colnames]
phot2
xcenter | ycenter | aperture_sum_r3 | aperture_sum_err_r3 | aperture_sum_r4 | aperture_sum_err_r4 | aperture_sum_r5 | aperture_sum_err_r5 |
---|---|---|---|---|---|---|---|
pix | pix | electron / s | electron / s | electron / s | electron / s | electron / s | electron / s |
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
2815.73 | 4209.43 | 0.0866436609693 | 0.00480744058841 | 0.107911591525 | 0.00641468014996 | 0.120391545825 | 0.00802117560768 |
2798.63 | 4289.41 | 0.393646538117 | 0.00483447616294 | 0.520937788127 | 0.00643953948754 | 0.615024363298 | 0.00805158065863 |
2768.62 | 4211.63 | 0.130109734904 | 0.00480517791676 | 0.175370895678 | 0.0064132830568 | 0.213526159949 | 0.00802964643742 |
# add the aperture radii to the table
phot2['aperture_radius_r3'] = radii[0] * u.pix
phot2['aperture_radius_r4'] = radii[1] * u.pix
phot2['aperture_radius_r5'] = radii[2] * u.pix
phot2
xcenter | ycenter | aperture_sum_r3 | aperture_sum_err_r3 | aperture_sum_r4 | aperture_sum_err_r4 | aperture_sum_r5 | aperture_sum_err_r5 | aperture_radius_r3 | aperture_radius_r4 | aperture_radius_r5 |
---|---|---|---|---|---|---|---|---|---|---|
pix | pix | electron / s | electron / s | electron / s | electron / s | electron / s | electron / s | pix | pix | pix |
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
2815.73 | 4209.43 | 0.0866436609693 | 0.00480744058841 | 0.107911591525 | 0.00641468014996 | 0.120391545825 | 0.00802117560768 | 3.0 | 4.0 | 5.0 |
2798.63 | 4289.41 | 0.393646538117 | 0.00483447616294 | 0.520937788127 | 0.00643953948754 | 0.615024363298 | 0.00805158065863 | 3.0 | 4.0 | 5.0 |
2768.62 | 4211.63 | 0.130109734904 | 0.00480517791676 | 0.175370895678 | 0.0064132830568 | 0.213526159949 | 0.00802964643742 | 3.0 | 4.0 | 5.0 |
# multiple apertures at each position
radii = np.linspace(0.1, 20, 100)
flux = []
for r in radii:
ap = CircularAperture(positions[0], r=r)
phot = aperture_photometry(data, ap)
flux.append(phot['aperture_sum'][0])
plt.plot(radii, flux)
plt.title('Encircled Flux')
plt.xlabel('Radius (pixels)')
plt.ylabel('Aperture Sum ($e^{-1}/s$)')
<matplotlib.text.Text at 0x10e3ef950>
# define the annular apertures
from photutils import CircularAnnulus
apertures = CircularAperture(positions, r=3)
annulus_apertures = CircularAnnulus(positions, r_in=10., r_out=15.)
# plot the apertures
from astropy.visualization import scale_image
plt.imshow(scale_image(data, scale='sqrt', percent=98.),
origin='lower', cmap='Greys_r')
apertures.plot(color='blue')
annulus_apertures.plot(color='green', hatch='//', alpha=0.8)
plt.xlim(2700, 2900)
plt.ylim(4150, 4350)
(4150, 4350)
# measure the aperture sums
phot = aperture_photometry(data, apertures)
bkg = aperture_photometry(data, annulus_apertures)
# calculate the mean background level (per pixel) in the annuli
bkg_mean = bkg['aperture_sum'] / annulus_apertures.area()
bkg_mean
5.07023247856e-05 |
9.13944855008e-05 |
4.24459249229e-05 |
# now calculate the total background in the circular aperture
bkg_sum = bkg_mean * apertures.area()
phot['bkg_sum'] = bkg_sum
phot
aperture_sum | xcenter | ycenter | bkg_sum |
---|---|---|---|
pix | pix | ||
float64 | float64 | float64 | float64 |
0.0866436609693 | 2815.73 | 4209.43 | 0.0014335744596 |
0.393646538117 | 2798.63 | 4289.41 | 0.00258411819805 |
0.130109734904 | 2768.62 | 4211.63 | 0.00120013025321 |
# subtract the background
flux_bkgsub = phot['aperture_sum'] - bkg_sum
phot['aperture_sum_bkgsub'] = flux_bkgsub
phot
aperture_sum | xcenter | ycenter | bkg_sum | aperture_sum_bkgsub |
---|---|---|---|---|
pix | pix | |||
float64 | float64 | float64 | float64 | float64 |
0.0866436609693 | 2815.73 | 4209.43 | 0.0014335744596 | 0.0852100865097 |
0.393646538117 | 2798.63 | 4289.41 | 0.00258411819805 | 0.391062419919 |
0.130109734904 | 2768.62 | 4211.63 | 0.00120013025321 | 0.128909604651 |
# define a threshold image
bkg = 0.
threshold = bkg + (3.0 * err) # 3-sigma (per pixel!) above the background
# detect sources
from photutils import detect_sources
segm = detect_sources(data, threshold, npixels=5)
# object minimum snr
obj_min_snr = np.sqrt(5) * 3.
obj_min_snr
6.7082039324993694
# plot the segmentation image
plt.imshow(segm, origin='lower')
plt.colorbar()
print('Found {0} sources'.format(np.max(segm)))
Found 2737 sources
# define the kernel
from astropy.convolution import Gaussian2DKernel
from astropy.stats import gaussian_fwhm_to_sigma
sigma = 2.0 * gaussian_fwhm_to_sigma # FWHM = 2 pixels
kernel = Gaussian2DKernel(sigma, x_size=5, y_size=5)
kernel.array
array([[ 0.00086186, 0.00689486, 0.01378973, 0.00689486, 0.00086186], [ 0.00689486, 0.0551589 , 0.1103178 , 0.0551589 , 0.00689486], [ 0.01378973, 0.1103178 , 0.2206356 , 0.1103178 , 0.01378973], [ 0.00689486, 0.0551589 , 0.1103178 , 0.0551589 , 0.00689486], [ 0.00086186, 0.00689486, 0.01378973, 0.00689486, 0.00086186]])
# detect the sources
segm = detect_sources(data, threshold, npixels=5, filter_kernel=kernel)
plt.imshow(segm)
plt.colorbar()
print('Found {0} sources'.format(np.max(segm)))
Found 2409 sources
from photutils import segment_properties, properties_table
props = segment_properties(data, segm, error=err)
tbl = properties_table(props)
tbl
id | xcentroid | ycentroid | ra_icrs_centroid | dec_icrs_centroid | segment_sum | segment_sum_err | background_sum | background_mean | background_atcentroid | xmin | xmax | ymin | ymax | min_value | max_value | minval_xpos | minval_ypos | maxval_xpos | maxval_ypos | area | equivalent_radius | perimeter | semimajor_axis_sigma | semiminor_axis_sigma | eccentricity | orientation | ellipticity | elongation | covar_sigx2 | covar_sigxy | covar_sigy2 | cxx | cxy | cyy |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
pix | pix | pix | pix | pix | pix | pix | pix | pix | pix | pix2 | pix | pix | pix | pix | rad | pix2 | pix2 | pix2 | 1 / pix2 | 1 / pix2 | 1 / pix2 | |||||||||||||
int64 | float64 | float64 | object | object | float64 | float64 | object | object | object | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
1 | 2465.87081094 | 399.726636976 | None | None | 19.3759366767 | 0.117206007737 | None | None | None | 2459.0 | 2472.0 | 372.0 | 426.0 | 0.0140414275229 | 0.16770106554 | 2460.0 | 402.0 | 2466.0 | 400.0 | 494.0 | 12.5397401797 | 122.183766184 | 11.4193004726 | 2.39595712088 | 0.977740791852 | -1.4520215309 | 0.790183547002 | 4.76607046641 | 7.4909884669 | -14.667582062 | 128.650045341 | 0.171859229177 | 0.039187850115 | 0.0100069572481 |
2 | 2525.42098782 | 393.588988417 | None | None | 1.98278233688 | 0.0393283794627 | None | None | None | 2522.0 | 2530.0 | 390.0 | 397.0 | 0.0138792432845 | 0.107574552298 | 2526.0 | 390.0 | 2525.0 | 393.0 | 53.0 | 4.10736216662 | 24.7279220614 | 1.96965166431 | 1.47909709675 | 0.660366469015 | 0.527680534087 | 0.249056509049 | 1.33165812348 | 3.45058343464 | 0.7359989638 | 2.61667246568 | 0.308302566847 | -0.173434293143 | 0.406555938495 |
3 | 2659.53286882 | 454.227861931 | None | None | 1.18773694895 | 0.0355890020075 | None | None | None | 2656.0 | 2663.0 | 451.0 | 458.0 | 0.0167008060962 | 0.0429179444909 | 2659.0 | 458.0 | 2660.0 | 454.0 | 44.0 | 3.74241031851 | 23.313708499 | 2.2863305038 | 1.30419071985 | 0.821346461013 | -0.702741514254 | 0.429570345282 | 1.75306453956 | 3.75426438872 | -1.73915893667 | 3.17395621763 | 0.356977381157 | 0.391209178739 | 0.422245102882 |
4 | 2553.39549224 | 458.68878555 | None | None | 1.98222968075 | 0.045835992135 | None | None | None | 2548.0 | 2559.0 | 455.0 | 462.0 | 0.0142367426306 | 0.0467581525445 | 2555.0 | 462.0 | 2553.0 | 459.0 | 75.0 | 4.88602511903 | 30.1421356237 | 2.83650150654 | 1.75330029915 | 0.786083201191 | 0.263646185028 | 0.381879299162 | 1.6178070054 | 7.70809594701 | 1.25086278848 | 3.41170678858 | 0.137940897925 | -0.1011488659 | 0.311650954233 |
5 | 2634.03819015 | 492.064461678 | None | None | 14.0897114333 | 0.0594858061938 | None | None | None | 2629.0 | 2639.0 | 487.0 | 498.0 | 0.0124509194866 | 0.937757313251 | 2631.0 | 487.0 | 2634.0 | 492.0 | 111.0 | 5.94410610323 | 36.1421356237 | 1.9692171003 | 1.92742100098 | 0.204936499671 | 1.43341001729 | 0.0212247290096 | 1.02168498698 | 3.7180064855 | 0.0220948262943 | 3.87476121764 | 0.268970452714 | -0.00306746924377 | 0.258089165094 |
6 | 2677.36215021 | 495.847499234 | None | None | 0.180360917002 | 0.0147379330151 | None | None | None | 2676.0 | 2678.0 | 495.0 | 497.0 | 0.0187366288155 | 0.0316804014146 | 2676.0 | 495.0 | 2678.0 | 496.0 | 7.0 | 1.49270533036 | 6.20710678119 | 0.850405435895 | 0.590477573242 | 0.71963909768 | 1.05814969591 | 0.305651694687 | 1.44019938171 | 0.438765605852 | 0.160084690283 | 0.633087564049 | 2.51075833176 | -1.2697579063 | 1.74009799456 |
7 | 2585.36607904 | 507.279808167 | None | None | 1.35954320803 | 0.0332199331556 | None | None | None | 2582.0 | 2589.0 | 504.0 | 510.0 | 0.013419944793 | 0.080189101398 | 2589.0 | 507.0 | 2585.0 | 507.0 | 40.0 | 3.56824823231 | 21.8994949366 | 1.5806193275 | 1.44046177895 | 0.411682257834 | -0.640956158222 | 0.0886725513931 | 1.09730042865 | 2.34695734778 | -0.202940754602 | 2.22633024728 | 0.429468714271 | 0.0782962950432 | 0.452738202624 |
8 | 2444.01454498 | 543.504578668 | None | None | 5.52228014916 | 0.0583988152155 | None | None | None | 2438.0 | 2450.0 | 538.0 | 550.0 | 0.0149926077574 | 0.150368094444 | 2445.0 | 550.0 | 2444.0 | 544.0 | 123.0 | 6.25716517287 | 39.7989898732 | 2.86791213588 | 2.23140955349 | 0.628189175028 | 0.728803760024 | 0.221939359447 | 1.28524686622 | 6.78535256101 | 1.61248096699 | 6.41875605354 | 0.156733059323 | -0.0787470572017 | 0.165684605647 |
9 | 2615.78157114 | 559.569130829 | None | None | 2.64729802497 | 0.0499982890498 | None | None | None | 2610.0 | 2622.0 | 555.0 | 564.0 | 0.0127965640277 | 0.0616609677672 | 2620.0 | 558.0 | 2615.0 | 559.0 | 89.0 | 5.32255388609 | 35.7989898732 | 3.11710753315 | 1.76953470132 | 0.823245964919 | 0.52380650105 | 0.43231515676 | 1.76154077726 | 8.06889802156 | 2.85211837829 | 4.77871441083 | 0.157068588055 | -0.187489005674 | 0.265211584214 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
2400 | 2954.50459075 | 4663.40108932 | None | None | 0.90555095952 | 0.024554066751 | None | None | None | 2952.0 | 2957.0 | 4660.0 | 4667.0 | 0.014181945473 | 0.0761082842946 | 2952.0 | 4663.0 | 2956.0 | 4663.0 | 29.0 | 3.03825388987 | 19.8994949366 | 2.13677091142 | 0.919719512154 | 0.902626369348 | -0.953635180475 | 0.56957505026 | 2.32328539646 | 2.09175319704 | -1.75565549076 | 3.31992071191 | 0.859609062631 | 0.909164707038 | 0.541606309606 |
2401 | 2796.88545675 | 4675.15644346 | None | None | 0.250964800827 | 0.0168687496741 | None | None | None | 2796.0 | 2798.0 | 4673.0 | 4677.0 | 0.0131820207462 | 0.0220986753702 | 2798.0 | 4674.0 | 2797.0 | 4675.0 | 14.0 | 2.11100412282 | 11.4142135624 | 1.33050410686 | 0.756443209504 | 0.822656448265 | 1.41410898484 | 0.431461199101 | 1.75889490465 | 0.601379272947 | 0.184659531598 | 1.74106823462 | 1.71882117723 | -0.364599970495 | 0.59369495654 |
2402 | 2630.47744623 | 4712.51646276 | None | None | 3.39533200581 | 0.0484386460895 | None | None | None | 2626.0 | 2635.0 | 4706.0 | 4719.0 | 0.0117573924363 | 0.0804541036487 | 2631.0 | 4718.0 | 2631.0 | 4712.0 | 102.0 | 5.69803548521 | 36.3847763109 | 2.76911541442 | 2.09244119699 | 0.654992656961 | -1.19357677139 | 0.244364757752 | 1.32338983691 | 4.82462931623 | -1.12652158255 | 7.22168102498 | 0.215104583106 | 0.0671090164564 | 0.143706136303 |
2403 | 2452.8236964 | 4709.96896329 | None | None | 1.25633653346 | 0.0286567356471 | None | None | None | 2449.0 | 2456.0 | 4707.0 | 4713.0 | 0.0134650403634 | 0.10125181824 | 2450.0 | 4707.0 | 2452.0 | 4711.0 | 36.0 | 3.38513750129 | 21.313708499 | 1.93691968624 | 1.17560048859 | 0.794745111018 | 0.662916838818 | 0.393056667791 | 1.64760027326 | 2.85418758189 | 1.14943979796 | 2.27950679783 | 0.439640961959 | -0.443377329666 | 0.550477750411 |
2404 | 2877.94634142 | 4736.98986036 | None | None | 1.36224744003 | 0.0328719461685 | None | None | None | 2874.0 | 2882.0 | 4733.0 | 4741.0 | 0.0111631425098 | 0.0540334023535 | 2880.0 | 4740.0 | 2878.0 | 4737.0 | 55.0 | 4.18414193594 | 25.5563491861 | 2.00382361671 | 1.68106812787 | 0.544239043059 | -1.34077728592 | 0.161069809811 | 1.19199429428 | 2.88781347616 | -0.264018258899 | 3.95348566125 | 0.348409956941 | 0.0465344246047 | 0.25449516075 |
2405 | 2549.41023647 | 4735.8273186 | None | None | 0.0924080368131 | 0.0105660669828 | None | None | None | 2549.0 | 2550.0 | 4735.0 | 4737.0 | 0.0154587738216 | 0.0201562419534 | 2549.0 | 4735.0 | 2549.0 | 4736.0 | 5.0 | 1.26156626101 | 5.20710678119 | 0.777657360408 | 0.434465643547 | 0.829379809236 | -1.20514942011 | 0.441314818496 | 1.78991681381 | 0.241942507859 | -0.138906156943 | 0.551568857761 | 4.83183711206 | 2.43368317398 | 2.11945756547 |
2406 | 2505.76276783 | 4750.36601243 | None | None | 4.6651768852 | 0.0621128438764 | None | None | None | 2498.0 | 2518.0 | 4744.0 | 4758.0 | 0.0128603102639 | 0.0700320526958 | 2512.0 | 4751.0 | 2505.0 | 4750.0 | 172.0 | 7.39927702033 | 58.941125497 | 3.82499273606 | 2.99642159785 | 0.621543459408 | -0.0985667623724 | 0.21662031679 | 1.27652021291 | 14.5758352896 | -0.553500695801 | 9.03327653337 | 0.0687667086514 | 0.00842715728803 | 0.1109599839 |
2407 | 2694.19677544 | 4778.05236349 | None | None | 1.29121355433 | 0.0279628272846 | None | None | None | 2691.0 | 2697.0 | 4775.0 | 4781.0 | 0.0092290667817 | 0.115526095033 | 2692.0 | 4776.0 | 2694.0 | 4778.0 | 35.0 | 3.33779058906 | 19.313708499 | 1.3560784345 | 1.25932997781 | 0.37094282203 | -0.769145525155 | 0.0713442926549 | 1.07682534236 | 1.71654214699 | -0.126451530434 | 1.70831856654 | 0.585760602976 | 0.0867172272961 | 0.588580363609 |
2408 | 2718.17875847 | 4776.7593573 | None | None | 0.228583076969 | 0.0165781042621 | None | None | None | 2716.0 | 2720.0 | 4775.0 | 4778.0 | 0.0140871666372 | 0.0240108538419 | 2716.0 | 4776.0 | 2719.0 | 4777.0 | 12.0 | 1.95441004761 | 9.65685424949 | 1.30597263094 | 0.662986230977 | 0.861559229401 | 0.564139569439 | 0.492342936389 | 1.96983371587 | 1.34362101477 | 0.57203352767 | 0.801494240445 | 1.06911207995 | -1.52606949297 | 1.79225424873 |
2409 | 2570.96400717 | 4797.48344787 | None | None | 0.137088250369 | 0.0118955943912 | None | None | None | 2570.0 | 2572.0 | 4797.0 | 4798.0 | 0.0190680138767 | 0.0281395856291 | 2572.0 | 4798.0 | 2571.0 | 4797.0 | 6.0 | 1.38197659789 | 6.0 | 0.778941355061 | 0.499722123831 | 0.767089489726 | -0.00327358238764 | 0.358459888432 | 1.55874898852 | 0.60674580861 | -0.0011687503686 | 0.24972602706 | 1.64815150107 | 0.0154271278575 | 4.0044244767 |
snr = tbl['segment_sum'] / tbl['segment_sum_err']
print('Minimum SNR: {0}'.format(np.min(snr)))
Minimum SNR: 7.3668086159
# a subset of sources can be specified,
# defined by their labels in the segmentation image
labels = [1, 5, 20, 50, 75, 80]
props2 = segment_properties(data, segm, error=err, wcs=wcs, labels=labels)
tbl = properties_table(props2)
tbl
id | xcentroid | ycentroid | ra_icrs_centroid | dec_icrs_centroid | segment_sum | segment_sum_err | background_sum | background_mean | background_atcentroid | xmin | xmax | ymin | ymax | min_value | max_value | minval_xpos | minval_ypos | maxval_xpos | maxval_ypos | area | equivalent_radius | perimeter | semimajor_axis_sigma | semiminor_axis_sigma | eccentricity | orientation | ellipticity | elongation | covar_sigx2 | covar_sigxy | covar_sigy2 | cxx | cxy | cyy |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
pix | pix | deg | deg | pix | pix | pix | pix | pix | pix | pix | pix | pix2 | pix | pix | pix | pix | rad | pix2 | pix2 | pix2 | 1 / pix2 | 1 / pix2 | 1 / pix2 | |||||||||||
int64 | float64 | float64 | float64 | float64 | float64 | float64 | object | object | object | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
1 | 2465.87081094 | 399.726636976 | 53.1655214682 | -27.8285090544 | 19.3759366767 | 0.117206007737 | None | None | None | 2459.0 | 2472.0 | 372.0 | 426.0 | 0.0140414275229 | 0.16770106554 | 2460.0 | 402.0 | 2466.0 | 400.0 | 494.0 | 12.5397401797 | 122.183766184 | 11.4193004726 | 2.39595712088 | 0.977740791852 | -1.4520215309 | 0.790183547002 | 4.76607046641 | 7.4909884669 | -14.667582062 | 128.650045341 | 0.171859229177 | 0.039187850115 | 0.0100069572481 |
5 | 2634.03819015 | 492.064461678 | 53.1623515826 | -27.8269710321 | 14.0897114333 | 0.0594858061938 | None | None | None | 2629.0 | 2639.0 | 487.0 | 498.0 | 0.0124509194866 | 0.937757313251 | 2631.0 | 487.0 | 2634.0 | 492.0 | 111.0 | 5.94410610323 | 36.1421356237 | 1.9692171003 | 1.92742100098 | 0.204936499671 | 1.43341001729 | 0.0212247290096 | 1.02168498698 | 3.7180064855 | 0.0220948262943 | 3.87476121764 | 0.268970452714 | -0.00306746924377 | 0.258089165094 |
20 | 2327.89478561 | 717.443551847 | 53.1681195861 | -27.8232129488 | 5.09327924345 | 0.0634265023838 | None | None | None | 2323.0 | 2333.0 | 709.0 | 725.0 | 0.0142875546589 | 0.101232789457 | 2332.0 | 722.0 | 2328.0 | 717.0 | 143.0 | 6.74672614861 | 43.7989898732 | 3.40785379096 | 2.40184286615 | 0.709409701436 | -1.38884042562 | 0.295203663809 | 1.41884960044 | 5.9602263893 | -1.04014501329 | 11.4220902249 | 0.170488266907 | 0.0310508001875 | 0.0889634600566 |
50 | 2654.19478477 | 964.965942774 | 53.161968865 | -27.8190894513 | 0.806716094725 | 0.0301505344762 | None | None | None | 2650.0 | 2658.0 | 962.0 | 968.0 | 0.0137104261667 | 0.0319027081132 | 2652.0 | 966.0 | 2654.0 | 965.0 | 37.0 | 3.43183125879 | 23.1066017178 | 2.48623246998 | 1.03722743727 | 0.908820008767 | 0.655381022095 | 0.582811563362 | 2.39699836376 | 4.28494470958 | 2.4669348987 | 2.97224794181 | 0.446944476139 | -0.741918539471 | 0.644338024956 |
75 | 2898.58539606 | 1116.73027789 | 53.1573626114 | -27.8165612716 | 0.122965287417 | 0.0132599867602 | None | None | None | 2898.0 | 2899.0 | 1115.0 | 1118.0 | 0.0160002671182 | 0.0187656655908 | 2898.0 | 1116.0 | 2899.0 | 1117.0 | 7.0 | 1.49270533036 | 7.20710678119 | 1.03257280072 | 0.473051042353 | 0.88888596853 | -1.42032587794 | 0.541871486426 | 2.18279362749 | 0.242707512418 | -0.124856004144 | 1.04727636505 | 4.38938760515 | 1.04660319913 | 1.01724566909 |
80 | 1755.70949679 | 1151.83352703 | 53.1788988607 | -27.8159691752 | 1.43870618287 | 0.034281153463 | None | None | None | 1753.0 | 1759.0 | 1148.0 | 1155.0 | 0.0126629555598 | 0.0654352232814 | 1759.0 | 1153.0 | 1756.0 | 1152.0 | 45.0 | 3.78469878303 | 22.4852813742 | 1.82258252039 | 1.49005603819 | 0.575854738034 | 1.36037865658 | 0.182447970656 | 1.22316374262 | 2.26832278138 | 0.225002250232 | 3.2737512592 | 0.443880506986 | -0.0610150894178 | 0.307556740415 |
# a subset of columns can be specified
columns = ['id', 'xcentroid', 'ycentroid', 'segment_sum', 'area']
tbl2 = properties_table(props2, columns=columns)
tbl2
id | xcentroid | ycentroid | segment_sum | area |
---|---|---|---|---|
pix | pix | pix2 | ||
int64 | float64 | float64 | float64 | float64 |
1 | 2465.87081094 | 399.726636976 | 19.3759366767 | 494.0 |
5 | 2634.03819015 | 492.064461678 | 14.0897114333 | 111.0 |
20 | 2327.89478561 | 717.443551847 | 5.09327924345 | 143.0 |
50 | 2654.19478477 | 964.965942774 | 0.806716094725 | 37.0 |
75 | 2898.58539606 | 1116.73027789 | 0.122965287417 | 7.0 |
80 | 1755.70949679 | 1151.83352703 | 1.43870618287 | 45.0 |
# props is a list of SegmentProperties objects
obj = props[0] # id = 1 (segmentation label)
print(obj.id)
print(obj.xcentroid)
print(obj.area)
print(obj.ellipticity)
print(obj.covariance)
print(obj.orientation)
print(obj.orientation.to(u.deg))
1 2465.87081094 pix 494.0 pix2 0.790183547002 [[ 7.49098847 -14.66758206] [ -14.66758206 128.65004534]] pix2 -1.4520215309 rad -83.1947054826 deg
# plot the first object
fig, ax = plt.subplots(figsize=(4, 6), ncols=2)
origin ='lower'
cmap = 'hot'
ax[0].imshow(obj.data_cutout_ma, origin=origin, cmap=cmap)
ax[1].imshow(obj.error_cutout_ma, origin=origin, cmap=cmap)
<matplotlib.image.AxesImage at 0x1461e4f50>
f160w_zpt = 25.9463
abmag = -2.5 * np.log10(obj.segment_sum) + f160w_zpt
print(abmag)
print(obj.segment_sum / obj.segment_sum_err)
22.7281432338 165.315217631
y0, y1 = 2700, 3000
x0, x1 = 3350, 3650
from photutils import EllipticalAperture
r = 3. # approximate isophotal extent
apertures = []
for prop in props:
position = (prop.xcentroid.value - x0, prop.ycentroid.value - y0)
a = prop.semimajor_axis_sigma.value * r
b = prop.semiminor_axis_sigma.value * r
theta = prop.orientation.value
apertures.append(EllipticalAperture(position, a, b, theta=theta))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8))
ax1.imshow(scale_image(data[y0:y1, x0:x1], scale='sqrt', percent=98.),
origin='lower', cmap='Greys_r')
ax2.imshow(segm[y0:y1, x0:x1], origin='lower', cmap='jet')
for aperture in apertures:
aperture.plot(color='red', lw=1.5, alpha=0.5, ax=ax1)
aperture.plot(color='white', lw=1.5, alpha=1.0, ax=ax2)
# read the F105W XDF data
path = '/Users/lbradley/Box Sync/workshop/data/glue'
fn = 'hlsp_xdf_hst_wfc3ir-60mas_hudf_f105w_v1_sci.fits'
data_f105w = fits.getdata(op.join(path, fn))
# extract and measure sources in the F105W data
props_f105w = segment_properties(data_f105w, segm, wcs=wcs, labels=labels)
tbl_f105w = properties_table(props_f105w, columns=columns)
tbl_f105w
id | xcentroid | ycentroid | segment_sum | area |
---|---|---|---|---|
pix | pix | pix2 | ||
int64 | float64 | float64 | float32 | float64 |
1 | 2465.84848452 | 399.740445382 | 22.7464 | 494.0 |
5 | 2633.99658613 | 492.0608717 | 15.9913 | 111.0 |
20 | 2327.88316059 | 717.51682841 | 5.6294 | 143.0 |
50 | 2654.19964617 | 964.971505227 | 1.15134 | 37.0 |
75 | 2898.56557568 | 1116.58094602 | 0.0895663 | 7.0 |
80 | 1755.71360509 | 1151.92670176 | 2.12652 | 45.0 |
# combine the F160W and F105W tables
from astropy.table import hstack
phot = hstack([tbl, tbl_f105w], table_names=['f160w', 'f105w'])
phot
id_f160w | xcentroid_f160w | ycentroid_f160w | ra_icrs_centroid | dec_icrs_centroid | segment_sum_f160w | segment_sum_err | background_sum | background_mean | background_atcentroid | xmin | xmax | ymin | ymax | min_value | max_value | minval_xpos | minval_ypos | maxval_xpos | maxval_ypos | area_f160w | equivalent_radius | perimeter | semimajor_axis_sigma | semiminor_axis_sigma | eccentricity | orientation | ellipticity | elongation | covar_sigx2 | covar_sigxy | covar_sigy2 | cxx | cxy | cyy | id_f105w | xcentroid_f105w | ycentroid_f105w | segment_sum_f105w | area_f105w |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
pix | pix | deg | deg | pix | pix | pix | pix | pix | pix | pix | pix | pix2 | pix | pix | pix | pix | rad | pix2 | pix2 | pix2 | 1 / pix2 | 1 / pix2 | 1 / pix2 | pix | pix | pix2 | |||||||||||||
int64 | float64 | float64 | float64 | float64 | float64 | float64 | object | object | object | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | int64 | float64 | float64 | float32 | float64 |
1 | 2465.87081094 | 399.726636976 | 53.1655214682 | -27.8285090544 | 19.3759366767 | 0.117206007737 | None | None | None | 2459.0 | 2472.0 | 372.0 | 426.0 | 0.0140414275229 | 0.16770106554 | 2460.0 | 402.0 | 2466.0 | 400.0 | 494.0 | 12.5397401797 | 122.183766184 | 11.4193004726 | 2.39595712088 | 0.977740791852 | -1.4520215309 | 0.790183547002 | 4.76607046641 | 7.4909884669 | -14.667582062 | 128.650045341 | 0.171859229177 | 0.039187850115 | 0.0100069572481 | 1 | 2465.84848452 | 399.740445382 | 22.7464 | 494.0 |
5 | 2634.03819015 | 492.064461678 | 53.1623515826 | -27.8269710321 | 14.0897114333 | 0.0594858061938 | None | None | None | 2629.0 | 2639.0 | 487.0 | 498.0 | 0.0124509194866 | 0.937757313251 | 2631.0 | 487.0 | 2634.0 | 492.0 | 111.0 | 5.94410610323 | 36.1421356237 | 1.9692171003 | 1.92742100098 | 0.204936499671 | 1.43341001729 | 0.0212247290096 | 1.02168498698 | 3.7180064855 | 0.0220948262943 | 3.87476121764 | 0.268970452714 | -0.00306746924377 | 0.258089165094 | 5 | 2633.99658613 | 492.0608717 | 15.9913 | 111.0 |
20 | 2327.89478561 | 717.443551847 | 53.1681195861 | -27.8232129488 | 5.09327924345 | 0.0634265023838 | None | None | None | 2323.0 | 2333.0 | 709.0 | 725.0 | 0.0142875546589 | 0.101232789457 | 2332.0 | 722.0 | 2328.0 | 717.0 | 143.0 | 6.74672614861 | 43.7989898732 | 3.40785379096 | 2.40184286615 | 0.709409701436 | -1.38884042562 | 0.295203663809 | 1.41884960044 | 5.9602263893 | -1.04014501329 | 11.4220902249 | 0.170488266907 | 0.0310508001875 | 0.0889634600566 | 20 | 2327.88316059 | 717.51682841 | 5.6294 | 143.0 |
50 | 2654.19478477 | 964.965942774 | 53.161968865 | -27.8190894513 | 0.806716094725 | 0.0301505344762 | None | None | None | 2650.0 | 2658.0 | 962.0 | 968.0 | 0.0137104261667 | 0.0319027081132 | 2652.0 | 966.0 | 2654.0 | 965.0 | 37.0 | 3.43183125879 | 23.1066017178 | 2.48623246998 | 1.03722743727 | 0.908820008767 | 0.655381022095 | 0.582811563362 | 2.39699836376 | 4.28494470958 | 2.4669348987 | 2.97224794181 | 0.446944476139 | -0.741918539471 | 0.644338024956 | 50 | 2654.19964617 | 964.971505227 | 1.15134 | 37.0 |
75 | 2898.58539606 | 1116.73027789 | 53.1573626114 | -27.8165612716 | 0.122965287417 | 0.0132599867602 | None | None | None | 2898.0 | 2899.0 | 1115.0 | 1118.0 | 0.0160002671182 | 0.0187656655908 | 2898.0 | 1116.0 | 2899.0 | 1117.0 | 7.0 | 1.49270533036 | 7.20710678119 | 1.03257280072 | 0.473051042353 | 0.88888596853 | -1.42032587794 | 0.541871486426 | 2.18279362749 | 0.242707512418 | -0.124856004144 | 1.04727636505 | 4.38938760515 | 1.04660319913 | 1.01724566909 | 75 | 2898.56557568 | 1116.58094602 | 0.0895663 | 7.0 |
80 | 1755.70949679 | 1151.83352703 | 53.1788988607 | -27.8159691752 | 1.43870618287 | 0.034281153463 | None | None | None | 1753.0 | 1759.0 | 1148.0 | 1155.0 | 0.0126629555598 | 0.0654352232814 | 1759.0 | 1153.0 | 1756.0 | 1152.0 | 45.0 | 3.78469878303 | 22.4852813742 | 1.82258252039 | 1.49005603819 | 0.575854738034 | 1.36037865658 | 0.182447970656 | 1.22316374262 | 2.26832278138 | 0.225002250232 | 3.2737512592 | 0.443880506986 | -0.0610150894178 | 0.307556740415 | 80 | 1755.71360509 | 1151.92670176 | 2.12652 | 45.0 |
See http://photutils.readthedocs.org/en/latest/photutils/segmentation.html
CHECKIMAGE_TYPE SEGMENTATION
CHECKIMAGE_NAME segmentation.fits
# load a star image
# https://github.com/astropy/photutils-datasets/blob/master/data/M6707HH.fits.gz
from photutils.datasets import load_star_image
hdu = load_star_image()
star_data = hdu.data[0:400, 0:400]
# estimate background and background rms from
# sigma-clipped statistics
from astropy.stats import sigma_clipped_stats
mean, median, std = sigma_clipped_stats(star_data, sigma=3.0)
# find stars using daofind
from photutils import daofind
sources = daofind(star_data - median, fwhm=3.0, threshold=5.*std)
sources
id | xcentroid | ycentroid | sharpness | roundness1 | roundness2 | npix | sky | peak | flux | mag |
---|---|---|---|---|---|---|---|---|---|---|
int64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
1 | 144.247567164 | 6.37979042704 | 0.581562572827 | 0.203512439213 | -0.00821880289245 | 25.0 | 0.0 | 6903.0 | 5.71451872131 | -1.89244914961 |
2 | 208.669068628 | 6.82058053777 | 0.483489662106 | -0.125851384513 | -0.0308111334993 | 25.0 | 0.0 | 7896.0 | 6.73850099891 | -2.07140824271 |
3 | 216.926136655 | 6.5775933198 | 0.693595253797 | -0.706646323127 | -0.0960886879515 | 25.0 | 0.0 | 2195.0 | 1.67120235604 | -0.557572598088 |
4 | 351.625190383 | 8.5459013233 | 0.485778336726 | -0.341513600178 | 0.0151249423679 | 25.0 | 0.0 | 6977.0 | 5.91447184003 | -1.92978992093 |
5 | 377.519909958 | 12.0655009987 | 0.520384876864 | 0.369715621436 | -0.0650824800673 | 25.0 | 0.0 | 1260.0 | 1.1211298456 | -0.124139785194 |
6 | 294.272840467 | 12.7371912508 | 0.680218915989 | 0.106312281745 | -0.344500901996 | 25.0 | 0.0 | 2059.0 | 1.48533945096 | -0.429564290413 |
7 | 85.2177283194 | 14.7103468345 | 0.623941337614 | -0.927179169514 | -0.274573940241 | 25.0 | 0.0 | 1458.0 | 1.24939080492 | -0.241745763414 |
8 | 137.941285661 | 17.4647184039 | 0.543405484171 | 0.0994376959054 | 0.00363728306813 | 25.0 | 0.0 | 5451.0 | 4.69715207557 | -1.67958655295 |
9 | 130.327972445 | 18.6780874278 | 0.627239868433 | 0.179234630855 | 0.00847608979707 | 25.0 | 0.0 | 3432.0 | 2.87258936961 | -1.14568387245 |
10 | 143.110499305 | 18.6235530714 | 0.521248764726 | 0.0759815635531 | -0.0342124559976 | 25.0 | 0.0 | 6722.0 | 5.77264892976 | -1.90343786547 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
276 | 92.3796198909 | 393.24931533 | 0.651596363908 | 0.336953560109 | -0.0029558463632 | 25.0 | 0.0 | 4580.0 | 3.64270928319 | -1.40356128128 |
277 | 256.498485334 | 394.563177212 | 0.356277116679 | -0.137873541286 | -0.157120921977 | 25.0 | 0.0 | 9381.0 | 4.67256968497 | -1.67388946752 |
278 | 302.411031757 | 395.275659094 | 0.605412961527 | -0.294933186934 | -0.108124784092 | 25.0 | 0.0 | 7257.0 | 5.26077425405 | -1.80262416529 |
279 | 351.479535039 | 394.67972307 | 0.641178433577 | 0.786780732556 | -0.231336085872 | 25.0 | 0.0 | 1545.0 | 1.06722678354 | -0.0706417898356 |
280 | 345.593064896 | 395.382215765 | 0.384077996336 | 0.0859556289744 | -0.0647132347379 | 25.0 | 0.0 | 9350.0 | 5.06751489395 | -1.76198758451 |
281 | 268.049236979 | 397.925371446 | 0.296507149185 | -0.888760213866 | -0.808085188561 | 25.0 | 0.0 | 9299.0 | 6.23450523045 | -1.98700498407 |
282 | 268.475068392 | 398.020998272 | 0.283257406393 | 0.334132805571 | -0.970759964067 | 25.0 | 0.0 | 8754.0 | 6.06468200209 | -1.95702008458 |
283 | 299.80943822 | 398.027911813 | 0.320113388974 | -0.0815605235252 | -0.906802614351 | 25.0 | 0.0 | 8890.0 | 6.13258007484 | -1.96910806905 |
284 | 315.689448343 | 398.70251891 | 0.295021382961 | -0.547103811972 | -0.653715309005 | 25.0 | 0.0 | 6485.0 | 5.56746265532 | -1.86414328153 |
285 | 360.437243037 | 398.698539555 | 0.811471438769 | -0.656797339623 | -0.969180784498 | 25.0 | 0.0 | 8079.0 | 5.27758085314 | -1.80608723868 |
# plot the results
from photutils import CircularAperture
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(star_data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
# lets add a large background gradient
ny, nx = star_data.shape
y, x = np.mgrid[:ny, :nx]
gradient = x * y * 1.e-1
data2 = star_data + gradient
plt.imshow(scale_image(data2, scale='sqrt', percent=99.), origin='lower', cmap='Greys')
<matplotlib.image.AxesImage at 0x18974f450>
from photutils.background import Background
bkg = Background(data2, (10, 10), filter_shape=(3, 3), method='median')
# plot the estimated background
plt.imshow(bkg.background, origin='lower', cmap='Greys')
<matplotlib.image.AxesImage at 0x15978f490>
# plot the background-subtracted data
plt.imshow(scale_image(data2 - bkg.background, scale='sqrt', percent=99.), origin='lower', cmap='Greys')
<matplotlib.image.AxesImage at 0x159a1d6d0>
threshold = median + 2.*std
segm_stars = detect_sources(data2 - bkg.background, threshold, npixels=5)
plt.imshow(segm_stars)
<matplotlib.image.AxesImage at 0x15a1058d0>
# create a source mask from the segmentation image
mask = (segm_stars > 0)
bkg2 = Background(data2, (10, 10), filter_shape=(3, 3),
method='sextractor', mask=mask)
plt.imshow(scale_image(data2 - bkg2.background, scale='sqrt', percent=99.),
origin='lower', cmap='Greys')
<matplotlib.image.AxesImage at 0x15ad2e710>
y0 = 200
data2_sub = data2 - bkg2.background
plt.plot(star_data[y0, :])
plt.plot(data2[y0, :])
plt.plot(data2_sub[y0, :], color='cyan')
plt.plot(star_data[y0, :] - median)
plt.axhline(0.0)
<matplotlib.lines.Line2D at 0x15b052990>
Experimental development code from agile sprints
Hope to move functionality to the astropy core or affiliated packages (or possibly upstream to numpy, scipy, or scikit-image)
(very preliminary) Documentation: http://imutils.readthedocs.org/en/latest/