The affiliated package photutils offers now some basic functions to perform PSF/PRF photometry on astronomical image data. In this notebook I'd like to give a short introduction on the basis of a Spitzer data set.
As an improvement, compared to aperture photometry, PSF photometry makes use of instrument related information. To estimate the the physical parameters of an astrophysical point source it uses a model how this point souce is imaged by the instrument. This model is called point spread function (PSF). The PSF can be either obtained from the data itself or modeled by some analytical function. Additionally there is a second concept, the so called point response function (PRF), which denotes the PSF after discretization with a CCD. We are going to deal mostly with PRFs.
The test image data can be downloaded from the Spitzer database. This notebook and the catalog is contained in the Gist PSFPhotometry on GitHub, which can be easily cloned, if you want to run and test the code yourself.
import numpy as np
from astropy.table import Table
from astropy.io import fits
from astropy.wcs import wcs
import astropy.units as u
At first we use astropy.table
to read the coordinates and fluxes for the sources from the catalog:
catalog = Table.read("catalog.xml")
l = catalog['l']
b = catalog['b']
fluxes_catalog = catalog['f4_5']
We read the image data using astropy.io.fits
. As the image data is given in MJy/sr
, but in the catalog the units are mJy
, we have to convert the units. Therfore we use that the image data has a resolution of (1.2 arcsec)^2 / pixel
. Furthermore we set up a mask to exclude NaN values.
data = fits.getdata('GLM_01800+0000_mosaic_I2.fits')
factor = (u.MJy / u.sr * (1.2 * u.arcsec) ** 2 / u.pixel).to(u.mJy / u.pixel)
data *= factor.value
mask = np.isfinite(data)
As the photometry function are currently only able to handle pixel coordinates, we define a WCS transformation to convert the galactic coordinates into pixel coordinates:
header = fits.getheader('GLM_01800+0000_mosaic_I2.fits')
wcs_transform = wcs.WCS(header)
coords = wcs_transform.wcs_world2pix(zip(l, b), 0)
# It seems that the coordinates are shifted. To have the PRF images properly aligned in the plot below, we correct for this.
coords.T[0] += 0.18
coords.T[1] += 0.3
We only consider a slice of the image, because the whole image is very large.
x = coords.T[0]
y = coords.T[1]
x_min, x_max = 3500, 4500
y_min, y_max= 4000, 4500
slice_mask = (x > x_min) * (x < x_max) * (y > y_min) * (y < y_max)
data_slice = data[y_min:y_max, x_min:x_max]
coords_slice = np.array([x[slice_mask] - x_min, y[slice_mask] - y_min]).T
fluxes_catalog_slice = fluxes_catalog[slice_mask]
Here is a plot with the source positions overlaied:
figure(figsize=(16, 12))
imshow(data_slice, cmap='hot', vmin=0, vmax=10, interpolation='None')
plot(coords_slice.T[0], coords_slice.T[1], marker="o", markerfacecolor='None', markeredgecolor='y', linestyle='None')
xlim(0, x_max - x_min)
ylim(0, y_max - y_min)
colorbar(orientation='horizontal')
<matplotlib.colorbar.Colorbar instance at 0xb8ffe6c>
Now we're ready to estimate the PRF from the data.
To estimate the PRF we work on the complete dataset. We choose a size of 7 pixels and a subsampling of 5 / pixel for the PRF.
from photutils.psf import create_prf
prf_discrete = create_prf(data, coords, 7, fluxes=fluxes_catalog, mask=np.logical_not(mask), subsampling=5)
Furthermore we measure the runtime:
#%timeit -n 3 create_prf(data, coords, 7, fluxes=fluxes_catalog, mask=np.logical_not(mask), subsampling=5)
The object returned by create_prf
is a DiscretePRF
object, which is basically a look-up table of the different PRFs at different subpixel positions. But it behaves and can be handled like a usual astropy.modeling.ParametricModel
. The PRFs are stored in the DiscretePRF._prf_array
attribute. We make a plot to illustrate this:
fig, axes = subplots(nrows=5, ncols=5)
fig.set_size_inches(12, 9)
# Plot kernels
for i in range(5):
for j in range(5):
prf_image = prf_discrete._prf_array[i, j]
im = axes[i, j].imshow(prf_image, interpolation='None')
cax = fig.add_axes([0.9, 0.1, 0.03, 0.8])
colorbar(im, cax=cax)
subplots_adjust(left=0.05, right=0.85, top=0.95, bottom=0.05)
show()
We perform photometry on the first 500 sources of the dataset:
from photutils.psf import psf_photometry
fluxes_photutils = psf_photometry(data_slice, coords_slice, prf_discrete)
To compare with the catalog results, we make a scatter plot:
scatter(fluxes_catalog_slice, fluxes_photutils)
loglog()
xlim(40, 500)
ylim(40, 500)
xlabel("Fluxes catalog [mJy]")
ylabel("Fluxes Photutils [mJy]")
<matplotlib.text.Text at 0xd2bd7cc>
We can also perform photometry using a GaussianPSF
:
from photutils.psf import GaussianPSF
psf_gaussian = GaussianPSF(1)
fluxes_gaussian = psf_photometry(data_slice, coords_slice, psf_gaussian)
# Measure runtime
%timeit -n 3 psf_photometry(data_slice, coords_slice, psf_gaussian)
3 loops, best of 3: 1.18 s per loop
Again we compare with the catalog:
scatter(fluxes_catalog_slice, fluxes_gaussian)
loglog()
xlim(40, 500)
ylim(40, 500)
xlabel("Fluxes catalog [mJy]")
ylabel("Fluxes Gaussian Photutils [mJy]")
<matplotlib.text.Text at 0xd4a3ecc>
It is also possible to remove PRFs from the data and make residual images:
from photutils.psf import remove_prf
residuals = remove_prf(data_slice.copy(), prf_discrete, coords_slice, fluxes_photutils)
# Measure runtime
%timeit -n 3 remove_prf(data_slice.copy(), prf_discrete, coords_slice, fluxes_photutils)
3 loops, best of 3: 54.1 ms per loop
We make a plot of the residual image:
figure(figsize=(16, 12))
imshow(residuals, cmap='hot', vmin=-1, vmax=10, interpolation='None')
plot(coords_slice.T[0], coords_slice.T[1], marker="o", markerfacecolor='None', markeredgecolor='y', linestyle='None')
xlim(0, x_max - x_min)
ylim(0, y_max - y_min)
colorbar(orientation='horizontal')
<matplotlib.colorbar.Colorbar instance at 0xdc7366c>