summaryrefslogtreecommitdiff
path: root/stwcs/updatewcs
diff options
context:
space:
mode:
Diffstat (limited to 'stwcs/updatewcs')
-rw-r--r--stwcs/updatewcs/__init__.py376
-rw-r--r--stwcs/updatewcs/apply_corrections.py248
-rw-r--r--stwcs/updatewcs/corrections.py326
-rw-r--r--stwcs/updatewcs/det2im.py299
-rw-r--r--stwcs/updatewcs/makewcs.py273
-rw-r--r--stwcs/updatewcs/npol.py343
-rw-r--r--stwcs/updatewcs/utils.py264
-rw-r--r--stwcs/updatewcs/wfpc2_dgeo.py123
8 files changed, 2252 insertions, 0 deletions
diff --git a/stwcs/updatewcs/__init__.py b/stwcs/updatewcs/__init__.py
new file mode 100644
index 0000000..eb5e850
--- /dev/null
+++ b/stwcs/updatewcs/__init__.py
@@ -0,0 +1,376 @@
+from __future__ import absolute_import, division, print_function # confidence high
+
+from astropy.io import fits
+from stwcs import wcsutil
+from stwcs.wcsutil import HSTWCS
+import stwcs
+
+from astropy import wcs as pywcs
+import astropy
+
+from . import utils, corrections, makewcs
+from . import npol, det2im
+from stsci.tools import parseinput, fileutil
+from . import apply_corrections
+
+import time
+import logging
+logger = logging.getLogger('stwcs.updatewcs')
+
+import atexit
+atexit.register(logging.shutdown)
+
+#Note: The order of corrections is important
+
+def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True,
+ checkfiles=True, verbose=False):
+ """
+
+ Updates HST science files with the best available calibration information.
+ This allows users to retrieve from the archive self contained science files
+ which do not require additional reference files.
+
+ Basic WCS keywords are updated in the process and new keywords (following WCS
+ Paper IV and the SIP convention) as well as new extensions are added to the science files.
+
+
+ Example
+ -------
+ >>>from stwcs import updatewcs
+ >>>updatewcs.updatewcs(filename)
+
+ Dependencies
+ ------------
+ `stsci.tools`
+ `astropy.io.fits`
+ `astropy.wcs`
+
+ Parameters
+ ----------
+ input: a python list of file names or a string (wild card characters allowed)
+ input files may be in fits, geis or waiver fits format
+ vacorr: boolean
+ If True, vecocity aberration correction will be applied
+ tddcorr: boolean
+ If True, time dependent distortion correction will be applied
+ npolcorr: boolean
+ If True, a Lookup table distortion will be applied
+ d2imcorr: boolean
+ If True, detector to image correction will be applied
+ checkfiles: boolean
+ If True, the format of the input files will be checked,
+ geis and waiver fits files will be converted to MEF format.
+ Default value is True for standalone mode.
+ """
+ if verbose == False:
+ logger.setLevel(100)
+ else:
+ formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
+ log_filename = 'stwcs.log'
+ fh = logging.FileHandler(log_filename, mode='w')
+ fh.setLevel(logging.DEBUG)
+ fh.setFormatter(formatter)
+ logger.addHandler(fh)
+ logger.setLevel(verbose)
+ args = "vacorr=%s, tddcorr=%s, npolcorr=%s, d2imcorr=%s, checkfiles=%s, \
+ " % (str(vacorr), str(tddcorr), str(npolcorr),
+ str(d2imcorr), str(checkfiles))
+ logger.info('\n\tStarting UPDATEWCS: %s', time.asctime())
+
+ files = parseinput.parseinput(input)[0]
+ logger.info("\n\tInput files: %s, " % [i for i in files])
+ logger.info("\n\tInput arguments: %s" %args)
+ if checkfiles:
+ files = checkFiles(files)
+ if not files:
+ print('No valid input, quitting ...\n')
+ return
+
+ for f in files:
+ acorr = apply_corrections.setCorrections(f, vacorr=vacorr, \
+ tddcorr=tddcorr,npolcorr=npolcorr, d2imcorr=d2imcorr)
+ if 'MakeWCS' in acorr and newIDCTAB(f):
+ logger.warning("\n\tNew IDCTAB file detected. All current WCSs will be deleted")
+ cleanWCS(f)
+
+ makecorr(f, acorr)
+
+ return files
+
+def makecorr(fname, allowed_corr):
+ """
+ Purpose
+ =======
+ Applies corrections to the WCS of a single file
+
+ :Parameters:
+ `fname`: string
+ file name
+ `acorr`: list
+ list of corrections to be applied
+ """
+ logger.info("Allowed corrections: {0}".format(allowed_corr))
+ f = fits.open(fname, mode='update')
+ #Determine the reference chip and create the reference HSTWCS object
+ nrefchip, nrefext = getNrefchip(f)
+ wcsutil.restoreWCS(f, nrefext, wcskey='O')
+ rwcs = HSTWCS(fobj=f, ext=nrefext)
+ rwcs.readModel(update=True,header=f[nrefext].header)
+
+ if 'DET2IMCorr' in allowed_corr:
+ kw2update = det2im.DET2IMCorr.updateWCS(f)
+ for kw in kw2update:
+ f[1].header[kw] = kw2update[kw]
+
+ for i in range(len(f))[1:]:
+ extn = f[i]
+
+ if 'extname' in extn.header:
+ extname = extn.header['extname'].lower()
+ if extname == 'sci':
+ wcsutil.restoreWCS(f, ext=i, wcskey='O')
+ sciextver = extn.header['extver']
+ ref_wcs = rwcs.deepcopy()
+ hdr = extn.header
+ ext_wcs = HSTWCS(fobj=f, ext=i)
+ ### check if it exists first!!!
+ # 'O ' can be safely archived again because it has been restored first.
+ wcsutil.archiveWCS(f, ext=i, wcskey="O", wcsname="OPUS", reusekey=True)
+ ext_wcs.readModel(update=True,header=hdr)
+ for c in allowed_corr:
+ if c != 'NPOLCorr' and c != 'DET2IMCorr':
+ corr_klass = corrections.__getattribute__(c)
+ kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs)
+ for kw in kw2update:
+ hdr[kw] = kw2update[kw]
+ # give the primary WCS a WCSNAME value
+ idcname = f[0].header.get('IDCTAB', " ")
+ if idcname.strip() and 'idc.fits' in idcname:
+ wname = ''.join(['IDC_',
+ utils.extract_rootname(idcname,suffix='_idc')])
+ else: wname = " "
+ hdr['WCSNAME'] = wname
+
+ elif extname in ['err', 'dq', 'sdq', 'samp', 'time']:
+ cextver = extn.header['extver']
+ if cextver == sciextver:
+ hdr = f[('SCI',sciextver)].header
+ w = pywcs.WCS(hdr, f)
+ copyWCS(w, extn.header)
+
+ else:
+ continue
+
+ if 'NPOLCorr' in allowed_corr:
+ kw2update = npol.NPOLCorr.updateWCS(f)
+ for kw in kw2update:
+ f[1].header[kw] = kw2update[kw]
+ # Finally record the version of the software which updated the WCS
+ if 'HISTORY' in f[0].header:
+ f[0].header.set('UPWCSVER', value=stwcs.__version__,
+ comment="Version of STWCS used to updated the WCS",
+ before='HISTORY')
+ f[0].header.set('PYWCSVER', value=astropy.__version__,
+ comment="Version of PYWCS used to updated the WCS",
+ before='HISTORY')
+ elif 'ASN_MTYP' in f[0].header:
+ f[0].header.set('UPWCSVER', value=stwcs.__version__,
+ comment="Version of STWCS used to updated the WCS",
+ after='ASN_MTYP')
+ f[0].header.set('PYWCSVER', value=astropy.__version__,
+ comment="Version of PYWCS used to updated the WCS",
+ after='ASN_MTYP')
+ else:
+ # Find index of last non-blank card, and insert this new keyword after that card
+ for i in range(len(f[0].header) - 1, 0, -1):
+ if f[0].header[i].strip() != '':
+ break
+ f[0].header.set('UPWCSVER', stwcs.__version__,
+ "Version of STWCS used to updated the WCS",
+ after=i)
+ f[0].header.set('PYWCSVER', astropy.__version__,
+ "Version of PYWCS used to updated the WCS",
+ after=i)
+ # add additional keywords to be used by headerlets
+ distdict = utils.construct_distname(f,rwcs)
+ f[0].header['DISTNAME'] = distdict['DISTNAME']
+ f[0].header['SIPNAME'] = distdict['SIPNAME']
+ # Make sure NEXTEND keyword remains accurate
+ f[0].header['NEXTEND'] = len(f)-1
+ f.close()
+
+def copyWCS(w, ehdr):
+ """
+ This is a convenience function to copy a WCS object
+ to a header as a primary WCS. It is used only to copy the
+ WCS of the 'SCI' extension to the headers of 'ERR', 'DQ', 'SDQ',
+ 'TIME' or 'SAMP' extensions.
+ """
+ hwcs = w.to_header()
+
+ if w.wcs.has_cd():
+ wcsutil.pc2cd(hwcs)
+ for k in hwcs.keys():
+ key = k[:7]
+ ehdr[key] = hwcs[k]
+
+def getNrefchip(fobj):
+ """
+
+ Finds which FITS extension holds the reference chip.
+
+ The reference chip has EXTNAME='SCI', can be in any extension and
+ is instrument specific. This functions provides mappings between
+ sci extensions, chips and fits extensions.
+ In the case of a subarray when the reference chip is missing, the
+ first 'SCI' extension is the reference chip.
+
+ Parameters
+ ----------
+ fobj: `astropy.io.fits.HDUList` object
+ """
+ nrefext = 1
+ nrefchip = 1
+ instrument = fobj[0].header['INSTRUME']
+
+ if instrument == 'WFPC2':
+ chipkw = 'DETECTOR'
+ extvers = [("SCI",img.header['EXTVER']) for img in
+ fobj[1:] if img.header['EXTNAME'].lower()=='sci']
+ detectors = [img.header[chipkw] for img in fobj[1:] if
+ img.header['EXTNAME'].lower()=='sci']
+ fitsext = [i for i in range(len(fobj))[1:] if
+ fobj[i].header['EXTNAME'].lower()=='sci']
+ det2ext=dict(list(zip(detectors, extvers)))
+ ext2det=dict(list(zip(extvers, detectors)))
+ ext2fitsext=dict(list(zip(extvers, fitsext)))
+
+ if 3 not in detectors:
+ nrefchip = ext2det.pop(extvers[0])
+ nrefext = ext2fitsext.pop(extvers[0])
+ else:
+ nrefchip = 3
+ extname = det2ext.pop(nrefchip)
+ nrefext = ext2fitsext.pop(extname)
+
+ elif (instrument == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC') or \
+ (instrument == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS'):
+ chipkw = 'CCDCHIP'
+ extvers = [("SCI",img.header['EXTVER']) for img in
+ fobj[1:] if img.header['EXTNAME'].lower()=='sci']
+ detectors = [img.header[chipkw] for img in fobj[1:] if
+ img.header['EXTNAME'].lower()=='sci']
+ fitsext = [i for i in range(len(fobj))[1:] if
+ fobj[i].header['EXTNAME'].lower()=='sci']
+ det2ext=dict(list(zip(detectors, extvers)))
+ ext2det=dict(list(zip(extvers, detectors)))
+ ext2fitsext=dict(list(zip(extvers, fitsext)))
+
+ if 2 not in detectors:
+ nrefchip = ext2det.pop(extvers[0])
+ nrefext = ext2fitsext.pop(extvers[0])
+ else:
+ nrefchip = 2
+ extname = det2ext.pop(nrefchip)
+ nrefext = ext2fitsext.pop(extname)
+ else:
+ for i in range(len(fobj)):
+ extname = fobj[i].header.get('EXTNAME', None)
+ if extname != None and extname.lower == 'sci':
+ nrefext = i
+ break
+
+ return nrefchip, nrefext
+
+def checkFiles(input):
+ """
+ Checks that input files are in the correct format.
+ Converts geis and waiver fits files to multiextension fits.
+ """
+ from stsci.tools.check_files import geis2mef, waiver2mef, checkFiles
+ logger.info("\n\tChecking files %s" % input)
+ removed_files = []
+ newfiles = []
+ if not isinstance(input, list):
+ input=[input]
+
+ for file in input:
+ try:
+ imgfits,imgtype = fileutil.isFits(file)
+ except IOError:
+ logger.warning( "\n\tFile %s could not be found, removing it from the input list.\n" %file)
+ removed_files.append(file)
+ continue
+ # Check for existence of waiver FITS input, and quit if found.
+ # Or should we print a warning and continue but not use that file
+ if imgfits:
+ if imgtype == 'waiver':
+ newfilename = waiver2mef(file, convert_dq=True)
+ if newfilename == None:
+ logger.warning("\n\tRemoving file %s from input list - could not convert waiver to mef" %file)
+ removed_files.append(file)
+ else:
+ newfiles.append(newfilename)
+ else:
+ newfiles.append(file)
+
+ # If a GEIS image is provided as input, create a new MEF file with
+ # a name generated using 'buildFITSName()'
+ # Convert the corresponding data quality file if present
+ if not imgfits:
+ newfilename = geis2mef(file, convert_dq=True)
+ if newfilename == None:
+ logger.warning("\n\tRemoving file %s from input list - could not convert geis to mef" %file)
+ removed_files.append(file)
+ else:
+ newfiles.append(newfilename)
+ if removed_files:
+ logger.warning('\n\tThe following files will be removed from the list of files to be processed %s' % removed_files)
+
+ newfiles = checkFiles(newfiles)[0]
+ logger.info("\n\tThese files passed the input check and will be processed: %s" % newfiles)
+ return newfiles
+
+def newIDCTAB(fname):
+ #When this is called we know there's a kw IDCTAB in the header
+ hdul = fits.open(fname)
+ idctab = fileutil.osfn(hdul[0].header['IDCTAB'])
+ try:
+ #check for the presence of IDCTAB in the first extension
+ oldidctab = fileutil.osfn(hdul[1].header['IDCTAB'])
+ except KeyError:
+ return False
+ if idctab == oldidctab:
+ return False
+ else:
+ return True
+
+def cleanWCS(fname):
+ # A new IDCTAB means all previously computed WCS's are invalid
+ # We are deleting all of them except the original OPUS WCS.
+ f = fits.open(fname, mode='update')
+ keys = wcsutil.wcskeys(f[1].header)
+ # Remove the primary WCS from the list
+ try:
+ keys.remove(' ')
+ except ValueError:
+ pass
+ fext = list(range(1, len(f)))
+ for key in keys:
+ try:
+ wcsutil.deleteWCS(fname, ext=fext, wcskey=key)
+ except KeyError:
+ # Some extensions don't have the alternate (or any) WCS keywords
+ continue
+
+def getCorrections(instrument):
+ """
+ Print corrections available for an instrument
+
+ :Parameters:
+ `instrument`: string, one of 'WFPC2', 'NICMOS', 'STIS', 'ACS', 'WFC3'
+ """
+ acorr = apply_corrections.allowed_corrections[instrument]
+
+ print("The following corrections will be performed for instrument %s\n" % instrument)
+ for c in acorr: print(c,': ' , apply_corrections.cnames[c])
diff --git a/stwcs/updatewcs/apply_corrections.py b/stwcs/updatewcs/apply_corrections.py
new file mode 100644
index 0000000..86b623a
--- /dev/null
+++ b/stwcs/updatewcs/apply_corrections.py
@@ -0,0 +1,248 @@
+from __future__ import division, print_function # confidence high
+
+import os
+from astropy.io import fits
+import time
+from stsci.tools import fileutil
+import os.path
+from stwcs.wcsutil import altwcs
+from . import utils
+from . import wfpc2_dgeo
+
+import logging
+logger = logging.getLogger("stwcs.updatewcs.apply_corrections")
+
+#Note: The order of corrections is important
+
+
+# A dictionary which lists the allowed corrections for each instrument.
+# These are the default corrections applied also in the pipeline.
+
+allowed_corrections={'WFPC2': ['DET2IMCorr', 'MakeWCS','CompSIP', 'VACorr'],
+ 'ACS': ['DET2IMCorr', 'TDDCorr', 'MakeWCS', 'CompSIP','VACorr', 'NPOLCorr'],
+ 'STIS': ['MakeWCS', 'CompSIP','VACorr'],
+ 'NICMOS': ['MakeWCS', 'CompSIP','VACorr'],
+ 'WFC3': ['DET2IMCorr', 'MakeWCS', 'CompSIP','VACorr', 'NPOLCorr'],
+ }
+
+cnames = {'DET2IMCorr': 'Detector to Image Correction',
+ 'TDDCorr': 'Time Dependent Distortion Correction',
+ 'MakeWCS': 'Recalculate basic WCS keywords based on the distortion model',
+ 'CompSIP': 'Given IDCTAB distortion model calculate the SIP coefficients',
+ 'VACorr': 'Velocity Aberration Correction',
+ 'NPOLCorr': 'Lookup Table Distortion'
+ }
+
+def setCorrections(fname, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True):
+ """
+ Creates a list of corrections to be applied to a file
+ based on user input paramters and allowed corrections
+ for the instrument.
+ """
+ instrument = fits.getval(fname, 'INSTRUME')
+ # make a copy of this list !
+ acorr = allowed_corrections[instrument][:]
+
+ # For WFPC2 images, the old-style DGEOFILE needs to be
+ # converted on-the-fly into a proper D2IMFILE here...
+ if instrument == 'WFPC2':
+ # check for DGEOFILE, and convert it to D2IMFILE if found
+ d2imfile = wfpc2_dgeo.update_wfpc2_d2geofile(fname)
+ # Check if idctab is present on disk
+ # If kw IDCTAB is present in the header but the file is
+ # not found on disk, do not run TDDCorr, MakeCWS and CompSIP
+ if not foundIDCTAB(fname):
+ if 'TDDCorr' in acorr: acorr.remove('TDDCorr')
+ if 'MakeWCS' in acorr: acorr.remove('MakeWCS')
+ if 'CompSIP' in acorr: acorr.remove('CompSIP')
+
+ if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr')
+ if 'TDDCorr' in acorr:
+ tddcorr = applyTDDCorr(fname, tddcorr)
+ if tddcorr == False: acorr.remove('TDDCorr')
+
+ if 'NPOLCorr' in acorr:
+ npolcorr = applyNpolCorr(fname, npolcorr)
+ if npolcorr == False: acorr.remove('NPOLCorr')
+ if 'DET2IMCorr' in acorr:
+ d2imcorr = applyD2ImCorr(fname, d2imcorr)
+ if d2imcorr == False: acorr.remove('DET2IMCorr')
+ logger.info("\n\tCorrections to be applied to %s: %s " % (fname, str(acorr)))
+ return acorr
+
+
+def foundIDCTAB(fname):
+ """
+ This functions looks for an "IDCTAB" keyword in the primary header.
+
+ Returns
+ -------
+ status : bool
+ If False : MakeWCS, CompSIP and TDDCorr should not be applied.
+ If True : there's no restriction on corrections, they all should be applied.
+
+ Raises
+ ------
+ IOError : If IDCTAB file not found on disk.
+ """
+
+ try:
+ idctab = fits.getval(fname, 'IDCTAB').strip()
+ if idctab == 'N/A' or idctab == "":
+ return False
+ except KeyError:
+ return False
+ idctab = fileutil.osfn(idctab)
+ if os.path.exists(idctab):
+ return True
+ else:
+ raise IOError("IDCTAB file {0} not found".format(idctab))
+
+
+def applyTDDCorr(fname, utddcorr):
+ """
+ The default value of tddcorr for all ACS images is True.
+ This correction will be performed if all conditions below are True:
+ - the user did not turn it off on the command line
+ - the detector is WFC
+ - the idc table specified in the primary header is available.
+ """
+
+ phdr = fits.getheader(fname)
+ instrument = phdr['INSTRUME']
+ try:
+ detector = phdr['DETECTOR']
+ except KeyError:
+ detector = None
+ try:
+ tddswitch = phdr['TDDCORR']
+ except KeyError:
+ tddswitch = 'PERFORM'
+
+ if instrument == 'ACS' and detector == 'WFC' and utddcorr == True and tddswitch == 'PERFORM':
+ tddcorr = True
+ try:
+ idctab = phdr['IDCTAB']
+ except KeyError:
+ tddcorr = False
+ #print "***IDCTAB keyword not found - not applying TDD correction***\n"
+ if os.path.exists(fileutil.osfn(idctab)):
+ tddcorr = True
+ else:
+ tddcorr = False
+ #print "***IDCTAB file not found - not applying TDD correction***\n"
+ else:
+ tddcorr = False
+ return tddcorr
+
+def applyNpolCorr(fname, unpolcorr):
+ """
+ Determines whether non-polynomial distortion lookup tables should be added
+ as extensions to the science file based on the 'NPOLFILE' keyword in the
+ primary header and NPOLEXT kw in the first extension.
+ This is a default correction and will always run in the pipeline.
+ The file used to generate the extensions is
+ recorded in the NPOLEXT keyword in the first science extension.
+ If 'NPOLFILE' in the primary header is different from 'NPOLEXT' in the
+ extension header and the file exists on disk and is a 'new type' npolfile,
+ then the lookup tables will be updated as 'WCSDVARR' extensions.
+ """
+ applyNPOLCorr = True
+ try:
+ # get NPOLFILE kw from primary header
+ fnpol0 = fits.getval(fname, 'NPOLFILE')
+ if fnpol0 == 'N/A':
+ utils.remove_distortion(fname, "NPOLFILE")
+ return False
+ fnpol0 = fileutil.osfn(fnpol0)
+ if not fileutil.findFile(fnpol0):
+ msg = """\n\tKw "NPOLFILE" exists in primary header but file %s not found
+ Non-polynomial distortion correction will not be applied\n
+ """ % fnpol0
+ logger.critical(msg)
+ raise IOError("NPOLFILE {0} not found".format(fnpol0))
+ try:
+ # get NPOLEXT kw from first extension header
+ fnpol1 = fits.getval(fname, 'NPOLEXT', ext=1)
+ fnpol1 = fileutil.osfn(fnpol1)
+ if fnpol1 and fileutil.findFile(fnpol1):
+ if fnpol0 != fnpol1:
+ applyNPOLCorr = True
+ else:
+ msg = """\n\tNPOLEXT with the same value as NPOLFILE found in first extension.
+ NPOL correction will not be applied."""
+ logger.info(msg)
+ applyNPOLCorr = False
+ else:
+ # npl file defined in first extension may not be found
+ # but if a valid kw exists in the primary header, non-polynomial
+ #distortion correction should be applied.
+ applyNPOLCorr = True
+ except KeyError:
+ # the case of "NPOLFILE" kw present in primary header but "NPOLEXT" missing
+ # in first extension header
+ applyNPOLCorr = True
+ except KeyError:
+ logger.info('\n\t"NPOLFILE" keyword not found in primary header')
+ applyNPOLCorr = False
+ return applyNPOLCorr
+
+ if isOldStyleDGEO(fname, fnpol0):
+ applyNPOLCorr = False
+ return (applyNPOLCorr and unpolcorr)
+
+def isOldStyleDGEO(fname, dgname):
+ # checks if the file defined in a NPOLFILE kw is a full size
+ # (old style) image
+
+ sci_hdr = fits.getheader(fname, ext=1)
+ dgeo_hdr = fits.getheader(dgname, ext=1)
+ sci_naxis1 = sci_hdr['NAXIS1']
+ sci_naxis2 = sci_hdr['NAXIS2']
+ dg_naxis1 = dgeo_hdr['NAXIS1']
+ dg_naxis2 = dgeo_hdr['NAXIS2']
+ if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2:
+ msg = """\n\tOnly full size (old style) DGEO file was found.\n
+ Non-polynomial distortion correction will not be applied."""
+ logger.critical(msg)
+ return True
+ else:
+ return False
+
+def applyD2ImCorr(fname, d2imcorr):
+ applyD2IMCorr = True
+ try:
+ # get D2IMFILE kw from primary header
+ fd2im0 = fits.getval(fname, 'D2IMFILE')
+ if fd2im0 == 'N/A':
+ utils.remove_distortion(fname, "D2IMFILE")
+ return False
+ fd2im0 = fileutil.osfn(fd2im0)
+ if not fileutil.findFile(fd2im0):
+ msg = """\n\tKw D2IMFILE exists in primary header but file %s not found\n
+ Detector to image correction will not be applied\n""" % fd2im0
+ logger.critical(msg)
+ print(msg)
+ raise IOError("D2IMFILE {0} not found".format(fd2im0))
+ try:
+ # get D2IMEXT kw from first extension header
+ fd2imext = fits.getval(fname, 'D2IMEXT', ext=1)
+ fd2imext = fileutil.osfn(fd2imext)
+ if fd2imext and fileutil.findFile(fd2imext):
+ if fd2im0 != fd2imext:
+ applyD2IMCorr = True
+ else:
+ applyD2IMCorr = False
+ else:
+ # D2IM file defined in first extension may not be found
+ # but if a valid kw exists in the primary header,
+ # detector to image correction should be applied.
+ applyD2IMCorr = True
+ except KeyError:
+ # the case of D2IMFILE kw present in primary header but D2IMEXT missing
+ # in first extension header
+ applyD2IMCorr = True
+ except KeyError:
+ print('D2IMFILE keyword not found in primary header')
+ applyD2IMCorr = False
+ return applyD2IMCorr
diff --git a/stwcs/updatewcs/corrections.py b/stwcs/updatewcs/corrections.py
new file mode 100644
index 0000000..d3641eb
--- /dev/null
+++ b/stwcs/updatewcs/corrections.py
@@ -0,0 +1,326 @@
+from __future__ import division, print_function # confidence high
+
+import copy
+import datetime
+import logging, time
+import numpy as np
+from numpy import linalg
+from stsci.tools import fileutil
+
+from . import npol
+from . import makewcs
+from .utils import diff_angles
+
+logger=logging.getLogger('stwcs.updatewcs.corrections')
+
+MakeWCS = makewcs.MakeWCS
+NPOLCorr = npol.NPOLCorr
+
+class TDDCorr(object):
+ """
+ Apply time dependent distortion correction to distortion coefficients and basic
+ WCS keywords. This correction **must** be done before any other WCS correction.
+
+ Parameters
+ ----------
+ ext_wcs: HSTWCS object
+ An HSTWCS object to be modified
+ ref_wcs: HSTWCS object
+ A reference HSTWCS object
+
+ Notes
+ -----
+ Compute the ACS/WFC time dependent distortion terms as described
+ in [1]_ and apply the correction to the WCS of the observation.
+
+ The model coefficients are stored in the primary header of the IDCTAB.
+ :math:`D_{ref}` is the reference date. The computed corrections are saved
+ in the science extension header as TDDALPHA and TDDBETA keywords.
+
+ .. math:: TDDALPHA = A_{0} + {A_{1}*(obsdate - D_{ref})}
+
+ .. math:: TDDBETA = B_{0} + B_{1}*(obsdate - D_{ref})
+
+
+ The time dependent distortion affects the IDCTAB coefficients, and the
+ relative location of the two chips. Because the linear order IDCTAB
+ coefficients ar eused in the computatuion of the NPOL extensions,
+ the TDD correction affects all components of the distortion model.
+
+ Application of TDD to the IDCTAB polynomial coefficients:
+ The TDD model is computed in Jay's frame, while the IDCTAB coefficients
+ are in the HST V2/V3 frame. The coefficients are transformed to Jay's frame,
+ TDD is applied and they are transformed back to the V2/V3 frame. This
+ correction is performed in this class.
+
+ Application of TDD to the relative location of the two chips is
+ done in `makewcs`.
+
+ References
+ ----------
+ .. [1] Jay Anderson, "Variation of the Distortion Solution of the WFC", ACS ISR 2007-08.
+
+ """
+ def updateWCS(cls, ext_wcs, ref_wcs):
+ """
+ - Calculates alpha and beta for ACS/WFC data.
+ - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA
+ """
+ logger.info("\n\tStarting TDDCorr: %s" % time.asctime())
+ ext_wcs.idcmodel.ocx = copy.deepcopy(ext_wcs.idcmodel.cx)
+ ext_wcs.idcmodel.ocy = copy.deepcopy(ext_wcs.idcmodel.cy)
+
+ newkw = {'TDDALPHA': None, 'TDDBETA':None,
+ 'OCX10':ext_wcs.idcmodel.ocx[1,0],
+ 'OCX11':ext_wcs.idcmodel.ocx[1,1],
+ 'OCY10':ext_wcs.idcmodel.ocy[1,0],
+ 'OCY11':ext_wcs.idcmodel.ocy[1,1],
+ 'TDD_CTA':None, 'TDD_CTB':None,
+ 'TDD_CYA':None, 'TDD_CYB':None,
+ 'TDD_CXA':None, 'TDD_CXB':None}
+
+ if ext_wcs.idcmodel.refpix['skew_coeffs']is not None and \
+ ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'] is not None:
+ cls.apply_tdd2idc2015(ref_wcs)
+ cls.apply_tdd2idc2015(ext_wcs)
+
+ newkw.update({
+ 'TDD_CTA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTA'],
+ 'TDD_CYA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYA'],
+ 'TDD_CXA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXA'],
+ 'TDD_CTB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'],
+ 'TDD_CYB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYB'],
+ 'TDD_CXB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXB']})
+
+ elif ext_wcs.idcmodel.refpix['skew_coeffs']is not None and \
+ ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'] is not None:
+ logger.info("\n\t Applying 2014-calibrated TDD: %s" % time.asctime())
+ # We have 2014-calibrated TDD, not J.A.-style TDD
+ cls.apply_tdd2idc2(ref_wcs)
+ cls.apply_tdd2idc2(ext_wcs)
+ newkw.update({'TDD_CYA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_ALPHA'],
+ 'TDD_CYB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'],
+ 'TDD_CXA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_ALPHA'],
+ 'TDD_CXB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_BETA']})
+
+ else:
+ alpha, beta = cls.compute_alpha_beta(ext_wcs)
+ cls.apply_tdd2idc(ref_wcs, alpha, beta)
+ cls.apply_tdd2idc(ext_wcs, alpha, beta)
+ ext_wcs.idcmodel.refpix['TDDALPHA'] = alpha
+ ext_wcs.idcmodel.refpix['TDDBETA'] = beta
+ ref_wcs.idcmodel.refpix['TDDALPHA'] = alpha
+ ref_wcs.idcmodel.refpix['TDDBETA'] = beta
+
+ newkw.update( {'TDDALPHA': alpha, 'TDDBETA':beta} )
+
+ return newkw
+ updateWCS = classmethod(updateWCS)
+
+ def apply_tdd2idc2015(cls, hwcs):
+ """ Applies 2015-calibrated TDD correction to a couple of IDCTAB
+ coefficients for ACS/WFC observations.
+ """
+ if not isinstance(hwcs.date_obs,float):
+ year,month,day = hwcs.date_obs.split('-')
+ rdate = datetime.datetime(int(year),int(month),int(day))
+ rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+ else:
+ rday = hwcs.date_obs
+
+ skew_coeffs = hwcs.idcmodel.refpix['skew_coeffs']
+ delta_date = rday - skew_coeffs['TDD_DATE']
+
+ if skew_coeffs['TDD_CXB'] is not None:
+ hwcs.idcmodel.cx[1,1] += skew_coeffs['TDD_CXB']*delta_date
+ if skew_coeffs['TDD_CTB'] is not None:
+ hwcs.idcmodel.cy[1,1] += skew_coeffs['TDD_CTB']*delta_date
+ if skew_coeffs['TDD_CYB'] is not None:
+ hwcs.idcmodel.cy[1,0] += skew_coeffs['TDD_CYB']*delta_date
+ #print("CX[1,1]_TDD={}, CY[1,1]_TDD={}, CY[1,0]_TDD={}".format(hwcs.idcmodel.cx[1,1],hwcs.idcmodel.cy[1,1],hwcs.idcmodel.cy[1,0]))
+
+ apply_tdd2idc2015 = classmethod(apply_tdd2idc2015)
+
+ def apply_tdd2idc2(cls, hwcs):
+ """ Applies 2014-calibrated TDD correction to single IDCTAB coefficient
+ of an ACS/WFC observation.
+ """
+ if not isinstance(hwcs.date_obs,float):
+ year,month,day = hwcs.date_obs.split('-')
+ rdate = datetime.datetime(int(year),int(month),int(day))
+ rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+ else:
+ rday = hwcs.date_obs
+
+ skew_coeffs = hwcs.idcmodel.refpix['skew_coeffs']
+ cy_beta = skew_coeffs['TDD_CY_BETA']
+ cy_alpha = skew_coeffs['TDD_CY_ALPHA']
+ delta_date = rday - skew_coeffs['TDD_DATE']
+ print("DELTA_DATE: {0} based on rday: {1}, TDD_DATE: {2}".format(delta_date,rday,skew_coeffs['TDD_DATE']))
+
+ if cy_alpha is None:
+ hwcs.idcmodel.cy[1,1] += cy_beta*delta_date
+ else:
+ new_beta = cy_alpha + cy_beta*delta_date
+ hwcs.idcmodel.cy[1,1] = new_beta
+ print("CY11: {0} based on alpha: {1}, beta: {2}".format(hwcs.idcmodel.cy[1,1],cy_alpha,cy_beta))
+
+ cx_beta = skew_coeffs['TDD_CX_BETA']
+ cx_alpha = skew_coeffs['TDD_CX_ALPHA']
+ if cx_alpha is not None:
+ new_beta = cx_alpha + cx_beta*delta_date
+ hwcs.idcmodel.cx[1,1] = new_beta
+ print("CX11: {0} based on alpha: {1}, beta: {2}".format(new_beta,cx_alpha,cx_beta))
+
+ apply_tdd2idc2 = classmethod(apply_tdd2idc2)
+
+ def apply_tdd2idc(cls, hwcs, alpha, beta):
+ """
+ Applies TDD to the idctab coefficients of a ACS/WFC observation.
+ This should be always the first correction.
+ """
+ theta_v2v3 = 2.234529
+ mrotp = fileutil.buildRotMatrix(theta_v2v3)
+ mrotn = fileutil.buildRotMatrix(-theta_v2v3)
+ tdd_mat = np.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],np.float64)
+ abmat1 = np.dot(tdd_mat, mrotn)
+ abmat2 = np.dot(mrotp,abmat1)
+ xshape, yshape = hwcs.idcmodel.cx.shape, hwcs.idcmodel.cy.shape
+ icxy = np.dot(abmat2,[hwcs.idcmodel.cx.ravel(), hwcs.idcmodel.cy.ravel()])
+ hwcs.idcmodel.cx = icxy[0]
+ hwcs.idcmodel.cy = icxy[1]
+ hwcs.idcmodel.cx.shape = xshape
+ hwcs.idcmodel.cy.shape = yshape
+
+ apply_tdd2idc = classmethod(apply_tdd2idc)
+
+ def compute_alpha_beta(cls, ext_wcs):
+ """
+ Compute the ACS time dependent distortion skew terms
+ as described in ACS ISR 07-08 by J. Anderson.
+
+ Jay's code only computes the alpha/beta values based on a decimal year
+ with only 3 digits, so this line reproduces that when needed for comparison
+ with his results.
+ rday = float(('%0.3f')%rday)
+
+ The zero-point terms account for the skew accumulated between
+ 2002.0 and 2004.5, when the latest IDCTAB was delivered.
+ alpha = 0.095 + 0.090*(rday-dday)/2.5
+ beta = -0.029 - 0.030*(rday-dday)/2.5
+ """
+ if not isinstance(ext_wcs.date_obs,float):
+ year,month,day = ext_wcs.date_obs.split('-')
+ rdate = datetime.datetime(int(year),int(month),int(day))
+ rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+ else:
+ rday = ext_wcs.date_obs
+
+ skew_coeffs = ext_wcs.idcmodel.refpix['skew_coeffs']
+ if skew_coeffs is None:
+ # Only print out warning for post-SM4 data where this may matter
+ if rday > 2009.0:
+ err_str = "------------------------------------------------------------------------ \n"
+ err_str += "WARNING: the IDCTAB geometric distortion file specified in the image \n"
+ err_str += " header did not have the time-dependent distortion coefficients. \n"
+ err_str += " The pre-SM4 time-dependent skew solution will be used by default.\n"
+ err_str += " Please update IDCTAB with new reference file from HST archive. \n"
+ err_str += "------------------------------------------------------------------------ \n"
+ print(err_str)
+ # Using default pre-SM4 coefficients
+ skew_coeffs = {'TDD_A':[0.095,0.090/2.5],
+ 'TDD_B':[-0.029,-0.030/2.5],
+ 'TDD_DATE':2004.5,'TDDORDER':1}
+
+ alpha = 0
+ beta = 0
+ # Compute skew terms, allowing for non-linear coefficients as well
+ for c in range(skew_coeffs['TDDORDER']+1):
+ alpha += skew_coeffs['TDD_A'][c]* np.power((rday-skew_coeffs['TDD_DATE']),c)
+ beta += skew_coeffs['TDD_B'][c]*np.power((rday-skew_coeffs['TDD_DATE']),c)
+
+ return alpha,beta
+ compute_alpha_beta = classmethod(compute_alpha_beta)
+
+
+class VACorr(object):
+ """
+ Apply velocity aberation correction to WCS keywords.
+
+ Notes
+ -----
+ Velocity Aberration is stored in the extension header keyword 'VAFACTOR'.
+ The correction is applied to the CD matrix and CRVALs.
+
+ """
+ def updateWCS(cls, ext_wcs, ref_wcs):
+ logger.info("\n\tStarting VACorr: %s" % time.asctime())
+ if ext_wcs.vafactor != 1:
+ ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor
+ crval0 = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0],
+ ref_wcs.wcs.crval[0])
+ crval1 = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1],
+ ref_wcs.wcs.crval[1])
+ crval = np.array([crval0,crval1])
+ ext_wcs.wcs.crval = crval
+ ext_wcs.wcs.set()
+ else:
+ pass
+ kw2update={'CD1_1': ext_wcs.wcs.cd[0,0], 'CD1_2':ext_wcs.wcs.cd[0,1],
+ 'CD2_1':ext_wcs.wcs.cd[1,0], 'CD2_2':ext_wcs.wcs.cd[1,1],
+ 'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]}
+ return kw2update
+
+ updateWCS = classmethod(updateWCS)
+
+
+class CompSIP(object):
+ """
+ Compute Simple Imaging Polynomial (SIP) coefficients as defined in [2]_
+ from IDC table coefficients.
+
+ This class transforms the TDD corrected IDCTAB coefficients into SIP format.
+ It also applies a binning factor to the coefficients if the observation was
+ binned.
+
+ References
+ ----------
+ .. [2] David Shupe, et al, "The SIP Convention of representing Distortion
+ in FITS Image headers", Astronomical Data Analysis Software And Systems, ASP
+ Conference Series, Vol. 347, 2005
+
+ """
+ def updateWCS(cls, ext_wcs, ref_wcs):
+ logger.info("\n\tStarting CompSIP: %s" %time.asctime())
+ kw2update = {}
+ if not ext_wcs.idcmodel:
+ logger.info("IDC model not found, SIP coefficient will not be computed")
+ return kw2update
+ order = ext_wcs.idcmodel.norder
+ kw2update['A_ORDER'] = order
+ kw2update['B_ORDER'] = order
+ pscale = ext_wcs.idcmodel.refpix['PSCALE']
+
+ cx = ext_wcs.idcmodel.cx
+ cy = ext_wcs.idcmodel.cy
+
+ matr = np.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=np.float64)
+ imatr = linalg.inv(matr)
+ akeys1 = np.zeros((order+1,order+1), dtype=np.float64)
+ bkeys1 = np.zeros((order+1,order+1), dtype=np.float64)
+ for n in range(order+1):
+ for m in range(order+1):
+ if n >= m and n>=2:
+ idcval = np.array([[cx[n,m]],[cy[n,m]]])
+ sipval = np.dot(imatr, idcval)
+ akeys1[m,n-m] = sipval[0]
+ bkeys1[m,n-m] = sipval[1]
+ Akey="A_%d_%d" % (m,n-m)
+ Bkey="B_%d_%d" % (m,n-m)
+ kw2update[Akey] = sipval[0,0] * ext_wcs.binned
+ kw2update[Bkey] = sipval[1,0] * ext_wcs.binned
+ kw2update['CTYPE1'] = 'RA---TAN-SIP'
+ kw2update['CTYPE2'] = 'DEC--TAN-SIP'
+ return kw2update
+
+ updateWCS = classmethod(updateWCS)
diff --git a/stwcs/updatewcs/det2im.py b/stwcs/updatewcs/det2im.py
new file mode 100644
index 0000000..bddb37b
--- /dev/null
+++ b/stwcs/updatewcs/det2im.py
@@ -0,0 +1,299 @@
+from __future__ import absolute_import, division # confidence high
+
+import numpy as np
+from astropy.io import fits
+from stsci.tools import fileutil
+
+from . import utils
+
+import logging, time
+logger = logging.getLogger('stwcs.updatewcs.d2im')
+
+class DET2IMCorr(object):
+ """
+ Defines a Lookup table prior distortion correction as per WCS paper IV.
+ It uses a reference file defined by the D2IMFILE (suffix 'd2im') keyword
+ in the primary header.
+
+ Notes
+ -----
+ - Using extensions in the reference file create a WCSDVARR extensions
+ and add them to the science file.
+ - Add record-valued keywords to the science extension header to describe
+ the lookup tables.
+ - Add a keyword 'D2IMEXT' to the science extension header to store
+ the name of the reference file used to create the WCSDVARR extensions.
+
+ If WCSDVARR extensions exist and `D2IMFILE` is different from `D2IMEXT`,
+ a subsequent update will overwrite the existing extensions.
+ If WCSDVARR extensions were not found in the science file, they will be added.
+
+ """
+
+ def updateWCS(cls, fobj):
+ """
+ Parameters
+ ----------
+ fobj: `astropy.io.fits.HDUList` object
+ Science file, for which a distortion correction in a NPOLFILE is available
+
+ """
+ logger.info("\n\tStarting DET2IM: %s" %time.asctime())
+ try:
+ assert isinstance(fobj, fits.HDUList)
+ except AssertionError:
+ logger.exception('\n\tInput must be a fits.HDUList object')
+ raise
+
+ cls.applyDet2ImCorr(fobj)
+ d2imfile = fobj[0].header['D2IMFILE']
+
+ new_kw = {'D2IMEXT': d2imfile}
+ return new_kw
+
+ updateWCS = classmethod(updateWCS)
+
+ def applyDet2ImCorr(cls, fobj):
+ """
+ For each science extension in a fits file object:
+ - create a WCSDVARR extension
+ - update science header
+ - add/update D2IMEXT keyword
+ """
+ d2imfile = fileutil.osfn(fobj[0].header['D2IMFILE'])
+ # Map D2IMARR EXTVER numbers to FITS extension numbers
+ wcsdvarr_ind = cls.getWCSIndex(fobj)
+ d2im_num_ext = 1
+ for ext in fobj:
+ try:
+ extname = ext.header['EXTNAME'].lower()
+ except KeyError:
+ continue
+ if extname == 'sci':
+ extversion = ext.header['EXTVER']
+ ccdchip = cls.get_ccdchip(fobj, extname='SCI', extver=extversion)
+ header = ext.header
+ # get the data arrays from the reference file
+ dx, dy = cls.getData(d2imfile, ccdchip)
+ # Determine EXTVER for the D2IMARR extension from the D2I file (EXTNAME, EXTVER) kw.
+ # This is used to populate DPj.EXTVER kw
+ for ename in zip(['DX', 'DY'], [dx, dy]):
+ if ename[1] is not None:
+ error_val = ename[1].max()
+ cls.addSciExtKw(header, wdvarr_ver=d2im_num_ext, d2im_extname=ename[0], error_val=error_val)
+ hdu = cls.createD2ImHDU(header, d2imfile=d2imfile,
+ wdvarr_ver=d2im_num_ext,
+ d2im_extname=ename[0],
+ data=ename[1],ccdchip=ccdchip)
+ if wcsdvarr_ind and d2im_num_ext in wcsdvarr_ind:
+ fobj[wcsdvarr_ind[d2im_num_ext]] = hdu
+ else:
+ fobj.append(hdu)
+ d2im_num_ext = d2im_num_ext + 1
+ applyDet2ImCorr = classmethod(applyDet2ImCorr)
+
+ def getWCSIndex(cls, fobj):
+
+ """
+ If fobj has WCSDVARR extensions:
+ returns a mapping of their EXTVER kw to file object extension numbers
+ if fobj does not have WCSDVARR extensions:
+ an empty dictionary is returned
+ """
+ wcsd = {}
+ for e in range(len(fobj)):
+ try:
+ ename = fobj[e].header['EXTNAME']
+ except KeyError:
+ continue
+ if ename == 'D2IMARR':
+ wcsd[fobj[e].header['EXTVER']] = e
+ logger.debug("A map of D2IMARR extensions %s" % wcsd)
+ return wcsd
+
+ getWCSIndex = classmethod(getWCSIndex)
+
+ def addSciExtKw(cls, hdr, wdvarr_ver=None, d2im_extname=None, error_val=0.0):
+ """
+ Adds kw to sci extension to define WCSDVARR lookup table extensions
+
+ """
+ if d2im_extname =='DX':
+ j=1
+ else:
+ j=2
+
+ d2imerror = 'D2IMERR%s' %j
+ d2imdis = 'D2IMDIS%s' %j
+ d2imext = 'D2IM%s.' %j + 'EXTVER'
+ d2imnaxes = 'D2IM%s.' %j +'NAXES'
+ d2imaxis1 = 'D2IM%s.' %j+'AXIS.1'
+ d2imaxis2 = 'D2IM%s.' %j+'AXIS.2'
+ keys = [d2imerror, d2imdis, d2imext, d2imnaxes, d2imaxis1, d2imaxis2]
+ values = {d2imerror: error_val,
+ d2imdis: 'Lookup',
+ d2imext: wdvarr_ver,
+ d2imnaxes: 2,
+ d2imaxis1: 1,
+ d2imaxis2: 2}
+
+ comments = {d2imerror: 'Maximum error of NPOL correction for axis %s' % j,
+ d2imdis: 'Detector to image correction type',
+ d2imext: 'Version number of WCSDVARR extension containing d2im lookup table',
+ d2imnaxes: 'Number of independent variables in d2im function',
+ d2imaxis1: 'Axis number of the jth independent variable in a d2im function',
+ d2imaxis2: 'Axis number of the jth independent variable in a d2im function'
+ }
+ # Look for HISTORY keywords. If present, insert new keywords before them
+ before_key = 'HISTORY'
+ if before_key not in hdr:
+ before_key = None
+
+ for key in keys:
+ hdr.set(key, value=values[key], comment=comments[key], before=before_key)
+
+ addSciExtKw = classmethod(addSciExtKw)
+
+ def getData(cls,d2imfile, ccdchip):
+ """
+ Get the data arrays from the reference D2I files
+ Make sure 'CCDCHIP' in the npolfile matches "CCDCHIP' in the science file.
+ """
+ xdata, ydata = (None, None)
+ d2im = fits.open(d2imfile)
+ for ext in d2im:
+ d2imextname = ext.header.get('EXTNAME',"")
+ d2imccdchip = ext.header.get('CCDCHIP',1)
+ if d2imextname == 'DX' and d2imccdchip == ccdchip:
+ xdata = ext.data.copy()
+ continue
+ elif d2imextname == 'DY' and d2imccdchip == ccdchip:
+ ydata = ext.data.copy()
+ continue
+ else:
+ continue
+ d2im.close()
+ return xdata, ydata
+ getData = classmethod(getData)
+
+ def createD2ImHDU(cls, sciheader, d2imfile=None, wdvarr_ver=1,
+ d2im_extname=None,data = None, ccdchip=1):
+ """
+ Creates an HDU to be added to the file object.
+ """
+ hdr = cls.createD2ImHdr(sciheader, d2imfile=d2imfile,
+ wdvarr_ver=wdvarr_ver, d2im_extname=d2im_extname,
+ ccdchip=ccdchip)
+ hdu = fits.ImageHDU(header=hdr, data=data)
+ return hdu
+
+ createD2ImHDU = classmethod(createD2ImHDU)
+
+ def createD2ImHdr(cls, sciheader, d2imfile, wdvarr_ver, d2im_extname, ccdchip):
+ """
+ Creates a header for the D2IMARR extension based on the D2I reference file
+ and sci extension header. The goal is to always work in image coordinates
+ (also for subarrays and binned images). The WCS for the D2IMARR extension
+ is such that a full size d2im table is created and then shifted or scaled
+ if the science image is a subarray or binned image.
+ """
+ d2im = fits.open(d2imfile)
+ d2im_phdr = d2im[0].header
+ for ext in d2im:
+ try:
+ d2imextname = ext.header['EXTNAME']
+ d2imextver = ext.header['EXTVER']
+ except KeyError:
+ continue
+ d2imccdchip = cls.get_ccdchip(d2im, extname=d2imextname, extver=d2imextver)
+ if d2imextname == d2im_extname and d2imccdchip == ccdchip:
+ d2im_header = ext.header
+ break
+ else:
+ continue
+ d2im.close()
+
+ naxis = d2im[1].header['NAXIS']
+ ccdchip = d2imextname
+
+ kw = { 'NAXIS': 'Size of the axis',
+ 'CDELT': 'Coordinate increment along axis',
+ 'CRPIX': 'Coordinate system reference pixel',
+ 'CRVAL': 'Coordinate system value at reference pixel',
+ }
+
+ kw_comm1 = {}
+ kw_val1 = {}
+ for key in kw.keys():
+ for i in range(1, naxis+1):
+ si = str(i)
+ kw_comm1[key+si] = kw[key]
+
+ for i in range(1, naxis+1):
+ si = str(i)
+ kw_val1['NAXIS'+si] = d2im_header.get('NAXIS'+si)
+ kw_val1['CDELT'+si] = d2im_header.get('CDELT'+si, 1.0) * sciheader.get('LTM'+si+'_'+si, 1)
+ kw_val1['CRPIX'+si] = d2im_header.get('CRPIX'+si, 0.0)
+ kw_val1['CRVAL'+si] = (d2im_header.get('CRVAL'+si, 0.0) - \
+ sciheader.get('LTV'+str(i), 0))
+ kw_comm0 = {'XTENSION': 'Image extension',
+ 'BITPIX': 'IEEE floating point',
+ 'NAXIS': 'Number of axes',
+ 'EXTNAME': 'WCS distortion array',
+ 'EXTVER': 'Distortion array version number',
+ 'PCOUNT': 'Special data area of size 0',
+ 'GCOUNT': 'One data group',
+ }
+
+ kw_val0 = { 'XTENSION': 'IMAGE',
+ 'BITPIX': -32,
+ 'NAXIS': naxis,
+ 'EXTNAME': 'D2IMARR',
+ 'EXTVER': wdvarr_ver,
+ 'PCOUNT': 0,
+ 'GCOUNT': 1,
+ 'CCDCHIP': ccdchip,
+ }
+ cdl = []
+ for key in kw_comm0.keys():
+ cdl.append((key, kw_val0[key], kw_comm0[key]))
+ for key in kw_comm1.keys():
+ cdl.append((key, kw_val1[key], kw_comm1[key]))
+ # Now add keywords from NPOLFILE header to document source of calibration
+ # include all keywords after and including 'FILENAME' from header
+ start_indx = -1
+ end_indx = 0
+ for i, c in enumerate(d2im_phdr):
+ if c == 'FILENAME':
+ start_indx = i
+ if c == '': # remove blanks from end of header
+ end_indx = i+1
+ break
+ if start_indx >= 0:
+ for card in d2im_phdr.cards[start_indx:end_indx]:
+ cdl.append(card)
+
+ hdr = fits.Header(cards=cdl)
+
+ return hdr
+
+ createD2ImHdr = classmethod(createD2ImHdr)
+
+ def get_ccdchip(cls, fobj, extname, extver):
+ """
+ Given a science file or npol file determine CCDCHIP
+ """
+ ccdchip = 1
+ if fobj[0].header['INSTRUME'] == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC':
+ ccdchip = fobj[extname, extver].header['CCDCHIP']
+ elif fobj[0].header['INSTRUME'] == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS':
+ ccdchip = fobj[extname, extver].header['CCDCHIP']
+ elif fobj[0].header['INSTRUME'] == 'WFPC2':
+ ccdchip = fobj[extname, extver].header['DETECTOR']
+ elif fobj[0].header['INSTRUME'] == 'STIS':
+ ccdchip = fobj[extname, extver].header['DETECTOR']
+ elif fobj[0].header['INSTRUME'] == 'NICMOS':
+ ccdchip = fobj[extname, extver].header['CAMERA']
+ return ccdchip
+
+ get_ccdchip = classmethod(get_ccdchip)
diff --git a/stwcs/updatewcs/makewcs.py b/stwcs/updatewcs/makewcs.py
new file mode 100644
index 0000000..06c6f9c
--- /dev/null
+++ b/stwcs/updatewcs/makewcs.py
@@ -0,0 +1,273 @@
+from __future__ import absolute_import, division # confidence high
+
+import datetime
+
+import numpy as np
+import logging, time
+from math import sin, sqrt, pow, cos, asin, atan2,pi
+
+from stsci.tools import fileutil
+from . import utils
+
+logger = logging.getLogger(__name__)
+
+class MakeWCS(object):
+ """
+ Recompute basic WCS keywords based on PA_V3 and distortion model.
+
+ Notes
+ -----
+ - Compute the reference chip WCS:
+
+ -- CRVAL: transform the model XREF/YREF to the sky
+ -- PA_V3 is calculated at the target position and adjusted
+ for each chip orientation
+ -- CD: PA_V3 and the model scale are used to cnstruct a CD matrix
+
+ - Compute the second chip WCS:
+
+ -- CRVAL: - the distance between the zero points of the two
+ chip models on the sky
+ -- CD matrix: first order coefficients are added to the components
+ of this distance and transfered on the sky. The difference
+ between CRVAL and these vectors is the new CD matrix for each chip.
+ -- CRPIX: chip's model zero point in pixel space (XREF/YREF)
+
+ - Time dependent distortion correction is applied to both chips' V2/V3 values.
+
+ """
+ tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]}
+ def updateWCS(cls, ext_wcs, ref_wcs):
+ """
+ recomputes the basic WCS kw
+ """
+ logger.info("\n\tStarting MakeWCS: %s" % time.asctime())
+ if not ext_wcs.idcmodel:
+ logger.info("IDC model not found, turning off Makewcs")
+ return {}
+ ltvoff, offshift = cls.getOffsets(ext_wcs)
+
+ v23_corr = cls.zero_point_corr(ext_wcs)
+ rv23_corr = cls.zero_point_corr(ref_wcs)
+
+ cls.uprefwcs(ext_wcs, ref_wcs, rv23_corr, ltvoff, offshift)
+ cls.upextwcs(ext_wcs, ref_wcs, v23_corr, rv23_corr, ltvoff, offshift)
+
+ kw2update = {'CD1_1': ext_wcs.wcs.cd[0,0],
+ 'CD1_2': ext_wcs.wcs.cd[0,1],
+ 'CD2_1': ext_wcs.wcs.cd[1,0],
+ 'CD2_2': ext_wcs.wcs.cd[1,1],
+ 'CRVAL1': ext_wcs.wcs.crval[0],
+ 'CRVAL2': ext_wcs.wcs.crval[1],
+ 'CRPIX1': ext_wcs.wcs.crpix[0],
+ 'CRPIX2': ext_wcs.wcs.crpix[1],
+ 'IDCTAB': ext_wcs.idctab,
+ 'OCX10' : ext_wcs.idcmodel.cx[1,0],
+ 'OCX11' : ext_wcs.idcmodel.cx[1,1],
+ 'OCY10' : ext_wcs.idcmodel.cy[1,0],
+ 'OCY11' : ext_wcs.idcmodel.cy[1,1]
+ }
+ return kw2update
+
+ updateWCS = classmethod(updateWCS)
+
+ def upextwcs(cls, ext_wcs, ref_wcs, v23_corr, rv23_corr, loff, offsh):
+ """
+ updates an extension wcs
+ """
+ ltvoffx, ltvoffy = loff[0], loff[1]
+ offshiftx, offshifty = offsh[0], offsh[1]
+ ltv1 = ext_wcs.ltv1
+ ltv2 = ext_wcs.ltv2
+ if ltv1 != 0. or ltv2 != 0.:
+ offsetx = ext_wcs.wcs.crpix[0] - ltv1 - ext_wcs.idcmodel.refpix['XREF']
+ offsety = ext_wcs.wcs.crpix[1] - ltv2 - ext_wcs.idcmodel.refpix['YREF']
+ ext_wcs.idcmodel.shift(offsetx,offsety)
+
+ fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy
+
+ ext_wcs.setPscale()
+ tddscale = (ref_wcs.pscale/fx[1,1])
+ v2 = ext_wcs.idcmodel.refpix['V2REF'] + v23_corr[0,0] * tddscale
+ v3 = ext_wcs.idcmodel.refpix['V3REF'] - v23_corr[1,0] * tddscale
+ v2ref = ref_wcs.idcmodel.refpix['V2REF'] + rv23_corr[0,0] * tddscale
+ v3ref = ref_wcs.idcmodel.refpix['V3REF'] - rv23_corr[1,0] * tddscale
+
+ R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0
+ off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0)
+
+ if v3 == v3ref:
+ theta=0.0
+ else:
+ theta = atan2(ext_wcs.parity[0][0]*(v2-v2ref), ext_wcs.parity[1][1]*(v3-v3ref))
+
+ if ref_wcs.idcmodel.refpix['THETA']: theta += ref_wcs.idcmodel.refpix['THETA']*pi/180.0
+
+ dX=(off*sin(theta)) + offshiftx
+ dY=(off*cos(theta)) + offshifty
+
+ px = np.array([[dX,dY]])
+ newcrval = ref_wcs.wcs.p2s(px, 1)['world'][0]
+ newcrpix = np.array([ext_wcs.idcmodel.refpix['XREF'] + ltvoffx,
+ ext_wcs.idcmodel.refpix['YREF'] + ltvoffy])
+ ext_wcs.wcs.crval = newcrval
+ ext_wcs.wcs.crpix = newcrpix
+ ext_wcs.wcs.set()
+
+ # Create a small vector, in reference image pixel scale
+ delmat = np.array([[fx[1,1], fy[1,1]], \
+ [fx[1,0], fy[1,0]]]) / R_scale/3600.
+
+ # Account for subarray offset
+ # Angle of chip relative to chip
+ if ext_wcs.idcmodel.refpix['THETA']:
+ dtheta = ext_wcs.idcmodel.refpix['THETA'] - ref_wcs.idcmodel.refpix['THETA']
+ else:
+ dtheta = 0.0
+
+ rrmat = fileutil.buildRotMatrix(dtheta)
+ # Rotate the vectors
+ dxy = np.dot(delmat, rrmat)
+ wc = ref_wcs.wcs.p2s((px + dxy), 1)['world']
+
+ # Calculate the new CDs and convert to degrees
+ cd11 = utils.diff_angles(wc[0,0],newcrval[0])*cos(newcrval[1]*pi/180.0)
+ cd12 = utils.diff_angles(wc[1,0],newcrval[0])*cos(newcrval[1]*pi/180.0)
+ cd21 = utils.diff_angles(wc[0,1],newcrval[1])
+ cd22 = utils.diff_angles(wc[1,1],newcrval[1])
+ cd = np.array([[cd11, cd12], [cd21, cd22]])
+ ext_wcs.wcs.cd = cd
+ ext_wcs.wcs.set()
+
+ upextwcs = classmethod(upextwcs)
+
+ def uprefwcs(cls, ext_wcs, ref_wcs, rv23_corr_tdd, loff, offsh):
+ """
+ Updates the reference chip
+ """
+ ltvoffx, ltvoffy = loff[0], loff[1]
+ ltv1 = ref_wcs.ltv1
+ ltv2 = ref_wcs.ltv2
+ if ref_wcs.ltv1 != 0. or ref_wcs.ltv2 != 0.:
+ offsetx = ref_wcs.wcs.crpix[0] - ltv1 - ref_wcs.idcmodel.refpix['XREF']
+ offsety = ref_wcs.wcs.crpix[1] - ltv2 - ref_wcs.idcmodel.refpix['YREF']
+ ref_wcs.idcmodel.shift(offsetx,offsety)
+
+ rfx, rfy = ref_wcs.idcmodel.cx, ref_wcs.idcmodel.cy
+
+ offshift = offsh
+ dec = ref_wcs.wcs.crval[1]
+ tddscale = (ref_wcs.pscale/ref_wcs.idcmodel.cx[1,1])
+ rv23 = [ref_wcs.idcmodel.refpix['V2REF'] + (rv23_corr_tdd[0,0] *tddscale),
+ ref_wcs.idcmodel.refpix['V3REF'] - (rv23_corr_tdd[1,0] * tddscale)]
+ # Get an approximate reference position on the sky
+ rref = np.array([[ref_wcs.idcmodel.refpix['XREF']+ltvoffx ,
+ ref_wcs.idcmodel.refpix['YREF']+ltvoffy]])
+
+ crval = ref_wcs.wcs.p2s(rref, 1)['world'][0]
+ # Convert the PA_V3 orientation to the orientation at the aperture
+ # This is for the reference chip only - we use this for the
+ # reference tangent plane definition
+ # It has the same orientation as the reference chip
+ pv = troll(ext_wcs.pav3,dec,rv23[0], rv23[1])
+ # Add the chip rotation angle
+ if ref_wcs.idcmodel.refpix['THETA']:
+ pv += ref_wcs.idcmodel.refpix['THETA']
+
+
+ # Set values for the rest of the reference WCS
+ ref_wcs.wcs.crval = crval
+ ref_wcs.wcs.crpix = np.array([0.0,0.0])+offsh
+ parity = ref_wcs.parity
+ R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0
+ cd11 = parity[0][0] * cos(pv*pi/180.0)*R_scale
+ cd12 = parity[0][0] * -sin(pv*pi/180.0)*R_scale
+ cd21 = parity[1][1] * sin(pv*pi/180.0)*R_scale
+ cd22 = parity[1][1] * cos(pv*pi/180.0)*R_scale
+
+ rcd = np.array([[cd11, cd12], [cd21, cd22]])
+ ref_wcs.wcs.cd = rcd
+ ref_wcs.wcs.set()
+
+ uprefwcs = classmethod(uprefwcs)
+
+ def zero_point_corr(cls,hwcs):
+ if hwcs.idcmodel.refpix['skew_coeffs'] is not None and 'TDD_CY_BETA' in hwcs.idcmodel.refpix['skew_coeffs']:
+ v23_corr = np.array([[0.],[0.]])
+ return v23_corr
+ else:
+ try:
+ alpha = hwcs.idcmodel.refpix['TDDALPHA']
+ beta = hwcs.idcmodel.refpix['TDDBETA']
+ except KeyError:
+ alpha = 0.0
+ beta = 0.0
+ v23_corr = np.array([[0.],[0.]])
+ logger.debug("\n\tTDD Zero point correction for chip %s defaulted to: %s" % (hwcs.chip, v23_corr))
+ return v23_corr
+
+ tdd = np.array([[beta, alpha], [alpha, -beta]])
+ mrotp = fileutil.buildRotMatrix(2.234529)/2048.
+ xy0 = np.array([[cls.tdd_xyref[hwcs.chip][0]-2048.], [cls.tdd_xyref[hwcs.chip][1]-2048.]])
+ v23_corr = np.dot(mrotp,np.dot(tdd,xy0)) * 0.05
+ logger.debug("\n\tTDD Zero point correction for chip %s: %s" % (hwcs.chip, v23_corr))
+ return v23_corr
+
+ zero_point_corr = classmethod(zero_point_corr)
+
+ def getOffsets(cls, ext_wcs):
+ ltv1 = ext_wcs.ltv1
+ ltv2 = ext_wcs.ltv2
+
+ offsetx = ext_wcs.wcs.crpix[0] - ltv1 - ext_wcs.idcmodel.refpix['XREF']
+ offsety = ext_wcs.wcs.crpix[1] - ltv2 - ext_wcs.idcmodel.refpix['YREF']
+
+ shiftx = ext_wcs.idcmodel.refpix['XREF'] + ltv1
+ shifty = ext_wcs.idcmodel.refpix['YREF'] + ltv2
+ if ltv1 != 0. or ltv2 != 0.:
+ ltvoffx = ltv1 + offsetx
+ ltvoffy = ltv2 + offsety
+ offshiftx = offsetx + shiftx
+ offshifty = offsety + shifty
+ else:
+ ltvoffx = 0.
+ ltvoffy = 0.
+ offshiftx = 0.
+ offshifty = 0.
+
+ ltvoff = np.array([ltvoffx, ltvoffy])
+ offshift = np.array([offshiftx, offshifty])
+ return ltvoff, offshift
+
+ getOffsets = classmethod(getOffsets)
+
+
+def troll(roll, dec, v2, v3):
+ """ Computes the roll angle at the target position based on:
+ the roll angle at the V1 axis(roll),
+ the dec of the target(dec), and
+ the V2/V3 position of the aperture (v2,v3) in arcseconds.
+
+ Based on algorithm provided by Colin Cox and used in
+ Generic Conversion at STScI.
+ """
+ # Convert all angles to radians
+ _roll = np.deg2rad(roll)
+ _dec = np.deg2rad(dec)
+ _v2 = np.deg2rad(v2 / 3600.)
+ _v3 = np.deg2rad(v3 / 3600.)
+
+ # compute components
+ sin_rho = sqrt((pow(sin(_v2),2)+pow(sin(_v3),2)) - (pow(sin(_v2),2)*pow(sin(_v3),2)))
+ rho = asin(sin_rho)
+ beta = asin(sin(_v3)/sin_rho)
+ if _v2 < 0: beta = pi - beta
+ gamma = asin(sin(_v2)/sin_rho)
+ if _v3 < 0: gamma = pi - gamma
+ A = pi/2. + _roll - beta
+ B = atan2( sin(A)*cos(_dec), (sin(_dec)*sin_rho - cos(_dec)*cos(rho)*cos(A)))
+
+ # compute final value
+ troll = np.rad2deg(pi - (gamma+B))
+
+ return troll
diff --git a/stwcs/updatewcs/npol.py b/stwcs/updatewcs/npol.py
new file mode 100644
index 0000000..33579f0
--- /dev/null
+++ b/stwcs/updatewcs/npol.py
@@ -0,0 +1,343 @@
+from __future__ import absolute_import, division # confidence high
+
+import logging, time
+
+import numpy as np
+from astropy.io import fits
+
+from stsci.tools import fileutil
+from . import utils
+
+logger = logging.getLogger('stwcs.updatewcs.npol')
+
+class NPOLCorr(object):
+ """
+ Defines a Lookup table prior distortion correction as per WCS paper IV.
+ It uses a reference file defined by the NPOLFILE (suffix 'NPL') keyword
+ in the primary header.
+
+ Notes
+ -----
+ - Using extensions in the reference file create a WCSDVARR extensions
+ and add them to the science file.
+ - Add record-valued keywords to the science extension header to describe
+ the lookup tables.
+ - Add a keyword 'NPOLEXT' to the science extension header to store
+ the name of the reference file used to create the WCSDVARR extensions.
+
+ If WCSDVARR extensions exist and `NPOLFILE` is different from `NPOLEXT`,
+ a subsequent update will overwrite the existing extensions.
+ If WCSDVARR extensions were not found in the science file, they will be added.
+
+ It is assumed that the NPL reference files were created to work with IDC tables
+ but will be applied with SIP coefficients. A transformation is applied to correct
+ for the fact that the lookup tables will be applied before the first order coefficients
+ which are in the CD matrix when the SIP convention is used.
+ """
+
+ def updateWCS(cls, fobj):
+ """
+ Parameters
+ ----------
+ fobj : `astropy.io.fits.HDUList` object
+ Science file, for which a distortion correction in a NPOLFILE is available
+
+ """
+ logger.info("\n\tStarting NPOL: %s" %time.asctime())
+ try:
+ assert isinstance(fobj, fits.HDUList)
+ except AssertionError:
+ logger.exception('\n\tInput must be a fits.HDUList object')
+ raise
+
+ cls.applyNPOLCorr(fobj)
+ nplfile = fobj[0].header['NPOLFILE']
+
+ new_kw = {'NPOLEXT': nplfile}
+ return new_kw
+
+ updateWCS = classmethod(updateWCS)
+
+ def applyNPOLCorr(cls, fobj):
+ """
+ For each science extension in a fits file object:
+ - create a WCSDVARR extension
+ - update science header
+ - add/update NPOLEXT keyword
+ """
+ nplfile = fileutil.osfn(fobj[0].header['NPOLFILE'])
+ # Map WCSDVARR EXTVER numbers to extension numbers
+ wcsdvarr_ind = cls.getWCSIndex(fobj)
+ for ext in fobj:
+ try:
+ extname = ext.header['EXTNAME'].lower()
+ except KeyError:
+ continue
+ if extname == 'sci':
+ extversion = ext.header['EXTVER']
+ ccdchip = cls.get_ccdchip(fobj, extname='SCI', extver=extversion)
+ header = ext.header
+ # get the data arrays from the reference file and transform
+ # them for use with SIP
+ dx,dy = cls.getData(nplfile, ccdchip)
+ idccoeffs = cls.getIDCCoeffs(header)
+
+ if idccoeffs is not None:
+ dx, dy = cls.transformData(dx,dy, idccoeffs)
+
+ # Determine EXTVER for the WCSDVARR extension from the
+ # NPL file (EXTNAME, EXTVER) kw.
+ # This is used to populate DPj.EXTVER kw
+ wcsdvarr_x_version = 2 * extversion -1
+ wcsdvarr_y_version = 2 * extversion
+ for ename in zip(['DX', 'DY'], [wcsdvarr_x_version,wcsdvarr_y_version],[dx, dy]):
+ error_val = ename[2].max()
+ cls.addSciExtKw(header, wdvarr_ver=ename[1], npol_extname=ename[0], error_val=error_val)
+ hdu = cls.createNpolHDU(header, npolfile=nplfile, \
+ wdvarr_ver=ename[1], npl_extname=ename[0], data=ename[2],ccdchip=ccdchip)
+ if wcsdvarr_ind:
+ fobj[wcsdvarr_ind[ename[1]]] = hdu
+ else:
+ fobj.append(hdu)
+
+
+ applyNPOLCorr = classmethod(applyNPOLCorr)
+
+ def getWCSIndex(cls, fobj):
+
+ """
+ If fobj has WCSDVARR extensions:
+ returns a mapping of their EXTVER kw to file object extension numbers
+ if fobj does not have WCSDVARR extensions:
+ an empty dictionary is returned
+ """
+ wcsd = {}
+ for e in range(len(fobj)):
+ try:
+ ename = fobj[e].header['EXTNAME']
+ except KeyError:
+ continue
+ if ename == 'WCSDVARR':
+ wcsd[fobj[e].header['EXTVER']] = e
+ logger.debug("A map of WSCDVARR externsions %s" % wcsd)
+ return wcsd
+
+ getWCSIndex = classmethod(getWCSIndex)
+
+ def addSciExtKw(cls, hdr, wdvarr_ver=None, npol_extname=None, error_val=0.0):
+ """
+ Adds kw to sci extension to define WCSDVARR lookup table extensions
+
+ """
+ if npol_extname =='DX':
+ j=1
+ else:
+ j=2
+
+ cperror = 'CPERR%s' %j
+ cpdis = 'CPDIS%s' %j
+ dpext = 'DP%s.' %j + 'EXTVER'
+ dpnaxes = 'DP%s.' %j +'NAXES'
+ dpaxis1 = 'DP%s.' %j+'AXIS.1'
+ dpaxis2 = 'DP%s.' %j+'AXIS.2'
+ keys = [cperror, cpdis, dpext, dpnaxes, dpaxis1, dpaxis2]
+ values = {cperror: error_val,
+ cpdis: 'Lookup',
+ dpext: wdvarr_ver,
+ dpnaxes: 2,
+ dpaxis1: 1,
+ dpaxis2: 2}
+
+ comments = {cperror: 'Maximum error of NPOL correction for axis %s' % j,
+ cpdis: 'Prior distortion function type',
+ dpext: 'Version number of WCSDVARR extension containing lookup distortion table',
+ dpnaxes: 'Number of independent variables in distortion function',
+ dpaxis1: 'Axis number of the jth independent variable in a distortion function',
+ dpaxis2: 'Axis number of the jth independent variable in a distortion function'
+ }
+ # Look for HISTORY keywords. If present, insert new keywords before them
+ before_key = 'HISTORY'
+ if before_key not in hdr:
+ before_key = None
+
+ for key in keys:
+ hdr.set(key, value=values[key], comment=comments[key], before=before_key)
+
+ addSciExtKw = classmethod(addSciExtKw)
+
+ def getData(cls,nplfile, ccdchip):
+ """
+ Get the data arrays from the reference NPOL files
+ Make sure 'CCDCHIP' in the npolfile matches "CCDCHIP' in the science file.
+ """
+ npl = fits.open(nplfile)
+ for ext in npl:
+ nplextname = ext.header.get('EXTNAME',"")
+ nplccdchip = ext.header.get('CCDCHIP',1)
+ if nplextname == 'DX' and nplccdchip == ccdchip:
+ xdata = ext.data.copy()
+ continue
+ elif nplextname == 'DY' and nplccdchip == ccdchip:
+ ydata = ext.data.copy()
+ continue
+ else:
+ continue
+ npl.close()
+ return xdata, ydata
+ getData = classmethod(getData)
+
+ def transformData(cls, dx, dy, coeffs):
+ """
+ Transform the NPOL data arrays for use with SIP
+ """
+ ndx, ndy = np.dot(coeffs, [dx.ravel(), dy.ravel()]).astype(np.float32)
+ ndx.shape = dx.shape
+ ndy.shape=dy.shape
+ return ndx, ndy
+
+ transformData = classmethod(transformData)
+
+ def getIDCCoeffs(cls, header):
+ """
+ Return a matrix of the scaled first order IDC coefficients.
+ """
+ try:
+ ocx10 = header['OCX10']
+ ocx11 = header['OCX11']
+ ocy10 = header['OCY10']
+ ocy11 = header['OCY11']
+ coeffs = np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32)
+ except KeyError:
+ logger.exception('\n\tFirst order IDCTAB coefficients are not available. \n\
+ Cannot convert SIP to IDC coefficients.')
+ return None
+ try:
+ idcscale = header['IDCSCALE']
+ except KeyError:
+ logger.exception("IDCSCALE not found in header - setting it to 1.")
+ idcscale = 1
+
+ return np.linalg.inv(coeffs/idcscale)
+
+ getIDCCoeffs = classmethod(getIDCCoeffs)
+
+ def createNpolHDU(cls, sciheader, npolfile=None, wdvarr_ver=1, npl_extname=None,data = None, ccdchip=1):
+ """
+ Creates an HDU to be added to the file object.
+ """
+ hdr = cls.createNpolHdr(sciheader, npolfile=npolfile, wdvarr_ver=wdvarr_ver, npl_extname=npl_extname, ccdchip=ccdchip)
+ hdu = fits.ImageHDU(header=hdr, data=data)
+ return hdu
+
+ createNpolHDU = classmethod(createNpolHDU)
+
+ def createNpolHdr(cls, sciheader, npolfile, wdvarr_ver, npl_extname, ccdchip):
+ """
+ Creates a header for the WCSDVARR extension based on the NPOL reference file
+ and sci extension header. The goal is to always work in image coordinates
+ (also for subarrays and binned images. The WCS for the WCSDVARR extension
+ i ssuch that a full size npol table is created and then shifted or scaled
+ if the science image is a subarray or binned image.
+ """
+ npl = fits.open(npolfile)
+ npol_phdr = npl[0].header
+ for ext in npl:
+ try:
+ nplextname = ext.header['EXTNAME']
+ nplextver = ext.header['EXTVER']
+ except KeyError:
+ continue
+ nplccdchip = cls.get_ccdchip(npl, extname=nplextname, extver=nplextver)
+ if nplextname == npl_extname and nplccdchip == ccdchip:
+ npol_header = ext.header
+ break
+ else:
+ continue
+ npl.close()
+
+ naxis = npl[1].header['NAXIS']
+ ccdchip = nplextname #npol_header['CCDCHIP']
+
+ kw = { 'NAXIS': 'Size of the axis',
+ 'CDELT': 'Coordinate increment along axis',
+ 'CRPIX': 'Coordinate system reference pixel',
+ 'CRVAL': 'Coordinate system value at reference pixel',
+ }
+
+ kw_comm1 = {}
+ kw_val1 = {}
+ for key in kw.keys():
+ for i in range(1, naxis+1):
+ si = str(i)
+ kw_comm1[key+si] = kw[key]
+
+ for i in range(1, naxis+1):
+ si = str(i)
+ kw_val1['NAXIS'+si] = npol_header.get('NAXIS'+si)
+ kw_val1['CDELT'+si] = npol_header.get('CDELT'+si, 1.0) * \
+ sciheader.get('LTM'+si+'_'+si, 1)
+ kw_val1['CRPIX'+si] = npol_header.get('CRPIX'+si, 0.0)
+ kw_val1['CRVAL'+si] = (npol_header.get('CRVAL'+si, 0.0) - \
+ sciheader.get('LTV'+str(i), 0))
+
+ kw_comm0 = {'XTENSION': 'Image extension',
+ 'BITPIX': 'IEEE floating point',
+ 'NAXIS': 'Number of axes',
+ 'EXTNAME': 'WCS distortion array',
+ 'EXTVER': 'Distortion array version number',
+ 'PCOUNT': 'Special data area of size 0',
+ 'GCOUNT': 'One data group',
+ }
+
+ kw_val0 = { 'XTENSION': 'IMAGE',
+ 'BITPIX': -32,
+ 'NAXIS': naxis,
+ 'EXTNAME': 'WCSDVARR',
+ 'EXTVER': wdvarr_ver,
+ 'PCOUNT': 0,
+ 'GCOUNT': 1,
+ 'CCDCHIP': ccdchip,
+ }
+ cdl = []
+ for key in kw_comm0.keys():
+ cdl.append((key, kw_val0[key], kw_comm0[key]))
+ for key in kw_comm1.keys():
+ cdl.append((key, kw_val1[key], kw_comm1[key]))
+ # Now add keywords from NPOLFILE header to document source of calibration
+ # include all keywords after and including 'FILENAME' from header
+ start_indx = -1
+ end_indx = 0
+ for i, c in enumerate(npol_phdr):
+ if c == 'FILENAME':
+ start_indx = i
+ if c == '': # remove blanks from end of header
+ end_indx = i+1
+ break
+ if start_indx >= 0:
+ for card in npol_phdr.cards[start_indx:end_indx]:
+ cdl.append(card)
+
+ hdr = fits.Header(cards=cdl)
+
+ return hdr
+
+ createNpolHdr = classmethod(createNpolHdr)
+
+ def get_ccdchip(cls, fobj, extname, extver):
+ """
+ Given a science file or npol file determine CCDCHIP
+ """
+ ccdchip = 1
+ if fobj[0].header['INSTRUME'] == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC':
+ ccdchip = fobj[extname, extver].header['CCDCHIP']
+ elif fobj[0].header['INSTRUME'] == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS':
+ ccdchip = fobj[extname, extver].header['CCDCHIP']
+ elif fobj[0].header['INSTRUME'] == 'WFPC2':
+ ccdchip = fobj[extname, extver].header['DETECTOR']
+ elif fobj[0].header['INSTRUME'] == 'STIS':
+ ccdchip = fobj[extname, extver].header['DETECTOR']
+ elif fobj[0].header['INSTRUME'] == 'NICMOS':
+ ccdchip = fobj[extname, extver].header['CAMERA']
+ return ccdchip
+
+ get_ccdchip = classmethod(get_ccdchip)
diff --git a/stwcs/updatewcs/utils.py b/stwcs/updatewcs/utils.py
new file mode 100644
index 0000000..ebfc701
--- /dev/null
+++ b/stwcs/updatewcs/utils.py
@@ -0,0 +1,264 @@
+from __future__ import division # confidence high
+import os
+from astropy.io import fits
+from stsci.tools import fileutil
+
+import logging
+logger = logging.getLogger("stwcs.updatewcs.utils")
+
+
+def diff_angles(a,b):
+ """
+ Perform angle subtraction a-b taking into account
+ small-angle differences across 360degree line.
+ """
+
+ diff = a - b
+
+ if diff > 180.0:
+ diff -= 360.0
+
+ if diff < -180.0:
+ diff += 360.0
+
+ return diff
+
+def getBinning(fobj, extver=1):
+ # Return the binning factor
+ binned = 1
+ if fobj[0].header['INSTRUME'] == 'WFPC2':
+ mode = fobj[0].header.get('MODE', "")
+ if mode == 'AREA': binned = 2
+ else:
+ binned = fobj['SCI', extver].header.get('BINAXIS',1)
+ return binned
+
+def updateNEXTENDKw(fobj):
+ """
+ Updates PRIMARY header with correct value for NEXTEND, if present.
+
+ Parameters
+ -----------
+ fobj : `astropy.io.fits.HDUList`
+ The FITS object for file opened in `update` mode
+
+ """
+ if 'nextend' in fobj[0].header:
+ fobj[0].header['nextend'] = len(fobj) -1
+
+def extract_rootname(kwvalue,suffix=""):
+ """ Returns the rootname from a full reference filename
+
+ If a non-valid value (any of ['','N/A','NONE','INDEF',None]) is input,
+ simply return a string value of 'NONE'
+
+ This function will also replace any 'suffix' specified with a blank.
+ """
+ # check to see whether a valid kwvalue has been provided as input
+ if kwvalue.strip() in ['','N/A','NONE','INDEF',None]:
+ return 'NONE' # no valid value, so return 'NONE'
+
+ # for a valid kwvalue, parse out the rootname
+ # strip off any environment variable from input filename, if any are given
+ if '$' in kwvalue:
+ fullval = kwvalue[kwvalue.find('$')+1:]
+ else:
+ fullval = kwvalue
+ # Extract filename without path from kwvalue
+ fname = os.path.basename(fullval).strip()
+
+ # Now, rip out just the rootname from the full filename
+ rootname = fileutil.buildNewRootname(fname)
+
+ # Now, remove any known suffix from rootname
+ rootname = rootname.replace(suffix,'')
+ return rootname.strip()
+
+def construct_distname(fobj,wcsobj):
+ """
+ This function constructs the value for the keyword 'DISTNAME'.
+ It relies on the reference files specified by the keywords 'IDCTAB',
+ 'NPOLFILE', and 'D2IMFILE'.
+
+ The final constructed value will be of the form:
+ <idctab rootname>-<npolfile rootname>-<d2imfile rootname>
+ and have a value of 'NONE' if no reference files are specified.
+ """
+ idcname = extract_rootname(fobj[0].header.get('IDCTAB', "NONE"),
+ suffix='_idc')
+ if (idcname is None or idcname=='NONE') and wcsobj.sip is not None:
+ idcname = 'UNKNOWN'
+
+ npolname, npolfile = build_npolname(fobj)
+ if npolname is None and wcsobj.cpdis1 is not None:
+ npolname = 'UNKNOWN'
+
+ d2imname, d2imfile = build_d2imname(fobj)
+ if d2imname is None and wcsobj.det2im is not None:
+ d2imname = 'UNKNOWN'
+
+ sipname, idctab = build_sipname(fobj)
+
+ distname = build_distname(sipname,npolname,d2imname)
+ return {'DISTNAME':distname,'SIPNAME':sipname}
+
+def build_distname(sipname,npolname,d2imname):
+ """
+ Core function to build DISTNAME keyword value without the HSTWCS input.
+ """
+
+ distname = sipname.strip()
+ if npolname != 'NONE' or d2imname != 'NONE':
+ if d2imname != 'NONE':
+ distname+= '-'+npolname.strip() + '-'+d2imname.strip()
+ else:
+ distname+='-'+npolname.strip()
+
+ return distname
+
+def build_default_wcsname(idctab):
+
+ idcname = extract_rootname(idctab,suffix='_idc')
+ wcsname = 'IDC_' + idcname
+ return wcsname
+
+def build_sipname(fobj, fname=None, sipname=None):
+ """
+ Build a SIPNAME from IDCTAB
+
+ Parameters
+ ----------
+ fobj : `astropy.io.fits.HDUList`
+ file object
+ fname : string
+ science file name (to be used if ROOTNAME is not present
+ sipname : string
+ user supplied SIPNAME keyword
+
+ Returns
+ -------
+ sipname, idctab
+ """
+ try:
+ idctab = fobj[0].header['IDCTAB']
+ idcname = extract_rootname(idctab,suffix='_idc')
+ except KeyError:
+ idctab = 'N/A'
+ idcname= 'N/A'
+ if not fname:
+ try:
+ fname = fobj.filename()
+ except:
+ fname = " "
+ if not sipname:
+ try:
+ sipname = fobj[0].header["SIPNAME"]
+ if idcname == 'N/A' or idcname not in sipname:
+ raise KeyError
+ except KeyError:
+ if idcname != 'N/A':
+ try:
+ rootname = fobj[0].header['rootname']
+ except KeyError:
+ rootname = fname
+ sipname = rootname +'_'+ idcname
+ else:
+ if 'A_ORDER' in fobj[1].header or 'B_ORDER' in fobj[1].header:
+ sipname = 'UNKNOWN'
+ else:
+ idcname = 'NOMODEL'
+
+ return sipname, idctab
+
+def build_npolname(fobj, npolfile=None):
+ """
+ Build a NPOLNAME from NPOLFILE
+
+ Parameters
+ ----------
+ fobj : `astropy.io.fits.HDUList`
+ file object
+ npolfile : string
+ user supplied NPOLFILE keyword
+
+ Returns
+ -------
+ npolname, npolfile
+ """
+ if not npolfile:
+ try:
+ npolfile = fobj[0].header["NPOLFILE"]
+ except KeyError:
+ npolfile = ' '
+ if fileutil.countExtn(fobj, 'WCSDVARR'):
+ npolname = 'UNKNOWN'
+ else:
+ npolname = 'NOMODEL'
+ npolname = extract_rootname(npolfile, suffix='_npl')
+ if npolname == 'NONE':
+ npolname = 'NOMODEL'
+ else:
+ npolname = extract_rootname(npolfile, suffix='_npl')
+ if npolname == 'NONE':
+ npolname = 'NOMODEL'
+ return npolname, npolfile
+
+def build_d2imname(fobj, d2imfile=None):
+ """
+ Build a D2IMNAME from D2IMFILE
+
+ Parameters
+ ----------
+ fobj : `astropy.io.fits.HDUList`
+ file object
+ d2imfile : string
+ user supplied NPOLFILE keyword
+
+ Returns
+ -------
+ d2imname, d2imfile
+ """
+ if not d2imfile:
+ try:
+ d2imfile = fobj[0].header["D2IMFILE"]
+ except KeyError:
+ d2imfile = 'N/A'
+ if fileutil.countExtn(fobj, 'D2IMARR'):
+ d2imname = 'UNKNOWN'
+ else:
+ d2imname = 'NOMODEL'
+ d2imname = extract_rootname(d2imfile,suffix='_d2i')
+ if d2imname == 'NONE':
+ d2imname = 'NOMODEL'
+ else:
+ d2imname = extract_rootname(d2imfile,suffix='_d2i')
+ if d2imname == 'NONE':
+ d2imname = 'NOMODEL'
+ return d2imname, d2imfile
+
+
+def remove_distortion(fname, dist_keyword):
+ logger.info("Removing distortion {0} from file {0}".format(dist_keyword, fname))
+ from ..wcsutil import altwcs
+ if dist_keyword == "NPOLFILE":
+ extname = "WCSDVARR"
+ keywords = ["CPERR*", "DP1.*", "DP2.*", "CPDIS*", "NPOLEXT"]
+ elif dist_keyword == "D2IMFILE":
+ extname = "D2IMARR"
+ keywords = ["D2IMERR*", "D2IM1.*", "D2IM2.*", "D2IMDIS*", "D2IMEXT"]
+ else:
+ raise AttributeError("Unrecognized distortion keyword "
+ "{0} when attempting to remove distortion".format(dist_keyword))
+ ext_mapping = altwcs.mapFitsExt2HDUListInd(fname, "SCI").values()
+ f = fits.open(fname, mode="update")
+ for hdu in ext_mapping:
+ for kw in keywords:
+ try:
+ del f[hdu].header[kw]
+ except KeyError:
+ pass
+ ext_mapping = list(altwcs.mapFitsExt2HDUListInd(fname, extname).values())
+ ext_mapping.sort()
+ for hdu in ext_mapping[::-1]:
+ del f[hdu]
+ f.close()
diff --git a/stwcs/updatewcs/wfpc2_dgeo.py b/stwcs/updatewcs/wfpc2_dgeo.py
new file mode 100644
index 0000000..e57bb5c
--- /dev/null
+++ b/stwcs/updatewcs/wfpc2_dgeo.py
@@ -0,0 +1,123 @@
+""" wfpc2_dgeo - Functions to convert WFPC2 DGEOFILE into D2IMFILE
+
+"""
+import os
+import datetime
+
+import astropy
+from astropy.io import fits
+import numpy as np
+
+from stsci.tools import fileutil
+
+import logging
+logger = logging.getLogger("stwcs.updatewcs.apply_corrections")
+
+def update_wfpc2_d2geofile(filename, fhdu=None):
+ """
+ Creates a D2IMFILE from the DGEOFILE for a WFPC2 image (input), and
+ modifies the header to reflect the new usage.
+
+ Parameters
+ ----------
+ filename: string
+ Name of WFPC2 file to be processed. This file will be updated
+ to delete any reference to a DGEOFILE and add a D2IMFILE to replace
+ that correction when running updatewcs.
+ fhdu: object
+ FITS object for WFPC2 image. If user has already opened the WFPC2
+ file, they can simply pass that FITS object in for direct processing.
+
+ Returns
+ -------
+ d2imfile: string
+ Name of D2IMFILE created from DGEOFILE. The D2IMFILE keyword in the
+ image header will be updated/added to point to this newly created file.
+
+ """
+
+ close_fhdu = False
+ if fhdu is None:
+ fhdu = fileutil.openImage(filename, mode='update')
+ close_fhdu = True
+
+ dgeofile = fhdu['PRIMARY'].header.get('DGEOFILE', None)
+ already_converted = dgeofile not in [None, "N/A", "", " "]
+ if already_converted or 'ODGEOFIL' in fhdu['PRIMARY'].header:
+ if not already_converted:
+ dgeofile = fhdu['PRIMARY'].header.get('ODGEOFIL', None)
+ logger.info('Converting DGEOFILE %s into D2IMFILE...' % dgeofile)
+ rootname = filename[:filename.find('.fits')]
+ d2imfile = convert_dgeo_to_d2im(dgeofile,rootname)
+ fhdu['PRIMARY'].header['ODGEOFIL'] = dgeofile
+ fhdu['PRIMARY'].header['DGEOFILE'] = 'N/A'
+ fhdu['PRIMARY'].header['D2IMFILE'] = d2imfile
+ else:
+ d2imfile = None
+ fhdu['PRIMARY'].header['DGEOFILE'] = 'N/A'
+ if 'D2IMFILE' not in fhdu['PRIMARY'].header:
+ fhdu['PRIMARY'].header['D2IMFILE'] = 'N/A'
+
+ # Only close the file handle if opened in this function
+ if close_fhdu:
+ fhdu.close()
+
+ # return the d2imfile name so that calling routine can keep
+ # track of the new file created and delete it later if necessary
+ # (multidrizzle clean=True mode of operation)
+ return d2imfile
+
+def convert_dgeo_to_d2im(dgeofile,output,clobber=True):
+ """ Routine that converts the WFPC2 DGEOFILE into a D2IMFILE.
+ """
+ dgeo = fileutil.openImage(dgeofile)
+ outname = output+'_d2im.fits'
+
+ removeFileSafely(outname)
+ data = np.array([dgeo['dy',1].data[:,0]])
+ scihdu = fits.ImageHDU(data=data)
+ dgeo.close()
+ # add required keywords for D2IM header
+ scihdu.header['EXTNAME'] = ('DY', 'Extension name')
+ scihdu.header['EXTVER'] = (1, 'Extension version')
+ fits_str = 'PYFITS Version '+str(astropy.__version__)
+ scihdu.header['ORIGIN'] = (fits_str, 'FITS file originator')
+ scihdu.header['INHERIT'] = (False, 'Inherits global header')
+
+ dnow = datetime.datetime.now()
+ scihdu.header['DATE'] = (str(dnow).replace(' ','T'),
+ 'Date FITS file was generated')
+
+ scihdu.header['CRPIX1'] = (0, 'Distortion array reference pixel')
+ scihdu.header['CDELT1'] = (1, 'Grid step size in first coordinate')
+ scihdu.header['CRVAL1'] = (0, 'Image array pixel coordinate')
+ scihdu.header['CRPIX2'] = (0, 'Distortion array reference pixel')
+ scihdu.header['CDELT2'] = (1, 'Grid step size in second coordinate')
+ scihdu.header['CRVAL2'] = (0, 'Image array pixel coordinate')
+
+ phdu = fits.PrimaryHDU()
+ phdu.header['INSTRUME'] = 'WFPC2'
+ d2imhdu = fits.HDUList()
+ d2imhdu.append(phdu)
+ scihdu.header['DETECTOR'] = (1, 'CCD number of the detector: PC 1, WFC 2-4 ')
+ d2imhdu.append(scihdu.copy())
+ scihdu.header['EXTVER'] = (2, 'Extension version')
+ scihdu.header['DETECTOR'] = (2, 'CCD number of the detector: PC 1, WFC 2-4 ')
+ d2imhdu.append(scihdu.copy())
+ scihdu.header['EXTVER'] = (3, 'Extension version')
+ scihdu.header['DETECTOR'] = (3, 'CCD number of the detector: PC 1, WFC 2-4 ')
+ d2imhdu.append(scihdu.copy())
+ scihdu.header['EXTVER'] = (4, 'Extension version')
+ scihdu.header['DETECTOR'] = (4, 'CCD number of the detector: PC 1, WFC 2-4 ')
+ d2imhdu.append(scihdu.copy())
+ d2imhdu.writeto(outname)
+ d2imhdu.close()
+
+ return outname
+
+
+def removeFileSafely(filename,clobber=True):
+ """ Delete the file specified, but only if it exists and clobber is True.
+ """
+ if filename is not None and filename.strip() != '':
+ if os.path.exists(filename) and clobber: os.remove(filename)
generated by cgit (git 2.34.1) at 2025-05-02 11:37:33 -0400