Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/cleanhelper.py: 6%
2054 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
2import glob
3import os
4import math
5import numpy
6import re
7import shutil
8import pwd
9from numpy import unique
11import subprocess
12from collections import OrderedDict as odict
14###some helper tools
15from casatasks import casalog as default_casalog
16from casatools import table, quanta, measures, regionmanager, image, imager, msmetadata
17from casatools import ms as mstool
18from casatools import ctsys
20ms = mstool( )
21tb = table( )
22qa = quanta( )
23me = measures( )
24rg = regionmanager( )
25ia = image( )
26im = imager( )
27msmd=msmetadata( )
29# trying to avoid sharing a single instance with all other functions in this module
30iatool = image
32def _casa_version_string():
33 """ produce a version string with the same format as the mstools.write_history. """
34 return 'version: ' + ctsys.version_string() + ' ' + ctsys.version_desc()
37class cleanhelper:
38 def __init__(self, imtool='', vis='', usescratch=False, casalog=default_casalog):
39 """
40 Contruct the cleanhelper object with an imager tool
41 like so:
42 a=cleanhelper(im, vis)
43 """
44 ###fix for assumption that if it is a list of 1 it is sngle ms mode
45 if((type(vis)==list) and (len(vis)==1)):
46 vis=vis[0]
47 ####
48 if((type(imtool) != str) and (len(vis) !=0)):
49 # for multi-mses input (not fully implemented yet)
50 if(type(vis)!=list):
51 vis=[vis]
52 self.sortedvisindx=[0]
53 self.initmultims(imtool, vis, usescratch)
54 # else:
55 # self.initsinglems(imtool, vis, usescratch)
56 #self.maskimages={}
57 self.maskimages=odict( )
58 self.finalimages={}
59 self.usescratch=usescratch
60 self.dataspecframe='LSRK'
61 self.usespecframe=''
62 self.inframe=False
63 # to use phasecenter parameter in initChannelizaiton stage
64 # this is a temporary fix need.
65 self.srcdir=''
66 # for multims handling
67 self.sortedvislist=[]
68 if not casalog: # Not good!
69 casalog = default_casalog
70 #casalog.setglobal(True)
71 self._casalog = casalog
73 @staticmethod
74 def getsubtable(visname='', subtab='SPECTRAL_WINDOW'):
75 """needed because mms has it somewhere else
76 """
77 tb.open(visname)
78 spectable=str.split(tb.getkeyword(subtab))
80 if(len(spectable) ==2):
81 spectable=spectable[1]
82 else:
83 spectable=visname+"/"+subtab
84 return spectable
86 @staticmethod
87 def checkimageusage(imagename=''):
88 """
89 check if images which will be written into is open elsewhere
90 """
91 retval=[]
92 if(imagename=='' or imagename==[]):
93 return retval
94 images=[]
95 if(type(imagename)==str):
96 images=[imagename]
97 elif(type(imagename)==list):
98 images=imagename
99 for ima in images:
100 for postfix in ['.model', '.image', '.residual', '.mask'] :
101 diskim=ima+postfix
102 if(os.path.exists(diskim)):
103 tb.open(diskim)
104 if(tb.ismultiused()):
105 retval.append(diskim)
106 tb.done()
107 return retval
110 def initsinglems(self, imtool, vis, usescratch):
111 """
112 initialization for single ms case
113 """
114 self.im=imtool
115 # modified for self.vis to be a list for handling multims
116 #self.vis=vis
117 self.vis=[vis]
118 self.sortedvisindx=[0]
120 if ((type(vis)==str) & (os.path.exists(vis))):
121 self.im.open(vis, usescratch=usescratch)
122 else:
123 raise Exception('Visibility data set not found - please verify the name')
124 self.phasecenter=''
125 self.spwindex=-1
126 self.fieldindex=-1
127 self.outputmask=''
128 self.csys={}
130 def initmultims(self, imtool, vislist, usescratch):
131 """
132 initialization for multi-mses
133 """
134 self.im=imtool
135 self.vis=vislist
136 if type(vislist)==list:
137 # self.im.selectvis(vis)
138 if ((type(self.vis[0])==str) & (os.path.exists(self.vis[0]))):
139 pass
140 else:
141 raise Exception('Visibility data set not found - please verify the name')
142 self.phasecenter=''
143 self.spwindex=-1
144 self.fieldindex=-1
145 self.outputmask=''
146 self.csys={}
148 def sortvislist(self,spw,mode,width):
149 """
150 sorting user input vis if it is multiple MSes based on
151 frequencies (spw). So that in sebsequent processing such
152 as setChannelization() correctly works.
153 Returns sorted vis. name list
154 """
155 import operator
156 reverse=False
157 if mode=='velocity':
158 if qa.quantity(width)['value']>0:
159 reverse=True
160 elif mode=='frequency':
161 if qa.quantity(width)['value']<0:
162 reverse=True
163 minmaxfs = []
164 # get selection
165 if len(self.vis) > 1:
166 for i in range(len(self.vis)):
167 visname = self.vis[i]
168 if type(spw)==list and len(spw) > 1:
169 inspw=spw[i]
170 else:
171 inspw=spw
172 if len(inspw)==0:
173 # empty string = select all (='*', for msselectindex)
174 inspw='*'
175 mssel=ms.msseltoindex(vis=visname,spw=inspw)
176 spectable=self.getsubtable(visname, "SPECTRAL_WINDOW")
177 tb.open(spectable)
178 chanfreqs=tb.getvarcol('CHAN_FREQ')
179 tb.close()
180 kys = list(chanfreqs.keys())
181 selspws=mssel['spw']
182 # find extreme freq in each selected spw
183 minmax0=0.
184 firstone=True
185 minmaxallspw=0.
186 for chansel in mssel['channel']:
187 if reverse:
188 minmaxf = max(chanfreqs[kys[chansel[0]]][chansel[1]:chansel[2]+1])
189 else:
190 minmaxf = min(chanfreqs[kys[chansel[0]]][chansel[1]:chansel[2]+1])
191 if firstone:
192 minmaxf0=minmaxf
193 firstone=False
194 if reverse:
195 minmaxallspw=max(minmaxf,minmaxf0)
196 else:
197 minmaxallspw=min(minmaxf,minmaxf0)
198 minmaxfs.append(minmaxallspw)
199 self.sortedvisindx = [x for x, y in sorted(enumerate(minmaxfs),
200 key=operator.itemgetter(1),reverse=reverse)]
201 self.sortedvislist = [self.vis[k] for k in self.sortedvisindx]
202 else:
203 self.sortedvisindx=[0]
204 self.sortedvislist=self.vis
205 #casalog.post("sortedvislist=" + self.sortedvislist)
206 #casalog.post("sortedvisindx=" + self.sortedvisindx)
207 return
210 def defineimages(self, imsize, cell, stokes, mode, spw, nchan, start,
211 width, restfreq, field, phasecenter, facets=1, outframe='',
212 veltype='radio'):
213 """
214 Define image parameters -calls im.defineimage.
215 for processing a single field
216 (also called by definemultiimages for multifield)
217 """
218 if((type(cell)==list) and (len(cell)==1)):
219 cell.append(cell[0])
220 elif ((type(cell)==str) or (type(cell)==int) or (type(cell)==float)):
221 cell=[cell, cell]
222 elif (type(cell) != list):
223 raise TypeError("parameter cell %s is not understood" % str(cell))
224 cellx=qa.quantity(cell[0], 'arcsec')
225 celly=qa.quantity(cell[1], 'arcsec')
226 if(cellx['unit']==''):
227 #string with no units given
228 cellx['unit']='arcsec'
229 if(celly['unit']==''):
230 #string with no units given
231 celly['unit']='arcsec'
232 if((type(imsize)==list) and (len(imsize)==1)):
233 imsize.append(imsize[0])
234 elif(type(imsize)==int):
235 imsize=[imsize, imsize]
236 elif(type(imsize) != list):
237 raise TypeError("parameter imsize %s is not understood" % str(imsize))
239 elstart=start
240 if(mode=='frequency'):
241 ##check that start and step have units
242 if(qa.quantity(start)['unit'].find('Hz') < 0):
243 raise TypeError("start parameter %s is not a valid frequency quantity " % str(start))
244 ###make sure we respect outframe
245 if(self.usespecframe != ''):
246 elstart=me.frequency(self.usespecframe, start)
247 if(qa.quantity(width)['unit'].find('Hz') < 0):
248 raise TypeError("width parameter %s is not a valid frequency quantity " % str(width))
249 elif(mode=='velocity'):
250 ##check that start and step have units
251 if(qa.quantity(start)['unit'].find('m/s') < 0):
252 raise TypeError("start parameter %s is not a valid velocity quantity " % str(start))
253 ###make sure we respect outframe
254 if(self.usespecframe != ''):
255 elstart=me.radialvelocity(self.usespecframe, start)
256 if(qa.quantity(width)['unit'].find('m/s') < 0):
257 raise TypeError("width parameter %s is not a valid velocity quantity " % str(width))
258 else:
259 if((type(width) != int) or
260 (type(start) != int)):
261 raise TypeError("start (%s), width (%s) have to be integers with mode %s" % (str(start),str(width),mode))
263 # changes related multims handling added below (TT)
264 # multi-mses are sorted internally (stored in self.sortedvislist and
265 # indices in self.sortedvisindx) in frequency-wise so that first vis
266 # contains lowest/highest frequency. Note: self.vis is in original user input order.
267 #####understand phasecenter
268 if(type(phasecenter)==str):
269 ### blank means take field[0]
270 if (phasecenter==''):
271 fieldoo=field
272 if(fieldoo==''):
273 msmd.open(self.vis[self.sortedvisindx[0]])
274 allfields=msmd.fieldsforintent('*')
275 fieldoo=str(allfields[0]) if(len(allfields) >0) else '0'
276 msmd.done()
277 #phasecenter=int(ms.msseltoindex(self.vis,field=fieldoo)['field'][0])
278 phasecenter=int(ms.msseltoindex(self.vis[self.sortedvisindx[0]],field=fieldoo)['field'][0])
279 msmd.open(self.vis[self.sortedvisindx[0]])
280 phasecenter=msmd.phasecenter(phasecenter)
281 msmd.done()
282 else:
283 tmppc=phasecenter
284 try:
285 #if(len(ms.msseltoindex(self.vis, field=phasecenter)['field']) > 0):
286 # tmppc = int(ms.msseltoindex(self.vis,
287 # field=phasecenter)['field'][0])
288 # to handle multims set to the first ms that matches
289 for i in self.sortedvisindx:
290 try:
291 if(len(ms.msseltoindex(self.vis[i], field=phasecenter)['field']) > 0):
292 tmppc = int(ms.msseltoindex(self.vis[i],
293 field=phasecenter)['field'][0])
294 # casalog.post("tmppc,i=" + tmppc, i)
295 except Exception:
296 pass
298 ##succesful must be string like '0' or 'NGC*'
299 except Exception:
300 ##failed must be a string 'J2000 18h00m00 10d00m00'
301 tmppc = phasecenter
302 phasecenter = tmppc
303 self.phasecenter = phasecenter
304 #casalog.post('cell' + cellx + celly + restfreq)
305 ####understand spw
306 if spw in (-1, '-1', '*', '', ' '):
307 spwindex = -1
308 else:
309 # old, single ms case
310 #spwindex=ms.msseltoindex(self.vis, spw=spw)['spw'].tolist()
311 #if(len(spwindex) == 0):
312 # spwindex = -1
313 #self.spwindex = spwindex
315 # modified for multims hanlding
316 # multims case
317 if len(self.vis)>1:
318 self.spwindex=[]
319 # will set spwindex in datsel
320 spwindex=-1
321 for i in self.sortedvisindx:
322 if type(spw)==list:
323 inspw=spw[i]
324 else:
325 inspw=spw
326 if len(inspw)==0:
327 inspw='*'
328 self.spwindex.append(ms.msseltoindex(self.vis[i],spw=inspw)['spw'].tolist())
329 # single ms
330 else:
331 spwindex=ms.msseltoindex(self.vis[0], spw=spw)['spw'].tolist()
332 if(len(spwindex) == 0):
333 spwindex = -1
334 self.spwindex = spwindex
336 ##end spwindex
338 if self.usespecframe=='':
339 useframe=self.dataspecframe
340 else:
341 useframe=self.usespecframe
343 self.im.defineimage(nx=imsize[0], ny=imsize[1],
344 cellx=cellx, celly=celly,
345 mode=mode, nchan=nchan,
346 start=elstart, step=width,
347 spw=spwindex, stokes=stokes,
348 restfreq=restfreq, outframe=useframe,
349 veltype=veltype, phasecenter=phasecenter,
350 facets=facets)
352 def definemultiimages(self, rootname, imsizes, cell, stokes, mode, spw,
353 nchan, start, width, restfreq, field, phasecenters,
354 names=[], facets=1, outframe='', veltype='radio',
355 makepbim=False, checkpsf=False):
356 """
357 Define image parameters - multiple field version
358 This fucntion set "private" variables (imagelist and imageids),
359 and then calls defineimages for ecah field.
360 """
361 #will do loop in reverse assuming first image is main field
362 if not hasattr(imsizes, '__len__'):
363 imsizes = [imsizes]
364 self.nimages=len(imsizes)
365 #if((len(imsizes)<=2) and ((type(imsizes[0])==int) or
366 # (type(imsizes[0])==long))):
367 if((len(imsizes)<=2) and (numpy.issubdtype(type(imsizes[0]), int))):
368 self.nimages=1
369 if(len(imsizes)==2):
370 imsizes=[(imsizes[0], imsizes[1])]
371 else:
372 imsizes=[(imsizes[0], imsizes[0])]
373 self._casalog.post('Number of images: ' + str(self.nimages), 'DEBUG1')
375 #imagelist is to have the list of image model names
376 self.imagelist={}
377 #imageids is to tag image to mask in aipsbox style file
378 self.imageids={}
379 if(type(phasecenters) == str):
380 phasecenters=[phasecenters]
381 if(type(phasecenters) == int or type(phasecenters) == float ):
382 phasecenters=[phasecenters]
383 self._casalog.post('Number of phase centers: ' + str(len(phasecenters)),
384 'DEBUG1')
386 if((self.nimages==1) and (type(names)==str)):
387 names=[names]
388 if((len(phasecenters)) != (len(imsizes))):
389 errmsg = "Mismatch between the number of phasecenters (%d), image sizes (%d) , and images (%d)" % (len(phasecenters), len(imsizes), self.nimages)
390 self._casalog.post(errmsg, 'SEVERE')
391 raise ValueError(errmsg)
392 self.skipclean=False
393 lerange=range(self.nimages)
394 lerange.reverse()
395 for n in lerange:
396 self.defineimages(list(imsizes[n]), cell, stokes, mode, spw, nchan,
397 start, width, restfreq, field, phasecenters[n],
398 facets,outframe,veltype)
399 if(len(names)==self.nimages):
400 self.imageids[n]=names[n]
401 if(rootname != ''):
402 self.imagelist[n]=rootname+'_'+names[n]
403 else:
404 self.imagelist[n]=names[n]
405 else:
406 self.imagelist[n]=rootname+'_'+str(n)
407 ###make the image only if it does not exits
408 ###otherwise i guess you want to continue a clean
409 if(not os.path.exists(self.imagelist[n])):
410 self.im.make(self.imagelist[n])
411 #if(makepbim and n==0):
412 if(makepbim):
413 ##make .flux image
414 # for now just make for a main field
415 ###need to get the pointing so select the fields
416 # single ms
417# if len(self.vis)==1:
418# self.im.selectvis(field=field,spw=spw)
419# # multi-ms
420# else:
421# if len(self.vis) > 1:
422# # multi-mses case: use first vis that has the specified field
423# # (use unsorted vis list)
424# nvis=len(self.vis)
425# for i in range(nvis):
426# #if type(field)!=list:
427# # field=[field]
428# try:
429# selparam=self._selectlistinputs(nvis,i,self.paramlist)
430# self.im.selectvis(vis=self.vis[i],field=selparam['field'],spw=selparam['spw'])
431# except:
432# pass
434 # set to default minpb(=0.1), should use input minpb?
435 self.im.setmfcontrol()
436 self.im.setvp(dovp=True)
437 self.im.makeimage(type='pb', image=self.imagelist[n]+'.flux',
438 compleximage="", verbose=False)
439 self.im.setvp(dovp=False, verbose=False)
442 def checkpsf(self,chan):
443 """
444 a check to make sure selected channel plane is not entirely flagged
445 (for chinter=T interactive clean)
446 """
447 #lerange=range(self.nimages)
448 #lerange.reverse()
449 #for n in lerange:
450 # self.getchanimage(self.finalimages[n]+'.psf',self.imagelist[n]+'.test.psf',chan)
451 # ia.open(self.imagelist[n]+'.test.psf')
452 # imdata=ia.getchunk()
453 # casalog.post("imagelist[", n, "]=" + self.imagelist[n] + " imdata.sum()=" + imdata.sum())
454 #if n==0 and imdata.sum()==0.0:
455 # self.skipclean=True
456 # if self.skipclean:
457 # pass
458 # elif imdata.sum()==0.0:
459 # self.skipclean=True
461 # need to check only for main field
462 self.getchanimage(self.finalimages[0]+'.psf',self.imagelist[0]+'.test.psf',chan)
463 ia.open(self.imagelist[0]+'.test.psf')
464 imdata=ia.getchunk()
465 if imdata.sum()==0.0:
466 self.skipclean=True
469 def makeEmptyimages(self):
470 """
471 Create empty images (0.0 pixel values) for
472 image, residual, psf
473 must run after definemultiimages()
474 and it is assumed that definemultiimages creates
475 empty images (self.imagelist).
476 """
477 lerange=range(self.nimages)
478 for n in lerange:
479 os.system('cp -r '+self.imagelist[n]+' '+self.imagelist[n]+'.image')
480 os.system('cp -r '+self.imagelist[n]+' '+self.imagelist[n]+'.residual')
481 os.system('cp -r '+self.imagelist[n]+' '+self.imagelist[n]+'.psf')
482 os.system('cp -r '+self.imagelist[n]+' '+self.imagelist[n]+'.model')
483 os.system('cp -r '+self.imagelist[n]+' '+self.imagelist[n]+'.mask')
486 def makemultifieldmask(self, maskobject=''):
487 """
488 This function assumes that the function definemultiimages has been run and thus
489 self.imagelist is defined
490 if single image use the single image version
491 (this is not up to date for current mask handling but used in task_widefield,
492 to be removed after task_widefield is gone)
493 """
494 if((len(self.maskimages)==(len(self.imagelist)))):
495 if(self.imagelist[0] not in self.maskimages):
496 self.maskimages={}
497 else:
498 self.maskimages={}
499 masktext=[]
500 if (not hasattr(maskobject, '__len__')) \
501 or (len(maskobject) == 0) or (maskobject == ['']):
502 return
503 if(type(maskobject)==str):
504 maskobject=[maskobject]
505 if(type(maskobject) != list):
506 ##don't know what to do with this
507 raise TypeError('Dont know how to deal with mask object')
508 n=0
509 for masklets in maskobject:
510 if(type(masklets)==str):
511 if(os.path.exists(masklets)):
512 if(subprocess.getoutput('file '+masklets).count('directory')):
513 self.maskimages[self.imagelist[n]]=masklets
514 n=n+1
515 elif(subprocess.getoutput('file '+masklets).count('text')):
516 masktext.append(masklets)
517 else:
518 raise TypeError('Can only read text mask files or mask images')
519 else:
520 raise TypeError(masklets+' seems to be non-existant')
521 if(len(masktext) > 0):
522 circles, boxes=self.readmultifieldboxfile(masktext)
523 if(len(self.maskimages)==0):
524 for k in range(len(self.imageids)):
525 if(self.imagelist[k] not in self.maskimages):
526 self.maskimages[self.imagelist[k]]=self.imagelist[k]+'.mask'
527 for k in range(len(self.imageids)):
528 ###initialize mask if its not there yet
529 if(not (os.path.exists(self.maskimages[self.imagelist[k]]))):
530 ia.fromimage(outfile=self.maskimages[self.imagelist[k]],
531 infile=self.imagelist[k])
532 ia.open(self.maskimages[self.imagelist[k]])
533 ia.set(pixels=0.0)
534 ia.done(verbose=False)
535 if(self.imageids[k] in circles and self.imageids[k] in boxes):
536 self.im.regiontoimagemask(mask=self.maskimages[self.imagelist[k]],
537 boxes=boxes[self.imageids[k]],
538 circles=circles[self.imageids[k]])
539 elif(self.imageids[k] in circles):
540 self.im.regiontoimagemask(mask=self.maskimages[self.imagelist[k]],
541 circles=circles[self.imageids[k]])
542 elif(self.imageids[k] in boxes):
543 self.im.regiontoimagemask(mask=self.maskimages[self.imagelist[k]],
544 boxes=boxes[self.imageids[k]])
545 else:
546 ###need to make masks that select that whole image
547 ia.open(self.maskimages[self.imagelist[k]])
548 ia.set(pixels=1.0)
549 ia.done(verbose=False)
552 def makemultifieldmask2(self, maskobject='',slice=-1, newformat=True, interactive=False):
554 """
555 Create mask images for multiple fields (flanking fields)
556 This new makemultifieldmask to accomodate different kinds of masks supported
557 in clean with flanking fields.
559 Keyword arguments:
560 maskobject -- input mask, a list (of string or of list(s)) or a string
561 slice -- channel slice (to handle chaniter mode): default = -1 (all)
562 newformat -- if mask is read from new format text file: default = True
564 Prerequiste: definemultiimages has already ran so that imagelist is defined
566 Notes:
567 * It makes empty mask images at begenning, calls makemaskimage, and if no
568 mask to be specified, the corresponding empty mask image is removed
570 * When clean is executed in commnad line style (e.g. clean(vis='..', ..))
571 it is possible to have mask parameter consists of a mix of strings and int lists
572 (i.e. mask=[['newreg.txt',[55,55,65,70]],['outliermask.rgn']]), and this function
573 should be able to parse these properly. - although this won't work for the task execution
574 by go() and tends to messed up inp() after such execution
576 * Currently it is made to handle old outlier text file format and boxfile-style
577 mask box specification for backward compartibility. But it is planned to
578 be removed for 3.4.
580 * This is a refactored version of the previous makemultifieldmask2
581 and calls makemaskimage() for each field.
582 this was called makemultifieldmask3 in CASA 3.3 release but now renamed
583 makemultifieldmask2 as the previous makemultifieldmask2 was removed.
584 """
585 #casalog.post("Inside makemultifieldmask2")
586 if((len(self.maskimages)==(len(self.imagelist)))):
587 if(self.imagelist[0] not in self.maskimages):
588 self.maskimages=odict()
589 else:
590 self.maskimages=odict()
591 # clean up temp mask image
592 if os.path.exists('__tmp_mask'):
593 shutil.rmtree('__tmp_mask')
595 if (not hasattr(maskobject, '__len__')) \
596 or (len(maskobject) == 0) or (maskobject == ['']):
597 return
598 # for empty maskobject list
599 if all([msk==[''] or msk==[] for msk in maskobject]):
600 return
601 # determine number of input elements
602 if (type(maskobject)==str):
603 maskobject=[maskobject]
604 if(type(maskobject) != list):
605 ##don't know what to do with this
606 raise TypeError('Dont know how to deal with maskobject with type: %s' % type(maskobject))
607 #if(type(maskobject[0])==int or type(maskobject[0])==float):
608 if(numpy.issubdtype(type(maskobject[0]),int) or numpy.issubdtype(type(maskobject[0]),float)):
609 maskobject=[maskobject]
610 if(type(maskobject[0][0])==list):
611 #if(type(maskobject[0][0][0])!=int and type(maskobject[0][0][0])!=float):
612 if not (numpy.issubdtype(type(maskobject[0][0][0]),int) or \
613 numpy.issubdtype(type(maskobject[0][0][0]),float)):
614 maskobject=maskobject[0]
616 # define maskimages
617 if(len(self.maskimages)==0):
618 for k in range(len(self.imageids)):
619 if(self.imagelist[k] not in self.maskimages):
620 self.maskimages[self.imagelist[k]]=self.imagelist[k]+'.mask'
621 # initialize maskimages - create empty maskimages
622 # --- use outframe or dataframe for mask creation
623 # * It appears somewhat duplicating with makemaskimage
624 # but it is necessary to create a maskimage for
625 # each field at this point...
626 if self.usespecframe=='':
627 maskframe=self.dataspecframe
628 else:
629 maskframe=self.usespecframe
630 #casalog.post("Frame : " + maskframe)
631 #casalog.post("dataframe : " + self.dataspecframe + " useframe : " + self.usespecframe)
632 for k in range(len(self.imagelist)):
633 if(not (os.path.exists(self.maskimages[self.imagelist[k]]))):
634 ia.fromimage(outfile=self.maskimages[self.imagelist[k]],
635 infile=self.imagelist[k])
636 ia.open(self.maskimages[self.imagelist[k]])
637 ia.set(pixels=0.0)
638 #mcsys=ia.coordsys().torecord()
639 #if mcsys['spectral2']['conversion']['system']!=maskframe:
640 # mcsys['spectral2']['conversion']['system']=maskframe
641 #ia.setcoordsys(mcsys)
642 #
643 ## This code to set the maskframe is copied from makemaskimages()
644# mycsys=ia.coordsys()
645# if mycsys.torecord()['spectral2']['conversion']['system']!=maskframe:
646# mycsys.setreferencecode(maskframe,'spectral',True)
647# self.csys=mycsys.torecord()
648# if self.csys['spectral2']['conversion']['system']!=maskframe:
649# self.csys['spectral2']['conversion']['system']=maskframe
650# ia.setcoordsys(self.csys)
652 ia.done(verbose=False)
654 self.setReferenceFrameLSRK(img =self.maskimages[self.imagelist[k]])
656 # take out extra []'s
657 maskobject=self.flatten(maskobject)
658 masktext=[]
659 # to keep backward compatibility for a mixed outlier file
660 # look for boxfiles contains multiple boxes with image ids
661 for maskid in range(len(maskobject)):
662 if type(maskobject[maskid])==str:
663 maskobject[maskid] = [maskobject[maskid]]
665 for masklets in maskobject[maskid]:
666 if type(masklets)==str:
667 if (os.path.exists(masklets) and
668 (not subprocess.getoutput('file '+masklets).count('directory')) and
669 (not subprocess.getoutput('file '+masklets).split(':')[-1].count('data'))):
670 # extract boxfile name
671 masktext.append(masklets)
673 # === boxfile handling ===
674 #extract boxfile mask info only for now the rest is
675 #processed by makemaskimage. - DEPRECATED and will be removed
676 #in 3.4
677 #
678 #group circles and boxes in dic for each image field
679 circles, boxes, oldfmts=self.readmultifieldboxfile(masktext)
681 # Loop over imagename
682 # take out text file names contain multifield boxes and field info
683 # from maskobject and create updated one (updatedmaskobject)
684 # by adding boxlists to it instead.
685 # Use readmultifieldboxfile to read old outlier/box text file format
686 # Note: self.imageids and boxes's key are self.imagelist for new outlier
687 # format while for the old format, they are 'index' in string.
689 #maskobject_tmp = maskobject
690 # need to do a deep copy
691 import copy
692 maskobject_tmp=copy.deepcopy(maskobject)
693 updatedmaskobject = []
694 for maskid in range(len(maskobject_tmp)):
695 if len(circles)!=0 or len(boxes)!=0:
696 # remove the boxfiles from maskobject list
697 for txtf in masktext:
698 if maskobject_tmp[maskid].count(txtf) and oldfmts[txtf]:
699 maskobject_tmp[maskid].remove(txtf)
700 updatedmaskobject = maskobject_tmp
701 else:
702 updatedmaskobject = maskobject
703 # adjust no. of elements of maskoject list with []
704 if len(updatedmaskobject)-len(self.imagelist)<0:
705 for k in range(len(self.imagelist)-len(updatedmaskobject)):
706 updatedmaskobject.append([])
707 #for maskid in range(len(self.maskimages)):
709 # === boxfile handling ====
710 for maskid in self.maskimages.keys():
711 # for handling old format
712 #if nmaskobj <= maskid:
713 # add circles,boxes back
714 maskindx = self.maskimages.keys().index(maskid)
715 if len(circles) != 0:
716 for key in circles:
717 if (newformat and maskid==key) or \
718 (not newformat and maskid.split('_')[-1]==key):
719 if len(circles[key])==1:
720 incircles=circles[key][0]
721 else:
722 incircles=circles[key]
723 # put in imagelist order
724 if len(incircles)>1 and isinstance(incircles[0],list):
725 updatedmaskobject[self.imagelist.values().index(maskid)].extend(incircles)
726 else:
727 updatedmaskobject[self.imagelist.values().index(maskid)].append(incircles)
728 if len(boxes) != 0:
729 for key in boxes:
730 #try:
731 # keyid=int(key)
732 #except:
733 # keyid=key
734 if (newformat and maskid==key) or \
735 (not newformat and maskid.split('_')[-1]==key):
736 if len(boxes[key])==1:
737 inboxes=boxes[key][0]
738 else:
739 inboxes=boxes[key]
740 # add to maskobject (extra list bracket taken out)
741 # put in imagelist order
742 # take out extra []
743 if len(inboxes)>1 and isinstance(inboxes[0],list):
744 updatedmaskobject[self.imagelist.values().index(maskid)].extend(inboxes)
745 else:
746 updatedmaskobject[self.imagelist.values().index(maskid)].append(inboxes)
747 # === boxfile handling ====
749 # do maskimage creation (call makemaskimage)
750 for maskid in range(len(self.maskimages)):
751 if maskid < len(updatedmaskobject):
752 self._casalog.post("Matched masks: maskid=%s mask=%s" % (maskid, updatedmaskobject[maskid]), 'DEBUG1')
753 self.outputmask=''
754 self.makemaskimage(outputmask=self.maskimages[self.imagelist[maskid]],
755 imagename=self.imagelist[maskid], maskobject=updatedmaskobject[maskid], slice=slice)
756# self._casalog.post("Matched masks: maskid=%s mask=%s" % (maskid, updatedmaskobject[maskid]), 'DEBUG1')
757# self.outputmask=''
758# self.makemaskimage(outputmask=self.maskimages[self.imagelist[maskid]],
759# imagename=self.imagelist[maskid], maskobject=updatedmaskobject[maskid], slice=slice)
761 for key in self.maskimages.keys():
762 if(os.path.exists(self.maskimages[key])):
763 ia.open(self.maskimages[key])
764 fsum=ia.statistics(verbose=False,list=False)['sum']
765 if(len(fsum)!=0 and fsum[0]==0.0):
766 # make an empty mask
767 ia.set(pixels=0.0)
768 # should not remove empty mask for multifield case
769 # interactive=F.
770 # Otherwise makemaskimage later does not work
771 # remove the empty mask
772 if not interactive:
773 ia.remove()
774 ia.done(verbose=False)
777 def make_mask_from_threshhold(self, imagename, thresh, outputmask=None):
778 """
779 Makes a mask image with the same coords as imagename where each
780 pixel is True if and only if the corresponding pixel in imagename
781 is >= thresh.
783 The mask will be named outputmask (if provided) or imagename +
784 '.thresh_mask'. The name is returned on success, or False on failure.
785 """
786 if not outputmask:
787 outputmask = imagename + '.thresh_mask'
789 # im.mask would be a lot shorter, but it (unnecessarily) needs im to be
790 # open with an MS.
791 # I am not convinced that im.mask should really be using Quantity.
792 # qa.quantity(quantity) = quantity.
793 self.im.mask(imagename, outputmask, qa.quantity(thresh))
795 ## # Copy imagename to a safe name to avoid problems with /, +, -, and ia.
796 ## ia.open(imagename)
797 ## shp = ia.shape()
798 ## ia.close()
799 ## self.copymaskimage(imagename, shp, '__temp_mask')
801 ## self.copymaskimage(imagename, shp, outputmask)
802 ## ia.open(outputmask)
803 ## ###getchunk is a mem hog
804 ## #arr=ia.getchunk()
805 ## #arr[arr>0.01]=1
806 ## #ia.putchunk(arr)
807 ## #inpix="iif("+"'"+outputmask.replace('/','\/')+"'"+">0.01, 1, 0)"
808 ## #ia.calc(pixels=inpix)
809 ## ia.calc(pixels="iif(__temp_mask>" + str(thresh) + ", 1, 0)")
810 ## ia.close()
811 ## ia.removefile('__temp_mask')
812 return outputmask
814 def makemaskimage(self, outputmask='', imagename='', maskobject=[], slice=-1):
815 """
816 This function is an attempt to convert all the kind of 'masks' that
817 people want to throw at it and convert it to a mask image to be used
818 by imager...For now 'masks' include
820 a)set of previous mask images
821 b)lists of blc trc's
822 c)record output from rg tool for e.g
824 * for a single field
825 """
826 if (not hasattr(maskobject, '__len__')) \
827 or (len(maskobject) == 0) or (maskobject == ['']):
828 return
829 maskimage=[]
830 masklist=[]
831 textreglist=[]
832 masktext=[]
833 maskrecord={}
834 tablerecord=[]
835 # clean up any left over temp files from previous clean runs
836 if os.path.exists("__temp_mask"):
837 shutil.rmtree("__temp_mask")
838 if os.path.exists("__temp_mask2"):
839 shutil.rmtree("__temp_mask2")
841 # relax to allow list input for imagename
842 if(type(imagename)==list):
843 imagename=imagename[0]
845 if(type(maskobject)==dict):
846 maskrecord=maskobject
847 maskobject=[]
848 if(type(maskobject)==str):
849 maskobject=[maskobject]
851 if(type(maskobject) != list):
852 ##don't know what to do with this
853 raise TypeError('Dont know how to deal with maskobject')
854 if(numpy.issubdtype(type(maskobject[0]),numpy.int) or \
855 numpy.issubdtype(type(maskobject[0]),numpy.float)):
856 # check and convert if list consist of python int or float
857 maskobject_tmp = convert_numpydtype(maskobject)
858 masklist.append(maskobject_tmp)
859 else:
860 for masklets in maskobject:
861 if(type(masklets)==str): ## Can be a file name, or an explicit region-string
862 if(os.path.exists(masklets)):
863 if(subprocess.getoutput('file '+masklets).count('directory')):
864 maskimage.append(masklets)
865 elif(subprocess.getoutput('file '+masklets).count('text')):
866 masktext.append(masklets)
867 else:
868 tablerecord.append(masklets)
869 else:
870 textreglist.append(masklets);
871 #raise TypeError, masklets+' seems to be non-existant'
872 if(type(masklets)==list):
873 masklets_tmp = convert_numpydtype(masklets)
874 masklist.append(masklets_tmp)
875 if(type(masklets)==dict):
876 maskrecord=masklets
877 if(len(outputmask)==0):
878 outputmask=imagename+'.mask'
879 if(os.path.exists(outputmask)):
880 # for multiple field
881 # outputmask is always already defined
882 # cannot use copymaskiamge since self.csys used in the code
883 # fixed to that of main field
884 if len(self.imagelist)>1:
885 ia.fromimage('__temp_mask',outputmask,overwrite=True)
886 ia.close()
887 else:
888 self.im.make('__temp_mask')
889 ia.open('__temp_mask')
890 shp=ia.shape()
891 self.csys=ia.coordsys().torecord()
892 ia.close()
893 ia.removefile('__temp_mask')
894 ia.open(outputmask)
895 outim = ia.regrid(outfile='__temp_mask',shape=shp,axes=[3,0,1], csys=self.csys,overwrite=True, asvelocity=False)
896 outim.done(verbose=False)
897 ia.done(verbose=False)
898 ia.removefile(outputmask)
899 os.rename('__temp_mask',outputmask)
900 else:
901 self.im.make(outputmask)
902 if len(self.imagelist)>1:
903 raise Exception("Multifield case - requires initial mask images but undefined.")
905 # respect dataframe or outframe
906 if self.usespecframe=='':
907 maskframe=self.dataspecframe
908 else:
909 maskframe=self.usespecframe
910 if len(self.vis)!=1:
911 if not self.inframe:
912 # for multi-ms case default output frame is default to LSRK
913 # (set by baseframe in imager_cmt.cc)
914 maskframe='LSRK'
916 ia.open(outputmask)
917 # make output mask template unit less to avoid imagenalysis warning...
918 if (ia.brightnessunit()!=""):
919 ia.setbrightnessunit("")
920 shp=ia.shape()
921 self.csys=ia.coordsys().torecord()
922 # keep this info for reading worldbox
923 self.csysorder=ia.coordsys().coordinatetype()
924# ia.close()
926# self.setReferenceFrameLSRK( outputmask )
928 mycsys=ia.coordsys()
929 if mycsys.torecord()['spectral2']['conversion']['system']!=maskframe:
930 mycsys.setreferencecode(maskframe,'spectral',True)
931 self.csys=mycsys.torecord()
932 if self.csys['spectral2']['conversion']['system']!=maskframe:
933 self.csys['spectral2']['conversion']['system']=maskframe
934 ia.setcoordsys(self.csys)
935 #ia.setcoordsys(mycsys.torecord())
936 ia.close() # close outputmask
938 if(len(maskimage) > 0):
939 for ima in maskimage :
940 ia.open(ima)
941 maskshape = ia.shape()
942 ia.close()
943 if shp[3]>1 and maskshape[3]==1:
944 self._casalog.post("Input mask,"+ima+" appears to be a continuum,"+
945 " it may not map to a cube correctly and will be"+
946 " ignored."+" Please use "+
947 "makemask to generate a proper mask.","WARN")
948 if slice>-1:
949 self.getchanimage(ima, ima+'chanim',slice)
950 self.copymaskimage(ima+'chanim',shp,'__temp_mask')
951 ia.removefile(ima+'chanim')
952 else:
953 self.copymaskimage(ima, shp, '__temp_mask')
954 #ia.open(ima)
955 #ia.regrid(outfile='__temp_mask',shape=shp,csys=self.csys,
956 # overwrite=True)
957 #ia.done(verbose=False)
958 os.rename(outputmask,'__temp_mask2')
959 outim = ia.imagecalc(outfile=outputmask,
960 pixels='__temp_mask + __temp_mask2',
961 overwrite=True)
962 outim.done(verbose=False)
963 ia.done(verbose=False)
964 ia.removefile('__temp_mask')
965 ia.removefile('__temp_mask2')
966 if(not os.path.exists(outputmask)):
967 outputmask = self.make_mask_from_threshhold(outputmask, 0.01,
968 outputmask)
969 #pdb.set_trace()
970 #### This goes when those tablerecord goes
971 ### Make masks from tablerecords
972 if(len(tablerecord) > 0):
973 reg={}
974 for tabl in tablerecord:
975 try:
976 reg.update({tabl:rg.fromfiletorecord(filename=tabl, verbose=False)})
977 except:
978 raise Exception('Region-file (binary) format not recognized. Please check. If box-file, please start the file with \'#boxfile\' on the first line')
979 if(len(reg)==1):
980 reg=reg[reg.keys()[0]]
981 else:
982 reg=rg.makeunion(reg)
983 self.im.regiontoimagemask(mask=outputmask, region=reg)
984 ###############
985 ### Make masks from region dictionaries
986 if((type(maskrecord)==dict) and (len(maskrecord) > 0)):
987 self.im.regiontoimagemask(mask=outputmask, region=maskrecord)
988 ### Make masks from text files
989 if(len(masktext) >0):
990 for textfile in masktext :
991 # Read a box file
992 polydic,listbox=self.readboxfile(textfile);
994 masklist.extend(listbox)
995 if(len(polydic) > 0):
996 ia.open(outputmask)
997 ia.close()
998 self.im.regiontoimagemask(mask=outputmask, region=polydic)
999 # If box lists are empty, it may be a region format
1000 if(len(polydic)==0 and len(listbox)==0):
1001 # Read in a region file
1002 try:
1003 ia.open(outputmask);
1004 mcsys = ia.coordsys();
1005 mshp = ia.shape();
1006 ia.close();
1007 mreg = rg.fromtextfile(filename=textfile,shape=mshp,csys=mcsys.torecord());
1008 self.im.regiontoimagemask(mask=outputmask, region=mreg);
1009 except:
1010 raise Exception('Region-file (text) format not recognized. Please check. If box-file, please start the file with \'#boxfile\' on the first line, and have at-least one valid box in it')
1011 ### Make masks from inline lists of pixel coordinates
1012 if((type(masklist)==list) and (len(masklist) > 0)):
1013 self.im.regiontoimagemask(mask=outputmask, boxes=masklist)
1014 ### Make masks from inline region-strings
1015 if((type(textreglist)==list) and (len(textreglist)>0)):
1016 ia.open(outputmask);
1017 mcsys = ia.coordsys();
1018 mshp = ia.shape();
1019 ia.close();
1020 for textlet in textreglist:
1021 try:
1022 mreg = rg.fromtext(text=textlet,shape=mshp,csys=mcsys.torecord());
1023 self.im.regiontoimagemask(mask=outputmask, region=mreg);
1024 except:
1025 raise Exception('\''+textlet+'\' is not recognized as a text file on disk or as a region string')
1026 ### Make mask from an image-mask
1027 if(os.path.exists(imagename) and (len(rg.namesintable(imagename)) !=0)):
1028 regs=rg.namesintable(imagename)
1029 if(type(regs)==str):
1030 regs=[regs]
1031 for reg in regs:
1032 elrec=rg.fromtabletorecord(imagename, reg)
1033 self.im.regiontoimagemask(mask=outputmask, region=elrec)
1035 self.outputmask=outputmask
1037 ## CAS-5227
1038 ia.open( outputmask )
1039 ia.calc('iif("'+outputmask+'"!=0.0,1.0,0.0)', False)
1040 ia.close()
1042 ## CAS-5221
1043 #### CAS-6676 : Force frame to LSRK only if it's a fresh image being made.
1044 #### If residual exists, then it's likely to not be in LSRK, so don't force the mask to be it.
1045 #### The call to setFrameConversionForMasks() from task_clean.py would have set the
1046 #### mask frame to match the ouput frame in the previous run.
1047 if not os.path.exists(imagename+'.residual'):
1048 self.setReferenceFrameLSRK( outputmask )
1049 #Done with making masks
1052 def datselweightfilter(self, field, spw, timerange, uvrange, antenna,scan,
1053 wgttype, robust, noise, npixels, mosweight,
1054 innertaper, outertaper, usescratch, nchan=-1, start=0, width=1):
1055 """
1056 Make data selection
1057 (not in use, split into datsel and datweightfileter)
1058 """
1059 rmode='none'
1060 weighting='natural';
1061 if(wgttype=='briggsabs'):
1062 weighting='briggs'
1063 rmode='abs'
1064 elif(wgttype=='briggs'):
1065 weighting='briggs'
1066 rmode='norm'
1067 else:
1068 weighting=wgttype
1070 self.fieldindex=ms.msseltoindex(self.vis,field=field)['field'].tolist()
1071 if(len(self.fieldindex)==0):
1072 fieldtab=self.getsubtable(self.vis, 'FIELD')
1073 tb.open(fieldtab)
1074 self.fieldindex=range(tb.nrows())
1075 tb.close()
1076 #weighting and tapering should be done together
1077 if(weighting=='natural'):
1078 mosweight=False
1079 self.im.selectvis(nchan=nchan,start=start,step=width,field=field,spw=spw,time=timerange,
1080 baseline=antenna, scan=scan, uvrange=uvrange, usescratch=usescratch)
1081 self.im.weight(type=weighting,rmode=rmode,robust=robust, npixels=npixels, noise=qa.quantity(noise,'Jy'), mosaic=mosweight)
1082 if((type(outertaper)==list) and (len(outertaper) > 0)):
1083 if(len(outertaper)==1):
1084 outertaper.append(outertaper[0])
1085 outertaper.append('0deg')
1086 if(qa.quantity(outertaper[0])['value'] > 0.0):
1087 self.im.filter(type='gaussian', bmaj=outertaper[0],
1088 bmin=outertaper[1], bpa=outertaper[2])
1091 # split version of datselweightfilter
1092 def datsel(self, field, spw, timerange, uvrange, antenna, scan, observation,intent,
1093 usescratch, nchan=-1, start=0, width=1):
1094 """
1095 Make selections in visibility data
1096 """
1098 # for multi-MSes, if field,spw,timerage,uvrange,antenna,scan are not
1099 # lists the same selection is applied to all the MSes.
1100 self.fieldindex=[]
1101 #nvislist=range(len(self.vis))
1102 vislist=self.sortedvisindx
1103 self.paramlist={'field':field,'spw':spw,'timerange':timerange,'antenna':antenna,
1104 'scan':scan, 'observation': observation, 'intent':intent, 'uvrange':uvrange}
1105 for i in vislist:
1106 selectedparams=self._selectlistinputs(len(vislist),i,self.paramlist)
1107 tempfield=selectedparams['field']
1108# if type(field)==list:
1109# if len(field)==nvislist:
1110# tempfield=field[i]
1111# else:
1112# if len(field)==1:
1113# tempfield=field[0]
1114# else:
1115# raise Exception, 'The number of field list does not match with the number of vis list'
1116# else:
1117# tempfield=field
1119 if len(tempfield)==0:
1120 tempfield='*'
1121 self.fieldindex.append(ms.msseltoindex(self.vis[i],field=tempfield)['field'].tolist())
1123 ############################################################
1124 # Not sure I need this now.... Nov 15, 2010
1125 vislist.reverse()
1126 writeaccess=True
1127 for i in vislist:
1128 writeaccess=writeaccess and os.access(self.vis[i], os.W_OK)
1129 #if any ms is readonly then no model will be stored, MSs will be in readmode only...but clean can proceed
1130 for i in vislist:
1131 # select apropriate parameters
1132 selectedparams=self._selectlistinputs(len(vislist),i,self.paramlist)
1133 inspw=selectedparams['spw']
1134 intimerange=selectedparams['timerange']
1135 inantenna=selectedparams['antenna']
1136 inscan=selectedparams['scan']
1137 inobs = selectedparams['observation']
1138 inintent = selectedparams['intent']
1139 inuvrange=selectedparams['uvrange']
1141 #if len(self.vis)==1:
1142 #casalog.post("single ms case")
1143 # self.im.selectvis(nchan=nchan,start=start,step=width,field=field,
1144 # spw=inspw,time=intimerange, baseline=inantenna,
1145 # scan=inscan, observation=inobs, intent=inintent, uvrange=inuvrange,
1146 # usescratch=usescratch)
1147 #else:
1148 #casalog.post("multims case: selectvis for vis[",i,"]: spw,field=", inspw, self.fieldindex[i])
1149 self.im.selectvis(vis=self.vis[i],nchan=nchan,start=start,step=width,
1150 field=self.fieldindex[i], spw=inspw,time=intimerange,
1151 baseline=inantenna, scan=inscan,
1152 observation=inobs, intent=inintent,
1153 uvrange=inuvrange, usescratch=usescratch, writeaccess=writeaccess)
1155 # private function for datsel and datweightfilter
1156 def _selectlistinputs(self,nvis,indx,params):
1157 """
1158 A little private function to do selection and checking for a parameter
1159 given in list of strings.
1160 It checks nelement in each param either match with nvis or nelement=1
1161 (or a string) otherwise exception is thrown.
1162 """
1163 outparams={}
1164 if type(params)==dict:
1165 for param,val in params.items():
1166 msg = 'The number of %s list given in list does not match with the number of vis list given.' % param
1167 if type(val)==list:
1168 if len(val)==nvis:
1169 outval=val[indx]
1170 else:
1171 if len(val)==1:
1172 outval=val[0]
1173 else:
1174 raise Exception(msg)
1175 outparams[param]=outval
1176 else:
1177 #has to be a string
1178 outparams[param]=val
1179 return outparams
1180 else:
1181 raise Exception('params must be a dictionary')
1183 # weighting/filtering part of datselweightfilter.
1184 # The scan parameter is not actually used, so observation is not included
1185 # as a parameter. Both are used via self._selectlistinputs().
1186 def datweightfilter(self, field, spw, timerange, uvrange, antenna,scan,
1187 wgttype, robust, noise, npixels, mosweight,
1188 uvtaper,innertaper, outertaper, usescratch, nchan=-1, start=0, width=1):
1189 """
1190 Apply weighting and tapering
1191 """
1192 rmode='none'
1193 weighting='natural';
1194 if(wgttype=='briggsabs'):
1195 weighting='briggs'
1196 rmode='abs'
1197 elif(wgttype=='briggs'):
1198 weighting='briggs'
1199 rmode='norm'
1200 else:
1201 weighting=wgttype
1202 #weighting and tapering should be done together
1203 if(weighting=='natural'):
1204 mosweight=False
1205# vislist=self.sortedvisindx
1206 #nvislist.reverse()
1207# for i in vislist:
1208# # select apropriate parameters
1209# selectedparams=self._selectlistinputs(len(vislist),i,self.paramlist)
1210# inspw=selectedparams['spw']
1211# intimerange=selectedparams['timerange']
1212# inantenna=selectedparams['antenna']
1213# inscan=selectedparams['scan']
1214# inobs=selectedparams['observation']
1215# inuvrange=selectedparams['uvrange']
1216#
1217# if len(self.vis) > 1:
1218# casalog.post('from datwtfilter - multi';)
1219# self.im.selectvis(vis=self.vis[i], field=self.fieldindex[i],spw=inspw,time=intimerange,
1220# baseline=inantenna, scan=inscan, observation=inobs,
1221# uvrange=inuvrange, usescratch=calready)
1222# else:
1223# casalog.post('from datwtfilter - single')
1224# self.im.selectvis(field=field,spw=inspw,time=intimerange,
1225# baseline=inantenna, scan=inscan, observation=inobs,
1226# uvrange=inuvrange, usescratch=calready)
1227# self.im.weight(type=weighting,rmode=rmode,robust=robust,
1228# npixels=npixels, noise=qa.quantity(noise,'Jy'), mosaic=mosweight)
1229 self.im.weight(type=weighting,rmode=rmode,robust=robust,
1230 npixels=npixels, noise=qa.quantity(noise,'Jy'), mosaic=mosweight)
1232 oktypes = (str,int,float)
1234 if((uvtaper==True) and (type(outertaper) in oktypes)):
1235 outertaper=[outertaper]
1236 if((uvtaper==True) and (type(outertaper)==list) and (len(outertaper) > 0)):
1237 if(len(outertaper)==1):
1238 outertaper.append(outertaper[0])
1239 if(len(outertaper)==2):
1240 outertaper.append('0deg')
1241 if(qa.quantity(outertaper[0])['unit']==''):
1242 outertaper[0]=qa.quantity(qa.quantity(outertaper[0])['value'],'lambda')
1243 if(qa.quantity(outertaper[1])['unit']==''):
1244 outertaper[1]=qa.quantity(qa.quantity(outertaper[1])['value'],'lambda')
1245 if(qa.quantity(outertaper[0])['value'] > 0.0):
1246 self.im.filter(type='gaussian', bmaj=outertaper[0],
1247 bmin=outertaper[1], bpa=outertaper[2])
1250 def setrestoringbeam(self, restoringbeam):
1251 """
1252 Set restoring beam
1253 """
1254 if((restoringbeam == ['']) or (len(restoringbeam) ==0)):
1255 return
1256 resbmaj=''
1257 resbmin=''
1258 resbpa='0deg'
1259 if((type(restoringbeam) == list) and len(restoringbeam)==1):
1260 restoringbeam=restoringbeam[0]
1261 if((type(restoringbeam)==str)):
1262 if(qa.quantity(restoringbeam)['unit'] == ''):
1263 restoringbeam=restoringbeam+'arcsec'
1264 resbmaj=qa.quantity(restoringbeam, 'arcsec')
1265 resbmin=qa.quantity(restoringbeam, 'arcsec')
1266 if(type(restoringbeam)==list):
1267 resbmaj=qa.quantity(restoringbeam[0], 'arcsec')
1268 resbmin=qa.quantity(restoringbeam[1], 'arcsec')
1269 if(resbmaj['unit']==''):
1270 resbmaj=restoringbeam[0]+'arcsec'
1271 if(resbmin['unit']==''):
1272 resbmin=restoringbeam[1]+'arcsec'
1273 if(len(restoringbeam)==3):
1274 resbpa=qa.quantity(restoringbeam[2], 'deg')
1275 if(resbpa['unit']==''):
1276 resbpa=restoringbeam[2]+'deg'
1277 if((resbmaj != '') and (resbmin != '')):
1278 self.im.setbeam(bmaj=resbmaj, bmin=resbmin, bpa=resbpa)
1280 def convertmodelimage(self, modelimages=[], outputmodel='',imindex=0):
1281 """
1282 Convert model inputs to a model image
1284 Keyword arguments:
1285 modleimages -- input model list
1286 outputmodel -- outout modelimage name
1287 imindex -- image name index (corresponding to imagelist)
1288 for multi field hanlding
1289 """
1290 modelos=[]
1291 maskelos=[]
1292 if((modelimages=='') or (modelimages==[]) or (modelimages==[''])):
1293 #if((modelimages=='') or (modelimages==[])):
1294 return
1295 if(type(modelimages)==str):
1296 modelimages=[modelimages]
1297 k=0
1298 for modim in modelimages:
1299 if not os.path.exists(modim):
1300 raise Exception("Model image file name="+modim+" does not exist.")
1302 ia.open(modim)
1303 modelosname='modelos_'+str(k)
1305 # clean up any temp files left from preveous incomplete run
1306 if os.path.exists(modelosname):
1307 ia.removefile(modelosname)
1308 if os.path.exists('__temp_model2'):
1309 ia.removefile('__temp_model2')
1311 modelos.append(modelosname)
1314 if( (ia.brightnessunit().count('/beam')) > 0):
1315 ##single dish-style model
1316 maskelos.append(modelos[k]+'.sdmask')
1317 self.im.makemodelfromsd(sdimage=modim,modelimage=modelos[k],maskimage=maskelos[k])
1318 ia.open(maskelos[k])
1319 ##sd mask cover whole image...delete it as it is not needed
1320 if((ia.statistics(verbose=False,list=False)['min']) >0):
1321 ia.remove(done=True, verbose=False)
1322 maskelos.remove(maskelos[k])
1323 ia.done()
1324 else:
1325 ##assuming its a model image already then just regrid it
1326 #self.im.make(modelos[k])
1327 shutil.copytree(self.imagelist[imindex],modelos[k])
1328 ia.open(modelos[k])
1329 newcsys=ia.coordsys()
1330 newshape=ia.shape()
1331 ia.open(modim)
1332 ib=ia.regrid(outfile=modelos[k], shape=newshape, axes=[0,1,3], csys=newcsys.torecord(), overwrite=True, asvelocity=False)
1333 ib.done(verbose=False)
1335 k=k+1
1336 if ia.isopen(): ia.close()
1337 #########
1338 if((len(maskelos)==1) and (self.outputmask == '')):
1339 self.outputmask=modelimages[0]+'.mask'
1340 if(os.path.exists(self.outputmask)):
1341 ia.removefile(self.outputmask)
1342 os.rename(maskelos[0],self.outputmask)
1343 elif(len(maskelos) > 0):
1344 if(self.outputmask == ''):
1345 self.outputmask=modelimages[0]+'.mask'
1347 else:
1348 outputmask=self.outputmask
1349 ##okay if outputmask exists then we need to do an "and" with
1350 ##the sdmask one
1351 doAnd=False;
1352 if(os.path.exists(outputmask)):
1353 ia.open(outputmask)
1354 if((ia.statistics(verbose=False, list=False)['max'].max()) > 0.00001):
1355 doAnd=True
1356 ia.close()
1357 if(doAnd):
1358 tmpomask='__temp_o_mask'
1359 self.makemaskimage(outputmask=tmpomask, maskobject=maskelos)
1360 os.rename(outputmask, '__temp_i_mask')
1361 outim = ia.imagecalc(outfile=outputmask, pixels='__temp_o_mask * __temp_i_mask', overwrite=True)
1362 outim.done(verbose=False)
1363 ia.removefile('__temp_o_mask')
1364 ia.removefile('__temp_i_mask')
1365 self.outputmask=outputmask
1366 else:
1367 self.makemaskimage(outputmask=outputmask, maskobject=maskelos)
1368 for ima in maskelos:
1369 if(os.path.exists(ima)):
1370 ia.removefile(ima)
1371 if(not (os.path.exists(outputmodel))):
1372 # im.make uses the main field coord. so it does
1373 # not make correct coord. for outlier fields
1374 if len(self.imagelist)>1:
1375 ia.fromimage(outputmodel, self.imagelist[imindex])
1376 else:
1377 self.im.make(outputmodel)
1378 #ia.open(outputmodel)
1379 #ia.close()
1380 for k in range(len(modelos)):
1381 # if os.rename() or shutil.move() is used here,
1382 # for k=1 and at image.imagecalc, it seems to cause
1383 # casapy to crash...
1384 #os.rename(outputmodel,'__temp_model2')
1385 shutil.copytree(outputmodel,'__temp_model2')
1387 outim = ia.imagecalc(outfile=outputmodel,
1388 pixels=modelos[k]+' + '+'__temp_model2',
1389 overwrite=True)
1390 outim.done(verbose=False)
1391 ia.removefile('__temp_model2')
1392 ia.removefile(modelos[k]);
1394 def readboxfile(self, boxfile):
1395 """ Read a file containing clean boxes (compliant with AIPS BOXFILE)
1397 Format is:
1398 #FIELDID BLC-X BLC-Y TRC-X TRC-Y
1399 0 110 110 150 150
1400 or
1401 0 hh:mm:ss.s dd.mm.ss.s hh:mm:ss.s dd.mm.ss.s
1403 Note all lines beginning with '#' are ignored.
1405 """
1406 union=[]
1407 polyg={}
1408 f=open(boxfile)
1409 temprec={}
1410 counter=0
1411 while 1:
1412 try:
1413 counter=counter+1
1414 line=f.readline()
1415 if(len(line)==0):
1416 raise Exception
1417 if (line.find('#')!=0):
1418 if(line.count('[')==2):
1419 ##its an output from qtclean
1420 line=line.replace('\n','')
1421 line=line.replace('\t',' ')
1422 line=line.replace('[',' ')
1423 line=line.replace(']',' ')
1424 line=line.replace(',',' ')
1425 splitline=line.split()
1426 if(len(splitline)==5):
1427 ##its box
1428 if(int(splitline[4]) > 0):
1429 ##it was a "mask" region not "erase"
1430 boxlist=[int(splitline[0]),int(splitline[1]),
1431 int(splitline[2]),int(splitline[3])]
1432 union.append(boxlist)
1433 else:
1434 #its a polygon
1435 x=[]
1436 y=[]
1437 if(int(splitline[len(splitline)-1]) > 0):
1438 ###ignore erase regions
1439 nnodes=(len(splitline)-1)/2
1440 for kk in range(nnodes):
1441 x.append(splitline[kk]+'pix')
1442 y.append(splitline[kk+nnodes]+'pix')
1443 elreg=rg.wpolygon(x=x, y=y, csys=self.csys)
1444 temprec.update({counter:elreg})
1446 elif(line.count('worldbox')==1):
1447 self._casalog.post('\'worldbox\' is deprecated please use CRTF format','WARN')
1448 #ascii box file from viewer or boxit
1449 # expected foramt: 'worldbox' pos_ref [lat
1450 line=line.replace('[',' ')
1451 line=line.replace(']',' ')
1452 line=line.replace(',',' ')
1453 line=line.replace('\'',' ')
1454 splitline=line.split()
1455 if len(splitline) != 13:
1456 raise TypeError('Error reading worldbox file')
1457 #
1458 refframe=self.csys['direction0']['conversionSystem']
1459 if refframe.find('_VLA')>0:
1460 refframe=refframe[0:refframe.find('_VLA')]
1461 ra =[splitline[2],splitline[3]]
1462 dec = [splitline[4],splitline[5]]
1463 #set frames
1464 obsdate=self.csys['obsdate']
1465 me.doframe(me.epoch(obsdate['refer'], str(obsdate['m0']['value'])+obsdate['m0']['unit']))
1466 me.doframe(me.observatory(self.csys['telescope']))
1467 #
1468 if splitline[1]!=refframe:
1469 # coversion between different epoch (and to/from AZEL also)
1470 radec0 = me.measure(me.direction(splitline[1],ra[0],dec[0]), refframe)
1471 radec1 = me.measure(me.direction(splitline[1],ra[1],dec[1]), refframe)
1472 ra=[str(radec0['m0']['value'])+radec0['m0']['unit'],\
1473 str(radec1['m0']['value'])+radec1['m0']['unit']]
1474 dec=[str(radec0['m1']['value'])+radec0['m1']['unit'],\
1475 str(radec1['m1']['value'])+radec1['m1']['unit']]
1476 # check for stokes
1477 stokes=[]
1478 imstokes = self.csys['stokes1']['stokes']
1479 for st in [splitline[10],splitline[11]]:
1480 prevlen = len(stokes)
1481 for i in range(len(imstokes)):
1482 if st==imstokes[i]:
1483 stokes.append(str(i)+'pix')
1484 elif st.count('pix') > 0:
1485 stokes.append(st)
1486 if len(stokes)<=prevlen:
1487 #raise TypeError, "Stokes %s for the box boundaries is outside image" % st
1488 self._casalog.post('Stokes %s for the box boundaries is outside image, -ignored' % st, 'WARN')
1489 # frequency
1490 freqs=[splitline[7].replace('s-1','Hz'), splitline[9].replace('s-1','Hz')]
1491 fframes=[splitline[6],splitline[8]]
1492 #imframe = self.csys['spectral2']['system']
1493 imframe = self.csys['spectral2']['conversion']['system']
1494 #casalog.post("imframe=" + imframe + " system frame =" + self.csys['spectral2']['system'] + " frame in boxfile=" + fframes[0])
1495 # the worldbox file created viewer's "box in file"
1496 # currently says TOPO in frequency axis but seems to
1497 # wrong (the freuencies look like in the image's base
1498 # frame).
1499 for k in [0,1]:
1500 if fframes[k]!=imframe and freqs[k].count('pix')==0:
1501 #do frame conversion
1502 #self._casalog.post('Ignoring the frequency frame of the box for now', 'WARN')
1503 # uncomment the following when box file correctly labeled the frequency frame
1504 me.doframe(me.direction(splitline[1],ra[k],dec[k]))
1505 mf=me.measure(me.frequency(fframes[k],freqs[k]),imframe)
1506 freqs[k]=str(mf['m0']['value'])+mf['m0']['unit']
1507 coordorder=self.csysorder
1508 wblc = []
1509 wtrc = []
1510 for type in coordorder:
1511 if type=='Direction':
1512 wblc.append(ra[0])
1513 wblc.append(dec[0])
1514 wtrc.append(ra[1])
1515 wtrc.append(dec[1])
1516 if type=='Stokes':
1517 wblc.append(stokes[0])
1518 wtrc.append(stokes[1])
1519 if type=='Spectral':
1520 wblc.append(freqs[0])
1521 wtrc.append(freqs[1])
1523 #wblc = [ra[0], dec[0], stokes[0], freqs[0]]
1524 #wtrc = [ra[1], dec[1], stokes[1], freqs[1]]
1525 #wblc = ra[0]+" "+dec[0]
1526 #wtrc = ra[1]+" "+dec[1]
1527 #casalog.post... "wblc=",wblc," wtrc=",wtrc," using system frame=",self.csys['spectral2']['system'], " convertion frame=",self.csys['spectral2']['conversion']['system']
1528 wboxreg = rg.wbox(blc=wblc,trc=wtrc,csys=self.csys)
1529 temprec.update({counter:wboxreg})
1531 else:
1532 ### its an AIPS boxfile
1533 splitline=line.split('\n')
1534 splitline2=splitline[0].split()
1535 if (len(splitline2)<6):
1536 if(int(splitline2[1])<0):
1537 ##circle
1538 #circlelist=[int(splitline2[2]),
1539 # int(splitline2[3]),int(splitline2[4])]
1540 #circles[splitline2[0]].append(circlelist)
1541 continue
1542 else:
1543 boxlist=[int(splitline2[1]),int(splitline2[2]),
1544 int(splitline2[3]),int(splitline2[4])]
1545 union.append(boxlist)
1546 else:
1547 ## Don't know what that is
1548 ## might be a facet definition
1549 continue
1551 except:
1552 break
1554 f.close()
1555 if(len(temprec)==1):
1556 polyg=temprec[temprec.keys()[0]]
1557 elif (len(temprec) > 1):
1558 #polyg=rg.dounion(temprec)
1559 polyg=rg.makeunion(temprec)
1561 return polyg,union
1563 def readmultifieldboxfile(self, boxfiles):
1564 """
1565 Read boxes and circles in text files in the
1566 AIPS clean boxfile format.
1568 Keyword arguments:
1569 boxfiles -- text files in boxfile format
1571 returns:
1572 circles -- dictionary containing circles
1573 boxes -- dictionary conatining boxes ([blc, trc])
1574 oldfileformats -- a list of boolean if the input textfiles
1575 are boxfile format.
1576 """
1577 circles={}
1578 boxes={}
1579 oldfilefmts={}
1580 for k in range(len(self.imageids)):
1581 circles[self.imageids[k]]=[]
1582 boxes[self.imageids[k]]=[]
1583 for boxfile in boxfiles:
1584 f=open(boxfile)
1585 setonce=False
1586 oldfilefmts[boxfile]=False
1587 while 1:
1588 try:
1589 line=f.readline()
1590 if(len(line)==0):
1591 raise Exception
1592 if line.find('#')==0:
1593 if not setonce and line.find('boxfile')>0:
1594 oldfilefmts[boxfile]=True
1595 setonce=True
1596 self._casalog.post(boxfile+" is in a deprecated boxfile format,"+\
1597 " will not be supported in the future releases","WARN")
1598 #raise Exception
1599 else:
1600 ### its an AIPS boxfile
1601 splitline=line.split('\n')
1602 splitline2=splitline[0].split()
1603 if (len(splitline2)<6):
1604 ##circles
1605 if(int(splitline2[1]) <0):
1606 circlelist=[int(splitline2[2]),
1607 int(splitline2[3]),int(splitline2[4])]
1608 #circles[splitline2[0]].append(circlelist)
1609 circles[self.imageids[int(splitline2[0])]].append(circlelist)
1610 else:
1611 #boxes
1612 if(len(splitline2)==5):
1613 boxlist=[int(splitline2[1]),int(splitline2[2]),
1614 int(splitline2[3]),int(splitline2[4])]
1615 #boxes[splitline2[0]].append(boxlist)
1616 boxes[self.imageids[int(splitline2[0])]].append(boxlist)
1617 else:
1618 ## Don't know what that is
1619 ## might be a facet definition
1620 continue
1624 except:
1625 break
1627 f.close()
1628 ###clean up the records
1629 for k in range(len(self.imageids)):
1630 if(circles[self.imageids[k]]==[]):
1631 circles.pop(self.imageids[k])
1632 if(boxes[self.imageids[k]]==[]):
1633 boxes.pop(self.imageids[k])
1635 return circles,boxes,oldfilefmts
1638 def readoutlier(self, outlierfile):
1639 """ Read a file containing clean boxes (kind of
1640 compliant with AIPS FACET FILE)
1642 Format is:
1643 col0 col1 col2 col3 col4 col5 col6 col7 col8 col9
1644 C FIELDID SIZEX SIZEY RAHH RAMM RASS DECDD DECMM DECSS
1645 why first column has to have C ... because its should
1646 not to be A or B ...now D would be a totally different thing.
1648 'C' as in AIPS BOXFILE format indicates the file specify the coordiates
1649 for field center(s).
1651 Note all lines beginning with '#' are ignored.
1652 (* Lines with first column other than C or c are also ignored)
1654 """
1655 imsizes=[]
1656 phasecenters=[]
1657 imageids=[]
1658 f=open(outlierfile)
1659 while 1:
1660 try:
1661 line=f.readline()
1663 if(len(line)==0):
1664 raise Exception
1665 if (line.find('#')!=0):
1666 splitline=line.split('\n')
1667 splitline2=splitline[0].split()
1668 if (len(splitline2)==10):
1669 if(splitline2[0].upper()=='C'):
1670 imageids.append(splitline2[1])
1671 imsizes.append((int(splitline2[2]),int(splitline2[3])))
1672 mydir='J2000 '+splitline2[4]+'h'+splitline2[5]+'m'+splitline2[6]+' '+splitline2[7]+'d'+splitline2[8]+'m'+splitline2[9]
1673 phasecenters.append(mydir)
1675 except:
1676 break
1678 f.close()
1679 return imsizes,phasecenters,imageids
1681 def newreadoutlier(self, outlierfile):
1682 """
1683 Read a outlier file (both old and new format)
1685 The new format consists of a set of task parameter inputs.
1686 imagename="outlier1" imsize=[128,128] phasecenter="J2000 hhmmss.s ddmmss.s"
1687 imagename="outlier2" imsize=[128,128] phasecenter="J2000 hhmmss.s ddmmss.s"
1688 mask=['viewermask.rgn','box[[50,60],[50,60]]'] ...
1690 Currently supported paramters are:
1691 imagename (required: used to delinate a paramater set for each field)
1692 imsize (required)
1693 phasecenter (required)
1694 mask (optional)
1695 modelimage (optional)
1696 * other parameters can be included in the file but not parsed
1698 For the old format, readoutlier() is called internally currently. But
1699 the support will be removed by CASA3.4.
1701 Returns:
1702 lists of imageids, imsizes, phasecenters, masks,
1703 and modelimages, and a dictionary contains all the parameters,
1704 and a boolean to indicate read file is in the new format.
1705 """
1706 import ast
1707 import re
1709 imsizes=[]
1710 phasecenters=[]
1711 imageids=[]
1712 masks=[]
1713 modelimages=[]
1714 keywd = "imagename"
1715 oldformat=False
1716 nchar= len(keywd)
1717 content0=''
1718 # set this to True to disallow old outlier file
1719 noOldOutlierFileSupport=False;
1721 with open(outlierfile) as f:
1722 for line in f:
1723 try:
1724 if len(line)!=0 and line.split()!=[]:
1725 if line.split()[0]=='C' or line.split()[0]=='c':
1726 oldformat=True
1727 elif line.split()[0]!='#':
1728 content0+=line
1729 except:
1730 casalog.post("Unkown error while reading the file: " + outlierfile)
1731 break
1732 if oldformat:
1733 if noOldOutlierFileSupport:
1734 self._casalog.post("You are using the old outlier file format no longer supported. Please use the new format (see help).","SEVERE")
1735 raise Exception
1736 else:
1737 self._casalog.post("This file format is deprecated. Use of a new format is encouraged.","WARN")
1738 # do old to new data format conversion....(watch out for different order of return parameters...)
1739 (imsizes,phasecenters,imageids)=self.readoutlier(outlierfile)
1740 for i in range(len(imageids)):
1741 modelimages.append('')
1742 #f.seek(0)
1743 #content0 = f.read()
1744 #f.close()
1745 content=content0.replace('\n',' ')
1746 last = len(content)
1747 # split the content using imagename as a key
1748 # and store each parameter set in pars dict
1749 pars={}
1750 initi = content.find(keywd)
1751 if initi > -1:
1752 i = 0
1753 prevstart=initi
1754 while True:
1755 #step = nchar*(i+1)
1756 step = nchar+1
1757 start = prevstart+step
1758 nexti = content[start:].find(keywd)
1759 #casalog.post... ("With start=",start, " found next one at(nexti)=",nexti, " step used =",step, " prevstart=",prevstart)
1761 if nexti == -1:
1762 pars[i]=content[prevstart:]
1763 # casalog.post("range=" + prevstart + " to the end")
1764 break
1765 pars[i]=content[prevstart:prevstart+nexti+step]
1766 #casalog.post("pars[",i,"]=" + pars[i])
1767 #casalog.post("range=" + prevstart + " to" + prevstart+nexti+step-1)
1768 prevstart=prevstart+nexti+step
1769 i+=1
1771 # for parsing of indiviual par (per field)
1772 #casalog.post("pars=" + pars)
1773 dparm ={}
1774 indx=0
1775 for key in pars.keys():
1776 # do parsing
1777 parstr = pars[key]
1778 # clean up extra white spaces
1779 parm=' '.join(parstr.split())
1780 # more clean up
1781 parm=parm.replace("[ ","[")
1782 parm=parm.replace(" ]","]")
1783 parm=parm.replace(" ,",",")
1784 parm=parm.replace(", ",",")
1785 parm=parm.replace(" =","=")
1786 parm=parm.replace("= ","=")
1787 #casalog.post("parm=" + parm)
1788 subdic={}
1789 # final parameter sets
1790 values=re.compile('\w+=').split(parm)
1791 values=values[1:len(values)]
1793 ipar = 0
1794 for pv in parm.split():
1795 if pv.find('=') != -1:
1796 if ipar >= len(values):
1797 raise TypeError("mismath in no. parameters in parsing outlier file.")
1798 (k,v) = pv.split('=')
1799 # fix a string to proper litral value
1800 # take out any commas at end which will
1801 # confuse literal_eval function
1802 pat = re.compile(',+$')
1803 subdic[k]=ast.literal_eval(pat.sub('',values[ipar]))
1804 ipar += 1
1805 dparm[indx]=subdic
1806 indx+=1
1807 #casalog.post("DONE for parsing parm for each field")
1808 # put into list of parameters
1809 # imsizes, phasecenters, imagenames(imageids)
1810 # mask is passed to other function for forther processing
1811 if not oldformat:
1812 #pack them by parameter name
1813 for fld in dparm.keys():
1814 # before process, check if it contains all required keys
1815 # namely, imagename, phasecenter, imsize
1816 #casalog.post("dparm[",fld,"]=" + dparm[fld])
1817 if not ("imagename" in dparm[fld] and \
1818 "phasecenter" in dparm[fld] and \
1819 "imsize" in dparm[fld]):
1820 raise TypeError("Missing one or more of the required parameters: \
1821 imagename, phasecenter, and imsize in the outlier file. Please check the input outlier file.")
1822 for key in dparm[fld].keys():
1823 if key == "imagename":
1824 imageids.append(dparm[fld][key])
1825 if key == "phasecenter":
1826 phasecenters.append(dparm[fld][key])
1827 if key == "imsize":
1828 imsizes.append(dparm[fld][key])
1829 if key == "mask":
1830 if type(dparm[fld][key])==str:
1831 masks.append([dparm[fld][key]])
1832 else:
1833 masks.append(dparm[fld][key])
1834 if key == "modelimage":
1835 if type(dparm[fld][key])==str:
1836 modelimages.append([dparm[fld][key]])
1837 else:
1838 modelimages.append(dparm[fld][key])
1839 if "mask" not in dparm[fld]:
1840 masks.append([])
1841 if "modelimage" not in dparm[fld]:
1842 modelimages.append('')
1845 return (imageids,imsizes,phasecenters,masks,modelimages,dparm, not oldformat)
1848 def copymaskimage(self, maskimage, shp, outfile):
1849 """
1850 Copy mask image
1852 Keyword arguments:
1853 maskimage -- input maskimage
1854 shp -- shape of output image
1855 outfile -- output image name
1856 """
1857 if outfile == maskimage: # Make it a no-op,
1858 return # this is more than just peace of mind.
1859 #pdb.set_trace()
1860 ia.open(maskimage)
1861 oldshp=ia.shape()
1862 if((len(oldshp) < 4) or (shp[2] != oldshp[2]) or (shp[3] != oldshp[3])):
1863 #take the first plane of mask
1864 tmpshp=oldshp
1865 tmpshp[0]=shp[0]
1866 tmpshp[1]=shp[1]
1867 if len(oldshp)==4: # include spectral axis for regrid
1868 tmpshp[3]=shp[3]
1869 ib=ia.regrid(outfile='__looloo', shape=tmpshp, axes=[3,0,1], csys=self.csys, overwrite=True, asvelocity=False)
1870 else:
1871 ib=ia.regrid(outfile='__looloo', shape=tmpshp, axes=[0,1], csys=self.csys, overwrite=True, asvelocity=False)
1873 #dat=ib.getchunk()
1874 ib.done(verbose=False)
1875 ia.fromshape(outfile=outfile, shape=shp, csys=self.csys, overwrite=True)
1876 ##getchunk is a massive memory hog
1877 ###so going round in a funny fashion
1878 #arr=ia.getchunk()
1879 #for k in range(shp[2]):
1880 # for j in range(shp[3]):
1881 # if(len(dat.shape)==2):
1882 # arr[:,:,k,j]=dat
1883 # elif(len(dat.shape)==3):
1884 # arr[:,:,k,j]=dat[:,:,0]
1885 # else:
1886 # arr[:,:,k,j]=dat[:,:,0,0]
1887 #ia.putchunk(arr)
1888 ia.calc(outfile+'[index3 in [0]]+__looloo')
1889 ia.calcmask('mask(__looloo)')
1890 ia.done(verbose=False)
1891 ia.removefile('__looloo')
1892 else:
1893 ib=ia.regrid(outfile=outfile ,shape=shp, axes=[0,1], csys=self.csys, overwrite=True, asvelocity=False)
1894 ia.done(verbose=False)
1895 ib.done(verbose=False)
1898 def flatten(self,l):
1899 """
1900 A utility function to flatten nested lists
1901 but allow nesting of [[elm1,elm2,elm3],[elm4,elm5],[elm6,elm7]]
1902 to handle multifield masks.
1903 This does not flatten if an element is a list of int or float.
1904 And also leave empty list as is.
1905 """
1906 retlist = []
1907 l = list(l)
1908 #casalog.post('l='+l)
1909 for i in range(len(l)):
1910 #casalog.post("ith l="+i+ l[i] )
1911 if isinstance(l[i],list) and l[i]:
1912 # and (not isinstance(l[i][0],(int,float))):
1913 #casalog.post("recursive l="+l)
1914 if isinstance(l[i][0],list) and isinstance(l[i][0][0],list):
1915 retlist.extend(self.flatten(l[i]))
1916 else:
1917 retlist.append(l[i])
1918 else:
1919 retlist.append(l[i])
1920 return retlist
1922 def getchanimage(self,cubeimage,outim,chan):
1923 """
1924 Create a slice of channel image from cubeimage
1926 Keyword arguments:
1927 cubeimage -- input image cube
1928 outim -- output sliced image
1929 chan -- nth channel
1930 """
1931 #pdb.set_trace()
1932 ia.open(cubeimage)
1933 modshape=ia.shape()
1934 if modshape[3]==1:
1935 return False
1936 if modshape[3]-1 < chan:
1937 return False
1938 blc=[0,0,modshape[2]-1,chan]
1939 trc=[modshape[0]-1,modshape[1]-1,modshape[2]-1,chan]
1940 sbim=ia.subimage(outfile=outim, region=rg.box(blc,trc), overwrite=True)
1941 sbim.close()
1942 ia.close()
1943 return True
1945 def putchanimage(self,cubimage,inim,chan):
1946 """
1947 Put channel image back to a pre-exisiting cubeimage
1949 Keyword arguments:
1950 cubimage -- image cube
1951 inim -- input channel image
1952 chan -- nth channel
1953 """
1954 ia.open(inim)
1955 inimshape=ia.shape()
1956 imdata=ia.getchunk()
1957 immask=ia.getchunk(getmask=True)
1958 ia.close()
1959 blc=[0,0,inimshape[2]-1,chan]
1960 trc=[inimshape[0]-1,inimshape[1]-1,inimshape[2]-1,chan]
1961 ia.open(cubimage)
1962 cubeshape=ia.shape()
1963 if not (cubeshape[3] > (chan+inimshape[3]-1)):
1964 return False
1965 #rg0=ia.setboxregion(blc=blc,trc=trc)
1966 rg0=rg.box(blc=blc,trc=trc)
1967 if inimshape[0:3]!=cubeshape[0:3]:
1968 return False
1969 #ia.putchunk(pixels=imdata,blc=blc)
1970 ia.putregion(pixels=imdata,pixelmask=immask, region=rg0)
1971 ia.close()
1972 return True
1974 def qatostring(self,q):
1975 """
1976 A utility function to return a quantity in string
1977 (currently only used in setChannelization which is deprecated)
1978 """
1979 if 'unit' not in q:
1980 raise TypeError("Does not seems to be quantity")
1981 return str(q['value'])+q['unit']
1983 def convertvf(self,vf,frame,field,restf,veltype='radio'):
1984 """
1985 returns doppler(velocity) or frequency in string
1986 currently use first rest frequency
1987 Assume input vf (velocity or fequency in a string) and
1988 output are the same 'frame'.
1989 """
1990 #pdb.set_trace()
1991 docalcf=False
1992 #if(frame==''): frame='LSRK'
1993 #Use datasepcframe, it is cleanhelper initialized to set
1994 #to LSRK
1995 if(frame==''): frame=self.dataspecframe
1996 if(qa.quantity(vf)['unit'].find('m/s') > -1):
1997 docalcf=True
1998 elif(qa.quantity(vf)['unit'].find('Hz') > -1):
1999 docalcf=False
2000 else:
2001 if vf !=0:
2002 raise TypeError("Unrecognized unit for the velocity or frequency parameter")
2003 ##fldinds=ms.msseltoindex(self.vis, field=field)['field'].tolist()
2004 fldinds=ms.msseltoindex(self.vis[self.sortedvisindx[0]], field=field)['field'].tolist()
2005 if(len(fldinds) == 0):
2006 fldid0=0
2007 else:
2008 fldid0=fldinds[0]
2009 if restf=='':
2010 #tb.open(self.vis+'/FIELD')
2011 fldtab=self.getsubtable(self.vis[self.sortedvisindx[0]],'FIELD')
2012 tb.open(fldtab)
2013 nfld = tb.nrows()
2014 if nfld >= fldid0:
2015 srcid=tb.getcell('SOURCE_ID',fldid0)
2016 else:
2017 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
2018 "no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
2019 tb.close()
2020 # SOUECE_ID in FIELD table = -1 if no SOURCE table
2021 if srcid==-1:
2022 raise TypeError("Rest frequency info is not supplied")
2023 #tb.open(self.vis+'/SOURCE')
2024 sourcetab=self.getsubtable(self.vis[self.sortedvisindx[0]], 'SOURCE')
2025 tb.open(sourcetab)
2026 tb2=tb.query('SOURCE_ID==%s' % srcid)
2027 tb.close()
2028 nsrc = tb2.nrows()
2029 if nsrc > 0:
2030 rfreq=tb2.getcell('REST_FREQUENCY',0)
2031 else:
2032 raise TypeError( "Cannot set REST_FREQUENCY from the data: "+
2033 " no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
2034 tb2.close()
2035 if(rfreq<=0):
2036 raise TypeError("Rest frequency does not seems to be properly set, check the data")
2037 else:
2038 if type(restf)==str: restf=[restf]
2039 if(qa.quantity(restf[0])['unit'].find('Hz') > -1):
2040 rfreq=[qa.convert(qa.quantity(restf[0]),'Hz')['value']]
2041 #casalog.post("using user input rest freq=" + rfreq)
2042 else:
2043 raise TypeError("Unrecognized unit or type for restfreq")
2044 if(vf==0):
2045 # assume just want to get a restfrequecy from the data
2046 ret=str(rfreq[0])+'Hz'
2047 else:
2048 if(docalcf):
2049 dop=me.doppler(veltype, qa.quantity(vf))
2050 rvf=me.tofrequency(frame, dop, qa.quantity(rfreq[0],'Hz'))
2051 else:
2052 frq=me.frequency(frame, qa.quantity(vf))
2053 rvf=me.todoppler(veltype, frq, qa.quantity(rfreq[0],'Hz'))
2054 ret=str(rvf['m0']['value'])+rvf['m0']['unit']
2055 return ret
2058 def getfreqs(self,nchan,spw,start,width, dummy=False):
2059 """
2060 (not in used - currently commented out in its caller, initChaniter())
2061 returns a list of frequencies to be used in output clean image
2062 if width = -1, start is actually end (max) freq
2063 """
2064 #pdb.set_trace()
2065 freqlist=[]
2066 finc=1
2067 loc_nchan=0
2069 if spw in (-1, '-1', '*', '', ' '):
2070 spwinds = -1
2071 else:
2072 #spwinds=ms.msseltoindex(self.vis, spw=spw)['spw'].tolist()
2073 spwinds=ms.msseltoindex(self.vis[self.sortedvisindx[0]], spw=spw)['spw'].tolist()
2074 if(len(spwinds) == 0):
2075 spwinds = -1
2077 if(spwinds==-1):
2078 # first row
2079 spw0=0
2080 else:
2081 spw0=spwinds[0]
2082 #tb.open(self.vis+'/SPECTRAL_WINDOW')
2083 spectable=self.getsubtable(self.vis[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
2084 tb.open(spectable)
2085 chanfreqscol=tb.getvarcol('CHAN_FREQ')
2086 chanwidcol=tb.getvarcol('CHAN_WIDTH')
2087 spwframe=tb.getcol('MEAS_FREQ_REF');
2088 tb.close()
2089 # assume spw[0]
2090 elspecframe=["REST",
2091 "LSRK",
2092 "LSRD",
2093 "BARY",
2094 "GEO",
2095 "TOPO",
2096 "GALACTO",
2097 "LGROUP",
2098 "CMB"]
2099 self.dataspecframe=elspecframe[spwframe[spw0]];
2100 if(dummy):
2101 return freqlist, finc
2102 #DP extract array from dictionary returned by getvarcol
2103 chanfreqs1dx = numpy.array([])
2104 chanfreqs=chanfreqscol['r'+str(spw0+1)].transpose()
2105 chanfreqs1dx = chanfreqs[0]
2106 if(spwinds!=-1):
2107 for ispw in range(1,len(spwinds)):
2108 chanfreqs=chanfreqscol['r'+str(spwinds[ispw]+1)].transpose()
2109 chanfreqs1dx = numpy.concatenate((chanfreqs1dx, chanfreqs[0]))
2110 chanfreqs1d = chanfreqs1dx.flatten()
2111 #RI this is woefully inadequate assuming the first chan's width
2112 #applies to everything selected, but we're going to replace all
2113 #this with MSSelect..
2114 chanwids=chanwidcol['r'+str(spw0+1)].transpose()
2115 chanfreqwidth=chanwids[0][0]
2117 if(type(start)==int or type(start)==float):
2118 if(start > len(chanfreqs1d)):
2119 raise TypeError("Start channel is outside the data range")
2120 startf = chanfreqs1d[start]
2121 elif(type(start)==str):
2122 if(qa.quantity(start)['unit'].find('Hz') > -1):
2123 startf=qa.convert(qa.quantity(start),'Hz')['value']
2124 else:
2125 raise TypeError("Unrecognized start parameter")
2126 if(type(width)==int or type(width)==float):
2127 if(type(start)==int or type(start)==float):
2128 #finc=(chanfreqs1d[start+1]-chanfreqs1d[start])*width
2129 finc=(chanfreqwidth)*width
2130 # still need to convert to target reference frame!
2131 elif(type(start)==str):
2132 if(qa.quantity(start)['unit'].find('Hz') > -1):
2133 # assume called from setChannelization with local width=1
2134 # for the default width(of clean task parameter)='' for
2135 # velocity and frequency modes. This case set width to
2136 # first channel width (for freq) and last one (for vel)
2137 if width==-1:
2138 finc=chanfreqs1d[-1]-chanfreqs1d[-2]
2139 else:
2140 finc=chanfreqs1d[1]-chanfreqs1d[0]
2142 # still need to convert to target reference frame!
2143 elif(type(width)==str):
2144 if(qa.quantity(width)['unit'].find('Hz') > -1):
2145 finc=qa.convert(qa.quantity(width),'Hz')['value']
2146 if(nchan ==-1):
2147 if(qa.quantity(start)['unit'].find('Hz') > -1):
2148 if width==-1: # must be in velocity order (i.e. startf is max)
2149 bw=startf-chanfreqs1d[0]
2150 else:
2151 bw=chanfreqs1d[-1]-startf
2152 else:
2153 bw=chanfreqs1d[-1]-chanfreqs1d[start]
2154 if(bw < 0):
2155 raise TypeError("Start parameter is outside the data range")
2156 if(qa.quantity(width)['unit'].find('Hz') > -1):
2157 qnchan=qa.convert(qa.div(qa.quantity(bw,'Hz'),qa.quantity(width)))
2158 #DP loc_nchan=int(math.ceil(qnchan['value']))+1
2159 loc_nchan=int(round(qnchan['value']))+1
2160 else:
2161 #DP loc_nchan=int(math.ceil(bw/finc))+1
2162 loc_nchan=int(round(bw/finc))+1
2163 else:
2164 loc_nchan=nchan
2165 for i in range(int(loc_nchan)):
2166 if(i==0):
2167 freqlist.append(startf)
2168 else:
2169 freqlist.append(freqlist[-1]+finc)
2170 return freqlist, finc
2173 def setChannelizeDefault(self,mode,spw,field,nchan,start,width,frame,veltype,phasec, restf,obstime=''):
2174 """
2175 Determine appropriate values for channelization
2176 parameters when default values are used
2177 for mode='velocity' or 'frequency' or 'channel'
2178 This makes use of ms.cvelfreqs.
2179 """
2180 ###############
2181 # for debugging
2182 ###############
2183 debug=False
2184 ###############
2185 spectable=self.getsubtable(self.vis[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
2186 tb.open(spectable)
2187 chanfreqscol=tb.getvarcol('CHAN_FREQ')
2188 chanwidcol=tb.getvarcol('CHAN_WIDTH')
2189 spwframe=tb.getcol('MEAS_FREQ_REF');
2190 tb.close()
2191 # first parse spw parameter:
2192 # use MSSelect if possible
2193 if len(self.sortedvislist) > 0:
2194 invis = self.sortedvislist[0]
2195 inspw = self.vis.index(self.sortedvislist[0])
2196 else:
2197 invis = self.vis[0]
2198 inspw = 0
2199 ms.open(invis)
2200 if type(spw)==list:
2201 spw=spw[inspw]
2202 if spw in ('-1', '*', '', ' '):
2203 spw='*'
2204 if field=='':
2205 field='*'
2206 mssel=ms.msseltoindex(vis=invis, spw=spw, field=field)
2207 selspw=mssel['spw']
2208 selfield=mssel['field']
2209 chaninds=mssel['channel'].tolist()
2210 chanst0 = chaninds[0][1]
2212 # frame
2213 spw0=selspw[0]
2214 chanfreqs=chanfreqscol['r'+str(spw0+1)].transpose()[0]
2215 chanres = chanwidcol['r'+str(spw0+1)].transpose()[0]
2217 # ascending or desending data frequencies?
2218 # based on selected first spw's first CHANNEL WIDTH
2219 # ==> some MS data may have positive chan width
2220 # so changed to look at first two channels of chanfreq (TT)
2221 #if chanres[0] < 0:
2222 descending = False
2223 if len(chanfreqs) > 1 :
2224 if chanfreqs[1]-chanfreqs[0] < 0:
2225 descending = True
2226 else:
2227 if chanres[0] < 0:
2228 descending = True
2230 # set dataspecframe:
2231 elspecframe=["REST",
2232 "LSRK",
2233 "LSRD",
2234 "BARY",
2235 "GEO",
2236 "TOPO",
2237 "GALACTO",
2238 "LGROUP",
2239 "CMB"]
2240 self.dataspecframe=elspecframe[spwframe[spw0]];
2242 # set usespecframe: user's frame if set, otherwise data's frame
2243 if(frame != ''):
2244 self.usespecframe=frame
2245 self.inframe=True
2246 else:
2247 self.usespecframe=self.dataspecframe
2249 # some start and width default handling
2250 if mode!='channel':
2251 if width==1:
2252 width=''
2253 if start==0:
2254 start=''
2256 #get restfreq
2257 if restf=='':
2258 fldtab=self.getsubtable(invis,'FIELD')
2259 tb.open(fldtab)
2260 nfld=tb.nrows()
2261 try:
2262 if nfld >= selfield[0]:
2263 srcid=tb.getcell('SOURCE_ID',selfield[0])
2264 else:
2265 if mode=='velocity':
2266 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
2267 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
2268 finally:
2269 tb.close()
2270 #SOUECE_ID in FIELD table = -1 if no SOURCE table
2271 if srcid==-1:
2272 if mode=='velocity':
2273 raise TypeError("Rest frequency info is not supplied")
2274 try:
2275 srctab=self.getsubtable(invis, 'SOURCE')
2276 tb.open(srctab)
2277 tb2=tb.query('SOURCE_ID==%s' % srcid)
2278 nsrc = tb2.nrows()
2279 if nsrc > 0 and tb2.iscelldefined('REST_FREQUENCY',0):
2280 rfqs = tb2.getcell('REST_FREQUENCY',0)
2281 if len(rfqs)>0:
2282 restf=str(rfqs[0])+'Hz'
2283 else:
2284 if mode=='velocity':
2285 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
2286 "REST_FREQUENCY entry for ID %s in SOURCE table is empty, please supply restfreq" % srcid )
2287 else:
2288 if mode=='velocity':
2289 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
2290 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
2291 finally:
2292 tb.close()
2293 tb2.close()
2295 if type(phasec)==list:
2296 inphasec=phasec[0]
2297 else:
2298 inphasec=phasec
2299 if type(inphasec)==str and inphasec.isdigit():
2300 inphasec=int(inphasec)
2301 #if nchan==1:
2302 # use data chan freqs
2303 # newfreqs=chanfreqs
2304 #else:
2305 # obstime not included here
2306 if debug: casalog.post("before ms.cvelfreqs (start,width,nchan)===> {} {} {}".
2307 format(start, width, nchan))
2308 try:
2309 newfreqs=ms.cvelfreqs(spwids=selspw,fieldids=selfield,mode=mode,nchan=nchan,
2310 start=start,width=width,phasec=inphasec, restfreq=restf,
2311 outframe=self.usespecframe,veltype=veltype).tolist()
2312 except:
2313 # ms must be closed here if ms.cvelfreqs failed with an exception
2314 ms.close()
2315 raise
2316 ms.close()
2318 #casalog.post(newfreqs)
2319 descendingnewfreqs=False
2320 if len(newfreqs)>1:
2321 if newfreqs[1]-newfreqs[0] < 0:
2322 descendingnewfreqs=True
2325 try:
2326 if((nchan in [-1, "-1", "", " "]) or (start in [-1, "-1", "", " "])):
2327 frange=im.advisechansel(msname=invis, spwselection=spw, fieldid=selfield[0], getfreqrange=True)
2328 startchan=0
2329 endchan=len(newfreqs)-1
2330 if(descendingnewfreqs):
2331 startchan=numpy.min(numpy.where(frange['freqend'] < numpy.array(newfreqs)))
2332 endchan=numpy.min(numpy.where(frange['freqstart'] < numpy.array(newfreqs)))
2333 else:
2334 startchan=numpy.max(numpy.where(frange['freqstart'] > numpy.array(newfreqs)))
2335 endchan=numpy.max(numpy.where(frange['freqend'] > numpy.array(newfreqs)))
2336 if((start not in [-1, "-1", "", " "]) and (mode=="channel")):
2337 startchan=start
2338 if(nchan not in [-1, "-1", "", " "]):
2339 endchan=startchan+nchan-1
2340 newfreqs=(numpy.array(newfreqs)[startchan:endchan]).tolist()
2341 except:
2342 pass
2343 if debug: casalog.post("Mode, Start, width after cvelfreqs = {} {} {}".
2344 format(mode, start, width))
2345 if type(newfreqs)==list and len(newfreqs) ==0:
2346 raise TypeError( "Output frequency grid cannot be calculated: " +
2347 " please check start and width parameters" )
2348 if debug:
2349 if len(newfreqs)>1:
2350 casalog.post("FRAME=" + self.usespecframe)
2351 casalog.post("newfreqs[0]===>" + newfreqs[0])
2352 casalog.post("newfreqs[1]===>" + newfreqs[1])
2353 casalog.post("newfreqs[-1]===>" + newfreqs[-1])
2354 casalog.post("len(newfreqs)===>" + len(newfreqs))
2355 else:
2356 casalog.post("newfreqs=" + newfreqs)
2358 # set output number of channels
2359 if nchan ==1:
2360 retnchan=1
2361 else:
2362 if len(newfreqs)>1:
2363 retnchan=len(newfreqs)
2364 else:
2365 retnchan=nchan
2366 newfreqs=chanfreqs.tolist()
2368 # set start parameter
2369 # first analyze data order etc
2370 reverse=False
2371 negativew=False
2372 if descending:
2373 # channel mode case (width always >0)
2374 if width!="" and (type(width)==int or type(width)==float):
2375 if descendingnewfreqs:
2376 reverse=False
2377 else:
2378 reverse=True
2379 elif width=="": #default width
2380 if descendingnewfreqs and mode=="frequency":
2381 reverse=False
2382 else:
2383 reverse=True
2385 elif type(width)==str:
2386 if width.lstrip().find('-')==0:
2387 negativew=True
2388 if descendingnewfreqs:
2389 if negativew:
2390 reverse=False
2391 else:
2392 reverse=True
2393 else:
2394 if negativew:
2395 reverse=True
2396 else:
2397 reverse=False
2398 else: #ascending data
2399 # depends on sign of width only
2400 # with CAS-3117 latest change(rev.15179), velocity start
2401 # means lowest velocity for default width
2402 if width=="" and mode=="velocity": #default width
2403 # ms.cvelfreqs returns correct order so no reversing
2404 reverse=False
2405 elif type(width)==str:
2406 if width.lstrip().find('-')==0:
2407 reverse=True
2408 else:
2409 reverse=False
2411 if reverse:
2412 newfreqs.reverse()
2413 #if (start!="" and mode=='channel') or \
2414 # (start!="" and type(start)!=int and mode!='channel'):
2415 # for now to avoid inconsistency later in imagecoordinates2 call
2416 # user's start parameter is preserved for channel mode only.
2417 # (i.e. the current code may adjust start parameter for other modes but
2418 # this probably needs to be changed, especially for multiple ms handling.)
2419 if ((start not in [-1, "", " "]) and mode=='channel'):
2420 retstart=start
2421 else:
2422 # default cases
2423 if mode=="frequency":
2424 retstart=str(newfreqs[0])+'Hz'
2425 elif mode=="velocity":
2426 #startfreq=str(newfreqs[-1])+'Hz'
2427 startfreq=(str(max(newfreqs))+'Hz') if(start=="") else (str(newfreqs[-1])+'Hz')
2428 retstart=self.convertvf(startfreq,frame,field,restf,veltype)
2429 elif mode=="channel":
2430 # default start case, use channel selection from spw
2431 retstart=chanst0
2433 # set width parameter
2434 if width!="":
2435 retwidth=width
2436 else:
2437 if nchan==1:
2438 finc = chanres[0]
2439 else:
2440 finc = newfreqs[1]-newfreqs[0]
2441 if debug: casalog.post("finc(newfreqs1-newfreqs0)=" + finc)
2442 if mode=="frequency":
2443 # It seems that this is no longer necessary... TT 2013-08-12
2444 #if descendingnewfreqs:
2445 # finc = -finc
2446 retwidth=str(finc)+'Hz'
2447 elif mode=="velocity":
2448 # for default width assume it is vel<0 (incresing in freq)
2449 if descendingnewfreqs:
2450 ind1=-2
2451 ind0=-1
2452 else:
2453 ind1=-1
2454 ind0=-2
2455 v1 = self.convertvf(str(newfreqs[ind1])+'Hz',frame,field,restf,veltype=veltype)
2456 v0 = self.convertvf(str(newfreqs[ind0])+'Hz',frame,field,restf,veltype=veltype)
2457 ##v1 = self.convertvf(str(newfreqs[-1])+'Hz',frame,field,restf,veltype=veltype)
2458 ##v0 = self.convertvf(str(newfreqs[-2])+'Hz',frame,field,restf,veltype=veltype)
2459 #v1 = self.convertvf(str(newfreqs[1])+'Hz',frame,field,restf,veltype=veltype)
2460 #v0 = self.convertvf(str(newfreqs[0])+'Hz',frame,field,restf,veltype=veltype)
2461 if(qa.lt(v0, v1) and start==""):
2462 ###user used "" as start make sure step is +ve in vel as start is min vel possible for freqs selected
2463 retwidth=qa.tos(qa.sub(v1, v0))
2464 else:
2465 retwidth = qa.tos(qa.sub(v0, v1))
2466 else:
2467 retwidth=1
2468 if debug: casalog.post("setChan retwidth=" + retwidth)
2469 return retnchan, retstart, retwidth
2471 def setChannelizeNonDefault(self, mode,spw,field,nchan,start,width,frame,
2472 veltype,phasec, restf):
2473 """
2474 Determine appropriate values for channelization
2475 parameters when default values are used
2476 for mode='velocity' or 'frequency' or 'channel'
2477 This does not replaces setChannelization and make no use of ms.cvelfreqs.
2478 """
2481 #spw='0:1~4^2;10~12, ,1~3:3~10^3,4~6,*:7'
2482 #vis='ngc5921/ngc5921.demo.ms'
2484 if type(spw)!=str:
2485 spw=''
2487 if spw.strip()=='':
2488 spw='*'
2490 freqs=set()
2491 wset=[]
2492 chunk=spw.split(',')
2493 for i in range(len(chunk)):
2494 #casalog.post(chunk[i] + '------')
2495 ck=chunk[i].strip()
2496 if len(ck)==0:
2497 continue
2499 wc=ck.split(':')
2500 window=wc[0].strip()
2502 if len(wc)==2:
2503 sec=wc[1].split(';')
2504 for k in range(len(sec)):
2505 chans=sec[k].strip()
2506 sep=chans.split('^')
2507 se=sep[0].strip()
2508 t=1
2509 if len(sep)==2:
2510 t=sep[1].strip()
2511 se=se.split('~')
2512 s=se[0].strip()
2513 if len(se)==2:
2514 e=se[1].strip()
2515 else:
2516 e=-1
2517 wd=window.split('~')
2518 if len(wd)==2:
2519 wds=int(wd[0])
2520 wde=int(wd[1])
2521 for l in range(wds, wde):
2522 #casalog.post... (l, s, e, t)
2523 wset.append([l, s, e, t])
2524 else:
2525 #casalog.post ...(wd[0], s, e, t)
2526 if e==-1:
2527 try:
2528 e=int(s)+1
2529 except:
2530 e=s
2531 wset.append([wd[0], s, e, t])
2532 else:
2533 win=window.split('~')
2534 if len(win)==2:
2535 wds=int(win[0])
2536 wde=int(win[1])
2537 for l in range(wds, wde):
2538 #casalog.post... (l, 0, -1, 1)
2539 wset.append([l, 0, -1, 1])
2540 else:
2541 #casalog.post...(win[0], 0, -1, 1)
2542 wset.append([win[0], 0, -1, 1])
2544 #casalog.post(wset)
2545 for i in range(len(wset)):
2546 for j in range(4):
2547 try:
2548 wset[i][j]=int(wset[i][j])
2549 except:
2550 wset[i][j]=-1
2551 #casalog.post(wset)
2552 spectable=self.getsubtable(self.vis[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
2553 tb.open(spectable)
2554 nr=tb.nrows()
2555 for i in range(len(wset)):
2556 if wset[i][0]==-1:
2557 w=range(nr)
2558 elif wset[i][0]<nr:
2559 w=[wset[i][0]]
2560 else:
2561 w=range(0)
2562 for j in w:
2563 chanfreqs=tb.getcell('CHAN_FREQ', j)
2564 if wset[i][2]==-1:
2565 wset[i][2]=len(chanfreqs)
2566 if wset[i][2]>len(chanfreqs):
2567 wset[i][2]=len(chanfreqs)
2568 #casalog.post...(wset[i][1], wset[i][2], len(chanfreqs), wset[i][3])
2569 for k in range(wset[i][1], wset[i][2], wset[i][3]):
2570 #casalog.post(k)
2571 freqs.add(chanfreqs[k])
2572 tb.close()
2573 freqs=list(freqs)
2574 freqs.sort()
2575 #casalog.post...(freqs[0], freqs[-1])
2577 if mode=='channel':
2578 star=0
2579 if type(start)==str:
2580 try:
2581 star=int(start)
2582 except:
2583 star=0
2584 if type(start)==int:
2585 star=start
2586 if star>len(freqs) or star<0:
2587 star=0
2589 if nchan==-1:
2590 nchan=len(freqs)
2592 widt=1
2593 if type(width)==str:
2594 try:
2595 widt=int(width)
2596 except:
2597 widt=1
2598 if type(width)==int:
2599 widt=width
2600 if widt==0:
2601 widt=1
2602 if widt>0:
2603 nchan=max(min(int((len(freqs)-star)/widt), nchan), 1)
2604 else:
2605 nchan=max(min(int((-star)/widt), nchan), 1)
2606 widt=-widt
2607 star=max(star-nchan*widt, 0)
2609 if mode=='frequency':
2610 star=freqs[0]
2611 if type(start)!=str:
2612 star=freqs[0]
2613 else:
2614 star=max(qa.quantity(start)['value'], freqs[0])
2616 if nchan==-1:
2617 nchan=len(freqs)
2619 widt=freqs[-1]
2620 if len(freqs)>1:
2621 for k in range(len(freqs)-1):
2622 widt=min(widt, freqs[k+1]-freqs[k])
2623 if type(width)==str and width.strip()!='':
2624 widt=qa.quantity(width)['value']
2626 if widt>0:
2627 #casalog.post...(star, widt, (freqs[-1]-star)//widt)
2628 nchan=max(min(int((freqs[-1]-star)//widt), nchan), 1)
2629 else:
2630 nchan=max(min(int(freqs[0]-star)//widt, nchan), 1)
2631 widt=-widt
2632 star=max(star-nchan*widt, freqs[0])
2634 widt=str(widt)+'Hz'
2635 star=str(star)+'Hz'
2637 if mode=='velocity':
2638 beg1=self.convertvf(str(freqs[0])+'Hz',frame,field,restf,veltype=veltype)
2639 beg1=qa.quantity(beg1)['value']
2640 end0=self.convertvf(str(freqs[-1])+'Hz',frame,field,restf,veltype=veltype)
2641 end0=qa.quantity(end0)['value']
2642 star=beg1
2643 if type(start)==str and start.strip()!='':
2644 star=min(qa.quantity(start)['value'], star)
2645 star=min(star, end0)
2647 #casalog.post...(beg1, star, end0)
2649 widt=-end0+beg1
2650 if len(freqs)>1:
2651 for k in range(len(freqs)-1):
2652 st=self.convertvf(str(freqs[k])+'Hz',frame,field,restf,veltype=veltype)
2653 en=self.convertvf(str(freqs[k+1])+'Hz',frame,field,restf,veltype=veltype)
2654 widt=min(widt, qa.quantity(en)['value']-qa.quantity(st)['value'])
2655 widt=-abs(widt)
2657 if type(width)==str and width.strip()!='':
2658 widt=qa.quantity(width)['value']
2660 #casalog.post(widt)
2661 if widt>0:
2662 nchan=max(min(int((beg1-star)/widt), nchan), 1)
2663 #star=0
2664 else:
2665 nchan=max(min(int((end0-star)/widt), nchan), 1)
2666 #widt=-widt
2668 widt=str(widt)+'m/s'
2669 star=str(star)+'m/s'
2671 return nchan, star, widt
2673 def convertframe(self,fin,frame,field):
2674 """
2675 (not in use: its caller is setChannelization...)
2676 convert freq frame in dataframe to specfied frame, assume fin in Hz
2677 retruns converted freq in Hz (value only)
2678 """
2679 # assume set to phasecenter before initChanelization is called
2680 pc=self.srcdir
2681 if(type(pc)==str):
2682 if (pc==''):
2683 fieldused = field
2684 if (fieldused ==''):
2685 fieldused ='0'
2686 #dir = int(ms.msseltoindex(self.vis,field=fieldused)['field'][0])
2687 dir = int(ms.msseltoindex(self.vis[self.sortedvisindx[0]],field=fieldused)['field'][0])
2688 else:
2689 tmpdir = phasecenter
2690 try:
2691 #if(len(ms.msseltoindex(self.vis, field=pc)['field']) > 0):
2692 if(len(ms.msseltoindex(self.vis[self.sortedvisindx[0]], field=pc)['field']) > 0):
2693 #tmpdir = int(ms.msseltoindex(self.vis,field=pc)['field'][0])
2694 tmpdir = int(ms.msseltoindex(self.vis[self.sortedvisindx[0]],field=pc)['field'][0])
2695 except Exception as instance:
2696 tmpdir = pc
2697 dir = tmpdir
2698 if type(dir)==str:
2699 try:
2700 mrf, ra, dec = dir.split()
2701 except Exception as instance:
2702 raise TypeError("Error in a string format for phasecenter")
2703 mdir = me.direction(mrf, ra, dec)
2704 else:
2705 #tb.open(self.vis+'/FIELD')
2706 fldtab=self.getsubtable(self.vis[self.sortedvisindx[0]],'FIELD')
2707 tb.open(fldtab)
2708 srcdir=tb.getcell('DELAY_DIR',dir)
2709 mrf=tb.getcolkeywords('DELAY_DIR')['MEASINFO']['Ref']
2710 tb.close()
2711 mdir = me.direction(mrf,str(srcdir[0][0])+'rad',str(srcdir[1][0])+'rad')
2712 #tb.open(self.vis+'/OBSERVATION')
2713 obstab=self.getsubtable(self.vis[self.sortedvisindx[0]],'OBSERVATION')
2714 tb.open(obstab)
2715 telname=tb.getcell('TELESCOPE_NAME',0)
2716 # use time in main table instead?
2717 tmr=tb.getcell('TIME_RANGE',0)
2718 tb.close()
2719 #casalog.post...("direction=", me.direction(mrf,str(srcdir[0][0])+'rad',str(srcdir[1][0])+'rad'))
2720 #casalog.post("tmr[1]=" + tmr[1])
2721 #caalog.post...("epoch=", me.epoch('utc',qa.convert(qa.quantity(str(tmr[1])+'s'),'d')))
2722 me.doframe(me.epoch('utc',qa.convert(qa.quantity(str(tmr[0])+'s'),'d')))
2723 me.doframe(me.observatory(telname))
2724 me.doframe(mdir)
2725 f0 = me.frequency(self.dataspecframe, str(fin)+'Hz')
2726 #casalog.post...("frame=", frame, ' f0=',f0)
2727 fout = me.measure(f0,frame)['m0']['value']
2728 return fout
2730 def setspecframe(self,spw):
2731 """
2732 set spectral frame for mfs to data frame based
2733 on spw selection
2734 (part copied from setChannelization)
2735 """
2736 #tb.open(self.vis+'/SPECTRAL_WINDOW')
2737 spectable=self.getsubtable(self.vis[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
2738 tb.open(spectable)
2739 spwframe=tb.getcol('MEAS_FREQ_REF');
2740 tb.close()
2742 # first parse spw parameter:
2744 # use MSSelect if possible
2745 if type(spw)==list:
2746 spw=spw[self.sortedvisindx[0]]
2748 if spw in (-1, '-1', '*', '', ' '):
2749 spw="*"
2751 sel=ms.msseltoindex(self.vis[self.sortedvisindx[0]], spw=spw)
2752 # spw returned by msseletoindex, spw='0:5~10;10~20'
2753 # will give spw=[0] and len(spw) not equal to len(chanids)
2754 # so get spwids from chaninds instead.
2755 chaninds=sel['channel'].tolist()
2756 spwinds=[]
2757 for k in range(len(chaninds)):
2758 spwinds.append(chaninds[k][0])
2759 if(len(spwinds) == 0):
2760 raise Exception('unable to parse spw parameter '+spw)
2762 # the first selected spw
2763 spw0=spwinds[0]
2765 # set dataspecframe:
2766 elspecframe=["REST",
2767 "LSRK",
2768 "LSRD",
2769 "BARY",
2770 "GEO",
2771 "TOPO",
2772 "GALACTO",
2773 "LGROUP",
2774 "CMB"]
2775 self.dataspecframe=elspecframe[spwframe[spw0]];
2776 return
2778 def initChaniter(self,nchan,spw,start,width,imagename,mode,tmpdir='_tmpimdir/'):
2779 """
2780 initialize for channel iteration in interactive clean
2781 --- create a temporary directory, get frequencies for
2782 mode='channel'
2784 Keyword arguments:
2785 nchan -- no. channels
2786 spw -- spw
2787 start -- start modified after channelization function
2788 width -- width modified after channelization function
2789 imagename -- from task input
2790 mode -- from task input
2791 tmpdir -- temporary directory name to store channel images
2793 returns:
2794 frequencies in a list
2795 frequency increment
2796 newmode -- force to set mode to frequency
2797 tmppath -- path for the temporary directory
2798 """
2799 # create a temporary directory to put channel images
2800 tmppath=[]
2801 freqs=[]
2802 finc=0
2803 newmode=mode
2804 for imname in imagename:
2805 if os.path.dirname(imname)=='':
2806 tmppath.append(tmpdir)
2807 else:
2808 tmppath.append(os.path.dirname(imname)+'/'+tmpdir)
2809 # clean up old directory
2810 if os.path.isdir(tmppath[-1]):
2811 shutil.rmtree(tmppath[-1])
2812 os.mkdir(tmppath[-1])
2813 #internally converted to frequency mode for mode='channel'
2814 #to ensure correct frequency axis for output image
2815 #if mode == 'channel':
2816 # freqs, finc = self.getfreqs(nchan, spw, start, width)
2817 # newmode = 'frequency'
2818 if mode == 'channel':
2819 # get spectral axis info from the dirty image
2820 ia.open(imagename[0]+'.image')
2821 imcsys=ia.coordsys().torecord()
2822 ia.close()
2823 # for optical velocity mode, the image will be in tabular form.
2824 if 'tabular' in imcsys['spectral2']:
2825 key='tabular'
2826 else:
2827 key='wcs'
2828 cdelt=imcsys['spectral2'][key]['cdelt']
2829 crval=imcsys['spectral2'][key]['crval']
2830 #cdelt=imcsys['spectral2']['wcs']['cdelt']
2831 #crval=imcsys['spectral2']['wcs']['crval']
2832 for i in range(nchan):
2833 if i==0: freqs.append(crval)
2834 freqs.append(freqs[-1]+cdelt)
2835 finc = cdelt
2836 newmode = 'frequency'
2837 return freqs,finc,newmode,tmppath
2840 def makeTemplateCubes(self, imagename,outlierfile, field, spw, selectdata, timerange,
2841 uvrange, antenna, scan, observation, intent, mode, facets, cfcache, interpolation,
2842 imagermode, localFTMachine, mosweight, locnchan, locstart, locwidth, outframe,
2843 veltype, imsize, cell, phasecenter, restfreq, stokes, weighting,
2844 robust, uvtaper, outertaper, innertaper, modelimage, restoringbeam,
2845 usescratch, noise, npixels, padding):
2846 """
2847 make template cubes to be used for chaniter=T interactive clean
2848 """
2849 imageids=[]
2850 imsizes=[]
2851 phasecenters=[]
2852 rootname=''
2853 multifield=False
2854 loc_modelimage=modelimage
2855 newformat=False
2857 if len(outlierfile) != 0:
2858 f_imageids,f_imsizes,f_phasecenters,f_masks,f_modelimages,parms,newformat=self.newreadoutlier(outlierfile)
2859 if type(imagename) == list or newformat:
2860 rootname = ''
2861 else:
2862 rootname = imagename
2864 # combine with the task parameter input
2865 if type(imagename) == str:
2866 if newformat:
2867 imageids.append(imagename)
2868 imsizes.append(imsize)
2869 phasecenters.append(phasecenter)
2870 else:
2871 imageids=imagename
2872 imsizes=imsize
2873 phasecenters=phasecenter
2875 #if type(mask) != list:
2876 # mask=[mask]
2877 #elif type(mask[0]) != list:
2878 # mask=[mask]
2879 if type(loc_modelimage) != list:
2880 loc_modelimage=[loc_modelimage]
2882 #elif type(loc_modelimage[0]) != list and type(imagename) != str:
2883 #if type(loc_modelimage[0]) != list and \
2884 # (type(imagename) != str or (type(imageids)==list and len(imageids)=1)):
2885 # loc_modelimage=[loc_modelimage]
2886 if type(loc_modelimage[0]) != list:
2887 loc_modelimage=[loc_modelimage]
2889 # now append readoutlier content
2890 for indx, name in enumerate(f_imageids):
2891 imageids.append(name)
2892 imsizes.append(f_imsizes[indx])
2893 phasecenters.append(f_phasecenters[indx])
2895 if newformat:
2896 #mask.append(f_masks[indx])
2897 loc_modelimage.append([f_modelimages[indx]])
2898 else:
2899 if indx!=0:
2900 loc_modelimage.append([f_modelimages[indx]])
2902 ##if len(imageids) > 1:
2903 # multifield=True
2904 else:
2905 imsizes=imsize
2906 phasecenters=phasecenter
2907 imageids=imagename
2909 if len(imageids) > 1:
2910 multifield=True
2912 self.imageids=imageids
2913 # readoutlier need to be run first....
2914 self.datsel(field=field, spw=spw, timerange=timerange, uvrange=uvrange,
2915 antenna=antenna,scan=scan, observation=observation, intent=intent,
2916 usescratch=usescratch,
2917 nchan=-1, start=0, width=1)
2919 self.definemultiimages(rootname=rootname,imsizes=imsizes,cell=cell,
2920 stokes=stokes,mode=mode,
2921 spw=spw, nchan=locnchan, start=locstart,
2922 width=locwidth, restfreq=restfreq,
2923 field=field, phasecenters=phasecenters,
2924 names=imageids, facets=facets,
2925 outframe=outframe, veltype=veltype,
2926 makepbim=False, checkpsf=False)
2928 self.datweightfilter(field=field, spw=spw, timerange=timerange,
2929 uvrange=uvrange, antenna=antenna,scan=scan,
2930 wgttype=weighting, robust=robust, noise=noise,
2931 npixels=npixels, mosweight=mosweight,
2932 uvtaper=uvtaper, innertaper=innertaper, outertaper=outertaper,
2933 usescratch=usescratch, nchan=-1, start=0, width=1)
2934 # split this
2935 #self.datselweightfilter(field=field, spw=spw,
2936 # timerange=timerange, uvrange=uvrange,
2937 # antenna=antenna, scan=scan,
2938 # wgttype=weighting, robust=robust,
2939 # noise=noise, npixels=npixels,
2940 # mosweight=mosweight,
2941 # innertaper=innertaper,
2942 # outertaper=outertaper,
2943 # calready=calready, nchan=-1,
2944 # start=0, width=1)
2946 #localAlgorithm = getAlgorithm(psfmode, imagermode, gridmode, mode,
2947 # multiscale, multifield, facets, nterms,
2948 # 'clark');
2950 #localAlgorithm = 'clark'
2951 #casalog.post("localAlogrithm=" + localAlgorithm)
2953 #self.im.setoptions(ftmachine=localFTMachine,
2954 # wprojplanes=wprojplanes,
2955 # freqinterp=interpolation, padding=padding,
2956 # cfcachedirname=cfcache, pastep=painc,
2957 # epjtablename=epjtable,
2958 # applypointingoffsets=False,
2959 # dopbgriddingcorrections=True)
2960 self.im.setoptions(ftmachine=localFTMachine,
2961 freqinterp=interpolation, padding=padding,
2962 cfcachedirname=cfcache)
2964 modelimages=[]
2965 restoredimage=[]
2966 residualimage=[]
2967 psfimage=[]
2968 fluximage=[]
2969 for k in range(len(self.imagelist)):
2970 ia.open(self.imagelist[k])
2971 #if (modelimage =='' or modelimage==[]) and multifield:
2972 # ia.rename(self.imagelist[k]+'.model',overwrite=True)
2973 #else:
2974 # ia.remove(verbose=False)
2975 if ((loc_modelimage =='' or loc_modelimage==[]) or \
2976 (type(loc_modelimage)==list and \
2977 (loc_modelimage[k]=='' or loc_modelimage[k]==[''] or loc_modelimage[k]==[]))) and multifield:
2978 ia.rename(self.imagelist[k]+'.model',overwrite=True)
2979 else:
2980 modlist=[]
2981 if type(modelimage)==str:
2982 modlist=[modelimage]
2983 # make sure input model image is not removed
2984 if (not any([inmodel == self.imagelist[k] for inmodel in modlist])) and \
2985 (not any([inmodel == self.imagelist[k]+'.model' for inmodel in modlist])):
2986 # ia.remove(verbose=False)
2987 ia.rename(self.imagelist[k]+'.model',overwrite=True)
2988 ia.close()
2990 modelimages.append(self.imagelist[k]+'.model')
2991 restoredimage.append(self.imagelist[k]+'.image')
2992 residualimage.append(self.imagelist[k]+'.residual')
2993 psfimage.append(self.imagelist[k]+'.psf')
2994 if(imagermode=='mosaic'):
2995 fluximage.append(self.imagelist[k]+'.flux')
2997 # make dirty image cube
2998 if multifield:
2999 alg='mfclark'
3000 else:
3001 alg='clark'
3003 self.im.clean(algorithm=alg, niter=0,
3004 model=modelimages, residual=residualimage,
3005 image=restoredimage, psfimage=psfimage,
3006 mask='', interactive=False)
3009 def setChaniterParms(self,finalimagename, spw,chan,start,width,freqs,finc,tmppath):
3010 """
3011 Set parameters for channel by channel iterations
3012 returns:
3013 start and width to define each channel image plane
3014 """
3015 retparms={}
3016 self.maskimages={}
3017 retparms['imagename']=[tmppath[indx]+os.path.basename(imn)+'.ch'+str(chan)
3018 for indx, imn in enumerate(finalimagename)]
3020 #casalog.post("Processing channel %s " % chan)
3021 #self._casalog.post("Processing channel %s "% chan)
3023 # Select only subset of vis data if possible.
3024 # It does not work well for multi-spw so need
3025 # to select with nchan=-1
3026 retparms['imnchan']=1
3027 retparms['chanslice']=chan
3029 q = qa.quantity
3031 # 2010-08-18 note: disable this. Has the problem
3032 # getting imaging weights correctly when the beginning
3033 # channels were flagged.
3034 #if type(spw)==int or len(spw)==1:
3035 # if width>1:
3036 # visnchan=width
3037 # else:
3038 # visnchan=1
3039 #else:
3040 # visnchan=-1
3042 visnchan=-1
3043 retparms['visnchan']=visnchan
3044 visstart=0
3046 if type(start)==int:
3047 # need to convert to frequencies
3048 # to ensure correct frequencies in
3049 # output images(especially for multi-spw)
3050 # Use freq list instead generated in initChaniter
3051 imstart=q(freqs[chan],'Hz')
3052 width=q(finc,'Hz')
3053 elif start.find('m/s')>0:
3054 imstart=qat.add(q(start),qat.mul(chan,q(width)))
3055 elif start.find('Hz')>0:
3056 imstart=qat.add(q(start),qat.mul(chan,q(width)))
3057 retparms['width']=width
3058 retparms['imstart']=imstart
3059 retparms['visstart']=visstart
3061 #
3062 return retparms
3064 def defineChaniterModelimages(self,modelimage,chan,tmppath):
3065 """
3066 chaniter=T specific function to convert input models
3067 to a model image
3068 """
3069 chanmodimg=[]
3070 if type(modelimage)==str:
3071 modelimage=[modelimage]
3072 indx=0
3073 for modimg in modelimage:
3074 if modimg=='':
3075 return
3076 if type(modimg)==list:
3077 chanmodimg=[]
3078 for img in modimg:
3079 if img!='':
3080 if os.path.dirname(img) != '':
3081 chanmodimg.append(tmppath[0] + '_tmp.' +
3082 os.path.basename(img))
3083 else:
3084 chanmodimg.append(tmppath[0] + '_tmp.' + img)
3085 self.getchanimage(cubeimage=img, outim=chanmodimg[-1], chan=chan)
3086 #self.convertmodelimage(modelimages=chanmodimg,
3087 # outputmodel=self.imagelist.values()[0]+'.model')
3088 self.convertmodelimage(modelimages=chanmodimg,
3089 outputmodel=self.imagelist.values()[indx]+'.model', imindex=indx)
3090 chanmodimg=[]
3091 indx+=1
3092 else:
3093 if os.path.dirname(modimg) != '':
3094 chanmodimg.append(tmppath[0] + '_tmp.' + os.path.basename(modimg))
3095 else:
3096 chanmodimg.append(tmppath[0] + '_tmp.' + modimg)
3097 self.getchanimage(cubeimage=modimg, outim=chanmodimg[-1],chan=chan)
3099 #self.convertmodelimage(modelimages=chanmodimg,
3100 # outputmodel=self.imagelist.values()[0]+'.model')
3101 self.convertmodelimage(modelimages=chanmodimg,
3102 outputmodel=self.imagelist.values()[indx]+'.model',imindex=indx)
3103 # clean up temporary channel model image
3104 self.cleanupTempFiles(chanmodimg)
3106 def convertAllModelImages_old(self,modelimage, mode, nterms, dochaniter, chan, tmppath):
3107 """
3108 wrapper function for convertmodelimage for all different cases
3109 """
3110 if (type(modelimage)!=str and type(modelimage)!=list):
3111 raise Exception('modelimage must be a string or a list of strings')
3112 #spectralline modes
3113 if (not mode=='mfs') or (mode=='mfs' and nterms==1):
3114 if (not all(img=='' or img==[] or img==[''] for img in modelimage)):
3115 if dochaniter:
3116 self.defineChaniterModelimages(modelimage,chan,tmppath)
3117 else:
3118 if type(modelimage)== str or \
3119 (type(modelimage)==list and len(self.imagelist)==1 and len(modelimage)>1):
3120 modelimage=[modelimage]
3122 #casalog.post("Run convertmodelimage for this list : " + self.imagelist + " with these models : " + modelimage;)
3123 for j in range(len(self.imagelist)):
3124 self._casalog.post("Use modelimages: "+str(modelimage[j])+" to create a combined modelimage: " \
3125 +self.imagelist.values()[j]+".model", 'DEBUG1')
3126 if modelimage[j] != '' and modelimage[j] != []:
3127 self.convertmodelimage(modelimages=modelimage[j],
3128 outputmodel=self.imagelist.values()[j]+'.model',imindex=j)
3130 # elif .......
3131 # put mfs with nterms>1 case here
3133 ##########################################################
3134 # Multiple models for one field : [ [ 'm0', 'm1' ] ]
3135 # Multiple taylor terms and one field : [ [ 't0','t1'] ]
3136 # Multiple models per field : [ [ 'm0f0', 'm1f0' ] , [ 'm0f1', 'm1f1' ] ]
3137 # Multiple taylor terms per field : [ [ 't0f0','t1f0' ] , [ 't0f1','t1f1' ] ]
3138 ##########################################################
3139 # Cannot do multiple models per taylor term and per field for now.
3140 # ....... later... [ [ ['m0t0f0','m1t0f0'],['m0t1f0','m1t1f0'] ] , [ [ ['t0f1'] ],[ ['t1f1'] ] ] ]
3141 ##########################################################
3142 def convertAllModelImages(self,modelimage, mode, nterms, dochaniter, chan, tmppath):
3143 """
3144 wrapper function for convertmodelimage for all different cases
3145 """
3146 if (type(modelimage)!=str and type(modelimage)!=list):
3147 raise Exception('modelimage must be a string or a list of strings')
3148 if (not all(img=='' or img==[] or img==[''] for img in modelimage)):
3149 if dochaniter:
3150 self.defineChaniterModelimages(modelimage,chan,tmppath)
3151 else:
3152 if type(modelimage)== str or \
3153 (type(modelimage)==list and len(self.imagelist)==1 and len(modelimage)>1):
3154 modelimage=[modelimage]
3156 #casalog.post convertmodelimage for this list : " + self.imagelist + " with these models : " + modelimage;)
3157 #spectralline modes + basic mfs
3158# if (not mode=='mfs') or (mode=='mfs' and nterms==1):
3159# for j in range(len(self.imagelist)): # = nfield
3160# self._casalog.post("Use modelimages: "+str(modelimage[j])+" to create a combined modelimage: " \
3161# +self.imagelist.values()[j]+".model", 'DEBUG1')
3162# if modelimage[j] != '' and modelimage[j] != []:
3163# self.convertmodelimage(modelimages=modelimage[j],
3164# outputmodel=self.imagelist.values()[j]+'.model',imindex=j)
3166# else: # mfs and nterms>1
3167 if 1:
3168 nfld = len(self.imagelist);
3169 # if only one field, then modelimage must be a list of strings. convert to list of list of str
3170 # if multiple fields, then model image : list of list of strings
3171 if nfld != len(modelimage):
3172 raise Exception('Model images must be same length as fields : '+ str(nfld) + str(modelimage))
3174 for fld in range(nfld):
3175 modsforfield = modelimage[fld]; # a list
3176 if type(modsforfield)==str:
3177 modsforfield = [modsforfield];
3178 if nterms==1:
3179 nimages = len(modsforfield);
3180 else:
3181 nimages = min( len(modsforfield), nterms ); ## one model per term
3182 for tt in range(0,nimages):
3183 if nterms==1:
3184 modname = self.imagelist[fld]+'.model';
3185 else:
3186 modname = self.imagelist[fld]+'.model.tt'+str(tt) ;
3187 if( os.path.exists(modsforfield[tt]) ):
3188# casalog.post("Found user-specified model image : "+modsforfield[tt]+" . Adding to starting model : "+modname;)
3189 self._casalog.post("Found user-specified model image : "+modsforfield[tt]+" . Adding to starting model : "+modname);
3190 self.convertmodelimage(modelimages=modsforfield[tt],outputmodel=modname, imindex=fld);
3191 else:
3192 self._casalog.post("Cannot find user-specified model image : "+modsforfield[tt]+" . Continuing with current model : "+modname);
3198 def storeCubeImages(self,cubeimageroot,chanimageroot,chan,imagermode):
3199 """
3200 Put channel images back into CubeImages for chaniter=T mode
3202 Keyword arguments:
3203 cubeimageroot -- root name for output cube image
3204 chanimageroot -- root name for channel image
3205 chan -- channel plane index
3206 imagermode -- imagermode
3207 """
3208 imagext = ['.image','.model','.flux','.residual','.psf','.mask']
3209 if imagermode=='mosaic':
3210 imagext.append('.flux.pbcoverage')
3211 lerange=range(self.nimages)
3212 for n in lerange:
3213 cubeimagerootname=cubeimageroot[n]
3214 chanimagerootname=chanimageroot[n]
3215 for ext in imagext:
3216 nomaskim=False
3217 cubeimage=cubeimagerootname+ext
3218 chanimage=chanimagerootname+ext
3219 if not os.path.exists(cubeimage):
3220 if os.path.exists(chanimage):
3221 outim=ia.newimagefromimage(cubeimagerootname+'.model',cubeimage)
3222 outim.done(verbose=False)
3223 elif ext=='.mask':
3224 # unless mask image is given or in interactive mode
3225 # there is no mask image
3226 nomaskim=True
3227 if not nomaskim:
3228 self.putchanimage(cubeimage, chanimage,chan)
3230 def cleanupTempFiles(self, tmppath):
3231 """
3232 Remove the directories listed by tmppath.
3233 """
3234 # Created to deal with temporary dirs created by chaniter=T clean,
3235 # now used elsewhere too.
3236 for dir in tmppath:
3237 if os.path.exists(dir):
3238 shutil.rmtree(dir)
3240 def convertImageFreqFrame(self,imlist):
3241 """
3242 Convert output images to proper output frame
3243 (after im.clean() executed)
3244 """
3245 if type(imlist)==str:
3246 imlist=[imlist]
3247 if self.usespecframe.lower() != 'lsrk':
3248 if self.usespecframe=='':
3249 inspectral=self.dataspecframe
3250 else:
3251 inspectral=self.usespecframe
3252 for img in imlist:
3253 if os.path.exists(img):
3254 ia.open(img)
3255 csys=ia.coordsys()
3256 csys.setconversiontype(spectral=inspectral)
3257 #casalog.post("csys.torecord spectral2=" + csys.torecord()['spectral2'])
3258 ia.setcoordsys(csys.torecord())
3259 ia.close()
3261 def setFrameConversionForMasks(self):
3262 ''' To be called at the end of clean, so that the output csys can be
3263 read and set for the mask. This will have the users desired
3264 conversion layer '''
3265 if self.usespecframe=='':
3266 useframe=self.dataspecframe
3267 else:
3268 useframe=self.usespecframe
3270 #casalog.post('maskimages.keys : ' + self.maskimages.keys())
3271 #casalog.post('imagelist : ' + self.imagelist)
3273 for key in self.imagelist.keys():
3274 imgmask = self.imagelist[key]+'.mask'
3275 img = self.imagelist[key]+'.image'
3276 if not os.path.exists(img):
3277 img = img+'.tt0'
3278 # casalog.post('Converting frame for ' + imgmask + ' to ' + useframe)
3279 if os.path.exists(imgmask) and os.path.exists(img):
3280 ia.open(img)
3281 imcsys = ia.coordsys()
3282 ia.close()
3283 ia.open(imgmask)
3284# csys=ia.coordsys()
3285# csys.setreferencecode('LSRK','spectral',True)
3286# val = csys.setconversiontype(spectral=useframe)
3287# casalog.post('Ret val : ' + val + csys.getconversiontype('spectral'))
3288# ia.setcoordsys(csys.torecord())
3289# casalog.post('conv type : '+ imcsys.getconversiontype('spectral'))
3290 ia.setcoordsys( imcsys.torecord() )
3291 ia.close()
3292 else:
3293 self._casalog.post('Not converting spectral reference frame for mask image','DEBUG1')
3296 def setReferenceFrameLSRK(self, img = ''):
3297 ''' To be called to reset reference code and conversion layer to LSRK '''
3298 if os.path.exists( img ):
3299 ia.open( img )
3300 mycsys=ia.coordsys()
3301 if (mycsys.torecord()['spectral2']['conversion']['system']=='REST') :
3302 ia.close()
3303 return
3304 if (mycsys.torecord()['spectral2']['conversion']['system']!='LSRK') :
3305 mycsys.setreferencecode('LSRK','spectral',True)
3306 mycsys.setconversiontype(spectral='LSRK')
3307 ia.setcoordsys( mycsys.torecord() )
3308 ia.close()
3310 def resmooth(self, model, residual, restored, minOrMax):
3311 if(minOrMax=="common"):
3312 ia.open(restored)
3313 beam=ia.restoringbeam();
3314 if('nChannels' not in beam):
3315 return
3316 combeam=ia.commonbeam()
3317 ia.done()
3318 ia.fromimage(outfile='__restored-copy', infile=restored, overwrite=True)
3319 ia.open('__restored-copy')
3320 ib=ia.convolve2d(outfile=restored, major='', minor='', pa='', targetres=True, beam=combeam, overwrite=True)
3321 ib.done()
3322 ia.remove()
3323 ia.done()
3324 ia.fromimage(outfile='__residual-copy', infile=residual, overwrite=True)
3325 ia.open('__residual-copy')
3326 ###need to set a beam first to go around CAS-5433 and then loop
3327 ia.setrestoringbeam(major=beam['beams']['*0']['*0']['major'], minor=beam['beams']['*0']['*0']['minor'], pa=beam['beams']['*0']['*0']['positionangle'], channel=0, polarization=0)
3328 nchan=beam['nChannels']
3329 for k in range(nchan):
3330 chstr='*'+str(k)
3331 ia.setrestoringbeam(beam=beam['beams'][chstr]['*0'], channel=k, polarization=0)
3332 ib=ia.convolve2d(outfile=residual, major='', minor='', pa='', targetres=True, beam=combeam, overwrite=True)
3333 ib.done()
3334 ia.remove()
3335 ia.done()
3336 return
3338 ###############for min or max
3339 ia.open(restored)
3340 beams=ia.restoringbeam()
3341 if('beams' not in beams):
3342 ########already has one beam only
3343 ia.done()
3344 return
3345 minArea=1e37
3346 maxArea=-1e37
3347 maxchan=-1
3348 minchan=-1
3349 theArea=numpy.zeros(beams['nChannels'])
3350 for k in range(beams['nChannels']):
3351 ##it must have been really hard to provide proper indices
3352 theArea[k]=qa.convert(beams['beams']['*'+str(k)]['*0']['major'], 'arcsec')['value'] * qa.convert(beams['beams']['*'+str(k)]['*0']['minor'], 'arcsec')['value']
3353 if(theArea[k] > maxArea):
3354 maxArea=theArea[k]
3355 maxchan=k
3356 if(theArea[k] < minArea):
3357 minArea=theArea[k]
3358 minchan=k
3359 maxbeam=[beams['beams']['*'+str(maxchan)]['*0']['major'], beams['beams']['*'+str(maxchan)]['*0']['minor'], beams['beams']['*'+str(maxchan)]['*0']['positionangle']]
3360 minbeam=[beams['beams']['*'+str(minchan)]['*0']['major'], beams['beams']['*'+str(minchan)]['*0']['minor'], beams['beams']['*'+str(minchan)]['*0']['positionangle']]
3361 thebeam=minbeam
3362 tobeDiv=theArea[minchan]
3363 if(minOrMax=='max'):
3364 thebeam=maxbeam
3365 tobeDiv=theArea[maxchan]
3366 ia.open(residual)
3367 shp=ia.shape()
3368 for k in range(beams['nChannels']):
3369 reg=rg.box(blc=[0,0,0,k], trc=[shp[0]-1, shp[1]-1, shp[2]-1, k])
3370 pix=ia.getregion(region=reg)
3371 pix=pix*theArea[k]/tobeDiv
3372 ia.putregion(pixels=pix, region=reg)
3373 ia.done()
3374 ia.open(model)
3375 ib=ia.convolve2d(outfile=restored, axes=[0,1], major=thebeam[0], minor=thebeam[1], pa=thebeam[2], overwrite=True)
3376 ib.calc('"'+restored+'" + '+'"'+residual+'"')
3377 ib.done()
3378 ia.done()
3384 @staticmethod
3385 def getOptimumSize(size):
3386 '''
3387 This returns the next largest even composite of 2, 3, 5, 7
3388 '''
3389 def prime_factors(n, douniq=True):
3390 """ Return the prime factors of the given number. """
3391 factors = []
3392 lastresult = n
3393 sqlast=int(numpy.sqrt(n))+1
3394 # 1 pixel must a single dish user
3395 if n == 1:
3396 return [1]
3397 c=2
3398 while 1:
3399 if (lastresult == 1) or (c > sqlast):
3400 break
3401 sqlast=int(numpy.sqrt(lastresult))+1
3402 while 1:
3403 if(c > sqlast):
3404 c=lastresult
3405 break
3406 if lastresult % c == 0:
3407 break
3408 c += 1
3410 factors.append(c)
3411 lastresult //= c
3412 if(factors==[]): factors=[n]
3413 return numpy.unique(factors).tolist() if douniq else factors
3414 n=size
3415 if(n%2 != 0):
3416 n+=1
3417 fac=prime_factors(n, False)
3418 for k in range(len(fac)):
3419 if(fac[k] > 7):
3420 val=fac[k]
3421 while(numpy.max(prime_factors(val)) > 7):
3422 val +=1
3423 fac[k]=val
3424 newlarge=numpy.product(fac)
3425 for k in range(n, newlarge, 2):
3426 if((numpy.max(prime_factors(k)) < 8)):
3427 return k
3428 return newlarge
3431def getFTMachine(gridmode, imagermode, mode, wprojplanes, userftm):
3432 """
3433 A utility function which implements the logic to determine the
3434 ftmachine name to be used in the under-laying tool.
3435 """
3436# ftm = userftm;
3437 ftm='ft';
3438 if((gridmode != '') and (imagermode=='mosaic')):
3439 raise Exception(pwd.getpwuid(os.getuid())[4].split()[0]+ " gridmode='"+gridmode + "' combined with " + "imagermode='"+imagermode + "' is not supported yet")
3440 if ((gridmode == 'widefield') and(wprojplanes > 1 or wprojplanes==-1)): ftm = 'wproject';
3441 elif (gridmode == 'aprojection'): ftm = 'awproject';
3442 elif (gridmode == 'advancedaprojection'): ftm = 'awproject';
3443 elif (imagermode == 'csclean'): ftm = 'ft';
3444 elif (imagermode == 'mosaic'): ftm = userftm;
3445 return ftm;
3447def getAlgorithm(psfmode, imagermode, gridmode, mode,
3448 multiscale, multifield, facets, nterms, useralg):
3449 """
3450 A utility function which implements the logic to determine the
3451 deconvolution algorithm to be used in the under-laying tool.
3452 """
3453 alg=useralg
3454 addMultiField=False;
3456 if((type(multiscale)==list) and
3457 (len(multiscale) > 0) and
3458 (sum(multiscale) > 0)): alg = 'multiscale';
3459 elif ((psfmode == 'clark') or (psfmode == 'clarkstokes') or (psfmode == 'hogbom')): alg=psfmode;
3461 if ((imagermode == '') and (multifield)): addMultiField=True;
3462 if (imagermode == 'mosaic'): addMultiField=True;
3463 if (imagermode == 'csclean'): addMultiField = True; #!!!!
3465 if ((mode == 'mfs') and (nterms > 1)):
3466 alg = 'msmfs';
3467 if(imagermode == 'mosaic'):
3468 ##casalog.post('Multi-Term MFS with a mosaic is experimental')
3469 raise Exception('msmfs (nterms>1) not allowed with imagermode=' + imagermode + '. For now, msmfs automatically performs cs-clean type iterations')
3470 if (multifield):
3471 addMultiField = True;
3472 if facets > 1:
3473 raise Exception('msmfs (nterms>1) with facets>1 is not yet available')
3474 if( (mode=='mfs') and (nterms<1) ):
3475 raise Exception('nterms must be > 0')
3477# if (gridmode == 'widefield'): alg='mfclark';
3479 if (gridmode == 'widefield'):
3480 addMultiField=True;
3481 if (facets > 1):
3482 if(alg.count('multiscale') > 0):
3483 raise Exception('multiscale with facets > 1 not allowed for now')
3484 if (psfmode==''): psfmode='clark';
3485 if((psfmode == 'clark') or (psfmode == 'hogbom')):
3486 alg='wf'+psfmode;
3487 addMultiField=False;
3488 else:
3489 addMultiField=True;
3490# addMultiField=False;
3492#
3493# if facets > 1 && mutliscale ==> fail
3496 if (addMultiField and (alg[0:2] != 'mf') and (alg != 'msmfs')): alg = 'mf' + alg;
3497 return alg;
3499def convert_numpydtype(listobj):
3500 """
3501 utility function to covert list with elements in numpy.int or
3502 numpy.float types to python int/float
3503 """
3504 import array as pyarr
3505 floatarr=False
3506 intarr=False
3507 for elm in listobj:
3508 if numpy.issubdtype(type(elm), numpy.float):
3509 floatarr = True
3510 elif numpy.issubdtype(type(elm), numpy.int):
3511 intarr = True
3512 if floatarr or (floatarr and intarr):
3513 temparr=pyarr.array('f', listobj)
3514 elif intarr:
3515 temparr=pyarr.array('i', listobj)
3516 else:
3517 temparr = listobj
3518 return temparr
3519 return temparr.tolist()
3521def get_func_params(func, loc_vars):
3522 """ returns a dictionary of parameter name:vale for all the parameters of a function
3524 :param func: function object (for example a task function)
3525 :param loc_vars: locals() from inside the function.
3526 """
3527 param_names = func.__code__.co_varnames[:func.__code__.co_argcount]
3528 params = [(name, loc_vars[name]) for name in param_names]
3529 return params
3531def write_task_history(images, tname, params, clog):
3532 """
3533 Update image/logtable with the taskname, CASA version and all task parameters
3534 CASR-571. Replicates the same format as mstools.write_history.
3536 :param images: list of images to write the history to
3537 :param tname: task name to use as origin of the history
3538 :param params: list of task parameters as a tuple (name, value)
3539 :param clog: casa log object
3540 """
3541 iat = iatool()
3543 history = ['taskname={0}'.format(tname)]
3544 history.append(_casa_version_string())
3545 # Add all task arguments.
3546 for name, val in params:
3547 msg = "%-11s = " % name
3548 if type(val) == str:
3549 msg += '"'
3550 msg += str(val)
3551 if type(val) == str:
3552 msg += '"'
3553 history.append(msg)
3555 for img in images:
3556 iat_open = False
3557 try:
3558 iat.open(img)
3559 iat_open = True
3560 iat.sethistory(origin=tname, history=history)
3561 except RuntimeError:
3562 clog.post('Could not open this directory as an image to write history: {}'.
3563 format(img), 'DEBUG')
3564 finally:
3565 if iat_open:
3566 iat.close()
3567 iat.done()
3569def write_tclean_history(imagename, tname, params, clog):
3570 """
3571 Update image/logtable with the taskname, CASA version and all task parameters
3572 CASR-571. Replicates the same format as mstools.write_history.
3574 :param imagename: imagename prefi as used in tclean
3575 :param tname: task name to use as origin of the history
3576 :param params: list of task parameters as a tuple (name, value)
3577 :param clog: casa log object
3578 """
3580 def filter_img_names(img_exts):
3581 """
3582 Applies name (extension) exclusions to not try to open tclean outputs files
3583 and/or leftovers that are not images (even if they share the same name as the
3584 proper images).
3586 Some of the files excluded here can cause spurious messages (image tool will
3587 print a SEVERE error to the log when open fails) if for example while running
3588 several test cases one of them fails or misbehaves leaving temporary files
3589 behind.
3591 :param img_exts: list of image names (different extensions)
3592 :returns: list of image names after filtering out undesired ones
3593 """
3594 accept = []
3595 regex = re.compile('^' + re.escape(imagename) + '[0-9]*_?[0-9]*\..+')
3596 for img in img_exts:
3597 if img.endswith(('.cf', '.cfcache', '.workdirectory', '.work.temp', '.txt')):
3598 continue
3600 # ensure only 'imgname' + optional integer + .*
3601 res = re.match(regex, img)
3602 if res:
3603 accept.append(img)
3604 return accept
3606 def filter_obvious_nonimages(img_exts):
3607 """
3608 Try to filter out files that are not images but have been placed in the same
3609 directory and share the same prefix name as the images.
3610 For example, images have to be directories. All non-dir files can be filtered
3611 out. It also checks for a logtable subdirectory with a table.info file, which
3612 is expected in CASA images.
3614 This is to not even try to open them (with the iatool or similar).
3615 See CAS-13464 for additional complications around tclean output image names.
3617 :param img_exts: list of image names (different extensions)
3618 :returns: list of image names after filtering out the ones that do not seem
3619 to be images.
3620 """
3621 path_check = ['logtable', 'table.info']
3622 accept = [img for img in img_exts if
3623 os.path.isdir(img) and os.path.isfile(os.path.join(img, *path_check))]
3624 return accept
3626 img_exts = glob.glob(imagename + '*.*')
3627 img_exts = filter_img_names(img_exts)
3628 img_exts = filter_obvious_nonimages(img_exts)
3629 clog.post("Searching for images with prefix '{}'... Found these, writing history "
3630 "into them: {}".format(imagename, img_exts))
3631 write_task_history(img_exts, tname, params, clog)