Photometry with Photutils

Photutils is an affiliated package of Astropy to provide tools for detecting and performing photometry of astronomical sources.

It contains tools for:

  • Background Estimation (photutils.background)
  • Source Detection (photutils.detection)
  • Grouping Algorithms
  • Aperture Photometry (photutils.aperture)
  • PSF Photometry (photutils.psf)
  • PSF Matching (photutils.psf.matching)
  • Image Segmentation (photutils.segmentation)
  • Centroids (photutils.centroids)
  • Morphological Properties (photutils.morphology)
  • Geometry Functions (photutils.geometry)
  • Datasets (photutils.datasets)
  • Utility Functions (photutils.utils)

In [ ]:
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
import seaborn
from astropy.io import fits
from astropy import units as u
from astropy.coordinates import SkyCoord
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 14
plt.rcParams['lines.linewidth'] = 2
plt.rcParams['xtick.labelsize'] = 13
plt.rcParams['ytick.labelsize'] = 13
plt.rcParams['axes.titlesize'] = 14
plt.rcParams['legend.fontsize'] = 13

As we will be using the image visualization several times, let's make a function that makes our code cleaner:


In [ ]:
from astropy.visualization import ZScaleInterval

def my_python_ds9(data):
    norm = ZScaleInterval()
    vmin, vmax = norm.get_limits(data)
    plt.imshow(data, vmin=vmin, vmax=vmax, interpolation='none', origin='lower')
    plt.colorbar()

Before performing any photometry we need to have a first guess of the image background properties. For this section we will use a simpler optical image from the Sloan Digital Sky Survey (SDSS). This image contains the Spindle Galaxy, also known as Messier 102 or NGC 5866.


In [ ]:
sdss_g_hdu_list = fits.open('../resources/sdss_g.fits')
sdss_g_hdu_list.info()

In [ ]:
my_python_ds9(sdss_g_hdu_list[0].data)

In [ ]:
plt.imshow(np.power(sdss_g_hdu_list[0].data[580:730, 500:750], 0.5), cmap='inferno', origin='lower')
plt.grid('off')

Let's get the basic statistics of it. For that we will need to remove the sources using a sigma clipping routine:


In [ ]:
from astropy.stats import sigma_clipped_stats

mean, median, std = sigma_clipped_stats(sdss_g_hdu_list[0].data, sigma=3.0, iters=5)    
print((mean, median, std))

We will need also the pixel scale that we can retrieve from the WCS:


In [ ]:
from astropy import wcs

sdss_g_image_wcs = wcs.WCS(sdss_g_hdu_list[0].header)
sdss_pixelscale = np.mean(wcs.utils.proj_plane_pixel_scales(sdss_g_image_wcs) * u.degree).to('arcsec')
sdss_pixelscale

Object detection

To detect the sources inside a astronomical image PhotUtils provides DAO Phot implementations to detect


In [ ]:
from photutils import DAOStarFinder
sigma_detection = 5.0
daofind = DAOStarFinder(fwhm=1.5 * u.arcsec / sdss_pixelscale, threshold=sigma_detection*std)    
sources = daofind(sdss_g_hdu_list[0].data - median)    
print(sources)

In [ ]:
my_python_ds9(sdss_g_hdu_list[0].data)

plt.scatter(sources['xcentroid'], sources['ycentroid'], alpha=0.5, color='limegreen')

PSF Modelling

We assumed that the image had a typical value of 1.5" per pixel. But we can make a more accurate estimation by fitting the pixels to a moffat profile around a detected star.

First, let's select the source we want to use for PSF modelling:


In [ ]:
isource = sources[50]

# High SN - 50
# Low SN - 33
# Non-uniform background - 78

print ("x pos: " + str(isource['xcentroid']))
print ("y pos: " + str(isource['ycentroid']))

stamp_radius = 50
my_python_ds9(sdss_g_hdu_list[0].data[int(isource['ycentroid'] - stamp_radius):
                                      int(isource['ycentroid'] + stamp_radius), 
                                      int(isource['xcentroid'] - stamp_radius): 
                                      int(isource['xcentroid'] + stamp_radius)])

As we intend to find the profile of the source, we need to remove the possible sky background that lies behind:


