Photutils is an affiliated package of Astropy to provide tools for detecting and performing photometry of astronomical sources.
It contains tools for:
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
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')
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()
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.
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)
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)
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.
In [ ]:
# %load -r 43-52 solutions/10_Astronomy.py
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)
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
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)
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')
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')
In [ ]:
def to_sdss_mag(fluxes):
return 22.5 - np.log(fluxes)
In [ ]:
mags_g = to_sdss_mag(phot_table_local_bkg['residual_aperture_sum'])
plt.scatter(mags_g, phot_table_local_bkg['residual_aperture_sum'])
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)
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)