NIRISS AMI Pipeline Notebook

stsci_logo

This notebook currently fails to execute, use as reference only

NIRISS AMI Pipeline Notebook#

Authors: R. Cooper
Last Updated: July 16, 2025
Pipeline Version: 1.19.1 (Build 12.0)

Purpose:
This notebook provides a framework for processing Near-Infrared Imager and Slitless Spectrograph (NIRISS) Aperture Masking Interferometry (AMI) data through all three James Webb Space Telescope (JWST) pipeline stages. Data is assumed to be located in one observation folder according to paths set up below. It should not be necessary to edit any cells other than in the Configuration section unless modifying the standard pipeline processing options.

Data: This notebook uses an example dataset from Program ID 1093 (PI: Thatte) which is the AMI commissioning program. For illustrative purposes, we will use a single target and reference star pair. Each exposure was taken in the F480W filter filter with the non-redundant mask (NRM) that enables AMI in the pupil. The observations used are observation 12 for the target and observation 15 for the reference star.

Example input data to use will be downloaded automatically unless disabled (i.e., to use local files instead).

JWST pipeline version and CRDS context:
This notebook was written for the above-specified pipeline version and associated build context for this version of the JWST Calibration Pipeline. Information about this and other contexts can be found in the JWST Calibration Reference Data System (CRDS server). If you use different pipeline versions, please refer to the table here to determine what context to use. To learn more about the differences for the pipeline, read the relevant documentation.

Please note that pipeline software development is a continuous process, so results in some cases may be slightly different if a subsequent version is used. For optimal results, users are strongly encouraged to reprocess their data using the most recent pipeline version and associated CRDS context, taking advantage of bug fixes and algorithm improvements. Any known issues for this build are noted in the notebook.

Updates:
This notebook is regularly updated as improvements are made to the pipeline. Find the most up to date version of this notebook at: https://github.com/spacetelescope/jwst-pipeline-notebooks/

Recent Changes:
March 31, 2025: original notebook released
July 16, 2025: Updated to jwst 1.19.1 (no significant changes)

Table of Contents#

  1. Configuration

  2. Package Imports

  3. Demo Mode Setup (ignore if not using demo data)

  4. Directory Setup

  5. Detector 1 Pipeline

  6. Image2 Pipeline

  7. AMI3 Pipeline

  8. Visualize the data

1. Configuration#

Install dependencies and parameters#

To make sure that the pipeline version is compatabile with the steps discussed below and the required dependencies and packages are installed, you can create a fresh conda environment and install the provided requirements.txt file:

conda create -n niriss_ami_pipeline python=3.11
conda activate niriss_ami_pipeline
pip install -r requirements.txt

Set the basic parameters to use with this notebook. These will affect what data is used, where data is located (if already in disk), and pipeline modules run in this data. The list of parameters are:

  • demo_mode

  • directories with data

  • pipeline modules

# Basic import necessary for configuration
import os
Note that demo_mode must be set appropriately below.

Set demo_mode = True to run in demonstration mode. In this mode this notebook will download example data from the Barbara A. Mikulski Archive for Space Telescopes (MAST) and process it through the pipeline. This will all happen in a local directory unless modified in Section 3 below.

Set demo_mode = False if you want to process your own data that has already been downloaded and provide the location of the data.

# Set parameters for demo_mode, channel, band, data mode directories, and 
# processing steps.

# -----------------------------Demo Mode---------------------------------
demo_mode = True

if demo_mode:
    print('Running in demonstration mode using online example data!')

# --------------------------User Mode Directories------------------------
# If demo_mode = False, look for user data in these paths
if not demo_mode:
    # Set directory paths for processing specific data; these will need
    # to be changed to your local directory setup (below are given as
    # examples)
    basedir = os.path.join(os.getcwd(), '')

    # Point to location of science observation data.
    # Assumes both science and PSF reference data are in the same directory
    # with uncalibrated data in sci_dir/uncal/ and results in stage1,
    # stage2, stage3 directories
    sci_dir = os.path.join(basedir, 'JWSTData/PID_1093/')

# Set which filter to process (empty will process all)
use_filter = '' # E.g., F480M

# --------------------------Set Processing Steps--------------------------
# Individual pipeline stages can be turned on/off here.  Note that a later
# stage won't be able to run unless data products have already been
# produced from the prior stage.

# Science processing
dodet1 = True  # calwebb_detector1
doimage2 = True  # calwebb_image2
doami3 = True  # calwebb_ami3
doviz = True  # Visualize calwebb_ami3 output

Set CRDS context and server#

Before importing CRDS and JWST modules, we need to configure our environment. This includes defining a CRDS cache directory in which to keep the reference files that will be used by the calibration pipeline.

