Photometry

On this page we’ll detect sources and perform simple aperture photometry (along other measurements) for each pre-reduced image. The same procedure could of course also be applied to coadded images.

As before, you could copy or write the code shown below in a script, or alternatively directly download this page as a jupyter notebook file.

To run the code, you’ll need the module dataredconfig.py, as explained here.

Note

The photometry presented here is sufficient for our purposes, but remains extremely simple. Notably, it ignores the effect of variable seeing (in time and/or accross filters). The next steps to improve upon the minimal approach shown here would be to perform a PSF homogeneization or to use PSF-fitting.

import dataredconfig
from pathlib import Path

import numpy as np
import astropy
import astropy.visualization
from astropy import units as u

%matplotlib widget
import matplotlib
from matplotlib import pyplot as plt

import ccdproc
import photutils.background
import photutils.aperture
import photutils.segmentation
import photutils.psf

from time import time

# We'll ignore some astropy warnings that get raised as our FITS headers (from NINA) are not 100% standards compliant.
import warnings
warnings.simplefilter('ignore', category=astropy.wcs.FITSFixedWarning)

# Setting the correct path to the pre-reduced science frames to work on:
#science_files_dir = dataredconfig.work_dir / "LIGHT_PRERED" # In case astrometric calibration was already performed before the pre-reduction
science_files_dir = dataredconfig.work_dir / "ASTROMETRY"
science_files = ccdproc.ImageFileCollection(science_files_dir, keywords=dataredconfig.ifc_header_keywords)

# Where to write the catalogs:
photometry_dir = dataredconfig.work_dir / "PHOTOMETRY"
photometry_dir.mkdir(exist_ok=True)

# Overview of all available files:
science_files.summary
# You may want to try this instead, for a better display:
#science_files.summary.show_in_notebook()

# Optional: selecting the science frames to run on : 
object_to_process = "HD92670"

science_files = science_files.filter(object=object_to_process)

#science_files.summary

We will measure photometry on all these files. But we start by picking one reference image, used to detect sources and to check the overall procedure.

ref_image_filepath = science_files.files[0]
ref_image = ccdproc.CCDData.read(ref_image_filepath)

if(ref_image.wcs is None): # Performing forced photometry requires each frame to have a good WCS. This test is just to avoid obvious mistakes.
    raise RuntimeError("The image has no WCS, make sure you specify the correct directory (i.e., with WCS) above!")

Estimating and subtracting the background

Some of the following tasks will be run on several images, and we therefore define them as very simple functions. Of course, depending on the data to be processed, you might want to adjust (hard-coded) parameters in these functions. They are only provided as examples.

def estimate_background(image):
    """Function to estimate the background (in form of an image),
    and the root-mean-square deviation (RMS) of the background.
    
    Parameters:
    image: a CCDData object

    Returns:
    background, background_rms: two numpy arrays
    """
    
    bkg = photutils.background.Background2D(
        image.data,
        box_size=(500, 500),
        filter_size=(3, 3),
        bkg_estimator=photutils.background.MedianBackground())
    
    return bkg.background, bkg.background_rms

# We run this on the reference image
background, background_rms = estimate_background(ref_image)

# Visualize the background image
plt.figure()
plt.imshow(background, origin='lower', cmap='Greys_r', interpolation='nearest')
plt.tight_layout()
plt.show()

# We subtract the background from the reference image
ref_image_noback = ref_image.subtract(ccdproc.CCDData(background, unit="adu"))

# And we define a "threshold" image for this reference frame,
# whose values are multiples of the estimated RMS of the background in each pixel.
# (RMS is an estimate of the standard deviation, so 3*RMS is "3 sigma" above the background)
threshold = 3.0 * background_rms # this is an array

Source detection and measurement

We now detect sources (as groups of connected pixels above the threshold) and measure their exact positions (among other parameters), on the selected reference image. Let’s nevertheless group all this in a function again.

def detect_and_measure(image, threshold):
    """Function to detect and measure sources.
    This is close to what Source Extractor typically does.
    
    Parameters:
    image: a CCDData object
    threshold: a 2D-array with threshold values

    Returns:
    An Astropy Table with source positions (from  windowed centroids, 
    both in pixel and sky coordinates) and other measurements.
    """
    data = image.data 

    # Filtering (convolving with a 2D Gaussian) for better source detection
    kernel = photutils.segmentation.make_2dgaussian_kernel(fwhm=4.0, size=21)  # FWHM in pixels
    convolved_data = astropy.convolution.convolve(data, kernel)

    # Segmenting the image
    finder = photutils.segmentation.SourceFinder(npixels=4, connectivity=4, progress_bar=False)
    segment_map = finder(convolved_data, threshold)
    #segment_map = finder(data, threshold) # One could also run this on the unfiltered image

    # And measuring sources
    source_catalog = photutils.segmentation.SourceCatalog(data, segment_map, convolved_data=convolved_data, wcs=image.wcs)
    source_table = source_catalog.to_table(
        columns=["xcentroid_win", "ycentroid_win", "sky_centroid_win",
                 "fwhm", "max_value", "kron_flux", "segment_flux", "segment_area"]
        )

    return astropy.table.Table(source_table) # We convert the QTable to a Table, better for later FITS writing.

