-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathinput_skymaker.py
191 lines (156 loc) · 7.38 KB
/
input_skymaker.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
from numpy import *
import numpy as np
import os
import csv
from astropy import wcs
import astropy.units as u
from astropy import coordinates as coord
from astroquery.sdss import SDSS
from astroquery.vizier import Vizier
#import matplotlib.pyplot as plt
def makelist(ra, dec, ccd, date, lname, variable=0., transients=0.):
vname = date+'/var/'+lname+'.var'
tname = date+'/var/'+lname+'.trans'
rname = date+'/reg/'+lname+'.reg'
radec = date+'/radec/'+lname+'.txt'
#Define size of camera in pixels:
xsi = 8176.
ysi = 6132.
#Corresponding size in degs:
wid = xsi*1.24/3600.
hei = ysi*1.24/3600.
#Define the columns to get from UCAC, magnitude filter and set unlimited rows:
v = Vizier(columns=['_RAJ2000', '_DEJ2000', 'Vmag'], column_filters={"Vmag":"<17."})
v.ROW_LIMIT=-1
#Query UCAC4
ucac = v.query_region(\
coord.SkyCoord(ra=ra, dec=dec,\
unit=(u.deg, u.deg),\
frame='icrs'),\
height=str(hei+1)+"d", width=str(wid+1)+"d",\
catalog=["UCAC4"])[0]
#For SDSS, need to calculate bounding box:
declim = dec+np.array([-hei,hei])/2.
declim_r = np.pi*declim/180.
wid_r = np.pi*wid/180.
d_lim_r = 2.*np.arcsin(np.sin(wid_r/4.)/np.cos(declim_r))
ralim = ra+np.array([-1,1])*np.max(180.*d_lim_r/np.pi)
#Query SDSS:
query = "SELECT p.objid, p.ra, p.dec, p.g, "+\
"p.deVRad_g, p.deVPhi_g, p.deVAB_g, "+\
"p.type, p.flags_g, (flags & dbo.fPhotoFlags('SATURATED')) as SAT "+\
"FROM PhotoPrimary AS p "+\
"WHERE p.ra BETWEEN "+str(ralim[0])+" AND "+str(ralim[1])+" AND "+\
"p.dec BETWEEN "+str(declim[0])+" AND "+str(declim[1])+" AND "+\
"p.g BETWEEN 16 AND 22"# AND "+\
#"p.htmid*37 & 0x000000000000FFFF < (650 * 10)"
sdss = SDSS.query_sql(query)
#Remove saturated sources and separate gals from stars:
o = np.where(sdss['SAT'] == 0)
sdss = sdss[o]
#Find and remove matching sources:
c = coord.SkyCoord(ra=sdss['ra']*u.degree, dec=sdss['dec']*u.degree)
cs = coord.SkyCoord(ra=np.array(ucac['_RAJ2000'])*u.degree, \
dec=np.array(ucac['_DEJ2000'])*u.degree)
idx, d2d, d3d = cs.match_to_catalog_sky(c)
match = np.where(d2d.value*3600. > 1.0)
ucac = ucac[match]
#Group stars from SDSS and UCAC together:
o = np.where(sdss['type'] == 6)
sdss_stars = sdss[o]
stars = np.append(np.array([ucac['_RAJ2000'], ucac['_DEJ2000'], ucac['Vmag']]),\
np.array([sdss_stars['ra'], sdss_stars['dec'], sdss_stars['g']]),\
axis=1)
#Extract gals from SDSS:
o = np.where(sdss['type'] != 6)
sdss_gals = sdss[o]
gals = np.array([sdss_gals['ra'], sdss_gals['dec'],\
sdss_gals['g'],\
sdss_gals['deVRad_g'], sdss_gals['deVAB_g'], sdss_gals['deVPhi_g']])
#Generate WCS to convert (RA,Dec) to (x,y)
w = wcs.WCS(naxis=2)
w.wcs.crpix = np.array([4088, 3066])+np.random.normal(scale=3.,size=2)
w.wcs.cdelt = np.array([3.444e-4, 3.444e-4])
w.wcs.crval = [ra, dec] #Pointing position of telescope.
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
pixcrds = np.column_stack((stars[0,:],stars[1,:]))
pixcrdg = np.column_stack((gals[0,:], gals[1,:]))
worldg = w.wcs_world2pix(pixcrdg, 1)
worlds = w.wcs_world2pix(pixcrds, 1)
gals[0,:] = worldg[:,0]
gals[1,:] = worldg[:,1]
stars[0,:] = worlds[:,0]
stars[1,:] = worlds[:,1]
ind = [-1]
if variable > 0:
#Select a set proportion of m<19 stars and add random variation:
bright = np.squeeze(np.where((stars[2,:]<19) &\
(stars[0,:]>0) & (stars[0,:]<xsi) &\
(stars[1,:]>0) & (stars[1,:]<ysi)))
n_vary = int(np.round(variable*(np.size(bright))))
if n_vary > 0:
ind = np.random.choice(bright, n_vary, replace=False)
orig = stars[2,ind]
stars[2,ind] = stars[2,ind] + np.random.normal(0, 1, ind.size)
np.savetxt(vname, np.c_[pixcrds[ind,0], pixcrds[ind,1], orig, stars[2,ind]], \
header = 'RA Dec Orig_Mag New_Mag', \
fmt='%9.5f %8.5f %5.2f %5.2f')
if transients > 0:
#Randomly distribute transients around image:
tx = np.random.uniform(0, xsi, transients)
ty = np.random.uniform(0, ysi, transients)
tm = np.random.uniform(14, 19, transients)
t = np.array([tx, ty, tm])
stars = np.append(stars, t, axis=1)
txy = (np.array([tx,ty])).transpose()
pixcrdt = w.wcs_pix2world(txy, 1)
#Append transients to stars:
pixcrds = np.append(pixcrds,pixcrdt, axis=0)
np.savetxt(tname, np.c_[pixcrdt[:,0], pixcrdt[:,1], tm], \
header = 'RA Dec Orig_Mag New_Mag', \
fmt='%9.5f %8.5f %5.2f')
#Write stars to file:
myfile = open('templist'+ccd+'.list','w')
regfile = open(rname,'w')
regfile.write('image\n')
obj = np.concatenate((100*np.ones(pixcrds[:,0].size), 200*np.ones(pixcrdg[:,0].size)))
ras = np.concatenate((pixcrds[:,0], pixcrdg[:,0]))
decs = np.concatenate((pixcrds[:,1], pixcrdg[:,1]))
mags = np.concatenate((stars[2,:], gals[2,:]))
np.savetxt(radec, \
np.c_[obj, ras, decs, mags], \
fmt='%3i %9.5f %8.5f %5.2f', \
header='Type RA DEC Mag')
for i in range(0,(stars.shape)[1]):
myfile.write((str(100)+' '+str(stars[0,i])+' '+str(stars[1,i])+' '+str(stars[2,i]))+'\n')
if np.any(ind == i):
regfile.write('circle('+str(stars[0,i])+','+str(stars[1,i])+',3) #color=blue\n')
else:
regfile.write('circle('+str(stars[0,i])+','+str(stars[1,i])+',3)\n')
#Write gals to file:
for i in range(0,(gals.shape)[1]):
myfile.write((str(200)+' '+str(gals[0,i])+' '+str(gals[1,i])+' '+str(gals[2,i])+' ' +\
str(0)+' ' +str(0)+' ' +str(0)+' ' +str(0)+' ' + str(gals[3,i])+' '+\
str(gals[4,i])+' '+str(gals[5,i])+'\n'))
regfile.write('circle('+str(gals[0,i])+','+str(gals[1,i])+',3) #color=red\n')
myfile.close()
regfile.close()
#What needs to be done now is:
# - Generate a list of pointings for the mount, covering ~100sq deg.
# Note that SDSS only covers certain parts of the sky, so you need
# to point only in these regions.
# - For each pointing, calculate the pointing for each telescope.
# There will need to be some overlap between pointings.
# - Generate a filename `fname` for each telescope; see below for example.
# The source list will have name `fname`.list.
# - Automatically edit the .conf file to include the correct output
# filename `fname`.fits and its own PSF. Each simulation should
# have its own .conf file called `fname`.conf
# - We also need to include telescope jitter.
# - Edit the header information to reflect the settings for each pointing.
# - Loop over the mount pointings, generating twelve 2-minute sim images
# (three for each telescope) per mount pointing. We can
# use the same PSF for each set of twelve. In fact, we can probably assume
# the PSf remains stable for about 30 minutes (so 5 mount pointings).
#fname="GOTO_01_20170323_00001"
#makelist(180.,30.,fname)