If the root directory for the local CRDS cache directory has not been set already, it will be set to create one in the home directory.

# ------------------------Set CRDS context and paths----------------------
# Each version of the calibration pipeline is associated with a specific CRDS
# context file. The pipeline will select the appropriate context file behind
# the scenes while running. However, if you wish to override the default context
# file and run the pipeline with a different context, you can set that using
# the CRDS_CONTEXT environment variable. Here we show how this is done,
# although we leave the line commented out in order to use the default context.
# If you wish to specify a different context, uncomment the line below.
#os.environ['CRDS_CONTEXT'] = 'jwst_1322.pmap'  # CRDS context for 1.16.0

# Check whether the local CRDS cache directory has been set.
# If not, set it to the user home directory
if (os.getenv('CRDS_PATH') is None):
    os.environ['CRDS_PATH'] = os.path.join(os.path.expanduser('~'), 'crds')
    
# Check whether the CRDS server URL has been set.  If not, set it.
if (os.getenv('CRDS_SERVER_URL') is None):
    os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'

# Echo CRDS path in use
print(f"CRDS local filepath: {os.environ['CRDS_PATH']}")
print(f"CRDS file server: {os.environ['CRDS_SERVER_URL']}")

2. Package Imports#

# Use the entire available screen width for this notebook
from IPython.display import display, HTML
display(HTML("<style>.container { width:95% !important; }</style>"))

# Basic system utilities for interacting with files
# ----------------------General Imports------------------------------------
import glob
import time
import json
from pathlib import Path
from collections import defaultdict

# Numpy for doing calculations
import numpy as np

# -----------------------Astroquery Imports--------------------------------
# ASCII files, and downloading demo files
from astroquery.mast import Observations

# For visualizing data
import matplotlib.pyplot as plt
from astropy.visualization import (MinMaxInterval, SqrtStretch,
                                   ImageNormalize)

# For file manipulation
from astropy.io import fits
import asdf

# for JWST calibration pipeline
import jwst
import crds

from jwst.pipeline import Detector1Pipeline
from jwst.pipeline import Image2Pipeline
from jwst.pipeline import Ami3Pipeline

# JWST pipeline utilities
from jwst import datamodels
from jwst.associations import asn_from_list  # Tools for creating association files
from jwst.associations.lib.rules_level3_base import DMS_Level3_Base  # Definition of a Lvl3 association file

# Echo pipeline version and CRDS context in use
print(f"JWST Calibration Pipeline Version: {jwst.__version__}")
print(f"Using CRDS Context: {crds.get_context_name('jwst')}")

Define convenience functions#

# Define a convenience function to select only files of a given filter from an input set
def select_filter_files(files, use_filter):
    files_culled = []
    
    if (use_filter != ''):
        for file in files:
            model = datamodels.open(file)
            filt = model.meta.instrument.filter
            if (filt == use_filter):
                files_culled.append(file)
            model.close()
    else:
        files_culled = files
        
    return files_culled

# Define a convenience function to separate a list of input files into science and PSF reference exposures
def split_scipsf_files(files):
    psffiles = []
    scifiles = []

    for file in files:
        model = datamodels.open(file)
        if model.meta.exposure.psf_reference is True:
            psffiles.append(file)
        else:
            scifiles.append(file)
        model.close()
    
    return scifiles, psffiles

# Start a timer to keep track of runtime
time0 = time.perf_counter()

3. Demo Mode Setup (ignore if not using demo data)#

If running in demonstration mode, set up the program information to retrieve the uncalibrated data automatically from MAST using astroquery. MAST has a dedicated service for JWST data retrieval, so the archive can be searched by instrument keywords rather than just filenames or proposal IDs.

The list of searchable keywords for filtered JWST MAST queries is here.

For illustrative purposes, we will use a single target and reference star pair. Each exposure was taken in the F480W filter filter with the non-redundant mask (NRM) that enables AMI in the pupil.

We will start with uncalibrated data products. The files are named jw010930nn001_03102_00001_nis_uncal.fits, where nn refers to the observation number: in this case, observation 12 for the target and observation 15 for the reference star.

More information about the JWST file naming conventions can be found at: https://jwst-pipeline.readthedocs.io/en/latest/jwst/data_products/file_naming.html

