Aperture Photometry

Aperture Photometry#

Kernel Information and Read-Only Status#

To run this notebook, please select the “Roman Calibration” kernel at the top right of your window.

This notebook is read-only. You can run cells and make edits, but you must save changes to a different location. We recommend saving the notebook within your home directory, or to a new folder within your home (e.g. file > save notebook as > my-nbs/nb.ipynb). Note that a directory must exist before you attempt to add a notebook to it.

Imports#

We are using standard libraries from the Numpy and Astropy stack.

  • astropy

  • numpy

  • copy

  • matplotlib

  • photutils is an Astropy-affiliated package for photometry

  • roman_datamodels opens and validates WFI data files

  • asdf opens WFI data files

  • os for checking if files exist

  • s3fs streams data from Simple Storage Service (S3) buckets on Amazon Web Services (AWS)

  • stpsf for retrieving the encircled energy profile of stars

from astropy.table import Table
from astropy.stats import SigmaClip, sigma_clipped_stats
import asdf
import crds
import copy
import matplotlib.pyplot as plt
import numpy as np
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry, ApertureStats
import roman_datamodels as rdm
import os
import s3fs
import stpsf

Introduction#

This notebook explains how to perform forced aperture photometry (also known as forced photometry or aperture photometry) on Roman WFI images. Aperture photometry is useful for measuring the integrated fluxes for a set of specified source positions and aperture sizes. This technique is often valuable for extracting fluxes of known sources when it is impractical to fit for their positions and light distributions, or when simplicity and speed are high priorities. For example, aperture photometry can be used for:

  • Faint sources. If a source is too faint, then it can be difficult to fit for the source centroid, light profile, and flux.

  • Sources detected at other wavelengths. If a source is detected in a given bandpass, and the source size is expected to be the same across wavelengths, then it can be useful to define one aperture and measure the source fluxes across multiple filters.

  • Time-series evolution of sources. If a source brightness is decreasing over time, then we can use earlier observations obtained when the source was bright to specify an aperture and extract the flux in a time series.

Here, we cover a simple example using the photutils package to extract fluxes from a catalog of known sources.

Please note that the WFI calibration is currently measured from ground test data, and changes are expected after launch.

Tutorial Data#

Image Data#

In this tutorial, we use a Level 2 (L2; calibrated rate image) WFI data file that is the result of RomanCal processing of a Level 1 (L1; uncalibrated ramp cube) simulated file created with Roman I-Sim. If you have already worked through the tutorials “Simulating WFI Imaging Data with Roman I-Sim” and “Calibrating WFI Exposures with RomanCal,” then you may already have these files saved locally. If not, then these files are also stored in the Nexus S3 bucket. For more information on how to access these data, see the Data Discovery and Access tutorial.

As a reminder, the file we are using is a L2 file meaning that the data were processed to flag and/or correct for detector-level effects (e.g., saturation, classic non-linearity, etc.), and that the ramp was fitted into a count rate image in units of Data Numbers (DN) per second.

#Stream the files from the S3 bucket if they are not in local storage

if os.path.exists('r0003201001001001004_0001_wfi01_f106_cal.asdf'):
    f = rdm.open('r0003201001001001004_0001_wfi01_f106_cal.asdf')
else:
    fs = s3fs.S3FileSystem(anon=True)
    asdf_dir_uri = 's3://stpubdata/roman/nexus/soc_simulations/tutorial_data/'
    asdf_file_uri = asdf_dir_uri + 'r0003201001001001004_0001_wfi01_f106_cal.asdf'
    f = rdm.open(fs.open(asdf_file_uri, 'rb'))

image = f.data
err = f.err
dq = np.bool(f.dq)
meta = f.meta
wcs = f.meta.wcs

Source Catalog#

We also have access to a source catalog that was used to simulate the WFI image. It contains stars and galaxies, which are labeled as “PSF” and “SER” under the column type. Source fluxes are available in all WFI filters (F062, F087, F106, F129, F146, F158, F184, and F213) and are sampled from a lognormal distribution. Note that fluxes are all given in maggies, which are defined as ${\rm maggie} \equiv 10^{-0.4 m_{AB}}$, for an AB apparent magnitude $m_{AB}$.

For galaxies, morphological parameters like n (Sersic index), half_light_radius, pa (position angle), and ba (axis ratio) are also provided in the catalog. These are sampled according to fiducial (and likely unrealistic) distributions.