# Run this on the reference image
ref_catalog = detect_and_measure(ref_image_noback, threshold)

print(f"Reference catalog: {len(ref_catalog)} sources detected.")

After running a detection, it’s almost mandatory to check the result with a visualization! We do this with the following interactive Figure.

# Figure that overplots the image and the positions of the detected sources

plt.figure(figsize=(10, 6))
ax = plt.subplot(projection=ref_image_noback.wcs)
ax.imshow(ref_image_noback.data, origin='lower', cmap='Greys_r', interpolation='nearest',
    norm=astropy.visualization.simple_norm(ref_image_noback.data, stretch="log", vmin=-20, vmax=2000))
ax.scatter(
    ref_catalog["xcentroid_win"],
    ref_catalog["ycentroid_win"],
    transform=ax.get_transform('pixel'),
    s=50, # The size of these markers is not related to any measurement apertures!
    edgecolor='red', facecolor='none'
    )
ax.grid(color='white', ls='solid')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
plt.tight_layout()
plt.show()

Zoom in and check if all relevant sources are detected. If there are issues, one might have to adjust the sensitivity of the detection.

# Make this catalog smaller, for testing:
ref_catalog = ref_catalog[ref_catalog["kron_flux"] > 10000] # pick the brightest sources
len(ref_catalog)

# Rerun the vizualisation above to see the effect...

Note

A large reference catalog will slow down the photometry. It’s helpful to adapt the flux threshold above to keep mainly those stars you actually care about, and not every tiny detection.

# For reference, we write the reference catalog in this directory.
# It contains the aperture positions.
ref_catalog.write(photometry_dir / f"ref_catalog_{object_to_process}.fits", overwrite=True)

Forced aperture photometry

We’ll now perform forced aperture photometry on all science images of our target, with a series of aperture sizes (radii given in arcseconds in the code below). This means that we use the sky coordinates of the sources detected in the reference frame to place the same apertures on all dithered exposures of the target.

# We group all steps of the photometric measurement into one function:

def measure_photometry(image, sky_positions, run_simple_fit=False):
    """Function for forced aperture photometry given the provided sky_positions.
    This function first estimates and subtracts the background from the image,
    and performs aperture photometry both on the background-subtracted image and on the background alone.

    Parameters:
    image: a CCDData object
    sky_positions: a list or column of aperture positions, in sky coordinates (i.e. as SkyCoord objects)
    run_fit: if True, a simple 2D Gaussian fit is also run, providing FWHM and flux measurements. Take significantly more time.

    Returns:
    An Astropy Table with measurements for several aperture radii at the given positions.
    """

    # Estimating and subtracting the background:
    tstart = time()    
    background, _ = estimate_background(image)
    image_noback = image.subtract(ccdproc.CCDData(background, unit="adu"))
    tend = time()
    print(f"Background subtraction took {tend-tstart:.2f} sec")
    
    
    # Aperture photometry
    aperture_radii = [4, 6, 8, 10] # In arcsec
    tstart = time()
    # We generate a separate catalog for each aperture radius, and combine these afterwards
    aperture_catalogs = []
    for r in aperture_radii:

        apertures = photutils.aperture.SkyCircularAperture(sky_positions, r=r*u.arcsec)

        # Photometry on the background-subtracted image:
        aperture_measurements = photutils.aperture.ApertureStats(image_noback.data, aperture=apertures, wcs=image_noback.wcs)
        # Same photometry on the background image:
        background_aperture_measurements = photutils.aperture.ApertureStats(background.data, aperture=apertures, wcs=image_noback.wcs)
        
        aperture_catalog = aperture_measurements.to_table(columns=["sum"])
        background_aperture_catalog = background_aperture_measurements.to_table(columns=["sum"])

        # Note: ApertureStats also provides "fwhm" and "elongation", among other stats.
        # While these measurements are fast and can be useful, they are computed from unweighted moments,
        # and will therefore be strongly influenced by the aperture size and the noise.
        # Just for reference, this is how to get them:
        #aperture_catalog = aperture_measurements.to_table(columns=["fwhm", "elongation", "sum"])        

        # Rename columns by adding the aperture radius, to get for example "sum_4" instead of sum
        for name in aperture_catalog.colnames:
            aperture_catalog.rename_column(name, new_name=f"{name}_{str(r)}")
        for name in background_aperture_catalog.colnames:
            background_aperture_catalog.rename_column(name, new_name=f"back_{name}_{str(r)}")
         
        aperture_catalogs.append(aperture_catalog)
        aperture_catalogs.append(background_aperture_catalog)
    
    # And now merge these catalogs into a single one:
    merged_catalog = astropy.table.hstack(aperture_catalogs,
        join_type="exact", metadata_conflicts="silent"
        )
    tend = time()
    print(f"Photometry took {tend-tstart:.2f} sec")


    # Optional simple 2D Gaussian fit for FWHM-measurement of sources:
    if run_simple_fit:
        tstart = time() 
        # Compute pixel-coordinate positions of the positions:
        image_pixel_positions = astropy.wcs.utils.skycoord_to_pixel(sky_positions, ref_image.wcs)
        xypos = np.array(image_pixel_positions).transpose() # Rearranging the data
        # Perform the fitting:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            fit_catalog = photutils.psf.fit_2dgaussian(image_noback.data, xypos=xypos, fit_shape=15, fwhm=6, fix_fwhm=False).results
        # Copy over some relevant columns into our final catalog:
        merged_catalog["fwhm_fit"] = fit_catalog["fwhm_fit"] # The FWHM obtained form the fit, in pixels
        merged_catalog["flux_fit"] = fit_catalog["flux_fit"] # The flux from the fit, in ADU
        merged_catalog["q_fit"] = fit_catalog["qfit"] # A quality-of-fit metric, see doc of photutils.psf for details.
        tend = time()
        print(f"2D Gaussian fitting took {tend-tstart:.2f} sec")

    return merged_catalog