# Set up the program information and paths for demo program
if demo_mode:
    print('Running in demonstration mode and will download example data from MAST!')
    
    # --------------Program and observation information--------------
    program = '01093'
    sci_observtn = ['012', '015']  # Obs 12 is the target, Obs 15 is the reference star
    visit = '001'
    visitgroup = '03'
    seq_id = "1"
    act_id = '02'
    expnum = '00001'

    # --------------Program and observation directories--------------
    data_dir = os.path.join('.', 'nis_ami_demo_data')
    sci_dir = os.path.join(data_dir, 'PID_1093')
    uncal_dir = os.path.join(sci_dir, 'uncal')  # Uncalibrated pipeline inputs should be here

    if not os.path.exists(uncal_dir):
        os.makedirs(uncal_dir)

    # Create directory if it does not exist
    if not os.path.isdir(data_dir):
        os.mkdir(data_dir)

Identify list of science (SCI) uncalibrated files associated with visits.

# Obtain a list of observation IDs for the specified demo program
if demo_mode:
    # Science data
    sci_obs_id_table = Observations.query_criteria(instrument_name=["NIRISS/AMI"],
                                                   proposal_id=[program],
                                                   filters=['F480M;NRM'],  # Data for Specific Filter
                                                   obs_id=['jw' + program + '*'])

# Turn the list of visits into a list of uncalibrated data files
if demo_mode:
    # Define types of files to select
    file_dict = {'uncal': {'product_type': 'SCIENCE',
                           'productSubGroupDescription': 'UNCAL',
                           'calib_level': [1]}}

    # Science files
    sci_files = []

    # Loop over visits identifying uncalibrated files that are associated
    # with them
    for exposure in (sci_obs_id_table):
        products = Observations.get_product_list(exposure)
        for filetype, query_dict in file_dict.items():
            filtered_products = Observations.filter_products(products, productType=query_dict['product_type'],
                                                             productSubGroupDescription=query_dict['productSubGroupDescription'],
                                                             calib_level=query_dict['calib_level'])
            sci_files.extend(filtered_products['dataURI'])

    # Select only the exposures we want to use based on filename
    # Construct the filenames and select files based on them
    filestrings = ['jw' + program + sciobs + visit + '_' + visitgroup + seq_id + act_id + '_' + expnum for sciobs in sci_observtn]
    sci_files_to_download = [scifile for scifile in sci_files if any(filestr in scifile for filestr in filestrings)]
    sci_files_to_download = sorted(set(sci_files_to_download))
    print(f"Science files selected for downloading: {len(sci_files_to_download)}")

Download all the uncal files and place them into the appropriate directories.

Warning: If this notebook is halted during this step the downloaded file may be incomplete, and cause crashes later on! 
if demo_mode:
    for filename in sci_files_to_download:
        sci_manifest = Observations.download_file(filename,
                                                  local_path=os.path.join(uncal_dir, Path(filename).name))

4. Directory Setup#

Set up detailed paths to input/output stages here.

# Define output subdirectories to keep science data products organized
# -----------------------------Science Directories------------------------------
uncal_dir = os.path.join(sci_dir, 'uncal')  # Uncalibrated pipeline inputs should be here
det1_dir = os.path.join(sci_dir, 'stage1')  # calwebb_detector1 pipeline outputs will go here
image2_dir = os.path.join(sci_dir, 'stage2')  # calwebb_image2 pipeline outputs will go here
ami3_dir = os.path.join(sci_dir, 'stage3')  # calwebb_ami3 pipeline outputs will go here

# We need to check that the desired output directories exist, and if not create them
# Ensure filepaths for input data exist
if not os.path.exists(uncal_dir):
    os.makedirs(uncal_dir)
 
if not os.path.exists(det1_dir):
    os.makedirs(det1_dir)
if not os.path.exists(image2_dir):
    os.makedirs(image2_dir)
if not os.path.exists(ami3_dir):
    os.makedirs(ami3_dir)

Print the exposure parameters of all potential input files:

uncal_files = sorted(glob.glob(os.path.join(uncal_dir, '*_uncal.fits')))
# Restrict to selected filter if applicable
uncal_files = select_filter_files(uncal_files, use_filter)

for file in uncal_files:
    model = datamodels.open(file)
    # print file name
    print(model.meta.filename)
    # Print out exposure info
    print(f"Instrument: {model.meta.instrument.name}")
    print(f"Filter: {model.meta.instrument.filter}")
    print(f"Pupil: {model.meta.instrument.pupil}")
    print(f"Number of integrations: {model.meta.exposure.nints}")
    print(f"Number of groups: {model.meta.exposure.ngroups}")
    print(f"Readout pattern: {model.meta.exposure.readpatt}")
    print(f"Dither position number: {model.meta.dither.position_number}")
    print("\n")
    model.close()

For the demo data, files should be for the NIRISS instrument using the F480M filter in the Filter Wheel and the NRM in the Pupil Wheel.

