-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmakeChi2Cat_tract.py
423 lines (383 loc) · 18.9 KB
/
makeChi2Cat_tract.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
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 4 11:53:21 2016
@author: anneya
.-.
.--.( ).--.
<-. .-.-.(.-> )_ .--.
`-`( )-' `) )
(o o ) `)`-'
( ) ,)
( () ) )
`---"\ , , ,/`
`--' `--' `--'
| | | |
| | | |
' | ' |
For a single HSC tract, create a chi^2 - selected multiband catalog
"""
from __future__ import print_function, division
import sys, os
import glob
import argparse
import numpy as np
from astropy import units as u
from astropy.io import fits
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def GoodPatches(drPath,tract,bands):
'''
Returns a list of patches where a calexp images exists for all bands
'''
patches = []
for band in bands:
os.chdir(drPath+'deepCoadd/'+band+'/'+tract)
bandPatches = np.array(glob.glob('*'))
if len(patches)==0:
patches = bandPatches
else:
patches = np.intersect1d(np.array(patches),bandPatches)
return(patches)
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def KillCommas(drPath,tract,bands,patches):
'''
Remove commas from paths and image names so SExtractor can parse
'''
for band in bands:
for patch in patches:
os.chdir(drPath+'deepCoadd/'+band+'/'+tract+'/'+patch)
if os.path.isfile('calexp-'+band+'-'+tract+'-'+patch+'.fits') and patch[1]==',':
os.system('mv calexp-'+band+'-'+tract+'-'+patch+'.fits calexp-'+band+'-'+tract+'-'+patch[0]+patch[-1]+'.fits')
if not os.path.isfile('calexp-'+band+'-'+tract+'-'+patch[0]+patch[-1]+'.fits'):
raise ValueError('calexp-'+band+'-'+tract+'-'+patch[0]+patch[-1]+'.fits does not exist.')
for patch in patches:
if os.path.exists(drPath+'deepCoadd/'+band+'/'+tract+'/'+patch) and patch[1]==',':
os.system('mv '+drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/ '+\
drPath+'deepCoadd/'+band+'/'+tract+'/'+patch[0]+patch[2]+'/')
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def MakeChi2Image(fluxData,varMaps):
'''
Calculates chi2 image from arrays of flux data and corresponding variance maps
Inputs:
fluxData = NumPy array of shape(nImages,nXpixels,nYpixels) containing fluxes
varMaps = NumPy array of shape(nImages,nXpixels,nYpixels) containing variances
Outputs:
Numpy array of shape(nXpixels,nYpixels), chi2
'''
medians = np.median(fluxData.reshape(fluxData.shape[0],np.product(fluxData.shape[1:])),axis=1)
chi2 = np.zeros(shape=fluxData[0].shape)
for i in range(len(fluxData)):
chi2 += ((fluxData[i]-medians[i])/np.sqrt(varMaps[i]))**2
return(np.sqrt(chi2))
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def MakeChi2HSC(images,baseIm,outputImage,clobber):
'''
Calls MakeChi2Image() to create a chi2 image from a list of HSC images.
Inputs:
images = list of filenames of HSC images with extensions [IMAGE] and [VARIANCE]
baseIm = filename of an image whose [IMAGE] header will be copied for WCS information
outputImage = filename of output image
Outputs:
None
'''
fluxData,varMaps = [],[]
for image in images:
f = fits.open(image)
fluxData.append(f['IMAGE'].data)
varMaps.append(f['VARIANCE'].data)
header = fits.getheader(baseIm,'IMAGE')
for i in range(len(images)):
header['chi2in_'+str(i)] = images[i]
fluxData = np.array(fluxData)
varMaps = np.array(varMaps)
chi2 = MakeChi2Image(fluxData,varMaps)
fits.writeto(outputImage, chi2, header,clobber=clobber)
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def MakeChi2Images(drPath,tract,bands,patches,clobber,prefix):
'''
Make chi2 image and variance map for each patch in patches
'''
# Create chi2 'filter' directory with all patches
if not os.path.exists(drPath+'deepCoadd/'+prefix+'chi2/'):
os.system('mkdir '+drPath+'deepCoadd/'+prefix+'chi2/')
if not os.path.exists(drPath+'deepCoadd/'+prefix+'chi2/'+tract):
os.system('mkdir '+drPath+'deepCoadd/'+prefix+'chi2/'+tract)
for patch in patches:
if not os.path.exists(drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch):
os.system('mkdir '+drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch)
# Make images
for patch in patches:
outImage = drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch+'/'+prefix+'chi2-'+tract+'-'+patch+'.fits'
if not (clobber==False and os.path.isfile(outImage)):
images = []
for band in bands:
images.append(drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/calexp-'+band+'-'+tract+'-'+patch[0]+patch[-1]+'.fits')
MakeChi2HSC(images,images[0],outImage,clobber)
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def SExtractChi2(drPath,tract,patches,prefix,sexdir,dotsex):
'''
Run SExtractor on chi2 images so photometry from detection image can be used for corrections
'''
os.chdir(sexdir)
for patch in patches:
inImage = drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch+'/'+prefix+'chi2-'+tract+'-'+patch+'.fits'
outCat = drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch+'/'+prefix+'chi2-'+tract+'-'+patch+'.cat'
os.system('sex '+inImage+' -c '+dotsex+' -CATALOG_NAME '+outCat)
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def SExtractorDualImage(drPath,tract,bands,patches,prefix,sexdir,dotsex,zps):
'''
Run SExtractor in dual image mode for each patch,band using chi2 image for detection
'''
os.chdir(sexdir)
for band in bands:
for patch in patches:
# Move Variance extension to new file for SExtractor
imname = drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/'+\
'calexp-'+band+'-'+tract+'-'+patch[0]+patch[-1]
varHead = fits.getheader(imname+'.fits','VARIANCE')
maskHead = fits.getheader(imname+'.fits','MASK')
im = fits.open(imname+'.fits')
fits.writeto(imname+'_var.fits',im['VARIANCE'].data,varHead,clobber=True)
fits.writeto(imname+'_mask.fits',im['MASK'].data.astype(float),maskHead,clobber=True)
inImage = imname+'.fits[1]'
outCat = drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/'+band+'-'+tract+'-'+patch[0]+patch[-1]+'-chi2.cat'
chi2Image = drPath+'deepCoadd/'+prefix+'chi2/'+tract+'/'+patch+'/'+prefix+'chi2-'+tract+'-'+patch+'.fits'
# Run SExtractor
os.system('sex '+chi2Image+','+inImage+' -c '+dotsex+' -CATALOG_NAME '+outCat+' -WEIGHT_IMAGE None,'+\
imname+'_var.fits -WEIGHT_TYPE NONE,MAP_VAR -MAG_ZEROPOINT '+str(zps[bands.index(band)]))
# Add flags to catalog
os.system('./venice -m '+imname+'_mask.fits -cat '+outCat+ ' -f all -xcol 2 -ycol 3 -o '+\
drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/'+band+'-'+tract+'-'+patch[0]+patch[-1]+'-chi2-flags.cat')
os.system('rm '+imname+'_var.fits')
os.system('rm '+imname+'_mask.fits')
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def CatPatch(drPath,tract,bands,patch,onceCols,apertures,chi2prefix,zps):
'''
Generate multiband catalog for a given patch
'''
# Make Header for combined catalog
columns,bigColumns = [],[]
fluxCols,aperFluxCols = [],[]
f = open(drPath+'deepCoadd/'+bands[0]+'/'+tract+'/'+patch+'/'+bands[0]+'-'+tract+'-'+patch[0]+patch[-1]+'-chi2-flags.cat','r')
i = 0
for line in f.readlines():
if line[0]=='#':
if 'APER' in line:
for aperture in apertures:
columns.append(aperture+'_'+line[6:-1])
else:
columns.append(line[6:-1])
if 'FLUX_' in line and not 'RADIUS' in line and not 'APER' in line:
fluxCols.append(i)
if 'FLUX_APER' in line:
aperFluxCols.append(i)
i+=1
f.close()
columns.append('flags')
for i in onceCols:
bigColumns.append(columns[i])
for band in bands:
keepCols = []
for i in range(len(columns)):
if not i in onceCols:
keepCols.append(i)
bigColumns.append(band+'_'+columns[i])
header = ''
i = 0
for h in bigColumns:
header+='('+str(i)+') '+h+'\n'
i+=1
print(fluxCols)
# combine Catalogs from each band
bigCat = None
for band in bands:
cat = np.genfromtxt(drPath+'deepCoadd/'+band+'/'+tract+'/'+patch+'/'+band+'-'+tract+'-'+patch[0]+patch[-1]+'-chi2-flags.cat')
# convert fluxes to Jy
for col in fluxCols:
diffrac = cat[:,col]/cat[:,col+1]
newFlux = cat[:,col] * 10**((8.9-zps[bands.index(band)])/2.5)
newFluxErr = diffrac * newFlux
cat[:,col+1] = newFluxErr
cat[:,col] = newFlux
for col in aperFluxCols:
diffrac = cat[:,col]/cat[:,col+len(apertures)]
newFlux = cat[:,col] * 10**((8.9-zps[bands.index(band)])/2.5)
newFluxErr = diffrac * newFlux
cat[:,col] = newFlux
cat[:,col+len(apertures)] = newFluxErr
if bigCat == None:
bigCat = np.c_[cat]
else:
bigCat = np.c_[bigCat,cat[:,keepCols]]
# add unique identifiers
ids = []
for j in bigCat[:,0]:
ids.append(tract+patch+str(j))
ids = np.array(ids)
bigCat[:,0] = ids
print('LEN(HEAD)=',len(bigColumns),'nColumns=',len(bigCat[0]))
# identify objects in regions that overlap with other patches
# get patch dimension
imname = drPath+'deepCoadd/'+bands[0]+'/'+tract+'/'+patch+'/'+\
'calexp-'+bands[0]+'-'+tract+'-'+patch[0]+patch[-1]
f = fits.open(imname+'.fits')
# Add binary flags for objects in overlap region and those to be removed in tract catalog
patchDims=[f[1].data.shape[0],f[1].data.shape[1]]
xmax,ymax=patchDims[1]-100,patchDims[0]-100
overlap = np.logical_not(np.sum(np.c_[bigCat[:,1]<100,bigCat[:,2]<100,bigCat[:,1]>xmax,bigCat[:,2]>ymax].astype(int),axis=1))
xmax,ymax=patchDims[1]-50,patchDims[0]-50
remove = np.logical_not(np.sum(np.c_[bigCat[:,1]<50,bigCat[:,2]<50,bigCat[:,1]>xmax,bigCat[:,2]>ymax].astype(int),axis=1))
header+='('+str(i)+') in_overlap_region \n'
i+=1
header+='('+str(i)+') remove_for_unique'
np.savetxt(drPath+'chicat/'+tract+'_'+patch+'_'+chi2prefix+'_chi2.cat',np.c_[bigCat,overlap,remove],fmt='%f4',header=header.replace('count','jansky'))
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def TractCat(drPath,tract,patches,chi2prefix):
'''
Combine multiband catalogs for each patch in tract
'''
flagCols,header = [],''
f = open(drPath+'chicat/'+tract+'_'+patches[0]+'_'+chi2prefix+'_chi2.cat','r')
i = 0
for line in f.readlines():
if line[0]=='#':
header+=line[2:]
if 'flags' in line:
flagCols.append(i)
i+=1
f.close()
bigCat = None
for patch in [patches[0]]:
cat = np.genfromtxt(drPath+'chicat/'+tract+'_'+patch+'_'+chi2prefix+'_chi2.cat')
flags = cat[:,flagCols].astype(int)
bad = flags>=256
bad = np.sum(np.c_[bad,np.logical_not(cat[:,-1].astype(int).astype(bool)).astype(int)],axis=1)
goodN = np.where(bad==0)
if bigCat == None:
bigCat = cat[goodN]
else:
bigCat = np.r_[bigCat,cat[goodN]]
np.savetxt(drPath+'chicat'+'/'+tract+'_'+chi2prefix+'_chi2.cat',bigCat,fmt='%d%06.d '+(len(bigCat[0])-2)*'%f ',header=header[:-1])
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
def CatToFits(inCat,outCat,deredCat,cs,bands,getDustLoc):
'''
Create .fits table from ascii catalog, deredden
'''
colNames,colUnits,colFormats = [],[],[]
f = open(inCat,'r')
for line in f.readlines():
if line[0]=='#':
lis = line[3:-1].split()
name = lis[1]
colNames.append(name.replace('-','_'))
if name=='NUMBER' or 'flag' in name:
colFormats.append('K')
else:
colFormats.append('D')
if '[' in line and ']' in line:
unit = line[line.index('[')+1:line.index(']')]
if unit =='deg':
unit = 'degree'
colUnits.append(unit)
else:
colUnits.append(None)
else:
break
f.close()
cat = np.genfromtxt(inCat)
cols = []
for k in range(len(cat[0])):
cols.append(fits.Column(name=colNames[k],format=colFormats[k],unit=colUnits[k],array=cat[:,k]))
cols = fits.ColDefs(cols)
tbhdu = fits.BinTableHDU.from_columns(cols)
tbhdu.writeto(outCat,clobber=True)
# deredden
print('Dereddening...')
coeffs,bs = '',''
for c in colNames:
if 'FLUX' in c and not 'ERR' in c and not 'RADIUS' in c and not 'CHI2' in c:
bs+=c+','
ind = None
for band in bands:
band=band.replace('-','_')
if band in c:
ind = c.index(band)
break
coeffs+=str(cs[ind])+','
os.system(getDustLoc+'./getDust.py '+outCat+' '+deredCat+' -band '+bs[:-1]+' -coef '+coeffs[:-1]+' -correct '+\
'-s '+getDustLoc+'SFD_dust_4096_sgp.fits -n '+getDustLoc+'SFD_dust_4096_ngp.fits' )
#------------------------------------------------------------------------------
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#------------------------------------------------------------------------------
# Set up parser
parser = argparse.ArgumentParser()
parser.add_argument("-v", "--verbosity", action="count", default=0, help="increase output verbosity")
parser.add_argument("-b", "--bands",nargs='+', help="bands to include in chi2 images",\
default=['MegaCam-u','HSC-G','HSC-R','HSC-I','HSC-Z','HSC-Y'])
parser.add_argument("-p","--dataPath",help="full path to data release (where deepCoadds is a subfolder)",\
default = '~/data/HSC/s15b_deep/')
parser.add_argument("--tract",help='HSC tract id', default = '8524')
parser.add_argument("--clobberChi2",help='if True, will overwrite existing chi2 images',\
default = True)
parser.add_argument("--chi2prefix",help="prefix to be added to chi2 images",default='ugrizy')
parser.add_argument("--sexdir",help="Path to directory where sextractor control files live, venice and getDust.py stuff assumed there",\
default = '~/local/source/HSCchi2Cat/')
parser.add_argument("--dotsex",help="SExtractor .param file",default='default.sex')
parser.add_argument("--zps",help="ZeroPoints for each band in bands",nargs='+',\
default=[30.,27.,27.,27.,27.,27.])
parser.add_argument("--onceCols",help="Columns in dual image SExtractor catalogs that should only appear once in combined catalog",\
nargs='+',default=[0,1,2,3,4])
parser.add_argument("--apertures",help="List of apertures",nargs='+',\
default=['12pix','24pix'])
parser.add_argument("--deredCoeffs",help="Corresponding Albda/E(B-V) (in mags)",nargs='+',\
default=[4.732,3.711,2.626,1.916,1.469,1.242])
args = parser.parse_args()
# Locate patches where all bands have coadds
print('Locating patches where all bands have coadds...')
patches = GoodPatches(args.dataPath,args.tract,args.bands)
# Remove commas
print('Removing commas...')
KillCommas(args.dataPath,args.tract,args.bands,patches)
# Get new patch names
patches = GoodPatches(args.dataPath,args.tract,args.bands)
# Make chi2 image for each good patch
print('Making chi^2 images...')
MakeChi2Images(args.dataPath,args.tract,args.bands,patches,args.clobberChi2,args.chi2prefix)
# Run SExtractor on chi2 images
print('SExtracting on chi2 images...')
SExtractChi2(args.dataPath,args.tract,patches,args.chi2prefix,args.sexdir,args.dotsex)
# Run SExtractor in dual image mode on other bands
print('SExtracting in dual image mode...')
SExtractorDualImage(args.dataPath,args.tract,args.bands,patches,args.chi2prefix,args.sexdir,args.dotsex,args.zps)
# Create a multiband catalog for each patch
print('Creating multiband catalogs for each patch...')
for patch in patches:
CatPatch(args.dataPath,args.tract,args.bands,patch,args.onceCols,args.apertures,args.chi2prefix,args.zps)
# Make multiband catalog for tract from catalogs from each patch
print('Creating catalog for tract '+args.tract+'...')
TractCat(args.dataPath,args.tract,patches,args.chi2prefix)
# Convert to fits and deredden
print('Creating .fits catalog...')
CatToFits(args.dataPath+'chicat'+'/'+args.tract+'_'+args.chi2prefix+'_chi2.cat',\
args.dataPath+'chicat'+'/'+args.tract+'_'+args.chi2prefix+'_chi2.fits',\
args.dataPath+'chicat'+'/'+args.tract+'_'+args.chi2prefix+'_chi2_dered.fits',\
args.deredCoeffs,args.bands,args.sexdir)