In [ ]:
# Median bkg subtracted image
bkg_subtracted_image = sdss_g_hdu_list[0].data - median

# 2D background subtracted - For later
#bkg_subtracted_image = sdss_g_hdu_list[0].data - bkg.background

To simplify the problem, we turn the 2D profile into a 1D distance array from the center of each pixel to the centroid of the source estimated by DAO Phot:


In [ ]:
flux_counts = []
pixel_distance = []

x_cen = int(isource['xcentroid'])
y_cen = int(isource['ycentroid'])

# Pixels around detection loop
analysis_radius = 10
for x in range(x_cen - analysis_radius, x_cen + analysis_radius):
    for y in range(y_cen - analysis_radius, y_cen + analysis_radius):
        flux_counts.append(((bkg_subtracted_image[y][x]) / isource['peak']))
        pixel_distance.append(np.linalg.norm((isource['xcentroid'] - x, isource['ycentroid'] - y)))

Here we present two possible models that can fit the PSF distribution. A Gaussian and a Moffat profile:


In [ ]:
from astropy.modeling import models, fitting

model = 'moffat'

if model == 'gaussian':
    # Fit the data using a Gaussian
    fwhm_best_guess = 1.5
    model_init = models.Gaussian1D(amplitude=1.0, mean=0, stddev=fwhm_best_guess)
elif model == 'moffat':
    # Fit the data using a Moffat
    model_init = models.Moffat1D(amplitude=1.0, x_0=0, gamma=2., alpha=3.5)
else:
    raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

fitter = fitting.SimplexLSQFitter()
fitted_model = fitter(model_init, pixel_distance, flux_counts)

print ("Fit value:",  fitter.fit_info['final_func_val'])
print ("SN:", isource['flux'] * sigma_detection)

Once fitted the models, we need to convert from the parameters to the actual FWHM estimate.


In [ ]:
# FHWM conversion
if model == 'gaussian':
    iFWHM = 2.355 * fitted_model.stddev * sdss_pixelscale.value
elif model == 'moffat':
    iFWHM = fitted_model.gamma * 2 * np.sqrt(2 ** (1. / fitted_model.alpha) - 1) * sdss_pixelscale.value
else:
    raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

print ("FWHM estimated: ", iFWHM)

We can finally plot and see how the model traces the pixel values traces our fitted model.


In [ ]:
# Check fitting
if fitter.fit_info['final_func_val'] < 5.0:
    color = 'green'
else:
    color = 'red'
    
# Plot the data with the best-fit model
plt.figure()
plt.plot(pixel_distance, flux_counts, 'ko')
rx = np.linspace(0, int(max(pixel_distance)), int(max(pixel_distance)) * 10)
plt.plot(rx,
         fitted_model(rx),
         color=color,
         lw=3.0,
         label='SN: %.2f, Fit: %.2f, FWHM: %.2f"' % (isource['flux'] * sigma_detection,
                                                       fitter.fit_info['final_func_val'],
                                                       iFWHM))
plt.xlabel('Distance (pixels)')
plt.ylabel('Normalized flux (ADU)')
plt.title('%s profile fit' % model)
plt.legend()
plt.show()

 Background modelling

As we have seen in the case with non-uniform background, the constant median can be insuficient. Here we produce a 2D model of the background that can be subtracted from the original image to improve the modelling of the stars close to a very large extended source (or when the backrgound is not flat for any other reason)


In [ ]:
from photutils import Background2D, SigmaClip, MedianBackground
sigma_clip = SigmaClip(sigma=3., iters=10)
bkg_estimator = MedianBackground()
bkg = Background2D(data=sdss_g_hdu_list[0].data, 
                   box_size=(100, 100), 
                   filter_size=(3, 3),
                   sigma_clip=sigma_clip, 
                   bkg_estimator=bkg_estimator)
my_python_ds9(bkg.background)

Now let's go back to where the background was subtracted to verify the difference.

Aperture photometry

We will work with the simple circular apertures. First we need to set the size we want and create the aperture objects:


In [ ]:
from photutils import CircularAperture

aperture_radius = 5.0

positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=aperture_radius)