Likewise, both demo exposures use the NISRAPID readout pattern. The target has 5 groups per integration, and 69 integrations per exposure. The reference star has 12 groups per integration, and 61 integrations per exposure. They were taken at the same dither position; primary dither pattern position 1.

For more information about how JWST exposures are defined by up-the-ramp sampling, see the Understanding Exposure Times JDox article.

# Print out the time benchmark
time1 = time.perf_counter()
print(f"Runtime so far: {time1 - time0:0.0f} seconds")

5. Detector1 Pipeline#

Run the datasets through the Detector1 stage of the pipeline to apply detector level calibrations and create a countrate data product where slopes are fitted to the integration ramps. These *_rateints.fits products are 3D (nintegrations x nrows x ncols) and contain the fitted ramp slopes for each integration. 2D countrate data products (*_rate.fits) are also created (nrows x ncols) which have been averaged over all integrations.

By default, all steps in the Detector1 stage of the pipeline are run for NIRISS except: the ipc correction step and the gain_scale step. Note that while the persistence step is set to run by default, this step does not automatically correct the science data for persistence. The persistence step creates a *_trapsfilled.fits file which is a model that records the number of traps filled at each pixel at the end of an exposure. This file would be used as an input to the persistence step, via the input_trapsfilled argument, to correct a science exposure for persistence. Since persistence is not well calibrated for NIRISS, we do not perform a persistence correction and thus turn off this step to speed up calibration and to not create files that will not be used in the subsequent analysis. This step can be turned off when running the pipeline in Python by doing:

rate_result = Detector1Pipeline.call(uncal,steps={'persistence': {'skip': True}})

or as indicated in the cell bellow using a dictionary.

The charge_migration step is particularly important for NIRISS images to mitigate apparent flux loss in resampled images due to the spilling of charge from a central pixel into its neighboring pixels (see Goudfrooij et al. 2023 for details). Charge migration occurs when the accumulated charge in a central pixel exceeds a certain signal limit, which is ~25,000 ADU. This step is turned on by default for NIRISS imaging mode when using CRDS contexts of jwst_1159.pmap or later. Different signal limits for each filter are provided by the pars-chargemigrationstep parameter files. Users can specify a different signal limit by running this step with the signal_threshold flag and entering another signal limit in units of ADU. The effect is stronger when there is high contrast between a bright pixel and neighboring faint pixel, as is the case for the strongly peaked AMI PSF.

For AMI mode, preliminary investigation shows that dark subtraction does not improve calibration, and may in fact have a detrimental effect, so we turn it off here.

# Set up a dictionary to define how the Detector1 pipeline should be configured

# Boilerplate dictionary setup
det1dict = defaultdict(dict)

# Step names are copied here for reference
det1_steps = ['group_scale', 'dq_init', 'saturation', 'ipc', 'superbias', 'refpix',
              'linearity', 'persistence', 'dark_current', 'charge_migration',
              'jump', 'ramp_fit', 'gain_scale']

# Overrides for whether or not certain steps should be skipped
# skipping the ipc, persistence, and dark steps
det1dict['ipc']['skip'] = True
det1dict['persistence']['skip'] = True
det1dict['dark_current']['skip'] = True

# Overrides for various reference files
# Files should be in the base local directory or provide full path
#det1dict['dq_init']['override_mask'] = 'myfile.fits'  # Bad pixel mask
#det1dict['saturation']['override_saturation'] = 'myfile.fits'  # Saturation
#det1dict['linearity']['override_linearity'] = 'myfile.fits'  # Linearity
#det1dict['dark_current']['override_dark'] = 'myfile.fits'  # Dark current subtraction
#det1dict['jump']['override_gain'] = 'myfile.fits'  # Gain used by jump step
#det1dict['ramp_fit']['override_gain'] = 'myfile.fits'  # Gain used by ramp fitting step
#det1dict['jump']['override_readnoise'] = 'myfile.fits'  # Read noise used by jump step
#det1dict['ramp_fit']['override_readnoise'] = 'myfile.fits'  # Read noise used by ramp fitting step

# Turn on multi-core processing (off by default). Choose what fraction of cores to use (quarter, half, or all)
det1dict['jump']['maximum_cores'] = 'half'

# Alter parameters of certain steps (example)
#det1dict['charge_migration']['signal_threshold'] = X

The clean_flicker_noise step removes 1/f noise from calibrated ramp images, after the jump step and prior to performing the ramp_fitting step. By default, this step is skipped in the calwebb_detector1 pipeline for all instruments and modes. Although available, this step has not been extensively tested for the NIRISS AMI subarray and is thus not recommended at the present time.

Run Detector1 stage of pipeline

