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
|
from __future__ import absolute_import, print_function
from astropy import wcs as pywcs
from collections import OrderedDict
from astropy.io import fits
from .headerlet import parse_filename
import numpy as np
def is_wcs_identical(scifile, file2, sciextlist, fextlist, scikey=" ",
file2key=" ", verbose=False):
"""
Compares the WCS solution of 2 files.
Parameters
----------
scifile: string
name of file1 (usually science file)
IRAF style extension syntax is accepted as well
for example scifile[1] or scifile[sci,1]
file2: string
name of second file (for example headerlet)
sciextlist - list
a list of int or tuple ('SCI', 1), extensions in the first file
fextlist - list
a list of int or tuple ('SIPWCS', 1), extensions in the second file
scikey: string
alternate WCS key in scifile
file2key: string
alternate WCS key in file2
verbose: boolean
True: print to stdout
Notes
-----
These can be 2 science observations or 2 headerlets
or a science observation and a headerlet. The two files
have the same WCS solution if the following are the same:
- rootname/destim
- primary WCS
- SIP coefficients
- NPOL distortion
- D2IM correction
"""
result = True
diff = OrderedDict()
fobj, fname, close_file = parse_filename(file2)
sciobj, sciname, close_scifile = parse_filename(scifile)
diff['file_names'] = [scifile, file2]
if get_rootname(scifile) != get_rootname(file2):
# logger.info('Rootnames do not match.')
diff['rootname'] = ("%s: %s", "%s: %s") % (sciname, get_rootname(scifile), file2,
get_rootname(file2))
result = False
for i, j in zip(sciextlist, fextlist):
w1 = pywcs.WCS(sciobj[i].header, sciobj, key=scikey)
w2 = pywcs.WCS(fobj[j].header, fobj, key=file2key)
diff['extension'] = [get_extname_extnum(sciobj[i]), get_extname_extnum(fobj[j])]
if not np.allclose(w1.wcs.crval, w2.wcs.crval, rtol=10**(-7)):
# logger.info('CRVALs do not match')
diff['CRVAL'] = w1.wcs.crval, w2.wcs.crval
result = False
if not np.allclose(w1.wcs.crpix, w2.wcs.crpix, rtol=10**(-7)):
# logger.info('CRPIX do not match')
diff['CRPIX'] = w1.wcs.crpix, w2.wcs.crpix
result = False
if not np.allclose(w1.wcs.cd, w2.wcs.cd, rtol=10**(-7)):
# logger.info('CDs do not match')
diff['CD'] = w1.wcs.cd, w2.wcs.cd
result = False
if not (np.array(w1.wcs.ctype) == np.array(w2.wcs.ctype)).all():
# logger.info('CTYPEs do not match')
diff['CTYPE'] = w1.wcs.ctype, w2.wcs.ctype
result = False
if w1.sip or w2.sip:
if (w2.sip and not w1.sip) or (w1.sip and not w2.sip):
diff['sip'] = 'one sip extension is missing'
result = False
if not np.allclose(w1.sip.a, w2.sip.a, rtol=10**(-7)):
diff['SIP_A'] = 'SIP_A differ'
result = False
if not np.allclose(w1.sip.b, w2.sip.b, rtol=10**(-7)):
# logger.info('SIP coefficients do not match')
diff['SIP_B'] = (w1.sip.b, w2.sip.b)
result = False
if w1.cpdis1 or w2.cpdis1:
if w1.cpdis1 and not w2.cpdis1 or w2.cpdis1 and not w1.cpdis1:
diff['CPDIS1'] = "CPDIS1 missing"
result = False
if w1.cpdis2 and not w2.cpdis2 or w2.cpdis2 and not w1.cpdis2:
diff['CPDIS2'] = "CPDIS2 missing"
result = False
if not np.allclose(w1.cpdis1.data, w2.cpdis1.data, rtol=10**(-7)):
# logger.info('NPOL distortions do not match')
diff['CPDIS1_data'] = (w1.cpdis1.data, w2.cpdis1.data)
result = False
if not np.allclose(w1.cpdis2.data, w2.cpdis2.data, rtol=10**(-7)):
# logger.info('NPOL distortions do not match')
diff['CPDIS2_data'] = (w1.cpdis2.data, w2.cpdis2.data)
result = False
if w1.det2im1 or w2.det2im1:
if w1.det2im1 and not w2.det2im1 or \
w2.det2im1 and not w1.det2im1:
diff['DET2IM'] = "Det2im1 missing"
result = False
if not np.allclose(w1.det2im1.data, w2.det2im1.data, rtol=10**(-7)):
# logger.info('Det2Im corrections do not match')
diff['D2IM1_data'] = (w1.det2im1.data, w2.det2im1.data)
result = False
if w1.det2im2 or w2.det2im2:
if w1.det2im2 and not w2.det2im2 or \
w2.det2im2 and not w1.det2im2:
diff['DET2IM2'] = "Det2im2 missing"
result = False
if not np.allclose(w1.det2im2.data, w2.det2im2.data, rtol=10**(-7)):
# logger.info('Det2Im corrections do not match')
diff['D2IM2_data'] = (w1.det2im2.data, w2.det2im2.data)
result = False
if not result and verbose:
for key in diff:
print(key, ":\t", diff[key][0], "\t", diff[key][1])
if close_file:
fobj.close()
if close_scifile:
sciobj.close()
return result, diff
def get_rootname(fname):
"""
Returns the value of ROOTNAME or DESTIM
"""
hdr = fits.getheader(fname)
try:
rootname = hdr['ROOTNAME']
except KeyError:
try:
rootname = hdr['DESTIM']
except KeyError:
rootname = fname
return rootname
def get_extname_extnum(ext):
"""
Return (EXTNAME, EXTNUM) of a FITS extension
"""
extname = ""
extnum = 1
extname = ext.header.get('EXTNAME', extname)
extnum = ext.header.get('EXTVER', extnum)
return (extname, extnum)
|
