summaryrefslogtreecommitdiff
path: root/lib/utils.py
blob: 6f63db031078190e2175c197f11cef7d60341a2f (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
from __future__ import division # confidence high

from pytools import parseinput, fileutil
import pyfits
import pywcs
import numpy as np
import string 

def restoreWCS(fnames, wcskey, clobber=False):
    """
    Purpose
    =======
    Reads in a WCS defined with wcskey and saves it as the primary WCS.
    If clobber is False, writes out new files whose names are the original 
    names with an attached 3 character string representing _'WCSKEY'_. 
    Otherwise overwrites the files. The WCS is restored from the 'SCI' 
    extension but the primary WCS of all extensions are updated.
    
    WCS keywords:
    'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 
    'CRVAL*',
    'CTYPE*',
    'CRPIX*', 
    'CDELT*',
    'CUNIT*',
    'ORIENTAT' - ?
    'TDDALPHA', 'TDDBETA'
    'A_x_x', B_x_x' - SIP coefficients
    'CPERROR*', 'CPDIS*', 'DP*', 
    
    `fnames`: a python list of file names, a string of comma separated file names, 
              an @file
    `wcskey`: a charater
              Used for one of 26 alternate WCS definitions.
    `clobber`: boolean
              A flag to define if the original files should be overwritten
    """
    files = parseinput.parseinput(fnames)[0]
    for f in files:
        isfits, ftype = fileutil.isFits(f)
        if not isfits or (isfits and ftype == 'waiver'):
            print "RestoreWCS works only with true fits files."
            return
        else:
            if clobber:
                print 'Overwriting original files\n'
                fobj = pyfits.open(f, mode='update')
                name = f
            else:
                fobj = pyfits.open(f)
                name = (f.split('.fits')[0] + '_%s_' + '.fits') %wcskey
            for e in range(len(fobj)):
                try:
                    extname = fobj[e].header['EXTNAME'].lower()
                except KeyError:
                    continue
                #Restore always from a 'SCI' extension but write it out to 'ERR' and 'DQ'
                if extname == 'sci':
                    sciver = fobj[e].header['extver']
                    try:
                        nwcs = pywcs.WCS(fobj[e].header, fobj=fobj, key=wcskey)
                    except:
                        print 'utils.restoreWCS: Could not read WCS with key %s in file %s,  \
                        extension %d' % (wcskey, f, e)
                        return #raise
                    hwcs = nwcs.to_header()
                    if nwcs.wcs.has_pc():
                        for c in ['1_1', '1_2', '2_1', '2_2']:
                            del hwcs['CD'+c+wcskey]
                    elif nwcs.wcs.has_cd():
                        for c in ['1_1', '1_2', '2_1', '2_2']:
                            hwcs.update(key='CD'+c+wcskey, value=hwcs['PC'+c+wcskey])
                            del hwcs['PC'+c]
                    for k in hwcs.keys():
                        key = k[:-1]
                        if key in fobj[e].header.keys():
                            fobj[e].header.update(key=key, value = hwcs[k])
                        else:
                            continue
                    if wcskey == 'O' and fobj[e].header.has_key('TDDALPHA'):
                        fobj[e].header['TDDALPHA'] = 0.0
                        fobj[e].header['TDDBETA'] = 0.0
                    if fobj[e].header.has_key('ORIENTAT'):
                        cd12 = 'CD1_2%s' % wcskey
                        cd22 = 'CD2_2%s' % wcskey
                        norient = np.rad2deg(np.arctan2(hwcs[cd12],hwcs[cd22]))
                        fobj[e].header.update(key='ORIENTAT', value=norient)
                elif extname in ['err', 'dq', 'sdq']:
                    cextver = fobj[e].header['extver']
                    if cextver == sciver:
                        for k in hwcs.keys():
                            key = k[:-1]
                            fobj[e].header.update(key=key, value = hwcs[k])
                        if fobj[e].header.has_key('ORIENTAT'):
                            cd12 = 'CD1_2%s' % wcskey
                            cd22 = 'CD2_2%s' % wcskey
                            norient = np.rad2deg(np.arctan2(hwcs[cd12],hwcs[cd22]))
                            fobj[e].header.update(key='ORIENTAT', value=norient)
                else:
                    continue
                
            if not clobber:
                fobj.writeto(name)
            fobj.close()

def archiveWCS(fname, ext, wcskey, wcsname=" ", clobber=False):
    """
    Copy the primary WCS to an alternate WCS 
    with wcskey and name WCSNAME.
    """
    assert len(wcskey) == 1
    if wcskey == " " and clobber==False: 
        print "Please provide a valid wcskey for this WCS."
        print 'Use "utils.next_wcskey" to obtain a valid wcskey.'
        print 'Or use utils.restoreWCS a specific WCS to the primary WCS.'
        print 'WCS was NOT archived.'
        return
    if (wcskey not in available_wcskeys(pyfits.getheader(fname, ext=ext))) and clobber==False:
        print 'Wcskey %s is already used in this header.' % wcskey
        print 'Use "utils.next_wcskey" to obtain a valid wcskey'
        print 'or use "clobber=True" to overwrite the values.'
        return
        
    f = pyfits.open(fname, mode='update')
    w = pywcs.WCS(f[ext].header, fobj=f)
    hwcs = w.to_header()
    wkey = 'WCSNAME' + wcskey
    f[ext].header.update(key=wkey, value=wcsname)
    if w.wcs.has_pc():
        for c in ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2']:
            del hwcs[c]
    elif w.wcs.has_cd():
        for c in ['1_1', '1_2', '2_1', '2_2']:
            hwcs.update(key='CD'+c, value=hwcs['PC'+c])
            del hwcs['PC'+c]
    for k in hwcs.keys():
        key = k+wcskey
        f[ext].header.update(key=key, value = hwcs[k])
    norient = np.rad2deg(np.arctan2(hwcs['CD1_2'],hwcs['CD2_2']))
    okey = 'ORIENT%s' % wcskey
    f[ext].header.update(key=okey, value=norient) 
    f.close()
    

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 wcsnames(header):
    """
    Purpose
    =======
    Return a dictionary of wcskey: WCSNAME pairs
    """
    names = header["WCSNAME*"]
    d = {}
    for c in names:
        d[c.key[-1]] = c.value
    return d
    

def wcskeys(header):
    """
    Purpose
    =======
    Returns a list of characters used in the header for alternate 
    WCS description via WCSNAME keyword
    
    `header`: pyfits.Header
    """
    names = header["WCSNAME*"]
    return [key.split('WCSNAME')[1] for key in names.keys()]

def available_wcskeys(header):
    """
    Purpose
    =======
    Returns a list of characters which are not used in the header
    with WCSNAME keyword. Any of them can be used to save a new
    WCS.
    
    `header`: pyfits.Header
    """
    all_keys = list(string.ascii_uppercase)
    used_keys = wcskeys(header)
    try: 
        used_keys.remove("")
    except ValueError: 
        pass
    [all_keys.remove(key) for key in used_keys]
    return all_keys

def next_wcskey(header):
    """
    Purpose
    =======
    Returns next available character to be used for an alternate WCS
    
    `header`: pyfits.Header    
    """
    allkeys = available_wcskeys(header)
    if allkeys != []:
        return allkeys[0]
    else:
        return None