# Run Detector1 stage of pipeline, specifying:
# output directory to save *_rateints.fits files
# save_results flag set to True so the *rateints.fits files are saved
# save_calibrated_ramp set to True so *ramp.fits files are saved

if dodet1:
    for uncal in uncal_files:
        rate_result = Detector1Pipeline.call(uncal,
                                             output_dir=det1_dir,
                                             steps=det1dict,
                                             save_results=True,
                                             save_calibrated_ramp=True)
else:
    print('Skipping Detector1 processing')

# Print out the time benchmark
time1 = time.perf_counter()
print(f"Runtime for Detector1: {time1 - time0:0.0f} seconds")

Exploring the data#

Identify *_rateints.fits files and verify which pipeline steps were run and which calibration reference files were applied.

The header contains information about which calibration steps were completed and skipped and which reference files were used to process the data.

if dodet1:
    # find rateints files
    rateints_files = sorted(glob.glob(os.path.join(det1_dir, '*_rateints.fits')))
    # Restrict to selected filter if applicable
    rateints_files = select_filter_files(rateints_files, use_filter)
    
    # Read in the first file as datamodel as an example
    rateints = datamodels.open(rateints_files[0])
    
    # Check which steps were run
    for step, status in rateints.meta.cal_step.instance.items():
        print(f"{step}: {status}")

Check which CRDS version and reference files were used to calibrate the dataset:

if dodet1:
    for key, val in rateints.meta.ref_file.instance.items():
        print(f"{key}:")
        for param in rateints.meta.ref_file.instance[key]:
            print(f"\t{rateints.meta.ref_file.instance[key][param]}")

6. Image2 Pipeline#

In the Image2 stage of the pipeline, calibrated data products are created (*_cal.fits or *_calints.fits files, depending on whether the input files are *_rate.fits or *_rateints.fits).

In this pipeline processing stage, the data are flat fielded. For AMI, we do not perform any of the other Image2 steps such as photometric calibration, so our output data products still have units of countrate (ADU/s).

time_image2 = time.perf_counter()

# Set up a dictionary to define how the Image2 pipeline should be configured.

# Boilerplate dictionary setup
image2dict = defaultdict(dict)

image2steps = ['assign_wcs', 'flat_field', 'photom', 'resample']

# Overrides for whether or not certain steps should be skipped (example)
image2dict['photom']['skip'] = True
image2dict['resample']['skip'] = True

# Overrides for various reference files
# Files should be in the base local directory or provide full path
#image2dict['flat_field']['override_flat'] = 'myfile.fits'  # Pixel flatfield

Find and sort the input files, ensuring use of absolute paths:

# Use files from the detector1 output folder
rateints_files = sorted(glob.glob(os.path.join(det1_dir, 'jw*rateints.fits')))
# Restrict to selected filter if applicable
rateints_files = select_filter_files(rateints_files, use_filter)

for ii in range(len(rateints_files)):
    rateints_files[ii] = os.path.abspath(rateints_files[ii])

print(f"Found {str(len(rateints_files))} science files")

# Run Image2 stage of pipeline, specifying:
# output directory to save *_calints.fits files
# save_results flag set to True so the calints files are saved

if doimage2:
    for rateints in rateints_files:
        calints_result = Image2Pipeline.call(rateints,
                                             output_dir=image2_dir,
                                             steps=image2dict,
                                             save_results=True)
else:
    print("Skipping Image2 processing.")

# Print out the time benchmark
time1 = time.perf_counter()
print(f"Runtime so far: {time1 - time0:0.0f} seconds")
print(f"Runtime for Image2: {time1 - time_image2:0.0f} seconds")

Verify which pipeline steps were run:

if doimage2:
    # Identify *_calints.fits files
    calints_files = sorted(glob.glob(os.path.join(image2_dir, '*_calints.fits')))
    # Restrict to selected filter if applicable
    calints_files = select_filter_files(calints_files, use_filter)

    # Read in the first file as datamodel as an example
    calints = datamodels.open(calints_files[0])
    
    # Check which steps were run
    for step, status in calints.meta.cal_step.instance.items():
        print(f"{step}: {status}")

Check which reference files were used to calibrate the dataset:

if doimage2:
    for key, val in calints.meta.ref_file.instance.items():
        print(f"{key}:")
        for param in calints.meta.ref_file.instance[key]:
            print(f"\t{calints.meta.ref_file.instance[key][param]}")

7. AMI3 Pipeline#

In the AMI3 stage of the pipeline, the target and reference star *_calints.fits files are analyzed to extract interferometric observables, and then the target’s observables are normalized by those of the reference star to produce a final set of calibrated observables.

