-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgal_trace.py
114 lines (93 loc) · 3.53 KB
/
gal_trace.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
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits as pyfits
from scipy.optimize import curve_fit
def _gaus(x,a,x0):
return a*np.exp(-(x-x0)**2/(2*5.0**2))
def _quadfit(x,a,b,c):
'''define quadratic galaxy fitting function'''
return a*(x-2032)**2 + b*(x-2032) + c
def gal_trace(cutout):
means = []
meansx = []
#define 40 'pins' to fit quadratic function to
for i in range(40):
i += 1
x = np.arange(0,len(cutout),1)
#stack the flux values for the nearest 100 xpixels and take the median along
#the xpixels to identify galaxy light
y = np.median(cutout[:,(100*i)-50:(100*i)+50],axis=1)
ymg = 25.0 #approx value (guess) for galaxy
try:
popt,pcov = curve_fit(_gaus,x,y,p0=[1.0,ymg]) #fit gaussian to points
except RuntimeError:
popt = [np.nan,np.nan]
pcov = [[np.inf,np.inf],[np.inf,np.inf]]
#if not a good fit, fill with nan values
if np.sqrt(np.diag(pcov))[1] < 5.0:
means.append(popt[1])
else:
means.append(np.nan)
meansx.append(100*i)
means.insert(0,means[0])
meansx.insert(0,0.0)
means.append(means[-1])
meansx.append(4064.0)
means = np.array(means)
meansx = np.array(meansx)
if means[np.isfinite(means)].size >= 10: #otherwise don't trust the fit and return the entire cutout summed
#fit quadratic equation to estimated galaxy position
popt,pcov = curve_fit(_quadfit,meansx[np.isfinite(means)],means[np.isfinite(means)],p0=[-1e-6,1e-6,20])
#plt.plot(meansx,means,'ko')
#plt.plot(meansx,_quadfit(np.array(meansx),*popt))
#plt.show()
gal_pos = _quadfit(np.arange(0,cutout.shape[1],1),*popt)
gal_posr = np.round(np.array(gal_pos))
gal_sum = []
for i in range(cutout.shape[1]):
gal_sum.append(np.sum(cutout[gal_posr[i]-5:gal_posr[i]+5,i]))
gal_sum = np.array(gal_sum)
else:
gal_sum = np.sum(cutout,axis=0)
return gal_sum
if __name__ == '__main__':
from astropy.io import fits as pyfits
fits = pyfits.open('C4_0199/science/C4_0199_science.reduced.fits')
img = fits[0].data
plt.imshow(img,vmin=-4,vmax=14)
plt.show()
mini = np.round(1478.0)
maxi = np.round(1525.0)
cutout = img[mini:maxi]
plt.imshow(cutout,vmin=-4,vmax=14,aspect='auto')
plt.show()
sigma = 5.0
def gaus(x,a,x0):
return a*np.exp(-(x-x0)**2/(2*sigma**2))
means = []
meansx = []
for i in range(40):
i += 1
x = np.arange(0,len(cutout),1)
y = np.median(cutout[:,(100*i)-50:(100*i)+50],axis=1)
ymg = 25.0
popt,pcov = curve_fit(gaus,x,y,p0=[1.0,ymg])
print np.sqrt(np.diag(pcov))
if i % 1.0 == 0:
plt.plot(np.arange(0,len(cutout),1),np.median(cutout[:,100*i-50:100*i+50],axis=1),'ko',markersize=4)
plt.plot(x,gaus(x,*popt),'r',ls='--',alpha=0.6)
plt.show()
if np.sqrt(np.diag(pcov))[1] < 5.0:
means.append(popt[1])
else:
means.append(np.nan)
meansx.append(100*i)
means.insert(0,means[0])
meansx.insert(0,0.0)
means.append(means[-1])
meansx.append(4064.0)
#define quadratic galaxy fitting function
def quadfit(x,a,b,c):
return a*(x-2032)**2 + b*(x-2032) + c
popt,pcov = curve_fit(quadfit,np.array(meansx)[np.isfinite(means)],np.array(means)[np.isfinite(means)],p0=[-1e-6,1e-6,35])
#gal_spec = gal_trace(cutout)