summaryrefslogtreecommitdiff
path: root/wcsutil/__init__.py
blob: b6820321ddb35a95fb356581d2dc83e78cd3a481 (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
229
230
231
232
233
234
235
236
237
238
239
240
#from .. pywcs.sip import SIP
from pywcs import WCS
import pyfits
import instruments
#from .. distortion import models
from updatewcs.distortion import models
import numpy as N
from pytools import fileutil
from pytools.fileutil import DEGTORAD, RADTODEG

#from .. mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG, basic_wcs
import mappings
from mappings import inst_mappings, ins_spec_kw
from mappings import basic_wcs, prim_hdr_kw

__docformat__ = 'restructuredtext'
__version__ = '0.3'


class HSTWCS(WCS):
    """
    Purpose
    =======
    Create a WCS object based on the instrument.
    It has all basic WCS kw as attribbutes (set by pywcs).
    It also uses the primary and extension header to define 
    instrument specific attributes needed by the correction classes.
    """
    
    def __init__(self, hdr0=None, ehdr=None, fobj=None, fname=None, instrument=None):
        """
        :Parameters:
        `fname`: string
                file/extension name
                filename[EXTNAME,EXTNUM]
                filename[extension_number]
        `hdr0`: Pyfits Header
                primary header
        `ehdr`: Pyfits Header
               extension header
        `fobj`: PyFITS HDUList object or None
                pyfits file object
        `instrument`: string
                one of 'ACS', 'NICMOS', 'WFPC2', 'STIS', 'WFC3' 
        """
        
        self.inst_kw = ins_spec_kw
        
        if fname == None and hdr0 == None and ehdr == None and instrument != None:
            # create a default HSTWCS object based on instrument only
            self.instrument = instrument
            self.setInstrSpecKw()
        elif fname != None:
            # create an HSTWCS object from a filename
            filename, extname = fileutil.parseFilename(fname)
            self.filename = filename
            self.extname = extname
            if extname == None:
                #data may be in the primary array
                ehdr = pyfits.getheader(filename)
                hdr0 = None
            else:
                ext = fileutil.parseExtn(extname)
                hdr0 = pyfits.getheader(filename)
                try:
                    ehdr = pyfits.getheader(filename, ext=ext)
                except:
                    print 'Unable to get extension header based on filename %s. /n' % fname
        elif fname == None and  instrument == None:
            # hdr0 may be None, a WCS object will still be created
            if hdr0 == None and ehdr == None:
                print 'Not enough information to create a WCS object\n'
                self.help()
                return
            if hdr0 != None:
                assert isinstance (hdr0, pyfits.Header)
            if ehdr != None:
                assert isinstance (ehdr, pyfits.Header)
                
        WCS.__init__(self, ehdr, fobj=fobj)    
        self.setHDR0kw(hdr0, ehdr)
        self.setInstrSpecKw(hdr0, ehdr)
        self.pscale = self.setPscale()
        self.orientat = self.setOrient()
        self.readIDCCoeffs(ehdr)
        # if input was not a file name, try to get it from the primary header
        self.filename = hdr0.get('FILENAME', "")
        extname = ehdr.get('EXTNAME', "")
        extnum = ehdr.get('EXTVER', None)
        self.extname = (extname, extnum)
            
        
        
    def setHDR0kw(self, primhdr, ehdr):
        # Set attributes from kw defined in the primary header.
        self.instrument = primhdr.get('INSTRUME', None)
        self.offtab = primhdr.get('OFFTAB', None) 
        self.idctab = primhdr.get('IDCTAB', None)
        self.date_obs = primhdr.get('DATE-OBS', None)
        self.pav3 = primhdr.get('PA_V3', None)
        self.ra_targ = primhdr.get('RA_TARG', None)
        self.dec_targ = primhdr.get('DEC_TARG', None)
        self.filename = primhdr.get('FILENAME', "")
        self.detector = primhdr.get('DETECTOR', None)
    
    
    def readIDCCoeffs(self, header):
        """
        Reads in first order IDCTAB coefficients if present in the header
        """
        coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale']
        for c in coeffs:
            self.__setattr__(c, header.get(c, None))
            
    def setInstrSpecKw(self, prim_hdr=None, ext_hdr=None):
        # Based on the instrument kw creates an instalnce of an instrument WCS class
        # and sets attributes from instrument specific kw 
        if self.instrument in inst_mappings.keys():
            inst_kl = inst_mappings[self.instrument]
            inst_kl = instruments.__dict__[inst_kl]
            insobj = inst_kl(prim_hdr, ext_hdr)
            for key in self.inst_kw:
                try:
                    self.__setattr__(key, insobj.__getattribute__(key))
                except AttributeError:
                    # Some of the instrument's attributes are recorded in the primary header and
                    # were already set, (e.g. 'DETECTOR'), the code below is a check for that case.
                    if not self.__getattribute__(key):
                        raise
                    else:
                        pass
        else:
            raise KeyError, "Unsupported instrument - %s" %self.instrument
    
    def setPscale(self):
        # Calculates the plate scale from the cd matrix
        # this may not be needed now and shoufl probably be done after all 
        # corrections
        cd11 = self.wcs.cd[0][0]
        cd21 = self.wcs.cd[1][0]
        return N.sqrt(N.power(cd11,2)+N.power(cd21,2)) * 3600.
    
    def setOrient(self):
        # Recompute ORIENTAT
        cd12 = self.wcs.cd[0][1]
        cd22 = self.wcs.cd[1][1]
        return RADTODEG(N.arctan2(cd12,cd22))
    
    def updatePscale(self, pscale):
        """Given a plate scale, update the CD matrix"""
        old_pscale = self.pscale
        self.pscale = pscale
        self.wcs.cd = self.wcs.cd * pscale/old_pscale
        self.naxis1 = self.naxis1 * old_pscale/ pscale 
        self.naxis2 = self.naxis2 * old_pscale/ pscale 
        self.wcs.crpix = self.wcs.crpix *old_pscale/pscale
        
    def updateOrient(self, orient):
        """Given n angle update the CD matrix"""
        if self.orientat == orient:
            return
        old_orient = self.orientat
        self.orientat = orient
        angle = fileutil.DEGTORAD(orient)
        cd11 = -N.cos(angle)
        cd12 = N.sin(angle)
        cd21 = cd12
        cd22 = -cd11
        cdmat = N.array([[cd11, cd12],[cd21,cd22]])
        self.wcs.cd = cdmat * self.pscale/3600
        
    def verifyPa_V3(self):
        """make sure PA_V3 is populated"""    
        
    def verifyKw(self):
        """verify that all required kw have meaningful values"""
            
    def readModel(self, update=False, header=None):
        """
        Purpose
        =======
        Read distortion model from idc table.
        Save some of the information as kw needed for interpreting the distortion
        If header is provided and update is true, some IDC model kw 
        will be recorded in the header.
        """
        if self.idctab == None or self.date_obs == None:
            print 'idctab or date_obs not available\n'
            self.idcmodel = None
            return
        if self.filter1 ==None and self.filter2 == None:
            'No filter information available\n'
            self.idcmodel = None
            return
        
        self.idcmodel = models.IDCModel(self.idctab,
                    chip=self.chip, direction='forward', date=self.date_obs,
                    filter1=self.filter1, filter2=self.filter2, 
                    offtab=self.offtab, binned=self.binned)
        
        if update:
            if header==None:
                print 'Update header with IDC model kw requested but header was not provided\n.'
                return
            else:
                self.updatehdr(header)
                   
    
    def updatehdr(self, ext_hdr, newkeywords=None):
        #kw2add : OCX10, OCX11, OCY10, OCY11 
        # record the model in the header for use by pydrizzle
        ext_hdr.update('OCX10', self.idcmodel.cx[1,0])
        ext_hdr.update('OCX11', self.idcmodel.cx[1,1])
        ext_hdr.update('OCY10', self.idcmodel.cy[1,0])
        ext_hdr.update('OCY11', self.idcmodel.cy[1,1])
        ext_hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE'])
        ext_hdr.update('IDCTHETA', self.idcmodel.refpix['THETA'])
        ext_hdr.update('IDCXREF', self.idcmodel.refpix['XREF'])
        ext_hdr.update('IDCYREF', self.idcmodel.refpix['YREF'])
        ext_hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF'])
        ext_hdr.update('IDCV3REF', self.idcmodel.refpix['V3REF'])
        
    
    def help(self):
        print 'How to create an HSTWCS object:\n\n'
        print """ \
        1. Create an HSTWCS object using pyfits.Header objects. \n
        Example:\n
        hdr0 = pyfits.getheader('j9irw4b1q_flt.fits', ext=0)\n
        hdr1 = pyfits.getheader('j9irw4b1q_flt.fits', ext=1)\n
        f = pyfits.open('j9irw4b1q_flt.fits')\n
        f is required only if a lookup table distortion is available.\n
        w = wcsutil.HSTWCS(hdr0, hdr1,f)\n\n
        
        2. Create an HSTWCS object using a qualified file name. \n
        Example:\n
        w = wcsutil.HSTWCS('j9irw4b1q_flt.fits[sci,1]')\n\n
        """