#Stream the files from the S3 bucket if they are not in local storage

if os.path.exists('full_catalog.ecsv'):
    cat = Table.read('full_catalog.ecsv')
else:
    fs = s3fs.S3FileSystem(anon=True)
    asdf_dir_uri = 's3://stpubdata/roman/nexus/soc_simulations/tutorial_data/'
    asdf_file_uri = asdf_dir_uri + 'full_catalog.ecsv'
    cat = Table.read(fs.open(asdf_file_uri, 'rb'), format='ascii.ecsv')

We can also display the first five rows and all columns of the catalog:

cat[:5]
Table length=5
radectypenhalf_light_radiuspabaF062F087F106F129F146F158F184F213
float64float64str3float64float64float64float64float64float64float64float64float64float64float64float64
271.6295510520961-0.03965028067096237SER3.5162152365994530.8056893348693848252.73039245605470.72357076406478883.6905281053378758e-091.1479291295302119e-081.2472877841673835e-086.86679868522333e-091.4913801305027619e-091.0117606308313043e-084.943988773931096e-092.402744225804554e-09
270.7539155586658-1.175434326974954SER1.03993985924195490.946480393409729336.083068847656250.48005193471908573.1115434673267828e-095.331648900153141e-099.049362148516593e-096.6658691899590394e-092.0099595321454444e-091.5802696040623232e-081.0280717610555712e-081.096286439405958e-08
270.86489243319850.5861049874846664SER2.82067304151365851.327717542648315473.264846801757810.69381815195083621.3153379718744418e-086.075517422488019e-093.1850986292880634e-097.20732051817663e-096.512435035688213e-096.66888944067523e-094.776757656088648e-096.596417634341378e-09
270.70544440611565-0.1905627389524337SER1.20807579181093570.573934018611908273.726013183593750.57724624872207643.2391787030405794e-091.9316861443741118e-093.5669256437387276e-091.699146068290247e-089.461913919039944e-097.503919152718197e-104.7414845383286774e-095.4636619672976394e-09
271.07632728190820.28776831750222365SER2.73018077847475875.05323600769043236.82214355468750.9073979854583747.582467986821939e-087.082852704343168e-083.547680762494565e-072.0032227610045084e-072.793574651605013e-072.3345998556578706e-073.1447200399270514e-071.8613071972595208e-07

We can convert (RA, Dec) to (x, y) positions on the WFI01 detector. Note that several of these positions may be NaN (Not a Number) as they fall outside the bounding box of the WCS.

x_cat, y_cat = wcs.invert(cat["ra"].data, cat["dec"].data)
cat['x'] = x_cat
cat['y'] = y_cat

It might be helpful to quantify the number of sources there are before we try to visualize them! We first create a mask of sources that actually fall on the detector, and then divide the remaining sources into stars and galaxies.

mask = (np.isfinite(x_cat) & np.isfinite(y_cat))
print(f"Number of sources on detector: {sum(mask)}")

stars = cat[mask & (cat["type"] == "PSF")]
galaxies = cat[mask & (cat["type"] == "SER")]
print(f"Number of stars: {len(stars)}\nNumber of galaxies: {len(galaxies)}")

Number of sources on detector: 1724
Number of stars: 1579
Number of galaxies: 145

We can now plot the distribution of source fluxes using a histogram. We will plot the flux from the input catalog in units of AB magnitudes. Recall that the input catalog is in units of maggies. The transformation from maggies to AB magnitudes is given by:

$m_{\mathrm{AB}} = -2.5 \times \log_{10}(f_{\mathrm{maggies}})$

stars['abmag'] = -2.5 * np.log10(stars['F106'])
galaxies['abmag'] = -2.5 * np.log10(galaxies['F106'])

fig, ax = plt.subplots()
ax.hist(stars['abmag'], bins=50, range=[16, 24], log=True, label='Stars')
ax.hist(galaxies['abmag'], bins=50, range=[16, 24], log=True, label='Galaxies')

ax.set_xlabel(r'$m_{\mathrm{AB, F106}}$')
ax.set_ylabel('Number of Sources')
ax.xaxis.set_inverted(True)
ax.legend()

ax.grid(alpha=0.3)
../../_images/4127df376ea370ff88f19a89f1ca7abb8a8ce19b9cf3cbe7e45c4f6ba4345440.png

Forced Aperture Photometry#

