In [1]:
%matplotlib inline

Preliminaries

Install dependencies for running on cloud notebook servers like Google Collab.

In [2]:
try:
    import grizli
except:
    with open('requirements.txt','w') as fp:
        fp.write("""
cython
tqdm
astropy==4.2.1
drizzlepac==3.2.1
photutils==1.0.2
grizli>=1.3""")

    !pip install -r requirements.txt
In [3]:
import os
import glob
import grizli
import grizli.utils

# Make directories
for sub in ['CONF','iref','jref','templates']:
    dir = os.path.join(grizli.GRIZLI_PATH, sub)
    if not os.path.exists(dir):
        print(f'mkdir {dir}')
        os.mkdir(dir)

for sub in ['iref','jref']:
    if not os.getenv(sub):
        print(f'set {sub}')
        os.environ[sub] = os.path.join(grizli.GRIZLI_PATH, sub)

if len(glob.glob(os.path.join(grizli.GRIZLI_PATH, 'CONF/*'))) == 0:
    print('Fetch')
    grizli.utils.fetch_default_calibs(ACS=False) # to iref/iref
    grizli.utils.fetch_config_files()            # to $GRIZLI/CONF

if len(glob.glob(os.path.join(grizli.GRIZLI_PATH, 'templates/*'))) == 0:
    print('link templates')
    grizli.utils.symlink_templates(force=False)
In [4]:
# Pysynphot
pysyn_cdbs = os.getenv('PYSYN_CDBS')
if not pysyn_cdbs:
    pysyn_cdbs = os.path.join(grizli.GRIZLI_PATH, 'cdbs')
    if not os.path.exists(pysyn_cdbs):
        os.mkdir(pysyn_cdbs)
    
    os.environ['PYSYN_CDBS'] = pysyn_cdbs
    
if not os.path.exists(os.path.join(pysyn_cdbs, 'comp')):
    pwd = os.getcwd()
    os.chdir(pysyn_cdbs)

    os.system('wget http://ssb.stsci.edu/trds/tarfiles/synphot1.tar.gz')
    os.system('tar xzf synphot1.tar.gz')
    os.system('ln -s grp/redcat/trds/* .')
    os.system('rm synphot1.tar.gz')
    
    os.chdir(pwd)

try:
    import pysynphot as S
except:
    ! pip install pysynphot
    
import pysynphot as S
bp = S.ObsBandpass('wfc3,ir,f140w')
WARNING: VerifyWarning: Invalid keyword for column 1: Column disp option (TDISPn) failed verification: Format 26A is not recognized. The invalid value will be ignored for the purpose of formatting the data in this column. [astropy.io.fits.column]
WARNING: VerifyWarning: Invalid keyword for column 2: Column disp option (TDISPn) failed verification: Format 18A is not recognized. The invalid value will be ignored for the purpose of formatting the data in this column. [astropy.io.fits.column]
WARNING: VerifyWarning: Invalid keyword for column 3: Column disp option (TDISPn) failed verification: Format 56A is not recognized. The invalid value will be ignored for the purpose of formatting the data in this column. [astropy.io.fits.column]
WARNING: VerifyWarning: Invalid keyword for column 4: Column disp option (TDISPn) failed verification: Format 68A is not recognized. The invalid value will be ignored for the purpose of formatting the data in this column. [astropy.io.fits.column]

Multimission Simulation

This notebook demonstrates two analysis capabilities of the grizli software:

  1. The capability of providing a reference image and SExtractor-like segmentation file created entirely independently from any given grism exposure. The reference image is typically much deeper than a single direct image taken accompanying a grism exposure. The code assumes that the grism file is astrometrically aligned to the reference image, but the reference image can have any pixel scale.

  2. Simulation tools for comparison of slitless spectroscopy from a number of different space-based missions and instruments, notably HST/WFC3-IR, JWST/NIRISS, and the WFIRST wide-field instrument.

Here we take as an example the extremely deep WFC3 F140W imaging from the Hubble Ultra-Deep Field and processed by the "eXtreme Deep Field" project. We use SExtractor (Bertin & Arnouts 1996) to detect objects in the deep image, creating a catlog and an accompanying segmentation image that defines which pixels are assigned to each object.

