In the previous lab, you extracted aperture photometry for all of the stars in a given image, and experimented with different parameters to understand how identifying stars and measuring their fluxes can change based on how you choose various parameters. This week, we want to use that information to do the following:
Calculate the instrumental magnitudes of the standard stars and M52 stars, based on your flux measurements.
Determine the zeropoint: the difference between the instrumental magnitude of a standard star and its known, calibated magnitude value in a given photometric system.
Use the zeropoint to calibrate the photometry of the cluster stars.
Estimate the errors on our measured fluxes, and propagate these errors into our final magnitude calculations.
Plot a V vs. V-R color-magnitude diagram for M52!
In this lab, you will work with the aligned V-band and R-band images of M52 that you created for Homework 10. Start by downloading your two final aligned M52 FITS files into this directory.
We've provided reduced, aligned frames of the standard star for the M52 dataset, SA 114637.
The first step is to identify the standard star in the provided reduced images. Use iObserve to look up the star and get a DSS image roughly the same size as the Smith detector. You may have to flip the image N-S or E-W to match the orientation of the Smith images.
(1) Include a screenshot of the DSS image from iObserve below.
(2) Decide which of the stars in the Smith images (either V or R band -- they are aligned to each other) is the standard star. Check your answer with one other neighboring group. Place a circle region around the star in the Smith image and put a screenshot below.
(3) What are the x-y coordinates of the standard star in the Smith image?
Answer:
In [ ]:
# The standard fare, plus a few extra packages:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import astropy.io.fits as fits
import os.path
%matplotlib inline
# Newer packages:
from astropy.stats import mad_std
from astropy.stats import sigma_clip
from photutils.utils import calc_total_error
import astropy.stats as stat
from photutils import aperture_photometry, CircularAperture, CircularAnnulus, DAOStarFinder
We would like our photometry module to have two functions:
(1) A function to extract the stars in a given image,
and
(2) A function to perform photometry on that image, with adjustable parameters, which provides errors on the measurements.
The first step is to use a new function (provided below) to estimate the background level. Comment this function and add a docstring, and we will apply it shortly.
In [ ]:
def bg_error_estimate(fitsfile):
"""
Edit this docstring
"""
fitsdata = fits.getdata(fitsfile)
hdr = fits.getheader(fitsfile)
# What is happening in the next step? Read the docstring for sigma_clip.
# Answer:
filtered_data = sigma_clip(fitsdata, sigma=3.,copy=False)
# Summarize the following steps:
#
#
#
bkg_values_nan = filtered_data.filled(fill_value=np.nan)
bkg_error = np.sqrt(bkg_values_nan)
bkg_error[np.isnan(bkg_error)] = np.nanmedian(bkg_error)
print("Writing the background-only error image: ", fitsfile.split('.')[0]+"_bgerror.fit")
fits.writeto(fitsfile.split('.')[0]+"_bgerror.fit", bkg_error, hdr, overwrite=True)
effective_gain = 1.4 # electrons per ADU
error_image = calc_total_error(fitsdata, bkg_error, effective_gain)
print("Writing the total error image: ", fitsfile.split('.')[0]+"_error.fit")
fits.writeto(fitsfile.split('.')[0]+"_error.fit", error_image, hdr, overwrite=True)
return error_image
(1): What is the purpose of the calc_total_error step? How does it relate to the CCD equation? (The docstring and the section on Error Estimation at https://photutils.readthedocs.io/en/stable/aperture.html are helpful resources here.)
Answer:
Part (a): Creating functions for a photometry module
In the following two cells are incomplete functions that you will edit and run. You may need to refer to Lab 9 to remember how to complete various steps. Anywhre you see asterisks are places you will need to edit. Comment as necessary -- some places have prompts for you to comment specifically.
In [ ]:
# Star extraction function -- this function can be to also return the x and y positions to the notebook to use later:
# target_filter_xpos, target_filter_ypos = starExtractor("image.fit", nsigma_value=#, fwhm_value=#)
def starExtractor(fitsfile, nsigma_value, fwhm_value):
"""
This is an incomplete function! Asterisks denote a step where you need to complete the code.
Also, replace this with your docstring, including how to use the function.
"""
# First, check if the region file exists yet, so it doesn't get overwritten
regionfile = fitsfile.split(".")[0] + ".reg"
if os.path.exists(regionfile) == True:
print(regionfile, "already exists in this directory. Rename or remove the .reg file and run again.")
return
# *** Read in the data from the fits file ***
image =
# *** Measure the median absolute standard deviation of the image: ***
bkg_sigma =
# *** Define the parameters for DAOStarFinder ***
daofind =
# Apply DAOStarFinder to the image
sources = daofind(image)
nstars = len(sources)
print("Number of stars found in ",fitsfile,":", nstars)
# Define arrays of x-position and y-position
xpos = np.array(sources['xcentroid'])
ypos = np.array(sources['ycentroid'])
# Write the positions to a .reg file based on the input file name
if os.path.exists(regionfile) == False:
f = open(regionfile, 'w')
for i in range(0,len(xpos)):
f.write('circle '+str(xpos[i])+' '+str(ypos[i])+' '+str(fwhm_value)+'\n')
f.close()
print("Wrote ", regionfile)
return xpos, ypos # Return the x and y positions of each star as variables
In [ ]:
# Photometry function, which returns a table of photometry values for a list of stars
# This function can be used as
# target_phot_table = measurePhotometry(file, star_xpos, star_ypos, aperture_radius, sky_inner, sky_outer, error_array)
def measurePhotometry(fitsfile, star_xpos, star_ypos, aperture_radius, sky_inner, sky_outer, error_array):
"""
Add a docstring here. Add comments at the # below
"""
# *** Read in the data from the fits file:
image =
starapertures = CircularAperture((star_xpos, star_ypos),r = aperture_radius)
skyannuli = CircularAnnulus((star_xpos, star_ypos), r_in = sky_inner, r_out = sky_outer)
phot_apers = [starapertures, skyannuli]
# What is new about the way we're calling aperture_photometry?
# *** Add descriptive comments here
phot_table = aperture_photometry(image, phot_apers, error=error_array)
# Calculate mean background in annulus and subtract from aperture flux
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['bg_subtracted_star_counts'] = final_sum
# *** Add descriptive comments here.
bkg_mean_err = phot_table['aperture_sum_err_1'] / skyannuli.area()
bkg_sum_err = bkg_mean_err * starapertures.area()
# *** Add descriptive comments here.
phot_table['bg_sub_star_cts_err'] = np.sqrt((phot_table['aperture_sum_err_0']**2)+(bkg_sum_err**2))
return phot_table
Part (b): Extracting photometry for the standard star images
(1) Take a look at the standard star images. The data for SA 114637 are not the highest quality. When examining the reduced images, also think back to the conditions during our first observing night at Smith. What factors might be affecting image quality and how?
Answer:
In the cells below, we will complete the following steps for the reduced standard star V-band image:
In [ ]:
# Measure the background of the image
std_V_bgerror = bg_error_estimate("std114637_V_stack.fit")
(2): Open the resulting image, std114637_Vstack_error.fit, in ds9. What do you notice about the image? Describe quantitatively.
Answer:
In [ ]:
# Extract the star positions. Replace ?? with values for the extraction parameters that capture the stars of interest.
std_V_xpos, std_V_ypos = starExtractor("std114637_V_stack.fit", nsigma_value=??, fwhm_value=??)
(3): Open the standard star image in ds9 and load the regions it created. Was your standard star identified correctly? If not, why/what did you have to change? (Pause and check at this point.)
Answer:
In [ ]:
# Measure photometry for the V band image. Replace ?? with reasonable values
std_V_phottable = measurePhotometry("std114637_V_stack.fit", star_xpos=std_V_xpos, star_ypos=std_V_ypos, \
aperture_radius=??, sky_inner=??, sky_outer=??, error_array=std_V_bgerror)
(4): Check the resulting photometry table below. What index is your standard star? (Refer back to the earlier portion of the lab where you measured the standard star position in the image)
Answer:
In [ ]:
# Print a single row from the array with the standard star only
std_V_phottable[??]
In [ ]:
# Measure the background of the R image
In [ ]:
# Measure photometry for the R band image.
# NOTE: Use std_V_xpos and std_V_ypos to extract photometry for the same stars in the same locations & order!
In [ ]:
# Print a single row from the array with the standard star only
std_V_phottable[??]
We can now stitch together the results into a streamlined pandas dataframe, by defining the dataframe index and column labels of interest:
In [ ]:
columns = ['id','xcenter', 'ycenter','Vflux','Vfluxerr','Rflux','Rfluxerr']
In [ ]:
std_fluxtable = pd.DataFrame(
{'id' : std_V_phottable['id'],
'xcenter' : std_V_phottable['xcenter'],
'ycenter' : std_V_phottable['ycenter'],
'Vflux' : std_V_phottable['bg_subtracted_star_counts'],
'Vfluxerr': std_V_phottable['bg_sub_star_cts_err'],
'Rflux' : std_R_phottable['bg_subtracted_star_counts'],
'Rfluxerr': std_R_phottable['bg_sub_star_cts_err']}, columns=columns)
In [ ]:
# Below, check the dataframe to ensure that the combination worked:
std_fluxtable.head()
Part (c): Extracting photometry for the M52 cluster stars
In the cells below, you will repeat the process with the actual cluster images you made for M52 in both V and R bands. Adjust the function inputs (e.g., star extraction) and the photometry parameters accordingly, and we will also estimate the combined errors in the final image and save these results as FITS files.
In the last cell, create a single pandas flux dataframe for the M52 cluster image (M52_fluxtable, just like the std_fluxtable). We will use these to calibrate the photometry in the next section.
In [ ]:
# M52 V-band analysis -- add any blank cells as needed
# Measure the background of the image
M52_V_bgerror =
In [ ]:
# Extract the star positions and save them to new variables.
# Do you have too few stars? Too many? Check the quality of the extraction.
M52_V_xpos, M52_V_ypos =
Once you have decided on extraction parameters, include a screenshot here of the ds9 image of M52 with the regions overlaid.
[screenshot]
In [ ]:
# Measure photometry
M52_V_phottable =
In [ ]:
# Check the M52 photometry table
M52_V_phottable
In [ ]:
# Follow the same approach for the the R band images, using the M52 V band star locations.
# Add blank cells below as needed.
M52_R_bgerror =
M52_R_phottable =
In [ ]:
# Finally, combine the photometry into a single pandas dataframe for M52.
M52_fluxtable =
In [ ]:
# Check the contents of the new dataframe here
Note that when calculating zeropoint, units are in flux per second, so it's critically important than we divide by exposure time of the image.
(1): What would happen if we simply compared the measured fluxes without performing this scaling for exposure time first? By what factor would our estimates of the cluster star magnitudes be incorrect?
Answer:
Recall the relation for instrumental magnitude:
$m_{inst} = -2.5 log_{10}(\textit{flux in 1-second})$
Part (a): Update the pandas dataframes for both the standard star and M52 images, by retrieving the exposure times from the FITS headers and adding new columns, "Vflux_1sec" and "Rflux_1sec", and "Vflux_1sec_err" and "Rflux_1sec_err":
In [ ]:
# Your code here for the standard star dataframe
In [ ]:
# Check that new columns were added
std_fluxtable.head()
In [ ]:
# Your code here for the M52 dataframe
In [ ]:
# Check that new columns were added
M52_fluxtable.head()
Part (b): Calculate the instrumental magnitudes (label them "V_inst" and "R_inst"), and also add these as new columns to your ever-growing standard star and M52 dataframes.
In [ ]:
# Your code here
Part (c): Propagate errors on the fluxes into errors on the instrumental magnitudes:
We discussed in lecture that the uncertainty propagation for log values goes as follows:
If $x = k * log(a)$, where is $k$ is a constant value, then we can evaluate the uncertainty on $x$, that is $\sigma_{x}$ as:
$\sigma_{x} = k * 0.434 \frac{\sigma_{a}}{a}$
Add these to your two pandas dataframes as new columns,
"Vinst_err", "Rinst_err"
In [ ]:
# Your code here
In [ ]:
# Check the tables again in the following cells
In [ ]:
In [ ]:
# Print the row in the standard star table corresponding to the standard star
Part(d): Look up the actual magnitude values for SA 114637 using the Simbad Astronomical Database:
http://simbad.u-strasbg.fr/simbad/sim-fid
(2) What are the V and R band magnitudes for this standard star, and what are their uncertainties? (bracketed values)
Answer:
Now we can use these values for the standard star to determine the zeropoints, as follows:
$m_{calibrated} = m_{inst} + zeropoint$
So to determine the zeropoints, we will need to propagate uncertainties. (Writing zeropoint as magzp):
$magzp \pm \sigma_{magzp}$ = ($m_{calib} \pm \sigma_{m_{calib}}$) - ($m_{inst} \pm \sigma_{m_{inst}}$ )
(3) How do you calculate the value of $\sigma_{magzp}$?
Answer:
First, you'll need to locate the row coresponding to the standard star in the standard star image dataframe.
(4) What are the V and R band instrumental magnitudes and uncertainties of the standard star?
Answer:
Finally, in the cell below calculate the zeropoints for each band. Before proceeding to 10.4, check your method of calculating the uncertainty on each zeropoint from the errors on the insturmental and calibrated magnitudes. After you've estimated the zeropoints, pause here and check your values with us and other groups.
In [ ]:
magzp_V =
magzp_V_error =
magzp_R =
magzp_R_error =
In [ ]:
print("Zeropoint in V: ", magzp_V, "+/-", magzp_V_error)
print("Zeropoint in R: ", magzp_R, "+/-", magzp_R_error)
Now that we have the zeropoint values in hand, we can efficiently calibrate all of the fluxes in the cluster dataframe. Write a quick function below that takes in a zeropoint value, its error, and its corresponding filter, then operates on the series in the panda dataframe and adds new columns for the final calibrated magnitude columns and their uncertainties: Vmag, Vmag_err, Rmag, and Rmag_err.
In [ ]:
# Define your function here
def zpcalc(magzp, magzp_err, filtername, dataframe):
Apply your function to finalize your M52 data table, and show a portion of the dataframe below:
In [ ]:
# Apply function to dataframes
zpcalc(magzp_V, magzp_V_error, "V", M52_fluxtable)
zpcalc(magzp_R, magzp_R_error, "R", M52_fluxtable)
In [ ]:
# Look at dataframe snippet
M52_fluxtable.head()
The final step for our photometry table is to calculate a color, namely V-R. Add this quantity as a final column with its error calculated from the relative errors on V and R. Then we will write the table to a file and save it!
In [ ]:
# Add V-R column and V-R error column
In [ ]:
# Finally, save both calibrated dataframes (standard and M52) here as .csv files;
# These can later be read into Excel, Google Sheets, back into pandas, etc. for future use
std_fluxtable.to_csv('SA114637_photometry.csv')
M52_fluxtable.to_csv('M52_photometry.csv')
Congratulations! You have converted your raw FITS datasets into usable science products, and have measured and calibrated the photometry to a magnitude system that is both universally recognized by other astronomers and includes reasonable error estimates.
One thing we have not done is compare our estimated calibrated photometry for M52 with the known literature values -- this is something you will explore shortly.
Let's see what the data look like!
To plot data with errorbars, we'll use a slightly different plotting method in matplotlib, called plt.scatter. It can be used as follows:
plt.errorbar(x_data, y_data, xerr = x_errors, yerr = y_errors, marker = 'o', linestyle='None')
In the cell below, plot the V vs. V-R color-magnitude digram for M52, adjusting the axes as necessary, and adding all relevant labels, etc.:
In [ ]:
# your CMD here:
plt.figure(figsize=(10,10))
1) What do you notice about the uncertainties in the dataset?
2) How do the median uncertainties compare for the V-band data vs. the R-band data?
3) How does this CMD compare to those you created back in Lab 2?
4) Are any parts of the H-R diagram missing? Which types of stars appear to be present? Which might be absent?
5) What are some potential sources of contamination in our datasets?