Forced aperture photometry uses predefined source positions to place apertures and measure the flux within them. We use the source catalog information provided above for this. If source positions are unavailable, you will need to perform source detection first to determine their locations (see Additional Resources).

Correcting for Pixel Area#

First, we need to correct our data for geometric distortion, which causes pixels to have different areas projected on the sky. The WFI distortion is described by a 2-D polynomial, and the area of each pixel is given by the determinant of the Jacobian matrix of this polynomial. A reference file, called the Pixel Area Map (PAM), provides the area per pixel of each WFI detector normalized to the nominal pixel area. The nominal pixel area is the pixel area at the tangent point about which the distortion polynomial is defined.

First, we can use the crds Python package to request the best pixel area map reference file that matches our observation. In this case, CRDS needs to know the instrument name (WFI), detector name (WFI01), and the start time of the exposure. We can pull all of this information out of the metadata of the L2 file. Note that we need to specify the start time as a string, so we use meta.exposure.start_time.isot to convert the astropy.time.Time object in the metadata to an ISOT string (i.e., it has the format YYYY-MM-DDTHH:MM:SS.SSS). We set the argument ignore_cache=True below to force CRDS to update the files in our local cache. If this cell is re-run, we can set ignore_cache=False as we are confident we have pulled the latest information from the CRDS server. We also specify that the observatory is “roman” and that we want the “area” reference file type:

crds_result = crds.getreferences({'ROMAN.META.INSTRUMENT.DETECTOR': meta.instrument.detector,
                                  'ROMAN.META.INSTRUMENT.NAME': meta.instrument.name,
                                  'ROMAN.META.EXPOSURE.START_TIME': meta.exposure.start_time.isot}, 
                                 observatory='roman', reftypes=['area'], ignore_cache=True)

