This notebook provides a demonstration of simulating NIRISS wide-field slitless spectroscopic observations of the Hubble WFC3/IR Ultra-Deep Field. We simulate two dispersion position angles and three of the broad-band blocking filters, for a total of six exposures.
For the example here, we use the Hubble H-band F140W image as the direct image, though in practice the NIRISS WFSS APT templates will enable/require users to obtain direct images with NIRISS itself in the same blocking filter as used for the spectroscopic exposures.
The deep Hubble images allows us to simulate objects with reasonable spatial morphologies. Of course the NIRISS / NIRCam images themselves will have better image quality and different PSFs than that of the hubble images used here.
To approximate the spectra of the simulated objects, we take the best fit galaxy template spectrum for each object as fit to the many-band photometric catalog of the entire GOODS-S field produced by the 3D-HST survey (Skelton et al. 2015).
Steps of the simulation are as follows:
Download the UDF F140W images
Download the NIRISS configuration files
Use astropy utilities (photutils, wcs, coordinates) to generate an object catalog from the Hubble UDF images
Generate fake images with dimensions and pixel scale appropriate for NIRISS but with calibrated WCS.
Simulate full 2D exposures with morphologies from the HST image and spectra from the photo-z fits
Once the simulation is complete, we then analyze the simulated data as if they were real observations. NB: The simulated images here are constructed to look like HST WFC3/IR FLT images and the actual data format of the calibrated exposures produced by the JWST pipeline will be different.
The initial analysis steps of the full-frame ismulated exposures are as follows:
Read in the simulated exposures and the reference image / catalog / segmentation image
Generate a first-pass contamination model assuming flat f_lambda SEDs normalized to the detection filter, in this case WFC3/IR F140W.
Generate a second pass refined contamination model fitting 2nd-order polynomials to the observed spectra
With the full contamination model computed as above, now take cutout spectra of a given object from each of the input simulations in the separate position angles / blocking filters. Use the grizli wrapper tools to fit the spectrum to determine the redshift and emission line fluxes and also extract drizzled continuum-subtracted emission line maps.
In [1]:
%matplotlib inline
In [2]:
import os
import glob
from collections import OrderedDict
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec
mpl.rcParams['figure.figsize'] = (10.0, 6.0)
mpl.rcParams['font.size'] = 14
mpl.rcParams['savefig.dpi'] = 72
import numpy as np
from astropy.io import fits
import astropy.wcs
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
import drizzlepac
import pysynphot as S
import grizli
import grizli.fake_image
print('\n grizli version: %s' %(grizli.__version__))
In [3]:
## Fetch current placeholder NIRISS configuration files
if not os.path.exists(os.path.join(os.getenv('GRIZLI'), 'CONF/NIRISS.F150W.conf')):
cwd = os.getcwd()
os.chdir(os.path.join(os.getenv('GRIZLI'), 'CONF'))
url = 'http://www.stsci.edu/jwst/instruments/niriss/software-tools'
os.system('wget {url}/wfss-simulations/niriss-wfss.tar'.format(url=url))
os.system('tar xvf niriss-wfss.tar cookbook/CONF')
os.system('ln -s cookbook/CONF/NIRISS* .')
In [4]:
# Work in /tmp for now
os.chdir('/tmp/')
In [5]:
## Fetch WFC3/IR F140W images of the UDF
if not os.path.exists('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits'):
url = 'https://archive.stsci.edu/missions/hlsp/xdf/'
wget='wget --no-check-certificate'
os.system('%s %s/hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits' %(wget, url))
os.system('%s %s/hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_wht.fits' %(wget, url))
## 3D-HST photometry around the HUDF field, and an extension with best-fit photo-z SED
if not os.path.exists('udf_3dhst_cat.fits'):
os.system('wget http://www.stsci.edu/~brammer/Grizli/Demos/udf_3dhst_cat.fits')
In [6]:
## Make an object catalog / segmentation image with photutils
sci = fits.open('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits')
wht = fits.open('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_wht.fits')
rms = 1/np.sqrt(wht[0].data)
rms[~np.isfinite(rms)] = 1000
dq = rms > 900
wcs = astropy.wcs.WCS(sci[0].header)
# Run the detection with the grizli / photutils wrapper
from grizli.utils import detect_with_photutils
cat, seg = detect_with_photutils(sci[0].data, err=rms, dq=dq, seg=None,
detect_thresh=1.4, npixels=6, grow_seg=3,
gauss_fwhm=2.0, gsize=3, wcs=wcs, save_detection=False,
root='udf_f140w_photutils', background=None, gain=None,
AB_zeropoint=26.452,
rename_columns={'ycentroid': 'y_flt', 'xcentroid': 'x_flt',
'dec_icrs_centroid': 'dec',
'ra_icrs_centroid': 'ra'},
clobber=True, verbose=True)
# code expects this columns later....
cat['NUMBER'] = cat['id']
In [7]:
# Get matches from 3D-HST catalog
ref_3dhst = fits.open('udf_3dhst_cat.fits')
ref_cat = Table.read(ref_3dhst[1])
gs = SkyCoord(ra=ref_cat['ra']*u.degree, dec=ref_cat['dec']*u.degree)
cat_rd = SkyCoord(ra=cat['ra'], dec=cat['dec'])
gs_idx, d2d, d3d = cat_rd.match_to_catalog_sky(gs)
has_gs_match = np.where(d2d < 2*u.arcsec)[0]
# Use 3D-HST mags because quick photutils catalog has some issues,
# perhaps related to background subtraction
gs_mag = 25-2.5*np.log10(ref_cat['F204'])
cat['MAG_AUTO'] = gs_mag[gs_idx]
cat.write('udf_f140w_photutils.cat', format='ascii.commented_header')
fits.writeto('udf_f140w_photutils_seg.fits', data=np.cast[int](seg),
header=sci[0].header, clobber=True)
In [8]:
## Setup fake images, cendered in the UDF/XDF
ra, dec = 53.1592277508136, -27.782056346146
pa_aper = 128.589
np.random.seed(1)
# Crude exposure parameters. Rough read noise and backgrounds coded in
# grizli.fake_image to make some reasonable noise estimate
EXPTIME = 1.e4 # 10 ks ~ 4 HST orbits
NEXP = 10 # divided between 10 exposures
# JWST NIRISS, three filters & two orients
for filt in ['F115W', 'F150W', 'F200W']:
for theta in [0,90]:
h, wcs = grizli.fake_image.niriss_header(filter=filt, ra=ra, dec=dec,
pa_aper=pa_aper+theta)
print('Filter: {filter}, Background: {bg} e/s/pix, RN: {RN} e/exp'.format(filter=filt,
bg=h['BACKGR'], RN=h['READN']))
output = 'niriss_{filt}_{theta:02d}_flt.fits'.format(filt=filt, theta=theta)
grizli.fake_image.make_fake_image(h, output=output, exptime=EXPTIME, nexp=NEXP)
In [9]:
# Load GroupFLT for simulation, NB: input files are just noise
sim = grizli.multifit.GroupFLT(grism_files=glob.glob('niriss_*flt.fits'), direct_files=[],
ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits',
seg_file='udf_f140w_photutils_seg.fits',
catalog='udf_f140w_photutils.cat',
cpu_count=0, pad=200)
In [10]:
# First pass, flat continuum
sim.compute_full_model(mag_limit=27)
In [11]:
# Compute model grism spectra for 3D-HST matches based on full photo-z templates
detection_bp = S.ObsBandpass('wfc3,ir,f140w')
for ix in has_gs_match:
templ = ref_3dhst['WAVE'].data*(1+ref_cat['zbest'][gs_idx[ix]])
tempf = ref_3dhst['FLAMBDA'].data[gs_idx[ix],:]
# Needs to be normalized to unity in the detection band
spec = S.ArraySpectrum(wave=templ, flux=tempf, waveunits='angstroms', fluxunits='flam')
spec = spec.renorm(1., 'flam', detection_bp)
id = cat['id'][ix]
#print(id)
sim.compute_single_model(id, mag=cat['MAG_AUTO'][ix], size=-1, store=False,
spectrum_1d=[spec.wave, spec.flux], get_beams=None,
in_place=True)
In [12]:
# Blotted reference image
# "grism" exposures are still empty, just noise
fig = plt.figure(figsize=[9,9*2./3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
ax.imshow(sim.FLTs[i].direct['REF'], vmin=-0.01, vmax=0.05, cmap='viridis',
origin='lower')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.grid(color='w', alpha=0.8)
ax.text(100,100,sim.FLTs[i].grism_file, color='w', size=10, ha='left', va='bottom')
fig.tight_layout(pad=0.1)
In [13]:
# Blotted segmentation image
# "grism" exposures are still empty, just noise
fig = plt.figure(figsize=[9,9*2./3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
ax.imshow(sim.FLTs[i].seg, vmin=-0.01, vmax=3000, cmap='viridis', origin='lower')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.grid(color='w', alpha=0.8)
ax.text(100,100,sim.FLTs[i].grism_file, color='w', size=10, ha='left', va='bottom')
fig.tight_layout(pad=0.1)
In [14]:
# "grism" exposures are still empty, just noise
fig = plt.figure(figsize=[9,9*2./3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
ax.imshow(sim.FLTs[i].grism['SCI'], vmin=-0.01, vmax=0.05, cmap='viridis', origin='lower')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.grid(color='w', alpha=0.8)
ax.text(100,100,sim.FLTs[i].grism_file, color='w', size=10, ha='left', va='bottom')
fig.tight_layout(pad=0.1)
In [15]:
# Model stored in `FLTs[i].model` attribute
fig = plt.figure(figsize=[9,9*2./3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
# show as if it were the rotated grism
if (i % 2) > 0:
ax.imshow(np.rot90(sim.FLTs[i].model,-1), vmin=-0.01, vmax=0.05, cmap='viridis',
origin='lower')
else:
ax.imshow(sim.FLTs[i].model, vmin=-0.01, vmax=0.05, cmap='viridis', origin='lower')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.grid(color='w', alpha=0.8)
ax.text(100,100,sim.FLTs[i].grism_file, color='w', size=10, ha='left', va='bottom')
fig.tight_layout(pad=0.1)
In [16]:
# Update SCI extension of the fake FLT images with the models just computed
for flt in sim.FLTs:
print('Update', flt.grism_file)
orig_flt = fits.open(flt.grism_file, mode='update')
orig_flt['SCI'].data += flt.model[flt.pad:-flt.pad, flt.pad:-flt.pad]
orig_flt.flush()
In [17]:
# Now reload for fitting
grp = grizli.multifit.GroupFLT(grism_files=glob.glob('niriss_*flt.fits'), direct_files=[],
ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits',
seg_file='udf_f140w_photutils_seg.fits',
catalog='udf_f140w_photutils.cat',
cpu_count=0, pad=200)
In [18]:
# First pass contamination model, flat continuum
grp.compute_full_model(mag_limit=26)
In [19]:
# Refine the (polynomial) continuum model for brighter objects
grp.refine_list(poly_order=2, mag_limits=[16, 24], verbose=False)
In [20]:
## Fit parameters
pzfit, pspec2, pline = grizli.multifit.get_redshift_fit_defaults()
# Redshift fit
pzfit ['zr'] = [0.5, 2.4]
pzfit['dz'] = [0.005, 0.0005]
# Drizzled line maps
pline = {'kernel': 'square', 'pixfrac': 0.8, 'pixscale': 0.06, 'size': 10}
# Full rectified 2D spectrum
pspec2 = {'NY': 20, 'dlam': 50, 'spatial_scale': 1}
In [21]:
from importlib import reload
reload(grizli.multifit)
## emission line object
id = 898 # ID in the detection catalog / segm image
# Extract spectrum cutouts from individual FLTs
beams = grp.get_beams(id, size=40)
# Put them in a `MultiBeam` object
mb = grizli.multifit.MultiBeam(beams, fcontam=1, group_name='niriss-udf')
# Run the redshift fit and generate the emission line map
out = mb.run_full_diagnostics(pzfit=pzfit, pspec2=pspec2, pline=pline,
GroupFLT=grp, prior=None, verbose=False)
fit, fig, fig2, hdu2, hdu_line = out
In [22]:
cmap = 'viridis_r'
In [23]:
# "Beams" are extracted spectra of a given order. Have attributes for contam, model etc.
fig = plt.figure(figsize=[9,9*1.2/3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
beam = mb.beams[i]
ax.imshow(beam.grism['SCI'], vmin=-0.01, vmax=0.05, cmap=cmap, origin='lower',
aspect='auto')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.text(0.1,0.1,beam.grism.parent_file, color='k', backgroundcolor='w',
transform=ax.transAxes, size=10, ha='left', va='bottom')
fig.axes[0].set_ylabel('Observed\nspectrum')
fig.tight_layout(pad=0.1)
In [24]:
# Each beam carries with it a static contamination model extracted from the full field
fig = plt.figure(figsize=[9,9*1.2/3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
beam = mb.beams[i]
ax.imshow(beam.contam, vmin=-0.01, vmax=0.05, cmap=cmap, origin='lower', aspect='auto')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.text(0.1,0.1,beam.grism.parent_file, color='k', backgroundcolor='w',
transform=ax.transAxes, size=10, ha='left', va='bottom')
fig.axes[0].set_ylabel('Contamination model')
fig.tight_layout(pad=0.1)
In [25]:
# Under the hood, the fitting is done by specifying a single 1D template, which
# is used to generate model 2D spectra for each beam
fig = plt.figure(figsize=[9,9*1.2/3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
beam = mb.beams[i]
ax.imshow(beam.model, vmin=-0.01, vmax=0.05, cmap=cmap, origin='lower', aspect='auto')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.text(0.1,0.1,beam.grism.parent_file, color='k', backgroundcolor='w',
transform=ax.transAxes, size=10, ha='left', va='bottom')
fig.axes[0].set_ylabel('Observed\nspectrum')
fig.tight_layout(pad=0.1)
In [26]:
# Goodness of fit is computed by comparing the models in the full 2D pixel space
fig = plt.figure(figsize=[9,9*1.2/3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
beam = mb.beams[i]
ax.imshow(beam.grism['SCI'] - beam.contam - beam.model, vmin=-0.01, vmax=0.05, cmap=cmap,
origin='lower', aspect='auto')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.text(0.1,0.1,beam.grism.parent_file, color='k', backgroundcolor='w',
transform=ax.transAxes, size=10, ha='left', va='bottom')
fig.axes[0].set_ylabel('Full residuals')
fig.tight_layout(pad=0.1)
In [27]:
# Trivial demo model, dropout in the middle of the F115W filter
xspec = np.arange(0.8,2.4,0.02)*1.e4
yspec = (xspec/1.4e4)**(-2) # Beta=-2
yspec[xspec < 1.2e4] = 0.
plt.plot(xspec, yspec)
mb.compute_model(spectrum_1d=[xspec, yspec])
In [28]:
fig = plt.figure(figsize=[9,9*1.2/3])
for ix, i in enumerate([0,2,4,1,3,5]):
ax = fig.add_subplot(2,3,ix+1)
beam = mb.beams[i]
ax.imshow(beam.beam.model, vmin=-0.01, vmax=0.05, cmap=cmap, origin='lower', aspect='auto')
ax.set_xticklabels([]); ax.set_yticklabels([])
ax.text(0.1,0.1,beam.grism.parent_file, color='k', backgroundcolor='w',
transform=ax.transAxes, size=10, ha='left', va='bottom')
fig.axes[0].set_ylabel('Observed\nspectrum')
fig.tight_layout(pad=0.1)
In [29]:
# Emission line map
line = fits.open('niriss-udf_zfit_00898.line.fits')
print(line[0].header['HASLINES'])
line.info()
In [30]:
cmap = 'cubehelix_r'
fig = plt.figure(figsize=[12,3])
ax = fig.add_subplot(141)
ax.imshow(line['DSCI'].data, vmin=-0.01, vmax=0.02, cmap=cmap, origin='lower')
ax.text(5,5,'F140W direct image', ha='left', va='bottom')
ax = fig.add_subplot(142)
ax.imshow(line['LINE', 'Ha'].data, vmin=-0.01, vmax=0.02, cmap=cmap, origin='lower')
ax.text(5,5,r'H$\alpha$', ha='left', va='bottom')
ax = fig.add_subplot(143)
ax.imshow(line['LINEWHT', 'Ha'].data, vmin=-0.01, vmax=20000, cmap='gray', origin='lower')
ax.text(5,5,r'H$\alpha$, weight', ha='left', va='bottom', color='w')
ax = fig.add_subplot(144)
ax.imshow(line['LINE', 'OIII'].data, vmin=-0.03, vmax=0.06, cmap=cmap, origin='lower')
ax.text(5,5,r'[OIII]$\lambda$4959,5007', ha='left', va='bottom')
for ax in fig.axes:
ax.set_xticklabels([]); ax.set_yticklabels([])
fig.tight_layout(pad=0.1)