Commits
Darrell Schiebel authored 8d246bd2b2a Merge
96 96 | |
97 97 | contextmanager | .
98 98 | def open_old_imager(self, vis): |
99 99 | try: |
100 100 | self.imager.open(vis) |
101 101 | yield self.imager |
102 102 | finally: |
103 103 | self.imager.close() |
104 104 | |
105 105 | contextmanager | .
106 - | def open_and_select_old_imager(self, vislist, field, spw, antenna, scan, intent): |
106 + | def open_and_select_old_imager(self, vislist, field, spw, antenna, scan, intent, timerange): |
107 107 | if isinstance(vislist, str): |
108 108 | with self.open_old_imager(vislist) as im: |
109 109 | im.selectvis(field=field, |
110 110 | spw=spw, |
111 111 | nchan=-1, |
112 112 | start=0, |
113 113 | step=1, |
114 114 | baseline=antenna, |
115 115 | scan=scan, |
116 - | intent=intent) |
116 + | intent=intent, |
117 + | time=timerange) |
117 118 | yield im |
118 119 | else: |
119 120 | fieldsel = SelectionHandler(field) |
120 121 | spwsel = SelectionHandler(spw) |
121 122 | antennasel = SelectionHandler(antenna) |
122 123 | scansel = SelectionHandler(scan) |
123 124 | intentsel = SelectionHandler(intent) |
125 + | timerangesel = SelectionHandler(timerange) |
124 126 | try: |
125 127 | for i in range(len(vislist)): |
126 128 | vis = vislist[i] |
127 129 | _field = fieldsel(i) |
128 130 | _spw = spwsel(i) |
129 131 | _antenna = antennasel(i) |
130 132 | _scan = scansel(i) |
131 133 | _intent = intentsel(i) |
134 + | _timerangesel = timerangesel(i) |
132 135 | if len(_antenna) == 0: |
133 136 | _baseline = _antenna |
134 137 | elif len(_antenna) < 4 or _antenna[:-3] != '&&&': |
135 138 | _baseline = _antenna + '&&&' |
136 139 | else: |
137 140 | _baseline = _antenna |
138 141 | self.imager.selectvis(vis, field=_field, spw=_spw, nchan=-1, start=0, step=1, |
139 - | baseline=_baseline, scan=_scan, intent=_intent) |
142 + | baseline=_baseline, scan=_scan, intent=_intent, time=_timerangesel) |
140 143 | yield self.imager |
141 144 | finally: |
142 145 | self.imager.close() |
143 146 | |
144 147 | def test(self, vis): |
145 148 | with self.open_old_imager(vis) as im: |
146 149 | casalog.post('test') |
147 150 | raise RuntimeError('ERROR!') |
148 151 | |
149 - | def get_pointing_sampling_params(self, vis, field, spw, baseline, scan, intent, outref, movingsource, pointingcolumntouse, antenna_name): |
152 + | def get_pointing_sampling_params(self, vis, field, spw, baseline, scan, intent, timerange, outref, movingsource, pointingcolumntouse, antenna_name): |
150 153 | with self.open_old_imager(vis) as im: |
151 154 | im.selectvis(field=field, |
152 155 | spw=spw, |
153 156 | nchan=-1, |
154 157 | start=0, |
155 158 | step=1, |
156 159 | baseline=baseline, |
157 160 | scan=scan, |
158 - | intent=intent) |
161 + | intent=intent, |
162 + | time=timerange) |
159 163 | sampling_params = im.pointingsampling(pattern='raster', |
160 164 | ref=outref, |
161 165 | movingsource=movingsource, |
162 166 | pointingcolumntouse=pointingcolumntouse, |
163 167 | antenna='{0}&&&'.format(antenna_name)) |
164 168 | return sampling_params |
165 169 | |
166 - | def get_map_extent(self, vislist, field, spw, antenna, scan, intent, |
170 + | def get_map_extent(self, vislist, field, spw, antenna, scan, intent, timerange, |
167 171 | ref, movingsource, pointingcolumntouse): |
168 172 | |
169 173 | with self.open_and_select_old_imager(vislist=vislist, field=field, |
170 174 | spw=spw, antenna=antenna, scan=scan, |
171 - | intent=intent) as im: |
175 + | intent=intent, timerange=timerange) as im: |
172 176 | map_param = im.mapextent(ref=ref, movingsource=movingsource, |
173 177 | pointingcolumntouse=pointingcolumntouse) |
174 178 | return map_param |
175 179 | |
176 - | def sort_vis(self, vislist, spw, mode, width, field, antenna, scan, intent): |
180 + | def sort_vis(self, vislist, spw, mode, width, field, antenna, scan, intent, timerange): |
177 181 | if isinstance(vislist, str) or len(vislist) == 1: |
178 - | return vislist, field, spw, antenna, scan, intent |
182 + | return vislist, field, spw, antenna, scan, intent, timerange |
179 183 | imhelper = cleanhelper(imtool=self.imager, vis=vislist, casalog=casalog) |
180 184 | imhelper.sortvislist(spw=spw, mode=mode, width=width) |
181 185 | sorted_idx = list(imhelper.sortedvisindx) |
182 186 | # reverse the order |
183 187 | sorted_idx.reverse() |
184 188 | sorted_vislist = [vislist[i] for i in sorted_idx] |
185 189 | fieldsel = SelectionHandler(field) |
186 190 | sorted_field = [fieldsel(i) for i in sorted_idx] |
187 191 | spwsel = SelectionHandler(spw) |
188 192 | sorted_spw = [spwsel(i) for i in sorted_idx] |
189 193 | antennasel = SelectionHandler(antenna) |
190 194 | sorted_antenna = [antennasel(i) for i in sorted_idx] |
191 195 | scansel = SelectionHandler(scan) |
192 196 | sorted_scan = [scansel(i) for i in sorted_idx] |
193 197 | intentsel = SelectionHandler(intent) |
194 198 | sorted_intent = [intentsel(i) for i in sorted_idx] |
195 - | return sorted_vislist, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent |
199 + | timerangesel = SelectionHandler(timerange) |
200 + | sorted_timerange = [timerangesel(i) for i in sorted_idx] |
201 + | return sorted_vislist, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent, sorted_timerange |
196 202 | |
197 203 | def _configure_spectral_axis(mode, nchan, start, width, restfreq): |
198 204 | # fix default |
199 205 | if mode == 'channel': |
200 206 | if start == '': start = 0 |
201 207 | if width == '': width = 1 |
202 208 | else: |
203 209 | if start == 0: start = '' |
204 210 | if width == 1: width = '' |
205 211 | # fix unit |
297 303 | ny = numpy.ceil( ( my_qa.convert(height, my_qa.getunit(dy))['value'] / \ |
298 304 | my_qa.getvalue(dy) ) ) |
299 305 | nx = numpy.ceil( ( my_qa.convert(width, my_qa.getunit(dx))['value'] / \ |
300 306 | my_qa.getvalue(dx) ) ) |
301 307 | casalog.post("- Map extent: [%s, %s]" % (my_qa.tos(width), my_qa.tos(height))) |
302 308 | casalog.post("- Cell size: [%s, %s]" % (my_qa.tos(dx), my_qa.tos(dy))) |
303 309 | casalog.post("Image pixel numbers to cover the extent: [%d, %d] (projected)" % \ |
304 310 | (nx+1, ny+1)) |
305 311 | return [int(nx+1), int(ny+1)] |
306 312 | |
307 - | def _get_pointing_extent(phasecenter, vislist, field, spw, antenna, scan, intent, |
313 + | def _get_pointing_extent(phasecenter, vislist, field, spw, antenna, scan, intent, timerange, |
308 314 | pointingcolumntouse, ephemsrcname): |
309 315 | ### MS selection is ignored. This is not quite right. |
310 316 | casalog.post("Calculating map extent from pointings.") |
311 317 | # CAS-5410 Use private tools inside task scripts |
312 318 | my_qa = quanta() |
313 319 | ret_dict = {} |
314 320 | |
315 321 | if isinstance(vislist, str): |
316 322 | vis = vislist |
317 323 | else: |
331 337 | pattern = '^([\-\+]?[0-9.]+([eE]?-?[0-9])?)|([\-\+]?[0-9][:h][0-9][:m][0-9.]s?)|([\-\+]?[0-9][.d][0-9][.d][0-9.]s?)$' |
332 338 | items = phasecenter.split() |
333 339 | base_mref = 'J2000' |
334 340 | for i in items: |
335 341 | s = i.strip() |
336 342 | if re.match(pattern, s) is None: |
337 343 | base_mref = s |
338 344 | break |
339 345 | |
340 346 | t = OldImagerBasedTools() |
341 - | mapextent = t.get_map_extent(vislist, field, spw, antenna, scan, intent, |
347 + | mapextent = t.get_map_extent(vislist, field, spw, antenna, scan, intent, timerange, |
342 348 | ref=base_mref, movingsource=ephemsrcname, |
343 349 | pointingcolumntouse=pointingcolumntouse) |
344 350 | #mapextent = self.imager.mapextent(ref=base_mref, movingsource=ephemsrcname, |
345 351 | # pointingcolumntouse=colname) |
346 352 | if mapextent['status'] is True: |
347 353 | qheight = my_qa.quantity(mapextent['extent'][1], 'rad') |
348 354 | qwidth = my_qa.quantity(mapextent['extent'][0], 'rad') |
349 355 | qcent0 = my_qa.quantity(mapextent['center'][0], 'rad') |
350 356 | qcent1 = my_qa.quantity(mapextent['center'][1], 'rad') |
351 357 | scenter = '%s %s %s'%(base_mref, my_qa.formxxx(qcent0, 'hms'), |
358 364 | ret_dict['width'] = qwidth |
359 365 | ret_dict['height'] = qheight |
360 366 | else: |
361 367 | casalog.post('Failed to derive map extent from the MSs registered to the imager probably due to mising valid data.', priority='SEVERE') |
362 368 | ret_dict['center'] = '' |
363 369 | ret_dict['width'] = my_qa.quantity(0.0, 'rad') |
364 370 | ret_dict['height'] = my_qa.quantity(0.0, 'rad') |
365 371 | return ret_dict |
366 372 | |
367 373 | def _handle_image_params(imsize, cell, phasecenter, |
368 - | vislist, field, spw, antenna, scan, intent, |
374 + | vislist, field, spw, antenna, scan, intent, timerange, |
369 375 | restfreq, pointingcolumntouse, ephemsrcname): |
370 376 | # round-up imsize |
371 377 | _imsize = sdutil._to_list(imsize, int) or sdutil._to_list(imsize, numpy.integer) |
372 378 | if _imsize is None: |
373 379 | _imsize = imsize if hasattr(imsize, '__iter__') else [ imsize ] |
374 380 | _imsize = [ int(numpy.ceil(v)) for v in _imsize ] |
375 381 | casalog.post("imsize is not integers. force converting to integer pixel numbers.", priority="WARN") |
376 382 | casalog.post("rounded-up imsize: %s --> %s" % (str(imsize), str(_imsize))) |
377 383 | |
378 384 | # calculate cell based on PB if it is not given |
401 407 | grid_factor = 3. |
402 408 | casalog.post("The cell size will be calculated using PB size of antennas in the first MS") |
403 409 | qpb = _calc_PB(vis, antenna_id, restfreq) |
404 410 | _cell = '%f%s' % (qpb['value']/grid_factor, qpb['unit']) |
405 411 | casalog.post("Using cell size = PB/%4.2F = %s" % (grid_factor, _cell)) |
406 412 | |
407 413 | # Calculate Pointing center and extent (if necessary) |
408 414 | _phasecenter = phasecenter |
409 415 | if _phasecenter == '' or len(_imsize) == 0 or _imsize[0] < 1: |
410 416 | # return a dictionary with keys 'center', 'width', 'height' |
411 - | map_param = _get_pointing_extent(_phasecenter, vislist, field, spw, antenna, scan, intent, |
417 + | map_param = _get_pointing_extent(_phasecenter, vislist, field, spw, antenna, scan, intent, timerange, |
412 418 | pointingcolumntouse, ephemsrcname) |
413 419 | # imsize |
414 420 | (cellx,celly) = sdutil.get_cellx_celly(_cell, unit='arcmin') |
415 421 | if len(_imsize) == 0 or _imsize[0] < 1: |
416 422 | _imsize = _get_imsize(map_param['width'], map_param['height'], cellx, celly) |
417 423 | if _phasecenter != "": |
418 424 | casalog.post("You defined phasecenter but not imsize. The image will cover as wide area as pointing in MS extends, but be centered at phasecenter. This could result in a strange image if your phasecenter is a part from the center of pointings", priority='WARN') |
419 425 | if _imsize[0] > 1024 or _imsize[1] > 1024: |
420 426 | casalog.post("The calculated image pixel number is larger than 1024. It could take time to generate the image depending on your computer resource. Please wait...", priority='WARN') |
421 427 | |
551 557 | # otherwise, return mean frequency of given spectral window |
552 558 | with open_table(os.path.join(vis, 'SPECTRAL_WINDOW')) as tb: |
553 559 | cf = tb.getcell('CHAN_FREQ', spwid) |
554 560 | rf = cf.mean() |
555 561 | |
556 562 | assert rf is not None |
557 563 | |
558 564 | return rf |
559 565 | |
560 566 | def set_beam_size(vis, imagename, |
561 - | field, spw, baseline, scan, intent, |
567 + | field, spw, baseline, scan, intent, timerange, |
562 568 | ephemsrcname, pointingcolumntouse, antenna_name, antenna_diameter, |
563 569 | restfreq, gridfunction, convsupport, truncate, gwidth, jwidth): |
564 570 | """ |
565 571 | Set estimated beam size to the image. |
566 572 | """ |
567 573 | is_alma = antenna_name[0:2] in ['PM', 'DV', 'DA', 'CM'] |
568 574 | blockage = '0.75m' if is_alma else '0.0m' |
569 575 | |
570 576 | with open_ia(imagename) as ia: |
571 577 | csys = ia.coordsys() |
572 578 | outref = csys.referencecode('direction')[0] |
573 579 | cell = list(csys.increment(type='direction', format='s')['string']) |
574 580 | csys.done() |
575 581 | |
576 582 | old_tool = OldImagerBasedTools() |
577 - | sampling_params = old_tool.get_pointing_sampling_params(vis, field, spw, baseline, scan, intent, |
583 + | sampling_params = old_tool.get_pointing_sampling_params(vis, field, spw, baseline, scan, intent, timerange, |
578 584 | outref=outref, |
579 585 | movingsource=ephemsrcname, |
580 586 | pointingcolumntouse=pointingcolumntouse, |
581 587 | antenna_name=antenna_name) |
582 588 | qa = quanta() |
583 589 | casalog.post('sampling_params={0}'.format(sampling_params)) |
584 590 | xsampling, ysampling = qa.getvalue(qa.convert(sampling_params['sampling'], 'arcsec')) |
585 591 | angle = qa.getvalue(qa.convert(sampling_params['angle'], 'deg'))[0] |
586 592 | |
587 593 | casalog.post('Detected raster sampling = [{0:f}, {1:f}] arcsec'.format(xsampling, ysampling)) |
696 702 | if image_unit == '': image_unit = get_ms_column_unit(tb, 'FLOAT_DATA') |
697 703 | if image_unit.upper() == 'K': |
698 704 | image_unit = 'K' |
699 705 | else: |
700 706 | image_unit = 'Jy/beam' |
701 707 | |
702 708 | return image_unit |
703 709 | |
704 710 | |
705 711 | |
706 - | def tsdimaging(infiles, outfile, overwrite, field, spw, antenna, scan, intent, mode, nchan, start, width, veltype, |
712 + | def tsdimaging(infiles, outfile, overwrite, field, spw, antenna, scan, intent, timerange, mode, nchan, start, width, veltype, |
707 713 | specmode, outframe, |
708 714 | gridfunction, convsupport, truncate, gwidth, jwidth, imsize, cell, phasecenter, projection, |
709 715 | pointingcolumn, restfreq, stokes, minweight, brightnessunit, clipminmax): |
710 716 | |
711 717 | origin = 'tsdimaging' |
712 718 | imager = None |
713 719 | |
714 720 | try: |
715 - | # if spw starts with ':', add '*' at the beginning |
721 + | # tweak input parameters |
722 + | # ---- if spw starts with ':', add '*' at the beginning |
716 723 | if isinstance(spw, str): |
717 724 | _spw = '*' + spw if spw.startswith(':') else spw |
718 725 | else: |
719 726 | _spw = ['*' + v if v.startswith(':') else v for v in spw] |
720 727 | |
721 - | # if antenna doesn't contain '&&&', append it |
728 + | # ---- if antenna doesn't contain '&&&', append it |
722 729 | def antenna_to_baseline(s): |
723 730 | if len(s) == 0: |
724 731 | return s |
725 732 | elif len(s) > 3 and s.endswith('&&&'): |
726 733 | return s |
727 734 | else: |
728 735 | return '{0}&&&'.format(s) |
729 736 | if isinstance(antenna, str): |
730 737 | baseline = antenna_to_baseline(antenna) |
731 738 | else: |
732 739 | baseline = [antenna_to_baseline(a) for a in antenna] |
733 740 | |
734 - | |
735 741 | # handle overwrite parameter |
736 742 | _outfile = outfile.rstrip('/') |
737 743 | presumed_imagename = _outfile + image_suffix |
738 744 | if os.path.exists(presumed_imagename): |
739 745 | if overwrite == False: |
740 746 | raise RuntimeError('Output file \'{0}\' exists.'.format(presumed_imagename)) |
741 747 | else: |
742 748 | # delete existing images |
743 749 | casalog.post('Removing \'{0}\''.format(presumed_imagename)) |
744 750 | _remove_image(presumed_imagename) |
745 751 | assert not os.path.exists(presumed_imagename) |
746 752 | for _suffix in associate_suffixes: |
747 753 | casalog.post('Removing \'{0}\''.format(_outfile + _suffix)) |
748 754 | _remove_image(_outfile + _suffix) |
749 755 | assert not os.path.exists(_outfile + _suffix) |
750 756 | |
751 - | # parse parameter for spectral axis |
757 + | # handle image spectral axis parameters |
752 758 | imnchan, imstart, imwidth = _configure_spectral_axis(mode, nchan, start, width, restfreq) |
759 + | |
760 + | # handle image restfreq parameter's default value |
753 761 | _restfreq = _get_restfreq_if_empty(infiles, _spw, field, restfreq) |
754 762 | |
755 - | # translate some default values into the ones that are consistent with the current framework |
763 + | # handle gridder parameters |
764 + | # ---- translate some default values into the ones that are consistent with the current framework |
756 765 | gtruncate = _handle_grid_defaults(truncate) |
757 766 | ggwidth = _handle_grid_defaults(gwidth) |
758 767 | gjwidth = _handle_grid_defaults(jwidth) |
759 768 | |
769 + | # handle infiles parameter |
770 + | # ---- sort input data using cleanhelper function to get results consistent with older sdimaging task |
771 + | old_way = OldImagerBasedTools() |
772 + | _sorted = old_way.sort_vis(infiles, _spw, mode, imwidth, field, antenna, scan, intent, timerange) |
773 + | sorted_vis, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent, sorted_timerange = _sorted |
774 + | |
775 + | # handle image geometric parameters |
760 776 | _ephemsrcname = '' |
761 - | if isinstance(phasecenter, str) and phasecenter.strip().upper() in ['MERCURY', 'VENUS', 'MARS', 'JUPITER', 'SATURN', 'URANUS', 'NEPTUNE', 'PLUTO', 'SUN', 'MOON', 'TRACKFIELD']: |
777 + | ephem_sources = ['MERCURY', 'VENUS', 'MARS', 'JUPITER', 'SATURN', 'URANUS', 'NEPTUNE', 'PLUTO', 'SUN', 'MOON', 'TRACKFIELD'] |
778 + | if isinstance(phasecenter, str) and phasecenter.strip().upper() in ephem_sources: |
762 779 | _ephemsrcname = phasecenter |
763 - | |
764 - | # handle image parameters |
765 - | if isinstance(infiles, str) or len(infiles) == 1: |
766 - | _imsize, _cell, _phasecenter = _handle_image_params(imsize, cell, phasecenter, infiles, |
767 - | field, _spw, antenna, scan, intent, |
768 - | _restfreq, pointingcolumn, _ephemsrcname) |
769 - | else: |
770 - | # sort input data using cleanhelper function to get consistent result with older sdimaging |
771 - | o = OldImagerBasedTools() |
772 - | _sorted = o.sort_vis(infiles, _spw, mode, imwidth, field, antenna, scan, intent) |
773 - | sorted_vis, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent = _sorted |
774 - | _imsize, _cell, _phasecenter = _handle_image_params(imsize, cell, phasecenter, sorted_vis, |
775 - | sorted_field, sorted_spw, sorted_antenna, |
776 - | sorted_scan, sorted_intent, |
777 - | _restfreq, pointingcolumn, _ephemsrcname) |
780 + | _imsize, _cell, _phasecenter = _handle_image_params(imsize, cell, phasecenter, sorted_vis, |
781 + | sorted_field, sorted_spw, sorted_antenna, |
782 + | sorted_scan, sorted_intent, sorted_timerange, |
783 + | _restfreq, pointingcolumn, _ephemsrcname) |
778 784 | |
779 785 | # calculate pblimit from minweight |
780 786 | pblimit = _calc_pblimit(minweight) |
781 787 | |
782 788 | ## (2) Set up Input Parameters |
783 789 | ## - List all parameters that you need here |
784 790 | ## - Defaults will be assumed for unspecified parameters |
785 791 | ## - Nearly all parameters are identical to that in the task. Please look at the |
786 792 | ## list of parameters under __init__ using " help ImagerParameters " ) |
787 793 | casalog.post('*** Creating paramList ***', origin=origin) |
788 794 | paramList = ImagerParameters( |
789 795 | # input file name |
790 796 | msname =infiles,#'sdimaging.ms', |
791 797 | # data selection |
792 798 | field=field,#'', |
793 799 | spw=_spw,#'0', |
800 + | timestr=timerange, |
794 801 | antenna=baseline, |
795 802 | scan=scan, |
796 803 | state=intent, |
797 804 | # image parameters |
798 805 | imagename=_outfile,#'try2', |
799 806 | nchan=imnchan,#1024, |
800 807 | start=imstart,#'0', |
801 808 | width=imwidth,#'1', |
802 809 | outframe=outframe, |
803 810 | veltype=veltype, |
871 878 | # set beam size |
872 879 | # TODO: re-define related functions in the new tool framework (sdms?) |
873 880 | imagename = outfile + image_suffix |
874 881 | ms_index = 0 |
875 882 | rep_ms = _get_param(0, infiles) |
876 883 | rep_field = _get_param(0, field) |
877 884 | rep_spw = _get_param(0, _spw) |
878 885 | rep_antenna = _get_param(0, antenna) |
879 886 | rep_scan = _get_param(0, scan) |
880 887 | rep_intent = _get_param(0, intent) |
888 + | rep_timerange = _get_param(0, timerange) |
881 889 | if len(rep_antenna) > 0: |
882 890 | baseline = '{0}&&&'.format(rep_antenna) |
883 891 | else: |
884 892 | baseline = '*&&&' |
885 893 | with open_ms(rep_ms) as ms: |
886 894 | ms.msselect({'baseline': baseline}) |
887 895 | ndx = ms.msselectedindices() |
888 896 | antenna_index = ndx['antenna1'][0] |
889 897 | with open_table(os.path.join(rep_ms, 'ANTENNA')) as tb: |
890 898 | antenna_name = tb.getcell('NAME', antenna_index) |
891 899 | antenna_diameter = tb.getcell('DISH_DIAMETER', antenna_index) |
892 900 | set_beam_size(rep_ms, imagename, |
893 - | rep_field, rep_spw, baseline, rep_scan, rep_intent, |
901 + | rep_field, rep_spw, baseline, rep_scan, rep_intent, rep_timerange, |
894 902 | _ephemsrcname, pointingcolumn, antenna_name, antenna_diameter, |
895 903 | _restfreq, gridfunction, convsupport, truncate, gwidth, jwidth) |
896 904 | |
897 905 | # set brightness unit (CAS-11503) |
898 906 | if len(image_unit) == 0: |
899 907 | image_unit = get_brightness_unit_from_ms(rep_ms) |
900 908 | if len(image_unit) > 0: |
901 909 | with open_ia(imagename) as ia: |
902 910 | casalog.post("Setting brightness unit '%s' to image." % image_unit) |
903 911 | ia.setbrightnessunit(image_unit) |