-
Notifications
You must be signed in to change notification settings - Fork 3
/
donut
148 lines (116 loc) · 4.22 KB
/
donut
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
#! /usr/bin/env python
import os,sys
import numpy as np
import pylab as py
from donut.don11 import Donut
from astropy.io import fits
from astropy.table import Table
import json
import logging
from tempfile import NamedTemporaryFile
from mpi4py import MPI
comm = MPI.COMM_WORLD
root=0
rs_string ='[rank:i/i]'%(comm.rank,comm.size)
logging.basicConfig(format='%(asctime)s[%(levelname)s:%(threadName)s]-%(name)s-(%(filename)s:%(lineno)d):: %(message)s',
level=logging.DEBUG,
filename=os.path.basename(__file__) '.log',
filemode='w')
log = logging.getLogger(__name__)
def main(argv):
from optparse import OptionParser
parser = OptionParser()
parser.add_option('-i','--image',
help = 'Fit image file with donuts to measure.'
,type='string')
parser.add_option('-c','--catalog',
help = 'Sextractor catalog with sources position.'
,type='string')
parser.add_option('--start',
help = 'Initial index of the catalog to process (default = 0).'
,type='int',default=0)
parser.add_option('--end',
help = 'Last index of the catalog to process (default = -1).'
,type='int',default=-1)
parser.add_option('-p','--parameters',
help = 'Base donut json parameters. Position will be overwritten.'
,type='string')
parser.add_option('-o','--output',
help = 'Output file name.'
,type='string')
opt, args = parser.parse_args(argv)
# Read sextractor catalog
cat = Table.read(opt.catalog,
format='ascii.sextractor')[opt.start:opt.end]
# Read basic parameters
with open(opt.parameters) as par_file:
donutpars = json.load(par_file)
# Write donut parameter catalog
basepar = donutpars['donpar'][0]
parcatalog = {'donpar':[]}
for src in cat:
basepar['XC'] = int(src['X_IMAGE'])
basepar['YC'] = int(src['Y_IMAGE'])
parcatalog['donpar'].append(dict(basepar))
tmpCat = NamedTemporaryFile(delete=False)
log.info('Donut catalog file: %s'%tmpCat.name)
with open(tmpCat.name,'w') as fp:
json.dump(parcatalog,fp)
# Read image
img = fits.getdata(opt.image)
def donutfit(index):
don = Donut()
don.readpar(tmpCat.name,index)
don.init()
piximg = don.extract(img.T)
zres = np.zeros(basepar['NZER'] 3)
try:
if piximg.shape == (don.ngrid/2,don.ngrid/2):
x2,immod,zer = don.fit(piximg.T)
zres[0] = 1 # True
zres[1] = don.xc
zres[2] = don.yc
zres[3:] = zer
else:
log.warning('Source %i to close to the border. Skipping...'%(index))
zres[0] = 0 # False
except AttributeError,e:
log.exception(e)
zres[0] = 0 # False
pass
except Exception,e:
zres[0] = 0 # False
pass
return zres
joblist = []
if comm.rank == 0:
nprocs = comm.size
njobs_perproc = len(cat)/nprocs
joblist = np.arange(njobs_perproc*nprocs).reshape(nprocs,njobs_perproc)
proclist = comm.scatter(joblist,root)
# procres = np.zeros(len(proclist))
# Store the results
procres = np.zeros((len(proclist),basepar['NZER'] 3))
fitted = np.zeros(len(cat)) == 0
for i,p in enumerate(proclist):
procres[i] = donutfit(p)
recvbuf=comm.gather(procres,root)
if comm.rank==0:
log.info('Writing output to %s '%opt.output)
zres = np.array(recvbuf)
# mask = recvbuf[:,:1].T[0] == 1
# zres = recvbuf[:,1:]
np.save(opt.output,zres)
# for i in range(len(cat)):
# log.info('Working on %4i/%4i'%(i 1,len(cat)))
# print 'Working on %4i/%4i'%(i 1,len(cat))
# donutfit(i)
# # exclude non-fitted sources
sendbuf=""
if MPI.COMM_WORLD.rank==0:
sendbuf="done"
recvbuf=MPI.COMM_WORLD.bcast(sendbuf,root)
log.info(recvbuf)
# print tmpCat.name
if __name__ == '__main__':
main(sys.argv)