my_python_ds9(sdss_g_hdu_list[0].data[0:400, 0:400])
apertures.plot(color='limegreen', lw=2)

Global Background subtraction

As we only want the flux from the sources we are interested in, we need to remove the contribution from the background. The simplest way is to remove the median value constant of the sigma clipped image that we calculated before to the entire array:


In [ ]:
from photutils import aperture_photometry

bkg_subtracted_image = sdss_g_hdu_list[0].data - median
phot_table_global_bkg = aperture_photometry(bkg_subtracted_image, apertures)
print (phot_table_global_bkg)

Local Background subtraction

We could also remove the 2D background model that we calculated before, but it is usually more precise to create an annulus around the source we are interested in, and estimate the background level there.


In [ ]:
from photutils import CircularAnnulus

r_in = 10
r_out = 15

annulus_apertures = CircularAnnulus(positions, r_in=r_in, r_out=r_out)

my_python_ds9(sdss_g_hdu_list[0].data[0:400, 0:400])
apertures.plot(color='limegreen', lw=2)
annulus_apertures.plot(color='purple', lw=2, alpha=1)

The circular apertures and the annulus can be joined for a common photometry processing.


In [ ]:
apers = [apertures, annulus_apertures]
phot_table_local_bkg = aperture_photometry(sdss_g_hdu_list[0].data, apers)
print(phot_table_local_bkg)

Now we use the aperture_sum_1 to estimate the level of background around the source. We need to know the area of the annulus for this estimation:


In [ ]:
bkg_mean = phot_table_local_bkg['aperture_sum_1'] / annulus_apertures.area()

And finally we can remove the background estimation to all pixels in the aperture:


In [ ]:
bkg_sum = bkg_mean * apertures.area()
final_sum = phot_table_local_bkg['aperture_sum_0'] - bkg_sum
phot_table_local_bkg['residual_aperture_sum'] = final_sum
print(phot_table_local_bkg['residual_aperture_sum'])

In this comparison we see that many sources have the same flux with both methods but a significant amount of sources (the ones in the galaxy halo) have significantly more flux in the global subtraction method, as the flux from M102 is added to the one of the stars.


In [ ]:
plt.scatter(phot_table_local_bkg['residual_aperture_sum'], phot_table_global_bkg['aperture_sum'])
plt.xlim(0,100)
plt.ylim(0,100)
plt.xlabel('Local background')
plt.ylabel('Global background')

Let's make a plot to verify this! Astropy qtables work similarly to Pandas, so we can make iloc searches.

Exercise PhotUtils: Identify problematic sources

Locate which sources differ their flux in more than 50% between the measurement with local and global background estimation. Plot the results on the image.

TIP: Using the "pandas-like" tools can facilitate the selection


In [ ]:
# %load -r 43-52 solutions/10_Astronomy.py

Accessing to annulus pixels

By default when calling aperture_photometry, we only receive the sum of pixels and therefore we can only access to the mean of the pixel values inside an aperture. To properly measure the local background, we would need to get the median of the pixels in the annulus and preferably, perform a sigma clipping over the values.

Recently photutils enabled masks so we can actually do that!


In [ ]:
# Flat annulus - 11
# Source in annulus - 12

source_num = 12

x = phot_table_local_bkg[source_num]['xcenter'].value
y = phot_table_local_bkg[source_num]['ycenter'].value

stamp_radius = 20
my_python_ds9(sdss_g_hdu_list[0].data[int(y - stamp_radius):
                                      int(y + stamp_radius), 
                                      int(x - stamp_radius): 
                                      int(x + stamp_radius)])
plt.grid('off')

We can access to the annulus mask:


In [ ]:
# Create annulus mask
annulus_apertures = CircularAnnulus((x, y), r_in=r_in, r_out=r_out)
masks = annulus_apertures.to_mask(method='center')
m0 = masks[0]

plt.imshow(m0)
plt.grid('off')

And apply it to the data:


In [ ]:
# Cut to data
cutout_data = m0.apply(sdss_g_hdu_list[0].data)

plt.imshow(cutout_data)
plt.grid('off')

Now that we have access to the pixels, we can perform the median and compare to the mean


