diff options
Diffstat (limited to 'stwcs/wcsutil/headerlet.py')
| -rw-r--r-- | stwcs/wcsutil/headerlet.py | 2754 |
1 files changed, 2754 insertions, 0 deletions
diff --git a/stwcs/wcsutil/headerlet.py b/stwcs/wcsutil/headerlet.py new file mode 100644 index 0000000..c0dd9b0 --- /dev/null +++ b/stwcs/wcsutil/headerlet.py @@ -0,0 +1,2754 @@ +""" +This module implements headerlets. + +A headerlet serves as a mechanism for encapsulating WCS information +which can be used to update the WCS solution of an image. The idea +came up first from the desire for passing improved astrometric +solutions for HST data and provide those solutions in a manner +that would not require getting entirely new images from the archive +when only the WCS information has been updated. + +""" + +from __future__ import absolute_import, division, print_function +import os +import sys +import functools +import logging +import textwrap +import copy +import time + +import numpy as np +from astropy.io import fits +#import pywcs +from astropy import wcs as pywcs +from astropy.utils import lazyproperty + +from stsci.tools.fileutil import countExtn +from stsci.tools import fileutil as fu +from stsci.tools import parseinput + +from stwcs.updatewcs import utils +from . import altwcs +from . import wcscorr +from .hstwcs import HSTWCS +from .mappings import basic_wcs + +#### Logging support functions +class FuncNameLoggingFormatter(logging.Formatter): + def __init__(self, fmt=None, datefmt=None): + if '%(funcName)s' not in fmt: + fmt = '%(funcName)s' + fmt + logging.Formatter.__init__(self, fmt=fmt, datefmt=datefmt) + + def format(self, record): + record = copy.copy(record) + if hasattr(record, 'funcName') and record.funcName == 'init_logging': + record.funcName = '' + else: + record.funcName += ' ' + return logging.Formatter.format(self, record) + + +logger = logging.getLogger(__name__) +formatter = FuncNameLoggingFormatter("%(levelname)s: %(message)s") +ch = logging.StreamHandler() +ch.setFormatter(formatter) +ch.setLevel(logging.CRITICAL) +logger.addHandler(ch) +logger.setLevel(logging.DEBUG) + +FITS_STD_KW = ['XTENSION', 'BITPIX', 'NAXIS', 'PCOUNT', + 'GCOUNT', 'EXTNAME', 'EXTVER', 'ORIGIN', + 'INHERIT', 'DATE', 'IRAF-TLM'] + +DEFAULT_SUMMARY_COLS = ['HDRNAME', 'WCSNAME', 'DISTNAME', 'AUTHOR', 'DATE', + 'SIPNAME', 'NPOLFILE', 'D2IMFILE', 'DESCRIP'] +COLUMN_DICT = {'vals': [], 'width': []} +COLUMN_FMT = '{:<{width}}' + + +def init_logging(funcname=None, level=100, mode='w', **kwargs): + """ + + Initialize logging for a function + + Parameters + ---------- + funcname: string + Name of function which will be recorded in log + level: int, or bool, or string + int or string : Logging level + bool: False - switch off logging + Text logging level for the message ("DEBUG", "INFO", + "WARNING", "ERROR", "CRITICAL") + mode: 'w' or 'a' + attach to logfile ('a' or start a new logfile ('w') + + """ + for hndl in logger.handlers: + if isinstance(hndl, logging.FileHandler): + has_file_handler = True + else: + has_file_handler = False + if level: + if not has_file_handler: + logname = 'headerlet.log' + fh = logging.FileHandler(logname, mode=mode) + fh.setFormatter(formatter) + fh.setLevel(logging.DEBUG) + logger.addHandler(fh) + logger.info("%s: Starting %s with arguments:\n\t %s" % + (time.asctime(), funcname, kwargs)) + + +def with_logging(func): + @functools.wraps(func) + def wrapped(*args, **kw): + level = kw.get('logging', 100) + mode = kw.get('logmode', 'w') + func_args = kw.copy() + if sys.version_info[0] >= 3: + argnames = func.__code__.co_varnames + else: + argnames = func.func_code.co_varnames + + for argname, arg in zip(argnames, args): + func_args[argname] = arg + + init_logging(func.__name__, level, mode, **func_args) + return func(*args, **kw) + return wrapped + +#### Utility functions +def is_par_blank(par): + return par in ['', ' ', 'INDEF', "None", None] + +def parse_filename(fname, mode='readonly'): + """ + Interprets the input as either a filename of a file that needs to be opened + or a PyFITS object. + + Parameters + ---------- + fname : str, `astropy.io.fits.HDUList` + Input pointing to a file or `astropy.io.fits.HDUList` object. + An input filename (str) will be expanded as necessary to + interpret any environmental variables + included in the filename. + + mode : string + Specifies what mode to use when opening the file, if it needs + to open the file at all [Default: 'readonly'] + + Returns + ------- + fobj : `astropy.io.fits.HDUList` + FITS file handle for input + + fname : str + Name of input file + + close_fobj : bool + Flag specifying whether or not fobj needs to be closed since it was + opened by this function. This allows a program to know whether they + need to worry about closing the FITS object as opposed to letting + the higher level interface close the object. + + """ + close_fobj = False + if not isinstance(fname, list): + if sys.version_info[0] >= 3: + is_string = isinstance(fname, str) + else: + is_string = isinstance(fname, basestring) + if is_string: + fname = fu.osfn(fname) + fobj = fits.open(fname, mode=mode) + close_fobj = True + else: + fobj = fname + if hasattr(fobj, 'filename'): + fname = fobj.filename() + else: + fname = '' + return fobj, fname, close_fobj + +def get_headerlet_kw_names(fobj, kw='HDRNAME'): + """ + Returns a list of specified keywords from all HeaderletHDU + extensions in a science file. + + Parameters + ---------- + fobj : str, `astropy.io.fits.HDUList` + kw : str + Name of keyword to be read and reported + """ + + fobj, fname, open_fobj = parse_filename(fobj) + + hdrnames = [] + for ext in fobj: + if isinstance(ext, fits.hdu.base.NonstandardExtHDU): + hdrnames.append(ext.header[kw]) + + if open_fobj: + fobj.close() + + return hdrnames + +def get_header_kw_vals(hdr, kwname, kwval, default=0): + if kwval is None: + if kwname in hdr: + kwval = hdr[kwname] + else: + kwval = default + return kwval + +@with_logging +def find_headerlet_HDUs(fobj, hdrext=None, hdrname=None, distname=None, + strict=True, logging=False, logmode='w'): + """ + Returns all HeaderletHDU extensions in a science file that matches + the inputs specified by the user. If no hdrext, hdrname or distname are + specified, this function will return a list of all HeaderletHDU objects. + + Parameters + ---------- + fobj : str, `astropy.io.fits.HDUList` + Name of FITS file or open fits object (`astropy.io.fits.HDUList` instance) + hdrext : int, tuple or None + index number(EXTVER) or extension tuple of HeaderletHDU to be returned + hdrname : string + value of HDRNAME for HeaderletHDU to be returned + distname : string + value of DISTNAME for HeaderletHDUs to be returned + strict : bool [Default: True] + Specifies whether or not at least one parameter needs to be provided + If False, all extension indices returned if hdrext, hdrname and distname + are all None. If True and hdrext, hdrname, and distname are all None, + raise an Exception requiring one to be specified. + logging : boolean + enable logging to a file called headerlet.log + logmode : 'w' or 'a' + log file open mode + + Returns + ------- + hdrlets : list + A list of all matching HeaderletHDU extension indices (could be just one) + + """ + + get_all = False + if hdrext is None and hdrname is None and distname is None: + if not strict: + get_all = True + else: + mess = """\n + ===================================================== + No valid Headerlet extension specified. + Either "hdrname", "hdrext", or "distname" needs to be specified. + ===================================================== + """ + logger.critical(mess) + raise ValueError + + fobj, fname, open_fobj = parse_filename(fobj) + + hdrlets = [] + if hdrext is not None and isinstance(hdrext, int): + if hdrext in range(len(fobj)): # insure specified hdrext is in fobj + if isinstance(fobj[hdrext], fits.hdu.base.NonstandardExtHDU) and \ + fobj[hdrext].header['EXTNAME'] == 'HDRLET': + hdrlets.append(hdrext) + else: + for ext in fobj: + if isinstance(ext, fits.hdu.base.NonstandardExtHDU): + if get_all: + hdrlets.append(fobj.index(ext)) + else: + if hdrext is not None: + if isinstance(hdrext, tuple): + hdrextname = hdrext[0] + hdrextnum = hdrext[1] + else: + hdrextname = 'HDRLET' + hdrextnum = hdrext + hdrext_match = ((hdrext is not None) and + (hdrextnum == ext.header['EXTVER']) and + (hdrextname == ext.header['EXTNAME'])) + hdrname_match = ((hdrname is not None) and + (hdrname == ext.header['HDRNAME'])) + distname_match = ((distname is not None) and + (distname == ext.header['DISTNAME'])) + if hdrext_match or hdrname_match or distname_match: + hdrlets.append(fobj.index(ext)) + + if open_fobj: + fobj.close() + + if len(hdrlets) == 0: + if hdrname: + kwerr = 'hdrname' + kwval = hdrname + elif hdrext: + kwerr = 'hdrext' + kwval = hdrext + else: + kwerr = 'distname' + kwval = distname + message = """\n + ===================================================== + No valid Headerlet extension found!' + "%s" = %s not found in %s.' % (kwerr, kwval, fname) + ===================================================== + """ + logger.critical(message) + raise ValueError + + return hdrlets + +def verify_hdrname_is_unique(fobj, hdrname): + """ + Verifies that no other HeaderletHDU extension has the specified hdrname. + + Parameters + ---------- + fobj : str, `astropy.io.fits.HDUList` + Name of FITS file or open fits file object + hdrname : str + value of HDRNAME for HeaderletHDU to be compared as unique + + Returns + ------- + unique: bool + If True, no other HeaderletHDU has the specified HDRNAME value + """ + hdrnames_list = get_headerlet_kw_names(fobj) + unique = not(hdrname in hdrnames_list) + + return unique + +def update_versions(sourcehdr, desthdr): + """ + Update keywords which store version numbers + """ + phdukw = {'PYWCSVER': 'Version of PYWCS used to updated the WCS', + 'UPWCSVER': 'Version of STWCS used to updated the WCS' + } + for key in phdukw: + try: + desthdr[key] = (sourcehdr[key], sourcehdr.comments[key]) + except KeyError: + desthdr[key] = (" ", phdukw[key]) + +def update_ref_files(source, dest): + """ + Update the reference files name in the primary header of 'dest' + using values from 'source' + + Parameters + ---------- + source : `astropy.io.fits.Header` + dest : `astropy.io.fits.Header` + """ + logger.info("Updating reference files") + phdukw = {'NPOLFILE': True, + 'IDCTAB': True, + 'D2IMFILE': True, + 'SIPNAME': True, + 'DISTNAME': True} + + for key in phdukw: + try: + try: + del dest[key] + except: + pass + dest.set(key, source[key], source.comments[key]) + except KeyError: + phdukw[key] = False + return phdukw + +def print_summary(summary_cols, summary_dict, pad=2, maxwidth=None, idcol=None, + output=None, clobber=True, quiet=False ): + """ + Print out summary dictionary to STDOUT, and possibly an output file + + """ + nrows = None + if idcol: + nrows = len(idcol['vals']) + + # Find max width of each column + column_widths = {} + for kw in summary_dict: + colwidth = np.array(summary_dict[kw]['width']).max() + if maxwidth: + colwidth = min(colwidth, maxwidth) + column_widths[kw] = colwidth + pad + if nrows is None: + nrows = len(summary_dict[kw]['vals']) + + # print rows now + outstr = '' + # Start with column names + if idcol: + outstr += COLUMN_FMT.format(idcol['name'], width=idcol['width'] + pad) + for kw in summary_cols: + outstr += COLUMN_FMT.format(kw, width=column_widths[kw]) + outstr += '\n' + # Now, add a row for each headerlet + for row in range(nrows): + if idcol: + outstr += COLUMN_FMT.format(idcol['vals'][row], + width=idcol['width']+pad) + for kw in summary_cols: + val = summary_dict[kw]['vals'][row][:(column_widths[kw]-pad)] + outstr += COLUMN_FMT.format(val, width=column_widths[kw]) + outstr += '\n' + if not quiet: + print(outstr) + + # If specified, write info to separate text file + write_file = False + if output: + output = fu.osfn(output) # Expand any environment variables in filename + write_file = True + if os.path.exists(output): + if clobber: + os.remove(output) + else: + print('WARNING: Not writing results to file!') + print(' Output text file ', output, ' already exists.') + print(' Set "clobber" to True or move file before trying again.') + write_file = False + if write_file: + fout = open(output, mode='w') + fout.write(outstr) + fout.close() + +#### Private utility functions +def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, + sipname, npolfile, d2imfile, + nmatch,catalog, wcskey, + author, descrip, history): + # convert input values into valid FITS kw values + if author is None: + author = '' + if descrip is None: + descrip = '' + + sipname, idctab = utils.build_sipname(fobj, fname, sipname) + logger.info("Setting sipname value to %s" % sipname) + + npolname, npolfile = utils.build_npolname(fobj, npolfile) + logger.info("Setting npolfile value to %s" % npolname) + + d2imname, d2imfile = utils.build_d2imname(fobj,d2imfile) + logger.info("Setting d2imfile value to %s" % d2imname) + + distname = utils.build_distname(sipname, npolname, d2imname) + logger.info("Setting distname to %s" % distname) + + # open file and parse comments + if history not in ['', ' ', None, 'INDEF'] and os.path.isfile(history): + f = open(fu.osfn(history)) + history = f.readlines() + f.close() + else: + history = '' + + rms_ra = fobj[wcsext].header.get("CRDER1"+wcskey, 0) + rms_dec = fobj[wcsext].header.get("CRDER2"+wcskey, 0) + if not nmatch: + nmatch = fobj[wcsext].header.get("NMATCH"+wcskey, 0) + if not catalog: + catalog = fobj[wcsext].header.get('CATALOG'+wcskey, "") + # get the version of STWCS used to create the WCS of the science file. + #try: + #upwcsver = fobj[0].header.cards[fobj[0].header.index('UPWCSVER')] + #except KeyError: + #upwcsver = pyfits.Card("UPWCSVER", " ", + #"Version of STWCS used to update the WCS") + #try: + #pywcsver = fobj[0].header.cards[fobj[0].header.index('PYWCSVER')] + #except KeyError: + #pywcsver = pyfits.Card("PYWCSVER", " ", + #"Version of PYWCS used to update the WCS") + upwcsver = fobj[0].header.get('UPWCSVER', "") + pywcsver = fobj[0].header.get('PYWCSVER', "") + # build Primary HDU + phdu = fits.PrimaryHDU() + phdu.header['DESTIM'] = (destim, 'Destination observation root name') + phdu.header['HDRNAME'] = (hdrname, 'Headerlet name') + fmt = "%Y-%m-%dT%H:%M:%S" + phdu.header['DATE'] = (time.strftime(fmt), 'Date FITS file was generated') + phdu.header['WCSNAME'] = (wcsname, 'WCS name') + phdu.header['DISTNAME'] = (distname, 'Distortion model name') + phdu.header['SIPNAME'] = (sipname, + 'origin of SIP polynomial distortion model') + phdu.header['NPOLFILE'] = (npolfile, + 'origin of non-polynmial distortion model') + phdu.header['D2IMFILE'] = (d2imfile, + 'origin of detector to image correction') + phdu.header['IDCTAB'] = (idctab, + 'origin of Polynomial Distortion') + phdu.header['AUTHOR'] = (author, 'headerlet created by this user') + phdu.header['DESCRIP'] = (descrip, + 'Short description of headerlet solution') + phdu.header['RMS_RA'] = (rms_ra, + 'RMS in RA at ref pix of headerlet solution') + phdu.header['RMS_DEC'] = (rms_dec, + 'RMS in Dec at ref pix of headerlet solution') + phdu.header['NMATCH'] = (nmatch, + 'Number of sources used for headerlet solution') + phdu.header['CATALOG'] = (catalog, + 'Astrometric catalog used for headerlet ' + 'solution') + phdu.header['UPWCSVER'] = (upwcsver, "Version of STWCS used to update the WCS") + phdu.header['PYWCSVER'] = (pywcsver, "Version of PYWCS used to update the WCS") + + # clean up history string in order to remove whitespace characters that + # would cause problems with FITS + if isinstance(history, list): + history_str = '' + for line in history: + history_str += line + else: + history_str = history + history_lines = textwrap.wrap(history_str, width=70) + for hline in history_lines: + phdu.header.add_history(hline) + + return phdu + + +#### Public Interface functions +@with_logging +def extract_headerlet(filename, output, extnum=None, hdrname=None, + clobber=False, logging=True): + """ + Finds a headerlet extension in a science file and writes it out as + a headerlet FITS file. + + If both hdrname and extnum are given they should match, if not + raise an Exception + + Parameters + ---------- + filename: string or HDUList or Python list + This specifies the name(s) of science file(s) from which headerlets + will be extracted. + + String input formats supported include use of wild-cards, IRAF-style + '@'-files (given as '@<filename>') and comma-separated list of names. + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + If a list of filenames has been specified, it will extract a + headerlet from the same extnum from all filenames. + output: string + Filename or just rootname of output headerlet FITS file + If string does not contain '.fits', it will create a filename with + '_hlet.fits' suffix + extnum: int + Extension number which contains the headerlet to be written out + hdrname: string + Unique name for headerlet, stored as the HDRNAME keyword + It stops if a value is not provided and no extnum has been specified + clobber: bool + If output file already exists, this parameter specifies whether or not + to overwrite that file [Default: False] + logging: boolean + enable logging to a file + + """ + + if isinstance(filename, fits.HDUList): + filename = [filename] + else: + filename, oname = parseinput.parseinput(filename) + + for f in filename: + fobj, fname, close_fobj = parse_filename(f) + frootname = fu.buildNewRootname(fname) + if hdrname in ['', ' ', None, 'INDEF'] and extnum is None: + if close_fobj: + fobj.close() + logger.critical("Expected a valid extnum or hdrname parameter") + raise ValueError + if hdrname is not None: + extn_from_hdrname = find_headerlet_HDUs(fobj, hdrname=hdrname)[0] + if extn_from_hdrname != extnum: + logger.critical("hdrname and extnmu should refer to the same FITS extension") + raise ValueError + else: + hdrhdu = fobj[extn_from_hdrname] + else: + hdrhdu = fobj[extnum] + + if not isinstance(hdrhdu, HeaderletHDU): + logger.critical("Specified extension is not a headerlet") + raise ValueError + + hdrlet = hdrhdu.headerlet + + if output is None: + output = frootname + + if '.fits' in output: + outname = output + else: + outname = '%s_hlet.fits' % output + + hdrlet.tofile(outname, clobber=clobber) + + if close_fobj: + fobj.close() + + +@with_logging +def write_headerlet(filename, hdrname, output=None, sciext='SCI', + wcsname=None, wcskey=None, destim=None, + sipname=None, npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + nmatch=None, catalog=None, + attach=True, clobber=False, logging=False): + + """ + Save a WCS as a headerlet FITS file. + + This function will create a headerlet, write out the headerlet to a + separate headerlet file, then, optionally, attach it as an extension + to the science image (if it has not already been archived) + + Either wcsname or wcskey must be provided; if both are given, they must + match a valid WCS. + + Updates wcscorr if necessary. + + Parameters + ---------- + filename: string or HDUList or Python list + This specifies the name(s) of science file(s) from which headerlets + will be created and written out. + String input formats supported include use of wild-cards, IRAF-style + '@'-files (given as '@<filename>') and comma-separated list of names. + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + hdrname: string + Unique name for this headerlet, stored as HDRNAME keyword + output: string or None + Filename or just rootname of output headerlet FITS file + If string does not contain '.fits', it will create a filename + starting with the science filename and ending with '_hlet.fits'. + If None, a default filename based on the input filename will be + generated for the headerlet FITS filename + sciext: string + name (EXTNAME) of extension that contains WCS to be saved + wcsname: string + name of WCS to be archived, if " ": stop + wcskey: one of A...Z or " " or "PRIMARY" + if " " or "PRIMARY" - archive the primary WCS + destim: string + DESTIM keyword + if NOne, use ROOTNAME or science file name + sipname: string or None (default) + Name of unique file where the polynomial distortion coefficients were + read from. If None, the behavior is: + The code looks for a keyword 'SIPNAME' in the science header + If not found, for HST it defaults to 'IDCTAB' + If there is no SIP model the value is 'NOMODEL' + If there is a SIP model but no SIPNAME, it is set to 'UNKNOWN' + npolfile: string or None (default) + Name of a unique file where the non-polynomial distortion was stored. + If None: + The code looks for 'NPOLFILE' in science header. + If 'NPOLFILE' was not found and there is no npol model, it is set to 'NOMODEL' + If npol model exists, it is set to 'UNKNOWN' + d2imfile: string + Name of a unique file where the detector to image correction was + stored. If None: + The code looks for 'D2IMFILE' in the science header. + If 'D2IMFILE' is not found and there is no d2im correction, + it is set to 'NOMODEL' + If d2im correction exists, but 'D2IMFILE' is missing from science + header, it is set to 'UNKNOWN' + author: string + Name of user who created the headerlet, added as 'AUTHOR' keyword + to headerlet PRIMARY header + descrip: string + Short description of the solution provided by the headerlet + This description will be added as the single 'DESCRIP' keyword + to the headerlet PRIMARY header + history: filename, string or list of strings + Long (possibly multi-line) description of the solution provided + by the headerlet. These comments will be added as 'HISTORY' cards + to the headerlet PRIMARY header + If filename is specified, it will format and attach all text from + that file as the history. + attach: bool + Specify whether or not to attach this headerlet as a new extension + It will verify that no other headerlet extension has been created with + the same 'hdrname' value. + clobber: bool + If output file already exists, this parameter specifies whether or not + to overwrite that file [Default: False] + logging: boolean + enable file logging + """ + + if isinstance(filename, fits.HDUList): + filename = [filename] + else: + filename, oname = parseinput.parseinput(filename) + + for f in filename: + if isinstance(f, str): + fname = f + else: + fname = f.filename() + + if wcsname in [None, ' ', '', 'INDEF'] and wcskey is None: + message = """\n + No valid WCS found found in %s. + A valid value for either "wcsname" or "wcskey" + needs to be specified. + """ % fname + logger.critical(message) + raise ValueError + + # Translate 'wcskey' value for PRIMARY WCS to valid altwcs value of ' ' + if wcskey == 'PRIMARY': + wcskey = ' ' + + if attach: + umode = 'update' + else: + umode = 'readonly' + + fobj, fname, close_fobj = parse_filename(f, mode=umode) + + # Interpret sciext input for this file + if isinstance(sciext, int): + sciextlist = [sciext] # allow for specification of simple FITS header + elif isinstance(sciext, str): + numsciext = countExtn(fobj, sciext) + if numsciext > 0: + sciextlist = [tuple((sciext,i)) for i in range(1, numsciext+1)] + else: + sciextlist = [0] + elif isinstance(sciext, list): + sciextlist = sciext + else: + errstr = "Expected sciext to be a list of FITS extensions with science data\n"+\ + " a valid EXTNAME string, or an integer." + logger.critical(errstr) + raise ValueError + + wnames = altwcs.wcsnames(fobj,ext=sciextlist[0]) + + # Insure that WCSCORR table has been created with all original + # WCS's recorded prior to adding the headerlet WCS + wcscorr.init_wcscorr(fobj) + + if wcsname is None: + scihdr = fobj[sciextlist[0]].header + wname = scihdr['wcsname'+wcskey] + else: + wname = wcsname + if hdrname in [None, ' ', '']: + hdrname = wcsname + + logger.critical('Creating the headerlet from image %s' % fname) + hdrletobj = create_headerlet(fobj, sciext=sciextlist, + wcsname=wname, wcskey=wcskey, + hdrname=hdrname, + sipname=sipname, npolfile=npolfile, + d2imfile=d2imfile, author=author, + descrip=descrip, history=history, + nmatch=nmatch, catalog=catalog, + logging=False) + + if attach: + # Check to see whether or not a HeaderletHDU with + #this hdrname already exists + hdrnames = get_headerlet_kw_names(fobj) + if hdrname not in hdrnames: + hdrlet_hdu = HeaderletHDU.fromheaderlet(hdrletobj) + + if destim is not None: + hdrlet_hdu.header['destim'] = destim + + fobj.append(hdrlet_hdu) + + # Update the WCSCORR table with new rows from the headerlet's WCSs + wcscorr.update_wcscorr(fobj, source=hdrletobj, + extname='SIPWCS', wcs_id=wname) + + utils.updateNEXTENDKw(fobj) + fobj.flush() + else: + message = """ + Headerlet with hdrname %s already archived for WCS %s. + No new headerlet appended to %s. + """ % (hdrname, wname, fname) + logger.critical(message) + + if close_fobj: + logger.info('Closing image in write_headerlet()...') + fobj.close() + + frootname = fu.buildNewRootname(fname) + + if output is None: + # Generate default filename for headerlet FITS file + outname = '{0}_hlet.fits'.format(frootname) + else: + outname = output + + if not outname.endswith('.fits'): + outname = '{0}_{1}_hlet.fits'.format(frootname,outname) + + # If user specifies an output filename for headerlet, write it out + hdrletobj.tofile(outname, clobber=clobber) + logger.critical( 'Created Headerlet file %s ' % outname) + + del hdrletobj + +@with_logging +def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, + wcskey=" ", wcsname=None, + sipname=None, npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + nmatch=None, catalog=None, + logging=False, logmode='w'): + """ + Create a headerlet from a WCS in a science file + If both wcskey and wcsname are given they should match, if not + raise an Exception + + Parameters + ---------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + sciext: string or python list (default: 'SCI') + Extension in which the science data with the linear WCS is. + The headerlet will be created from these extensions. + If string - a valid EXTNAME is expected + If int - specifies an extension with a valid WCS, such as 0 for a + simple FITS file + If list - a list of FITS extension numbers or strings representing + extension tuples, e.g. ('SCI, 1') is expected. + hdrname: string + value of HDRNAME keyword + Takes the value from the HDRNAME<wcskey> keyword, if not available from WCSNAME<wcskey> + It stops if neither is found in the science file and a value is not provided + destim: string or None + name of file this headerlet can be applied to + if None, use ROOTNAME keyword + wcskey: char (A...Z) or " " or "PRIMARY" or None + a char representing an alternate WCS to be used for the headerlet + if " ", use the primary (default) + if None use wcsname + wcsname: string or None + if wcskey is None use wcsname specified here to choose an alternate WCS for the headerlet + sipname: string or None (default) + Name of unique file where the polynomial distortion coefficients were + read from. If None, the behavior is: + The code looks for a keyword 'SIPNAME' in the science header + If not found, for HST it defaults to 'IDCTAB' + If there is no SIP model the value is 'NOMODEL' + If there is a SIP model but no SIPNAME, it is set to 'UNKNOWN' + npolfile: string or None (default) + Name of a unique file where the non-polynomial distortion was stored. + If None: + The code looks for 'NPOLFILE' in science header. + If 'NPOLFILE' was not found and there is no npol model, it is set to 'NOMODEL' + If npol model exists, it is set to 'UNKNOWN' + d2imfile: string + Name of a unique file where the detector to image correction was + If None: + The code looks for 'D2IMFILE' in the science header. + If 'D2IMFILE' is not found and there is no d2im correction, + it is set to 'NOMODEL' + If d2im correction exists, but 'D2IMFILE' is missing from science + header, it is set to 'UNKNOWN' + author: string + Name of user who created the headerlet, added as 'AUTHOR' keyword + to headerlet PRIMARY header + descrip: string + Short description of the solution provided by the headerlet + This description will be added as the single 'DESCRIP' keyword + to the headerlet PRIMARY header + history: filename, string or list of strings + Long (possibly multi-line) description of the solution provided + by the headerlet. These comments will be added as 'HISTORY' cards + to the headerlet PRIMARY header + If filename is specified, it will format and attach all text from + that file as the history. + nmatch: int (optional) + Number of sources used in the new solution fit + catalog: string (optional) + Astrometric catalog used for headerlet solution + logging: boolean + enable file logging + logmode: 'w' or 'a' + log file open mode + + Returns + ------- + Headerlet object + + """ + if wcskey == 'O': + message = "Warning: 'O' is a reserved key for the original WCS. Quitting..." + logger.info(message) + return None + + fobj, fname, close_file = parse_filename(filename) + # based on `sciext` create a list of (EXTNAME, EXTVER) tuples + # of extensions with WCS to be saved in a headerlet + sciext = get_extname_extver_list(fobj, sciext) + logger.debug("Data extensions from which to create headerlet:\n\t %s" + % (str(sciext))) + if not sciext: + logger.critical("No valid target extensions found in file %s" % fname) + raise ValueError + + # Define extension to evaluate for verification of input parameters + wcsext = sciext[0] + logger.debug("sciext in create_headerlet is %s" % str(sciext)) + # Translate 'wcskey' value for PRIMARY WCS to valid altwcs value of ' ' + if wcskey == 'PRIMARY': + wcskey = ' ' + logger.info("wcskey reset from 'PRIMARY' to ' '") + wcskey = wcskey.upper() + wcsnamekw = "".join(["WCSNAME", wcskey.upper()]).rstrip() + hdrnamekw = "".join(["HDRNAME", wcskey.upper()]).rstrip() + + wnames = altwcs.wcsnames(fobj, ext=wcsext) + if not wcsname: + # User did not specify a value for 'wcsname' + if wcsnamekw in fobj[wcsext].header: + #check if there's a WCSNAME for this wcskey in the header + wcsname = fobj[wcsext].header[wcsnamekw] + logger.info("Setting wcsname from header[%s] to %s" % (wcsnamekw, wcsname)) + else: + if hdrname not in ['', ' ', None, "INDEF"]: + """ + If wcsname for this wcskey was not provided + and WCSNAME<wcskey> does not exist in the header + and hdrname is provided, then + use hdrname as WCSNAME for the headerlet. + """ + wcsname = hdrname + logger.debug("Setting wcsname from hdrname to %s" % hdrname) + else: + if hdrnamekw in fobj[wcsext].header: + wcsname = fobj[wcsext].header[hdrnamekw] + logger.debug("Setting wcsname from header[%s] to %s" % (hdrnamekw, wcsname)) + else: + message = """ + Required keywords 'HDRNAME' or 'WCSNAME' not found! + Please specify a value for parameter 'hdrname' + or update header with 'WCSNAME' keyword. + """ + logger.critical(message) + raise KeyError + else: + # Verify that 'wcsname' and 'wcskey' values specified by user reference + # the same WCS + wname = fobj[wcsext].header[wcsnamekw] + if wcsname != wname: + message = "\tInconsistent values for 'wcskey' and 'wcsname' specified!\n" + message += " 'wcskey' = %s and 'wcsname' = %s. \n" % (wcskey, wcsname) + message += "Actual value of %s found to be %s. \n" % (wcsnamekw, wname) + logger.critical(message) + raise KeyError + wkeys = altwcs.wcskeys(fobj, ext=wcsext) + if wcskey != ' ': + if wcskey not in wkeys: + logger.critical('No WCS with wcskey=%s found in extension %s. Skipping...' % (wcskey, str(wcsext))) + raise ValueError("No WCS with wcskey=%s found in extension %s. Skipping...' % (wcskey, str(wcsext))") + + # get remaining required keywords + if destim is None: + if 'ROOTNAME' in fobj[0].header: + destim = fobj[0].header['ROOTNAME'] + logger.info("Setting destim to rootname of the science file") + else: + destim = fname + logger.info('DESTIM not provided') + logger.info('Keyword "ROOTNAME" not found') + logger.info('Using file name as DESTIM') + + if not hdrname: + # check if HDRNAME<wcskey> is in header + if hdrnamekw in fobj[wcsext].header: + hdrname = fobj[wcsext].header[hdrnamekw] + else: + if wcsnamekw in fobj[wcsext].header: + hdrname = fobj[wcsext].header[wcsnamekw] + message = """ + Using default value for HDRNAME of "%s" derived from %s. + """ % (hdrname, wcsnamekw) + logger.info(message) + logger.info("Setting hdrname to %s from header[%s]" + % (hdrname, wcsnamekw)) + else: + message = "Required keywords 'HDRNAME' or 'WCSNAME' not found" + logger.critical(message) + raise KeyError + + + + hdul = [] + phdu = _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, + sipname, npolfile, d2imfile, + nmatch, catalog, wcskey, + author, descrip, history) + hdul.append(phdu) + wcsdvarr_extns = [] + """ + nd2i is a counter for d2i extensions to be used when the science file + has an old d2i correction format. The old format did not write EXTVER + kw for the d2i correction in the science header bu tthe new format expects + them. + """ + nd2i_extver = 1 + for ext in sciext: + wkeys = altwcs.wcskeys(fobj, ext=ext) + if wcskey != ' ': + if wcskey not in wkeys: + logger.debug( + 'No WCS with wcskey=%s found in extension %s. ' + 'Skipping...' % (wcskey, str(ext))) + raise ValueError("") + + hwcs = HSTWCS(fobj, ext=ext, wcskey=wcskey) + + whdul = hwcs.to_fits(relax=True, key=" ") + if hasattr(hwcs, 'orientat'): + orient_comment = "positions angle of image y axis (deg. e of n)" + whdul[0].header['ORIENTAT'] = (hwcs.orientat, orient_comment) + + whdul[0].header.append(('TG_ENAME', ext[0], 'Target science data extname')) + whdul[0].header.append(('TG_EVER', ext[1], 'Target science data extver')) + + if hwcs.wcs.has_cd(): + whdul[0].header = altwcs.pc2cd(whdul[0].header) + + idckw = hwcs._idc2hdr() + whdul[0].header.extend(idckw) + + if hwcs.det2im1 or hwcs.det2im2: + try: + whdul[0].header.append(fobj[ext].header.cards['D2IMEXT']) + except KeyError: + pass + whdul[0].header.extend(fobj[ext].header.cards['D2IMERR*']) + if 'D2IM1.EXTVER' in whdul[0].header: + try: + whdul[0].header['D2IM1.EXTVER'] = fobj[ext].header['D2IM1.EXTVER'] + except KeyError: + whdul[0].header['D2IM1.EXTVER'] = nd2i_extver + nd2i_extver += 1 + if 'D2IM2.EXTVER' in whdul[0].header: + try: + whdul[0].header['D2IM2.EXTVER'] = fobj[ext].header['D2IM2.EXTVER'] + except KeyError: + whdul[0].header['D2IM2.EXTVER'] = nd2i_extver + nd2i_extver += 1 + + if hwcs.cpdis1 or hwcs.cpdis2: + whdul[0].header.extend(fobj[ext].header.cards['CPERR*']) + try: + whdul[0].header.append(fobj[ext].header.cards['NPOLEXT']) + except KeyError: + pass + if 'DP1.EXTVER' in whdul[0].header: + whdul[0].header['DP1.EXTVER'] = fobj[ext].header['DP1.EXTVER'] + if 'DP2.EXTVER' in whdul[0].header: + whdul[0].header['DP2.EXTVER'] = fobj[ext].header['DP2.EXTVER'] + ihdu = fits.ImageHDU(header=whdul[0].header, name='SIPWCS') + + if ext[0] != "PRIMARY": + ihdu.update_ext_version(fobj[ext].header['EXTVER'], comment='Extension version') + + hdul.append(ihdu) + + if hwcs.cpdis1: + whdu = whdul[('WCSDVARR', 1)].copy() + whdu.update_ext_version(fobj[ext].header['DP1.EXTVER']) + hdul.append(whdu) + if hwcs.cpdis2: + whdu = whdul[('WCSDVARR', 2)].copy() + whdu.update_ext_version(fobj[ext].header['DP2.EXTVER']) + hdul.append(whdu) + + if hwcs.det2im1: + whdu = whdul[('D2IMARR', 1)].copy() + whdu.update_ext_version(ihdu.header['D2IM1.EXTVER']) + hdul.append(whdu) + if hwcs.det2im2: + whdu = whdul[('D2IMARR', 2)].copy() + whdu.update_ext_version(ihdu.header['D2IM2.EXTVER']) + hdul.append(whdu) + + + #if hwcs.det2im1 or hwcs.det2im2: + #try: + #darr = hdul[('D2IMARR', 1)] + #except KeyError: + #whdu = whdul[('D2IMARR')] + #whdu.update_ext_version(1) + #hdul.append(whdu) + if close_file: + fobj.close() + + hlet = Headerlet(hdul, logging=logging, logmode='a') + hlet.init_attrs() + return hlet + +@with_logging +def apply_headerlet_as_primary(filename, hdrlet, attach=True, archive=True, + force=False, logging=False, logmode='a'): + """ + Apply headerlet 'hdrfile' to a science observation 'destfile' as the primary WCS + + Parameters + ---------- + filename: string or list of strings + File name(s) of science observation whose WCS solution will be updated + hdrlet: string or list of strings + Headerlet file(s), must match 1-to-1 with input filename(s) + attach: boolean + True (default): append headerlet to FITS file as a new extension. + archive: boolean + True (default): before updating, create a headerlet with the + WCS old solution. + force: boolean + If True, this will cause the headerlet to replace the current PRIMARY + WCS even if it has a different distortion model. [Default: False] + logging: boolean + enable file logging + logmode: 'w' or 'a' + log file open mode + """ + if not isinstance(filename, list): + filename = [filename] + if not isinstance(hdrlet, list): + hdrlet = [hdrlet] + if len(hdrlet) != len(filename): + logger.critical("Filenames must have matching headerlets. " + "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + + for fname,h in zip(filename,hdrlet): + print("Applying {0} as Primary WCS to {1}".format(h,fname)) + hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) + hlet.apply_as_primary(fname, attach=attach, archive=archive, + force=force) + + +@with_logging +def apply_headerlet_as_alternate(filename, hdrlet, attach=True, wcskey=None, + wcsname=None, logging=False, logmode='w'): + """ + Apply headerlet to a science observation as an alternate WCS + + Parameters + ---------- + filename: string or list of strings + File name(s) of science observation whose WCS solution will be updated + hdrlet: string or list of strings + Headerlet file(s), must match 1-to-1 with input filename(s) + attach: boolean + flag indicating if the headerlet should be attached as a + HeaderletHDU to fobj. If True checks that HDRNAME is unique + in the fobj and stops if not. + wcskey: string + Key value (A-Z, except O) for this alternate WCS + If None, the next available key will be used + wcsname: string + Name to be assigned to this alternate WCS + WCSNAME is a required keyword in a Headerlet but this allows the + user to change it as desired. + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + if not isinstance(filename, list): + filename = [filename] + if not isinstance(hdrlet, list): + hdrlet = [hdrlet] + if len(hdrlet) != len(filename): + logger.critical("Filenames must have matching headerlets. " + "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + + for fname,h in zip(filename,hdrlet): + print('Applying {0} as an alternate WCS to {1}'.format(h,fname)) + hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) + hlet.apply_as_alternate(fname, attach=attach, + wcsname=wcsname, wcskey=wcskey) + + +@with_logging +def attach_headerlet(filename, hdrlet, logging=False, logmode='a'): + """ + Attach Headerlet as an HeaderletHDU to a science file + + Parameters + ---------- + filename: HDUList or list of HDULists + science file(s) to which the headerlet should be applied + hdrlet: string, Headerlet object or list of strings or Headerlet objects + string representing a headerlet file(s), must match 1-to-1 input filename(s) + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + if not isinstance(filename, list): + filename = [filename] + if not isinstance(hdrlet, list): + hdrlet = [hdrlet] + if len(hdrlet) != len(filename): + logger.critical("Filenames must have matching headerlets. " + "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + + for fname,h in zip(filename,hdrlet): + print('Attaching {0} as Headerlet extension to {1}'.format(h,fname)) + hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) + hlet.attach_to_file(fname,archive=True) + + +@with_logging +def delete_headerlet(filename, hdrname=None, hdrext=None, distname=None, + logging=False, logmode='w'): + """ + Deletes HeaderletHDU(s) with same HDRNAME from science files + + Notes + ----- + One of hdrname, hdrext or distname should be given. + If hdrname is given - delete a HeaderletHDU with a name HDRNAME from fobj. + If hdrext is given - delete HeaderletHDU in extension. + If distname is given - deletes all HeaderletHDUs with a specific distortion model from fobj. + Updates wcscorr + + Parameters + ---------- + filename: string, HDUList or list of strings + Filename can be specified as a single filename or HDUList, or + a list of filenames + Each input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + hdrname: string or None + HeaderletHDU primary header keyword HDRNAME + hdrext: int, tuple or None + HeaderletHDU FITS extension number + tuple has the form ('HDRLET', 1) + distname: string or None + distortion model as specified in the DISTNAME keyword + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + if not isinstance(filename, list): + filename = [filename] + + for f in filename: + print("Deleting Headerlet from ",f) + _delete_single_headerlet(f, hdrname=hdrname, hdrext=hdrext, + distname=distname, logging=logging, logmode='a') + +def _delete_single_headerlet(filename, hdrname=None, hdrext=None, distname=None, + logging=False, logmode='w'): + """ + Deletes HeaderletHDU(s) from a SINGLE science file + + Notes + ----- + One of hdrname, hdrext or distname should be given. + If hdrname is given - delete a HeaderletHDU with a name HDRNAME from fobj. + If hdrext is given - delete HeaderletHDU in extension. + If distname is given - deletes all HeaderletHDUs with a specific distortion model from fobj. + Updates wcscorr + + Parameters + ---------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + hdrname: string or None + HeaderletHDU primary header keyword HDRNAME + hdrext: int, tuple or None + HeaderletHDU FITS extension number + tuple has the form ('HDRLET', 1) + distname: string or None + distortion model as specified in the DISTNAME keyword + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + hdrlet_ind = find_headerlet_HDUs(filename, hdrname=hdrname, hdrext=hdrext, + distname=distname, logging=logging, logmode='a') + if len(hdrlet_ind) == 0: + message = """ + No HDUs deleted... No Headerlet HDUs found with ' + hdrname = %s + hdrext = %s + distname = %s + Please review input parameters and try again. + """ % (hdrname, str(hdrext), distname) + logger.critical(message) + return + + fobj, fname, close_fobj = parse_filename(filename, mode='update') + + # delete row(s) from WCSCORR table now... + # + # + if hdrname not in ['', ' ', None, 'INDEF']: + selections = {'hdrname': hdrname} + elif hdrname in ['', ' ', None, 'INDEF'] and hdrext is not None: + selections = {'hdrname': fobj[hdrext].header['hdrname']} + else: + selections = {'distname': distname} + wcscorr.delete_wcscorr_row(fobj['WCSCORR'].data, selections) + + # delete the headerlet extension now + for hdrind in hdrlet_ind: + del fobj[hdrind] + + utils.updateNEXTENDKw(fobj) + # Update file object with changes + fobj.flush() + # close file, if was opened by this function + if close_fobj: + fobj.close() + logger.critical('Deleted headerlet from extension(s) %s ' % str(hdrlet_ind)) + + +def headerlet_summary(filename, columns=None, pad=2, maxwidth=None, + output=None, clobber=True, quiet=False): + """ + + Print a summary of all HeaderletHDUs in a science file to STDOUT, and + optionally to a text file + The summary includes: + HDRLET_ext_number HDRNAME WCSNAME DISTNAME SIPNAME NPOLFILE D2IMFILE + + Parameters + ---------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + columns: list + List of headerlet PRIMARY header keywords to report in summary + By default (set to None), it will use the default set of keywords + defined as the global list DEFAULT_SUMMARY_COLS + pad: int + Number of padding spaces to put between printed columns + [Default: 2] + maxwidth: int + Maximum column width(not counting padding) for any column in summary + By default (set to None), each column's full width will be used + output: string (optional) + Name of optional output file to record summary. This filename + can contain environment variables. + [Default: None] + clobber: bool + If True, will overwrite any previous output file of same name + quiet: bool + If True, will NOT report info to STDOUT + + """ + if columns is None: + summary_cols = DEFAULT_SUMMARY_COLS + else: + summary_cols = columns + + summary_dict = {} + for kw in summary_cols: + summary_dict[kw] = copy.deepcopy(COLUMN_DICT) + + # Define Extension number column + extnums_col = copy.deepcopy(COLUMN_DICT) + extnums_col['name'] = 'EXTN' + extnums_col['width'] = 6 + + fobj, fname, close_fobj = parse_filename(filename) + # find all HDRLET extensions and combine info into a single summary + for extn in fobj: + if 'extname' in extn.header and extn.header['extname'] == 'HDRLET': + hdrlet_indx = fobj.index_of(('hdrlet', extn.header['extver'])) + try: + ext_cols, ext_summary = extn.headerlet.summary(columns=summary_cols) + extnums_col['vals'].append(hdrlet_indx) + for kw in summary_cols: + for key in COLUMN_DICT: + summary_dict[kw][key].extend(ext_summary[kw][key]) + except: + print("Skipping headerlet") + print("Could not read Headerlet from extension ", hdrlet_indx) + + if close_fobj: + fobj.close() + + # Print out the summary dictionary + print_summary(summary_cols, summary_dict, pad=pad, maxwidth=maxwidth, + idcol=extnums_col, output=output, + clobber=clobber, quiet=quiet) + + +@with_logging +def restore_from_headerlet(filename, hdrname=None, hdrext=None, archive=True, + force=False, logging=False, logmode='w'): + """ + Restores a headerlet as a primary WCS + + Parameters + ---------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + hdrname: string + HDRNAME keyword of HeaderletHDU + hdrext: int or tuple + Headerlet extension number of tuple ('HDRLET',2) + archive: boolean (default: True) + When the distortion model in the headerlet is the same as the distortion model of + the science file, this flag indicates if the primary WCS should be saved as an alternate + nd a headerlet extension. + When the distortion models do not match this flag indicates if the current primary and + alternate WCSs should be archived as headerlet extensions and alternate WCS. + force: boolean (default:False) + When the distortion models of the headerlet and the primary do not match, and archive + is False, this flag forces an update of the primary. + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + + hdrlet_ind = find_headerlet_HDUs(filename, hdrext=hdrext, hdrname=hdrname) + + fobj, fname, close_fobj = parse_filename(filename, mode='update') + + if len(hdrlet_ind) > 1: + if hdrext: + kwerr = 'hdrext' + kwval = hdrext + else: + kwerr = 'hdrname' + kwval = hdrname + message = """ + Multiple Headerlet extensions found with the same name. + %d Headerlets with "%s" = %s found in %s. + """% (len(hdrlet_ind), kwerr, kwval, fname) + if close_fobj: + fobj.close() + logger.critical(message) + raise ValueError + + hdrlet_indx = hdrlet_ind[0] + + # read headerlet from HeaderletHDU into memory + if hasattr(fobj[hdrlet_ind[0]], 'hdulist'): + hdrlet = fobj[hdrlet_indx].hdulist + else: + hdrlet = fobj[hdrlet_indx].headerlet # older convention in PyFITS + + # read in the names of the extensions which HeaderletHDU updates + extlist = [] + for ext in hdrlet: + if 'extname' in ext.header and ext.header['extname'] == 'SIPWCS': + # convert from string to tuple or int + sciext = eval(ext.header['sciext']) + extlist.append(fobj[sciext]) + # determine whether distortion is the same + current_distname = hdrlet[0].header['distname'] + same_dist = True + if current_distname != fobj[0].header['distname']: + same_dist = False + if not archive and not force: + if close_fobj: + fobj.close() + message = """ + Headerlet does not have the same distortion as image! + Set "archive"=True to save old distortion model, or + set "force"=True to overwrite old model with new. + """ + logger.critical(message) + raise ValueError + + # check whether primary WCS has been archived already + # Use information from first 'SCI' extension + priwcs_name = None + + scihdr = extlist[0].header + sci_wcsnames = altwcs.wcsnames(scihdr).values() + if 'hdrname' in scihdr: + priwcs_hdrname = scihdr['hdrname'] + else: + if 'wcsname' in scihdr: + priwcs_hdrname = priwcs_name = scihdr['wcsname'] + else: + if 'idctab' in scihdr: + priwcs_hdrname = ''.join(['IDC_', + utils.extract_rootname(scihdr['idctab'], suffix='_idc')]) + else: + priwcs_hdrname = 'UNKNOWN' + priwcs_name = priwcs_hdrname + scihdr['WCSNAME'] = priwcs_name + + priwcs_unique = verify_hdrname_is_unique(fobj, priwcs_hdrname) + if archive and priwcs_unique: + if priwcs_unique: + newhdrlet = create_headerlet(fobj, sciext=scihdr['extname'], + hdrname=priwcs_hdrname) + newhdrlet.attach_to_file(fobj) + # + # copy hdrlet as a primary + # + hdrlet.apply_as_primary(fobj, attach=False, archive=archive, force=force) + + utils.updateNEXTENDKw(fobj) + fobj.flush() + if close_fobj: + fobj.close() + + +@with_logging +def restore_all_with_distname(filename, distname, primary, archive=True, + sciext='SCI', logging=False, logmode='w'): + """ + Restores all HeaderletHDUs with a given distortion model as alternate WCSs and a primary + + Parameters + -------------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + distname: string + distortion model as represented by a DISTNAME keyword + primary: int or string or None + HeaderletHDU to be restored as primary + if int - a fits extension + if string - HDRNAME + if None - use first HeaderletHDU + archive: boolean (default True) + flag indicating if HeaderletHDUs should be created from the + primary and alternate WCSs in fname before restoring all matching + headerlet extensions + logging: boolean + enable file logging + logmode: 'a' or 'w' + """ + + fobj, fname, close_fobj = parse_filename(filename, mode='update') + + hdrlet_ind = find_headerlet_HDUs(fobj, distname=distname) + if len(hdrlet_ind) == 0: + message = """ + No Headerlet extensions found with + + DISTNAME = %s in %s. + + For a full list of DISTNAMEs found in all headerlet extensions: + + get_headerlet_kw_names(fobj, kw='DISTNAME') + """ % (distname, fname) + if close_fobj: + fobj.close() + logger.critical(message) + raise ValueError + + # Interpret 'primary' parameter input into extension number + if primary is None: + primary_ind = hdrlet_ind[0] + elif isinstance(primary, int): + primary_ind = primary + else: + primary_ind = None + for ind in hdrlet_ind: + if fobj[ind].header['hdrname'] == primary: + primary_ind = ind + break + if primary_ind is None: + if close_fobj: + fobj.close() + message = """ + No Headerlet extensions found with DISTNAME = %s in %s. + """ % (primary, fname) + logger.critical(message) + raise ValueError + # Check to see whether 'primary' HeaderletHDU has same distname as user + # specified on input + + # read headerlet from HeaderletHDU into memory + if hasattr(fobj[primary_ind], 'hdulist'): + primary_hdrlet = fobj[primary_ind].hdulist + else: + primary_hdrlet = fobj[primary_ind].headerlet # older convention in PyFITS + pri_distname = primary_hdrlet[0].header['distname'] + if pri_distname != distname: + if close_fobj: + fobj.close() + message = """ + Headerlet extension to be used as PRIMARY WCS + has "DISTNAME" = %s + "DISTNAME" = %s was specified on input. + All updated WCSs must have same DISTNAME. Quitting...' + """ % (pri_distname, distname) + logger.critical(message) + raise ValueError + + # read in the names of the WCSs which the HeaderletHDUs will update + wnames = altwcs.wcsnames(fobj[sciext, 1].header) + + # work out how many HeaderletHDUs will be used to update the WCSs + numhlt = len(hdrlet_ind) + hdrnames = get_headerlet_kw_names(fobj, kw='wcsname') + + # read in headerletHDUs and update WCS keywords + for hlet in hdrlet_ind: + if fobj[hlet].header['distname'] == distname: + if hasattr(fobj[hlet], 'hdulist'): + hdrlet = fobj[hlet].hdulist + else: + hdrlet = fobj[hlet].headerlet # older convention in PyFITS + if hlet == primary_ind: + hdrlet.apply_as_primary(fobj, attach=False, + archive=archive, force=True) + else: + hdrlet.apply_as_alternate(fobj, attach=False, + wcsname=hdrlet[0].header['wcsname']) + + utils.updateNEXTENDKw(fobj) + fobj.flush() + if close_fobj: + fobj.close() + + +@with_logging +def archive_as_headerlet(filename, hdrname, sciext='SCI', + wcsname=None, wcskey=None, destim=None, + sipname=None, npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + nmatch=None, catalog=None, + logging=False, logmode='w'): + """ + Save a WCS as a headerlet extension and write it out to a file. + + This function will create a headerlet, attach it as an extension to the + science image (if it has not already been archived) then, optionally, + write out the headerlet to a separate headerlet file. + + Either wcsname or wcskey must be provided, if both are given, they must match a valid WCS + Updates wcscorr if necessary. + + Parameters + ---------- + filename: string or HDUList + Either a filename or PyFITS HDUList object for the input science file + An input filename (str) will be expanded as necessary to interpret + any environmental variables included in the filename. + hdrname: string + Unique name for this headerlet, stored as HDRNAME keyword + sciext: string + name (EXTNAME) of extension that contains WCS to be saved + wcsname: string + name of WCS to be archived, if " ": stop + wcskey: one of A...Z or " " or "PRIMARY" + if " " or "PRIMARY" - archive the primary WCS + destim: string + DESTIM keyword + if NOne, use ROOTNAME or science file name + sipname: string or None (default) + Name of unique file where the polynomial distortion coefficients were + read from. If None, the behavior is: + The code looks for a keyword 'SIPNAME' in the science header + If not found, for HST it defaults to 'IDCTAB' + If there is no SIP model the value is 'NOMODEL' + If there is a SIP model but no SIPNAME, it is set to 'UNKNOWN' + npolfile: string or None (default) + Name of a unique file where the non-polynomial distortion was stored. + If None: + The code looks for 'NPOLFILE' in science header. + If 'NPOLFILE' was not found and there is no npol model, it is set to 'NOMODEL' + If npol model exists, it is set to 'UNKNOWN' + d2imfile: string + Name of a unique file where the detector to image correction was + stored. If None: + The code looks for 'D2IMFILE' in the science header. + If 'D2IMFILE' is not found and there is no d2im correction, + it is set to 'NOMODEL' + If d2im correction exists, but 'D2IMFILE' is missing from science + header, it is set to 'UNKNOWN' + author: string + Name of user who created the headerlet, added as 'AUTHOR' keyword + to headerlet PRIMARY header + descrip: string + Short description of the solution provided by the headerlet + This description will be added as the single 'DESCRIP' keyword + to the headerlet PRIMARY header + history: filename, string or list of strings + Long (possibly multi-line) description of the solution provided + by the headerlet. These comments will be added as 'HISTORY' cards + to the headerlet PRIMARY header + If filename is specified, it will format and attach all text from + that file as the history. + logging: boolean + enable file folling + logmode: 'w' or 'a' + log file open mode + """ + + fobj, fname, close_fobj = parse_filename(filename, mode='update') + + if wcsname in [None, ' ', '', 'INDEF'] and wcskey is None: + message = """ + No valid WCS found found in %s. + A valid value for either "wcsname" or "wcskey" + needs to be specified. + """ % fname + if close_fobj: + fobj.close() + logger.critical(message) + raise ValueError + + # Translate 'wcskey' value for PRIMARY WCS to valid altwcs value of ' ' + if wcskey == 'PRIMARY': + wcskey = ' ' + wcskey = wcskey.upper() + + numhlt = countExtn(fobj, 'HDRLET') + + if wcsname is None: + scihdr = fobj[sciext, 1].header + wcsname = scihdr['wcsname'+wcskey] + + if hdrname in [None, ' ', '']: + hdrname = wcsname + + # Check to see whether or not a HeaderletHDU with this hdrname already + # exists + hdrnames = get_headerlet_kw_names(fobj) + if hdrname not in hdrnames: + hdrletobj = create_headerlet(fobj, sciext=sciext, + wcsname=wcsname, wcskey=wcskey, + hdrname=hdrname, + sipname=sipname, npolfile=npolfile, + d2imfile=d2imfile, author=author, + descrip=descrip, history=history, + nmatch=nmatch, catalog=catalog, + logging=False) + hlt_hdu = HeaderletHDU.fromheaderlet(hdrletobj) + + if destim is not None: + hlt_hdu[0].header['destim'] = destim + + fobj.append(hlt_hdu) + + utils.updateNEXTENDKw(fobj) + fobj.flush() + else: + message = """ + Headerlet with hdrname %s already archived for WCS %s + No new headerlet appended to %s . + """ % (hdrname, wcsname, fname) + logger.critical(message) + + if close_fobj: + fobj.close() + +#### Headerlet Class definitions +class Headerlet(fits.HDUList): + """ + A Headerlet class + Ref: http://mediawiki.stsci.edu/mediawiki/index.php/Telescopedia:Headerlets + """ + + def __init__(self, hdus=[], file=None, logging=False, logmode='w'): + """ + Parameters + ---------- + hdus : list + List of HDUs to be used to create the headerlet object itself + file: string + File-like object from which HDUs should be read + logging: boolean + enable file logging + logmode: 'w' or 'a' + for internal use only, indicates whether the log file + should be open in attach or write mode + """ + self.logging = logging + init_logging('class Headerlet', level=logging, mode=logmode) + + super(Headerlet, self).__init__(hdus, file=file) + + def init_attrs(self): + self.fname = self.filename() + self.hdrname = self[0].header["HDRNAME"] + self.wcsname = self[0].header["WCSNAME"] + self.upwcsver = self[0].header.get("UPWCSVER", "") + self.pywcsver = self[0].header.get("PYWCSVER", "") + self.destim = self[0].header["DESTIM"] + self.sipname = self[0].header["SIPNAME"] + self.idctab = self[0].header["IDCTAB"] + self.npolfile = self[0].header["NPOLFILE"] + self.d2imfile = self[0].header["D2IMFILE"] + self.distname = self[0].header["DISTNAME"] + + try: + self.vafactor = self[("SIPWCS", 1)].header.get("VAFACTOR", 1) #None instead of 1? + except (IndexError, KeyError): + self.vafactor = self[0].header.get("VAFACTOR", 1) #None instead of 1? + self.author = self[0].header["AUTHOR"] + self.descrip = self[0].header["DESCRIP"] + + self.fit_kws = ['HDRNAME', 'NMATCH', 'CATALOG'] + self.history = '' + # header['HISTORY'] returns an iterable of all HISTORY values + if 'HISTORY' in self[0].header: + for hist in self[0].header['HISTORY']: + self.history += hist + '\n' + + self.d2imerr = 0 + self.axiscorr = 1 + + # Overridden to support the Headerlet logging features + @classmethod + def fromfile(cls, fileobj, mode='readonly', memmap=False, + save_backup=False, logging=False, logmode='w', **kwargs): + hlet = super(cls, cls).fromfile(fileobj, mode, memmap, save_backup, + **kwargs) + if len(hlet) > 0: + hlet.init_attrs() + hlet.logging = logging + init_logging('class Headerlet', level=logging, mode=logmode) + return hlet + + @classmethod + def fromstring(cls, data, **kwargs): + hlet = super(cls, cls).fromstring(data, **kwargs) + hlet.logging = logging + init_logging('class Headerlet', level=logging, mode=logmode) + return hlet + + def apply_as_primary(self, fobj, attach=True, archive=True, force=False): + """ + Copy this headerlet as a primary WCS to fobj + + Parameters + ---------- + fobj: string, HDUList + science file to which the headerlet should be applied + attach: boolean + flag indicating if the headerlet should be attached as a + HeaderletHDU to fobj. If True checks that HDRNAME is unique + in the fobj and stops if not. + archive: boolean (default is True) + When the distortion model in the headerlet is the same as the + distortion model of the science file, this flag indicates if + the primary WCS should be saved as an alternate and a headerlet + extension. + When the distortion models do not match this flag indicates if + the current primary and alternate WCSs should be archived as + headerlet extensions and alternate WCS. + force: boolean (default is False) + When the distortion models of the headerlet and the primary do + not match, and archive is False this flag forces an update + of the primary + """ + self.hverify() + fobj, fname, close_dest = parse_filename(fobj, mode='update') + if not self.verify_dest(fobj, fname): + if close_dest: + fobj.close() + raise ValueError("Destination name does not match headerlet" + "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + + # Check to see whether the distortion model in the destination + # matches the distortion model in the headerlet being applied + + dname = self.get_destination_model(fobj) + dist_models_equal = self.equal_distmodel(dname) + if not dist_models_equal and not force: + raise ValueError("Distortion models do not match" + " To overwrite the distortion model, set force=True") + + orig_hlt_hdu = None + numhlt = countExtn(fobj, 'HDRLET') + hdrlet_extnames = get_headerlet_kw_names(fobj) + + # Insure that WCSCORR table has been created with all original + # WCS's recorded prior to adding the headerlet WCS + wcscorr.init_wcscorr(fobj) + + + ### start archive + # If archive has been specified + # regardless of whether or not the distortion models are equal... + + numsip = countExtn(self, 'SIPWCS') + sciext_list = [] + alt_hlethdu = [] + for i in range(1, numsip+1): + sipheader = self[('SIPWCS', i)].header + sciext_list.append((sipheader['TG_ENAME'], sipheader['TG_EVER'])) + target_ext = sciext_list[0] + if archive: + if 'wcsname' in fobj[target_ext].header: + hdrname = fobj[target_ext].header['WCSNAME'] + wcsname = hdrname + else: + hdrname = fobj[0].header['ROOTNAME'] + '_orig' + wcsname = None + if hdrname not in hdrlet_extnames: + # - if WCS has not been saved, write out WCS as headerlet extension + # Create a headerlet for the original Primary WCS data in the file, + # create an HDU from the original headerlet, and append it to + # the file + orig_hlt = create_headerlet(fobj, sciext=sciext_list, #[target_ext], + wcsname=wcsname, + hdrname=hdrname, + logging=self.logging) + orig_hlt_hdu = HeaderletHDU.fromheaderlet(orig_hlt) + numhlt += 1 + orig_hlt_hdu.header['EXTVER'] = numhlt + logger.info("Created headerlet %s to be attached to file" % hdrname) + else: + logger.info("Headerlet with name %s is already attached" % hdrname) + + + if dist_models_equal: + # Use the WCSNAME to determine whether or not to archive + # Primary WCS as altwcs + # wcsname = hwcs.wcs.name + scihdr = fobj[target_ext].header + if 'hdrname' in scihdr: + priwcs_name = scihdr['hdrname'] + else: + if 'wcsname' in scihdr: + priwcs_name = scihdr['wcsname'] + else: + if 'idctab' in scihdr: + priwcs_name = ''.join(['IDC_', + utils.extract_rootname(scihdr['idctab'], + suffix='_idc')]) + else: + priwcs_name = 'UNKNOWN' + nextkey = altwcs.next_wcskey(fobj, ext=target_ext) + altwcs.archiveWCS(fobj, ext=sciext_list, wcskey=nextkey, + wcsname=priwcs_name) + else: + + for hname in altwcs.wcsnames(fobj, ext=target_ext).values(): + if hname != 'OPUS' and hname not in hdrlet_extnames: + # get HeaderletHDU for alternate WCS as well + alt_hlet = create_headerlet(fobj, sciext=sciext_list, + wcsname=hname, wcskey=wcskey, + hdrname=hname, sipname=None, + npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + logging=self.logging) + numhlt += 1 + alt_hlet_hdu = HeaderletHDU.fromheaderlet(alt_hlet) + alt_hlet_hdu.header['EXTVER'] = numhlt + alt_hlethdu.append(alt_hlet_hdu) + hdrlet_extnames.append(hname) + + self._del_dest_WCS_ext(fobj) + for i in range(1, numsip+1): + target_ext = sciext_list[i-1] + self._del_dest_WCS(fobj, target_ext) + sipwcs = HSTWCS(self, ('SIPWCS', i)) + idckw = sipwcs._idc2hdr() + priwcs = sipwcs.to_fits(relax=True) + numnpol = 1 + numd2im = 1 + if sipwcs.wcs.has_cd(): + priwcs[0].header = altwcs.pc2cd(priwcs[0].header) + priwcs[0].header.extend(idckw) + if 'crder1' in sipheader: + for card in sipheader['crder*'].cards: + priwcs[0].header.set(card.keyword, card.value, card.comment, + after='WCSNAME') + # Update WCS with HDRNAME as well + + for kw in ['SIMPLE', 'BITPIX', 'NAXIS', 'EXTEND']: + try: + priwcs[0].header.remove(kw) + except ValueError: + pass + + priwcs[0].header.set('WCSNAME', self[0].header['WCSNAME'], "") + priwcs[0].header.set('WCSAXES', self[('SIPWCS', i)].header['WCSAXES'], "") + priwcs[0].header.set('HDRNAME', self[0].header['HDRNAME'], "") + if sipwcs.det2im1 or sipwcs.det2im2: + try: + d2imerr = self[('SIPWCS', i)].header['D2IMERR*'] + priwcs[0].header.extend(d2imerr) + except KeyError: + pass + try: + priwcs[0].header.append(self[('SIPWCS', i)].header.cards['D2IMEXT']) + except KeyError: + pass + if 'D2IM1.EXTVER' in priwcs[0].header: + priwcs[0].header['D2IM1.EXTVER'] = self[('SIPWCS', i)].header['D2IM1.EXTVER'] + priwcs[('D2IMARR', 1)].header['EXTVER'] = self[('SIPWCS', i)].header['D2IM1.EXTVER'] + if 'D2IM2.EXTVER' in priwcs[0].header: + priwcs[0].header['D2IM2.EXTVER'] = self[('SIPWCS', i)].header['D2IM2.EXTVER'] + priwcs[('D2IMARR', 2)].header['EXTVER'] = self[('SIPWCS', i)].header['D2IM2.EXTVER'] + # D2IM1 will NOT exist for WFPC2 data... + if 'D2IM1.EXTVER' in priwcs[0].header: + # only set number of D2IM extensions to 2 if D2IM1 exists + numd2im = 2 + + if sipwcs.cpdis1 or sipwcs.cpdis2: + try: + cperr = self[('SIPWCS', i)].header['CPERR*'] + priwcs[0].header.extend(cperr) + except KeyError: + pass + try: + priwcs[0].header.append(self[('SIPWCS', i)].header.cards['NPOLEXT']) + except KeyError: + pass + if 'DP1.EXTVER' in priwcs[0].header: + priwcs[0].header['DP1.EXTVER'] = self[('SIPWCS', i)].header['DP1.EXTVER'] + priwcs[('WCSDVARR', 1)].header['EXTVER'] = self[('SIPWCS', i)].header['DP1.EXTVER'] + if 'DP2.EXTVER' in priwcs[0].header: + priwcs[0].header['DP2.EXTVER'] = self[('SIPWCS', i)].header['DP2.EXTVER'] + priwcs[('WCSDVARR', 2)].header['EXTVER'] = self[('SIPWCS', i)].header['DP2.EXTVER'] + numnpol = 2 + + fobj[target_ext].header.extend(priwcs[0].header) + if sipwcs.cpdis1: + whdu = priwcs[('WCSDVARR', (i-1)*numnpol+1)].copy() + whdu.update_ext_version(self[('SIPWCS', i)].header['DP1.EXTVER']) + fobj.append(whdu) + if sipwcs.cpdis2: + whdu = priwcs[('WCSDVARR', i*numnpol)].copy() + whdu.update_ext_version(self[('SIPWCS', i)].header['DP2.EXTVER']) + fobj.append(whdu) + if sipwcs.det2im1: #or sipwcs.det2im2: + whdu = priwcs[('D2IMARR', (i-1)*numd2im+1)].copy() + whdu.update_ext_version(self[('SIPWCS', i)].header['D2IM1.EXTVER']) + fobj.append(whdu) + if sipwcs.det2im2: + whdu = priwcs[('D2IMARR', i*numd2im)].copy() + whdu.update_ext_version(self[('SIPWCS', i)].header['D2IM2.EXTVER']) + fobj.append(whdu) + + update_versions(self[0].header, fobj[0].header) + refs = update_ref_files(self[0].header, fobj[0].header) + # Update the WCSCORR table with new rows from the headerlet's WCSs + wcscorr.update_wcscorr(fobj, self, 'SIPWCS') + + # Append the original headerlet + if archive and orig_hlt_hdu: + fobj.append(orig_hlt_hdu) + # Append any alternate WCS Headerlets + if len(alt_hlethdu) > 0: + for ahdu in alt_hlethdu: + fobj.append(ahdu) + if attach: + # Finally, append an HDU for this headerlet + self.attach_to_file(fobj) + utils.updateNEXTENDKw(fobj) + if close_dest: + fobj.close() + + + def apply_as_alternate(self, fobj, attach=True, wcskey=None, wcsname=None): + """ + Copy this headerlet as an alternate WCS to fobj + + Parameters + ---------- + fobj: string, HDUList + science file/HDUList to which the headerlet should be applied + attach: boolean + flag indicating if the headerlet should be attached as a + HeaderletHDU to fobj. If True checks that HDRNAME is unique + in the fobj and stops if not. + wcskey: string + Key value (A-Z, except O) for this alternate WCS + If None, the next available key will be used + wcsname: string + Name to be assigned to this alternate WCS + WCSNAME is a required keyword in a Headerlet but this allows the + user to change it as desired. + + """ + self.hverify() + fobj, fname, close_dest = parse_filename(fobj, mode='update') + if not self.verify_dest(fobj, fname): + if close_dest: + fobj.close() + raise ValueError("Destination name does not match headerlet" + "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + + # Verify whether this headerlet has the same distortion + #found in the image being updated + dname = self.get_destination_model(fobj) + dist_models_equal = self.equal_distmodel(dname) + if not dist_models_equal: + raise ValueError("Distortion models do not match \n" + "Headerlet: %s \n" + "Destination file: %s\n" + "attach_to_file() can be used to append this headerlet" %(self.distname, dname)) + + # Insure that WCSCORR table has been created with all original + # WCS's recorded prior to adding the headerlet WCS + wcscorr.init_wcscorr(fobj) + + # determine value of WCSNAME to be used + if wcsname is not None: + wname = wcsname + else: + wname = self[0].header['WCSNAME'] + tg_ename = self[('SIPWCS', 1)].header['TG_ENAME'] + tg_ever = self[('SIPWCS', 1)].header['TG_EVER'] + # determine what alternate WCS this headerlet will be assigned to + if wcskey is None: + wkey = altwcs.next_wcskey(fobj[(tg_ename, tg_ever)].header) + else: + wcskey = wcskey.upper() + available_keys = altwcs.available_wcskeys(fobj[(tg_ename, tg_ever)].header) + if wcskey in available_keys: + wkey = wcskey + else: + mess = "Observation %s already contains alternate WCS with key %s" % (fname, wcskey) + logger.critical(mess) + if close_dest: + fobj.close() + raise ValueError(mess) + numsip = countExtn(self, 'SIPWCS') + for idx in range(1, numsip + 1): + siphdr = self[('SIPWCS', idx)].header + tg_ext = (siphdr['TG_ENAME'], siphdr['TG_EVER']) + + fhdr = fobj[tg_ext].header + hwcs = pywcs.WCS(siphdr, self) + hwcs_header = hwcs.to_header(key=wkey) + _idc2hdr(siphdr, fhdr, towkey=wkey) + if hwcs.wcs.has_cd(): + hwcs_header = altwcs.pc2cd(hwcs_header, key=wkey) + + fhdr.extend(hwcs_header) + fhdr['WCSNAME' + wkey] = wname + # also update with HDRNAME (a non-WCS-standard kw) + for kw in self.fit_kws: + #fhdr.insert(wind, pyfits.Card(kw + wkey, + # self[0].header[kw])) + fhdr.append(fits.Card(kw + wkey, self[0].header[kw])) + # Update the WCSCORR table with new rows from the headerlet's WCSs + wcscorr.update_wcscorr(fobj, self, 'SIPWCS') + + if attach: + self.attach_to_file(fobj) + utils.updateNEXTENDKw(fobj) + + if close_dest: + fobj.close() + + def attach_to_file(self, fobj, archive=False): + """ + Attach Headerlet as an HeaderletHDU to a science file + + Parameters + ---------- + fobj: string, HDUList + science file/HDUList to which the headerlet should be applied + archive: string + Specifies whether or not to update WCSCORR table when attaching + + Notes + ----- + The algorithm used by this method: + - verify headerlet can be applied to this file (based on DESTIM) + - verify that HDRNAME is unique for this file + - attach as HeaderletHDU to fobj + + """ + self.hverify() + fobj, fname, close_dest = parse_filename(fobj, mode='update') + destver = self.verify_dest(fobj, fname) + hdrver = self.verify_hdrname(fobj) + if destver and hdrver: + + numhlt = countExtn(fobj, 'HDRLET') + new_hlt = HeaderletHDU.fromheaderlet(self) + new_hlt.header['extver'] = numhlt + 1 + fobj.append(new_hlt) + utils.updateNEXTENDKw(fobj) + else: + message = "Headerlet %s cannot be attached to" % (self.hdrname) + message += "observation %s" % (fname) + if not destver: + message += " * Image %s keyword ROOTNAME not equal to " % (fname) + message += " DESTIM = '%s'\n" % (self.destim) + if not hdrver: + message += " * Image %s already has headerlet " % (fname) + message += "with HDRNAME='%s'\n" % (self.hdrname) + logger.critical(message) + if close_dest: + fobj.close() + + def info(self, columns=None, pad=2, maxwidth=None, + output=None, clobber=True, quiet=False): + """ + Prints a summary of this headerlet + The summary includes: + HDRNAME WCSNAME DISTNAME SIPNAME NPOLFILE D2IMFILE + + Parameters + ---------- + columns: list + List of headerlet PRIMARY header keywords to report in summary + By default (set to None), it will use the default set of keywords + defined as the global list DEFAULT_SUMMARY_COLS + pad: int + Number of padding spaces to put between printed columns + [Default: 2] + maxwidth: int + Maximum column width(not counting padding) for any column in summary + By default (set to None), each column's full width will be used + output: string (optional) + Name of optional output file to record summary. This filename + can contain environment variables. + [Default: None] + clobber: bool + If True, will overwrite any previous output file of same name + quiet: bool + If True, will NOT report info to STDOUT + + """ + summary_cols, summary_dict = self.summary(columns=columns) + print_summary(summary_cols, summary_dict, pad=pad, maxwidth=maxwidth, + idcol=None, output=output, clobber=clobber, quiet=quiet) + + def summary(self, columns=None): + """ + Returns a summary of this headerlet as a dictionary + + The summary includes a summary of the distortion model as : + HDRNAME WCSNAME DISTNAME SIPNAME NPOLFILE D2IMFILE + + Parameters + ---------- + columns: list + List of headerlet PRIMARY header keywords to report in summary + By default(set to None), it will use the default set of keywords + defined as the global list DEFAULT_SUMMARY_COLS + + Returns + ------- + summary: dict + Dictionary of values for summary + """ + if columns is None: + summary_cols = DEFAULT_SUMMARY_COLS + else: + summary_cols = columns + + # Initialize summary dict based on requested columns + summary = {} + for kw in summary_cols: + summary[kw] = copy.deepcopy(COLUMN_DICT) + + # Populate the summary with headerlet values + for kw in summary_cols: + if kw in self[0].header: + val = self[0].header[kw] + else: + val = 'INDEF' + summary[kw]['vals'].append(val) + summary[kw]['width'].append(max(len(val), len(kw))) + + return summary_cols, summary + + def hverify(self): + """ + Verify the headerlet file is a valid fits file and has + the required Primary Header keywords + """ + self.verify() + header = self[0].header + assert('DESTIM' in header and header['DESTIM'].strip()) + assert('HDRNAME' in header and header['HDRNAME'].strip()) + assert('UPWCSVER' in header) + + def verify_hdrname(self, dest): + """ + Verifies that the headerlet can be applied to the observation + + Reports whether or not this file already has a headerlet with this + HDRNAME. + """ + unique = verify_hdrname_is_unique(dest, self.hdrname) + logger.debug("verify_hdrname() returned %s"%unique) + return unique + + def get_destination_model(self, dest): + """ + Verifies that the headerlet can be applied to the observation + + Determines whether or not the file specifies the same distortion + model/reference files. + """ + destim_opened = False + if not isinstance(dest, fits.HDUList): + destim = fits.open(dest) + destim_opened = True + else: + destim = dest + dname = destim[0].header['DISTNAME'] if 'distname' in destim[0].header \ + else self.build_distname(dest) + if destim_opened: + destim.close() + return dname + + def equal_distmodel(self, dmodel): + if dmodel != self[0].header['DISTNAME']: + if self.logging: + message = """ + Distortion model in headerlet not the same as destination model + Headerlet model : %s + Destination model: %s + """ % (self[0].header['DISTNAME'], dmodel) + logger.critical(message) + return False + else: + return True + + def verify_dest(self, dest, fname): + """ + verifies that the headerlet can be applied to the observation + + DESTIM in the primary header of the headerlet must match ROOTNAME + of the science file (or the name of the destination file) + """ + try: + if not isinstance(dest, fits.HDUList): + droot = fits.getval(dest, 'ROOTNAME') + else: + droot = dest[0].header['ROOTNAME'] + except KeyError: + logger.debug("Keyword 'ROOTNAME' not found in destination file") + droot = dest.split('.fits')[0] + if droot == self.destim: + logger.debug("verify_destim() returned True") + return True + else: + logger.debug("verify_destim() returned False") + logger.critical("Destination name does not match headerlet. " + "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + return False + + def build_distname(self, dest): + """ + Builds the DISTNAME for dest based on reference file names. + """ + + try: + npolfile = dest[0].header['NPOLFILE'] + except KeyError: + npolfile = None + try: + d2imfile = dest[0].header['D2IMFILE'] + except KeyError: + d2imfile = None + + sipname, idctab = utils.build_sipname(dest, dest, None) + npolname, npolfile = utils.build_npolname(dest, npolfile) + d2imname, d2imfile = utils.build_d2imname(dest, d2imfile) + dname = utils.build_distname(sipname,npolname,d2imname) + return dname + + def tofile(self, fname, destim=None, hdrname=None, clobber=False): + """ + Write this headerlet to a file + + Parameters + ---------- + fname: string + file name + destim: string (optional) + provide a value for DESTIM keyword + hdrname: string (optional) + provide a value for HDRNAME keyword + clobber: boolean + a flag which allows to overwrte an existing file + """ + if not destim or not hdrname: + self.hverify() + self.writeto(fname, clobber=clobber) + + def _del_dest_WCS(self, dest, ext=None): + """ + Delete the WCS of a science file extension + """ + + logger.info("Deleting all WCSs of file %s" % dest.filename()) + numext = len(dest) + + if ext: + fext = dest[ext] + self._remove_d2im(fext) + self._remove_sip(fext) + self._remove_lut(fext) + self._remove_primary_WCS(fext) + self._remove_idc_coeffs(fext) + self._remove_fit_values(fext) + else: + for idx in range(numext): + # Only delete WCS from extensions which may have WCS keywords + if ('XTENSION' in dest[idx].header and + dest[idx].header['XTENSION'] == 'IMAGE'): + self._remove_d2im(dest[idx]) + self._remove_sip(dest[idx]) + self._remove_lut(dest[idx]) + self._remove_primary_WCS(dest[idx]) + self._remove_idc_coeffs(dest[idx]) + self._remove_fit_values(dest[idx]) + self._remove_ref_files(dest[0]) + """ + if not ext: + self._remove_alt_WCS(dest, ext=range(numext)) + else: + self._remove_alt_WCS(dest, ext=ext) + """ + def _del_dest_WCS_ext(self, dest): + numwdvarr = countExtn(dest, 'WCSDVARR') + numd2im = countExtn(dest, 'D2IMARR') + if numwdvarr > 0: + for idx in range(1, numwdvarr + 1): + del dest[('WCSDVARR', idx)] + if numd2im > 0: + for idx in range(1, numd2im + 1): + del dest[('D2IMARR', idx)] + + def _remove_ref_files(self, phdu): + """ + phdu: Primary HDU + """ + refkw = ['IDCTAB', 'NPOLFILE', 'D2IMFILE', 'SIPNAME', 'DISTNAME'] + for kw in refkw: + try: + del phdu.header[kw] + except KeyError: + pass + + def _remove_fit_values(self, ext): + """ + Remove the any existing astrometric fit values from a FITS extension + """ + + logger.debug("Removing astrometric fit values from (%s, %s)"% + (ext.name, ext.ver)) + dkeys = altwcs.wcskeys(ext.header) + if 'O' in dkeys: dkeys.remove('O') # Do not remove wcskey='O' values + for fitkw in ['NMATCH', 'CATALOG']: + for k in dkeys: + fkw = (fitkw+k).rstrip() + if fkw in ext.header: + del ext.header[fkw] + + def _remove_sip(self, ext): + """ + Remove the SIP distortion of a FITS extension + """ + + logger.debug("Removing SIP distortion from (%s, %s)" + % (ext.name, ext.ver)) + for prefix in ['A', 'B', 'AP', 'BP']: + try: + order = ext.header[prefix + '_ORDER'] + del ext.header[prefix + '_ORDER'] + except KeyError: + continue + for i in range(order + 1): + for j in range(order + 1): + key = prefix + '_%d_%d' % (i, j) + try: + del ext.header[key] + except KeyError: + pass + try: + del ext.header['IDCTAB'] + except KeyError: + pass + + def _remove_lut(self, ext): + """ + Remove the Lookup Table distortion of a FITS extension + """ + + logger.debug("Removing LUT distortion from (%s, %s)" + % (ext.name, ext.ver)) + try: + cpdis = ext.header['CPDIS*'] + except KeyError: + return + try: + for c in range(1, len(cpdis) + 1): + del ext.header['DP%s*...' % c] + del ext.header[cpdis.cards[c - 1].keyword] + del ext.header['CPERR*'] + del ext.header['NPOLFILE'] + del ext.header['NPOLEXT'] + except KeyError: + pass + + def _remove_d2im(self, ext): + """ + Remove the Detector to Image correction of a FITS extension + """ + + logger.debug("Removing D2IM correction from (%s, %s)" + % (ext.name, ext.ver)) + try: + d2imdis = ext.header['D2IMDIS*'] + except KeyError: + return + try: + for c in range(1, len(d2imdis) + 1): + del ext.header['D2IM%s*...' % c] + del ext.header[d2imdis.cards[c - 1].keyword] + del ext.header['D2IMERR*'] + del ext.header['D2IMFILE'] + del ext.header['D2IMEXT'] + except KeyError: + pass + + def _remove_alt_WCS(self, dest, ext): + """ + Remove Alternate WCSs of a FITS extension. + A WCS with wcskey 'O' is never deleted. + """ + dkeys = altwcs.wcskeys(dest[('SCI', 1)].header) + for val in ['O', '', ' ']: + if val in dkeys: + dkeys.remove(val) # Never delete WCS with wcskey='O' + + logger.debug("Removing alternate WCSs with keys %s from %s" + % (dkeys, dest.filename())) + for k in dkeys: + altwcs.deleteWCS(dest, ext=ext, wcskey=k) + + def _remove_primary_WCS(self, ext): + """ + Remove the primary WCS of a FITS extension + """ + + logger.debug("Removing Primary WCS from (%s, %s)" + % (ext.name, ext.ver)) + naxis = ext.header['NAXIS'] + for key in basic_wcs: + for i in range(1, naxis + 1): + try: + del ext.header[key + str(i)] + except KeyError: + pass + try: + del ext.header['WCSAXES'] + except KeyError: + pass + try: + del ext.header['WCSNAME'] + except KeyError: + pass + + def _remove_idc_coeffs(self, ext): + """ + Remove IDC coefficients of a FITS extension + """ + + logger.debug("Removing IDC coefficient from (%s, %s)" + % (ext.name, ext.ver)) + coeffs = ['OCX10', 'OCX11', 'OCY10', 'OCY11', 'IDCSCALE'] + for k in coeffs: + try: + del ext.header[k] + except KeyError: + pass + +@with_logging +def _idc2hdr(fromhdr, tohdr, towkey=' '): + """ + Copy the IDC (HST specific) keywords from one header to another + + """ + # save some of the idc coefficients + coeffs = ['OCX10', 'OCX11', 'OCY10', 'OCY11', 'IDCSCALE'] + for c in coeffs: + try: + tohdr[c+towkey] = fromhdr[c] + logger.debug("Copied %s to header") + except KeyError: + continue + + +def get_extname_extver_list(fobj, sciext): + """ + Create a list of (EXTNAME, EXTVER) tuples + + Based on sciext keyword (see docstring for create_headerlet) + walk throughh the file and convert extensions in `sciext` to + valid (EXTNAME, EXTVER) tuples. + """ + extlist = [] + if isinstance(sciext, int): + if sciext == 0: + extname = 'PRIMARY' + extver = 1 + else: + try: + extname = fobj[sciext].header['EXTNAME'] + except KeyError: + extname = "" + try: + extver = fobj[sciext].header['EXTVER'] + except KeyError: + extver = 1 + extlist.append((extname, extver)) + elif isinstance(sciext, str): + if sciext == 'PRIMARY': + extname = "PRIMARY" + extver = 1 + extlist.append((extname, extver)) + else: + for ext in fobj: + try: + extname = ext.header['EXTNAME'] + except KeyError: + continue + if extname.upper() == sciext.upper(): + try: + extver = ext.header['EXTVER'] + except KeyError: + extver = 1 + extlist.append((extname, extver)) + elif isinstance(sciext, list): + if isinstance(sciext[0], int): + for i in sciext: + try: + extname = fobj[i].header['EXTNAME'] + except KeyError: + if i == 0: + extname = "PRIMARY" + extver = 1 + else: + extname = "" + try: + extver = fobj[i].header['EXTVER'] + except KeyError: + extver = 1 + extlist.append((extname, extver)) + else: + extlist = sciext[:] + else: + errstr = "Expected sciext to be a list of FITS extensions with science data\n"+\ + " a valid EXTNAME string, or an integer." + logger.critical(errstr) + raise ValueError + return extlist + + +class HeaderletHDU(fits.hdu.nonstandard.FitsHDU): + """ + A non-standard extension HDU for encapsulating Headerlets in a file. These + HDUs have an extension type of HDRLET and their EXTNAME is derived from the + Headerlet's HDRNAME. + + The data itself is a FITS file embedded within the HDU data. The file name + is derived from the HDRNAME keyword, and should be in the form + `<HDRNAME>_hdr.fits`. If the COMPRESS keyword evaluates to `True`, the tar + file is compressed with gzip compression. + + The structure of this HDU is the same as that proposed for the 'FITS' + extension type proposed here: + http://listmgr.cv.nrao.edu/pipermail/fitsbits/2002-April/thread.html + + The Headerlet contained in the HDU's data can be accessed by the + `headerlet` attribute. + """ + + _extension = 'HDRLET' + + @lazyproperty + def headerlet(self): + """Return the encapsulated headerlet as a Headerlet object. + + This is similar to the hdulist property inherited from the FitsHDU + class, though the hdulist property returns a normal HDUList object. + """ + + return Headerlet(self.hdulist) + + @classmethod + def fromheaderlet(cls, headerlet, compress=False): + """ + Creates a new HeaderletHDU from a given Headerlet object. + + Parameters + ---------- + headerlet : `Headerlet` + A valid Headerlet object. + + compress : bool, optional + Gzip compress the headerlet data. + + Returns + ------- + hlet : `HeaderletHDU` + A `HeaderletHDU` object for the given `Headerlet` that can be + attached as an extension to an existing `HDUList`. + """ + + # TODO: Perhaps check that the given object is in fact a valid + # Headerlet + hlet = cls.fromhdulist(headerlet, compress) + + # Add some more headerlet-specific keywords to the header + phdu = headerlet[0] + + if 'SIPNAME' in phdu.header: + sipname = phdu.header['SIPNAME'] + else: + sipname = phdu.header['WCSNAME'] + + hlet.header['HDRNAME'] = (phdu.header['HDRNAME'], + phdu.header.comments['HDRNAME']) + hlet.header['DATE'] = (phdu.header['DATE'], + phdu.header.comments['DATE']) + hlet.header['SIPNAME'] = (sipname, 'SIP distortion model name') + hlet.header['WCSNAME'] = (phdu.header['WCSNAME'], 'WCS name') + hlet.header['DISTNAME'] = (phdu.header['DISTNAME'], + 'Distortion model name') + hlet.header['NPOLFILE'] = (phdu.header['NPOLFILE'], + phdu.header.comments['NPOLFILE']) + hlet.header['D2IMFILE'] = (phdu.header['D2IMFILE'], + phdu.header.comments['D2IMFILE']) + hlet.header['EXTNAME'] = (cls._extension, 'Extension name') + + return hlet + + +fits.register_hdu(HeaderletHDU) |