print(crds_result)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_skycells_0002.rmap      712 bytes  (1 / 38 files) (0 / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_saturation_0003.rmap    2.4 K bytes  (2 / 38 files) (712 / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_refpix_0003.rmap    2.3 K bytes  (3 / 38 files) (3.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_readnoise_0004.rmap    3.1 K bytes  (4 / 38 files) (5.5 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_photom_0004.rmap    2.3 K bytes  (5 / 38 files) (8.6 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-tweakregstep_0001.rmap      568 bytes  (6 / 38 files) (10.9 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-sourcedetectionstep_0001.rmap      643 bytes  (7 / 38 files) (11.4 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-skymatchstep_0001.rmap      568 bytes  (8 / 38 files) (12.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-saturationstep_0001.rmap      580 bytes  (9 / 38 files) (12.6 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-resamplestep_0001.rmap      568 bytes  (10 / 38 files) (13.2 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-refpixstep_0001.rmap      556 bytes  (11 / 38 files) (13.8 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-rampfitstep_0001.rmap      574 bytes  (12 / 38 files) (14.3 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-photomstep_0001.rmap      556 bytes  (13 / 38 files) (14.9 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-outlierdetectionstep_0001.rmap      649 bytes  (14 / 38 files) (15.5 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-linearitystep_0001.rmap      574 bytes  (15 / 38 files) (16.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-jumpstep_0001.rmap      712 bytes  (16 / 38 files) (16.7 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-highlevelpipeline_0001.rmap      807 bytes  (17 / 38 files) (17.4 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-fluxstep_0001.rmap      544 bytes  (18 / 38 files) (18.2 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-flatfieldstep_0001.rmap      574 bytes  (19 / 38 files) (18.8 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-exposurepipeline_0001.rmap      800 bytes  (20 / 38 files) (19.3 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-dqinitstep_0001.rmap      556 bytes  (21 / 38 files) (20.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-darkcurrentstep_0001.rmap      586 bytes  (22 / 38 files) (20.7 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_pars-assignwcsstep_0001.rmap      574 bytes  (23 / 38 files) (21.3 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_matable_0003.rmap      690 bytes  (24 / 38 files) (21.9 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_mask_0003.rmap    2.3 K bytes  (25 / 38 files) (22.5 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_linearity_0005.rmap    2.4 K bytes  (26 / 38 files) (24.8 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_ipc_0002.rmap    2.3 K bytes  (27 / 38 files) (27.2 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_inverselinearity_0005.rmap    2.6 K bytes  (28 / 38 files) (29.5 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_gain_0003.rmap    2.3 K bytes  (29 / 38 files) (32.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_flat_0006.rmap    8.5 K bytes  (30 / 38 files) (34.3 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_epsf_0003.rmap   16.3 K bytes  (31 / 38 files) (42.8 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_distortion_0002.rmap    3.6 K bytes  (32 / 38 files) (59.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_dark_0006.rmap    5.2 K bytes  (33 / 38 files) (62.7 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_area_0002.rmap    2.3 K bytes  (34 / 38 files) (67.9 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_apcorr_0002.rmap    2.3 K bytes  (35 / 38 files) (70.1 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_abvegaoffset_0002.rmap    2.5 K bytes  (36 / 38 files) (72.5 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_wfi_0039.imap    2.4 K bytes  (37 / 38 files) (74.9 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/mappings/roman/roman_0041.pmap       344 bytes  (38 / 38 files) (77.3 K / 77.7 K bytes)

CRDS - INFO -  Fetching  /home/runner/crds_cache/references/roman/wfi/roman_wfi_area_0009.asdf   66.8 M bytes  (1 / 1 files) (0 / 66.8 M bytes)

{'area': '/home/runner/crds_cache/references/roman/wfi/roman_wfi_area_0009.asdf'}

We see that a dictionary is returned. The keys of the dictionary are the reference file types (e.g., “area”), and the values are the path to the correct file in your local CRDS cache. If the file is not in your cache, then the crds.getreferences() function will download the file to your local cache. At this time, CRDS does not support streaming reference files into memory from the CRDS server, but work is ongoing to adapt CRDS to cloud applications.

Next, let’s take a look at the pixel area map that we now have:

pam = rdm.open(crds_result['area'])

fig, ax = plt.subplots()
img = ax.imshow(pam.data, origin='lower')
ax.set_xlabel('X science coordinate (pixels)')
ax.set_ylabel('Y science coordinate (pixels)')
ax.set_title(f'Pixel Area Map for {pam.meta.instrument.detector}')
plt.colorbar(img);
../../_images/8656eb0e3166f713ef0d609f15ea6c353e36006186c4486a40d655574c188448.png

Now let’s multiply the science data and the uncertainty array by the PAM to adjust for the relative sizes of the pixels on the sky:

image *= pam.data
err *= pam.data

Create Apertures#

For simplicity, we will focus on the stars from this point. From the catalogs, we now know the positions of every star. We can define set aperture radii in units of pixels; we choose radii of 3 pixels:

rad = 3  # aperture radius in pixels

star_positions = [(x, y) for y, x in zip(stars['y'].data, stars['x'].data)]
star_apertures = CircularAperture(positions=star_positions, r=rad)

Let’s create annulus apertures with $r_{in} = 20$ pixels and $r_{out} = 25$ pixels to measure local sky background values:

star_bkg_apertures = CircularAnnulus(positions=star_positions, r_in=20, r_out=25)

Source Visualization#

Let’s plot our apertures on a subsection of the image to visually confirm that we have placed our apertures correctly. There may be a few small, fuzzy galaxies in the image that are not marked since we are focusing on stars. Similarly, there may be several small features that are detector imperfections that are also not marked, but those will be masked in the data quality array.

fig, ax = plt.subplots(figsize=(9, 9))

# show the simulated image
ax.imshow(image, origin='lower', vmin=0, vmax=12, cmap='gray_r')

# plot circles over bright galaxies and stars
star_apertures.plot(color='cornflowerblue')
star_bkg_apertures.plot(color='cornflowerblue', ls=':')

# zoom in on 1/16th of the image
ax.set_xlim(1024, 1536)
ax.set_ylim(1024, 1536)

plt.axis('off')
plt.show()
../../_images/2bcdf3542b60037db76874824e99d09e237856ae01474544b1f50203aaf0b31b.png

Aperture Photometry with Photutils#

We can now perform aperture photometry on the selected sources. First we measure the photometry in a circular aperture on each source:

columns = ('id', 'xcentroid', 'ycentroid', 'sum', 'sum_err')
star_phot = ApertureStats(image, star_apertures, error=err)
phot_tab = star_phot.to_table(columns)
phot_tab[:5]
QTable length=5
idxcentroidycentroidsumsum_err
int64float64float64float64float64
14047.4972909152084640.29531539665862642.6815666570215.488672086430786
23983.703847283155878.6457742523635107.334819756334131.3541310292849316
34059.70511434428861198.8515033373212191.722647160864231.645956120742491
44066.514995718392858.5475960110696135.54745984489011.4682498573763318
53976.6420803446194716.3110121010073410.52864339682592.2933240160979875

And now let’s do the same for the background annuli using the ApertureStats class. We will use the median of the pixels within the annulus as our estimate of the sky background per pixel. We will need to adjust this to the area of the star circular apertures (after accounting for the masked pixels).

columns = ('id', 'xcentroid', 'ycentroid', 'median', 'std')
sigclip = SigmaClip(sigma=3.0, maxiters=5)
star_bkg = ApertureStats(image, star_bkg_apertures, sigma_clip=sigclip, mask=dq)
star_bkg.to_table(columns)[:5]
QTable length=5
idxcentroidycentroidmedianstd
int64float64float64float64float64
14047.778004866061640.74063694625070.63857176899909970.1630904260689358
23983.420897141212878.79133013385380.62351185083389280.16736493527375262
34059.499977124421198.92318992945570.65144473314285280.163246024548572
44063.877398350844857.92975068799630.62546834349632260.1680694038503727
53976.522889173548716.13572806560730.63738602399826050.17689512336474497

For our next steps, we will need the nominal pixel area in steradians (sr) and the zeropoint in megaJanskys per steradian (MJy/sr), both of which can be retrieved from the L2 metadata. Note that the pre-launch zeropoint is based on the unit response of the throughput curve of the integrated optical system (i.e., including all of the mirror reflectivities, optical element transmission, and detector quantum efficiency). The zeropoint also includes a factor of the median measured gain of the detector. The pre-launch zeropoint uncertainty is based solely on the standard deviation of the measured gain. After launch, these values will be updated using on-orbit observations of spectrophotometric standard stars and the photometric touchstone fields.

zpt = meta.photometry.conversion_megajanskys
pixel_area = meta.photometry.pixel_area

print(f'Zeropoint = {zpt} (MJy/sr) / (DN/s)')
print(f'Nominal pixel area = {pixel_area} sr')

Zeropoint = 0.7417487930165987 (MJy/sr) / (DN/s)
Nominal pixel area = 2.8083389953727505e-13 sr

Next, we subtract the sky background values from the source fluxes and propagate the uncertainty in the background into the total uncertainty:

flux_bkgsub = star_phot.sum - (star_bkg.median * star_apertures.area)
phot_tab['flux_jy'] = flux_bkgsub * zpt * 1e6 * pixel_area

We also need to correct our fluxes for the fraction of the starlight that falls outside of the aperture. The amount of light within the aperture as a function of the total light within a hypothetical infinite aperture is called the encircled energy and can be either measured from the image data using a series of apertures or can be determined from STPSF. In the future, a reference file in CRDS will also provide this information. For now, let’s use STPSF:

wfi = stpsf.WFI()
wfi.filter = meta.instrument.optical_element
wfi.detector = f'SCA{meta.instrument.detector[-2:]}'
wfi.detector_position = (2043.5, 2043.5)
wfi_psf = wfi.calc_psf()
ee_func = stpsf.measure_ee(wfi_psf)

We can also plot the encircled energy as a function of aperture radius:

stpsf.display_ee(wfi_psf)
../../_images/a0d2119840d3db7be2c467dc2e837f8f20a8d73b5915de53b1ee553e70555b92.png

With these settings, our encircled energy measurement is an average one for the whole detector, measured at the center of the detector, and for the default STPSF spectral energy distribution (SED), which is a 5700 K blackbody (similar to a G2V star). The variable ee_func above is a function that takes as its single argument the aperture radius in arcseconds and returns the percentage of the PSF contained within that aperture. The inverse of that encircled energy fraction is the correction factor that we must apply to our measurements. This is called the encircled energy correction or aperture correction. Let’s calculate that now for a circular aperture with a 3-pixel radius and assuming a pixel scale of 0.11 arcseconds per pixel (i.e., radius = 0.33 arcseconds):

ee_correction = 1/ee_func(rad * 0.11)
print(f'Encircled energy correction = {ee_correction}')

Encircled energy correction = 1.215690698785258

Now we can apply the aperture correction to our flux, and then convert the flux in Jy to maggies and AB magnitudes. For any source with a flux in maggies less than 0, we will assign a value of NaN for the AB magnitude. Recall that maggies are defined such that 1 maggie is equal to the reference flux for the AB magnitude system, which is 3,631 Jy. Thus, to convert our measurements from Jankys to maggies, we divide our fluxes by 3,631 Jy. We will also mask out any sources that have a flux of less than or equal to zero as those will not convert to magnitudes. Such sources exist where the background flux was statistically higher than the measured source flux.

phot_tab['flux_jy'] *= ee_correction
phot_tab['flux_maggies'] = phot_tab['flux_jy'] / 3631

# Mask out things with flux <= 0
mask = phot_tab['flux_maggies'] <= 0
phot_tab['flux_maggies'][mask] = np.nan
phot_tab['flux_jy'][mask] = np.nan

# Convert the flux to AB magnitudes
phot_tab['abmag'] = -2.5 * np.log10(phot_tab['flux_maggies'])

phot_tab[:5]
QTable length=5
idxcentroidycentroidsumsum_errflux_jyflux_maggiesabmag
int64float64float64float64float64float64float64float64
14047.4972909152084640.29531539665862642.6815666570215.4886720864307860.00066465604880100381.8305041277912526e-0716.843573218615585
23983.703847283155878.6457742523635107.334819756334131.35413102928493162.2716856849064223e-056.25636377005349e-0920.509195018982602
34059.70511434428861198.8515033373212191.722647160864231.6459561207424914.3887087392912794e-051.2086777029169043e-0819.79422372366485
44066.514995718392858.5475960110696135.54745984489011.46824985737633182.9847370460652654e-058.220151600289907e-0920.21280043225625
53976.6420803446194716.3110121010073410.52864339682592.29332401609798759.939781893483308e-052.7374778004635938e-0818.906623485105822

Let’s evaluate our results by plotting the difference between the simulated catalog and measured AB magnitudes. The blue points show the photometry of the stars, while the orange line shows the median offset between the input catalog and the measurements. The labeled uncertainty on the median is the standard deviation of the distribution. We also include a black line at an offset of 0 AB magnitudes as a reference.

fig, ax = plt.subplots(figsize=(7, 7))

# Compute the offset between the input catalogs (stars) and measured (phot_tab) AB magnitudes and plot them
offset = stars['abmag'] - phot_tab['abmag']
ax.scatter(stars['abmag'], offset, s=5, color='cornflowerblue', alpha=0.3)

# Compute the median offset ignoring NaNs and plot the median and a reference line at 0
_, med, std = sigma_clipped_stats(offset[~np.isnan(offset)], sigma=5, maxiters=2)
ax.axhline(med, c='darkorange', label=fr'Median $\Delta m = {med:.3f}$ AB mag, Std. $\Delta m = {std:.3f}$ AB mag', zorder=1e4)
ax.axhline(0, color='black', zorder=-100)

# Restrict the axis limits for visibility. You may need to change these if you alter the input
# data for this tutorial.
ax.set_xlim(21, 16.5)
ax.set_ylim(-0.5, 0.5)

# Label the axes
ax.set_xlabel(r'$m_{\mathrm{AB, input}}$')
ax.set_ylabel(r'$m_{\mathrm{AB, input}} - m_{\mathrm{AB, measured}}$')
ax.legend(loc='lower right', fontsize=10)

# Overplot a reference grid
ax.grid(alpha=0.1)
plt.savefig('photometry.png')
../../_images/249d7ba787fbab8b428b2187c5bafb9af7f7e24965f1a0adcfc5babf3a5df9df.png

Our measurements of stars in the image agree well with the input catalog fluxes at the level of approximately 0.9%. There could be several reasons why we see a small offset between the input and output magnitudes. First, we use a single value of the gain for our zeropoint, but we could do better and use a gain per pixel to convert the image from DN/s to photo-electrons per second and apply the zeropoint (with the gain removed). This would likely decrease the scatter in the plot. Secondly, we used a single aperture correction for the whole image based on the model of a star with an average SED and located at the center of the detector. A more detailed handling of the aperture correction may also decrease the scatter, and may also help bring the measurements into better alignment with the input catalog. Finally, a more careful treatment of other systematics (e.g., interpixel capacitance) may yield even better results.

A similar exercise in photometry can be repeated for the galaxies in the input catalog, but care should be taken when selecting apertures for galaxies as fixed aperture sizes will correspond to different physical sizes at the distances of galaxies. Other choices for galaxy apertures include, e.g., Kron apertures or isophotal apertures.

Additional Resources#

The Photutils documentation has additional tutorials for detecting and fitting sources in images:

About This Notebook#

Author: John F. Wu, Tyler Desjardins
Updated On: 2025-09-30

↑ Top of page
Contents