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 scrollinline
– 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)')