Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_tsdimaging.py: 88%
637 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1# sd task for imaging
3import collections
4import contextlib
5import os
6import re
7import shutil
8import time
10import numpy
12from functools import partial
14from casatasks import casalog
15from casatools import image, imager
16from casatools import ms as mstool
17from casatools import quanta
19from . import mslisthelper, sdbeamutil, sdutil
20# (1) Import the python application layer
21from .imagerhelpers.imager_base import PySynthesisImager
22from .imagerhelpers.input_parameters import ImagerParameters
24image_suffix = '.image'
25residual_suffix = '.residual'
26weight_suffix = '.weight'
27associate_suffixes = ['.psf', '.sumwt', weight_suffix, residual_suffix]
30@contextlib.contextmanager
31def open_ia(imagename):
32 ia = image()
33 ia.open(imagename)
34 try:
35 yield ia
36 finally:
37 ia.close()
40@contextlib.contextmanager
41def open_ms(vis):
42 ms = mstool()
43 ms.open(vis)
44 try:
45 yield ms
46 finally:
47 ms.close()
50class SelectionHandler(object):
51 def __init__(self, sel):
52 self.sel = sel
53 if isinstance(self.sel, str):
54 self.selector = self._select0
55 elif len(self.sel) == 1:
56 self.selector = self._select1
57 else:
58 self.selector = self._select2
60 def __call__(self, i):
61 return self.selector(i)
63 def _select0(self, i):
64 return self.sel
66 def _select1(self, i):
67 return self.sel[0]
69 def _select2(self, i):
70 return self.sel[i]
73class OldImagerBasedTools(object):
74 def __init__(self):
75 self.imager = imager()
77 @contextlib.contextmanager
78 def open_old_imager(self, vis):
79 try:
80 self.imager.open(vis)
81 yield self.imager
82 finally:
83 self.imager.close()
85 @contextlib.contextmanager
86 def open_and_select_old_imager(self, vislist, field, spw, antenna, scan, intent, timerange):
87 if isinstance(vislist, str):
88 with self.open_old_imager(vislist) as im:
89 im.selectvis(field=field,
90 spw=spw,
91 nchan=-1,
92 start=0,
93 step=1,
94 baseline=antenna,
95 scan=scan,
96 intent=intent,
97 time=timerange)
98 yield im
99 else:
100 fieldsel = SelectionHandler(field)
101 spwsel = SelectionHandler(spw)
102 antennasel = SelectionHandler(antenna)
103 scansel = SelectionHandler(scan)
104 intentsel = SelectionHandler(intent)
105 timerangesel = SelectionHandler(timerange)
106 try:
107 for i in range(len(vislist)):
108 vis = vislist[i]
109 _field = fieldsel(i)
110 _spw = spwsel(i)
111 _antenna = antennasel(i)
112 _scan = scansel(i)
113 _intent = intentsel(i)
114 _timerangesel = timerangesel(i)
115 if len(_antenna) == 0:
116 _baseline = _antenna
117 elif len(_antenna) < 4 or _antenna[:-3] != '&&&':
118 _baseline = _antenna + '&&&'
119 else:
120 _baseline = _antenna
121 self.imager.selectvis(
122 vis, field=_field, spw=_spw, nchan=-1, start=0, step=1,
123 baseline=_baseline, scan=_scan, intent=_intent, time=_timerangesel)
124 yield self.imager
125 finally:
126 self.imager.close()
128 def test(self, vis):
129 with self.open_old_imager(vis):
130 casalog.post('test')
131 raise RuntimeError('ERROR!')
133 def get_pointing_sampling_params(self, vis, field, spw, baseline, scan, intent, timerange,
134 outref, movingsource, pointingcolumntouse, antenna_name):
135 with self.open_old_imager(vis) as im:
136 im.selectvis(field=field,
137 spw=spw,
138 nchan=-1,
139 start=0,
140 step=1,
141 baseline=baseline,
142 scan=scan,
143 intent=intent,
144 time=timerange)
145 sampling_params = im.pointingsampling(pattern='raster',
146 ref=outref,
147 movingsource=movingsource,
148 pointingcolumntouse=pointingcolumntouse,
149 antenna='{0}&&&'.format(antenna_name))
150 return sampling_params
152 def get_map_extent(self, vislist, field, spw, antenna, scan, intent, timerange,
153 ref, movingsource, pointingcolumntouse):
155 with self.open_and_select_old_imager(vislist=vislist, field=field,
156 spw=spw, antenna=antenna, scan=scan,
157 intent=intent, timerange=timerange) as im:
158 map_param = im.mapextent(ref=ref, movingsource=movingsource,
159 pointingcolumntouse=pointingcolumntouse)
160 return map_param
163def check_conformance(mslist, check_result):
164 """Check conformance of input MS list
166 Check conformance of input MS list, particularlly existence of
167 WEIGHT_SPECTRUM column.
169 Args:
170 mslist (list): list of names for input MS
171 check_result (dict): result of conformance check.
172 see mslisthelper.check_mslist for detail about
173 the structure of check_result.
175 Returns:
176 dict: Per-column set of names for MS that needs to be
177 edited to resolve the conformance. Top level dict
178 has two keys, "remove" and "add", which indicate
179 the operation to be applied to the columns.
180 """
181 process_dict = {
182 'remove': collections.defaultdict(set),
183 'add': collections.defaultdict(set)
184 }
185 for name, summary in check_result.items():
186 if 'Main' in summary:
187 missingcol_list = [summary['Main']['missingcol_{}'.format(x)] for x in ['a', 'b']]
188 ms_a = mslist[0]
189 ms_b = name
191 # "remove" operation:
192 # - MS list is opposite order to missingcol_list
193 # - if WEIGHT_SPECTRUM column is missing in one MS,
194 # the column should be removed from another MS.
195 for c, t in zip(missingcol_list, [ms_b, ms_a]):
196 if 'WEIGHT_SPECTRUM' in c:
197 process_dict['remove']['WEIGHT_SPECTRUM'].add(t)
198 # "add" operation:
199 # - MS list is same order as missingcol_list
200 # - if CORRECTED_DATA column is missing in one MS,
201 # the column should be added to that MS.
202 for c, t in zip(missingcol_list, [ms_a, ms_b]):
203 if 'CORRECTED_DATA' in c:
204 process_dict['add']['CORRECTED_DATA'].add(t)
205 return process_dict
208def report_conformance(mslist, column_name, process_set):
209 """Report conformance of input MS
211 Report conformance of input MS, particularlly on the existence
212 of the column given by column_name.
214 Args:
215 mslist (list): list of names for input MS
216 column_name (str): name of the column
217 process_set (set): set of names of MS that need to be edited
218 """
219 if len(process_set) > 0:
220 casalog.post('Detected non-conformance of {} column in input list of MSes.'.format(column_name), priority='WARN')
221 cols = ['exists?', 'MS name']
222 header = ' '.join(cols)
223 casalog.post('', priority='WARN')
224 casalog.post('Summary of existence of {}:'.format(column_name), priority='WARN')
225 casalog.post(header, priority='WARN')
226 casalog.post('-' * len(header), priority='WARN')
227 for name in mslist:
228 basename = os.path.basename(name.rstrip('/'))
229 exists = 'YES' if name in process_set else 'NO'
230 row = '{:^7s} {:<s}'.format(exists, basename)
231 casalog.post(row, priority='WARN')
234def fix_conformance(process_dict):
235 """Resolve non-conformance by removing WEIGHT_SPECTRUM
237 Two non-conformances are fixed, WEIGHT_SPECTRUM and
238 CORRECTED_DATA. Remove WEIGHT_SPECTRUM column from, or
239 add CORRECTED_DATA to the MS provided by process_dict.
240 Backup is created with the name:
242 <original_name>.sdimaging.backup-<timestamp>
244 Args:
245 process_dict (dict): per-operation ("remove" and "add")
246 key-value pair of column name and
247 list of names for MS to be edited
249 Returns:
250 dict: mapping of original MS name and the name of backup
251 """
252 backup_list = {}
253 process_list = set()
254 for v in process_dict.values():
255 for w in v.values():
256 process_list = process_list.union(w)
257 for name in process_list:
258 basename = os.path.basename(name.rstrip('/'))
259 timestamp = time.strftime('%Y%m%dT%H%M%S', time.gmtime())
260 backup_name = basename + '.sdimaging.backup-{}'.format(timestamp)
261 with sdutil.table_manager(name) as tb:
262 tb.copy(backup_name, deep=True, returnobject=True).close()
263 backup_list[name] = backup_name
264 casalog.post('Copy of "{}" has been saved to "{}"'.format(name, backup_name), priority='WARN')
266 for colname, msnames in process_dict['remove'].items():
267 for name in msnames:
268 casalog.post('{} will be removed from "{}"'.format(colname, name), priority='WARN')
269 with sdutil.table_manager(name, nomodify=False) as tb:
270 if colname in tb.colnames():
271 tb.removecols(colname)
273 for colname, msnames in process_dict['add'].items():
274 for name in msnames:
275 casalog.post('{} will be added to "{}"'.format(colname, name), priority='WARN')
276 with sdutil.calibrater_manager(name, addmodel=False, addcorr=True):
277 pass
278 return backup_list
281def conform_mslist(mslist, ignore_columns=['CORRECTED_DATA']):
282 """Make given set of MS data conform
284 Here, only conformance on the existence of WEIGHT_SPECTRUM
285 is checked and resolved because non-conformance of WEIGHT_SPECTRUM,
286 i.e. some MS have the column while others don't, could cause
287 the task to crash. If non-conformance is detected, all existing
288 WEIGHT_SPECTRUM columns are removed. This opration modifies input
289 MS so data will be backed up with the name:
291 <original_name>.sdimaging.backup-<timestamp>
293 Args:
294 mslist (list): list of names for input MS
295 """
296 check_result = mslisthelper.check_mslist(mslist, testcontent=False)
297 process_dict = check_conformance(mslist, check_result)
298 fix_dict = dict()
299 for op, op_dict in process_dict.items():
300 fix_dict[op] = dict()
301 for col, process_set in op_dict.items():
302 if col not in ignore_columns:
303 report_conformance(mslist, col, process_set)
304 fix_dict[op][col] = process_set
305 fix_conformance(fix_dict)
308def sort_vis(vislist, spw, mode, width, field, antenna, scan, intent, timerange):
309 """Sort the given MeasurementSet path(s) by their earliest data-taking time, in increasing order.
311 Return a 7-tuple where
312 * the first entry holds the re-ordered paths
313 * the remaining entries hold their corresponding data selection parameters
314 * FIXME: input parameters mode and width are not used
315 """
316 if isinstance(vislist, str) or len(vislist) == 1:
317 return vislist, field, spw, antenna, scan, intent, timerange
318 # chronological sort
319 sorted_vislist, sorted_timelist = mslisthelper.sort_mslist(vislist)
320 _vislist = list(vislist)
321 sorted_idx = [_vislist.index(vis) for vis in sorted_vislist]
322 mslisthelper.report_sort_result(sorted_vislist, sorted_timelist, sorted_idx, mycasalog=casalog)
323 # conform MS
324 conform_mslist(sorted_vislist)
325 fieldsel = SelectionHandler(field)
326 sorted_field = [fieldsel(i) for i in sorted_idx]
327 spwsel = SelectionHandler(spw)
328 sorted_spw = [spwsel(i) for i in sorted_idx]
329 antennasel = SelectionHandler(antenna)
330 sorted_antenna = [antennasel(i) for i in sorted_idx]
331 scansel = SelectionHandler(scan)
332 sorted_scan = [scansel(i) for i in sorted_idx]
333 intentsel = SelectionHandler(intent)
334 sorted_intent = [intentsel(i) for i in sorted_idx]
335 timerangesel = SelectionHandler(timerange)
336 sorted_timerange = [timerangesel(i) for i in sorted_idx]
337 return sorted_vislist, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent, sorted_timerange
340def _configure_spectral_axis(mode, nchan, start, width, restfreq):
341 # fix default
342 if mode == 'channel':
343 if start == '':
344 start = 0
345 if width == '':
346 width = 1
347 else:
348 if start == 0:
349 start = ''
350 if width == 1:
351 width = ''
352 # fix unit
353 if mode == 'frequency':
354 myunit = 'Hz'
355 elif mode == 'velocity':
356 myunit = 'km/s'
357 else: # channel
358 myunit = ''
360 tmp_start = _format_quantum_unit(start, myunit)
361 if tmp_start is None:
362 raise ValueError("Invalid unit for %s in mode %s: %s" % ('start', mode, start))
363 start = tmp_start
364 if mode == 'channel':
365 start = int(start)
366 tmp_width = _format_quantum_unit(width, myunit)
367 if tmp_width is None:
368 raise ValueError("Invalid unit for %s in mode %s: %s" % ('width', mode, width))
369 width = tmp_width
370 if mode == 'channel':
371 width = int(width)
373 # TODO: work for nchan
374 imnchan = nchan
375 imstart = start
376 imwidth = width
377 return imnchan, imstart, imwidth
380def _format_quantum_unit(data, unit):
381 """Format quantity data.
383 Returns False if data has an unit which in not a variation of
384 input unit.
385 Otherwise, returns input data as a quantum string. The input
386 unit is added to the return value if no unit is in data.
387 """
388 my_qa = quanta()
389 if data == '' or my_qa.compare(data, unit):
390 return data
391 if my_qa.getunit(data) == '':
392 casalog.post("No unit specified. Using '%s'" % unit)
393 return '%f%s' % (data, unit)
394 return None
397def _handle_grid_defaults(value):
398 ret = ''
399 if isinstance(value, int) or isinstance(value, float):
400 ret = str(value)
401 elif isinstance(value, str):
402 ret = value
403 return ret
406def _calc_PB(vis, antenna_id, restfreq):
407 """Calculate the primary beam size of antenna.
409 Calculate the primary beam size of antenna, using dish diamenter
410 and rest frequency
411 Average antenna diamter and reference frequency are adopted for
412 calculation.
413 The input argument should be a list of antenna IDs.
414 """
415 logger = sdutil.Casalog(origin="_calc_PB")
416 logger.post("Calculating Primary beam size:")
417 # CAS-5410 Use private tools inside task scripts
418 my_qa = quanta()
420 pb_factor = 1.175
421 # Reference frequency
422 ref_freq = restfreq
423 if type(ref_freq) in [float, numpy.float64]:
424 ref_freq = my_qa.tos(my_qa.quantity(ref_freq, 'Hz'))
425 if not my_qa.compare(ref_freq, 'Hz'):
426 msg = "Could not get the reference frequency. " + \
427 "Your data does not seem to have valid one in selected field.\n" + \
428 "PB is not calculated.\n" + \
429 "Please set restreq or cell manually to generate an image."
430 raise RuntimeError(msg)
431 # Antenna diameter
432 with sdutil.table_manager(os.path.join(vis, 'ANTENNA')) as tb:
433 antdiam_ave = tb.getcell('DISH_DIAMETER', antenna_id)
434 # antdiam_ave = self._get_average_antenna_diameter(antenna)
435 # Calculate PB
436 wave_length = 0.2997924 / my_qa.convert(my_qa.quantity(ref_freq), 'GHz')['value']
437 D_m = my_qa.convert(antdiam_ave, 'm')['value']
438 lambda_D = wave_length / D_m * 3600. * 180 / numpy.pi
439 PB = my_qa.quantity(pb_factor * lambda_D, 'arcsec')
440 # Summary
441 logger.post(f"- Antenna diameter: {D_m} m")
442 logger.post(f"- Reference Frequency: {ref_freq}")
443 logger.post(f"PB size = {pb_factor:5.3f} * lambda/D = {my_qa.tos(PB)}")
444 return PB
447def _get_imsize(width, height, dx, dy):
448 casalog.post("Calculating pixel size.")
449 # CAS-5410 Use private tools inside task scripts
450 my_qa = quanta()
451 ny = numpy.ceil((my_qa.convert(height, my_qa.getunit(dy))['value'] /
452 my_qa.getvalue(dy)))
453 nx = numpy.ceil((my_qa.convert(width, my_qa.getunit(dx))['value'] /
454 my_qa.getvalue(dx)))
455 casalog.post("- Map extent: [%s, %s]" % (my_qa.tos(width), my_qa.tos(height)))
456 casalog.post("- Cell size: [%s, %s]" % (my_qa.tos(dx), my_qa.tos(dy)))
457 casalog.post("Image pixel numbers to cover the extent: [%d, %d] (projected)" %
458 (nx + 1, ny + 1))
459 return [int(nx + 1), int(ny + 1)]
462def _get_pointing_extent(phasecenter, vislist, field, spw, antenna, scan, intent, timerange,
463 pointingcolumntouse, ephemsrcname):
464 # MS selection is ignored. This is not quite right.
465 casalog.post("Calculating map extent from pointings.")
466 # CAS-5410 Use private tools inside task scripts
467 my_qa = quanta()
468 ret_dict = {}
470 if isinstance(vislist, str):
471 vis = vislist
472 else:
473 vis = vislist[0]
475 # colname = pointingcolumntouse.upper()
477 if phasecenter == "":
478 # defaut is J2000
479 base_mref = 'J2000'
480 elif isinstance(phasecenter, int) or phasecenter.isdigit():
481 # may be field id
482 with sdutil.table_manager(os.path.join(vis, 'FIELD')) as tb:
483 base_mref = tb.getcolkeyword('PHASE_DIR', 'MEASINFO')['Ref']
484 else:
485 # may be phasecenter is explicitly specified
486 # numeric value: 3.14, -.3e1, etc.
487 numeric_pattern = r'[-+]?([0-9]+(.[0-9]*)?|\.[0-9]+)([eE]-?[0-9])?'
488 # HMS string: 9:15:29, -9h15m29
489 hms_pattern = r'[-+]?[0-9]+[:h][0-9]+[:m][0-9.]+s?'
490 # DMS string: 9.15.29, -9d15m29s
491 dms_pattern = r'[-+]?[0-9]+[.d][0-9]+[.m][0-9.]+s?'
492 # composite pattern
493 pattern = fr'^({numeric_pattern}|{hms_pattern}|{dms_pattern})$'
494 items = phasecenter.split()
495 base_mref = 'J2000'
496 for i in items:
497 s = i.strip()
498 if re.match(pattern, s) is None:
499 base_mref = s
500 break
502 t = OldImagerBasedTools()
503 mapextent = t.get_map_extent(vislist, field, spw, antenna, scan, intent, timerange,
504 ref=base_mref, movingsource=ephemsrcname,
505 pointingcolumntouse=pointingcolumntouse)
506 # mapextent = self.imager.mapextent(ref=base_mref, movingsource=ephemsrcname,
507 # pointingcolumntouse=colname)
508 if mapextent['status']:
509 qheight = my_qa.quantity(mapextent['extent'][1], 'rad')
510 qwidth = my_qa.quantity(mapextent['extent'][0], 'rad')
511 qcent0 = my_qa.quantity(mapextent['center'][0], 'rad')
512 qcent1 = my_qa.quantity(mapextent['center'][1], 'rad')
513 scenter = '%s %s %s' % (base_mref, my_qa.formxxx(qcent0, 'hms'),
514 my_qa.formxxx(qcent1, 'dms'))
516 casalog.post("- Pointing center: %s" % scenter)
517 casalog.post("- Pointing extent: [%s, %s] (projected)" % (my_qa.tos(qwidth),
518 my_qa.tos(qheight)))
519 ret_dict['center'] = scenter
520 ret_dict['width'] = qwidth
521 ret_dict['height'] = qheight
522 else:
523 casalog.post(
524 'Failed to derive map extent from the MSs registered to the imager probably '
525 'due to mising valid data.',
526 priority='SEVERE')
527 ret_dict['center'] = ''
528 ret_dict['width'] = my_qa.quantity(0.0, 'rad')
529 ret_dict['height'] = my_qa.quantity(0.0, 'rad')
530 return ret_dict
533def _handle_image_params(imsize, cell, phasecenter,
534 vislist, field, spw, antenna, scan, intent, timerange,
535 restfreq, pointingcolumntouse, ephemsrcname):
536 logger = sdutil.Casalog(origin="_handle_image_params")
537 # round-up imsize
538 _imsize = sdutil.to_list(imsize, int) or sdutil.to_list(imsize, numpy.integer)
539 if _imsize is None:
540 _imsize = imsize if hasattr(imsize, '__iter__') else [imsize]
541 _imsize = [int(numpy.ceil(v)) for v in _imsize]
542 logger.post(
543 "imsize is not integers. force converting to integer pixel numbers.",
544 priority="WARN")
545 logger.post("rounded-up imsize: %s --> %s" % (str(imsize), str(_imsize)))
547 # calculate cell based on PB if it is not given
548 _cell = cell
549 if _cell == '' or _cell[0] == '':
550 # calc PB
551 if isinstance(vislist, str):
552 vis = vislist
553 else:
554 vis = vislist[0]
555 if isinstance(antenna, str):
556 antsel = antenna
557 else:
558 antsel = antenna[0]
559 if antsel == '':
560 antenna_id = 0
561 else:
562 if len(antsel) > 3 and antsel[:-3] == '&&&':
563 baseline = antsel
564 else:
565 baseline = antsel + '&&&'
566 with open_ms(vis) as ms:
567 ms.msselect({'baseline': baseline})
568 ndx = ms.msselectedindices()
569 antenna_id = ndx['antenna1'][0]
570 grid_factor = 3.
571 logger.post("The cell size will be calculated using PB size of antennas in the first MS")
572 qpb = _calc_PB(vis, antenna_id, restfreq)
573 _cell = '%f%s' % (qpb['value'] / grid_factor, qpb['unit'])
574 logger.post("Using cell size = PB/%4.2F = %s" % (grid_factor, _cell))
576 # Calculate Pointing center and extent (if necessary)
577 _phasecenter = phasecenter
578 if _phasecenter == '' or len(_imsize) == 0 or _imsize[0] < 1:
579 # return a dictionary with keys 'center', 'width', 'height'
580 map_param = _get_pointing_extent(_phasecenter, vislist, field, spw, antenna, scan, intent,
581 timerange, pointingcolumntouse, ephemsrcname)
582 # imsize
583 (cellx, celly) = sdutil.get_cellx_celly(_cell, unit='arcmin')
584 if len(_imsize) == 0 or _imsize[0] < 1:
585 _imsize = _get_imsize(map_param['width'], map_param['height'], cellx, celly)
586 if _phasecenter != "":
587 logger.post(
588 "You defined phasecenter but not imsize. "
589 "The image will cover as wide area as pointing in MS extends, "
590 "but be centered at phasecenter. "
591 "This could result in a strange image if your phasecenter is "
592 "apart from the center of pointings",
593 priority='WARN')
594 if _imsize[0] > 1024 or _imsize[1] > 1024:
595 logger.post(
596 "The calculated image pixel number is larger than 1024. "
597 "It could take time to generate the image depending on your computer resource. "
598 "Please wait...",
599 priority='WARN')
601 # phasecenter
602 # if empty, it should be determined here...
603 if _phasecenter == "":
604 _phasecenter = map_param['center']
606 return _imsize, _cell, _phasecenter
609def _get_param(ms_index, param):
610 if isinstance(param, str):
611 return param
612 elif hasattr(param, '__iter__'):
613 if len(param) == 1:
614 return param[0]
615 else:
616 return param[ms_index]
617 else:
618 raise RuntimeError('Invalid parameter')
621def _remove_image(imagename):
622 if os.path.exists(imagename):
623 if os.path.isdir(imagename):
624 shutil.rmtree(imagename)
625 elif os.path.isfile(imagename):
626 os.remove(imagename)
627 else:
628 # could be a symlink
629 os.remove(imagename)
632def _get_restfreq_if_empty(vislist, spw, field, restfreq):
633 qa = quanta()
634 rf = None
635 # if restfreq is nonzero float value, return it
636 if isinstance(restfreq, float):
637 if restfreq != 0.0:
638 rf = restfreq
639 # if restfreq is valid frequency string, return it
640 # numeric string is interpreted as a value in the unit of Hz
641 elif isinstance(restfreq, str):
642 q = qa.convert(qa.quantity(restfreq), 'Hz')
643 if q['unit'] == 'Hz' and q['value'] > 0.0:
644 rf = restfreq
645 # if restfreq is valid quantity, return it
646 elif isinstance(restfreq, dict):
647 q = qa.convert(restfreq, 'Hz')
648 if q['unit'] == 'Hz' and q['value'] > 0.0:
649 rf = restfreq
651 if isinstance(vislist, str):
652 vis = vislist
653 elif hasattr(vislist, '__iter__'):
654 vis = vislist[0]
655 else:
656 raise RuntimeError('Internal Error: invalid vislist \'{0}\''.format(vislist))
658 if isinstance(spw, str):
659 spwsel = spw
660 elif hasattr(spw, '__iter__'):
661 spwsel = spw[0]
662 else:
663 raise RuntimeError('Internal Error: invalid spw selection \'{0}\''.format(spw))
665 if isinstance(field, str):
666 fieldsel = field
667 elif hasattr(field, '__iter__'):
668 fieldsel = field[0]
669 else:
670 raise RuntimeError('Internal Error: invalid field selection \'{0}\''.format(field))
672 with open_ms(vis) as ms:
673 ms.msselect({'spw': spwsel, 'field': fieldsel})
674 ndx = ms.msselectedindices()
675 if len(ndx['spw']) > 0:
676 spwid = ndx['spw'][0]
677 else:
678 spwid = None
679 if len(ndx['field']) > 0:
680 fieldid = ndx['field'][0]
681 else:
682 fieldid = None
683 sourceid = None
684 if fieldid is not None:
685 with sdutil.table_manager(os.path.join(vis, 'FIELD')) as tb:
686 sourceid = tb.getcell('SOURCE_ID', fieldid)
687 if sourceid < 0:
688 sourceid = None
689 if rf is None:
690 # if restfrequency is defined in SOURCE table, return it
691 with sdutil.table_manager(os.path.join(vis, 'SOURCE')) as tb:
692 if 'REST_FREQUENCY' in tb.colnames():
693 tsel = None
694 taql = ''
695 if spwid is not None:
696 taql = 'SPECTRAL_WINDOW_ID == {0}'.format(spwid)
697 if sourceid is not None:
698 delimiter = '&&' if len(taql) > 0 else ''
699 taql += '{0}SOURCE_ID == {1}'.format(delimiter, sourceid)
700 if len(taql) > 0:
701 tsel = tb.query(taql)
702 t = tsel
703 else:
704 t = tb
705 try:
706 nrow = t.nrows()
707 if nrow > 0:
708 for irow in range(nrow):
709 if t.iscelldefined('REST_FREQUENCY', irow):
710 rfs = t.getcell('REST_FREQUENCY', irow)
711 if len(rfs) > 0:
712 rf = rfs[0]
713 break
714 finally:
715 if tsel is not None:
716 tsel.close()
718 if rf is None:
719 if spwid is None:
720 spwid = 0
721 # otherwise, return mean frequency of given spectral window
722 with sdutil.table_manager(os.path.join(vis, 'SPECTRAL_WINDOW')) as tb:
723 cf = tb.getcell('CHAN_FREQ', spwid)
724 rf = cf.mean()
726 assert rf is not None
728 return rf
731def set_beam_size(vis, imagename,
732 field, spw, baseline, scan, intent, timerange,
733 ephemsrcname, pointingcolumntouse, antenna_name, antenna_diameter,
734 restfreq, gridfunction, convsupport, truncate, gwidth, jwidth):
735 """Set estimated beam size to the image."""
736 is_alma = antenna_name[0:2] in ['PM', 'DV', 'DA', 'CM']
737 blockage = '0.75m' if is_alma else '0.0m'
738 log_origin = 'set_beam_size'
740 with open_ia(imagename) as ia:
741 csys = ia.coordsys()
742 outref = csys.referencecode('direction')[0]
743 cell = list(csys.increment(type='direction', format='s')['string'])
744 csys.done()
746 old_tool = OldImagerBasedTools()
747 sampling_params = old_tool.get_pointing_sampling_params(vis, field, spw, baseline,
748 scan, intent, timerange,
749 outref=outref,
750 movingsource=ephemsrcname,
751 pointingcolumntouse=pointingcolumntouse,
752 antenna_name=antenna_name)
753 qa = quanta()
754 casalog.post(
755 f'sampling_params={sampling_params}',
756 origin=log_origin
757 )
758 xsampling, ysampling = qa.getvalue(qa.convert(sampling_params['sampling'], 'arcsec'))
759 angle = qa.getvalue(qa.convert(sampling_params['angle'], 'deg'))[0]
761 casalog.post(
762 f'Detected raster sampling = [{xsampling:f}, {ysampling:f}] arcsec',
763 origin=log_origin
764 )
766 # handling of failed sampling detection
767 valid_sampling = True
768 # TODO: copy from sdimaging implementation
769 sampling = [xsampling, ysampling]
770 if abs(xsampling) < 2.2e-3 or not numpy.isfinite(xsampling):
771 casalog.post(
772 f"Invalid sampling={xsampling} arcsec. "
773 f"Using the value of orthogonal direction={ysampling} arcsec",
774 priority="WARN",
775 origin=log_origin
776 )
777 sampling = [ysampling]
778 angle = 0.0
779 valid_sampling = False
780 if abs(ysampling) < 1.0e-3 or not numpy.isfinite(ysampling):
781 if valid_sampling:
782 casalog.post(
783 f"Invalid sampling={ysampling} arcsec. "
784 f"Using the value of orthogonal direction={xsampling} arcsec",
785 priority="WARN",
786 origin=log_origin
787 )
788 sampling = [xsampling]
789 angle = 0.0
790 valid_sampling = True
791 # reduce sampling and cell if it's possible
792 if len(sampling) > 1 and abs(sampling[0] - sampling[1]) <= 0.01 * abs(sampling[0]):
793 sampling = [sampling[0]]
794 angle = 0.0
795 if cell[0] == cell[1]:
796 cell = [cell[0]]
797 if valid_sampling:
798 # actual calculation of beam size
799 bu = sdbeamutil.TheoreticalBeam()
800 bu.set_antenna(antenna_diameter, blockage)
801 bu.set_sampling(sampling, "%fdeg" % angle)
802 bu.set_image_param(cell, restfreq, gridfunction,
803 convsupport, truncate, gwidth,
804 jwidth, is_alma)
805 bu.summary()
806 imbeam_dict = bu.get_beamsize_image()
807 casalog.post(
808 f"Setting image beam: "
809 f"major={imbeam_dict['major']}, "
810 f"minor={imbeam_dict['minor']}, "
811 f"pa={imbeam_dict['pa']}",
812 origin=log_origin
813 )
814 # set beam size to image
815 with open_ia(imagename) as ia:
816 ia.setrestoringbeam(**imbeam_dict)
817 else:
818 # BOTH sampling was invalid
819 casalog.post(
820 "Could not detect valid raster sampling. "
821 "Exiting without setting beam size to image",
822 priority='WARN',
823 origin=log_origin
824 )
827def do_weight_mask(imagename, weightimage, minweight):
828 # Mask image pixels whose weight are smaller than minweight.
829 # Weight image should have 0 weight for pixels below < minweight
830 logger = sdutil.Casalog(origin="do_weight_mask")
831 logger.post(f"Start masking the map using minweight = {minweight:f}",
832 priority="INFO")
833 with open_ia(weightimage) as ia:
834 try:
835 stat = ia.statistics(mask="'" + weightimage + "' > 0.0", robust=True)
836 valid_pixels = stat['npts']
837 except RuntimeError as e:
838 if 'No valid data found.' in str(e):
839 valid_pixels = [0]
840 else:
841 raise e
843 if len(valid_pixels) == 0 or valid_pixels[0] == 0:
844 logger.post(
845 "All pixels weight zero. "
846 "This indicates no data in MS is in image area. "
847 "Mask will not be set. Please check your image parameters.",
848 priority="WARN")
849 return
850 median_weight = stat['median'][0]
851 weight_threshold = median_weight * minweight
852 logger.post(f"Median of weight in the map is {median_weight:f}",
853 priority="INFO")
854 logger.post(f"Pixels in map with weight <= median(weight)*minweight = "
855 "{weight_threshold:f} will be masked.",
856 priority="INFO")
857 # Leaving the original logic to calculate the number of masked pixels via
858 # product of median of and min_weight (which i don't understand the logic)
860 # Modify default mask
861 with open_ia(imagename) as ia:
862 ia.calcmask("'%s'>%f" % (weightimage, weight_threshold), asdefault=True)
864 ndim = len(ia.shape())
865 _axes = numpy.arange(start=0 if ndim <= 2 else 2, stop=ndim)
867 try:
868 collapsed = ia.collapse('npts', axes=_axes)
869 valid_pixels_after = collapsed.getchunk().sum()
870 collapsed.close()
871 except RuntimeError as e:
872 if 'All selected pixels are masked' in str(e):
873 valid_pixels_after = 0
874 else:
875 raise
877 masked_fraction = 100. * (1. - valid_pixels_after / float(valid_pixels[0]))
879 logger.post(f"This amounts to {masked_fraction:5.1f} % "
880 "of the area with nonzero weight.",
881 priority="INFO")
882 logger.post(
883 f"The weight image '{weightimage}' is returned by this task, "
884 "if the user wishes to assess the results in detail.",
885 priority="INFO")
888def get_ms_column_unit(tb, colname):
889 col_unit = ''
890 if colname in tb.colnames():
891 cdkw = tb.getcoldesc(colname)['keywords']
892 if 'QuantumUnits' in cdkw:
893 u = cdkw['QuantumUnits']
894 if isinstance(u, str):
895 col_unit = u.strip()
896 elif isinstance(u, list):
897 col_unit = u[0].strip()
898 return col_unit
901def get_brightness_unit_from_ms(msname):
902 image_unit = ''
903 with sdutil.table_manager(msname) as tb:
904 image_unit = get_ms_column_unit(tb, 'DATA')
905 if image_unit == '':
906 image_unit = get_ms_column_unit(tb, 'FLOAT_DATA')
907 if image_unit.upper() == 'K':
908 image_unit = 'K'
909 else:
910 image_unit = 'Jy/beam'
912 return image_unit
915@sdutil.sdtask_decorator
916def tsdimaging(
917 # Input data: list of MeasurementSets
918 infiles,
919 # Output data: CASA images: their path prefix, overwrite control
920 outfile, overwrite,
921 # Select data from input MeasurementSets, by
922 field, spw, antenna, scan, intent, timerange,
923 # Output images definition: frequency axis
924 outframe, # velocity frame
925 mode, nchan, start, width, veltype, # gridding type
926 specmode, # Doppler handling
927 interpolation, # interpolation mode
928 # Output images definition: spatial axes
929 pointingcolumn, convertfirst,
930 projection,
931 imsize, cell, phasecenter,
932 # Output images definition: stokes axis
933 stokes,
934 # Gridder parameters
935 gridfunction, convsupport, truncate, gwidth, jwidth,
936 clipminmax,
937 # Single-dish image: mask control
938 minweight,
939 # Single-dish image: metadata
940 brightnessunit,
941 restfreq # rest frequency to assign to image
942 ):
944 origin = 'tsdimaging'
945 imager = None
947 try: # Create the Single-Dish Image
948 # Validate brightnessunit parameter CAS-11503
949 image_unit = brightnessunit.lower().capitalize()
950 if not image_unit in ['', 'K', 'Jy/beam']:
951 raise ValueError(f"Invalid brightness unit: {brightnessunit}")
953 # Handle outfile and overwrite parameters
954 output_path_prefix = outfile.rstrip('/')
955 singledish_image_path = output_path_prefix + image_suffix
956 if os.path.exists(singledish_image_path):
957 if overwrite == False:
958 raise RuntimeError(
959 f"Output file exists: '{singledish_image_path}'"
960 )
961 else:
962 # delete existing images
963 if os.path.exists(singledish_image_path):
964 casalog.post(f"Removing '{singledish_image_path}'")
965 _remove_image(singledish_image_path)
966 assert not os.path.exists(singledish_image_path)
967 for _suffix in associate_suffixes:
968 path_to_remove = output_path_prefix + _suffix
969 if os.path.exists(path_to_remove):
970 casalog.post(f"Removing '{path_to_remove}'")
971 _remove_image(path_to_remove)
972 assert not os.path.exists(path_to_remove)
974 # Tweak spw parameter into _spw
975 if isinstance(spw, str):
976 _spw = '*' + spw if spw.startswith(':') else spw
977 else:
978 _spw = ['*' + v if v.startswith(':') else v for v in spw]
980 # Handle image spectral axis parameters
981 imnchan, imstart, imwidth = _configure_spectral_axis(
982 mode, nchan, start, width, restfreq
983 )
985 # Handle image restfreq parameter's default value
986 _restfreq = _get_restfreq_if_empty(
987 infiles, _spw, field, restfreq
988 )
990 # Handle gridder parameters
991 # convert type of task's default values
992 # to ones supported by the Synthesis Imager framework
993 gtruncate = _handle_grid_defaults(truncate)
994 ggwidth = _handle_grid_defaults(gwidth)
995 gjwidth = _handle_grid_defaults(jwidth)
997 # handle infiles parameter
998 # sort input data to get results consistent with old sdimaging task
999 _sorted = sort_vis(
1000 infiles, _spw, mode, imwidth, field,
1001 antenna, scan, intent, timerange
1002 )
1003 sorted_vis = _sorted[0]
1004 sorted_field = _sorted[1]
1005 sorted_spw = _sorted[2]
1006 sorted_antenna = _sorted[3]
1007 sorted_scan = _sorted[4]
1008 sorted_intent = _sorted[5]
1009 sorted_timerange = _sorted[6]
1011 def antenna_to_baseline(s):
1012 if len(s) == 0:
1013 return s
1014 elif len(s) > 3 and s.endswith('&&&'):
1015 return s
1016 else:
1017 return '{0}&&&'.format(s)
1019 if isinstance(sorted_antenna, str):
1020 sorted_baseline = antenna_to_baseline(sorted_antenna)
1021 else:
1022 sorted_baseline = [antenna_to_baseline(a) for a in sorted_antenna]
1024 # Handle image geometric parameters
1025 _ephemsrcname = ''
1026 ephem_sources = ['MERCURY', 'VENUS', 'MARS', 'JUPITER', 'SATURN',
1027 'URANUS', 'NEPTUNE', 'PLUTO', 'SUN', 'MOON', 'TRACKFIELD']
1028 if isinstance(phasecenter, str) and phasecenter.strip().upper() in ephem_sources:
1029 _ephemsrcname = phasecenter
1030 _imsize, _cell, _phasecenter = _handle_image_params(
1031 imsize, cell, phasecenter, sorted_vis,
1032 sorted_field, sorted_spw, sorted_antenna,
1033 sorted_scan, sorted_intent, sorted_timerange,
1034 _restfreq, pointingcolumn, _ephemsrcname
1035 )
1037 # Set up PySynthesisImager input parameters
1038 # - List all parameters that you need here
1039 # - Defaults will be assumed for unspecified parameters
1040 # - Nearly all parameters are identical to that in the task. Please look at the
1041 # list of parameters under __init__ using " help ImagerParameters " )
1042 casalog.post('*** Creating paramList ***', origin=origin)
1043 paramList = ImagerParameters(
1044 # input file name
1045 msname=sorted_vis, # 'sdimaging.ms',
1046 # data selection
1047 field=sorted_field, # '',
1048 spw=sorted_spw, # '0',
1049 timestr=sorted_timerange,
1050 antenna=sorted_baseline,
1051 scan=sorted_scan,
1052 state=sorted_intent,
1053 # image parameters
1054 imagename=output_path_prefix, # 'try2',
1055 nchan=imnchan, # 1024,
1056 start=imstart, # '0',
1057 width=imwidth, # '1',
1058 outframe=outframe,
1059 veltype=veltype,
1060 restfreq=_restfreq,
1061 phasecenter=_phasecenter, # 'J2000 17:18:29 +59.31.23',
1062 imsize=_imsize, # [75,75],
1063 cell=_cell, # ['3arcmin', '3arcmin'],
1064 projection=projection,
1065 stokes=stokes,
1066 specmode=specmode,
1067 gridder='singledish',
1068 # single dish specific parameters
1069 gridfunction=gridfunction,
1070 convsupport=convsupport,
1071 truncate=gtruncate,
1072 gwidth=ggwidth,
1073 jwidth=gjwidth,
1074 pointingcolumntouse=pointingcolumn,
1075 convertfirst=convertfirst,
1076 minweight=minweight,
1077 clipminmax=clipminmax,
1078 # normalizer
1079 normtype='flatsky',
1080 pblimit=1e-16, # TODO: explain why 1e-16 ?
1081 interpolation=interpolation
1082 )
1084 # Construct the PySynthesisImager object, with all input parameters
1085 casalog.post('*** Creating imager object ***', origin=origin)
1086 imager = PySynthesisImager(params=paramList)
1088 # Initialize PySynthesisImager "modules" required for single-dish imaging
1089 # - Pick only the modules you will need later on. For example, to only make
1090 # the PSF, there is no need for the deconvolver or iteration control modules.
1091 casalog.post('*** Initializing imagers ***', origin=origin)
1092 # This is where the underlying C++ synthesis imager is created
1093 imager.initializeImagers()
1094 casalog.post('*** Initializing normalizers ***', origin=origin)
1095 imager.initializeNormalizers()
1097 # Create Single-Dish images
1098 casalog.post('*** Creating single-dish images ***', origin=origin)
1099 imager.makeSdImage()
1100 casalog.post('*** Created single-dish images ***', origin=origin)
1102 finally: # Close tools and rename Synthesis Imager's residual image
1103 casalog.post('*** Cleaning up tools ***', origin=origin)
1104 if imager is not None:
1105 imager.deleteTools()
1106 # Change image suffix from .residual to .image
1107 residual_image_path = output_path_prefix + residual_suffix
1108 if os.path.exists(residual_image_path):
1109 os.rename(residual_image_path, singledish_image_path)
1111 # Set single-dish image's beam size
1112 # TODO: re-define related functions in the new tool framework (sdms?)
1113 # ms_index = 0
1114 rep_ms = _get_param(0, infiles)
1115 rep_field = _get_param(0, field)
1116 rep_spw = _get_param(0, _spw)
1117 rep_antenna = _get_param(0, antenna)
1118 rep_scan = _get_param(0, scan)
1119 rep_intent = _get_param(0, intent)
1120 rep_timerange = _get_param(0, timerange)
1121 if len(rep_antenna) > 0:
1122 baseline = '{0}&&&'.format(rep_antenna)
1123 else:
1124 baseline = '*&&&'
1125 with open_ms(rep_ms) as ms:
1126 ms.msselect({'baseline': baseline})
1127 ndx = ms.msselectedindices()
1128 antenna_index = ndx['antenna1'][0]
1129 with sdutil.table_manager(os.path.join(rep_ms, 'ANTENNA')) as tb:
1130 antenna_name = tb.getcell('NAME', antenna_index)
1131 antenna_diameter = tb.getcell('DISH_DIAMETER', antenna_index)
1132 casalog.post(f"Setting single-dish image's beam")
1133 set_beam_size(
1134 rep_ms, singledish_image_path,
1135 rep_field, rep_spw, baseline, rep_scan, rep_intent, rep_timerange,
1136 _ephemsrcname, pointingcolumn, antenna_name, antenna_diameter,
1137 _restfreq, gridfunction, convsupport, truncate, gwidth, jwidth
1138 )
1140 # Set single-dish image's brightness unit (CAS-11503)
1141 if len(image_unit) == 0:
1142 image_unit = get_brightness_unit_from_ms(rep_ms)
1143 if len(image_unit) > 0:
1144 with open_ia(singledish_image_path) as ia:
1145 casalog.post(f"Setting single-dish image's brightness unit to '{image_unit}'")
1146 ia.setbrightnessunit(image_unit)
1148 # Update single-dish image's mask: mask low weight pixels
1149 weight_image_path = output_path_prefix + weight_suffix
1150 casalog.post(f"Creating weight image mask")
1151 do_weight_mask(singledish_image_path, weight_image_path, minweight)
1153 # Delete images systematically generated by the Synthesis Imager
1154 # which are either not required or currently useless
1155 # in the context of single-dish imaging
1156 # CAS-10891
1157 _remove_image(output_path_prefix + '.sumwt')
1158 # CAS-10893
1159 # TODO: remove the following line once the 'correct' SD
1160 # PSF image based on primary beam can be generated
1161 _remove_image(output_path_prefix + '.psf')