In [6]:
# The standard fare:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
from astropy.io import fits
import astropy.stats as stat
from photutils import aperture_photometry, CircularAperture, CircularAnnulus, DAOStarFinder
The goal of this lab is to teach you how to extract aperture photometry for all of the stars in a given image. Although the extraction itself is automated, there is some art to setting the parameters for extraction, and this lab is designed to teach you to choose these in an intelligent, data-driven way.
To start we need to read in a raw image to work with.
In [57]:
#as an example, let's load in a randomly-selected raw Smith R band image
image = fits.getdata('M52-007_R.fit')
We will be using a python package called photutils, which is based on an old IRAF package of the same name. One of the key functions within this package is the DAOStarFinder function, which we'll just call "the star finder" here.
The star finder will extract sources, defined as some multiple (which you provide) of the background/sky level, so first we need to get a reasonable estimate of the background level. This is done using the function mad_std from the astropy.stats package, as below.
In [25]:
from astropy.stats import mad_std
bkg_sigma = mad_std(image)
print(bkg_sigma)
In [58]:
print(np.std(image))
Your answers go here
The star finder requires two parameters to extract stars in the image:
1) The threshhold, which we will define as some number of background levels, above which we will call something a star
2) An estimate for the FWHM of point-sources (stars) in the image
To start, estimate the FWHM of the stars in your image using pyraf's imexam functions, as you did in Lab 8. Measure the FWHM for at least 10 stars and average them. In the cell below, paste the imexam output and calculate the average of the measurements in the cell below that. Insert your calculated average FWHM in the third cell below.
insert imexam output here
In [30]:
#insert calculation of average FWHM here
In [32]:
#FWHM= this is a placeholder. INSERT YOUR VALUE IN PLACE OF 10 BELOW
FWHM=10
We also need to set the threshhold (described above) for star finder, which we define as some multiple (nsigma) of the background level. To start, let's set nsigma to 10, meaning that in order for somehting to be considered a star, it must have at least 10x the detector counts of the background level.
The next several lines below set up the parameters for the star finder (by specifying the FWHM and threshhold parameters), apply the star finder to the image, and then extract and print the number of sources.
In [59]:
nsigma=10
daofind = DAOStarFinder(fwhm=FWHM, threshold=nsigma*bkg_sigma)
sources = daofind(image)
nstars=len(sources)
nstars
Out[59]:
To check how well we're doing here, we need to be able to see how the sources that were automatically extracted line up with apparent sources in the image. To do so, we are going to write the information that star finder extracted from the sources it found into a DS9 region file, so that we can load it with the image. DS9 region files are text files where each line contains the following basic info:
regiontype xcen ycen FWHM
The code below writes the relevant outpt from daofind into a text file with this format.
In [36]:
xpos = np.array(sources['xcentroid'])
ypos = np.array(sources['ycentroid'])
f = open('M52_R.reg', 'w') #you will need to change the first input if you want to write to a different filename later
for i in range(0,len(xpos)):
f.write('circle '+str(xpos[i])+' '+str(ypos[i])+' '+str(FWHM)+'\n')
f.close()
To display this region file, you should open the science image in DS9, then click Region --> Load Regions and navigate to the .reg file above. When you load it, you will see green circles appear on top of all of the stars that the Star Finder has extracted. Place a screenshot of this overlay in the cell below.
DS9 screenshot goes here
Using the median-combined V and R images that you generated for Homework 9, answer the following questions. Include code, screenshots, etc. to support your argument, and add cells below as needed to do calculations, generate new region files, etc.
1) How many sources can you extract at V and R for nsigma=10 from your median-combined images? How much does the number of sources vary between the two wavelengths and why, do you think?
2) How different are the number of extracted stars in the raw R image you used in the example vs. the reduced and median combined R image that you generated for your homework? Name at least one potential reason why they are different, and find an example in the images that shows it.
Hint: An example really means a source that was identified in one image and not the other. Remember you can load multiple images in DS9, and can load a separate region file in each. Zoom in on your discrpant source. To match the zoom level and location between the two images, select Frame --> Match --> Frame --> Physical.
3) For one of your images (V or R), discuss how the number of extracted sources changes when you change the nsigma threshhold. Make an argument based on the images for what you think the most reasonable limit is for this data.
In [60]:
#problem 1 code goes here
In [61]:
#problem 2 code goes here
In [ ]:
#problem 3 code goes here
The next step is to actually extract the photometry, and here too there is some art to choosing parameters. Although photutils will extract the photometry for each star in an automated way, you need to intelligently choose the parameters based on your data.
The tunable parameters are:
We'll start with some potentially reasonable values for these parameters.
In [62]:
aprad = 8
skyin=10
skyout=15
In [55]:
#add comments to this cell
starapertures = CircularAperture((xpos,ypos),r=aprad)
skyannuli = CircularAnnulus((xpos,ypos),r_in=skyin,r_out=skyout)
phot_apers = [starapertures, skyannuli]
phot_table = aperture_photometry(image,phot_apers)
phot_table
Out[55]:
In [56]:
#add comments to this cell
bkg_mean = phot_table['aperture_sum_1']/skyannuli.area()
bkg_starap_sum = bkg_mean * starapertures.area()
final_sum = phot_table['aperture_sum_0']-bkg_starap_sum
phot_table['background subtracted star counts'] = final_sum
phot_table
Out[56]:
Spend the rest of the lab time investigating what changes about the photometric measurement (background subtracted sky counts caluclated column) when you adjust the tunable parameters (aperture radius and inner/outer sky annulus) and report your findings below. You may wish to examine only a handful of stars in the table and avoid pr (to print just a few rows of a table, see the example below), but make sure the rows you select include stars with a range of brightnesses. You may also wish to make separate versions/copies of the table with different aperture parameters so that you can compare without overwriting. Think in particular about crowded fields and see if you can derive the best parameters for this case as well (identify things in the table with very close values for xcenter and ycenter to find these crowded regions).
In [65]:
#example of printing just a few rows in a table - chosen to be one bright star, one faint
phot_table[6:8]
Out[65]: