From 749102aed2e2c258f00aeef6c78fadd704c14fe7 Mon Sep 17 00:00:00 2001 From: dencheva Date: Wed, 29 Apr 2009 12:13:55 +0000 Subject: Renamed hstwcs to udpatewcs git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/development/trunk/stwcs@7976 fe389314-cf27-0410-b35b-8c050e845b92 --- hstwcs/__init__.py | 196 -------------------------------- hstwcs/alignwcs.py | 95 ---------------- hstwcs/apply_corrections.py | 181 ----------------------------- hstwcs/corrections.py | 149 ------------------------ hstwcs/det2im.py | 155 ------------------------- hstwcs/dgeo.py | 250 ----------------------------------------- hstwcs/makewcs.py | 240 --------------------------------------- updatewcs/__init__.py | 196 ++++++++++++++++++++++++++++++++ updatewcs/alignwcs.py | 95 ++++++++++++++++ updatewcs/apply_corrections.py | 181 +++++++++++++++++++++++++++++ updatewcs/corrections.py | 149 ++++++++++++++++++++++++ updatewcs/det2im.py | 155 +++++++++++++++++++++++++ updatewcs/dgeo.py | 250 +++++++++++++++++++++++++++++++++++++++++ updatewcs/makewcs.py | 240 +++++++++++++++++++++++++++++++++++++++ 14 files changed, 1266 insertions(+), 1266 deletions(-) delete mode 100644 hstwcs/__init__.py delete mode 100755 hstwcs/alignwcs.py delete mode 100644 hstwcs/apply_corrections.py delete mode 100644 hstwcs/corrections.py delete mode 100644 hstwcs/det2im.py delete mode 100644 hstwcs/dgeo.py delete mode 100644 hstwcs/makewcs.py create mode 100644 updatewcs/__init__.py create mode 100755 updatewcs/alignwcs.py create mode 100644 updatewcs/apply_corrections.py create mode 100644 updatewcs/corrections.py create mode 100644 updatewcs/det2im.py create mode 100644 updatewcs/dgeo.py create mode 100644 updatewcs/makewcs.py diff --git a/hstwcs/__init__.py b/hstwcs/__init__.py deleted file mode 100644 index bf47cdc..0000000 --- a/hstwcs/__init__.py +++ /dev/null @@ -1,196 +0,0 @@ -import os -import pyfits -#from .. wcsutil import HSTWCS -from updatewcs.wcsutil import HSTWCS - -#from .. mappings import allowed_corrections -from updatewcs import utils -import corrections, makewcs -import dgeo, det2im -from pytools import parseinput, fileutil -import apply_corrections - -#Note: The order of corrections is important - -__docformat__ = 'restructuredtext' - -__version__ = '0.3' - -def updatewcs(input, vacorr=True, tddcorr=True, dgeocorr=True, checkfiles=True, d2imcorr=True): - """ - Purpose - ======= - Applies corrections to the WCS keywords. - - Example - ======= - >>>from hstwcs import updatewcs - >>>updatewcs.updatewcs(filename) - - Dependencies - ============ - `pytools` - `pyfits` - `pywcs` - - :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 - `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. - """ - files = parseinput.parseinput(input)[0] - 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,dgeocorr=dgeocorr, d2imcorr=d2imcorr) - - #restore the original WCS keywords - utils.restoreWCS(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 - - """ - f = pyfits.open(fname, mode='update') - #Determine the reference chip and create the reference HSTWCS object - nrefchip, nrefext = getNrefchip(f) - ref_wcs = HSTWCS(fobj=f, ext=nrefext) - ref_wcs.readModel(update=True,header=f[nrefext].header) - utils.write_archive(f[nrefext].header) - - if 'DET2IMCorr' in allowed_corr: - det2im.DET2IMCorr.updateWCS(f) - - for i in range(len(f))[1:]: - # Perhaps all ext headers should be corrected (to be consistent) - extn = f[i] - if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': - ref_wcs.restore(f[nrefext].header) - hdr = extn.header - utils.write_archive(hdr) - ext_wcs = HSTWCS(fobj=f, ext=i) - ext_wcs.readModel(update=True,header=hdr) - for c in allowed_corr: - if c != 'DGEOCorr' and c != 'DET2IMCorr': - corr_klass = corrections.__getattribute__(c) - kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) - for kw in kw2update: - hdr.update(kw, kw2update[kw]) - - if 'DGEOCorr' in allowed_corr: - kw2update = dgeo.DGEOCorr.updateWCS(f) - for kw in kw2update: - f[1].header.update(kw, kw2update[kw]) - - - - - f.close() -def getNrefchip(fobj): - """ - This handles the fact that WFPC2 subarray observations - may not include chip 3 which is the default reference chip for - full observations. Also for subarrays chip 3 may not be the third - extension in a MEF file. - """ - Nrefext = 1 - instrument = fobj[0].header['INSTRUME'] - if instrument == 'WFPC2': - detectors = [img.header['DETECTOR'] for img in fobj[1:]] - - if 3 not in detectors: - Nrefchip=detectors[0] - Nrefext = 1 - else: - Nrefchip = 3 - Nrefext = detectors.index(3) + 1 - elif instrument == 'ACS': - detector = fobj[0].header['DETECTOR'] - if detector == 'WCS': - Nrefchip =2 - else: - Nrefchip = 1 - elif instrument == 'NICMOS': - Nrefchip = fobj[0].header['CAMERA'] - elif instrument == 'WFC3': - detector = fobj[0].header['DETECTOR'] - if detector == 'UVIS': - Nrefchip =2 - else: - Nrefchip = 1 - else: - Nrefchip = 1 - return Nrefchip, Nrefext - -def checkFiles(input): - """ - Purpose - ======= - Checks that input files are in the correct format. - Converts geis and waiver fits files to multiextension fits. - """ - from pytools.check_files import geis2mef, waiver2mef, checkFiles - removed_files = [] - newfiles = [] - for file in input: - try: - imgfits,imgtype = fileutil.isFits(file) - except IOError: - print "Warning: File %s could not be found\n" %file - print "Removing file %s from input list" %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: - print "Removing 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: - print "Removing file %s from input list - could not convert geis to mef" %file - removed_files.append(file) - else: - newfiles.append(newfilename) - if removed_files: - print 'The following files will be removed from the list of files to be processed :\n' - for f in removed_files: - print f - - newfiles = checkFiles(newfiles)[0] - - return newfiles - diff --git a/hstwcs/alignwcs.py b/hstwcs/alignwcs.py deleted file mode 100755 index db96ef6..0000000 --- a/hstwcs/alignwcs.py +++ /dev/null @@ -1,95 +0,0 @@ -#!/usr/bin/env python - -from pydrizzle import process_input -import pywcs, pyfits, numpy as np -from updatewcs.distortion import utils -from pytools import fileutil as fu -from updatewcs import wcsutil - - -def align_wcs(input, shiftfile=None, writeasn=False, asnname=None): - """ - Purpose - ======= - Correct the WCS of a fits file, so that multidrizzle aligns the images. - To view the format of the shift file: - >>>from pytools.asnutil import ShiftFile - >>>print ShiftFile.__doc__ - - Example - ======= - >>>import alignwcs - >>>alignwcs.align_wcs('*flt.fits', shiftfile='shifts.txt') - It works also on the command line: - %./alignwcs.py *flt.fits shifts.txt - - Dependencies - ============ - `pytools` - `pyfits` - `pywcs` - `numpy` - `pydrizzle` - - :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 - `shiftfile`: a shift file, as the output of tweakshifts - `writeasn`: boolean - defines whether to write out the asn table updated with shift information - `asnname`: string - name for the output asn table - """ - if shiftfile == None: - print 'A shift file is required but was not provided\n' - return - - asndict,b,c = process_input.process_input(input, shiftfile=shiftfile, updatewcs=False) - if writeasn: - if asnname == None: - asnname = c.split('.fits')[0] + '_shifts' + '.fits' - a.write(asnname) - outwcs = get_outwcs(asndict) - apply_shifts(asndict, outwcs) - -def apply_shifts(asndict, outwcs): - for mem in asndict['members']: - filename = fu.buildRootname(mem) - xsh = asndict['members'][mem]['xshift'] - ysh = asndict['members'][mem]['yshift'] - shifts = np.array([[xsh, ysh]]) - f = pyfits.open(filename) - for extn in f: - if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': - extver = extn.header['extver'] - owcs = pywcs.WCS(extn.header, f) - crvals = np.array([owcs.wcs.crval]) - px = outwcs.wcs.s2p(crvals, 1)['pixcrd'] + shifts - ncrvals = outwcs.all_pix2sky(px, 1) - pyfits.setval(filename, 'CRVAL1', value=ncrvals[0,0], ext=extver) - pyfits.setval(filename, 'CRVAL2', value=ncrvals[0,1], ext=extver) - print 'Updated %s with shifts ' % filename, shifts - f.close() - - -def get_outwcs(asndict): - wcslist = [] - for mem in asndict['members']: - filename = fu.buildRootname(mem) - f = pyfits.open(filename) - hdr0 = f[0].header - for extn in f: - if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': - owcs = wcsutil.HSTWCS(hdr0, extn.header,f) - wcslist.append(owcs) - f.close() - outwcs = utils.output_wcs(wcslist) - return outwcs - -if __name__ == '__main__': - import sys - args = sys.argv[1:] - input = args[:-1] - shifts = args[-1] - - align_wcs(input, shiftfile=shifts) diff --git a/hstwcs/apply_corrections.py b/hstwcs/apply_corrections.py deleted file mode 100644 index 88d997e..0000000 --- a/hstwcs/apply_corrections.py +++ /dev/null @@ -1,181 +0,0 @@ -import os -import pyfits -#import allowed_corrections -import time -from pytools import fileutil -import os.path -#Note: The order of corrections is important - -__docformat__ = 'restructuredtext' - -# A dictionary which lists the allowed corrections for each instrument. -# These are the default corrections applied also in the pipeline. - -allowed_corrections={'WFPC2': ['MakeWCS','CompSIP', 'VACorr', 'DGEOCorr'], - 'ACS': ['DET2IMCorr', 'TDDCorr', 'MakeWCS', 'CompSIP','VACorr', 'DGEOCorr'], - 'STIS': ['MakeWCS', 'CompSIP','VACorr'], - 'NICMOS': ['MakeWCS', 'CompSIP','VACorr'], - 'WFC3': ['MakeWCS', 'CompSIP','VACorr'], - } - -def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=True): - """ - Purpose - ======= - Creates a list of corrections to be applied to a file. - based on user input paramters and allowed corrections - for the instrument, which are defined in mappings.py. - """ - instrument = pyfits.getval(fname, 'INSTRUME') - tddcorr = applyTDDCorr(fname, tddcorr) - dgeocorr = applyDgeoCorr(fname, dgeocorr) - d2imcorr = applyD2ImCorr(fname, d2imcorr) - # make a copy of this list ! - acorr = allowed_corrections[instrument][:] - if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') - if 'TDDCorr' in acorr and tddcorr==False: acorr.remove('TDDCorr') - if 'DGEOCorr' in acorr and dgeocorr==False: acorr.remove('DGEOCorr') - if 'DET2IMCorr' in acorr and d2imcorr==False: acorr.remove('DET2IMCorr') - - return acorr - -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. - """ - - instrument = pyfits.getval(fname, 'INSTRUME') - try: - detector = pyfits.getval(fname, 'DETECTOR') - except KeyError: - detector = None - try: - tddcorr = pyfits.getval(fname, 'TDDCORR') - if tddcorr == 'PERFORM': - tddcorr = True - else: - tddcorr = False - except KeyError: - tddcorr = None - - if instrument == 'ACS' and detector == 'WFC' and tddcorr== True and utddcorr == True: - tddcorr = True - try: - idctab = pyfits.getval(fname, '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 applyDgeoCorr(fname, udgeocorr): - """ - Purpose - ======= - Adds dgeo extensions to files based on the DGEOFILE keyword in the primary - header. This is a default correction and will always run in the pipeline. - The file used to generate the extensions is - recorded in the DGEOEXT keyword in each science extension. - If 'DGEOFILE' in the primary header is different from 'DGEOEXT' in the - extension header and the file exists on disk and is a 'new type' dgeofile, - then the dgeo extensions will be updated. - """ - applyDGEOCorr = True - try: - # get DGEOFILE kw from primary header - fdgeo0 = pyfits.getval(fname, 'DGEOFILE') - if fdgeo0 == 'N/A': - return False - fdgeo0 = fileutil.osfn(fdgeo0) - if not fileutil.findFile(fdgeo0): - print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0 - print 'DGEO correction will not be applied\n' - applyDGEOCorr = False - return applyDGEOCorr - try: - # get DGEOEXT kw from first extension header - fdgeo1 = pyfits.getval(fname, 'DGEOEXT', ext=1) - fdgeo1 = fileutil.osfn(fdgeo1) - if fdgeo1 and fileutil.findFile(fdgeo1): - if fdgeo0 != fdgeo1: - applyDGEOCorr = True - else: - applyDGEOCorr = False - else: - # dgeo file defined in first extension may not be found - # but if a valid kw exists in the primary header, dgeo should be applied. - applyDGEOCorr = True - except KeyError: - # the case of DGEOFILE kw present in primary header but DGEOEXT missing - # in first extension header - applyDGEOCorr = True - except KeyError: - print 'DGEOFILE keyword not found in primary header' - applyDGEOCorr = False - return applyDGEOCorr - - if isOldStyleDGEO(fname, fdgeo0): - applyDGEOCorr = False - return (applyDGEOCorr and udgeocorr) - -def isOldStyleDGEO(fname, dgname): - # checks if the file defined in a DGEOFILE kw is a full size - # (old style) image - - sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1) - sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1) - dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1) - dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1) - if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2: - print 'Only full size (old style) XY file was found.' - print 'DGEO correction will not be applied.\n' - return True - else: - return False - -def applyD2ImCorr(fname, d2imcorr): - applyD2IMCorr = True - try: - # get DGEOFILE kw from primary header - fd2im0 = pyfits.getval(fname, 'D2IMFILE') - if fd2im0 == 'N/A': - return False - fd2im0 = fileutil.osfn(fd2im0) - if not fileutil.findFile(fd2im0): - print 'Kw D2IMFILE exists in primary header but file %s not found\n' % fd2im0 - print 'DGEO correction will not be applied\n' - applyD2IMCorr = False - return applyD2IMCorr - try: - # get DGEOEXT kw from first extension header - fd2imext = pyfits.getval(fname, 'D2IMEXT', ext=1) - fd2imext = fileutil.osfn(fd2imext) - if fd2imext and fileutil.findFile(fd2imext): - if fd2im0 != fd2imext: - applyD2IMCorr = True - else: - applyD2IMCorr = False - else: - # dgeo file defined in first extension may not be found - # but if a valid kw exists in the primary header, dgeo should be applied. - applyD2IMCorr = True - except KeyError: - # the case of DGEOFILE kw present in primary header but DGEOEXT missing - # in first extension header - applyD2IMCorr = True - except KeyError: - print 'D2IMFILE keyword not found in primary header' - applyD2IMCorr = False - return applyD2IMCorr - \ No newline at end of file diff --git a/hstwcs/corrections.py b/hstwcs/corrections.py deleted file mode 100644 index e76a260..0000000 --- a/hstwcs/corrections.py +++ /dev/null @@ -1,149 +0,0 @@ -import datetime -import numpy -from numpy import linalg -from pytools import fileutil -from updatewcs.utils import diff_angles -import makewcs, dgeo - -MakeWCS = makewcs.MakeWCS -DGEOCorr = dgeo.DGEOCorr - -class TDDCorr(object): - """ - Purpose - ======= - Apply time dependent distortion correction to SIP coefficients and basic - WCS keywords. It is applicable only to ACS/WFC data. - - Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion - Solution of the WFC - - :Parameters: - `ext_wcs`: HSTWCS object - An extension HSTWCS object to be modified - `ref_wcs`: HSTWCS object - A reference HSTWCS object - """ - - 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 - """ - - 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 = {'TDDALPHA': alpha, 'TDDBETA':beta, '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 newkw - updateWCS = classmethod(updateWCS) - - 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 = numpy.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],numpy.float64) - abmat1 = numpy.dot(tdd_mat, mrotn) - abmat2 = numpy.dot(mrotp,abmat1) - xshape, yshape = hwcs.idcmodel.cx.shape, hwcs.idcmodel.cy.shape - icxy = numpy.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 time dependent distortion skew terms - default date of 2004.5 = 2004-7-1 - """ - dday = 2004.5 - 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 - alpha = 0.095 + 0.090*(rday-dday)/2.5 - beta = -0.029 - 0.030*(rday-dday)/2.5 - - return alpha, beta - - compute_alpha_beta = classmethod(compute_alpha_beta) - - -class VACorr(object): - """ - Purpose - ======= - Apply velocity aberation correction to WCS keywords. - Modifies the CD matrix and CRVAL1/2 - - """ - - def updateWCS(cls, ext_wcs, ref_wcs): - if ext_wcs.vafactor != 1: - ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor - ext_wcs.wcs.crval[0] = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], ref_wcs.wcs.crval[0]) - ext_wcs.wcs.crval[1] = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], ref_wcs.wcs.crval[1]) - 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): - """ - Purpose - ======= - Compute SIP coefficients from idc table coefficients. - """ - - def updateWCS(cls, ext_wcs, ref_wcs): - 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 = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float) - imatr = linalg.inv(matr) - akeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) - bkeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) - for n in range(order+1): - for m in range(order+1): - if n >= m and n>=2: - idcval = numpy.array([[cx[n,m]],[cy[n,m]]]) - sipval = numpy.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] - kw2update[Bkey] = sipval[1,0] - kw2update['CTYPE1'] = 'RA---TAN-SIP' - kw2update['CTYPE2'] = 'DEC--TAN-SIP' - return kw2update - - updateWCS = classmethod(updateWCS) - - diff --git a/hstwcs/det2im.py b/hstwcs/det2im.py deleted file mode 100644 index 90c9f28..0000000 --- a/hstwcs/det2im.py +++ /dev/null @@ -1,155 +0,0 @@ -import pyfits -from pytools import fileutil - -class DET2IMCorr(object): - def updateWCS(cls, fobj): - """ - :Parameters: - `fobj`: pyfits object - Science file, for which a detector to image correction - is available - - """ - assert isinstance(fobj, pyfits.NP_pyfits.HDUList) - - d2imfile = fobj[0].header['D2IMFILE'] - axiscorr = cls.getAxisCorr(d2imfile) - d2imerr = pyfits.getdata(d2imfile, ext=1).max() - if axiscorr == None: - new_kw = {} - else: - new_kw = {'D2IMEXT': d2imfile, 'AXISCORR': axiscorr, 'D2IMERR': d2imerr} - cls.applyDet2ImCorr(fobj, axiscorr) - cls.updatehdr(fobj, new_kw) - - updateWCS = classmethod(updateWCS) - - def getAxisCorr(cls, refname): - try: - direction = pyfits.getval(refname, ext=1, key='EXTNAME') - if direction == 'DX': return 1 - elif direction == 'DY': return 2 - else: - print '\tDET2IM correction expects the reference file to have' - print '\tan EXTNAME keyword of value "DX" or "DY"' - return None - except AttributeError: - print "\tAxis to which to apply the detector to image " - print "\tcorrection cannot be determined because the reference " - print "\tfile %s is missing a keyword EXTNAME" % refname - direction = None - return direction - getAxisCorr = classmethod(getAxisCorr) - - def applyDet2ImCorr(cls,fobj, axiscorr): - d2imfile = fileutil.osfn(fobj[0].header['D2IMFILE']) - hdu = cls.createDgeoHDU(d2imfile, axiscorr) - d2imarr_ind = cls.getD2imIndex(fobj) - if d2imarr_ind: - fobj[d2imarr_ind] = hdu - else: - fobj.append(hdu) - applyDet2ImCorr = classmethod(applyDet2ImCorr) - - def getD2imIndex(cls,fobj): - index = None - for e in range(len(fobj)): - try: - ename = fobj[e].header['EXTNAME'] - except KeyError: - continue - if ename == 'D2IMARR': - index = e - return index - getD2imIndex = classmethod(getD2imIndex) - - def createDgeoHDU(cls, d2imfile, axiscorr): - - d2im_data = pyfits.getdata(d2imfile, ext=1) - d2im_hdr = cls.createDet2ImHdr(d2im_data.shape, axiscorr) - hdu = pyfits.ImageHDU(header=d2im_hdr, data=d2im_data) - - return hdu - - createDgeoHDU = classmethod(createDgeoHDU) - - def createDet2ImHdr(cls, data_shape, axiscorr): - """ - Creates a header for the D2IMARR extension based on the - reference file recorded in D2IMFILE keyword in the primary header. - """ - - naxis1 = data_shape[0] - naxis2 = 0 - crpix1 = 0.0 - crpix2 = 0.0 - cdelt1 = 1.0 - cdelt2 = 1.0 - crval1 = 0.0 - crval2 = 0.0 - keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', - 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', - 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2', 'AXISCORR'] - - comments = {'XTENSION': 'Image extension', - 'BITPIX': 'IEEE floating point', - 'NAXIS': 'Number of axes', - 'NAXIS1': 'Number of image columns', - 'NAXIS2': 'Number of image rows', - 'EXTNAME': 'WCS distortion array', - 'EXTVER': 'Distortion array version number', - 'PCOUNT': 'Special data area of size 0', - 'GCOUNT': 'One data group', - 'CRPIX1': 'Distortion array reference pixel', - 'CDELT1': 'Grid step size in first coordinate', - 'CRVAL1': 'Image array pixel coordinate', - 'CRPIX2': 'Distortion array reference pixel', - 'CDELT2': 'Grid step size in second coordinate', - 'CRVAL2': 'Image array pixel coordinate', - 'AXISCORR': 'Direction in which the det2im correction is applied'} - - values = {'XTENSION': 'IMAGE', - 'BITPIX': -32, - 'NAXIS': 1, - 'NAXIS1': naxis1, - 'NAXIS2': naxis2, - 'EXTNAME': 'D2IMARR', - 'EXTVER': 1, - 'PCOUNT': 0, - 'GCOUNT': 1, - 'CRPIX1': crpix1, - 'CDELT1': cdelt1, - 'CRVAL1': crval1, - 'CRPIX2': crpix2, - 'CDELT2': cdelt2, - 'CRVAL2': crval2, - 'AXISCORR': axiscorr - } - - - cdl = pyfits.CardList() - for c in keys: - cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c])) - - hdr = pyfits.Header(cards=cdl) - return hdr - - createDet2ImHdr = classmethod(createDet2ImHdr) - - def updatehdr(cls, fobj, kwdict): - """ - Update extension headers to keep record of the files used for the - detector to image correction. - """ - for ext in fobj: - try: - extname = ext.header['EXTNAME'].lower() - except KeyError: - continue - if extname == 'sci': - for kw in kwdict: - ext.header.update(kw, kwdict[kw]) - else: - continue - updatehdr = classmethod(updatehdr) - \ No newline at end of file diff --git a/hstwcs/dgeo.py b/hstwcs/dgeo.py deleted file mode 100644 index 00178d9..0000000 --- a/hstwcs/dgeo.py +++ /dev/null @@ -1,250 +0,0 @@ -import pyfits -from pytools import fileutil -import numpy as np - -class DGEOCorr(object): - """ - Purpose - ======= - Defines a Lookup table prior distortion correction as per WCS paper IV. - It uses a reference file defined by the DGEOFILE keyword in the primary header. - - Algorithm - ========= - - Using extensions in the reference file create a WCSDVARR extension - and add it to the file object. - - Add record-valued keywords which describe the lookup tables to the - science extension header - - Add a keyword 'DGEOFILE' to the science extension header, whose - value is the reference file used to create the WCSVARR extension - - If WCSDVARR extensions exist, subsequent updates will overwrite them. - If not, they will be added to the file object. - - It is assumed that the DGEO 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`: pyfits object - Science file, for which a distortion correction in a DGEOFILE is available - - """ - assert isinstance(fobj, pyfits.NP_pyfits.HDUList) - cls.applyDgeoCorr(fobj) - dgfile = fobj[0].header['DGEOFILE'] - - new_kw = {'DGEOEXT': dgfile} - return new_kw - - updateWCS = classmethod(updateWCS) - - def applyDgeoCorr(cls, fobj): - """ - For each science extension in a pyfits file object: - - create a WCSDVARR extension - - update science header - - add/update DGEOEXT keyword - """ - dgfile = fileutil.osfn(fobj[0].header['DGEOFILE']) - # 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'] - header = ext.header - # get the data arrays from the reference file and transform them for use with SIP - dx,dy = cls.getData(dgfile, extversion) - ndx, ndy = cls.transformData(header, dx,dy) - # determine EXTVER for the WCSDVARR extension from the DGEO file (EXTNAME, 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],[ndx, ndy]): - cls.addSciExtKw(header, wdvarr_ver=ename[1], dgeo_name=ename[0]) - hdu = cls.createDgeoHDU(header, dgeofile=dgfile, wdvarr_ver=ename[1],dgeo_name=ename[0], data=ename[2],extver=extversion) - if wcsdvarr_ind: - fobj[wcsdvarr_ind[ename[1]]] = hdu - else: - fobj.append(hdu) - - - applyDgeoCorr = classmethod(applyDgeoCorr) - - def getWCSIndex(cls, fobj): - """ - If fobj has WCSDVARR extensions: - returns a mapping of their EXTVER kw are mapped to 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 - - return wcsd - - getWCSIndex = classmethod(getWCSIndex) - - def addSciExtKw(cls, hdr, wdvarr_ver=None, dgeo_name=None): - """ - Adds kw to sci extension to define WCSDVARR lookup table extensions - - """ - if dgeo_name =='DX': - j=1 - else: - j=2 - - cperror = 'CPERROR%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: 0.0, cpdis: 'Lookup', dpext: wdvarr_ver, dpnaxes: 2, - dpaxis1: 1, dpaxis2: 2} - - comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (wdvarr_ver/2), - cpdis: 'Prior distortion funcion 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' - } - - for key in keys: - hdr.update(key=key, value=values[key], comment=comments[key], before='HISTORY') - - addSciExtKw = classmethod(addSciExtKw) - - def getData(cls,dgfile, extver): - """ - Get the data arrays from the reference DGEO files - """ - xdata = pyfits.getdata(dgfile, ext=('DX',extver)) - ydata = pyfits.getdata(dgfile, ext=('DY',extver)) - return xdata, ydata - getData = classmethod(getData) - - def transformData(cls, header, dx, dy): - """ - Transform the DGEO data arrays for use with SIP - """ - coeffs = cls.getCoeffs(header) - idcscale = header['IDCSCALE'] - sclcoeffs = np.linalg.inv(coeffs/idcscale) - ndx, ndy = np.dot(sclcoeffs, [dx.ravel(), dy.ravel()]) - ndx.shape = dx.shape - ndy.shape=dy.shape - return ndx, ndy - transformData = classmethod(transformData) - - def getCoeffs(cls, header): - """ - Return a matrix of the first order IDC coefficients. - """ - try: - ocx10 = header['OCX10'] - ocx11 = header['OCX11'] - ocy10 = header['OCY10'] - ocy11 = header['OCY11'] - except AttributeError: - print 'First order IDCTAB coefficients are not available.\n' - print 'Cannot convert SIP to IDC coefficients.\n' - return None - return np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32) - - getCoeffs = classmethod(getCoeffs) - - def createDgeoHDU(cls, sciheader, dgeofile=None, wdvarr_ver=1, dgeo_name=None,data = None, extver=1): - """ - Creates an HDU to be added to the file object. - """ - hdr = cls.createDgeoHdr(sciheader, dgeofile=dgeofile, wdvarr_ver=wdvarr_ver, dgeoname=dgeo_name, extver=extver) - hdu=pyfits.ImageHDU(header=hdr, data=data) - return hdu - - createDgeoHDU = classmethod(createDgeoHDU) - - def createDgeoHdr(cls, sciheader, dgeofile, wdvarr_ver, dgeoname, extver): - """ - Creates a header for the WCSDVARR extension based on the DGEO reference file - and sci extension header. - """ - dgeo_header = pyfits.getheader(dgeofile, ext=(dgeoname,extver)) - sci_naxis1 = sciheader['NAXIS1'] - sci_naxis2 = sciheader['NAXIS2'] - sci_crpix1 = sciheader['CRPIX1'] - sci_crpix2 = sciheader['CRPIX2'] - - naxis1 = dgeo_header['naxis1'] - naxis2 = dgeo_header['naxis2'] - extver = dgeo_header['extver'] - crpix1 = naxis1/2. - crpix2 = naxis2/2. - cdelt1 = sci_naxis1/naxis1 - cdelt2 = sci_naxis2/naxis2 - crval1 = sci_crpix1 - crval2 = sci_crpix2 - keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', - 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', - 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2'] - - comments = {'XTENSION': 'Image extension', - 'BITPIX': 'IEEE floating point', - 'NAXIS': 'Number of axes', - 'NAXIS1': 'Number of image columns', - 'NAXIS2': 'Number of image rows', - 'EXTNAME': 'WCS distortion array', - 'EXTVER': 'Distortion array version number', - 'PCOUNT': 'Special data area of size 0', - 'GCOUNT': 'One data group', - 'CRPIX1': 'Distortion array reference pixel', - 'CDELT1': 'Grid step size in first coordinate', - 'CRVAL1': 'Image array pixel coordinate', - 'CRPIX2': 'Distortion array reference pixel', - 'CDELT2': 'Grid step size in second coordinate', - 'CRVAL2': 'Image array pixel coordinate'} - - values = {'XTENSION': 'IMAGE', - 'BITPIX': -32, - 'NAXIS': 2, - 'NAXIS1': naxis1, - 'NAXIS2': naxis2, - 'EXTNAME': 'WCSDVARR', - 'EXTVER': wdvarr_ver, - 'PCOUNT': 0, - 'GCOUNT': 1, - 'CRPIX1': crpix1, - 'CDELT1': cdelt1, - 'CRVAL1': crval1, - 'CRPIX2': crpix2, - 'CDELT2': cdelt2, - 'CRVAL2': crval2 - } - - - cdl = pyfits.CardList() - for c in keys: - cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c])) - - hdr = pyfits.Header(cards=cdl) - return hdr - - createDgeoHdr = classmethod(createDgeoHdr) - diff --git a/hstwcs/makewcs.py b/hstwcs/makewcs.py deleted file mode 100644 index 49ed824..0000000 --- a/hstwcs/makewcs.py +++ /dev/null @@ -1,240 +0,0 @@ -from updatewcs import DEGTORAD, RADTODEG -import numpy as np -from math import sin, sqrt, pow, cos, asin, atan2,pi -from updatewcs import utils -from pytools import fileutil - -class MakeWCS(object): - """ - Purpose - ======= - Recompute basic WCS keywords based on PA_V3 and distortion model. - - Algorithm - ========= - - update reference chip wcs - - -- CRVAL: reference chip model zero point (XREF/YREF) on 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 - - - update extension 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) - - """ - tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]} - def updateWCS(cls, ext_wcs, ref_wcs): - """ - recomputes the basic WCS kw - """ - 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], - } - 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'] - fx,fy = ext_wcs.idcmodel.shift(ext_wcs.idcmodel.cx,ext_wcs.idcmodel.cy,offsetx,offsety) - else: - fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy - - 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] - offshift = offsh - dec = ref_wcs.wcs.crval[1] - tddscale = (ref_wcs.pscale/ext_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): - try: - alpha = hwcs.idcmodel.refpix['TDDALPHA'] - beta = hwcs.idcmodel.refpix['TDDBETA'] - except KeyError: - return np.array([[0., 0.],[0.,0.]]) - - 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 - - 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 the algorithm provided by Colin Cox that is used in - Generic Conversion at STScI. - """ - # Convert all angles to radians - _roll = DEGTORAD(roll) - _dec = DEGTORAD(dec) - _v2 = DEGTORAD(v2 / 3600.) - _v3 = DEGTORAD(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 = RADTODEG(pi - (gamma+B)) - - return troll - diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py new file mode 100644 index 0000000..bf47cdc --- /dev/null +++ b/updatewcs/__init__.py @@ -0,0 +1,196 @@ +import os +import pyfits +#from .. wcsutil import HSTWCS +from updatewcs.wcsutil import HSTWCS + +#from .. mappings import allowed_corrections +from updatewcs import utils +import corrections, makewcs +import dgeo, det2im +from pytools import parseinput, fileutil +import apply_corrections + +#Note: The order of corrections is important + +__docformat__ = 'restructuredtext' + +__version__ = '0.3' + +def updatewcs(input, vacorr=True, tddcorr=True, dgeocorr=True, checkfiles=True, d2imcorr=True): + """ + Purpose + ======= + Applies corrections to the WCS keywords. + + Example + ======= + >>>from hstwcs import updatewcs + >>>updatewcs.updatewcs(filename) + + Dependencies + ============ + `pytools` + `pyfits` + `pywcs` + + :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 + `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. + """ + files = parseinput.parseinput(input)[0] + 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,dgeocorr=dgeocorr, d2imcorr=d2imcorr) + + #restore the original WCS keywords + utils.restoreWCS(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 + + """ + f = pyfits.open(fname, mode='update') + #Determine the reference chip and create the reference HSTWCS object + nrefchip, nrefext = getNrefchip(f) + ref_wcs = HSTWCS(fobj=f, ext=nrefext) + ref_wcs.readModel(update=True,header=f[nrefext].header) + utils.write_archive(f[nrefext].header) + + if 'DET2IMCorr' in allowed_corr: + det2im.DET2IMCorr.updateWCS(f) + + for i in range(len(f))[1:]: + # Perhaps all ext headers should be corrected (to be consistent) + extn = f[i] + if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': + ref_wcs.restore(f[nrefext].header) + hdr = extn.header + utils.write_archive(hdr) + ext_wcs = HSTWCS(fobj=f, ext=i) + ext_wcs.readModel(update=True,header=hdr) + for c in allowed_corr: + if c != 'DGEOCorr' and c != 'DET2IMCorr': + corr_klass = corrections.__getattribute__(c) + kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) + for kw in kw2update: + hdr.update(kw, kw2update[kw]) + + if 'DGEOCorr' in allowed_corr: + kw2update = dgeo.DGEOCorr.updateWCS(f) + for kw in kw2update: + f[1].header.update(kw, kw2update[kw]) + + + + + f.close() +def getNrefchip(fobj): + """ + This handles the fact that WFPC2 subarray observations + may not include chip 3 which is the default reference chip for + full observations. Also for subarrays chip 3 may not be the third + extension in a MEF file. + """ + Nrefext = 1 + instrument = fobj[0].header['INSTRUME'] + if instrument == 'WFPC2': + detectors = [img.header['DETECTOR'] for img in fobj[1:]] + + if 3 not in detectors: + Nrefchip=detectors[0] + Nrefext = 1 + else: + Nrefchip = 3 + Nrefext = detectors.index(3) + 1 + elif instrument == 'ACS': + detector = fobj[0].header['DETECTOR'] + if detector == 'WCS': + Nrefchip =2 + else: + Nrefchip = 1 + elif instrument == 'NICMOS': + Nrefchip = fobj[0].header['CAMERA'] + elif instrument == 'WFC3': + detector = fobj[0].header['DETECTOR'] + if detector == 'UVIS': + Nrefchip =2 + else: + Nrefchip = 1 + else: + Nrefchip = 1 + return Nrefchip, Nrefext + +def checkFiles(input): + """ + Purpose + ======= + Checks that input files are in the correct format. + Converts geis and waiver fits files to multiextension fits. + """ + from pytools.check_files import geis2mef, waiver2mef, checkFiles + removed_files = [] + newfiles = [] + for file in input: + try: + imgfits,imgtype = fileutil.isFits(file) + except IOError: + print "Warning: File %s could not be found\n" %file + print "Removing file %s from input list" %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: + print "Removing 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: + print "Removing file %s from input list - could not convert geis to mef" %file + removed_files.append(file) + else: + newfiles.append(newfilename) + if removed_files: + print 'The following files will be removed from the list of files to be processed :\n' + for f in removed_files: + print f + + newfiles = checkFiles(newfiles)[0] + + return newfiles + diff --git a/updatewcs/alignwcs.py b/updatewcs/alignwcs.py new file mode 100755 index 0000000..db96ef6 --- /dev/null +++ b/updatewcs/alignwcs.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python + +from pydrizzle import process_input +import pywcs, pyfits, numpy as np +from updatewcs.distortion import utils +from pytools import fileutil as fu +from updatewcs import wcsutil + + +def align_wcs(input, shiftfile=None, writeasn=False, asnname=None): + """ + Purpose + ======= + Correct the WCS of a fits file, so that multidrizzle aligns the images. + To view the format of the shift file: + >>>from pytools.asnutil import ShiftFile + >>>print ShiftFile.__doc__ + + Example + ======= + >>>import alignwcs + >>>alignwcs.align_wcs('*flt.fits', shiftfile='shifts.txt') + It works also on the command line: + %./alignwcs.py *flt.fits shifts.txt + + Dependencies + ============ + `pytools` + `pyfits` + `pywcs` + `numpy` + `pydrizzle` + + :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 + `shiftfile`: a shift file, as the output of tweakshifts + `writeasn`: boolean + defines whether to write out the asn table updated with shift information + `asnname`: string + name for the output asn table + """ + if shiftfile == None: + print 'A shift file is required but was not provided\n' + return + + asndict,b,c = process_input.process_input(input, shiftfile=shiftfile, updatewcs=False) + if writeasn: + if asnname == None: + asnname = c.split('.fits')[0] + '_shifts' + '.fits' + a.write(asnname) + outwcs = get_outwcs(asndict) + apply_shifts(asndict, outwcs) + +def apply_shifts(asndict, outwcs): + for mem in asndict['members']: + filename = fu.buildRootname(mem) + xsh = asndict['members'][mem]['xshift'] + ysh = asndict['members'][mem]['yshift'] + shifts = np.array([[xsh, ysh]]) + f = pyfits.open(filename) + for extn in f: + if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': + extver = extn.header['extver'] + owcs = pywcs.WCS(extn.header, f) + crvals = np.array([owcs.wcs.crval]) + px = outwcs.wcs.s2p(crvals, 1)['pixcrd'] + shifts + ncrvals = outwcs.all_pix2sky(px, 1) + pyfits.setval(filename, 'CRVAL1', value=ncrvals[0,0], ext=extver) + pyfits.setval(filename, 'CRVAL2', value=ncrvals[0,1], ext=extver) + print 'Updated %s with shifts ' % filename, shifts + f.close() + + +def get_outwcs(asndict): + wcslist = [] + for mem in asndict['members']: + filename = fu.buildRootname(mem) + f = pyfits.open(filename) + hdr0 = f[0].header + for extn in f: + if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': + owcs = wcsutil.HSTWCS(hdr0, extn.header,f) + wcslist.append(owcs) + f.close() + outwcs = utils.output_wcs(wcslist) + return outwcs + +if __name__ == '__main__': + import sys + args = sys.argv[1:] + input = args[:-1] + shifts = args[-1] + + align_wcs(input, shiftfile=shifts) diff --git a/updatewcs/apply_corrections.py b/updatewcs/apply_corrections.py new file mode 100644 index 0000000..88d997e --- /dev/null +++ b/updatewcs/apply_corrections.py @@ -0,0 +1,181 @@ +import os +import pyfits +#import allowed_corrections +import time +from pytools import fileutil +import os.path +#Note: The order of corrections is important + +__docformat__ = 'restructuredtext' + +# A dictionary which lists the allowed corrections for each instrument. +# These are the default corrections applied also in the pipeline. + +allowed_corrections={'WFPC2': ['MakeWCS','CompSIP', 'VACorr', 'DGEOCorr'], + 'ACS': ['DET2IMCorr', 'TDDCorr', 'MakeWCS', 'CompSIP','VACorr', 'DGEOCorr'], + 'STIS': ['MakeWCS', 'CompSIP','VACorr'], + 'NICMOS': ['MakeWCS', 'CompSIP','VACorr'], + 'WFC3': ['MakeWCS', 'CompSIP','VACorr'], + } + +def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=True): + """ + Purpose + ======= + Creates a list of corrections to be applied to a file. + based on user input paramters and allowed corrections + for the instrument, which are defined in mappings.py. + """ + instrument = pyfits.getval(fname, 'INSTRUME') + tddcorr = applyTDDCorr(fname, tddcorr) + dgeocorr = applyDgeoCorr(fname, dgeocorr) + d2imcorr = applyD2ImCorr(fname, d2imcorr) + # make a copy of this list ! + acorr = allowed_corrections[instrument][:] + if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') + if 'TDDCorr' in acorr and tddcorr==False: acorr.remove('TDDCorr') + if 'DGEOCorr' in acorr and dgeocorr==False: acorr.remove('DGEOCorr') + if 'DET2IMCorr' in acorr and d2imcorr==False: acorr.remove('DET2IMCorr') + + return acorr + +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. + """ + + instrument = pyfits.getval(fname, 'INSTRUME') + try: + detector = pyfits.getval(fname, 'DETECTOR') + except KeyError: + detector = None + try: + tddcorr = pyfits.getval(fname, 'TDDCORR') + if tddcorr == 'PERFORM': + tddcorr = True + else: + tddcorr = False + except KeyError: + tddcorr = None + + if instrument == 'ACS' and detector == 'WFC' and tddcorr== True and utddcorr == True: + tddcorr = True + try: + idctab = pyfits.getval(fname, '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 applyDgeoCorr(fname, udgeocorr): + """ + Purpose + ======= + Adds dgeo extensions to files based on the DGEOFILE keyword in the primary + header. This is a default correction and will always run in the pipeline. + The file used to generate the extensions is + recorded in the DGEOEXT keyword in each science extension. + If 'DGEOFILE' in the primary header is different from 'DGEOEXT' in the + extension header and the file exists on disk and is a 'new type' dgeofile, + then the dgeo extensions will be updated. + """ + applyDGEOCorr = True + try: + # get DGEOFILE kw from primary header + fdgeo0 = pyfits.getval(fname, 'DGEOFILE') + if fdgeo0 == 'N/A': + return False + fdgeo0 = fileutil.osfn(fdgeo0) + if not fileutil.findFile(fdgeo0): + print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0 + print 'DGEO correction will not be applied\n' + applyDGEOCorr = False + return applyDGEOCorr + try: + # get DGEOEXT kw from first extension header + fdgeo1 = pyfits.getval(fname, 'DGEOEXT', ext=1) + fdgeo1 = fileutil.osfn(fdgeo1) + if fdgeo1 and fileutil.findFile(fdgeo1): + if fdgeo0 != fdgeo1: + applyDGEOCorr = True + else: + applyDGEOCorr = False + else: + # dgeo file defined in first extension may not be found + # but if a valid kw exists in the primary header, dgeo should be applied. + applyDGEOCorr = True + except KeyError: + # the case of DGEOFILE kw present in primary header but DGEOEXT missing + # in first extension header + applyDGEOCorr = True + except KeyError: + print 'DGEOFILE keyword not found in primary header' + applyDGEOCorr = False + return applyDGEOCorr + + if isOldStyleDGEO(fname, fdgeo0): + applyDGEOCorr = False + return (applyDGEOCorr and udgeocorr) + +def isOldStyleDGEO(fname, dgname): + # checks if the file defined in a DGEOFILE kw is a full size + # (old style) image + + sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1) + sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1) + dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1) + dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1) + if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2: + print 'Only full size (old style) XY file was found.' + print 'DGEO correction will not be applied.\n' + return True + else: + return False + +def applyD2ImCorr(fname, d2imcorr): + applyD2IMCorr = True + try: + # get DGEOFILE kw from primary header + fd2im0 = pyfits.getval(fname, 'D2IMFILE') + if fd2im0 == 'N/A': + return False + fd2im0 = fileutil.osfn(fd2im0) + if not fileutil.findFile(fd2im0): + print 'Kw D2IMFILE exists in primary header but file %s not found\n' % fd2im0 + print 'DGEO correction will not be applied\n' + applyD2IMCorr = False + return applyD2IMCorr + try: + # get DGEOEXT kw from first extension header + fd2imext = pyfits.getval(fname, 'D2IMEXT', ext=1) + fd2imext = fileutil.osfn(fd2imext) + if fd2imext and fileutil.findFile(fd2imext): + if fd2im0 != fd2imext: + applyD2IMCorr = True + else: + applyD2IMCorr = False + else: + # dgeo file defined in first extension may not be found + # but if a valid kw exists in the primary header, dgeo should be applied. + applyD2IMCorr = True + except KeyError: + # the case of DGEOFILE kw present in primary header but DGEOEXT missing + # in first extension header + applyD2IMCorr = True + except KeyError: + print 'D2IMFILE keyword not found in primary header' + applyD2IMCorr = False + return applyD2IMCorr + \ No newline at end of file diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py new file mode 100644 index 0000000..e76a260 --- /dev/null +++ b/updatewcs/corrections.py @@ -0,0 +1,149 @@ +import datetime +import numpy +from numpy import linalg +from pytools import fileutil +from updatewcs.utils import diff_angles +import makewcs, dgeo + +MakeWCS = makewcs.MakeWCS +DGEOCorr = dgeo.DGEOCorr + +class TDDCorr(object): + """ + Purpose + ======= + Apply time dependent distortion correction to SIP coefficients and basic + WCS keywords. It is applicable only to ACS/WFC data. + + Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion + Solution of the WFC + + :Parameters: + `ext_wcs`: HSTWCS object + An extension HSTWCS object to be modified + `ref_wcs`: HSTWCS object + A reference HSTWCS object + """ + + 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 + """ + + 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 = {'TDDALPHA': alpha, 'TDDBETA':beta, '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 newkw + updateWCS = classmethod(updateWCS) + + 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 = numpy.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],numpy.float64) + abmat1 = numpy.dot(tdd_mat, mrotn) + abmat2 = numpy.dot(mrotp,abmat1) + xshape, yshape = hwcs.idcmodel.cx.shape, hwcs.idcmodel.cy.shape + icxy = numpy.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 time dependent distortion skew terms + default date of 2004.5 = 2004-7-1 + """ + dday = 2004.5 + 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 + alpha = 0.095 + 0.090*(rday-dday)/2.5 + beta = -0.029 - 0.030*(rday-dday)/2.5 + + return alpha, beta + + compute_alpha_beta = classmethod(compute_alpha_beta) + + +class VACorr(object): + """ + Purpose + ======= + Apply velocity aberation correction to WCS keywords. + Modifies the CD matrix and CRVAL1/2 + + """ + + def updateWCS(cls, ext_wcs, ref_wcs): + if ext_wcs.vafactor != 1: + ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor + ext_wcs.wcs.crval[0] = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], ref_wcs.wcs.crval[0]) + ext_wcs.wcs.crval[1] = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], ref_wcs.wcs.crval[1]) + 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): + """ + Purpose + ======= + Compute SIP coefficients from idc table coefficients. + """ + + def updateWCS(cls, ext_wcs, ref_wcs): + 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 = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float) + imatr = linalg.inv(matr) + akeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) + bkeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) + for n in range(order+1): + for m in range(order+1): + if n >= m and n>=2: + idcval = numpy.array([[cx[n,m]],[cy[n,m]]]) + sipval = numpy.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] + kw2update[Bkey] = sipval[1,0] + kw2update['CTYPE1'] = 'RA---TAN-SIP' + kw2update['CTYPE2'] = 'DEC--TAN-SIP' + return kw2update + + updateWCS = classmethod(updateWCS) + + diff --git a/updatewcs/det2im.py b/updatewcs/det2im.py new file mode 100644 index 0000000..90c9f28 --- /dev/null +++ b/updatewcs/det2im.py @@ -0,0 +1,155 @@ +import pyfits +from pytools import fileutil + +class DET2IMCorr(object): + def updateWCS(cls, fobj): + """ + :Parameters: + `fobj`: pyfits object + Science file, for which a detector to image correction + is available + + """ + assert isinstance(fobj, pyfits.NP_pyfits.HDUList) + + d2imfile = fobj[0].header['D2IMFILE'] + axiscorr = cls.getAxisCorr(d2imfile) + d2imerr = pyfits.getdata(d2imfile, ext=1).max() + if axiscorr == None: + new_kw = {} + else: + new_kw = {'D2IMEXT': d2imfile, 'AXISCORR': axiscorr, 'D2IMERR': d2imerr} + cls.applyDet2ImCorr(fobj, axiscorr) + cls.updatehdr(fobj, new_kw) + + updateWCS = classmethod(updateWCS) + + def getAxisCorr(cls, refname): + try: + direction = pyfits.getval(refname, ext=1, key='EXTNAME') + if direction == 'DX': return 1 + elif direction == 'DY': return 2 + else: + print '\tDET2IM correction expects the reference file to have' + print '\tan EXTNAME keyword of value "DX" or "DY"' + return None + except AttributeError: + print "\tAxis to which to apply the detector to image " + print "\tcorrection cannot be determined because the reference " + print "\tfile %s is missing a keyword EXTNAME" % refname + direction = None + return direction + getAxisCorr = classmethod(getAxisCorr) + + def applyDet2ImCorr(cls,fobj, axiscorr): + d2imfile = fileutil.osfn(fobj[0].header['D2IMFILE']) + hdu = cls.createDgeoHDU(d2imfile, axiscorr) + d2imarr_ind = cls.getD2imIndex(fobj) + if d2imarr_ind: + fobj[d2imarr_ind] = hdu + else: + fobj.append(hdu) + applyDet2ImCorr = classmethod(applyDet2ImCorr) + + def getD2imIndex(cls,fobj): + index = None + for e in range(len(fobj)): + try: + ename = fobj[e].header['EXTNAME'] + except KeyError: + continue + if ename == 'D2IMARR': + index = e + return index + getD2imIndex = classmethod(getD2imIndex) + + def createDgeoHDU(cls, d2imfile, axiscorr): + + d2im_data = pyfits.getdata(d2imfile, ext=1) + d2im_hdr = cls.createDet2ImHdr(d2im_data.shape, axiscorr) + hdu = pyfits.ImageHDU(header=d2im_hdr, data=d2im_data) + + return hdu + + createDgeoHDU = classmethod(createDgeoHDU) + + def createDet2ImHdr(cls, data_shape, axiscorr): + """ + Creates a header for the D2IMARR extension based on the + reference file recorded in D2IMFILE keyword in the primary header. + """ + + naxis1 = data_shape[0] + naxis2 = 0 + crpix1 = 0.0 + crpix2 = 0.0 + cdelt1 = 1.0 + cdelt2 = 1.0 + crval1 = 0.0 + crval2 = 0.0 + keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', + 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', + 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2', 'AXISCORR'] + + comments = {'XTENSION': 'Image extension', + 'BITPIX': 'IEEE floating point', + 'NAXIS': 'Number of axes', + 'NAXIS1': 'Number of image columns', + 'NAXIS2': 'Number of image rows', + 'EXTNAME': 'WCS distortion array', + 'EXTVER': 'Distortion array version number', + 'PCOUNT': 'Special data area of size 0', + 'GCOUNT': 'One data group', + 'CRPIX1': 'Distortion array reference pixel', + 'CDELT1': 'Grid step size in first coordinate', + 'CRVAL1': 'Image array pixel coordinate', + 'CRPIX2': 'Distortion array reference pixel', + 'CDELT2': 'Grid step size in second coordinate', + 'CRVAL2': 'Image array pixel coordinate', + 'AXISCORR': 'Direction in which the det2im correction is applied'} + + values = {'XTENSION': 'IMAGE', + 'BITPIX': -32, + 'NAXIS': 1, + 'NAXIS1': naxis1, + 'NAXIS2': naxis2, + 'EXTNAME': 'D2IMARR', + 'EXTVER': 1, + 'PCOUNT': 0, + 'GCOUNT': 1, + 'CRPIX1': crpix1, + 'CDELT1': cdelt1, + 'CRVAL1': crval1, + 'CRPIX2': crpix2, + 'CDELT2': cdelt2, + 'CRVAL2': crval2, + 'AXISCORR': axiscorr + } + + + cdl = pyfits.CardList() + for c in keys: + cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c])) + + hdr = pyfits.Header(cards=cdl) + return hdr + + createDet2ImHdr = classmethod(createDet2ImHdr) + + def updatehdr(cls, fobj, kwdict): + """ + Update extension headers to keep record of the files used for the + detector to image correction. + """ + for ext in fobj: + try: + extname = ext.header['EXTNAME'].lower() + except KeyError: + continue + if extname == 'sci': + for kw in kwdict: + ext.header.update(kw, kwdict[kw]) + else: + continue + updatehdr = classmethod(updatehdr) + \ No newline at end of file diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py new file mode 100644 index 0000000..00178d9 --- /dev/null +++ b/updatewcs/dgeo.py @@ -0,0 +1,250 @@ +import pyfits +from pytools import fileutil +import numpy as np + +class DGEOCorr(object): + """ + Purpose + ======= + Defines a Lookup table prior distortion correction as per WCS paper IV. + It uses a reference file defined by the DGEOFILE keyword in the primary header. + + Algorithm + ========= + - Using extensions in the reference file create a WCSDVARR extension + and add it to the file object. + - Add record-valued keywords which describe the lookup tables to the + science extension header + - Add a keyword 'DGEOFILE' to the science extension header, whose + value is the reference file used to create the WCSVARR extension + + If WCSDVARR extensions exist, subsequent updates will overwrite them. + If not, they will be added to the file object. + + It is assumed that the DGEO 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`: pyfits object + Science file, for which a distortion correction in a DGEOFILE is available + + """ + assert isinstance(fobj, pyfits.NP_pyfits.HDUList) + cls.applyDgeoCorr(fobj) + dgfile = fobj[0].header['DGEOFILE'] + + new_kw = {'DGEOEXT': dgfile} + return new_kw + + updateWCS = classmethod(updateWCS) + + def applyDgeoCorr(cls, fobj): + """ + For each science extension in a pyfits file object: + - create a WCSDVARR extension + - update science header + - add/update DGEOEXT keyword + """ + dgfile = fileutil.osfn(fobj[0].header['DGEOFILE']) + # 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'] + header = ext.header + # get the data arrays from the reference file and transform them for use with SIP + dx,dy = cls.getData(dgfile, extversion) + ndx, ndy = cls.transformData(header, dx,dy) + # determine EXTVER for the WCSDVARR extension from the DGEO file (EXTNAME, 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],[ndx, ndy]): + cls.addSciExtKw(header, wdvarr_ver=ename[1], dgeo_name=ename[0]) + hdu = cls.createDgeoHDU(header, dgeofile=dgfile, wdvarr_ver=ename[1],dgeo_name=ename[0], data=ename[2],extver=extversion) + if wcsdvarr_ind: + fobj[wcsdvarr_ind[ename[1]]] = hdu + else: + fobj.append(hdu) + + + applyDgeoCorr = classmethod(applyDgeoCorr) + + def getWCSIndex(cls, fobj): + """ + If fobj has WCSDVARR extensions: + returns a mapping of their EXTVER kw are mapped to 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 + + return wcsd + + getWCSIndex = classmethod(getWCSIndex) + + def addSciExtKw(cls, hdr, wdvarr_ver=None, dgeo_name=None): + """ + Adds kw to sci extension to define WCSDVARR lookup table extensions + + """ + if dgeo_name =='DX': + j=1 + else: + j=2 + + cperror = 'CPERROR%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: 0.0, cpdis: 'Lookup', dpext: wdvarr_ver, dpnaxes: 2, + dpaxis1: 1, dpaxis2: 2} + + comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (wdvarr_ver/2), + cpdis: 'Prior distortion funcion 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' + } + + for key in keys: + hdr.update(key=key, value=values[key], comment=comments[key], before='HISTORY') + + addSciExtKw = classmethod(addSciExtKw) + + def getData(cls,dgfile, extver): + """ + Get the data arrays from the reference DGEO files + """ + xdata = pyfits.getdata(dgfile, ext=('DX',extver)) + ydata = pyfits.getdata(dgfile, ext=('DY',extver)) + return xdata, ydata + getData = classmethod(getData) + + def transformData(cls, header, dx, dy): + """ + Transform the DGEO data arrays for use with SIP + """ + coeffs = cls.getCoeffs(header) + idcscale = header['IDCSCALE'] + sclcoeffs = np.linalg.inv(coeffs/idcscale) + ndx, ndy = np.dot(sclcoeffs, [dx.ravel(), dy.ravel()]) + ndx.shape = dx.shape + ndy.shape=dy.shape + return ndx, ndy + transformData = classmethod(transformData) + + def getCoeffs(cls, header): + """ + Return a matrix of the first order IDC coefficients. + """ + try: + ocx10 = header['OCX10'] + ocx11 = header['OCX11'] + ocy10 = header['OCY10'] + ocy11 = header['OCY11'] + except AttributeError: + print 'First order IDCTAB coefficients are not available.\n' + print 'Cannot convert SIP to IDC coefficients.\n' + return None + return np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32) + + getCoeffs = classmethod(getCoeffs) + + def createDgeoHDU(cls, sciheader, dgeofile=None, wdvarr_ver=1, dgeo_name=None,data = None, extver=1): + """ + Creates an HDU to be added to the file object. + """ + hdr = cls.createDgeoHdr(sciheader, dgeofile=dgeofile, wdvarr_ver=wdvarr_ver, dgeoname=dgeo_name, extver=extver) + hdu=pyfits.ImageHDU(header=hdr, data=data) + return hdu + + createDgeoHDU = classmethod(createDgeoHDU) + + def createDgeoHdr(cls, sciheader, dgeofile, wdvarr_ver, dgeoname, extver): + """ + Creates a header for the WCSDVARR extension based on the DGEO reference file + and sci extension header. + """ + dgeo_header = pyfits.getheader(dgeofile, ext=(dgeoname,extver)) + sci_naxis1 = sciheader['NAXIS1'] + sci_naxis2 = sciheader['NAXIS2'] + sci_crpix1 = sciheader['CRPIX1'] + sci_crpix2 = sciheader['CRPIX2'] + + naxis1 = dgeo_header['naxis1'] + naxis2 = dgeo_header['naxis2'] + extver = dgeo_header['extver'] + crpix1 = naxis1/2. + crpix2 = naxis2/2. + cdelt1 = sci_naxis1/naxis1 + cdelt2 = sci_naxis2/naxis2 + crval1 = sci_crpix1 + crval2 = sci_crpix2 + keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', + 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', + 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2'] + + comments = {'XTENSION': 'Image extension', + 'BITPIX': 'IEEE floating point', + 'NAXIS': 'Number of axes', + 'NAXIS1': 'Number of image columns', + 'NAXIS2': 'Number of image rows', + 'EXTNAME': 'WCS distortion array', + 'EXTVER': 'Distortion array version number', + 'PCOUNT': 'Special data area of size 0', + 'GCOUNT': 'One data group', + 'CRPIX1': 'Distortion array reference pixel', + 'CDELT1': 'Grid step size in first coordinate', + 'CRVAL1': 'Image array pixel coordinate', + 'CRPIX2': 'Distortion array reference pixel', + 'CDELT2': 'Grid step size in second coordinate', + 'CRVAL2': 'Image array pixel coordinate'} + + values = {'XTENSION': 'IMAGE', + 'BITPIX': -32, + 'NAXIS': 2, + 'NAXIS1': naxis1, + 'NAXIS2': naxis2, + 'EXTNAME': 'WCSDVARR', + 'EXTVER': wdvarr_ver, + 'PCOUNT': 0, + 'GCOUNT': 1, + 'CRPIX1': crpix1, + 'CDELT1': cdelt1, + 'CRVAL1': crval1, + 'CRPIX2': crpix2, + 'CDELT2': cdelt2, + 'CRVAL2': crval2 + } + + + cdl = pyfits.CardList() + for c in keys: + cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c])) + + hdr = pyfits.Header(cards=cdl) + return hdr + + createDgeoHdr = classmethod(createDgeoHdr) + diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py new file mode 100644 index 0000000..49ed824 --- /dev/null +++ b/updatewcs/makewcs.py @@ -0,0 +1,240 @@ +from updatewcs import DEGTORAD, RADTODEG +import numpy as np +from math import sin, sqrt, pow, cos, asin, atan2,pi +from updatewcs import utils +from pytools import fileutil + +class MakeWCS(object): + """ + Purpose + ======= + Recompute basic WCS keywords based on PA_V3 and distortion model. + + Algorithm + ========= + - update reference chip wcs + + -- CRVAL: reference chip model zero point (XREF/YREF) on 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 + + - update extension 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) + + """ + tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]} + def updateWCS(cls, ext_wcs, ref_wcs): + """ + recomputes the basic WCS kw + """ + 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], + } + 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'] + fx,fy = ext_wcs.idcmodel.shift(ext_wcs.idcmodel.cx,ext_wcs.idcmodel.cy,offsetx,offsety) + else: + fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy + + 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] + offshift = offsh + dec = ref_wcs.wcs.crval[1] + tddscale = (ref_wcs.pscale/ext_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): + try: + alpha = hwcs.idcmodel.refpix['TDDALPHA'] + beta = hwcs.idcmodel.refpix['TDDBETA'] + except KeyError: + return np.array([[0., 0.],[0.,0.]]) + + 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 + + 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 the algorithm provided by Colin Cox that is used in + Generic Conversion at STScI. + """ + # Convert all angles to radians + _roll = DEGTORAD(roll) + _dec = DEGTORAD(dec) + _v2 = DEGTORAD(v2 / 3600.) + _v3 = DEGTORAD(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 = RADTODEG(pi - (gamma+B)) + + return troll + -- cgit