forked from abostroem/13287
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cross_correlate_09ip.py
182 lines (159 loc) · 7.27 KB
/
cross_correlate_09ip.py
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
from astropy.io import fits
import numpy as np
from matplotlib import pyplot
from scipy.ndimage.interpolation import rotate
import time
from matplotlib.backends.backend_pdf import PdfPages
import pdb
import scipy.interpolate
import scipy.ndimage
def congrid(a, newdims, method='linear', centre=False, minusone=False):
if not a.dtype in [np.float64, np.float32]:
a = np.cast[float](a)
m1 = np.cast[int](1)
print m1
ofs = np.cast[int](centre) * 0.5
old = np.array( a.shape )
ndims = len( a.shape)
#if len( newdims ) != ndims:
# print "[congrid] dimensions error. "\
# "This routine currently only supports "\
# "rebinning to the same number of dimensions."
# return None
newdims = np.asarray( newdims, dtype=float )
#newdims = newdims.reshape(newdims,1)
dimlist = []
if method == 'neighbor':
for i in range( ndims ):
base = np.indices(newdims)
dimlist.append( (old - m1) / (newdims - m1) * (base + ofs) - ofs )
cd = np.array( dimlist ).round().astype(int)
newa = a[list( cd )]
return newa
elif method in ['nearest','linear']:
for i in range( ndims ):
base = np.arange( newdims )
dimlist.append( (old - m1) / (newdims - m1) * (base + ofs) - ofs)
olddims = [np.arange(i, dtype = np.float) for i in list( a.shape )]
mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method )
newa = mint
newa = mint( dimlist[-1] )
trorder = [ndims - 1] + range (ndims - 1 )
for i in range( ndims -2, -1, -1 ):
print 'here'
newa = newa.transpose( trorder)
mint = scipy.interpolate.interp1d( olddims[i], newa, kind=method )
newa = mint( dimlist[i] )
if ndims > 1:
newa = newa.transpose( trorder)
return newa
else:
print "Congrid error:"
return None
def cross_correlate_05ip(input_stis_file):
'''
A WFC3 image of 2010jl's host galaxy is rotated to the STIS orientation.
This code steps across a small region of a WFC3 image on the 2010jl galaxy and collapsed
every 2 columns to create a cross-dispersion profile the same width as the STIS
slit (0.2 arcseconds). It then cross correlates the STIS FUV cross dispersion profile
with the WFC3 cross-dispersion profile to find the slit position against the background.
The best cross-correlation for each WFC3 x-location is plotted in the file
xd_profile_match.pdf.
'''
pp = PdfPages('xd_profile_match_09ip.pdf')
orientat = 4.17551#+45.11698363500195 #(WFPC2+STIS/CCD)
#read in and rotate the WFC3 image
wfpc2_img = fits.getdata('2009ip_otfr/wfpc2/hst_06359_40_wfpc2_f606w_wf_drz.fits', 1) #F814W
rot_wfpc2_img = rotate(wfpc2_img, orientat)
#Make a STIS cross-dispersion profile
img = fits.getdata(input_stis_file, 1)
#xd_profile = np.sum(img, axis = 1)
xd_profile = np.sum(img[:, 200:900], axis = 1)
#Remove occulation bar (set values to the median of the profile)
#xd_profile[854:910] = np.median(xd_profile)
#Bin STIS to WFC3 resolution
stis_plate_scale = 0.05071 #/arcsec/pix (Make sure this is STIS MAMA plate scale)
wfc3_platescale = 0.0394 #arcsec/pix (Make sure this is wfpc2 platescale)
# binning_factor = 2.#*wfc3_platescale/stis_plate_scale
# binned_xd_profile = xd_profile[::int(binning_factor)]
# for indx in range(1,int(binning_factor)):
# binned_xd_profile = binned_xd_profile + xd_profile[indx::int(binning_factor)]
#Calculate the STIS slit height in WFC3 profile
stis_slit_height_pix = 25.0/wfc3_platescale
#0.2 arcsec slit = 2 WFPC2 pixels
stis_slit_width_pix = 0.2/wfc3_platescale
#columns = np.arange(1210, 1260, 2) #Determine which WFC3 columns to cross-correlate
wfpc2_x_start = 1200 #2480, 2520
wfpc2_x_end = 1400
columns = np.arange(wfpc2_x_start, wfpc2_x_end, 200) #Determine which WFC3 columns to cross-correlate stis_slit_width_pix
corr = []
offset = []
wfpc2_y_start = 1000 #bottom of the STIS slit
multfactor = 1.
#xd_profile = xd_profile.reshape(len(xd_profile),1)
binned_xd_profile = congrid(xd_profile, multfactor*int(stis_slit_height_pix))
fig = pyplot.figure(figsize = [20, 15])
ax1 = fig.add_subplot(1, 3, 1)
ax1.set_title('STIS Slit position on WFC3 Image')
im1 = ax1.imshow(rot_wfpc2_img, interpolation = 'nearest', vmin = 0, vmax = 0.1)
ax1.set_ylim(wfpc2_y_start+0, wfpc2_y_start+0+1*stis_slit_height_pix)
ax1.set_xlim(wfpc2_x_start, wfpc2_x_end)
ax1.set_yticks(np.arange(ax1.get_ylim()[0], ax1.get_ylim()[1], 50))
ax1.grid(color = 'w')
#make WFC3 xd profile 2 times actual profile to account for offsetting
for column in columns:
#0.2 arcsec slit = 2 WFPC2 pixels
wfpc2_profile = np.sum(rot_wfpc2_img[wfpc2_y_start: wfpc2_y_start + 1*stis_slit_height_pix, column:column+stis_slit_width_pix], axis = 1)
binned_wfpc2_profile=congrid(wfpc2_profile, multfactor*int(stis_slit_height_pix))
wfpc2_profile=binned_wfpc2_profile
#Normalize the 2 XD profiles
#Correlate bigger array, smaller array apply offset to smaller array
correlation_array = np.correlate(wfpc2_profile/np.max(wfpc2_profile), binned_xd_profile/np.max(binned_xd_profile))
max_corr_indx = 160
max_corr = 0
#max_corr_indx = np.argmax(correlation_array)
#max_corr = correlation_array[max_corr_indx]
corr.append(max_corr)
offset.append(max_corr_indx)
#Plot the XD profile for each column tested
fig = pyplot.figure(figsize = [7, 25])
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_title('STIS Slit position on WFC3 Image')
im1 = ax1.imshow(rot_wfpc2_img, interpolation = 'nearest', vmin = 0, vmax = 5)
ax1.set_ylim(wfpc2_y_start+max_corr_indx, wfpc2_y_start+max_corr_indx+1*stis_slit_height_pix)
ax1.set_xlim(wfpc2_x_start, wfpc2_x_end)
ax1.axvspan(wfpc2_x_start, column, color = 'k', alpha = 0.5)
ax1.axvspan(column+stis_slit_width_pix, wfpc2_x_end, color = 'k', alpha = 0.1)
ax1.set_yticks(np.arange(ax1.get_ylim()[0], ax1.get_ylim()[1], 50))
ax1.grid(color = 'w')
#fig_temp = pyplot.figure(figsize = [7, 25])
ax_temp = fig.add_subplot(1, 2, 2)
ax_temp.plot(binned_wfpc2_profile/max(binned_wfpc2_profile), np.arange(len(binned_wfpc2_profile)))
ax_temp.plot(binned_xd_profile/max(binned_xd_profile), np.arange(len(binned_xd_profile))+max_corr_indx)
print len(binned_xd_profile)
ax_temp.legend(['Norm WFC3 profile', 'Norm, binned, STIS profile'], loc = 'best')
ax_temp.set_title('WFC3 Column {}'.format(column))
#ax_temp.set_xlim(0, 0.5)
ax_temp.set_ylim(max_corr_indx, len(binned_xd_profile)+max_corr_indx)
pp.savefig()
pyplot.close()
pp.close()
corr = np.array(corr)
fig2 = pyplot.figure()
ax2 = fig2.add_subplot(1,1,1)
pyplot.plot(columns, corr)
final_max_corr_indx = np.argmax(corr)
fig3 = pyplot.figure()
ax3 = fig3.add_subplot(1,1,1)
ax3.plot(binned_xd_profile/binned_xd_profile.max(), np.arange(len(binned_xd_profile))+offset[final_max_corr_indx])
best_col = rot_wfpc2_img[wfpc2_y_start : wfpc2_y_start+ 2.0*stis_slit_height_pix, columns[final_max_corr_indx]]
ax3.plot(best_col/best_col.max(), range(len(best_col)))
#pdb.set_trace()
print 'highest correlation value = ', np.max(corr)
print 'Best offset value = ', offset[final_max_corr_indx]
print 'Best WFPC2 column = ', columns[final_max_corr_indx]
if __name__ == "__main__":
#cross_correlate_wfc3_stis_fuv('otfr_data2/ocdd03090_flt.fits') #718, 132
#cross_correlate_wfc3_stis_fuv('otfr_data2/ocdd030d0_flt.fits') #718, 61
#cross_correlate_wfc3_stis_ccd('2006gy/ocdd04010_flt.fits')
cross_correlate_05ip('2009ip_otfr/ocdd01070_flt.fits')