The second half of the notebook demonstrates how to use the grizli.fake_image scripts to create WFC3, NIRISS, and WFIRST-sized cutouts extracted from the deep reference image. Those cutouts are then used as the reference to simulate slitless spectra from those instruments.

In [5]:
import os
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

import astropy.io.fits as pyfits
import astropy.wcs as pywcs
from astropy.table import Table
import astropy.time

import drizzlepac
import photutils

import grizli
import grizli.model
import grizli.fake_image

print('\ngrizli version: %s' %(grizli.__version__))
print('Now: ', astropy.time.Time.now().iso)

The following task in the stsci.skypac package can be run with TEAL:
                                    skymatch                                    
The following tasks in the drizzlepac package can be run with TEAL:
    astrodrizzle       config_testbed      imagefindpars           mapreg       
       photeq            pixreplace           pixtopix            pixtosky      
  refimagefindpars       resetbits          runastrodriz          skytopix      
     tweakback            tweakreg           updatenpol

grizli version: 1.3.2.dev1
Now:  2021-08-12 09:20:44.281
In [6]:
## Fetch the UDF images from the HLSP pages
workdir = '/Users/gbrammer/Research/Roman'
if os.path.exists(workdir):
    os.chdir(workdir)

if not os.path.exists('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits'):
    print('Fetch XDF images')
    url = 'https://archive.stsci.edu/missions/hlsp/xdf/'
    os.system(f'wget {url}/hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits')
    os.system(f'wget {url}/hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_wht.fits')

## SourceExtractor products
if not os.path.exists('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits'):
    print('Fetch demo XDF catalog files')
    os.system('wget https://s3.amazonaws.com/grizli-v1/Demo/xdf-demo-f140w.tar.gz')
    os.system('tar xzvf xdf-demo-f140w.tar.gz')

## Need a WFC3/IR FLT file later
if not os.path.exists('ibhj34h6q_flt.fits'):
    print('Fetch grizli demo files')
    os.system('wget http://www.stsci.edu/~brammer/grism/grizli_demo_data.tar.gz')
    os.system('tar xzvf grizli_demo_data.tar.gz')

WFC3 simulation and observed image

