|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os 00002 import shutil 00003 from taskinit import * 00004 00005 def cvel(vis, outputvis, 00006 passall, field, spw, selectdata, antenna, timerange, scan, array, 00007 mode, nchan, start, width, interpolation, 00008 phasecenter, restfreq, outframe, veltype, hanning): 00009 00010 """ regrid an MS to a new spectral window / channel structure or frame 00011 00012 vis -- Name of input visibility file 00013 default: none; example: vis='ngc5921.ms' 00014 00015 outputvis -- Name of output measurement set (required) 00016 default: none; example: vis='ngc5921-regridded.ms' 00017 00018 passall -- if False, data not meeting the selection on field or spw 00019 is omitted; if True, data not meeting the selection 00020 is passed through without modification 00021 default: False; example: 00022 field='NGC5921' 00023 passall=False : only data from NGC5921 is included in output MS, 00024 no data from other fields (e.g. 1331+305) is included 00025 passall=True : data from NGC5921 is transformed by cvel, all other 00026 fields are passed through unchanged 00027 00028 field -- Select fields in MS. Use field id(s) or field name(s). 00029 ['go listobs' to obtain the list id's or names] 00030 default: ''= all fields 00031 If field string is a non-negative integer, it is assumed to 00032 be a field index otherwise, it is assumed to be a 00033 field name 00034 field='0~2'; field ids 0,1,2 00035 field='0,4,5~7'; field ids 0,4,5,6,7 00036 field='3C286,3C295'; field named 3C286 and 3C295 00037 field = '3,4C*'; field id 3, all names starting with 4C 00038 00039 spw --Select spectral window/channels 00040 NOTE: This selects the data passed as the INPUT to mode 00041 default: ''=all spectral windows and channels 00042 spw='0~2,4'; spectral windows 0,1,2,4 (all channels) 00043 spw='0:5~61'; spw 0, channels 5 to 61 00044 spw='<2'; spectral windows less than 2 (i.e. 0,1) 00045 spw='0,10,3:3~45'; spw 0,10 all channels, spw 3, 00046 channels 3 to 45. 00047 spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each. 00048 spw='0:0~10;15~60'; spectral window 0 with channels 00049 0-10,15-60 00050 spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10, 00051 spw 1, channels 20-30, and spw 2, channels, 1,2 and 3 00052 00053 selectdata -- Other data selection parameters 00054 default: True 00055 00056 selectdata=True expandable parameters 00057 00058 antenna -- Select data based on antenna/baseline 00059 default: '' (all) 00060 If antenna string is a non-negative integer, it is 00061 assumed to be an antenna index, otherwise, it is 00062 considered an antenna name. 00063 antenna='5&6'; baseline between antenna index 5 and 00064 index 6. 00065 antenna='VA05&VA06'; baseline between VLA antenna 5 00066 and 6. 00067 antenna='5&6;7&8'; baselines 5-6 and 7-8 00068 antenna='5'; all baselines with antenna index 5 00069 antenna='05'; all baselines with antenna number 05 00070 (VLA old name) 00071 antenna='5,6,9'; all baselines with antennas 5,6,9 00072 index numbers 00073 00074 timerange -- Select data based on time range: 00075 default = '' (all); examples, 00076 timerange = 'YYYY/MM/DD/hh:mm:ss~YYYY/MM/DD/hh:mm:ss' 00077 Note: if YYYY/MM/DD is missing date defaults to first 00078 day in data set 00079 timerange='09:14:0~09:54:0' picks 40 min on first day 00080 timerange= '25:00:00~27:30:00' picks 1 hr to 3 hr 00081 30min on NEXT day 00082 timerange='09:44:00' pick data within one integration 00083 of time 00084 timerange='>10:24:00' data after this time 00085 00086 scan -- Scan number range. 00087 default: '' (all) 00088 example: scan='1~5' 00089 Check 'go listobs' to insure the scan numbers are in 00090 order. 00091 00092 array -- Select data by (sub)array indices 00093 default: '' (all); example: 00094 array='0~2'; arrays 0 to 2 00095 00096 mode -- Frequency Specification: 00097 NOTE: See examples below: 00098 default: 'channel' 00099 mode = 'channel'; Use with nchan, start, width to specify 00100 output spw. See examples below 00101 mode = 'velocity', means channels are specified in 00102 velocity. 00103 mode = 'frequency', means channels are specified in 00104 frequency. 00105 00106 mode expandable parameters 00107 Start, width are given in units of channels, frequency 00108 or velocity as indicated by mode 00109 nchan -- Number of channels in output spw 00110 default: -1 = all channels; example: nchan=3 00111 start -- Start input channel (relative-0) 00112 default=0; example: start=5 00113 width -- Output channel width in units of the input 00114 channel width (>1 indicates channel averaging) 00115 default=1; example: width=4 00116 interpolation -- Interpolation method 00117 default = 'linear' 00118 00119 examples: 00120 spw = '0,1'; mode = 'channel' 00121 will produce a single spw containing all channels in spw 00122 0 and 1 00123 spw='0:5~28^2'; mode = 'channel' 00124 will produce a single spw made with channels 00125 (5,7,9,...,25,27) 00126 spw = '0'; mode = 'channel': nchan=3; start=5; width=4 00127 will produce an spw with 3 output channels 00128 new channel 1 contains data from channels (5+6+7+8) 00129 new channel 2 contains data from channels (9+10+11+12) 00130 new channel 3 contains data from channels (13+14+15+16) 00131 spw = '0:0~63^3'; mode='channel'; nchan=21; start = 0; 00132 width = 1 00133 will produce an spw with 21 channels 00134 new channel 1 contains data from channel 0 00135 new channel 2 contains data from channel 2 00136 new channel 21 contains data from channel 61 00137 spw = '0:0~40^2'; mode = 'channel'; nchan = 3; start = 00138 5; width = 4 00139 will produce an spw with three output channels 00140 new channel 1 contains channels (5,7) 00141 new channel 2 contains channels (13,15) 00142 new channel 3 contains channels (21,23) 00143 00144 phasecenter -- direction measure or fieldid for the mosaic center 00145 default: '' => first field selected ; example: phasecenter=6 00146 or phasecenter='J2000 19h30m00 -40d00m00' 00147 00148 restfreq -- Specify rest frequency to use for output image 00149 default='' Occasionally it is necessary to set this (for 00150 example some VLA spectral line data). For example for 00151 NH_3 (1,1) put restfreq='23.694496GHz' 00152 00153 outframe -- output reference frame 00154 default='' (keep original reference frame) ; example: outframe='bary' 00155 00156 veltype -- definition of velocity (in mode) 00157 default = 'radio' 00158 00159 hanning -- if true, Hanning smooth frequency channel data before regridding 00160 to remove Gibbs ringing 00161 default = False 00162 00163 """ 00164 00165 00166 #Python script 00167 00168 # make ms and tb tool local 00169 ms = casac.ms() 00170 tb = casac.table() 00171 00172 try: 00173 casalog.origin('cvel') 00174 00175 parsummary = 'vis="'+str(vis)+'", outputvis="'+str(outputvis)+'", passall="'+str(passall)+'", ' 00176 parsummary += 'field="'+str(field)+'", spw='+str(spw)+', antenna="'+str(antenna)+'", timerange="' 00177 parsummary += str(timerange)+'",' 00178 casalog.post(parsummary,'INFO') 00179 parsummary = 'mode='+str(mode)+', nchan='+str(nchan)+', start='+str(start)+', width='+str(width) 00180 parsummary += ', interpolation = '+str(interpolation)+', outframe="'+str(outframe)+'",' 00181 casalog.post(parsummary,'INFO') 00182 parsummary = 'phasecenter="'+str(phasecenter)+'", restfreq="'+str(restfreq)+'", veltype="'+str(veltype) 00183 parsummary += '", hanning="'+str(hanning)+'"' 00184 casalog.post(parsummary,'INFO') 00185 00186 00187 if not ((type(vis)==str) & (os.path.exists(vis))): 00188 raise Exception, 'Visibility data set not found - please verify the name' 00189 00190 if (outputvis == ""): 00191 raise Exception, "Must provide output data set name in parameter outputvis." 00192 00193 if os.path.exists(outputvis): 00194 raise Exception, "Output MS %s already exists - will not overwrite." % outputvis 00195 00196 # Handle selectdata explicitly 00197 # (avoid hidden globals) 00198 if (selectdata==False): 00199 timerange='' 00200 array='' 00201 antenna='' 00202 scan='' 00203 00204 if(type(antenna) == list): 00205 antenna = ', '.join([str(ant) for ant in antenna]) 00206 00207 if (field == ''): 00208 field = '*' 00209 00210 if (spw == ''): 00211 spw = '*' 00212 00213 if(passall and spw=='*' and field=='*'): 00214 # all spws and fields selected, nothing to pass through 00215 passall = False 00216 00217 doselection = True 00218 if(field=='*' and spw=='*' and 00219 (not selectdata or (selectdata and antenna=='' and timerange=='' and scan=='' and array=='')) 00220 ): 00221 doselection = False 00222 00223 # open input MS 00224 ms.open(vis) 00225 00226 if(type(phasecenter)==str): 00227 ### blank means take field 0 00228 if (phasecenter==''): 00229 phasecenter=ms.msseltoindex(vis,field=field)['field'][0] 00230 else: 00231 tmppc=phasecenter 00232 try: 00233 if(len(ms.msseltoindex(vis, field=phasecenter)['field']) > 0): 00234 tmppc=ms.msseltoindex(vis, field=phasecenter)['field'][0] 00235 ##succesful must be string like '0' or 'NGC*' 00236 except Exception, instance: 00237 ##failed must be a string 'J2000 18h00m00 10d00m00' 00238 tmppc=phasecenter 00239 phasecenter=tmppc 00240 00241 newphasecenter = phasecenter 00242 phasecentername = phasecenter 00243 if not (type(phasecenter)==str): 00244 # first check that this field will be in the output MS 00245 if not (phasecenter in ms.msseltoindex(vis,field=field)['field']): 00246 message = "Field id " + str(phasecenter) 00247 message += " was selected as phasecenter but is not among the fields selected for output: " 00248 message += str(ms.msseltoindex(vis,field=field)['field']) 00249 ms.close() 00250 raise Exception, message 00251 00252 tb.open(vis+"/FIELD") 00253 try: 00254 # get the name for later display 00255 phasecentername = tb.getcell('NAME', phasecenter) + " (original field " + str(phasecenter) + ")" 00256 tb.close() 00257 # phase center id will change if we select on fields, 00258 # the name column is only optionally filled and not necessarily unique. 00259 # But ms.msseltoindex(vis,field=field)['field'] returns the old field ids 00260 # in the order in which they will occur in the new field table. 00261 # => need to get index from there as new phase center ID 00262 newphasecenter = (ms.msseltoindex(vis,field=field)['field']).tolist().index(phasecenter) 00263 phasecentername += ", new field " + str(newphasecenter) + ")" 00264 except: 00265 tb.close() 00266 message = "phasecenter field id " + str(phasecenter) + " is not valid" 00267 ms.close() 00268 raise Exception, message 00269 00270 if(mode=='frequency'): 00271 ## reset the default values 00272 if(start==0): 00273 start = "" 00274 if(width==1): 00275 width = "" 00276 ##check that start and width have units if they are non-empty 00277 if not(start==""): 00278 if (qa.quantity(start)['unit'].find('Hz') < 0): 00279 ms.close() 00280 raise TypeError, "start parameter is not a valid frequency quantity " %start 00281 if(not(width=="") and (qa.quantity(width)['unit'].find('Hz') < 0)): 00282 ms.close() 00283 raise TypeError, "width parameter %s is not a valid frequency quantity " %width 00284 elif(mode=='velocity'): 00285 ## reset the default values 00286 if(start==0): 00287 start = "" 00288 if(width==1): 00289 width = "" 00290 ##check that start and width have units if they are non-empty 00291 if not(start==""): 00292 if (qa.quantity(start)['unit'].find('m/s') < 0): 00293 ms.close() 00294 raise TypeError, "start parameter %s is not a valid velocity quantity " %start 00295 if(not(width=="") and (qa.quantity(width)['unit'].find('m/s') < 0)): 00296 ms.close() 00297 raise TypeError, "width parameter %s is not a valid velocity quantity " %width 00298 elif(mode=='channel' or mode=='channel_b'): 00299 if((type(width) != int) or (type(start) != int)): 00300 ms.close() 00301 raise TypeError, "start and width have to be integers with mode = %s" %mode 00302 00303 ms.close() 00304 00305 # determine parameter "datacolumn" 00306 tb.open(vis) 00307 allcols = tb.colnames() 00308 tb.close() 00309 dpresent = ('DATA' in allcols) 00310 mpresent = ('MODEL_DATA' in allcols) 00311 cpresent = ('CORRECTED_DATA' in allcols) 00312 if (dpresent and cpresent): 00313 datacolumn = 'all' 00314 elif (dpresent and not cpresent): 00315 datacolumn = 'data' 00316 elif (cpresent and not dpresent): 00317 datacolumn = 'corrected_data' 00318 elif (mpresent and not cpresent and not dpresent): 00319 datacolumn = 'model_data' 00320 else: 00321 raise Exception, "Neither DATA nor CORRECTED_DATA nor MODEL_DATA column present. Cannot proceed." 00322 00323 if(doselection): 00324 casalog.post("Creating selected SubMS ...", 'INFO') 00325 ms.open(vis) 00326 ms.split(outputms=outputvis, field=field, 00327 spw=spw, step=[1], 00328 baseline=antenna, subarray=array, 00329 timebin='-1s', time=timerange, 00330 whichcol=datacolumn, 00331 scan=scan, uvrange="") 00332 ms.close() 00333 00334 else: 00335 # no selection necessary, just copy 00336 casalog.post("Creating working copy ...", 'INFO') 00337 shutil.rmtree(outputvis,ignore_errors=True) 00338 shutil.copytree(vis,outputvis) 00339 00340 # Combine and if necessary regrid it 00341 ms.open(outputvis, nomodify=False) 00342 00343 message = "Using " + phasecentername + " as phase center." 00344 casalog.post(message, 'INFO') 00345 00346 if not ms.cvel(mode=mode, nchan=nchan, start=start, width=width, 00347 interp=interpolation, 00348 phasec=newphasecenter, restfreq=restfreq, 00349 outframe=outframe, veltype=veltype, hanning=hanning): 00350 ms.close() 00351 raise Exception, "Error in regridding step ..." 00352 ms.close() 00353 00354 # deal with the passall option 00355 temp_suffix = ".deselected" 00356 if (passall): 00357 # determine range of fields 00358 fieldsel = ms.msseltoindex(vis=vis, field=field)['field'] 00359 tb.open(vis+"/FIELD") 00360 nfields = tb.nrows() 00361 tb.close() 00362 fielddesel = "" 00363 for i in xrange(0,nfields): 00364 if not (i in fieldsel): 00365 if not (fielddesel == ""): 00366 fielddesel += "," 00367 fielddesel += str(i) 00368 00369 # determine range of SPWs 00370 spwsel = ms.msseltoindex(vis=vis, spw=spw)['spw'] 00371 tb.open(vis+"/SPECTRAL_WINDOW") 00372 nspws = tb.nrows() 00373 tb.close() 00374 spwdesel = "" 00375 for i in xrange(0,nspws): 00376 if not (i in spwsel): 00377 if not (spwdesel == ""): 00378 spwdesel += "," 00379 spwdesel += str(i) 00380 00381 if not (fielddesel == "" and spwdesel == ""): 00382 # split out the data not selected by the conditions on field and spw 00383 # from the original MS and join it to the output MS 00384 00385 # need to do this in two steps 00386 # I) field == "*" and deselected spws 00387 if not (spwdesel == ""): 00388 ms.open(vis) 00389 casalog.post("Passing through data with", 'INFO') 00390 casalog.post(" spws: " + spwdesel, 'INFO') 00391 00392 ms.split(outputms=outputvis+temp_suffix, field='*', 00393 spw=spwdesel, step=[1], 00394 baseline=antenna, subarray=array, 00395 timebin='-1s', time=timerange, 00396 whichcol=datacolumn, 00397 scan=scan, uvrange="") 00398 ms.close() 00399 00400 # join with the deselected part 00401 ms.open(outputvis, nomodify=False) 00402 rval = ms.concatenate(msfile = outputvis+temp_suffix) 00403 ms.close() 00404 shutil.rmtree(outputvis+temp_suffix) 00405 if not rval: 00406 raise Exception, "Error in attaching passed-through data ..." 00407 00408 # II) deselected fields and selected spws 00409 if not (fielddesel == ""): 00410 ms.open(vis) 00411 casalog.post("Passing through data with", 'INFO') 00412 casalog.post(" fields: " + fielddesel, 'INFO') 00413 casalog.post(" spws: " + spw, 'INFO') 00414 00415 ms.split(outputms=outputvis+temp_suffix, field=fielddesel, 00416 spw=spw, step=[1], 00417 baseline=antenna, subarray=array, 00418 timebin='-1s', time=timerange, 00419 whichcol=datacolumn, 00420 scan=scan, uvrange="") 00421 ms.close() 00422 00423 # join with the deselected part 00424 ms.open(outputvis, nomodify=False) 00425 rval = ms.concatenate(msfile = outputvis+temp_suffix) 00426 ms.close() 00427 shutil.rmtree(outputvis+temp_suffix) 00428 if not rval: 00429 raise Exception, "Error in attaching passed-through data ..." 00430 00431 # Write history to output MS 00432 ms.open(outputvis, nomodify=False) 00433 ms.writehistory(message='taskname=cvel', origin='cvel') 00434 ms.writehistory(message='vis = "'+str(vis)+'"', 00435 origin='cvel') 00436 ms.writehistory(message='outputvis = "'+str(outputvis)+'"', 00437 origin='cvel') 00438 ms.writehistory(message='passall = "'+str(passall)+'"', 00439 origin='cvel') 00440 ms.writehistory(message='field = "'+str(field)+'"', 00441 origin='cvel') 00442 ms.writehistory(message='spw = '+str(spw), origin='cvel') 00443 ms.writehistory(message='antenna = "'+str(antenna)+'"', 00444 origin='cvel') 00445 ms.writehistory(message='timerange = "'+str(timerange)+'"', 00446 origin='cvel') 00447 ms.writehistory(message='mode = '+str(mode), origin='cvel') 00448 ms.writehistory(message='nchan = '+str(nchan), origin='cvel') 00449 ms.writehistory(message='start = '+str(start), origin='cvel') 00450 ms.writehistory(message='width = '+str(width), origin='cvel') 00451 ms.writehistory(message='interpolation = '+str(interpolation), origin='cvel') 00452 ms.writehistory(message='outframe = "'+str(outframe)+'"', 00453 origin='cvel') 00454 ms.writehistory(message='phasecenter = "'+ str(phasecentername) +'"', 00455 origin='cvel') 00456 ms.writehistory(message='hanning = "'+str(hanning)+'"', 00457 origin='cvel') 00458 ms.close() 00459 00460 return True 00461 00462 except Exception, instance: 00463 print '*** Error *** ',instance 00464 # delete temp output (comment out for debugging) 00465 if os.path.exists(outputvis+".spwCombined"): 00466 casalog.post("Deleting temporary output files ...", 'INFO') 00467 shutil.rmtree(outputvis+".spwCombined") 00468 raise Exception, instance 00469 00470
1.8.0