In [ ]:
annulus_array = cutout_data[cutout_data != 0]

# Apply statistics to masked data
mean = np.mean(annulus_array)
median = np.median(annulus_array)
std = np.std(annulus_array)
print (mean, median, std)

And even sigma-clip the sources that may appear on top of the background. This creates a numpy.masked array:


In [ ]:
from astropy.stats import sigma_clip

clip_annulus_array = sigma_clip(annulus_array, sigma=3, iters=2)

mean_clipped = np.ma.mean(clip_annulus_array)
median_clipped = np.ma.median(clip_annulus_array)
std_clipped = np.ma.std(clip_annulus_array)
print (mean_clipped, median_clipped, std_clipped)

Exercise PhotUtils 2: Estimate aperture error

Compute what would be the measurement error for the last aperture (the one with sigma clipping)

$$ \sigma^2_{F} = \sum_i{\sigma^2_{i}+F/g} $$

TIP: Gain is the exposure time when pixels are in e/s


In [ ]:
# %load -r 56-61 solutions/10_Astronomy.py

In this link you can find details on how to estimate the error: http://photutils.readthedocs.io/en/stable/photutils/aperture.html#error-estimation

Matching catalogues

astropy.coordinates contains commonly-used tools for comparing or matching coordinate objects. It supports leverages the coordinate framework to make it straightforward to find the closest coordinates in a catalog to a desired set of other coordinates.

Let's do a new catalogue with the SDSS r band image on the same location in the sky


In [ ]:
sdss_r_hdu_list = fits.open('../resources/sdss_r.fits')
sdss_r_hdu_list.info()

If we plot the previous apertures in the sky, we realize that there is a slight offset (either due to pointing inaccuracy or dithering observarions)


In [ ]:
my_python_ds9(sdss_r_hdu_list[0].data[0:200, 0:300])
apertures.plot(color='limegreen', lw=2)

In order to make a r band catalogue, we follow the same process as before