In order to run the AMI3 stage, an Association file needs to be created to inform the pipeline which exposures should be treated as targets and which as reference stars.

By default, the AMI3 stage of the pipeline performs the following steps on NIRISS data:

  • ami_analyze - fits a model to each integration of the input image and computes interferometric observables (fringe phases, amplitudes, and derived quantities).

  • ami_normalize - normalizes the target’s observables by the reference star observables.

While a previous version of the AMI3 pipeline included an intermediate step to average together the observables from multiple exposures (ami_average), the step is not currently supported.

time_ami3 = time.perf_counter()

Create ASDF File#

The ami_analyze step has several optional arguments, some can be provided in an ASDF file containing a particular data tree. In here we will modify this file to provide an specific affine distortion matrix (affine2d).

For affine2d the default parameters from commissioning are accessed with a special string: ‘commissioning.’ If affine2d = None, it will perform a search for the best-fit rotation only. To use a different affine distortion, such as one measured directly from the data or an ideal distortion, it must be specified in an ASDF file as shown here.

# Example of writing an ASDF file specifying an affine distortion matrix to use

# Create an ASDF file
asdf_affine2d = os.path.join(sci_dir, 'affine2d_ideal.asdf')

aff_tree = {'mx': 1.,  # dimensionless x-magnification
            'my': 1.,  # dimensionless y-magnification
            'sx': 0.,  # dimensionless x shear
            'sy': 0.,  # dimensionless y shear
            'xo': 0.,  # x-offset in pupil space
            'yo': 0.,  # y-offset in pupil space
            'rotradccw': None}

with open(asdf_affine2d, 'wb') as fh:
    af = asdf.AsdfFile(aff_tree)
    af.write_to(fh)
    
af.close()

An example of using the affine distortion ASDF file created above is shown below.

# Set up a dictionary to define how the AMI3 pipeline should be configured
# Boilerplate dictionary setup
ami3dict = defaultdict(dict)

# Options for ami_analyze step
#ami3dict['ami_analyze']['firstfew'] = 5  # Analyze only the first 5 integrations to speed up demo
#ami3dict['ami_analyze']['affine2d'] = 11  # Increase oversampling of image plane fit (increases runtime)
#ami3dict['ami_analyze']['bandpass'] = 'myfile.asdf'  # Provide a custom bandpass (e.g., from synphot)
ami3dict['ami_analyze']['affine2d'] = asdf_affine2d  # Use the affine distortion ASDF file we created
ami3dict['ami_analyze']['save_results'] = True  # Turn on optional output results for display purposes

# Overrides for whether or not certain steps should be skipped (example)
#ami3dict['ami_normalize']['skip'] = True

# Overrides for various reference files
# Files should be in the base local directory or provide full path
#ami3dict['ami_analyze']['override_nrm'] = 'mynrmfile.fits'  # NRM mask geometry file

Find and sort all of the input files, ensuring use of absolute paths

# AMI3 takes the calints.fits files
calints_files = sorted(glob.glob(os.path.join(image2_dir, 'jw*calints.fits')))
# Restrict to selected filter if applicable
calints_files = select_filter_files(calints_files, use_filter)
    
calints_files = [os.path.abspath(calints) for calints in calints_files]

print(f'Found {str(len(calints_files))} science files to process')

Create Association Files#

An association file lists the exposures to calibrate together in Stage 3 of the pipeline. Note that an association file is available for download from MAST, with a filename of *_asn.json, though it may require additional manipulation for AMI.

Here we create association files for each pairing of input science and reference PSF exposures. Note that the final output products will have a rootname that is specified by the product_name in the association file.

# Create Level 3 Associations
if doami3:
    # Separate inputs into science and PSF reference exposures
    scifiles, psffiles = split_scipsf_files(calints_files)
    
    # Make an association file for every valid science/psf file pair
    for sci in scifiles:
        hdr = fits.getheader(sci)
        thisfilter = hdr['FILTER']
        # Potential PSF calibration files are those for which filter matches
        psfoptions = select_filter_files(psffiles, thisfilter)
        
        for psf in psfoptions:
            # Construct product name from the input headers
            prodname = 'ami3_' + hdr['TARGPROP'] + '_' + hdr['FILTER'] + '_' + hdr['ACT_ID']
            prodname = prodname.replace(" ", "")
    
            asn = asn_from_list.asn_from_list([sci],
                                              rule=DMS_Level3_Base,
                                              product_name=prodname)
    
            # Set the second observation as the psf reference star
            asn['products'][0]['members'].append({'expname': psf, 'exptype': 'psf'})
            asn.data['asn_type'] = 'ami3'
            asn.data['program'] = hdr['PROGRAM']    
    
            # Format association as .json file
            asn_filename = prodname + '_asn.json'
            _, serialized = asn.dump(format="json")

            # Write out association file
            association_ami3 = os.path.join(sci_dir, asn_filename)
            with open(association_ami3, "w") as fd:
                fd.write(serialized)

    # List all associations            
    all_asn = sorted(glob.glob(os.path.join(sci_dir, '*_asn.json')))

