Source code for jwql.utils.instrument_properties

"""Collection of functions dealing with retrieving/calculating various
instrument properties

Authors
-------

    - Bryan Hilbert

Uses
----

    This module can be imported and used as such:

    ::

        from jwql.utils import instrument_properties as inst
        amps = inst.amplifier_info('my_files.fits')
"""

from copy import deepcopy
import datetime

from astropy.io import fits
from jwst.datamodels import dqflags
import numpy as np

from jwql.utils.constants import AMPLIFIER_BOUNDARIES, FOUR_AMP_SUBARRAYS, NIRCAM_SUBARRAYS_ONE_OR_FOUR_AMPS


[docs]def amplifier_info(filename, omit_reference_pixels=True): """Calculate the number of amplifiers used to collect the data in a given file using the array size and exposure time of a single frame (This is needed because there is no header keyword specifying how many amps were used.) Parameters ---------- filename : str Name of fits file to investigate omit_reference_pixels : bool If ``True``, return the amp boundary coordinates excluding reference pixels Returns ------- num_amps : int Number of amplifiers used to read out the data amp_bounds : dict Dictionary of amplifier boundary coordinates. Keys are strings of the amp number (1-4). Each value is a list composed of two tuples. The first tuple gives the coordinates of the (minimum x, maximum x, and x step value), and the second tuple gives the (minimum y, maximum y, and y step value). These are set up such that a list of indexes for each amplifier can be generated by using ``np.mgrid[x_min: x_max: x_step, y_min: y_max: y_step]`` """ # First get necessary metadata header = fits.getheader(filename) instrument = header['INSTRUME'].lower() detector = header['DETECTOR'] x_dim = header['SUBSIZE1'] y_dim = header['SUBSIZE2'] sample_time = header['TSAMPLE'] * 1.e-6 frame_time = header['TFRAME'] subarray_name = header['SUBARRAY'] aperture = "{}_{}".format(detector, subarray_name) # Full frame data will be 2048x2048 for all instruments if instrument.lower() == 'miri' or ((x_dim == 2048) and (y_dim == 2048)) or \ subarray_name in FOUR_AMP_SUBARRAYS: num_amps = 4 amp_bounds = deepcopy(AMPLIFIER_BOUNDARIES[instrument]) else: if subarray_name not in NIRCAM_SUBARRAYS_ONE_OR_FOUR_AMPS: num_amps = 1 amp_bounds = {'1': [(0, x_dim, 1), (0, y_dim, 1)]} else: # These are the tougher cases. Subarrays that can be # used with multiple amp combinations # Compare the given frametime with the calculated frametimes # using 4 amps or 1 amp. # Right now this is used only for the NIRCam grism stripe # subarrays, so we don't need this to be a general case that # can handle any subarray orientation relative to any amp # orientation amp4_time = calc_frame_time(instrument, aperture, x_dim, y_dim, 4, sample_time=sample_time) amp1_time = calc_frame_time(instrument, aperture, x_dim, y_dim, 1, sample_time=sample_time) if np.isclose(amp4_time, frame_time, atol=0.001, rtol=0): num_amps = 4 # In this case, keep the full frame amp boundaries in # the x direction, and set the boundaries in the y # direction equal to the height of the subarray amp_bounds = deepcopy(AMPLIFIER_BOUNDARIES[instrument]) for amp_num in ['1', '2', '3', '4']: newdims = (amp_bounds[amp_num][1][0], y_dim, 1) amp_bounds[amp_num][1] = newdims elif np.isclose(amp1_time, frame_time, atol=0.001, rtol=0): num_amps = 1 amp_bounds = {'1': [(0, x_dim, 1), (0, y_dim, 1)]} else: raise ValueError(('Unable to determine number of amps used for exposure. 4-amp frametime' 'is {}. 1-amp frametime is {}. Reported frametime is {}.') .format(amp4_time, amp1_time, frame_time)) if omit_reference_pixels: # If requested, ignore reference pixels by adjusting the indexes of # the amp boundaries. with fits.open(filename) as hdu: try: data_quality = hdu['DQ'].data except KeyError: try: data_quality = hdu['PIXELDQ'].data except KeyError: raise KeyError('DQ extension not found.') # Reference pixels should be flagged in the DQ array with the # REFERENCE_PIXEL flag. Find the science pixels by looping for # pixels that don't have that bit set. scipix = np.where(data_quality & dqflags.pixel['REFERENCE_PIXEL'] == 0) ymin = np.min(scipix[0]) xmin = np.min(scipix[1]) ymax = np.max(scipix[0]) + 1 xmax = np.max(scipix[1]) + 1 # Adjust the minimum and maximum x and y values if they are within # the reference pixels for i, key in enumerate(amp_bounds): bounds = amp_bounds[key] prev_xmin, prev_xmax, prev_xstep = bounds[0] prev_ymin, prev_ymax, prev_ystep = bounds[1] if prev_xmin < xmin: # Ensure the x start values remain interleaved for MIRI if instrument.lower() == 'miri': new_xmin = xmin + i else: new_xmin = xmin else: new_xmin = prev_xmin if prev_ymin < ymin: new_ymin = ymin else: new_ymin = prev_ymin if prev_xmax > xmax: new_xmax = xmax else: new_xmax = prev_xmax if prev_ymax > ymax: new_ymax = ymax else: new_ymax = prev_ymax amp_bounds[key] = [(new_xmin, new_xmax, prev_xstep), (new_ymin, new_ymax, prev_ystep)] return num_amps, amp_bounds
[docs]def calc_frame_time(instrument, aperture, xdim, ydim, amps, sample_time=1.e-5): """Calculate the readout time for a single frame of a given size and number of amplifiers. Note that for NIRISS and FGS, the fast readout direction is opposite to that in NIRCam, so we switch ``xdim`` and ``ydim`` so that we can keep a single equation. Parameters ---------- instrument : str Name of the instrument being simulated aperture : str Name of aperture being simulated (e.g ``NRCA1_FULL``). Currently this is only used to check for the FGS ``ACQ1`` aperture, which uses a unique value of ``colpad`` below. xdim : int Number of columns in the frame ydim : int Number of rows in the frame amps : int Number of amplifiers used to read out the frame sample_time : float Time to sample a pixel, in seconds. For NIRCam/NIRISS/FGS this is 10 microseconds = 1e-5 seconds Returns ------- frametime : float Readout time in seconds for the frame """ instrument = instrument.lower() if instrument == "nircam": colpad = 12 # Fullframe if amps == 4: rowpad = 1 fullpad = 1 else: # All subarrays rowpad = 2 fullpad = 0 if ((xdim <= 8) & (ydim <= 8)): # The smallest subarray rowpad = 3 elif instrument == "niriss": colpad = 12 # Fullframe if amps == 4: rowpad = 1 fullpad = 1 else: rowpad = 2 fullpad = 0 elif instrument == 'fgs': colpad = 6 if 'acq1' in aperture.lower(): colpad = 12 rowpad = 1 if amps == 4: fullpad = 1 else: fullpad = 0 return ((1.0 * xdim / amps + colpad) * (ydim + rowpad) + fullpad) * sample_time
[docs]def get_obstime(filename): """Extract the observation date and time from a fits file Parameters ---------- filename : str Name of fits file Returns ------- obs_time : datetime.datetime Observation date and time """ with fits.open(filename) as h: date = h[0].header['DATE-OBS'] time = h[0].header['TIME-OBS'] year, month, day = [int(element) for element in date.split('-')] hour, minute, second = [float(element) for element in time.split(':')] return datetime.datetime(year, month, day, int(hour), int(minute), int(second))
[docs]def mean_time(times): """Given a list of datetime objects, calculate the mean time Parameters ---------_ times : list List of datetime objects Returns ------- meantime ; datetime.datetime Mean time of the input ``times`` """ seconds_per_day = 24. * 3600. min_time = np.min(times) delta_times = [] for time in times: delta_time = time - min_time total_seconds = delta_time.days * seconds_per_day + delta_time.seconds delta_times.append(total_seconds) mean_delta = np.mean(delta_times) mean_time_delta = datetime.timedelta(0, mean_delta, 0) return min_time + mean_time_delta