Note - the data for all the interactive sessions is available for download.

Hands-On Exercise 1:

Constructing a Light Curve From PTF Photometric Files

Version 0.1

The following iPython notebook details a few different methods for constructing light curves from PTF images using photometric measurements produced by SExtractor (Bertin & Arnouts 1996). It details the steps to go from PTF SExtractor catalogs, which are in a fits table format, to an actual light curve. The methods presented here are for illustration only and are not reccommended for use in scientific publications. There are many caveats when using PTF data (see Mansi's talk or the notes throughout this document).


By AA Miller (c) 2014 Jul 24

Collecting Catalog Files

PTF SExtractor catalogs can be obtained from the IRSA PTF image archive. Images and source catalogs can be searched for via position or PTF field ID, more details are available on the PTF Data Access website. The example below uses source catalogs from images of the open cluster M44.

We begin with the assumption that the images have been downloaded and unzipped in the present working directory.

Making a Light Curve: I. The Basic Solution

For PTF a nightly photometric solution is determined for each image (for further details see Ofek et al. 2012), thus, the easiest way to create a light curve for a given star is to grab its measured brightness on each of the PTF images. This is simple to implement in python using a for loop.

In [1]:
import numpy as np                  # array manipulation
import astropy.io.fits as fits      # fits file manipulation
import glob                         # to grab source files

files = glob.glob("PTF*ctlg")       # list of all source files
targ_ra = 128.86346                 # right ascension of target
targ_dec = 19.75875                 # declination of target
match_r = 2.0                       # source matching radius in arcesc

Nims = len(files)
HJD = np.empty(Nims)                # store HJD of observation
mag = np.empty(Nims)                # store the magnitude
magunc = np.empty(Nims)             # store the mag uncertainty

# Now loop over the files to create the light curve
for imnum, cat in enumerate(files):
    HJD[imnum] = fits.getval(cat, "HJD", 2)     # HJD stored in cat header
    Sdat = fits.getdata(cat)
    Sdat_RA = Sdat["ALPHAWIN_J2000"]            # RA list for cat
    Sdat_Dec = Sdat["DELTAWIN_J2000"]           # Dec list for cat
    source_sep = np.hypot((Sdat_RA - targ_ra) * np.cos(np.radians(targ_dec)), 
        Sdat_Dec - targ_dec)

    # cat source is a match if it is closest to targ and < 2" away
    targ_match = ( (source_sep < match_r/3600.) & 
        (source_sep == np.min(source_sep)) )
    
    if sum(targ_match) == 1:
        ZP = Sdat["ZEROPOINT"][targ_match][0]
        mag[imnum] = Sdat["MAG_AUTO"][targ_match][0] + ZP
        magunc[imnum] = Sdat["MAGERR_AUTO"][targ_match][0]
    else:
        mag[imnum] = -99.           # non-detection
        magunc[imnum] = -99.        # non-detection

Recall that the zeropoint determined for PTF images is positionally dependent, and that SExtractor mag measurements are instrumental. We therefore must add the zeropoint, which is different for every source, to the instrumental mag to get a quasi-calibrated magnitude. Also note that the source of interest may not be detected on every image, here we elect to report both the mag and mag uncertainty as -99 to indicate the source was not detected. There are multiple reasons this may be the case, for example: it may be too faint to be detected, or it may not be within the positional limits of the CCD.

Now that we have a light curve, we can plot the results.

In [3]:
targ_det = mag != -99.          # identify epochs where targ is detected

# plot mag vs HJD along with uncertainties
figure(figsize(13,8))
errorbar(HJD[targ_det], mag[targ_det], magunc[targ_det], fmt='ok')
ylim(np.max(mag[targ_det]) + 0.3, np.min(mag[targ_det]) - 0.3)
xlim(np.min(HJD[targ_det]) - 15, np.max(HJD[targ_det]) + 15)
xlabel(r'${\rm HJD}$', fontsize=18)
ylabel(r'$R \; {\rm mag}$', fontsize=18)

print "The LC scatter = %5.3f mag" % np.std(mag[targ_det])
The LC scatter = 0.325 mag

Success!

Or... at least a partial success. From Ofek et al. 2012, we know that a typical star of this brightness should have a root-mean-square (RMS) scatter of $$0.03 mag, and, instead, we see that the scatter is more than 0.3 mag. This means we have either found a variable star, or our calibration method is not perfectly correct.

Important Caveat

An essential caveat to remember when dealing with PTF data is that the photometric calibration routine only works well on nights when the observing conditions are photometric, which is not the case for the majority of nights at Palomar. Thus, without an optimal calibration we find that the scatter in the light curve is larger than we would otherwise hope.

Making a Light Curve: II. Photometric Conditions Only

The easiest way to rectify the aforementioned problems is to identify which images were taken under photometric conditions, and only use those while constructing our light curve. This will require a very small tweak to the python commands that we used above.

We use the IPAC keywords > PHTCALFL - Was phot.-cal. module executed?

PCALRMSE - RMSE from (zeropoint, extinction) data fit

to determine whether or not the night was photometric. First we require that the photometric calibration module was executed, then we require that the RMSE from the fit to the data is less than 0.04 as described in Ofek et al. 2012.

In [4]:
Nims = len(files)
HJD = np.empty(Nims)                # store HJD of observation
mag = np.empty(Nims)                # store the magnitude
magunc = np.empty(Nims)             # store the mag uncertainty

# Now loop over the files to create the light curve
for imnum, cat in enumerate(files):
    HJD[imnum] = fits.getval(cat, "HJD", 2)     # HJD stored in cat header
    # vvv New condition is added here vvv
    if ( fits.getval(cat, "PHTCALFL", 2) == 1 and 
        fits.getval(cat, "PCALRMSE", 2) < 0.04 ):  # phot. conditions?
        Sdat = fits.getdata(cat)
        Sdat_RA = Sdat["ALPHAWIN_J2000"]            # RA list for cat
        Sdat_Dec = Sdat["DELTAWIN_J2000"]           # Dec list for cat
        source_sep = np.hypot( 
        (Sdat_RA - targ_ra) * np.cos(np.radians(targ_dec)), 
        Sdat_Dec - targ_dec)

        # cat source is a match if it is closest to targ and < 2" away
        targ_match = ( (source_sep < match_r/3600.) & 
            (source_sep == np.min(source_sep)) )
    
        if sum(targ_match) == 1:
            ZP = Sdat["ZEROPOINT"][targ_match][0]
            mag[imnum] = Sdat["MAG_AUTO"][targ_match][0] + ZP
            magunc[imnum] = Sdat["MAGERR_AUTO"][targ_match][0]
        else:
            mag[imnum] = -99.           # non-detection
            magunc[imnum] = -99.        # non-detection
    else:
        mag[imnum] = -99.           # non-detection
        magunc[imnum] = -99.        # non-detection

Note that in addition to the examples listed earlier, a new condition that can lead to a "non-detection" of a source is that it was observed on a non-photometric night.

Once again, we can plot the light curve that we have generated.

In [5]:
targ_det = mag != -99.          # identify epochs where targ is detected

# plot mag vs HJD along with uncertainties
figure(figsize(13,8))
errorbar(HJD[targ_det], mag[targ_det], magunc[targ_det], fmt='ok')
ylim(np.max(mag[targ_det]) + 0.3, np.min(mag[targ_det]) - 0.3)
xlim(np.min(HJD[targ_det]) - 15, np.max(HJD[targ_det]) + 15)
xlabel(r'${\rm HJD}$', fontsize=18)
ylabel(r'$R \; {\rm mag}$', fontsize=18)

print "The LC scatter = %5.3f mag" % np.std(mag[targ_det])
The LC scatter = 0.041 mag

Success!

And, unlike our previous efforts, the light curve has a scatter that is similar to what we would expect for a star of this brightness.

That being said, the above result leaves much to be desired. For starters, deciding that nights with a zeropoint RMSE < 0.04 are photometric requires the use of a somewhat arbitrary hard cut for the data (ex. are all nights with RMSE = 0.39 clearly photometric, while all nights with RMSE = 0.41 clearly not photometric?).

  • Please allow me a quick aside here -- while hard cuts are sometimes a necessary evil when analyzing astromical data, it is best to avoid them whenever possible. In the cases when they cannot be avoided, I would recommend searching for a physical criterion to make the cut before adopting a (semi-)arbitrary cut based on the distribution in a given dataset.

Even more importantly, however, the above cuts remove well over half of the observations from our light curve. Surely those osbervations are not useless? Otherwise, there would be no point in opening a telescope on a night when the conditions were not going to be photometric. Fortunately, we can use a technique called differential photometry to measure brightness changes even on nights with non-photometric conditions.

Making a Light Curve: III. Differential Photometry

The basic principle of differential photometry is the following: the vast majority of stars are non-variable (as measured from ground-based instruments, this is not true for Kepler data), and thus, relative brightness variations for a star of interest can be determined via direct comparison to nearby non-variable stars. For most astronomical detectors, the assumption that non-stellar brightness variations will be uniform for all sources on the detector is typically safe (though caution must be used with detectors as large as those used by PTF).

To proceed with the construction of a differential light curve for our source of interest, we must identify nearby non-variable stars in order to produce the differential correction to be applied to each photometric measurement. We will do this by identifying all bright (\(R < 18\) mag) stars within 5 arcmin of our target star.

In [6]:
diffphot_r = 5.0                    # diff. star match radius in arcmin

ZPscat = np.empty(Nims)             # ZP scatter for each image

for imnum, cat in enumerate(files):
    ZPscat[imnum] = fits.getval(cat, "PCALRMSE", 2)

# vvvvv arbitrary assumption time! vvvvv
# image with smallest ZP scatter is adopted as ref
refim = np.array(files)[ZPscat == min(ZPscat)][0]

refdat = fits.getdata(refim)
refdat_RA = refdat["ALPHAWIN_J2000"]            # RA list for cat
refdat_Dec = refdat["DELTAWIN_J2000"]           # Dec list for cat

source_sep = np.hypot((refdat_RA - targ_ra)*np.cos(np.radians(targ_dec)), 
    refdat_Dec - targ_dec)
# for diff. phot. refs only select sources that: (i) are not the targ; 
# (ii) are < 5' from the targ; (iii) have no SExtractor flags; 
# (iv) are classified as stars by SExtractor; and (v) are R < 18 mag
ref_stars = ((source_sep > match_r/3600.) & (source_sep < diffphot_r/60.) & 
    (refdat["FLAGS"] == 0) & (refdat["CLASS_STAR"] > 0.9) & 
    (refdat["MAG_AUTO"] + refdat["ZEROPOINT"] < 18))

ref_star_ra = refdat["ALPHAWIN_J2000"][ref_stars] # ref star RA
ref_star_dec = refdat["DELTAWIN_J2000"][ref_stars] # ref star Dec

We identify photometric reference stars using the "best" image that contains our target. In this case, we have defined "best" as the image taken on the night with the smallest zeropoint scatter. This definition is somewhat arbitrary, but also sufficient for our purposes. We further require that our reference stars are less than 5 arcmin from our target, have no SExtractor flags, are brighter than 18 mag, and have a SExtractor star/galaxy classification > 0.9. Warning: this particular cut does not ensure that our list of reference stars is galaxy free. The inclusion of galaxies in a list of reference stars for relative photometry will seriously compromise the efficacy of a light curve. For simplicity we ignore that possibility here, as the vast majority of sources with CLASS_STAR > 0.9 will be stars.

Now that we have identified photometric reference stars, we must construct "raw" light curves from the source catalogs. We will also construct the "raw" light curve for the target, which we will simply store for now until we have determined the differential corrections.

In [7]:
Nref_stars = sum(ref_stars)    # number of phot. ref. stars

# make an array to hold mag for all ref stars
ref_star_rawmag = np.zeros((Nref_stars,Nims)) 
targ_rawmag = np.zeros(Nims)   # targ "raw" mag
targ_Sunc = np.zeros(Nims)     # targ SExtractor mag uncertainties
HJD_ims = np.zeros(Nims)       # HJD of all images

for imnum, cat in enumerate(files):
    Sdat = fits.getdata(cat)  # grab data for image
    Sdat_RA = Sdat["ALPHAWIN_J2000"]            # RA list for cat
    Sdat_Dec = Sdat["DELTAWIN_J2000"]           # Dec list for cat
    HJD_ims[imnum] = fits.getval(cat, "HJD", 2)
    for refnum, ra, dec in zip(range(Nref_stars), ref_star_ra, ref_star_dec):
        source_sep = np.hypot((Sdat_RA - ra)*np.cos(np.radians(dec)), 
            Sdat_Dec - dec)
        match = ((source_sep < match_r/3600.) & (source_sep == np.min(source_sep)))
        if sum(match) == 1:
            imzp = Sdat["ZEROPOINT"][match][0]
            immag = Sdat["MAG_AUTO"][match][0] + imzp
            ref_star_rawmag[ refnum, imnum ] = immag
        else:
            ref_star_rawmag[ refnum, imnum ] = -99.

    targ_sep = np.hypot((Sdat_RA - targ_ra)*np.cos(np.radians(targ_dec)), 
        Sdat_Dec - targ_dec)
    targ_match = ((targ_sep < match_r/3600.) & (targ_sep == np.min(targ_sep)))
    if sum(targ_match) == 1:
        tzp = Sdat["ZEROPOINT"][targ_match][0]
        tmag = Sdat["MAG_AUTO"][targ_match][0] + tzp
        targ_rawmag[imnum] = tmag
        targ_Sunc[imnum] = Sdat["MAGERR_AUTO"][targ_match][0]
    else:
        targ_rawmag[imnum] = -99.
        targ_Sunc[imnum] = -99.
       

Using essentially the same for loop that we adopted perviously, we construct "raw" light curves for each of the reference stars, as well as our target of interest. Note that once again we record the mag and mag uncertainty to be -99 in lieu of "non-detections."

To determine relative photometric corrections for each image, we must measure how much each star deviates from its normal (i.e. mean) brightness. Thus, we must determine the mean mag for each of the reference stars.

In [39]:
# determine the mean magnitude for each ref star
ref_star_meanmag = np.empty(Nref_stars)
for refnum, ref_rawlc in enumerate(ref_star_rawmag):
    ref_star_meanmag[refnum] = np.mean(ref_rawlc[ref_rawlc != -99.])

Using the mean mag for each star, we will now determine how much each star deviates from its mean brightness in each image.

In [40]:
# determine differential offsets
ref_star_dm = 1.*ref_star_rawmag
for refnum, mean_ref_mag in enumerate(ref_star_meanmag):
    ref_star_det = ref_star_dm[refnum] != -99.   # use detections only 
    ref_star_dm[refnum][ref_star_det] -= mean_ref_mag

To determine the differential corrections that must be applied to the light curve of our target, we average the deviations from the mean for the reference stars on an image by image basis.

In [41]:
# determine the mean offset for each image
im_delmag = np.empty(Nims)
for imnum in range(Nims):
    ref_star_im_det = ref_star_dm[:,imnum] != -99.
    im_delmag[imnum] = np.mean(ref_star_dm[:,imnum][ref_star_im_det])

We can now apply these offsets to the "raw" light curve of the target in order to get a differential light curve.

In [42]:
targ_det = targ_rawmag != -99.    # targ detection epochs
targ_mag = targ_rawmag[targ_det] - im_delmag[targ_det]

We can plot results of our differential light curve, and for comparison we will also plot the lightcurve measured in part I of this exercise, where we assumed perfect photometric calibration for all of the PTF images.

In [43]:
figure(figsize=(14,8))
errorbar(HJD_ims[targ_det], targ_mag, targ_Sunc[targ_det], 
         fmt='o', label=r"${\rm Diff. \; LC}$")
errorbar(HJD_ims[targ_det], targ_rawmag[targ_det]-1, targ_Sunc[targ_det], 
         fmt='o', label=r"${\rm Raw \; LC} - 1$")
ylim(np.max(targ_mag) + 0.3, np.min(targ_mag) - 0.3 - 1)
xlim(np.min(HJD_ims[targ_det]) - 15, np.max(HJD_ims[targ_det]) + 15)
xlabel(r'${\rm HJD}$', fontsize=18)
ylabel(r'$R \; {\rm mag}$', fontsize=18)
legend(numpoints=1, fontsize=15)

print "The LC scatter = %5.3f mag" % np.std(targ_mag)
The LC scatter = 0.055 mag

Success!!!

We now have a light curve that includes every epoch on which our target was detected, and it has a photometric scatter similar to what is expected for a star of this brightness.

There are some downsides to using this differential method. For instance, all of the small corrections that we applied to the raw light curve enable us to measure relative changes in the brightness of our target, but they remove any information about the calibrated brightness of this source. This is important if, for instance, one wants to measure the absolute magnitude of the source in question. We also adopted some shortcuts in the interest of time and simplicity, shortcuts that, in most cases, would render the light curves generated above insufficient for inclusion in a scientific journal.

Ultimately, the decisions made when constructing a light curve (e.g., PSF vs. Aperture, relative vs. absolute, etc.) must be driven by the science. If \(10\%\) accuracy will do, don't waste time trying to achieve a \(\sim1\%\) calibration.

For a differential photometry method that fully accounts for missing detections and the uncertainties associated with the reference stars we refer the curious reader to Honeycutt 1992.

Can you think of any poor assumptions of possible corrections to the methods outlined above?

References

In []: