Extended Aperture Photometry#

Use case: Measure galaxy photometry in a field. Related to JDox Science Use Case #22.
Data: WST simulated NIRCam images from JADES JWST GTO extragalactic blank field.
(Williams et al. 2018) https://ui.adsabs.harvard.edu/abs/2018ApJS..236…33W.
Tools: photutils, matplotlib, scipy, scikit.
Cross-intrument: potentially MIRI.
Documentation: This notebook is part of a STScI’s larger post-pipeline Data Analysis Tools Ecosystem.

Introduction#

This notebook uses photutils to detect objects/galaxies in NIRCam deep imaging. Detections are first made in a F200W image, then isophotal photometry is obtained in all 9 filters (F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, F444W). The catalog is loaded back in and some simple analysis is performed on the full catalog and on an individual galaxy.

The notebook analyzes only the central 1000 x 1000 pixels (30” x 30”) of the full JADES simulation. These cutouts have been staged at STScI with permission from the authors (Williams et al.).

NOTE: The photometry is aperture matched, but no PSF corrections are made. For more accurate color measurements, PSF corrections should be implemented, given the large range of wavelengths (and thus PSF FWHM) spanning a factor of >4.

NOTE: The simulated JADES images have different units (e-/s) than JWST pipeline products (MJy/sr).

NOTE: An exposure map is missing but required to calculate flux uncertainties.

To Do#

  • PSF corrections

  • Check accuracy of photometry against simulated JADES catalog

  • Exposure map required for input to error calculation

  • ABmag units cannot be written to ecsv file (astropy update coming soon)

  • plot with text labels looks horrible (I wish cursor hover would show id number instead)

  • Fix plot secondary axis: mag vs. flux

  • requirements.txt file – but I don’t know what versions are “required”

  • rest of Robel’s comments: https://github.com/spacetelescope/dat_pyinthesky/pull/82#pullrequestreview-355206337

Download WEBBPSF Data#

import os
import tarfile
import urllib.request

boxlink = 'https://stsci.box.com/shared/static/34o0keicz2iujyilg4uz617va46ks6u9.gz'                                                           
boxfile = './webbpsf-data/webbpsf-data-1.0.0.tar.gz'

webbpsf_folder = './webbpsf-data'

# Check whether the specified path exists or not
isExist = os.path.exists(webbpsf_folder)
if not isExist:
  # Create a new directory because it does not exist 
  os.makedirs(webbpsf_folder)
    
urllib.request.urlretrieve(boxlink, boxfile)
gzf = tarfile.open(boxfile)
gzf.extractall(webbpsf_folder)

Import packages#

import os

import numpy as np

from astropy.convolution import convolve
from astropy.io import fits
from astropy.stats import gaussian_fwhm_to_sigma
from astropy.table import QTable
import astropy.units as u
from astropy.visualization import make_lupton_rgb, SqrtStretch, ImageNormalize, simple_norm
import astropy.wcs as wcs

from photutils.background import Background2D, MedianBackground
from photutils.segmentation import (detect_sources, deblend_sources, SourceCatalog,
                                    make_2dgaussian_kernel, SegmentationImage)
from photutils.utils import calc_total_error

Matplotlib setup for plotting#

There are two versions

  • notebook – gives interactive plots, but makes the overall notebook a bit harder to scroll

  • inline – gives non-interactive plots for better overall scrolling

import matplotlib.pyplot as plt
import matplotlib as mpl

# Use this version for non-interactive plots (easier scrolling of the notebook)
%matplotlib inline

# Use this version if you want interactive plots
# %matplotlib notebook

# These gymnastics are needed to make the sizes of the figures
# be the same in both the inline and notebook versions
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}

mpl.rcParams['savefig.dpi'] = 80
mpl.rcParams['figure.dpi'] = 80

Create list of images to be loaded and analyzed#

baseurl = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/nircam_photometry/'

filters = 'F090W F115W F150W F200W F277W F335M F356W F410M F444W'.split()

# Data images [e-/s]
imagefiles = {}
for filt in filters:
    filename = f'jades_jwst_nircam_goods_s_crop_{filt}.fits'
    imagefiles[filt] = os.path.join(baseurl, filename)

# Weight images (Inverse Variance Maps; IVM)
weightfiles = {}
for filt in filters:
    filename = f'jades_jwst_nircam_goods_s_crop_{filt}_wht.fits'
    weightfiles[filt] = os.path.join(baseurl, filename)

Load detection image: F200W#

filt = 'F200W'
infile = imagefiles[filt]
hdu = fits.open(infile)
data = hdu[0].data
imwcs = wcs.WCS(hdu[0].header, hdu)

weight = fits.open(weightfiles[filt])[0].data

Report image size and field of view#

ny, nx = data.shape
pixscale = wcs.utils.proj_plane_pixel_scales(imwcs)[0] 
pixscale *= imwcs.wcs.cunit[0].to('arcsec')
outline = '%d x %d pixels' % (ny, nx)
outline += ' = %g" x %g"' % (ny * pixscale, nx * pixscale)
outline += ' (%.2f" / pixel)' % pixscale
print(outline)

Create color image (optional)#

# 3 NIRCam short wavelength channel images
r = fits.open(imagefiles['F200W'])[0].data
g = fits.open(imagefiles['F150W'])[0].data
b = fits.open(imagefiles['F090W'])[0].data

rgb = make_lupton_rgb(r, g, b, Q=5, stretch=0.02)  # , minimum=-0.001

fig = plt.figure(figsize=(8, 8))
ax = plt.subplot(projection=imwcs)
plt.imshow(rgb, origin='lower')
plt.xlabel('Right Ascension')
plt.ylabel('Declination')
fig.tight_layout()
plt.subplots_adjust(left=0.15)

Detect Sources and Deblend using photutils#

https://photutils.readthedocs.io/en/latest/segmentation.html

# For detection, requiring 5 connected pixels 2-sigma above background

# Measure background and set detection threshold
bkg_estimator = MedianBackground()
bkg = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator)
threshold = bkg.background + (2. * bkg.background_rms)

# Before detection, smooth image with Gaussian FWHM = 3 pixels
smooth_kernel = make_2dgaussian_kernel(3.0, size=3)
convolved_data = convolve(data, smooth_kernel)

# Detect and deblend
segm_detect = detect_sources(convolved_data, threshold, npixels=5)

segm_deblend = deblend_sources(convolved_data, segm_detect, npixels=5, nlevels=32, contrast=0.001)

# Save segmentation map of detected objects
segm_hdu = fits.PrimaryHDU(segm_deblend.data.astype(np.uint32), header=imwcs.to_header())
segm_hdu.writeto('JADES_detections_segm.fits', overwrite=True)

Measure photometry (and more) in detection image#

https://photutils.readthedocs.io/en/latest/segmentation.html#centroids-photometry-and-morphological-properties

#error = bkg.background_rms
# Input weight should be exposure map. Fudging for now.
error = calc_total_error(data, bkg.background_rms, weight/500)
#cat = source_properties(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)
cat = SourceCatalog(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)

Show detections alongside images (optional)#

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(9.5, 6))
# For RA,Dec axes instead of pixels, add: , subplot_kw={'projection': imwcs})

# Color image
ax[0, 0].imshow(rgb, origin='lower')
ax[0, 0].set_title('Color Image')

# Data
norm = simple_norm(data, 'sqrt', percent=99.)
ax[0, 1].imshow(data, origin='lower', cmap='Greys_r', norm=norm)
ax[0, 1].set_title('Detection Image F200W')

# Segmentation map
cmap = segm_deblend.make_cmap(seed=12345)
ax[0, 2].imshow(segm_deblend, origin='lower', cmap=cmap, interpolation='nearest')
ax[0, 2].set_title('Detections (Segmentation Image)')

# Weight
ax[1, 0].imshow(weight, origin='lower', cmap='Greys_r', vmin=0)
ax[1, 0].set_title('Weight Image F200W')

# RMS
ax[1, 1].imshow(bkg.background_rms, origin='lower', norm=None)
ax[1, 1].set_title('Background RMS')

# Total error including Poisson noise
norm = simple_norm(error, 'sqrt', percent=99.)
ax[1, 2].imshow(error, origin='lower', norm=norm)
ax[1, 2].set_title('RMS + Poisson noise')

fig.tight_layout()

View all measured quantities in detection image (optional)#

cat.to_table()

Only keep some quantities#

columns = 'label xcentroid ycentroid sky_centroid area semimajor_sigma semiminor_sigma ellipticity orientation gini'.split()
tbl = cat.to_table(columns=columns)
tbl.rename_column('semimajor_sigma', 'a')
tbl.rename_column('semiminor_sigma', 'b')

tbl

Convert measured fluxes (data units) to magnitudes#

https://docs.astropy.org/en/stable/units/

https://docs.astropy.org/en/stable/units/equivalencies.html#photometric-zero-point-equivalency

https://docs.astropy.org/en/stable/units/logarithmic_units.html#logarithmic-units

