Source code for sofia_redux.instruments.forcast.imgshift_header

# Licensed under a 3-clause BSD style license - see LICENSE.rst

from astropy import log
from astropy.io import fits
import numpy as np

from sofia_redux.instruments.forcast.getpar import getpar

__all__ = ['imgshift_header']


[docs] def imgshift_header(header, chop=True, nod=True, dither=True, dripconf=True): """ Calculates the shift_image in the pixel frame for merging an image Reads the chop/nod/dither amplitudes, angles, and transform from the sky to the detector frame depending on the coordinate system. The chop/nod angle was defined in the opposite sense in flights earlier than 103. This distinction is recorded in the ANGLCONV keyword value in the header. If "positive", as for 99 and earlier, use the angle as written in the header; otherwise, angles will need to be converted to negative values. The chop, nod, and dither coordinate systems (CHPCRSYS, NODCRSYS, DTHCRSYS in the header) may take values of SIRF (science instrument reference frame), TARF (telescope assembly reference frame), or ERF (equatorial reference frame). If these keywords are not available in the header, we check the CHPCOORD, NODCOORD, and DTHRCS) keywords. The COORD/CS keyword values are integers where 0=SIRF, 1=TARF, and 2=ERF. If the coordinate system cannot be be determined, SIRF is used by default. TARF is seen as equivalent to SIRF for the purposes of this algorithm. Parameters ---------- header : astropy.io.fits.header.Header A FITS header chop : bool, optional Indicates whether to perform the calculation of chop shift_image nod : bool, optional Indicates whether to perform the calculation of nod shift_image dither : bool, optional Indicates whether to perform the calculation of dither shift_image dripconf : bool, optional Indicates whether configuration values can overwrite header keyword values. Returns ------- dict A dict containing the X and Y shifts necessary to merge or coadd an image based on the header keywords. """ shift = {'chopx': 0., 'chopy': 0., 'chopcoord': '', 'nodx': 0., 'nody': 0., 'nodcoord': '', 'ditherx': 0., 'dithery': 0., 'dithercoord': '', 'sky_angle': 0.} if not isinstance(header, fits.header.Header): log.error("header is not %s" % fits.header.Header) return shift # Get plate scale - assume square pixels # may be in degrees, in pltscale pixsize = header.get('PLTSCALE', None) if pixsize is None: telescope = header.get('TELESCOP', 'SOFIA').strip().upper() if telescope == 'PIXELS': plate_scale = [1.0, 1.0] else: # assume forcast default plate_scale = [0.768, 0.768] else: pixsize *= 3600. plate_scale = [pixsize, pixsize] sky_angle = getpar(header, 'SKY_ANGL', dripconf=dripconf, dtype=float, default=0.0) shift['sky_angle'] = sky_angle sky_angle = np.radians(sky_angle) cosa, sina = np.cos(sky_angle), np.sin(sky_angle) skyrot = np.array(((cosa, -sina), (sina, cosa))).T # Determine angle convention +1/-1 angle_convention = getpar( header, 'ANGLCONV', dripconf=dripconf, dtype=str, default='UNKNOWN', comment="Chopping/Nodding angle conventions").strip().lower() angle_convention = 1 if angle_convention == 'positive' else -1 # Get chop distances if int(header.get('CHOPPING', 0)) and chop: dchop = getpar(header, 'CHPAMP1', dtype=float, default=0) * 2 chop_angle = getpar(header, 'CHPANGLE', dtype=float, default=0, dripconf=dripconf) chop_angle = np.radians(chop_angle) * angle_convention shift['chopx'] = dchop * np.sin(chop_angle) / plate_scale[0] shift['chopy'] = dchop * np.cos(chop_angle) / plate_scale[1] coordsys = getpar(header, 'CHPCRSYS', dripconf=dripconf, default='unknown', dtype=str).upper().strip() if coordsys not in ['ERF', 'SIRF']: val = getpar(header, 'CHPCOORD', default=0, dtype=int, dripconf=dripconf) coordsys = 'ERF' if val == 2 else 'SIRF' if coordsys == 'ERF': shift['chopx'], shift['chopy'] = \ skyrot @ (shift['chopx'], shift['chopy']) shift['chopcoord'] = coordsys # Get nod distances if int(header.get('NODDING', 0)) and nod: dnod = getpar(header, 'NODAMP', default=0, dtype=float, dripconf=dripconf) nod_angle = getpar(header, 'NODANGLE', default=0, dtype=float, dripconf=dripconf) nod_angle = np.radians(nod_angle) * angle_convention shift['nodx'] = dnod * np.sin(nod_angle) / plate_scale[0] shift['nody'] = dnod * np.cos(nod_angle) / plate_scale[1] coordsys = getpar(header, 'NODCRSYS', dripconf=dripconf, default='unknown', dtype=str).upper().strip() if coordsys not in ['ERF', 'SIRF']: val = getpar(header, 'NODCOORD', default=0, dtype=int, dripconf=dripconf) coordsys = 'ERF' if val == 2 else 'SIRF' if coordsys == 'ERF': shift['nodx'], shift['nody'] = \ skyrot @ (shift['nodx'], shift['nody']) shift['nodcoord'] = coordsys # Get dither distances if int(header.get('DITHER', 0)) and dither: x = getpar(header, 'DITHERX', dtype=float, default=0, dripconf=dripconf) y = getpar(header, 'DITHERY', dtype=float, default=0, dripconf=dripconf) # from newest coordinate system header keyword to oldest for dthcrsys in ['DTHCRSYS', 'DTHRCS', 'DITHERCS']: if dthcrsys in header: coordsys = str(header.get(dthcrsys)).strip().upper() break else: coordsys = 'unknown' if coordsys not in ['ERF', 'SIRF']: coordsys = 'ERF' if coordsys == '2' else 'SIRF' if coordsys == 'SIRF': rot_angle = np.pi - sky_angle cosa, sina = np.cos(rot_angle), np.sin(rot_angle) # clockwise rotation sirf_x = (x * cosa) + (y * sina) sirf_y = (-x * sina) + (y * cosa) # x dither in the SIRF coorsys has the opposite sign x, y = sirf_x, sirf_y else: # ERF coordsys as dither values in arcseconds x /= -plate_scale[0] y /= plate_scale[1] shift['ditherx'], shift['dithery'] = x, y shift['dithercoord'] = coordsys return shift