# As a test, we run this on the ref_image.
# For this we read the ref image again, as the one in memory already has the background subtracted.
ref_image = ccdproc.CCDData.read(ref_image_filepath)
ref_photo_cat = measure_photometry(ref_image, ref_catalog["sky_centroid_win"], run_simple_fit=True)
print(ref_photo_cat.colnames)

Chimney Plot

As an illustration of how this photometric catalog can be used, we create a basic and very common scatter plot: the so-called “chimney plot”, named after its shape. The plot shows instrumental magnitude against FWHM (or some other size measurement) of all detected sources.

Note

Potentially, your image contains both unresolved sources (stars) and resolved sources (e.g., galaxies). This plot should show you a clear “stellar locus” (i.e., where the stars are), and directly tell you the seeing of this particular image. Can you measure it? A potential fine “smoke plume” coming from your chimney is related to the saturation of the brightest stars.

# Just as an example, let's use the flux within the apertures of 6 arcseconds: 
ref_photo_cat["mag_6"] = -2.5 * np.log10(ref_photo_cat["sum_6"].value)

# We create a scatter plot of FWHM versus mag, and color the datapoints according to the q_fit parameter (see legend).
# High values of q_fit indicate that the source was not well fitted by the simple 2D-Gaussian.
plt.figure()
plt.scatter(
    ref_photo_cat["fwhm_fit"].value, 
    ref_photo_cat["mag_6"], 
    c=ref_photo_cat["q_fit"].value,
    vmin=0.0, vmax=0.5,
    cmap="jet",
    s=10
)
plt.gca().invert_yaxis() # Bright is at the top!
plt.colorbar(label="Fit quality metric: sum(abs(fit residuals)) / (fitted flux)")
plt.xlabel("FWHM [pix]")
plt.xlim(0, 15)
plt.ylabel("Instrumental magnitude")
plt.tight_layout()
plt.show()

Performing the measurements on all images

Finally, we apply this measurement function to all images of the object.

If you don’t need the Gaussian fitting, you can run this faster by setting run_simple_fit=False below, in the call of measure_photometry().

# We read again the ref catalog from file:
ref_catalog = astropy.table.Table.read(photometry_dir / f"ref_catalog_{object_to_process}.fits")


n_files = len(science_files.summary)
# And now we loop over all selected exposures:
for (i, science_file) in enumerate(science_files.files):

    #if i < 104: # could use something like that to relaunch on selected files
    #    continue
    
    image_filename = Path(science_file).stem
    print(f"=== Processing image {i+1}/{n_files}: {image_filename} ===")

    image = ccdproc.CCDData.read(science_file)

    # We build a dict containing some meta-information of the image taken from the FITS header.
    # We'll later write this as meta-information of the catalogue.
    meta_from_image = {key: image.header[key] for key in [
        "DATE-OBS", "EXPTIME", "IMAGETYP", "AIRMASS", "PIERSIDE", "FILTER", "OBJECT",
        "FOCUSPOS", "FOCTEMP", "RA", "DEC", "SET-TEMP", "CCD-TEMP", "GAIN", "OFFSET"
        ]}
    
    # Now run the measurements:
    catalog = measure_photometry(image, ref_catalog["sky_centroid_win"], run_simple_fit=True)

    # Add the information from the FITS header
    catalog.meta.update(meta_from_image)

    # Could copy a few columns from the reference catalog
    # But no need to do this here: we keep the reference catalog, and will use it later. 
    #catalog["sky_centroid_win"] = ref_catalog["sky_centroid_win"]

    # And write all this as a FITS table
    catalog.write(photometry_dir / f"{image_filename}.fits", overwrite=True)