-
Notifications
You must be signed in to change notification settings - Fork 0
/
starfinder.py
executable file
·217 lines (208 loc) · 8.35 KB
/
starfinder.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
210
211
212
213
214
215
216
217
#!/usr/bin/python
'''
Program:
This is a program to find stars with IRAFstarfinder
Then save a target on frames table
Usage:
starfinder.py [image_list]
Editor:
Jacob975
20180626
#################################
update log
20180626 version alpha 1
1. The code works
20181029 version alpha 2
1. Make the program simplier, leave the func for finding source only.
Anything else(IDP, EP) is independant from here.
'''
from photutils.detection import IRAFStarFinder, DAOStarFinder
from photutils import aperture_photometry, CircularAperture
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.io import fits as pyfits
from astropy import wcs
from astropy.coordinates import SkyCoord
from astropy import units as u
from fit_lib import get_star, get_peak_filter, get_flux
from reduction_lib import image_info
from mysqlio_lib import save2sql, find_fileID, find_source_match_coords
import matplotlib.pyplot as plt
import numpy as np
import time
import os
from sys import argv
import TAT_env
from test_EP import flux2mag
# find a star through iraf finder
def starfinder(image_name):
infos = image_info(image_name)
mean_bkg = infos.bkg
std_bkg = infos.std
u_sigma = infos.u_sigma
sigma = u_sigma.n
print 3.0*std_bkg + mean_bkg
iraffind = IRAFStarFinder(threshold = 3.0*std_bkg + mean_bkg, \
# fwhm = sigma*gaussian_sigma_to_fwhm, \
fwhm = 5.0,
roundhi = 1.0,
roundlo = -1.0,
sharphi = 2.0,
sharplo = 0.5)
iraf_table = iraffind.find_stars(infos.data)
return iraf_table, infos
def star_extend(image_name, SE_table):
#------------------------------------------------------------
# Add infomations into the extend star array
# Get WCS infos with astrometry
try:
header_wcs = pyfits.getheader("stacked_image.wcs")
except:
print "WCS not found"
return 1, None, None
w = wcs.WCS(header_wcs)
# Convert pixel coord to RA and DEC
pixcrd = np.transpose(np.array([SE_table[:,2], SE_table[:,3]]))
world = w.wcs_pix2world(pixcrd, 1)
extend_star_array = np.full((len(SE_table), len(TAT_env.obs_data_titles)), None, dtype = object)
# RA and DEC
extend_star_array[:,14] = world[:,0]
extend_star_array[:,15] = world[:,1]
# Name targets with RA and DEC, and insert into table
table_length = len(world)
target_names_list = np.array(["target_{0:.4f}_{1:.4f}".format(world[i,0], world[i,1]) for i in range(table_length)])
# NAME
extend_star_array[:,1] = target_names_list
# FLUX
extend_star_array[:,4] = SE_table[:,6]
extend_star_array[:,5] = SE_table[:,7]
# INST MAG
extend_star_array[:,6] = SE_table[:,4]
extend_star_array[:,7] = SE_table[:,5]
# CENTER
extend_star_array[:,19] = SE_table[:,2]
extend_star_array[:,20] = SE_table[:,3]
# Load infos from the header of the image
header = pyfits.getheader(image_name)
# get the fileID from mysql.
fileID = find_fileID(image_name)
# The else
extend_star_array[:, 3] = header['BJD']
extend_star_array[:, 30] = header['JD']
extend_star_array[:, 31] = header['MJD-OBS']
extend_star_array[:, 32] = header['HJD']
extend_star_array[:, 33] = fileID
return 0, extend_star_array
def SExtractor(image_name):
# Initilized
config_name = "{0}.config".format(image_name[:-5])
catalog_name = "{0}.cat".format(image_name[:-5])
index_flag = TAT_env.SE_table_titles.index('FLAGS')
# Copy parameter files to working directory.
os.system('cp {0}/SE_workshop/SE.config {1}'.format(TAT_env.path_of_code, config_name))
os.system("sed -i 's/stack_image/{0}/g' {1}".format(image_name[:-5], config_name))
# Execute SExtrctor
os.system('sex {0} -c {1}'.format(image_name, config_name))
# Load the result table
table = np.loadtxt(catalog_name)
try:
table = table[table[:,index_flag ] == 0 ]
except:
table = np.array([])
return table
# check if there is new sources.
def check_new_sources(extend_star_array):
# Initialize
new_source_list = []
new = False
index_of_name = TAT_env.obs_data_titles.index('NAME')
index_of_ra = TAT_env.obs_data_titles.index('RA')
index_of_dec = TAT_env.obs_data_titles.index('`DEC`')
# Setup the spatial tolerance.
tolerance = TAT_env.pix1/3600.0 * 3.0
# Match positions
for i in xrange(len(extend_star_array)):
src_coord_list = find_source_match_coords( extend_star_array[i, index_of_ra],
extend_star_array[i, index_of_dec],
tolerance)
src_coord_array = np.asarray(src_coord_list)
if src_coord_list == None:
continue
# No nearby sources found in databases imply that's a new source.
if len(src_coord_list) == 0:
new_source_list.append([extend_star_array[i, index_of_name],
extend_star_array[i, index_of_ra],
extend_star_array[i, index_of_dec]])
continue
stu = find_sources(src_coord_array, tolerance)
failure, min_distance, jndex = stu.find([extend_star_array[i,index_of_ra],
extend_star_array[i,index_of_dec]])
# Update the name if we found a existed observation from database.
if not failure:
extend_star_array[i, index_of_name] = src_coord_array[jndex, index_of_name]
# Update the database if we don't.
else :
new_source_list.append([extend_star_array[i, index_of_name],
extend_star_array[i, index_of_ra],
extend_star_array[i, index_of_dec]])
# Check if we find new sources or not in this observations.
if len(new_source_list) > 0:
new = True
return extend_star_array, new_source_list, new
# This is a class for match the coordinates efficiently.
class find_sources():
def __init__(self, coord_table, tolerance = 0.0):
self.coord_table = coord_table
self.tolerance = tolerance
self.ref_coords = SkyCoord(self.coord_table[:,2], self.coord_table[:,3], unit = 'deg')
return
def find(self, coord):
source_coord = SkyCoord(coord[0], coord[1], unit = 'deg')
# Calculate the distance
distance_object_array = self.ref_coords.separation(source_coord)
distance_array = distance_object_array.deg
# Pick the nearest one
min_distance = np.min(distance_array)
index_min_distance = np.argmin(distance_array)
if min_distance < self.tolerance:
return False, min_distance, index_min_distance
else:
return True, 0.0, 0
# The def find the match name of target names.
def make_coord(src_list):
index_RA = TAT_env.src_titles.index("RA")
index_DEC = TAT_env.src_titles.index("`DEC`")
ra_array = src_list[index_RA]
dec_array= src_list[index_DEC]
ans_array = np.transpose([ra_array, dec_array])
return ans_array
#--------------------------------------------
# main code
if __name__ == "__main__":
# Measure time
start_time = time.time()
#----------------------------------------
# Initialize
if len(argv) != 2:
print "Error!\n The number of arguments is wrong."
print "Usage: starfinder.py [images list]"
exit(1)
name_image_list = argv[1]
image_list = np.loadtxt(name_image_list, dtype = str)
#----------------------------------------
# PSF register
for image_name in image_list:
# Find all stars with IRAFstarfinder
print ("{0}, start!".format(image_name))
print ('SExtractor!')
SE_table = SExtractor(image_name)
failure, extend_star_array = star_extend(image_name, SE_table)
# Rename if source is already named in previous observations.
print ('--- repeatness check ---')
extend_star_array, new_sources, new = check_new_sources(extend_star_array)
# Save and upload the result
save2sql(extend_star_array, new_sources, new)
print ("done.")
#---------------------------------------
# Measure time
elapsed_time = time.time() - start_time
print "Exiting Main Program, spending ", elapsed_time, "seconds."