Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_cvel.py: 3%
247 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 19:53 +0000
1import os
2import shutil
4from casatools import table, quanta
5from casatools import ms as mstool
6from casatasks import casalog
7from .mstools import write_history
9qa = quanta()
11def cvel(vis, outputvis,
12 passall, field, spw, selectdata, antenna, timerange, scan, array,
13 mode, nchan, start, width, interpolation,
14 phasecenter, restfreq, outframe, veltype, hanning):
16 """ regrid an MS to a new spectral window / channel structure or frame
18 vis -- Name of input visibility file
19 default: none; example: vis='ngc5921.ms'
21 outputvis -- Name of output measurement set (required)
22 default: none; example: vis='ngc5921-regridded.ms'
24 passall -- if False, data not meeting the selection on field or spw
25 is omitted; if True, data not meeting the selection
26 is passed through without modification
27 default: False; example:
28 field='NGC5921'
29 passall=False : only data from NGC5921 is included in output MS,
30 no data from other fields (e.g. 1331+305) is included
31 passall=True : data from NGC5921 is transformed by cvel, all other
32 fields are passed through unchanged
34 field -- Select fields in MS. Use field id(s) or field name(s).
35 ['go listobs' to obtain the list id's or names]
36 default: ''= all fields
37 If field string is a non-negative integer, it is assumed to
38 be a field index otherwise, it is assumed to be a
39 field name
40 field='0~2'; field ids 0,1,2
41 field='0,4,5~7'; field ids 0,4,5,6,7
42 field='3C286,3C295'; field named 3C286 and 3C295
43 field = '3,4C*'; field id 3, all names starting with 4C
45 spw --Select spectral window/channels
46 NOTE: This selects the data passed as the INPUT to mode
47 default: ''=all spectral windows and channels
48 spw='0~2,4'; spectral windows 0,1,2,4 (all channels)
49 spw='0:5~61'; spw 0, channels 5 to 61
50 spw='<2'; spectral windows less than 2 (i.e. 0,1)
51 spw='0,10,3:3~45'; spw 0,10 all channels, spw 3,
52 channels 3 to 45.
53 spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each.
54 spw='0:0~10;15~60'; spectral window 0 with channels
55 0-10,15-60
56 spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10,
57 spw 1, channels 20-30, and spw 2, channels, 1,2 and 3
59 selectdata -- Other data selection parameters
60 default: True
62 selectdata=True expandable parameters
64 antenna -- Select data based on antenna/baseline
65 default: '' (all)
66 If antenna string is a non-negative integer, it is
67 assumed to be an antenna index, otherwise, it is
68 considered an antenna name.
69 antenna='5&6'; baseline between antenna index 5 and
70 index 6.
71 antenna='VA05&VA06'; baseline between VLA antenna 5
72 and 6.
73 antenna='5&6;7&8'; baselines 5-6 and 7-8
74 antenna='5'; all baselines with antenna index 5
75 antenna='05'; all baselines with antenna number 05
76 (VLA old name)
77 antenna='5,6,9'; all baselines with antennas 5,6,9
78 index numbers
80 timerange -- Select data based on time range:
81 default = '' (all); examples,
82 timerange = 'YYYY/MM/DD/hh:mm:ss~YYYY/MM/DD/hh:mm:ss'
83 Note: if YYYY/MM/DD is missing date defaults to first
84 day in data set
85 timerange='09:14:0~09:54:0' picks 40 min on first day
86 timerange= '25:00:00~27:30:00' picks 1 hr to 3 hr
87 30min on NEXT day
88 timerange='09:44:00' pick data within one integration
89 of time
90 timerange='>10:24:00' data after this time
92 scan -- Scan number range.
93 default: '' (all)
94 example: scan='1~5'
95 Check 'go listobs' to insure the scan numbers are in
96 order.
98 array -- Select data by (sub)array indices
99 default: '' (all); example:
100 array='0~2'; arrays 0 to 2
102 mode -- Frequency Specification:
103 NOTE: See examples below:
104 default: 'channel'
105 mode = 'channel'; Use with nchan, start, width to specify
106 output spw. See examples below
107 mode = 'velocity', means channels are specified in
108 velocity.
109 mode = 'frequency', means channels are specified in
110 frequency.
112 mode expandable parameters
113 Start, width are given in units of channels, frequency
114 or velocity as indicated by mode
115 nchan -- Number of channels in output spw
116 default: -1 = all channels; example: nchan=3
117 start -- Start input channel (relative-0)
118 default=0; example: start=5
119 width -- Output channel width in units of the input
120 channel width (>1 indicates channel averaging)
121 default=1; example: width=4
122 interpolation -- Interpolation method
123 default = 'linear'
125 examples:
126 spw = '0,1'; mode = 'channel'
127 will produce a single spw containing all channels in spw
128 0 and 1
129 spw='0:5~28^2'; mode = 'channel'
130 will produce a single spw made with channels
131 (5,7,9,...,25,27)
132 spw = '0'; mode = 'channel': nchan=3; start=5; width=4
133 will produce an spw with 3 output channels
134 new channel 1 contains data from channels (5+6+7+8)
135 new channel 2 contains data from channels (9+10+11+12)
136 new channel 3 contains data from channels (13+14+15+16)
137 spw = '0:0~63^3'; mode='channel'; nchan=21; start = 0;
138 width = 1
139 will produce an spw with 21 channels
140 new channel 1 contains data from channel 0
141 new channel 2 contains data from channel 2
142 new channel 21 contains data from channel 61
143 spw = '0:0~40^2'; mode = 'channel'; nchan = 3; start =
144 5; width = 4
145 will produce an spw with three output channels
146 new channel 1 contains channels (5,7)
147 new channel 2 contains channels (13,15)
148 new channel 3 contains channels (21,23)
150 phasecenter -- direction measure or fieldid for the mosaic center
151 default: '' => first field selected ; example: phasecenter=6
152 or phasecenter='J2000 19h30m00 -40d00m00'
154 restfreq -- Specify rest frequency to use for output image
155 default='' Occasionally it is necessary to set this (for
156 example some VLA spectral line data). For example for
157 NH_3 (1,1) put restfreq='23.694496GHz'
159 outframe -- output reference frame (not case-sensitive)
160 possible values: LSRK, LSRD, BARY, GALACTO, LGROUP, CMB, GEO, TOPO or SOURCE
161 default='' (keep original reference frame) ; example: outframe='BARY'
163 veltype -- definition of velocity (in mode)
164 default = 'radio'
166 hanning -- if true, Hanning smooth frequency channel data before regridding
167 to remove Gibbs ringing
168 default = False
170 """
172 # Note: this is duplicated in task_cvel, and really needing CASA-wide harmonization
173 # (CAS-12871)
174 def copy_ms(src, dest):
175 """ This is a MMS-safe copy of an MMS tree directory.
176 :param src: path to the source MS
177 :param dest: path to the destination MS
178 """
179 shutil.copytree(src, dest, symlinks=True)
181 #Python script
183 # make ms, table tools local
184 _ms = mstool()
185 _tb = table()
187 try:
188 casalog.origin('cvel')
190 if not ((type(vis)==str) & (os.path.exists(vis))):
191 raise ValueError('Visibility data set not found - please verify the name')
193 assert outputvis != '', "Must provide output data set name in parameter outputvis."
194 assert not os.path.exists(outputvis), "Output MS %s already exists - will not overwrite." % outputvis
195 assert not os.path.exists(outputvis+".flagversions"), \
196 "The flagversions \"%s.flagversions\" for the output MS already exist. Please delete." % outputvis
198 # Handle selectdata explicitly
199 # (avoid hidden globals)
200 if (selectdata==False):
201 timerange=''
202 array=''
203 antenna=''
204 scan=''
206 if(type(antenna) == list):
207 antenna = ', '.join([str(ant) for ant in antenna])
209 if (field == ''):
210 field = '*'
212 if (spw == ''):
213 spw = '*'
215 if(passall and spw=='*' and field=='*'):
216 # all spws and fields selected, nothing to pass through
217 passall = False
219 doselection = True
220 if(field=='*' and spw=='*' and
221 (not selectdata or (selectdata and antenna=='' and timerange=='' and scan=='' and array==''))
222 ):
223 doselection = False
225 # open input MS
226 _ms.open(vis)
228 if(type(phasecenter)==str):
229 ### blank means take field 0
230 if (phasecenter==''):
231 phasecenter=_ms.msseltoindex(vis,field=field)['field'][0]
232 else:
233 tmppc=phasecenter
234 try:
235 if(len(_ms.msseltoindex(vis, field=phasecenter)['field']) > 0):
236 tmppc=_ms.msseltoindex(vis, field=phasecenter)['field'][0]
237 ##succesful must be string like '0' or 'NGC*'
238 except Exception:
239 ##failed must be a string 'J2000 18h00m00 10d00m00'
240 tmppc=phasecenter
241 phasecenter=tmppc
243 newphasecenter = phasecenter
244 phasecentername = phasecenter
245 if not (type(phasecenter)==str):
246 # first check that this field will be in the output MS
247 if not (phasecenter in _ms.msseltoindex(vis,field=field)['field']):
248 message = "Field id " + str(phasecenter)
249 message += " was selected as phasecenter but is not among the fields selected for output: "
250 message += str(_ms.msseltoindex(vis,field=field)['field'])
251 _ms.close()
252 raise RuntimeError(message)
254 _tb.open(vis+"/FIELD")
255 try:
256 # get the name for later display
257 phasecentername = _tb.getcell('NAME', phasecenter) + " (original field " + str(phasecenter)
258 _tb.close()
259 # phase center id will change if we select on fields,
260 # the name column is only optionally filled and not necessarily unique.
261 # But _ms.msseltoindex(vis,field=field)['field'] returns the old field ids
262 # in the order in which they will occur in the new field table.
263 # => need to get index from there as new phase center ID
264 newphasecenter = (_ms.msseltoindex(vis,field=field)['field']).tolist().index(phasecenter)
265 phasecentername += ", new field " + str(newphasecenter) + ")"
266 except:
267 _tb.close()
268 message = "phasecenter field id " + str(phasecenter) + " is not valid"
269 _ms.close()
270 raise RuntimeError(message)
272 if(mode=='frequency'):
273 ## reset the default values
274 if(start==0):
275 start = ""
276 if(width==1):
277 width = ""
278 ##check that start and width have units if they are non-empty
279 if not(start==""):
280 if (qa.quantity(start)['unit'].find('Hz') < 0):
281 _ms.close()
282 raise TypeError("start parameter is not a valid frequency quantity " %start)
283 if(not(width=="") and (qa.quantity(width)['unit'].find('Hz') < 0)):
284 _ms.close()
285 raise TypeError("width parameter %s is not a valid frequency quantity " %width)
287 elif(mode=='velocity'):
288 ## reset the default values
289 if(start==0):
290 start = ""
291 if(width==1):
292 width = ""
293 ##check that start and width have units if they are non-empty
294 if not(start==""):
295 if (qa.quantity(start)['unit'].find('m/s') < 0):
296 _ms.close()
297 raise TypeError("start parameter %s is not a valid velocity quantity " %start)
298 if(not(width=="") and (qa.quantity(width)['unit'].find('m/s') < 0)):
299 _ms.close()
300 raise TypeError("width parameter %s is not a valid velocity quantity " %width)
302 elif(mode=='channel' or mode=='channel_b'):
303 if((type(width) != int) or (type(start) != int)):
304 _ms.close()
305 raise TypeError("start and width have to be integers with mode = %s" %mode)
308 ## check if preaveraging is necessary
309 dopreaverage=False
310 preavwidth = [1]
312 thespwsel = _ms.msseltoindex(vis=vis, spw=spw)['spw']
313 thefieldsel = _ms.msseltoindex(vis=vis, field=field)['field']
314 outgrid = _ms.cvelfreqs(spwids=thespwsel, fieldids=thefieldsel,
315 mode=mode, nchan=nchan, start=start, width=width,
316 phasec=newphasecenter, restfreq=restfreq,
317 outframe=outframe, veltype=veltype, verbose=False)
318 if(len(outgrid)>1):
319 tmpavwidth = []
320 for thespw in thespwsel:
321 outgridw1 = _ms.cvelfreqs(spwids=[thespw], fieldids=thefieldsel,
322 mode='channel', nchan=-1, start=0, width=1, # native width
323 phasec=newphasecenter, restfreq=restfreq,
324 outframe=outframe, veltype=veltype, verbose=False)
325 if(len(outgridw1)>1):
326 widthratio = abs((outgrid[1]-outgrid[0])/(outgridw1[1]-outgridw1[0]))
327 if(widthratio>=2.0): # do preaverage
328 tmpavwidth.append( int(widthratio+0.001) )
329 dopreaverage=True
330 else:
331 tmpavwidth.append(1)
333 if dopreaverage:
334 preavwidth = tmpavwidth
336 _ms.close()
338 # if in channel mode and preaveraging,
339 if(dopreaverage and (mode=='channel' or mode=='channel_b')):
340 if(max(preavwidth)==1):
341 dopreaverage = False
342 else:
343 # switch to frequency mode
344 mode = 'frequency'
345 if(width>0):
346 start = str(outgrid[0]/1E6)+'MHz'
347 width = str((outgrid[1]-outgrid[0]))+'Hz'
348 else:
349 start = str(outgrid[len(outgrid)-1]/1E6)+'MHz'
350 width = str(-(outgrid[1]-outgrid[0]))+'Hz'
351 casalog.post("After pre-averaging, will switch to frequency mode with start="+start+", width = "+width, 'INFO')
353 if(dopreaverage and hanning and max(preavwidth)>2):
354 casalog.post("NOTE: You have requested Hanning smoothing and at the same time you have chosen\n"
355 +"a large width parameter which makes pre-averaging necessary.\n"
356 +"The Hanning-smoothing may be redundant in this context.\n", 'WARN')
358 if dopreaverage:
359 # Past this point we know we are going to 'dopreaverage'
360 # CAS-9798
361 raise RuntimeError('ERROR: cvel does not regrid properly for channel '
362 'widths > or = 2x the native channel width, please use '
363 'mstransform, clean, or tclean for larger regridding. '
364 'All earlier versions of CASA also have this issue.')
366 # determine parameter "datacolumn"
367 _tb.open(vis)
368 allcols = _tb.colnames()
369 _tb.close()
370 dpresent = ('DATA' in allcols)
371 mpresent = ('MODEL_DATA' in allcols)
372 cpresent = ('CORRECTED_DATA' in allcols)
373 if (dpresent and cpresent):
374 datacolumn = 'all'
375 elif (dpresent and not cpresent):
376 datacolumn = 'data'
377 elif (cpresent and not dpresent):
378 datacolumn = 'corrected_data'
379 elif (mpresent and not cpresent and not dpresent):
380 datacolumn = 'model_data'
381 else:
382 raise RuntimeError("Neither DATA nor CORRECTED_DATA nor MODEL_DATA column present. Cannot proceed.")
384 if(doselection and not dopreaverage):
385 casalog.post("Creating selected SubMS ...", 'INFO')
386 _ms.open(vis)
387 _ms.split(outputms=outputvis, field=field,
388 spw=spw, step=[1],
389 baseline=antenna, subarray=array,
390 timebin='-1s', time=timerange,
391 whichcol=datacolumn,
392 scan=scan, uvrange="")
393 _ms.close()
395 elif(dopreaverage and not doselection):
396 if(hanning):
397 casalog.post("Creating working copy for Hanning-smoothing ...", 'INFO')
398 shutil.rmtree(outputvis+'TMP',ignore_errors=True)
399 copy_ms(vis, outputvis+'TMP')
400 _ms.open(outputvis+'TMP', nomodify=False)
401 _ms.hanningsmooth(datacolumn=datacolumn)
402 _ms.close()
403 hanning = False
404 _ms.open(outputvis+'TMP')
405 else:
406 _ms.open(vis)
408 casalog.post("Creating preaveraged SubMS using widths "+str(preavwidth), 'INFO')
409 _ms.split(outputms=outputvis,
410 whichcol=datacolumn, step=preavwidth)
411 _ms.close()
413 elif(doselection and dopreaverage):
414 if(hanning):
415 casalog.post("Creating selected working copy for Hanning-smoothing ...", 'INFO')
416 shutil.rmtree(outputvis+'TMP',ignore_errors=True)
417 _ms.open(vis)
418 _ms.split(outputms=outputvis+'TMP', field=field,
419 spw=spw, step=[1],
420 baseline=antenna, subarray=array,
421 timebin='-1s', time=timerange,
422 whichcol=datacolumn,
423 scan=scan, uvrange="")
424 _ms.close()
425 _ms.open(outputvis+'TMP', nomodify=False)
426 _ms.hanningsmooth(datacolumn=datacolumn)
427 _ms.close()
428 hanning = False
429 _ms.open(outputvis+'TMP')
430 casalog.post("Creating preaveraged SubMS using widths "+str(preavwidth), 'INFO')
431 _ms.split(outputms=outputvis,
432 whichcol=datacolumn, step=preavwidth)
433 _ms.close()
434 else:
435 casalog.post("Creating selected, preaveraged SubMS using widths "+str(preavwidth), 'INFO')
436 _ms.open(vis)
437 _ms.split(outputms=outputvis, field=field,
438 spw=spw, step=preavwidth,
439 baseline=antenna, subarray=array,
440 timebin='-1s', time=timerange,
441 whichcol=datacolumn,
442 scan=scan, uvrange="")
443 _ms.close()
445 else:
446 # no selection or preaveraging necessary, just copy
447 casalog.post("Creating working copy ...", 'INFO')
448 shutil.rmtree(outputvis,ignore_errors=True)
449 copy_ms(vis, outputvis)
451 # Combine and if necessary regrid it
452 _ms.open(outputvis, nomodify=False)
454 message = "Using " + phasecentername + " as common direction for the output reference frame."
455 casalog.post(message, 'INFO')
457 if not _ms.cvel(mode=mode, nchan=nchan, start=start, width=width,
458 interp=interpolation,
459 phasec=newphasecenter, restfreq=restfreq,
460 outframe=outframe, veltype=veltype, hanning=hanning):
461 _ms.close()
462 raise RuntimeError("Error in regridding step ...")
463 _ms.close()
465 # deal with the passall option
466 temp_suffix = ".deselected"
467 if (passall):
468 # determine range of fields
469 fieldsel = _ms.msseltoindex(vis=vis, field=field)['field']
470 _tb.open(vis+"/FIELD")
471 nfields = _tb.nrows()
472 _tb.close()
473 fielddesel = ""
474 for i in range(0,nfields):
475 if not (i in fieldsel):
476 if not (fielddesel == ""):
477 fielddesel += ","
478 fielddesel += str(i)
480 # determine range of SPWs
481 spwsel = _ms.msseltoindex(vis=vis, spw=spw)['spw']
482 _tb.open(vis+"/SPECTRAL_WINDOW")
483 nspws = _tb.nrows()
484 _tb.close()
485 spwdesel = ""
486 for i in range(0,nspws):
487 if not (i in spwsel):
488 if not (spwdesel == ""):
489 spwdesel += ","
490 spwdesel += str(i)
492 if not (fielddesel == "" and spwdesel == ""):
493 # split out the data not selected by the conditions on field and spw
494 # from the original MS and join it to the output MS
496 # need to do this in two steps
497 # I) field == "*" and deselected spws
498 if not (spwdesel == ""):
499 _ms.open(vis)
500 casalog.post("Passing through data with", 'INFO')
501 casalog.post(" spws: " + spwdesel, 'INFO')
503 _ms.split(outputms=outputvis+temp_suffix, field='*',
504 spw=spwdesel, step=[1],
505 baseline=antenna, subarray=array,
506 timebin='-1s', time=timerange,
507 whichcol=datacolumn,
508 scan=scan, uvrange="")
509 _ms.close()
511 # join with the deselected part
512 _ms.open(outputvis, nomodify=False)
513 rval = _ms.concatenate(msfile = outputvis+temp_suffix)
514 _ms.close()
515 shutil.rmtree(outputvis+temp_suffix)
516 if not rval:
517 raise RuntimeError("Error in attaching passed-through data ...")
519 # II) deselected fields and selected spws
520 if not (fielddesel == ""):
521 _ms.open(vis)
522 casalog.post("Passing through data with", 'INFO')
523 casalog.post(" fields: " + fielddesel, 'INFO')
524 casalog.post(" spws: " + spw, 'INFO')
526 _ms.split(outputms=outputvis+temp_suffix, field=fielddesel,
527 spw=spw, step=[1],
528 baseline=antenna, subarray=array,
529 timebin='-1s', time=timerange,
530 whichcol=datacolumn,
531 scan=scan, uvrange="")
532 _ms.close()
534 # join with the deselected part
535 _ms.open(outputvis, nomodify=False)
536 rval = _ms.concatenate(msfile = outputvis+temp_suffix)
537 _ms.close()
538 shutil.rmtree(outputvis+temp_suffix)
539 if not rval:
540 raise RuntimeError("Error in attaching passed-through data ...")
542 # Write history to output MS
543 try:
544 param_names = cvel.__code__.co_varnames[:cvel.__code__.co_argcount]
545 local_vars = locals()
546 param_vals = [local_vars[p] for p in param_names]
548 write_history(mstool(), outputvis, 'cvel', param_names,
549 param_vals, casalog)
550 except Exception as instance:
551 casalog.post("*** Error \'%s\' updating HISTORY" % (instance),
552 'WARN')
554 finally:
555 # delete temp output (comment out for debugging)
556 if os.path.exists(outputvis+".spwCombined"):
557 casalog.post("Deleting temporary output files ...", 'INFO')
558 shutil.rmtree(outputvis+".spwCombined")