Check that file paths have been correctly updated.

# Open an ASN file as an example.
if doami3:
    if isinstance(all_asn, str):
        with open(all_asn, 'r') as f_obj:
            asnfile_data = json.load(f_obj)
    elif isinstance(all_asn, list):
        with open(all_asn[0], 'r') as f_obj:
            asnfile_data = json.load(f_obj)
    expanded_json = json.dumps(asnfile_data, indent=2)  # 'expanded=True' maps to 'indent'
    print(expanded_json)

Run AMI3 stage of the pipeline#

For the target and reference star exposures listed in the association file, the AMI3 stage of the pipeline will produce:

  • *ami-oi.fits files (one for the target, one for the reference star) from the ami_analyze step containing averaged interferometric observables

  • *aminorm-oi.fits file from the ami_normalize step containing normalized interferometric observables

The *ami-oi.fits and *aminorm-oi.fits files adhere to the OIFITS2 format, which is a registered FITS standard for optical and infrared interferometry data.

# Run Stage 3
if doami3:
    for asn in all_asn:
        ami3_result = Ami3Pipeline.call(asn,
                                        output_dir=ami3_dir,
                                        steps=ami3dict)
else:
    print('Skipping AMI3 processing')

# Print out the time benchmark
time1 = time.perf_counter()
print(f"Runtime so far: {time1 - time0:0.0f} seconds")
print(f"Runtime for AMI3: {time1 - time_ami3:0.0f} seconds")

Verify which pipeline steps were run#

if doami3:
    # Identify *ami-oi.fits file and open as datamodel
    ami_oi = glob.glob(os.path.join(ami3_dir, "*_ami-oi.fits"))[0]

    with datamodels.open(ami_oi) as oi_f:
        # Check which steps were run
        for step, status in oi_f.meta.cal_step.instance.items():
            print(f"{step}: {status}")

Check which reference files were used to calibrate the dataset

if doami3:
    for key, val in oi_f.meta.ref_file.instance.items():
        print(f"{key}:")
        for param in oi_f.meta.ref_file.instance[key]:
            print(f"\t{oi_f.meta.ref_file.instance[key][param]}")

8. Visualize the results#

Plot interferometric observables#

We will now plot the interferometric observables for the target, reference star, and normalized target.

# Define a function for plotting observables

def plot_observables(ami_oimodel):
    # Read the observables from the datamodel
    # Squared visibilities and uncertainties
    vis2 = ami_oimodel.vis2["VIS2DATA"]
    vis2_err = ami_oimodel.vis2["VIS2ERR"]
    # Closure phases and uncertainties
    t3phi = ami_oimodel.t3["T3PHI"]
    t3phi_err = ami_oimodel.t3["T3PHIERR"]

    # Construct baselines between the U and V coordinates of sub-apertures
    baselines = (ami_oimodel.vis2['UCOORD']**2 + ami_oimodel.vis2['VCOORD']**2)**0.5

    # Construct baselines between combinations of three sub-apertures
    u1 = ami_oimodel.t3['U1COORD']
    u2 = ami_oimodel.t3['U2COORD']
    v1 = ami_oimodel.t3['V1COORD']
    v2 = ami_oimodel.t3['V2COORD']
    u3 = -(u1 + u2)
    v3 = -(v1 + v2)
    baselines_t3 = []
    for k in range(len(u1)):
        B1 = np.sqrt(u1[k]**2 + v1[k]**2)
        B2 = np.sqrt(u2[k]**2 + v2[k]**2)
        B3 = np.sqrt(u3[k]**2 + v3[k]**2)
        # Use longest baseline of the three for plotting
        baselines_t3.append(np.max([B1, B2, B3])) 
    baselines_t3 = np.array(baselines_t3)

    # Plot closure phases, squared visibilities against their baselines 
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
    ax1.errorbar(baselines_t3, t3phi, yerr=t3phi_err, fmt="go")
    ax2.errorbar(baselines, vis2, yerr=vis2_err, fmt="go")
    ax1.set_xlabel(r"$B_{max}$ [m]", size=12)
    ax1.set_ylabel("Closure phase [deg]", size=12)
    ax1.set_title("Closure Phase", size=14)
    ax2.set_title("Squared Visibility", size=14)
    ax2.set_xlabel(r"$B_{max}$ [m]", size=12)
    ax2.set_ylabel("Squared Visibility", size=12)
    plt.suptitle(ami_oimodel.meta.filename, fontsize=16)
    ax1.set_ylim([-3.5, 3.5])
    ax2.set_ylim([0.5, 1.1])

