Here is an example of a mosaic script, using pyramid procedures.
1 #!/usr/bin/env python
2 #
3 # History:
4 # june 02 mchw. ALMA script.
5 # 15aug02 mchw. CARMA version edited from ALMA script.
6 # 23aug02 mchw. calculate region from source size.
7 # 20sep02 mchw. Re-import CARMA improvements for ALMA.
8 # 25sep02 mchw. Re-import improvements from hex7.csh to hex19.csh
9 # 26sep02 mchw. Increase imsize from 129 to 257.
10 # 12mar03 mchw. convert to PYTHON.
11 # 13mar03 pjt more conversion to PYTHON, now at 200ft, renamed to mosaic.py
12
13 import sys, os, time, string, math
14 from Miriad import *
15
16 version='2003-03-14'
17
18 print " --- ALMA Mosaicing (Cas A model) --- "
19
20 # command line arguments that can be changed...
21 keyval = {
22 "config" : "config1", # antenna config file (without the .ant extension)
23 "dec" : "-30", # declination (can be a real number)
24 "image" : "casc.vla", # image to test (nice Cas-A VLA image as default)
25 "cell" : "0.04", # scale size (should be calculated from image)
26 "nchan" : "1", # number of channels
27 "method" : "mosmem", # mosmem, joint, or default
28 "flux" : "732.063", # expected flux in the image (for mosmem)
29 "nring" : "3", # number of rings in the mosaic
30 "grid" : "12.0", # gridsize (in arcsec) for the mosaic
31 "center" : "", # optional center file that overrides (nring,grid)
32 "VERSION" : "1.0 mchw" # VERSION id for the user interface
33 }
34
35 help = """
36 The minimum amount of information you need to run this task is:
37 a miriad image (image=) for the model.
38 an antenna configuration file (<config>.ant) for uvgen
39 """
40
41 keyini(keyval,help,0)
42 # report current defaults, exit if --help given
43 setlogger('mosaic.log')
44 #
45 # -----------------------------------------------------------------------------
46
47 #
48 # define all variables, now in their proper type, for this script
49 #
50
51 config = keya('config')
52 dec = keyr('dec')
53 cell = keyr('cell')
54 nchan = keyi('nchan')
55 method = keya('method')
56 center = keya('center')
57 flux = keyr('flux')
58 image = keya('image')
59 nring = keyi('nring')
60 grid = keyr('grid')
61
62 harange = '-1,1,0.013'
63 select = '-shadow\(12\)'
64 freq = 230.0
65 imsize = 257 # avoid 2**N, image size 2**N + 1 is good. [or calculate from image]
66
67 mir = os.environ['MIR']
68
69 # -----------------------------------------------------------------------------
70
71 # returns a list of strings that are the ascii centers as uvgen wants them
72 # (in a file) via the center= keyword
73 def hex(nring,grid):
74 center=""
75 npoint=0
76 for row in range(-nring+1,nring,1):
77 y = 0.866025403 * grid * row
78 lo = 2-2*nring+abs(row)
79 hi = 2*nring-abs(row)-1
80 for k in range(lo,hi,2):
81 x = 0.5*grid*k
82 npoint = npoint + 1
83 if center=="":
84 center = center + "%.2f,%.2f" % (x,y)
85 else:
86 center = center + ",%.2f,%.2f" % (x,y)
87 return (npoint,center)
88
89 # get the (as a string) value of an item in a dataset
90 def itemize(data,item):
91 log = 'tmp.log'
92 cmd = [
93 'itemize',
94 'in=%s/%s' % (data, item),
95 'log=%s' % log
96 ]
97 miriad(cmd)
98 f = open(log,"r")
99 v = string.split(f.readline())
100 f.close()
101 return v[2]
102
103 # copy a file from source (s) to destination (d)
104 def copy_data(s,d):
105 zap(d)
106 os.system("cp -r %s %s" % (s,d))
107
108
109 # should this be
110 # units=None
111 # if units is None:
112 # bla
113 # else:
114 # bla
115
116 def puthd(map,item,value,units=0):
117 cmd = [
118 'puthd',
119 'in=%s/%s' % (map,item),
120 ]
121 if (units == 0):
122 cmd.append("value=%s" % value)
123 else:
124 cmd.append("value=%s,%s" % (value,units))
125 return cmd
126
127 def demos(map,vis,out):
128 cmd = [
129 'demos',
130 'map=%s' % map,
131 'vis=%s' % vis,
132 'out=%s' % out,
133 ]
134 zap_all(out+"*")
135 return cmd;
136
137 def invert(vis,map,beam,imsize,select):
138 cmd = [
139 'invert',
140 'vis=%s' % vis,
141 'map=%s' % map,
142 'beam=%s' % beam,
143 'imsize=%d' % imsize,
144 'select=%s' % select,
145 'sup=0',
146 'options=mosaic,double',
147 ]
148 zap(map)
149 zap(beam)
150 return cmd
151
152 def cgdisp(map):
153 cmd = [
154 'cgdisp',
155 'in=%s' % map,
156 'device=/xs',
157 'labtyp=arcsec',
158 'range=0,0,lin,8'
159 ]
160 return cmd;
161
162 def uvmodel(vis,model,out):
163 cmd = [
164 'uvmodel',
165 'vis=%s' % vis,
166 'model=%s' % model,
167 'out=%s' % out,
168 'options=add,selradec'
169 ]
170 zap(out)
171 return cmd;
172
173 def implot(map,region='quarter'):
174 cmd = [
175 'implot',
176 'in=' + map,
177 'device=/xs',
178 'units=s',
179 'conflag=l',
180 'conargs=1.4',
181 'region=%s' % region
182 ]
183 return cmd
184
185 def imlist(map):
186 cmd = [
187 'imlist',
188 'in=' + map,
189 'options=mosaic'
190 ]
191 return cmd
192
193
194 def imgen(map,out,pbfwhm):
195 cmd = [
196 'imgen',
197 'in=%s' % map,
198 'out=%s' % out,
199 'object=gaussian',
200 'factor=0',
201 'spar=1,0,0,%g,%g' % (pbfwhm,pbfwhm)
202 ]
203 zap(out)
204 return cmd
205
206
207 def mosmem(map,beam,out,region,flux=0,default=0):
208 cmd = [
209 'mosmem',
210 'map=%s' % map,
211 'beam=%s' % beam,
212 'out=%s' % out,
213 'region=%s' % region,
214 'rmsfac=200,1',
215 # 'niters=200'
216 'niters=2'
217 ]
218 if flux != 0:
219 cmd.append('flux=%g' % flux)
220 if default != 0:
221 cmd.appenx('default=%s' % default)
222 zap(out)
223 return cmd
224
225 def restor(map,beam,model,out):
226 cmd = [
227 'restor',
228 'model=%s' % model,
229 'map=%s' % map,
230 'beam=%s' % beam,
231 'out=%s' % out
232 ]
233 zap(out)
234 return cmd
235
236 def regrid(map,tin,out):
237 cmd = [
238 'regrid',
239 'in=%s' % map,
240 'tin=%s' % tin,
241 'out=%s' % out,
242 'axes=1,2'
243 ]
244 zap(out)
245 return cmd
246
247 def convol(map,out,b1,b2,pa):
248 cmd = [
249 'convol',
250 'map=%s' % map,
251 'out=%s' % out,
252 'fwhm=%g,%g' % (b1,b2),
253 'pa=%g' % pa
254 ]
255 zap(out)
256 return cmd
257
258 def imframe(map,out):
259 cmd = [
260 'imframe',
261 'in=%s' % map,
262 'out=%s' % out,
263 'frame=-1024,1024,-1024,1024' # TODO:: the 1024 here depends on the input image size
264 ]
265 zap(out)
266 return cmd
267
268 def uvgen(ant,dec,harange,freq,nchan,out,center):
269 cmd = [
270 'uvgen',
271 'ant=%s' % ant,
272 'baseunit=-3.33564',
273 'radec=23:23:25.803,%g' % dec,
274 'lat=-23.02',
275 'harange=%s' % harange,
276 'source=$MIRCAT/point.source',
277 'telescop=alma',
278 'systemp=40',
279 # 'pnoise=30',
280 'jyperk=40',
281 'freq=%g' % freq,
282 'corr=%d,1,0,8000' % nchan,
283 'out=%s' % out,
284 'center=%s' % center # notice we don't use a file, but a string of numbers
285 ]
286 zap(out)
287 return cmd
288
289 def imdiff(in1,in2,resid):
290 cmd = [
291 'imdiff',
292 'in1=%s' % in1,
293 'in2=%s' % in2,
294 'resid=%s' % resid,
295 'options=nox,noy,noex'
296 ]
297 zap(resid)
298 return cmd
299
300 def histo(map,region):
301 cmd = [
302 'histo',
303 'in=%s' % map,
304 'region=%s' % region
305 ]
306 return cmd
307
308 # ================================================================================
309 #
310 # start of the actual script
311 # -----------------------------------------------------------------------------
312 # Nyquist sample rate for each pointing.
313 # calc '6/(pi*250)*12'
314 cells = 500*cell
315 region = "arcsec,box\(%.2f,-%.2f,-%.2f,%.2f\)" % (cells,cells,cells,cells)
316
317 ant = config + '.ant' # antenna file for uvgen
318 uv = config + '.uv' # dataset for visibilities
319
320 demos1 = "%s.cas.%g.demos" % (config,cell)
321 base1 = "%s.%g.cas.%g" % (config,dec,cell)
322 base2 = "single.%g.cas.%g" % (dec,cell)
323
324 map1 = base1 + ".mp"
325 beam1 = base1 + ".bm"
326 map2 = base2 + ".map"
327 beam2 = base2 + ".beam"
328 mem = base1 + ".mem"
329 cm = base1 + ".cm"
330 mp = base1 + '.mp'
331 res = base1 + '.resid'
332 conv = base1 + '.conv'
333
334 if center == "":
335 (npoint,center) = hex(nring,grid)
336 print "MOSAIC FIELD, using hexagonal field with nring=%d and grid=%g (%d pointings) " % (nring,grid,npoint)
337 else:
338 centerfile = center
339 f = open(centerfile,"r")
340 center=f.read()
341 f.close()
342 npoint = len(string.split(center,","))-1
343 center=string.replace(center,'\n',',')
344 print "MOSAIC FIELD, using center file %s (%d pointings) " % (centerfile,npoint)
345
346 print " --- ALMA Mosaicing (Cas A model) --- "
347
348 print " config = %s" % config
349 print " dec = %g" % dec
350 print " scale = %g" % cell
351 print " harange = %s hours" % harange
352 print " select = %s" % select
353 print " freq = %g" % freq
354 print " nchan = %d" % nchan
355 print " imsize = %d" % imsize
356 print " region = %s" % region
357 print " method = %s" % method
358 print " "
359 print " --- TIMING --- "
360
361 if method == "mosmem":
362 print "Generate mosaic grid"
363 # lambda/2*antdiam (arcsec)
364 print 300/freq/2/12e3*2e5
365
366 print "Generate uv-data. Tsys=40K, bandwidth=8 GHz "
367 miriad(uvgen(ant,dec,harange,freq,nchan,uv,center))
368 os.system('uvindex vis=%s' % uv)
369
370 print "Scale model size from pixel 0.4 to %g arcsec" % cell
371 # with 0.4 arcsec pixel size Cas A is about 320 arcsec diameter; image size 1024 == 409.6 arcsec
372 # scale model size. eg. cell=0.1 arcsec -> 80 arcsec cell=.01 -> 8 arcsec diameter
373
374 copy_data(image,base2)
375
376 miriad(puthd(base2,'crval2',dec,units='dms'))
377 miriad(puthd(base2,'crval3',freq))
378 miriad(puthd(base2,'cdelt1',-cell,units='arcsec'))
379 miriad(puthd(base2,'cdelt2',cell,units='arcsec'))
380
381 print "Make model images for each pointing center"
382 miriad(demos(base2,uv,demos1))
383
384 print "Make model uv-data using VLA image of Cas A as a model (the model has the VLA primary beam)"
385 for i in range(1,npoint+1):
386 miriad(cgdisp(demos1+"%d"%i))
387 vis_i = base1+".uv%d"%i
388 demos_i = demos1+"%d"%i
389 miriad(uvmodel(uv,demos_i,vis_i))
390 if i==1:
391 vis_all=vis_i
392 else:
393 vis_all=vis_all + ',' + vis_i
394 print "UVMODEL: add the model to the noisy sampled uv-data"
395
396 miriad(invert(vis_all, base1+".mp", base1+".bm", imsize, select))
397
398 print "INVERT: "
399
400 miriad(implot(base1+'.mp',region=region))
401 miriad(imlist(base1+'.mp'))
402
403 print "Make single dish image and beam"
404
405 pbfwhm = string.atof(grepcmd("pbplot telescop=alma freq=%g" % freq, "FWHM", 2)) * 60.0
406
407 print "Single dish FWHM = %g arcsec at %g GHz" % (pbfwhm,freq)
408
409 miriad(imframe(base2,base2+".bigger"))
410 miriad(convol(base2+".bigger",base2+".bigger.map",pbfwhm,pbfwhm,0.0))
411 miriad(regrid(base2+".bigger.map",base1+'.mp',base2+".map"))
412 miriad(imgen(base2+".map",base2+".beam",pbfwhm))
413 miriad(implot(base2+".map"))
414 miriad(puthd(base2+".map",'rms','7.32')) # is that 1/100 of the flux ???
415
416 if method=='mosmem':
417 print " MOSMEM Interferometer only"
418 print " MOSMEM Interferometer only with niters=200 flux=%g rmsfac=200." % flux
419 miriad(mosmem(map1,beam1,mem,region,flux=flux))
420 elif method=='joint':
421 print "Joint deconvolution of interferometer and single dish data"
422 print "Joint deconvolution of interferometer and single dish data ; niters=200 rmsfac=200,1"
423 miriad(mosmem(map1+','+map2,beam1+','+beam2,mem,region))
424 elif method=='default':
425 print "MOSMEM with default single dish image"
426 print "MOSMEM with default single dish image; niters=200 rmsfac=200"
427 miriad(mosmem(map1,beam1,mem,region))
428 else:
429 print "Unknown method " + method
430
431 miriad(restor(map1,beam1,mem,cm))
432 miriad(implot(map1,region=region))
433
434 print "convolve the model by the beam and subtract from the deconvolved image"
435 b1 = string.atof(grepcmd('prthd in=%s' % cm, 'Beam', 2))
436 b2 = string.atof(grepcmd('prthd in=%s' % cm, 'Beam', 4))
437 b3 = string.atof(grepcmd('prthd in=%s' % cm, 'Position', 2))
438
439 miriad(convol(base2,base1+'.conv',b1,b2,b3))
440 miriad(implot(base1+'.conv',region=region))
441
442 print "regrid the convolved model to the deconvolved image template"
443
444 miriad(regrid(base1+".conv",base1+".cm",base1+'.regrid'))
445 miriad(implot(base1+'.regrid',region=region))
446
447 # skipping cgdisp /gif production
448
449 miriad(imdiff(base1+'.cm',base1+'.regrid',base1+'.resid'))
450 miriad(implot(base1+'.resid',region=region))
451 miriad(histo(base1+'.resid',region=region))
452
453 # ================================================================================
454
455 print "print out results - summarize rms and beam sidelobe levels"
456 print " --- RESULTS --- "
457
458 # extract information, the hard way
459
460 # BUG: doesn't look like 'mp' has rms???
461 #rms = string.atof(itemize(mp,'rms')) * 1000
462 rms = -1
463 srms = string.atof(grepcmd('histo in=%s' % res, 'Rms', 3))
464 smax = string.atof(grepcmd('histo in=%s' % res, 'Maximum', 2))
465 smin = string.atof(grepcmd('histo in=%s' % res, 'Minimum', 2))
466 Model_Flux = string.atof(grepcmd('histo in=%s region=%s' % (conv,region),'Flux',5))
467 Model_Peak = string.atof(grepcmd('histo in=%s region=%s' % (conv,region),'Maximum',2))
468 Flux = string.atof(grepcmd('histo in=%s region=%s' % (cm,region),'Flux',5))
469 Peak = string.atof(grepcmd('histo in=%s region=%s' % (cm,region),'Maximum',2))
470 Fidelity = Peak/srms
471
472 print " Config DEC HA[hrs] Beam[arcsec] scale Model_Flux,Peak Image_Flux,Peak Residual:Rms,Max,Min[Jy] Fidelity"
473 print " %s %g %s %.3f %g %g %g %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f" % (config,dec,harange,rms,b1,b2,
474 cell,Model_Flux,Model_Peak,Flux,Peak,srms,smax,smin,Fidelity)
475
476 #mv timing hex19.$config.$dec.$harange.$nchan.$imsize
477 #cat $config.$dec.$harange.$nchan.$imsize
478 #cat casa.results
479 #enscript -r casa.results
480
481 #print "DEBUGGING"
482 #string.atof(itemize(mp,'rms'))
C.-0