In [7]:
## Initialize the Grizli object, flt_file is the G141 exposure
pad=120 # allow simulation of objects at the edges
flt = grizli.model.GrismFLT(grism_file='ibhj34h8q_flt.fits', verbose=True, pad=pad,  
                            ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits', 
                            seg_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits')
Image cutout: x=slice(509, 4977, None), y=slice(882, 5250, None) [Out of range]
ibhj34h8q_flt.fits / blot reference hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits[0]
Using default C-based coordinate transformation...
Image cutout: x=slice(509, 4977, None), y=slice(882, 5250, None) [Out of range]
ibhj34h8q_flt.fits / Pad ref HDU with 559 pixels
ibhj34h8q_flt.fits / blot segmentation hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits[0]
Using default C-based coordinate transformation...
Using default C-based coordinate transformation...

The pad keyword above is used to extract a cutout from the reference image that is larger than the direct image itself so that objects can be accounted for and modeled that would otherwise "fall off" of the direct image but whose dispersed spectra fall on the grism images. The figure below compares the direct image itself to the "blotted" reference image that contains an extra border of pixels around the edges.

In [8]:
## Show the blotted reference direct image
f140w = pyfits.open('ibhj34h6q_flt.fits')
fig = plt.figure()

ax = fig.add_subplot(131) # FLT exposure
ax.imshow(f140w['SCI'].data, interpolation='Nearest', 
           origin='lower', vmin=-0.2, vmax=0.3, cmap='viridis')
ax.set_title('F140W FLT')

ax = fig.add_subplot(132) # Blotted reference image
blotted = flt.direct.data['REF']/flt.direct.ref_photflam
ax.imshow(blotted, interpolation='Nearest', 
           origin='lower', vmin=-0.2, vmax=0.3, cmap='viridis')
ax.set_title('"blotted" XDF reference')

ax = fig.add_subplot(133) # Grism
ax.imshow(f140w['SCI'].data - blotted[pad:-pad, pad:-pad],
           interpolation='Nearest', 
           origin='lower', vmin=-0.2, vmax=0.3, cmap='viridis')
ax.set_title('Difference')

for ax in fig.axes[1:]:
    ax.set_yticklabels([])

for i, ax in enumerate(fig.axes):
    ax.grid()
    offset = (i != 1)*pad # First is (1014,1014), others are (1014+2*pad, 1014+2*pad)
    ax.set_xlim(60-offset,360-offset)
    ax.set_ylim(60-offset,360-offset)

fig.tight_layout()

Since the grism configuration files are all defined in detector coordinates, the blot_catalog function is used to compute the detector positions of the objects detected in the rectified reference image.

In [9]:
## "blotted" SExtractor catalog, with catalog sky coordinates put into FLT frame
ref_cat = Table.read(pyfits.open('hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.cat')[2])
flt_cat = flt.blot_catalog(ref_cat, sextractor=True) # also stored in flt.catalog
In [10]:
## Compute flat-spectrum model of bright sources
mask = flt_cat['MAG_AUTO'] < 25
print('N=%d' %(mask.sum()))
#flt.compute_full_model(compute_beams=['A','B','C','D','E','F'], mask=mask, verbose=False)
flt.compute_full_model(ids=flt_cat['NUMBER'][mask], mags=flt_cat['MAG_AUTO'][mask], verbose=True)
N=515
515it [00:32, 15.96it/s]
In [11]:
# Compare to the actual G141 exposure
fig = plt.figure()

ax = fig.add_subplot(131)
ax.imshow(flt.grism.data['SCI'], interpolation='Nearest', 
           origin='lower', vmin=-0.01, vmax=0.1, cmap='gray_r')
ax.set_title('G141 exposure (1.1 ks)')

ax = fig.add_subplot(132)
ax.imshow(flt.model, interpolation='Nearest', 
           origin='lower', vmin=-0.01, vmax=0.1, cmap='gray_r')
ax.set_title(r'Flat $f_\lambda$ model')

ax = fig.add_subplot(133)
ax.imshow(flt.grism.data['SCI'] - flt.model, interpolation='Nearest', 
           origin='lower', vmin=-0.01, vmax=0.1, cmap='gray_r')
ax.set_title('Difference')

for ax in fig.axes[1:]:
    ax.set_yticklabels([])

# show lower corner and how objects are modeled right to the edge
for ax in fig.axes:
    ax.set_xlim(60,360)
    ax.set_ylim(60,360)

fig.tight_layout()

The bright spectrum at lower left and the fainter spectrum at the top would not be modeled if just the single direct image were available. However, they are included in the "padded" reference image and even the simple flat-spectrum model offers a reasonable first guess at the grism exposure model.

Slitless spectra from multi-missions JWST/NIRISS and WFIRST

The examples below show how to first generate a "dummy" cutout image from the reference image, where the image itself is not cutout but rather the dummy images have appropriate WCS so that the grizli.model.GrismFLT will then be able to blot the reference and segmentation images.

The image headers specify default values of the per-pixel sky backgrounds that are adopted with make_fake_image(background='auto'), that parameter can also be set to any desired float value. The make_fake_image generating script has parameters exptime and nexp that are the used to define a simple noise model for the resulting simulation.

$\sigma^2 = \mathrm{sky}\cdot\mathrm{exptime} + N_\mathrm{exp}\cdot\mathrm{readnoise}$

The value of $\sigma$ is stored in the ERR extension and random Gaussian deviates are put into the SCI extension of the output FITS files that are useful later for adding noise to the noiseless simulated spectra.

Again, these examples simulate grism exposures from the different missions/instruments with simple flat-spectrum SEDs with normalization set by the UDF F140W image (at 1.4 ┬Ám). However, it is straightforward to use any spectrum (i.e., stellar or galaxy templates) for any given object model.

In [12]:
### Fake images, cendered in the UDF/XDF
ra, dec = 53.1592277508136, -27.782056346146
pa_aper = 128.589

# allow simulation of objects at the edges
pad=0 # pixels

mag_limit = 25 # faint limit for the simulation

EXPTIME = 1.e4 # 10 ks ~ 4 HST orbits
NEXP = 10      # divided between 10 exposures

WFC3/IR G141

In [13]:
### WFC3/IR G141
h, wcs = grizli.fake_image.wfc3ir_header(filter='G141', ra=ra, dec=dec, pa_aper=pa_aper,
                                         flt='ibhj34h6q_flt.fits')

grizli.fake_image.make_fake_image(h, output='wfc3ir.fits', exptime=EXPTIME, nexp=NEXP, seed=1)

wfc3 = grizli.model.GrismFLT(grism_file='wfc3ir.fits', verbose=True, pad=pad,  
                           ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits', 
                           seg_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits')

wfc3_cat = wfc3.blot_catalog(ref_cat, sextractor=True) # also stored in flt.catalog
wfc3.catalog = wfc3_cat

mask = wfc3_cat['MAG_AUTO'] < mag_limit
print('N=%d' %(mask.sum()))

#wfc3.compute_full_model(compute_beams=['A','B','C','D','E','F'], mask=mask, verbose=False)
wfc3.compute_full_model(ids=wfc3_cat['NUMBER'][mask], mags=wfc3_cat['MAG_AUTO'][mask], verbose=True)
Image cutout: x=slice(1255, 4287, None), y=slice(1633, 4693, None)
wfc3ir.fits / blot reference hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits[0]
Using default C-based coordinate transformation...
Image cutout: x=slice(1255, 4287, None), y=slice(1633, 4693, None)
wfc3ir.fits / Pad ref HDU with 50 pixels
wfc3ir.fits / blot segmentation hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits[0]
Using default C-based coordinate transformation...
Using default C-based coordinate transformation...
N=432
2it [00:00, 13.23it/s]
ID 228 not found in segmentation image
78it [00:04, 15.42it/s]
ID 1504 not found in segmentation image
432it [00:26, 16.21it/s]

JWST NIRISS

(Seems to be broken in 2021 as the header doesn't work with the jwst pipeline data model)

In [14]:
try:
    import jwst
    h, wcs = grizli.fake_image.niriss_header(filter='F150W', ra=ra, dec=dec, pa_aper=pa_aper)
    grizli.fake_image.make_fake_image(h, output='niriss.fits', exptime=EXPTIME, nexp=NEXP)

    niriss = grizli.model.GrismFLT(grism_file='niriss.fits', verbose=True, pad=pad,  
                               ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits', 
                               seg_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits')

    niriss_cat = niriss.blot_catalog(ref_cat, sextractor=True) # also stored in flt.catalog
    niriss.catalog = niriss_cat

    mask = niriss_cat['MAG_AUTO'] < mag_limit
    print('N=%d' %(mask.sum()))
    #niriss.compute_full_model(compute_beams=['A','B','C','D','E'], mask=mask, verbose=False)
    niriss.compute_full_model(ids=niriss_cat['NUMBER'][mask], mags=niriss_cat['MAG_AUTO'][mask])
except:
    print('Skip JWST')
Skip JWST

Roman/G150 Grism

Note: The default simulation here is significantly deeper than the baseline High Latitude Survey

In [15]:
# Generate crude configuration files from Pandeia outputs
# Requires pandia.engine and the roman data files
from grizli import GRIZLI_PATH
conf_file = os.path.join(GRIZLI_PATH, 'CONF', 'Roman.G150.conf')

if not os.path.exists(conf_file):
    try:
        import pandeia
        %env pandeia_refdata=/Users/gbrammer/Research/Roman/pandeia_data-1.6_roman/
    
        grizli.fake_image.make_roman_config(save_to_conf=True)
    except:
        for file in ['Roman.G150.conf','Roman.G150.v1.6.sens.fits']:
            conf = os.path.join(grizli.GRIZLI_PATH, 'CONF')
            if not os.path.exists(os.path.join(conf, file)):
                print(file)
                os.system(f'wget https://s3.amazonaws.com/grizli-v1/Demo/{file} -O {conf}/{file}')
In [16]:
### Roman/G150 grism.  
h, wcs = grizli.fake_image.roman_header(ra=ra, dec=dec, pa_aper=pa_aper, naxis=(1180,1180))
grizli.fake_image.make_fake_image(h, output='roman.fits', exptime=EXPTIME, nexp=NEXP)

roman = grizli.model.GrismFLT(grism_file='roman.fits', verbose=True, pad=pad,  
                               ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits', 
                               seg_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits')

roman_cat = roman.blot_catalog(ref_cat, sextractor=True) 
roman.catalog = roman_cat

mask = roman_cat['MAG_AUTO'] < mag_limit
print('N=%d' %(mask.sum()))
#wfirst.compute_full_model(compute_beams=['A'], mask=mask, verbose=False)
roman.compute_full_model(ids=roman_cat['NUMBER'][mask], mags=roman_cat['MAG_AUTO'][mask], verbose=True)
Image cutout: x=slice(1275, 4315, None), y=slice(1666, 4706, None)
roman.fits / blot reference hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits[0]
Using default C-based coordinate transformation...
Image cutout: x=slice(1275, 4315, None), y=slice(1666, 4706, None)
roman.fits / Pad ref HDU with 50 pixels
roman.fits / blot segmentation hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits[0]
Using default C-based coordinate transformation...
Using default C-based coordinate transformation...
N=455
0it [00:00, ?it/s]
ID 85 not found in segmentation image
ID 133 not found in segmentation image
ID 198 not found in segmentation image
ID 301 not found in segmentation image
ID 791 not found in segmentation image
455it [00:02, 222.49it/s]

Roman HLS depth

In [17]:
# HLS exposure times

hdu, wcs = grizli.fake_image.roman_hls_image(ra=ra, dec=dec, pa_aper=pa_aper, naxis=(1180,1180), output='roman.fits')


hls = grizli.model.GrismFLT(grism_file='roman.fits', verbose=True, pad=pad,  
                               ref_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits', 
                               seg_file='hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits')

hls_cat = hls.blot_catalog(ref_cat, sextractor=True) 
hls.catalog = hls_cat

mask = hls_cat['MAG_AUTO'] < mag_limit
print('N=%d' %(mask.sum()))
#wfirst.compute_full_model(compute_beams=['A'], mask=mask, verbose=False)
hls.compute_full_model(ids=hls_cat['NUMBER'][mask], mags=hls_cat['MAG_AUTO'][mask], verbose=True)
Image cutout: x=slice(1275, 4315, None), y=slice(1666, 4706, None)
roman.fits / blot reference hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_sci.fits[0]
Using default C-based coordinate transformation...
Image cutout: x=slice(1275, 4315, None), y=slice(1666, 4706, None)
roman.fits / Pad ref HDU with 50 pixels
roman.fits / blot segmentation hlsp_xdf_hst_wfc3ir-60mas_hudf_f140w_v1_seg.fits[0]
Using default C-based coordinate transformation...
Using default C-based coordinate transformation...
N=455
0it [00:00, ?it/s]
ID 85 not found in segmentation image
ID 133 not found in segmentation image
ID 198 not found in segmentation image
ID 301 not found in segmentation image
29it [00:00, 280.21it/s]
ID 791 not found in segmentation image
455it [00:02, 219.32it/s]
In [18]:
### Show them!
# Compare to the actual G141 exposure
fig = plt.figure(figsize=[10,10.*2/3])

for i, sim, key in zip([0,1,2], [wfc3, roman, hls], ['WFC3','WFI','HLS']):
    # Direct
    axt = fig.add_subplot(231+i)
    axt.imshow(sim.direct.data['REF']/sim.direct.ref_photflam, interpolation='Nearest', 
           origin='lower', vmin=-0.1, vmax=0.2, cmap='viridis')
    axt.set_xticklabels([])
    
    # Grism
    axb = fig.add_subplot(234+i)
    axb.imshow(sim.model + sim.grism.data['SCI'], interpolation='Nearest', 
           origin='lower', vmin=-0.01, vmax=0.1, cmap='gray_r')
    axb.set_title('%s - %s' %(key, sim.grism.filter))
    
    if i > 0:
        axt.set_yticklabels([])
        axb.set_yticklabels([])

fig.tight_layout(pad=0.5)