[ don't replicate code like we do here... you should create a function or Scott will get mad!! ]


In [ ]:
print ("Image statistics...")
mean_r, median_r, std_r = sigma_clipped_stats(sdss_r_hdu_list[0].data, sigma=3.0, iters=5)    

print ("Source detection...")
sources_r = daofind(sdss_r_hdu_list[0].data - median_r)    

print ("Aperture definition...")
positions_r = (sources_r['xcentroid'], sources_r['ycentroid'])
apertures_r = CircularAperture(positions_r, r=aperture_radius)
annulus_apertures_r = CircularAnnulus(positions_r, r_in=10, r_out=15)

print ("Photometry calculation...")
apers_r = [apertures_r, annulus_apertures_r]
phot_table_local_bkg_r = aperture_photometry(sdss_r_hdu_list[0].data, apers_r)
bkg_mean_r = phot_table_local_bkg_r['aperture_sum_1'] / annulus_apertures_r.area()
bkg_sum_r = bkg_mean_r * apertures_r.area()
final_sum_r = phot_table_local_bkg_r['aperture_sum_0'] - bkg_sum_r
phot_table_local_bkg_r['residual_aperture_sum'] = final_sum_r

print ("..ready!")

In addition we se that more/different objects have been detected due to the different spectral emission of the sources


In [ ]:
my_python_ds9(sdss_r_hdu_list[0].data[0:200, 0:300])
apertures.plot(color='limegreen', lw=2)
apertures_r.plot(color='red', lw=2)

First we create the two catalogues in sky coordinates using the WCS information as matching needs to be happen the same frame


In [ ]:
origin = 0
### g band catalogue ###
# Create wcs object
sdss_g_image_wcs = wcs.WCS(sdss_g_hdu_list[0].header)

# Pixels to Sky
lon_g, lat_g = sdss_g_image_wcs.all_pix2world(phot_table_local_bkg['xcenter'], 
                                              phot_table_local_bkg['ycenter'], 
                                              origin)
cat_g = SkyCoord(ra=lon_g*u.degree, dec=lat_g*u.degree)  


### r band catalogue ###
# Create wcs object
sdss_r_image_wcs = wcs.WCS(sdss_r_hdu_list[0].header)

# Pixels to Sky
lon_r, lat_r = sdss_r_image_wcs.all_pix2world(phot_table_local_bkg_r['xcenter'], 
                                          phot_table_local_bkg_r['ycenter'], 
                                          origin)
cat_r = SkyCoord(ra=lon_r*u.degree, dec=lat_r*u.degree)

Closest catalogue matching

The following function performs the matching between the two catalogues, and returns what is the closest "r" source to each of the "g" sources. Therefore the resulting arrays contain the same number of indices than the catalogue "g".


In [ ]:
id_r, d2d, d3d = cat_g.match_to_catalog_sky(cat_r)

print (len(id_r))

hist_data = plt.hist(d2d.to('arcsec'), bins=100)
plt.xlabel('distance (arcsec)')
plt.ylabel('counts')

Comparing the fluxes, some matches show a very large difference, due to source mismatch and not due to color.


In [ ]:
plt.scatter(phot_table_local_bkg['residual_aperture_sum'], phot_table_local_bkg_r[id_r]['residual_aperture_sum'])
plt.xlabel('g flux')
plt.ylabel('r flux')

Search around coordinates matching

The following function performs the matching between the two catalogues, returning only those sources that match within certain tolerance in the sky. Therefore


In [ ]:
id_r_around, id_g_around, d2d_around, d3d_around = cat_g.search_around_sky(cat_r, 1*u.arcsec)  

print (len(id_r_around))

hist_data = plt.hist(d2d_around.to('arcsec'), bins=100)
plt.xlabel('distance (arcsec)')
plt.ylabel('counts')

Now we can see how the biggest outliers have disappeared, remaining only the difference in flux due to color.


In [ ]:
plt.scatter(phot_table_local_bkg[id_g_around]['residual_aperture_sum'], phot_table_local_bkg_r[id_r_around]['residual_aperture_sum'])
plt.xlabel('g flux')
plt.ylabel('r flux')

From fluxes to magnitudes (in SDSS)

The images we use are already calibrated in photometry, and the units used are nanomaggies. Therefore

m = [22.5 mag] – 2.5 log10 f.

http://www.sdss.org/dr12/algorithms/magnitudes/#nmgy


In [ ]:
def to_sdss_mag(fluxes):
    return 22.5 - np.log(fluxes)

Checking relation between measured flux and calibrated magnitude


In [ ]:
mags_g = to_sdss_mag(phot_table_local_bkg['residual_aperture_sum'])
plt.scatter(mags_g, phot_table_local_bkg['residual_aperture_sum'])

Matching to an external catalogue

Load the SDSS catalogue with astroquery that overlaps the image


In [ ]:
from astroquery.sdss import SDSS
from astropy import coordinates as coords

region = np.mean(sdss_r_image_wcs.calc_footprint(), axis=0)
c = SkyCoord(region[0] * u.degree, region[1] * u.degree, frame='icrs')

sdss_query = SDSS.query_region(c,
                               spectro=False,
                               radius=20*u.arcmin,
                               photoobj_fields=['ra','dec','u','g','r','i','z'])
print(sdss_query)

Checking the catalogue: Color - Color plot


In [ ]:
from matplotlib.colors import LogNorm

plt.hist2d(sdss_query['g'] - sdss_query['r'], 
           sdss_query['r'] - sdss_query['i'],
           range=[(-0.5,2), (-0.5,2)], 
           bins=50, 
           norm=LogNorm(),
          cmap='hot')
plt.colorbar()
plt.xlabel('g-r')
plt.ylabel('r-i')

Finally we can check that calibrated photometry matches within some margin (probably due to different extraction method - Aperture vs Model)


In [ ]:
cat_sdss_query = SkyCoord(ra=sdss_query['ra'] * u.degree, dec=sdss_query['dec'] * u.degree)  

id_query_around, id_g_around, d2d_around, d3d_around = cat_g.search_around_sky(cat_sdss_query, 0.5*u.arcsec)  

print (len(id_query_around))

plt.scatter(sdss_query[id_query_around]['g'], to_sdss_mag(phot_table_local_bkg[id_g_around]['residual_aperture_sum']))
plt.xlim(15, 25)
plt.ylim(15,25)