-
Notifications
You must be signed in to change notification settings - Fork 40
/
pixelCorrection.py
209 lines (176 loc) · 6.34 KB
/
pixelCorrection.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
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
import sys, os
import numpy as np
import matplotlib.pyplot as plt
from lightkurve import KeplerTargetPixelFile as ktpf
from scipy import ndimage
import collections
from astropy.wcs import WCS
from astropy.io import fits
from gaiaTIC import coneSearch, jsonTable
from gaiaTIC import ticSearchByContam as tsbc
from scipy.optimize import minimize
def openFITS(dir, fn):
""" Opens FITS file to get mast, mheader """
return fits.getdata(dir+fn, header=True)
def radec2pixel(header, r, contam):
""" Completes a cone search around center of FITS file """
pos = [header['CRVAL1'], header['CRVAL2']]
print("Getting sources")
data = tsbc(pos, r, contam)
dataTable = jsonTable(data)
ra, dec = dataTable['ra'], dataTable['dec']
return WCS(header).all_world2pix(ra, dec, 1), dataTable['ID']
#################
##### ADDED #####
#################
def sortByDate(fns, dir):
"""
Sorts FITS files by start date of observation
Parameters
----------
dir: directory the files are in
Returns
----------
fns: sorted filenames by date
"""
dates = []
for f in fns:
mast, header = fits.getdata(dir+f, header=True)
dates.append(header['DATE-OBS'])
dates, fns = np.sort(np.array([dates, fns]))
return fns
#################
##### ADDED #####
#################
def nearest(x, y, x_list, y_list):
"""
Calculates the distance to the nearest source
Parameters
----------
x: x-coordinate of a source
y: y-coordinate of a source
x_list: x-coordinates for all sources in file
y_list: y-coordinates for all sources in file
Returns
----------
dist: distance to closest source
closest: index of closest source
"""
x_list = np.delete(x_list, np.where(x_list==x))
y_list = np.delete(y_list, np.where(y_list==y))
closest = np.sqrt( (x-x_list)**2 + (y-y_list)**2 ).argmin()
dist = np.sqrt( (x-x_list[closest])**2 + (y-y_list[closest])**2 )
return dist, closest
#################
##### ADDED #####
#################
def findIsolated(x, y):
"""
Finds isolated sources in the image where an isolated source is >= 5 pixels away from
any other source
Parameters
----------
x: a list of x-coordinates for all sources in file
y: a list of y-coordinates for all sources in file
Returns
----------
isolated: a list of isolated source indices
"""
isolated = []
for i in range(len(x)):
dist, dist_ind = nearest(x[i], y[i], x, y)
if dist >= 8.0:
isolated.append(i)
return isolated
def calcShift(dir, fns, x, y, corrFile, mast):
"""
Calculates the deltaX, deltaY, and rotation for isolated sources in order to
put together a pointing model for the entire chip
Parameters
----------
fns: all FITS file names
x: a list of x-coordinates for isolated sources
y: a list of y-coordinates for isolated sources
corrFile: name of file to write corrections to for each cadence
Returns
----------
"""
def model(params):
nonlocal xy, centroids
theta, xT, yT = params
theta = np.radians(theta)
xRot = xy[:,0]*np.cos(theta) - xy[:,1]*np.sin(theta)
yRot = xy[:,0]*np.sin(theta) + xy[:,1]*np.cos(theta)
xRot += xT
yRot += yT
dist = np.sqrt((xRot-centroids[:,0])**2 + (yRot-centroids[:,1])**2)
dist = np.square(dist)
return np.sum(dist)
# rand = np.random.randint(0, len(x), size=200)
# x, y = x[rand], y[rand]
print(len(x))
fns = [dir+i for i in fns]
x_cen, y_cen = len(mast)/2, len(mast[0])/2
xy = np.zeros((len(x),2))
for i in range(len(x)):
xy[i][0] = x[i]-x_cen
xy[i][1] = y[i]-y_cen
with open(corrFile, 'w') as tf:
tf.write('cadence medT medX medY\n')
matrix = np.zeros((len(x), len(fns), 2))
for i in range(len(x)):
print(x[i], y[i])
tpf = ktpf.from_fits_images(images=fns, position=(x[i], y[i]), size=(6,6))
#tpf.to_fits(output_fn='test{}.fits'.format(i))
for j in range(len(fns)):
com = ndimage.measurements.center_of_mass(tpf.flux[j].T-np.median(tpf.flux[j])) # subtracts background
matrix[i][j][0] = com[0]+xy[i][0]
matrix[i][j][1] = com[1]+xy[i][1]
for i in range(len(fns)):
centroids = matrix[:,i]
if i == 0:
initGuess = [0.001, 0.1, -0.1]
else:
initGuess = solution.x
bnds = ((-0.08, 0.08), (-5.0, 5.0), (-5.0, 5.0))
solution = minimize(model, initGuess, method='L-BFGS-B', bounds=bnds, options={'ftol':5e-11,
'gtol':5e-05})
sol = solution.x
with open(corrFile, 'a') as tf:
if i == 0:
theta, delX, delY = sol[0], sol[1], sol[2]
else:
theta = oldSol[0] - sol[0]
delX = oldSol[1] - sol[1]
delY = oldSol[2] - sol[2]
tf.write('{}\n'.format(str(i) + ' ' + str(theta) + ' ' + str(delX) + ' ' + str(delY)))
oldSol = sol
return
def correctionFactors(camera, chip, dir):
"""
Creates a pointing model for a given camera and chip
Parameters
----------
camera: which camera to create a pointing model for
chip: which chip on said camera to create a pointing model for
dir: directory in which the FITS files can be found
Returns
----------
"""
filenames = np.array(os.listdir(dir))
fitsInds = np.array([i for i,item in enumerate(filenames) if "fits" in item])
filenames = filenames[fitsInds]
filenames = sortByDate(filenames, dir)
mast, header = openFITS(dir, filenames[0])
xy, id = radec2pixel(header, 6*np.sqrt(2), [0.0, 5e-5])
x, y = xy[0], xy[1]
inds = np.where((y>=44.) & (y<len(mast)-45.) & (x>=0.) & (x<len(mast[0])-41.))[0]
x, y = x[inds], y[inds]
isolated = findIsolated(x, y)
print(len(isolated))
x, y = x[isolated], y[isolated]
calcShift(dir, filenames, x, y, 'pointingModel_{}-{}.txt'.format(camera, chip), mast)
#correctionFactors(3, 1, './2019/2019_1_3-1/ffis/')
#correctionFactors(4, 4, './2019/2019_1_4-4/ffis/')
#correctionFactors(3, 3, './2019/2019_1_3-3/ffis/')
#correctionFactors(1, 3, './2019/2019_1_1-3/ffis/')