First we’ll look at the target and reference star’s squared visibilities and closure phases. We plot the squared visibilities against their corresponding baselines (units of meters projected back to the primary meter). We plot the closure phases against the longest of the three baselines that forms the “closure triangle” whose phases were summed to calculate the closure phase.

if doviz:
    # Get the first target ami-oi.fits file
    oi_scifile = sorted(glob.glob(os.path.join(ami3_dir, "jw*_ami-oi.fits")))
    if isinstance(oi_scifile, list):
        oi_scifile = oi_scifile[0]
        
    # Get the first PSF reference ami-oi fits file
    oi_psffile = sorted(glob.glob(os.path.join(ami3_dir, "jw*psf-ami-oi.fits")))
    if isinstance(oi_psffile, list):
        oi_psffile = oi_psffile[0]

    # Open them as datamodels
    amioi_targ = datamodels.open(oi_scifile)
    amioi_ref = datamodels.open(oi_psffile)

    # Plot the observables
    plot_observables(amioi_targ)
    plot_observables(amioi_ref)

The top two scatter plots show the observables from the target observation, and the lower two show the observables from the reference star. For a perfect point source observation at single wavelength, we would expect to recover closure phases of zero and squared visibilites of unity.

Since the closure phases are sensitive to asymmetries in the source brightness distribution, we can tell from the larger scatter of the target that it is likely not a point source; i.e, there is a faint companion. On the other hand, the reference star has closure phases with a much smaller scatter around zero. The squared visibilities of both decrease at longer wavelengths due to an effect called bandpass smearing. We expect that calibrating the target by the reference star should correct for this, as well as other systematics.

Now we will plot the final calibrated data product; the target normalized by the reference star. The reference star closure phases are subtracted from the target closure phases, and the target squared visibilities are divided by the reference star squared visibilities.

Further scientific analysis on these calibrated OIFITS files can be done with community-developed analysis software like CANDID (Gallenne et al. 2015) or Fouriever to extract binary parameters, or an image reconstruction code like SQUEEZE (Baron et al. 2010) or BSMEM (Skilling & Bryan 1984, Buscher 1994, Baron & Young 2008).

if doviz:
    # Identify calibrated *_aminorm-oi.fits file and open as datamodel
    abdor_oifits = glob.glob(os.path.join(ami3_dir, "*_aminorm-oi.fits"))[0]
    amioi_norm = datamodels.open(abdor_oifits)

    # Plot the observables
    plot_observables(amioi_norm)

Display the best-fit model#

We can also look at the cleaned data, model, and residual images that are saved in the auxiliary *amilg.fits data products:

if doviz:
    # Find the data files
    amilg = sorted(glob.glob(os.path.join(ami3_dir, '*amilg.fits')))

    # Open the first one as an AmiLGModel
    firstfile = amilg[0]
    amilgmodel = datamodels.open(firstfile)

    # Plot the data, model, residual
    norm = ImageNormalize(amilgmodel.norm_centered_image[0], 
                          interval=MinMaxInterval(), 
                          stretch=SqrtStretch())
    fig, axs = plt.subplots(1, 3, figsize=(12, 5))
    axs[0].set_title('Normalized Data')
    im1 = axs[0].imshow(amilgmodel.norm_centered_image[0], norm=norm)
    axs[1].set_title('Normalized Model')
    im2 = axs[1].imshow(amilgmodel.norm_fit_image[0], norm=norm)
    axs[2].set_title('Normalized Residual (Data-Model)')
    im3 = axs[2].imshow(amilgmodel.norm_resid_image[0])

    for im in [im1, im2, im3]:
        plt.colorbar(im, shrink=.95, location='bottom', pad=.05)
    for ax in axs:
        ax.axis('off')
        
    plt.suptitle(os.path.basename(firstfile))
    plt.tight_layout()

Each image has been normalized by the peak pixel value of the data, and the data and model are displayed on a square root stretch to emphasize the fainter features. By looking at the residual (data - model) image, we can see that the model is a good fit to the data. This model achieves better contrast than ground-based NRM, but has not reached the binary contrast science requirements of AMI. The faint vertical striping in the background of the residual image is 1/f noise (flicker noise), which is an active area of improvement for the NIRISS/AMI team, as is the best method of correcting for charge migration.

stsci_logo
Contents