# not detected: mag =  99; magerr = 1-sigma upper limit assuming zero flux
# not observed: mag = -99; magerr = 0
def fluxes2mags(flux, fluxerr):
    nondet = flux < 0  # Non-detection if flux is negative
    unobs = (fluxerr <= 0) + (fluxerr == np.inf)  # Unobserved if flux uncertainty is negative or infinity

    mag = flux.to(u.ABmag)
    magupperlimit = fluxerr.to(u.ABmag)  # 1-sigma upper limit if flux=0

    mag = np.where(nondet, 99 * u.ABmag, mag)
    mag = np.where(unobs, -99 * u.ABmag, mag)

    magerr = 2.5 * np.log10(1 + fluxerr/flux) 
    magerr = magerr.value * u.ABmag

    magerr = np.where(nondet, magupperlimit, magerr)
    magerr = np.where(unobs, 0*u.ABmag, magerr)
    
    return mag, magerr

# Includes features I couldn't find in astropy:
# mag = 99 / -99 for non-detections / unobserved
# flux uncertainties -> mag uncertainties

Multiband photometry using isophotal apertures defined in detection image#

(Similar to running SourceExtractor in double-image mode)

filters = 'F090W F115W F150W F200W F277W F335M F356W F410M F444W'.split()
for filt in filters:
    infile = imagefiles[filt]
    print(filt)
    print(infile)
    print(weightfiles[filt])
    hdu = fits.open(infile)
    data = hdu[0].data
    zp = hdu[0].header['ABMAG'] * u.ABmag  # zeropoint
    weight = fits.open(weightfiles[filt])[0].data
    
    # Measure background
    bkg = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator)
    #error = bkg.background_rms
    error = calc_total_error(data, bkg.background_rms, weight/500)
                             
    # Measure properties in each image of previously detected objects 
    filtcat = SourceCatalog(data-bkg.background, segm_deblend, wcs=imwcs, background=bkg.background, error=error)

    # Convert measured fluxes to fluxes in nJy and to AB magnitudes
    filttbl = filtcat.to_table()
    tbl[filt+'_flux']    = flux    = filttbl['segment_flux']     * zp.to(u.nJy)
    tbl[filt+'_fluxerr'] = fluxerr = filttbl['segment_fluxerr'] * zp.to(u.nJy)

    mag, magerr = fluxes2mags(flux, fluxerr)
    #mag = mag * u.ABmag  # incompatible with file writing
    tbl[filt+'_mag']    = mag.value
    tbl[filt+'_magerr'] = magerr.value

View complete results (optional)#

tbl

Save photometry as output catalog#

tbl.write('JADESphotometry.ecsv', overwrite=True)

!head -175 JADESphotometry.ecsv  # show the first 175 lines

Reformat output catalog for readability (optional)#

# Remove units (pixels) from area
tbl['area'] = tbl['area'].value.astype(int)

# Replace sky_centroid with ra, dec
tbl['ra'] = tbl['sky_centroid'].ra.degree
tbl['dec'] = tbl['sky_centroid'].dec.degree

columns = list(tbl.columns)
columns = columns[:3] + ['ra', 'dec'] + columns[4:-2]

tbl = tbl[columns]

for column in columns:
    tbl[column].info.format = '.4f'

tbl['ra'].info.format = '11.7f'
tbl['dec'].info.format = '11.7f'

tbl['label'].info.format = 'd'
tbl['area'].info.format = 'd'

tbl.write('JADESphotometry.cat', format='ascii.fixed_width_two_line', delimiter=' ', overwrite=True)

!head -10 JADESphotometry.cat  # show the first 10 lines

Start new session and analyze results#

Load catalog and segmentation map#

# Catalog: ecsv format preserves units for loading in Python notebooks
tbl = QTable.read('JADESphotometry.ecsv')

# Reconstitute filter list
filters = []
for param in tbl.columns:
    if param[-4:] == '_mag':
        filters.append(param[:-4])

# Segmentation map
segmfile = 'JADES_detections_segm.fits'
segm = fits.open(segmfile)[0].data
segm = SegmentationImage(segm)

Plot number counts vs. magnitude#

fig = plt.figure(figsize=(8, 4))

filt = 'F200W'
mag1 = tbl[filt + '_mag']

mag1 = mag1[(0 < mag1) & (mag1 < 90)]  # detections only
n = plt.hist(mag1, histtype='step', label=filt)

plt.xlabel('AB magnitude')
plt.ylabel('Number counts')
plt.legend()

Plot F200W vs. F090W magnitudes and look for dropouts#

#import mplcursors
# Would love a better solution here!

mag1 = tbl['F090W_mag']
mag2 = tbl['F200W_mag']

# Only plot detections in F200W
det2 = (0 < mag2) & (mag2 < 90)

mag1 = mag1[det2]
mag2 = mag2[det2]
ids = tbl['label'][det2]

plt.figure(figsize=(8, 4))

plt.plot(mag1, mag2, '.')

for i in range(len(mag1)):
    plt.text(mag1[i], mag2[i], ids[i])

plt.xlabel('F090W AB magnitude')
plt.ylabel('F200W AB magnitude')

Look at one object#

# Could select object by position
#x, y = 905, 276
#id = segm.data[y,x]

# Select by ID number
id = 261  # F090W dropout
obj = tbl[id-1]

obj

obj['ellipticity']

segmobj = segm.segments[segm.get_index(id)]
segmobj

Show the object in all the images#

fig, ax = plt.subplots(2, len(filters)+1, figsize=(9.5, 3.5), sharex=True, sharey=True)

ax[0, 0].imshow(rgb[segmobj.slices], origin='lower', extent=segmobj.bbox.extent)
ax[0, 0].set_title('Color')

cmap = segm.make_cmap(seed=12345)  # ERROR
ax[1, 0].imshow(segm.data[segmobj.slices], origin='lower', extent=segmobj.bbox.extent, cmap=cmap,
               interpolation='nearest')
ax[1, 0].set_title('Segment')

for i in range(1, len(filters)+1):
    filt = filters[i-1]

    # Show data on top row
    data = fits.open(imagefiles[filt])[0].data
    stamp = data[segmobj.slices]
    norm = ImageNormalize(stretch=SqrtStretch())  # scale each filter individually
    ax[0, i].imshow(stamp, extent=segmobj.bbox.extent, cmap='Greys_r', norm=norm, origin='lower')
    ax[0, i].set_title(filt.upper())

    # Show weights on bottom row
    weight = fits.open(weightfiles[filt])[0].data
    stamp = weight[segmobj.slices]
    # set black to zero weight (no exposure time / bad pixel)
    ax[1, i].imshow(stamp, extent=segmobj.bbox.extent, vmin=0, cmap='Greys_r', origin='lower')

ax[0, 0].set_ylabel('Data')
ax[1, 0].set_ylabel('Weight')

Plot SED (Spectral Energy Distribution)#

fig, ax = plt.subplots(figsize=(8, 6))

for filt in filters:
    lam = int(filt[1:4]) / 100
    plt.errorbar(lam, obj[filt+'_flux'].value, obj[filt+'_fluxerr'].value, marker='.', c='b')

plt.axhline(0, c='k', ls=':')
plt.xlim(0, 5)
plt.xlabel('Wavelength ($\mu$m)')
plt.ylabel('Flux (nJy)')

mlim = 31.4
flim = mlim * u.ABmag
flim = flim.to(u.nJy).value

# Add AB magnitudes as secondary x-axis at right
# https://matplotlib.org/gallery/subplots_axes_and_figures/secondary_axis.html#sphx-glr-gallery-subplots-axes-and-figures-secondary-axis-py
def AB2nJy(mAB):
    m = mAB * u.ABmag
    f = m.to(u.nJy)
    f = f.value
    f = np.where(f > flim, f, flim)
    return f

def nJy2AB(F_nJy):
    f = F_nJy * u.nJy
    m = f.to(u.ABmag)
    m = m.value
    m = np.where(m < mlim, m, mlim)
    return m
    
plt.ylim(flim, plt.ylim()[1])

secax = ax.secondary_yaxis('right', functions=(nJy2AB, AB2nJy))
secax.set_ylabel('magnitude (AB)')

Magnitude conversion fails for flux <= 0#

fig, ax = plt.subplots(figsize=(8, 6))

for filt in filters:
    lam = int(filt[1:4]) / 100
    plt.errorbar(lam, obj[filt+'_flux'].value, obj[filt+'_fluxerr'].value, marker='.', c='b')

plt.axhline(0, c='k', ls=':')
plt.xlim(0, 5)
plt.xlabel('Wavelength ($\mu$m)')
plt.ylabel('Flux (nJy)')

f0 = 10**(0.4 * 31.4)  # flux [nJy] at zero magnitude
b0 = 1.e-12  # this should be filter dependent

# Add AB magnitudes as secondary x-axis at right
# https://matplotlib.org/gallery/subplots_axes_and_figures/secondary_axis.html#sphx-glr-gallery-subplots-axes-and-figures-secondary-axis-py
def AB2nJy(m):
    f = np.sinh(-0.4 * m * np.log(10) - np.log(b0)) * 2 * b0 * f0
    return f

# Luptitudes
# https://www.sdss.org/dr12/algorithms/magnitudes/
def nJy2AB(f):
    m = -2.5 / np.log(10) * (np.arcsinh((f / f0) / (2 * b0)) + np.log(b0))
    return m

#plt.ylim(flim, plt.ylim()[1])

secax = ax.secondary_yaxis('right', functions=(nJy2AB, AB2nJy))
secax.set_ylabel('asinh magnitude (AB)')