Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/task_sdimaging.py: 5%
801 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-21 07:43 +0000
1# sd task for imaging
2import os
3import re
4import shutil
5import string
7import numpy
9from casatasks import casalog
10from casatools import image, measures, ms, msmetadata, quanta, table
12from . import mslisthelper, sdbeamutil, sdutil
13from .task_tsdimaging import conform_mslist
16@sdutil.sdtask_decorator
17def sdimaging(infiles, outfile, overwrite, field, spw, antenna, scan, intent,
18 mode, nchan, start, width, veltype, outframe,
19 gridfunction, convsupport, truncate, gwidth, jwidth,
20 imsize, cell, phasecenter, projection, ephemsrcname,
21 pointingcolumn, restfreq, stokes, minweight, brightnessunit, clipminmax,
22 # Performances optimization options
23 enablecache, convertfirst,
24 interpolation):
26 casalog.post(
27 'The sdimaging task is deprecated and will be removed in future releases. '
28 'Please use tsdimaging instead.',
29 "WARN")
30 with sdimaging_worker(**locals()) as worker:
31 worker.initialize()
32 worker.execute()
33 worker.finalize()
36def is_string_type(val):
37 """Return True if the argument is string type."""
38 return type(val) in [str, numpy.str_]
41def smart_remove(path):
42 if os.path.isdir(path):
43 shutil.rmtree(path)
44 else:
45 os.remove(path)
48def get_subtable_path(vis, name):
49 return os.path.join(vis, name)
52# functions from cleanhelper
53class cleanhelper_minimal(object):
54 def __init__(self, imtool='', vis='', sort_index=None, usescratch=False, casalog=casalog):
55 if((type(vis)==list) and (len(vis)==1)):
56 vis=vis[0]
57 ####
58 if isinstance(vis, str):
59 vislist = [vis]
60 else:
61 vislist = vis
62 assert not isinstance(imtool, str)
63 self.im = imtool
64 self.vislist = vislist
65 assert sort_index is not None
66 self.sortedvisindx = sort_index
67 self.sortedvislist = [self.vislist[i] for i in sort_index]
69 self.dataspecframe='LSRK'
70 self.usespecframe=''
71 self.inframe=False
72 self._casalog = casalog
74 def convertvf(self,vf,frame,field,restf,veltype='radio'):
75 """
76 returns doppler(velocity) or frequency in string
77 currently use first rest frequency
78 Assume input vf (velocity or fequency in a string) and
79 output are the same 'frame'.
80 """
81 #pdb.set_trace()
82 docalcf=False
83 #if(frame==''): frame='LSRK'
84 #Use datasepcframe, it is cleanhelper initialized to set
85 #to LSRK
86 if(frame==''): frame=self.dataspecframe
87 qa = quanta()
88 if(qa.quantity(vf)['unit'].find('m/s') > -1):
89 docalcf=True
90 elif(qa.quantity(vf)['unit'].find('Hz') > -1):
91 docalcf=False
92 else:
93 if vf !=0:
94 raise TypeError("Unrecognized unit for the velocity or frequency parameter")
95 ##fldinds=ms.msseltoindex(self.vis, field=field)['field'].tolist()
96 myms = ms()
97 fldinds=myms.msseltoindex(self.vislist[self.sortedvisindx[0]], field=field)['field'].tolist()
98 if(len(fldinds) == 0):
99 fldid0=0
100 else:
101 fldid0=fldinds[0]
102 if restf=='':
103 #tb.open(self.vis+'/FIELD')
104 fldtab=get_subtable_path(self.vislist[self.sortedvisindx[0]],'FIELD')
105 tb = table()
106 tb.open(fldtab)
107 nfld = tb.nrows()
108 if nfld >= fldid0:
109 srcid=tb.getcell('SOURCE_ID',fldid0)
110 else:
111 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
112 "no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
113 tb.close()
114 # SOUECE_ID in FIELD table = -1 if no SOURCE table
115 if srcid==-1:
116 raise TypeError("Rest frequency info is not supplied")
117 #tb.open(self.vis+'/SOURCE')
118 sourcetab=get_subtable_path(self.vislist[self.sortedvisindx[0]], 'SOURCE')
119 tb.open(sourcetab)
120 tb2=tb.query('SOURCE_ID==%s' % srcid)
121 tb.close()
122 nsrc = tb2.nrows()
123 if nsrc > 0:
124 rfreq=tb2.getcell('REST_FREQUENCY',0)
125 else:
126 raise TypeError( "Cannot set REST_FREQUENCY from the data: "+
127 " no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
128 tb2.close()
129 if(rfreq<=0):
130 raise TypeError("Rest frequency does not seems to be properly set, check the data")
131 else:
132 if type(restf)==str: restf=[restf]
133 if(qa.quantity(restf[0])['unit'].find('Hz') > -1):
134 rfreq=[qa.convert(qa.quantity(restf[0]),'Hz')['value']]
135 #print("using user input rest freq=",rfreq)
136 else:
137 raise TypeError("Unrecognized unit or type for restfreq")
138 if(vf==0):
139 # assume just want to get a restfrequecy from the data
140 ret=str(rfreq[0])+'Hz'
141 else:
142 me = measures()
143 if(docalcf):
144 dop=me.doppler(veltype, qa.quantity(vf))
145 rvf=me.tofrequency(frame, dop, qa.quantity(rfreq[0],'Hz'))
146 else:
147 frq=me.frequency(frame, qa.quantity(vf))
148 rvf=me.todoppler(veltype, frq, qa.quantity(rfreq[0],'Hz'))
149 ret=str(rvf['m0']['value'])+rvf['m0']['unit']
150 return ret
153 def setChannelizeDefault(self,mode,spw,field,nchan,start,width,frame,veltype,phasec, restf,obstime=''):
154 """
155 Determine appropriate values for channelization
156 parameters when default values are used
157 for mode='velocity' or 'frequency' or 'channel'
158 This makes use of ms.cvelfreqs.
159 """
160 ###############
161 # for debugging
162 ###############
163 debug=False
164 ###############
165 spectable=get_subtable_path(self.vislist[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
166 tb = table()
167 tb.open(spectable)
168 chanfreqscol=tb.getvarcol('CHAN_FREQ')
169 chanwidcol=tb.getvarcol('CHAN_WIDTH')
170 spwframe=tb.getcol('MEAS_FREQ_REF');
171 tb.close()
172 # first parse spw parameter:
173 # use MSSelect if possible
174 if len(self.sortedvislist) > 0:
175 invis = self.sortedvislist[0]
176 inspw = self.vislist.index(self.sortedvislist[0])
177 else:
178 invis = self.vislist[0]
179 inspw = 0
180 myms = ms()
181 myms.open(invis)
182 if type(spw)==list:
183 spw=spw[inspw]
184 if spw in ('-1', '*', '', ' '):
185 spw='*'
186 if field=='':
187 field='*'
188 mssel=myms.msseltoindex(vis=invis, spw=spw, field=field)
189 selspw=mssel['spw']
190 selfield=mssel['field']
191 chaninds=mssel['channel'].tolist()
192 chanst0 = chaninds[0][1]
194 # frame
195 spw0=selspw[0]
196 chanfreqs=chanfreqscol['r'+str(spw0+1)].transpose()[0]
197 chanres = chanwidcol['r'+str(spw0+1)].transpose()[0]
199 # ascending or desending data frequencies?
200 # based on selected first spw's first CHANNEL WIDTH
201 # ==> some MS data may have positive chan width
202 # so changed to look at first two channels of chanfreq (TT)
203 #if chanres[0] < 0:
204 descending = False
205 if len(chanfreqs) > 1 :
206 if chanfreqs[1]-chanfreqs[0] < 0:
207 descending = True
208 else:
209 if chanres[0] < 0:
210 descending = True
212 # set dataspecframe:
213 elspecframe=["REST",
214 "LSRK",
215 "LSRD",
216 "BARY",
217 "GEO",
218 "TOPO",
219 "GALACTO",
220 "LGROUP",
221 "CMB"]
222 self.dataspecframe=elspecframe[spwframe[spw0]];
224 # set usespecframe: user's frame if set, otherwise data's frame
225 if(frame != ''):
226 self.usespecframe=frame
227 self.inframe=True
228 else:
229 self.usespecframe=self.dataspecframe
231 # some start and width default handling
232 if mode!='channel':
233 if width==1:
234 width=''
235 if start==0:
236 start=''
238 #get restfreq
239 if restf=='':
240 fldtab=get_subtable_path(invis,'FIELD')
241 tb.open(fldtab)
242 nfld=tb.nrows()
243 try:
244 if nfld >= selfield[0]:
245 srcid=tb.getcell('SOURCE_ID',selfield[0])
246 else:
247 if mode=='velocity':
248 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
249 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
250 finally:
251 tb.close()
252 #SOUECE_ID in FIELD table = -1 if no SOURCE table
253 if srcid==-1:
254 if mode=='velocity':
255 raise TypeError("Rest frequency info is not supplied")
256 try:
257 srctab=get_subtable_path(invis, 'SOURCE')
258 tb.open(srctab)
259 tb2=tb.query('SOURCE_ID==%s' % srcid)
260 nsrc = tb2.nrows()
261 if nsrc > 0 and tb2.iscelldefined('REST_FREQUENCY',0):
262 rfqs = tb2.getcell('REST_FREQUENCY',0)
263 if len(rfqs)>0:
264 restf=str(rfqs[0])+'Hz'
265 else:
266 if mode=='velocity':
267 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
268 "REST_FREQUENCY entry for ID %s in SOURCE table is empty, please supply restfreq" % srcid )
269 else:
270 if mode=='velocity':
271 raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
272 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
273 finally:
274 tb.close()
275 tb2.close()
277 if type(phasec)==list:
278 inphasec=phasec[0]
279 else:
280 inphasec=phasec
281 if type(inphasec)==str and inphasec.isdigit():
282 inphasec=int(inphasec)
283 #if nchan==1:
284 # use data chan freqs
285 # newfreqs=chanfreqs
286 #else:
287 # obstime not included here
288 if debug: print("before ms.cvelfreqs (start,width,nchan)===>",start, width, nchan)
289 try:
290 newfreqs=myms.cvelfreqs(spwids=selspw,fieldids=selfield,mode=mode,nchan=nchan,
291 start=start,width=width,phasec=inphasec, restfreq=restf,
292 outframe=self.usespecframe,veltype=veltype).tolist()
293 except:
294 # ms must be closed here if ms.cvelfreqs failed with an exception
295 myms.close()
296 raise
297 myms.close()
299 #print(newfreqs)
300 descendingnewfreqs=False
301 if len(newfreqs)>1:
302 if newfreqs[1]-newfreqs[0] < 0:
303 descendingnewfreqs=True
306 try:
307 if((nchan in [-1, "-1", "", " "]) or (start in [-1, "-1", "", " "])):
308 frange=im.advisechansel(msname=invis, spwselection=spw, fieldid=selfield[0], getfreqrange=True)
309 startchan=0
310 endchan=len(newfreqs)-1
311 if(descendingnewfreqs):
312 startchan=numpy.min(numpy.where(frange['freqend'] < numpy.array(newfreqs)))
313 endchan=numpy.min(numpy.where(frange['freqstart'] < numpy.array(newfreqs)))
314 else:
315 startchan=numpy.max(numpy.where(frange['freqstart'] > numpy.array(newfreqs)))
316 endchan=numpy.max(numpy.where(frange['freqend'] > numpy.array(newfreqs)))
317 if((start not in [-1, "-1", "", " "]) and (mode=="channel")):
318 startchan=start
319 if(nchan not in [-1, "-1", "", " "]):
320 endchan=startchan+nchan-1
321 newfreqs=(numpy.array(newfreqs)[startchan:endchan]).tolist()
322 except:
323 pass
324 if debug: print("Mode, Start, width after cvelfreqs =",mode, start,width )
325 if type(newfreqs)==list and len(newfreqs) ==0:
326 raise TypeError( "Output frequency grid cannot be calculated: " +
327 " please check start and width parameters" )
328 if debug:
329 if len(newfreqs)>1:
330 print("FRAME=",self.usespecframe)
331 print("newfreqs[0]===>",newfreqs[0])
332 print("newfreqs[1]===>",newfreqs[1])
333 print("newfreqs[-1]===>",newfreqs[-1])
334 print("len(newfreqs)===>",len(newfreqs))
335 else:
336 print("newfreqs=",newfreqs)
338 # set output number of channels
339 if nchan ==1:
340 retnchan=1
341 else:
342 if len(newfreqs)>1:
343 retnchan=len(newfreqs)
344 else:
345 retnchan=nchan
346 newfreqs=chanfreqs.tolist()
348 # set start parameter
349 # first analyze data order etc
350 reverse=False
351 negativew=False
352 if descending:
353 # channel mode case (width always >0)
354 if width!="" and (type(width)==int or type(width)==float):
355 if descendingnewfreqs:
356 reverse=False
357 else:
358 reverse=True
359 elif width=="": #default width
360 if descendingnewfreqs and mode=="frequency":
361 reverse=False
362 else:
363 reverse=True
365 elif type(width)==str:
366 if width.lstrip().find('-')==0:
367 negativew=True
368 if descendingnewfreqs:
369 if negativew:
370 reverse=False
371 else:
372 reverse=True
373 else:
374 if negativew:
375 reverse=True
376 else:
377 reverse=False
378 else: #ascending data
379 # depends on sign of width only
380 # with CAS-3117 latest change(rev.15179), velocity start
381 # means lowest velocity for default width
382 if width=="" and mode=="velocity": #default width
383 # ms.cvelfreqs returns correct order so no reversing
384 reverse=False
385 elif type(width)==str:
386 if width.lstrip().find('-')==0:
387 reverse=True
388 else:
389 reverse=False
391 if reverse:
392 newfreqs.reverse()
393 #if (start!="" and mode=='channel') or \
394 # (start!="" and type(start)!=int and mode!='channel'):
395 # for now to avoid inconsistency later in imagecoordinates2 call
396 # user's start parameter is preserved for channel mode only.
397 # (i.e. the current code may adjust start parameter for other modes but
398 # this probably needs to be changed, especially for multiple ms handling.)
399 if ((start not in [-1, "", " "]) and mode=='channel'):
400 retstart=start
401 else:
402 # default cases
403 if mode=="frequency":
404 retstart=str(newfreqs[0])+'Hz'
405 elif mode=="velocity":
406 #startfreq=str(newfreqs[-1])+'Hz'
407 startfreq=(str(max(newfreqs))+'Hz') if(start=="") else (str(newfreqs[-1])+'Hz')
408 retstart=self.convertvf(startfreq,frame,field,restf,veltype)
409 elif mode=="channel":
410 # default start case, use channel selection from spw
411 retstart=chanst0
413 # set width parameter
414 if width!="":
415 retwidth=width
416 else:
417 if nchan==1:
418 finc = chanres[0]
419 else:
420 finc = newfreqs[1]-newfreqs[0]
421 if debug: print("finc(newfreqs1-newfreqs0)=",finc)
422 if mode=="frequency":
423 # It seems that this is no longer necessary... TT 2013-08-12
424 #if descendingnewfreqs:
425 # finc = -finc
426 retwidth=str(finc)+'Hz'
427 elif mode=="velocity":
428 # for default width assume it is vel<0 (incresing in freq)
429 if descendingnewfreqs:
430 ind1=-2
431 ind0=-1
432 else:
433 ind1=-1
434 ind0=-2
435 v1 = self.convertvf(str(newfreqs[ind1])+'Hz',frame,field,restf,veltype=veltype)
436 v0 = self.convertvf(str(newfreqs[ind0])+'Hz',frame,field,restf,veltype=veltype)
437 ##v1 = self.convertvf(str(newfreqs[-1])+'Hz',frame,field,restf,veltype=veltype)
438 ##v0 = self.convertvf(str(newfreqs[-2])+'Hz',frame,field,restf,veltype=veltype)
439 #v1 = self.convertvf(str(newfreqs[1])+'Hz',frame,field,restf,veltype=veltype)
440 #v0 = self.convertvf(str(newfreqs[0])+'Hz',frame,field,restf,veltype=veltype)
441 qa = quanta()
442 if(qa.lt(v0, v1) and start==""):
443 ###user used "" as start make sure step is +ve in vel as start is min vel possible for freqs selected
444 retwidth=qa.tos(qa.sub(v1, v0))
445 else:
446 retwidth = qa.tos(qa.sub(v0, v1))
447 else:
448 retwidth=1
449 if debug: print("setChan retwidth=",retwidth)
450 return retnchan, retstart, retwidth
453class sdimaging_worker(sdutil.sdtask_template_imaging):
454 def __init__(self, **kwargs):
455 super(sdimaging_worker, self).__init__(**kwargs)
456 self.imager_param = {}
457 self.sorted_idx = []
458 self.image_unit = ""
460 def __exit__(self, exc_type, exc_val, exc_tb):
461 self.__del__()
462 if exc_type:
463 return False
464 else:
465 return True
467 def parameter_check(self):
468 # outfile check
469 sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile,
470 'im',
471 self.overwrite)
472 sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile + '.weight',
473 'im',
474 self.overwrite)
475 # fix spw
476 if type(self.spw) == str:
477 self.spw = self.__format_spw_string(self.spw)
478 # check unit of start and width
479 # fix default
480 if self.mode == 'channel':
481 if self.start == '':
482 self.start = 0
483 if self.width == '':
484 self.width = 1
485 else:
486 if self.start == 0:
487 self.start = ''
488 if self.width == 1:
489 self.width = ''
490 # fix unit
491 if self.mode == 'frequency':
492 myunit = 'Hz'
493 elif self.mode == 'velocity':
494 myunit = 'km/s'
495 else: # channel
496 myunit = ''
498 for name in ['start', 'width']:
499 param = getattr(self, name)
500 new_param = self.__format_quantum_unit(param, myunit)
501 if new_param is None:
502 raise ValueError("Invalid unit for %s in mode %s: %s" %
503 (name, self.mode, param))
504 setattr(self, name, new_param)
506 casalog.post("mode='%s': start=%s, width=%s, nchan=%d" %
507 (self.mode, self.start, self.width, self.nchan))
509 # check length of selection parameters
510 if is_string_type(self.infiles):
511 nfile = 1
512 self.infiles = [self.infiles]
513 else:
514 nfile = len(self.infiles)
516 for name in ['field', 'spw', 'antenna', 'scanno']:
517 param = getattr(self, name)
518 if not self.__check_selection_length(param, nfile):
519 raise ValueError("Length of %s != infiles." % (name))
520 # set convsupport default
521 if self.convsupport >= 0 and self.gridfunction.upper() != 'SF':
522 casalog.post(
523 f"user defined convsupport is ignored for {self.gridfunction} kernel",
524 priority='WARN')
525 self.convsupport = -1
527 def __format_spw_string(self, spw):
528 """Return formatted spw selection string which is accepted by imager."""
529 if type(spw) != str:
530 raise ValueError("The parameter should be string.")
531 if spw.strip() == '*':
532 spw = ''
533 # WORKAROUND for CAS-6422, i.e., ":X~Y" fails while "*:X~Y" works.
534 if spw.startswith(":"):
535 spw = '*' + spw
536 return spw
538 def __format_quantum_unit(self, data, unit):
539 """Format quantity data.
541 Returns False if data has an unit which in not a variation of
542 input unit.
543 Otherwise, returns input data as a quantum string. The input
544 unit is added to the return value if no unit is in data.
545 """
546 my_qa = quanta()
547 if data == '' or my_qa.compare(data, unit):
548 return data
549 if my_qa.getunit(data) == '':
550 casalog.post("No unit specified. Using '%s'" % unit)
551 return '%f%s' % (data, unit)
552 return None
554 def __check_selection_length(self, data, nfile):
555 """Check length of the data.
557 Returns true if data is either a string, an array with length
558 1 or nfile
559 """
560 if not is_string_type(data) and len(data) not in [1, nfile]:
561 return False
562 return True
564 def get_selection_param_for_ms(self, fileid, param):
565 """Return valid selection string for a certain ms.
567 Arguments
568 fileid : file idx in infiles list
569 param : string (array) selection value
570 """
571 if is_string_type(param):
572 return param
573 elif len(param) == 1:
574 return param[0]
575 else:
576 return param[fileid]
578 def get_selection_idx_for_ms(self, file_idx):
579 """Return a dictionary of selection indices for i-th MS in infiles.
581 Argument: file idx in infiles list
582 """
583 if file_idx < len(self.infiles) and file_idx > -1:
584 vis = self.infiles[file_idx]
585 field = self.get_selection_param_for_ms(file_idx, self.field)
586 spw = self.get_selection_param_for_ms(file_idx, self.spw)
587 spw = self.__format_spw_string(spw)
588 antenna = self.get_selection_param_for_ms(file_idx, self.antenna)
589 if antenna == -1:
590 antenna = ''
591 scan = self.get_selection_param_for_ms(file_idx, self.scanno)
592 intent = self.get_selection_param_for_ms(file_idx, self.intent)
593 my_ms = ms()
594 sel_ids = my_ms.msseltoindex(vis=vis, spw=spw, field=field,
595 baseline=antenna, scan=scan)
596 fieldid = list(sel_ids['field']) if len(sel_ids['field']) > 0 else -1
597 baseline = self.format_ac_baseline(sel_ids['antenna1'])
598 scanid = list(sel_ids['scan']) if len(sel_ids['scan']) > 0 else ""
599 # SPW (need to get a list of valid spws instead of -1)
600 if len(sel_ids['channel']) > 0:
601 spwid = [chanarr[0] for chanarr in sel_ids['channel']]
602 elif spw == "": # No spw selection
603 my_ms.open(vis)
604 try:
605 spwinfo = my_ms.getspectralwindowinfo()
606 except Exception:
607 raise
608 finally:
609 my_ms.close()
611 spwid = [int(idx) for idx in spwinfo.keys()]
612 else:
613 raise RuntimeError("Invalid spw selction, %s ,for MS %d" % (str(spw), file_idx))
615 return {'field': fieldid, 'spw': spwid, 'baseline': baseline, 'scan': scanid,
616 'intent': intent, 'antenna1': sel_ids['antenna1']}
617 else:
618 raise ValueError("Invalid file index, %d" % file_idx)
620 def format_ac_baseline(self, in_antenna):
621 """Format auto-correlation baseline string from antenna idx list."""
622 # exact match string
623 if is_string_type(in_antenna):
624 # return sdutil.convert_antenna_spec_autocorr(in_antenna)
625 if (len(in_antenna) != 0) and (in_antenna.find('&') == -1) \
626 and (in_antenna.find(';') == -1):
627 in_antenna = + '&&&'
628 return in_antenna
629 # single integer -> list of int
630 if type(in_antenna) == int:
631 if in_antenna >= 0:
632 in_antenna = [in_antenna]
633 else:
634 return -1
635 # format auto-corr string from antenna idices.
636 baseline = ''
637 for idx in in_antenna:
638 if len(baseline) > 0:
639 baseline += ';'
640 if idx >= 0:
641 baseline += (str(idx) + '&&&')
642 return baseline
644 def compile(self):
645 # imaging mode
646 self.imager_param['mode'] = self.mode
648 # Work on selection of the first table in sorted list
649 # to get default restfreq and outframe
650 # chronological sort
651 sorted_vislist, sorted_timelist = mslisthelper.sort_mslist(self.infiles)
652 self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist]
653 mslisthelper.report_sort_result(sorted_vislist, sorted_timelist, self.sorted_idx, mycasalog=casalog)
654 # conform MS
655 conform_mslist(sorted_vislist, ignore_columns=[])
656 selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0])
657 self.__update_subtable_name(self.infiles[self.sorted_idx[0]])
658 # field
659 fieldid = selection_ids['field'][0] if type(
660 selection_ids['field']) != int else selection_ids['field']
661 sourceid = -1
662 self.open_table(self.field_table)
663 source_ids = self.table.getcol('SOURCE_ID')
664 self.close_table()
665 if self.field == '' or fieldid == -1:
666 sourceid = source_ids[0]
667 elif fieldid >= 0 and fieldid < len(source_ids):
668 sourceid = source_ids[fieldid]
669 else:
670 raise ValueError("No valid field in the first MS.")
672 # restfreq
673 if self.restfreq == '' and self.source_table != '':
674 self.open_table(self.source_table)
675 source_ids = self.table.getcol('SOURCE_ID')
676 for i in range(self.table.nrows()):
677 if sourceid == source_ids[i] \
678 and self.table.iscelldefined('REST_FREQUENCY', i) \
679 and (selection_ids['spw'] == -1 or
680 self.table.getcell('SPECTRAL_WINDOW_ID', i) in selection_ids['spw']):
681 rf = self.table.getcell('REST_FREQUENCY', i)
682 if len(rf) > 0:
683 self.restfreq = self.table.getcell('REST_FREQUENCY', i)[0]
684 break
685 self.close_table()
686 casalog.post("restfreq set to %s" % self.restfreq, "INFO")
687 # REST_FREQUENCY column is optional (need retry if not exists)
688 self.imager_param['restfreq'] = self.restfreq
690 #
691 # spw (define representative spw id = spwid_ref)
692 spwid_ref = selection_ids['spw'][0] \
693 if type(selection_ids['spw']) != int else selection_ids['spw']
694 # Invalid spw selection should have handled at msselectiontoindex().
695 # -1 means all spw are selected.
696 self.open_table(self.spw_table)
697 if spwid_ref < 0:
698 for id in range(self.table.nrows()):
699 if self.table.getcell('NUM_CHAN', id) > 0:
700 spwid_ref = id
701 break
702 if spwid_ref < 0:
703 self.close_table()
704 msg = 'No valid spw id exists in the first table'
705 raise ValueError(msg)
706 self.allchannels = self.table.getcell('NUM_CHAN', spwid_ref)
707 # freq_chan0 = self.table.getcell('CHAN_FREQ', spwid_ref)[0]
708 # freq_inc0 = self.table.getcell('CHAN_WIDTH', spwid_ref)[0]
709 # in case rest frequency is not defined yet.
710 if self.restfreq == '':
711 self.restfreq = '%fHz' % self.table.getcell('CHAN_FREQ', spwid_ref).mean()
712 self.imager_param['restfreq'] = self.restfreq
713 casalog.post("Using mean freq of spw %d as restfreq: %s" %
714 (spwid_ref, self.restfreq), "INFO")
715 self.close_table()
716 self.imager_param['spw'] = -1 # spwid_ref
718 # outframe (force using the current frame)
719 self.imager_param['outframe'] = self.outframe
720 if self.outframe == '':
721 if len(self.infiles) > 1:
722 # The default will be 'LSRK'
723 casalog.post("Multiple MS inputs. The default outframe is set to 'LSRK'")
724 self.imager_param['outframe'] = 'LSRK'
725 else:
726 # get from MS
727 my_ms = ms()
728 my_ms.open(self.infiles[0])
729 spwinfo = my_ms.getspectralwindowinfo()
730 my_ms.close()
731 del my_ms
732 for key, spwval in spwinfo.items():
733 if spwval['SpectralWindowId'] == spwid_ref:
734 self.imager_param['outframe'] = spwval['Frame']
735 casalog.post(
736 f"Using frequency frame of MS, '{self.imager_param['outframe']}'")
737 break
738 if self.imager_param['outframe'] == '':
739 raise Exception("Internal error of getting frequency frame of spw=%d." % spwid_ref)
740 else:
741 casalog.post(
742 f"Using frequency frame defined by user, '{self.imager_param['outframe']}'")
744 # brightness unit
745 if len(self.brightnessunit) > 0:
746 if self.brightnessunit.lower() == 'k':
747 self.image_unit = 'K'
748 elif self.brightnessunit.lower() == 'jy/beam':
749 self.image_unit = 'Jy/beam'
750 else:
751 raise ValueError("Invalid brightness unit, %s" % self.brightnessunit)
753# # antenna
754# in_antenna = self.antenna # backup for future use
755# if type(self.antenna)==int:
756# if self.antenna >= 0:
757# self.antenna=str(self.antenna)+'&&&'
758# else:
759# if (len(self.antenna) != 0) and (self.antenna.find('&') == -1) \
760# and (self.antenna.find(';')==-1):
761# self.antenna = self.antenna + '&&&'
763 def _configure_map_property(self):
764 selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0])
766 # stokes
767 if self.stokes == '':
768 self.stokes = 'I'
769 self.imager_param['stokes'] = self.stokes
771 # gridfunction
773 # outfile
774 if os.path.exists(self.outfile) and self.overwrite:
775 smart_remove(self.outfile)
776 if os.path.exists(self.outfile + '.weight') and self.overwrite:
777 smart_remove(self.outfile + '.weight')
779 # cell
780 cell = self.cell
781 if cell == '' or cell[0] == '':
782 # Calc PB
783 grid_factor = 3.
784 casalog.post(
785 "The cell size will be calculated using PB size of antennas in the first MS")
786 qPB = self._calc_PB(selection_ids['antenna1'])
787 cell = '%f%s' % (qPB['value'] / grid_factor, qPB['unit'])
788 casalog.post("Using cell size = PB/%4.2F = %s" % (grid_factor, cell))
790 (cellx, celly) = sdutil.get_cellx_celly(cell, unit='arcmin')
791 self.imager_param['cellx'] = cellx
792 self.imager_param['celly'] = celly
794 # Calculate Pointing center and extent (if necessary)
795 # return a dictionary with keys 'center', 'width', 'height'
796 imsize = sdutil.to_list(self.imsize, int) or \
797 sdutil.to_list(self.imsize, numpy.integer)
798 if imsize is None:
799 imsize = self.imsize if hasattr(self.imsize, '__iter__') else [self.imsize]
800 imsize = [int(numpy.ceil(v)) for v in imsize]
801 casalog.post(
802 "imsize is not integers. force converting to integer pixel numbers.",
803 priority="WARN")
804 casalog.post("rounded-up imsize: %s --> %s" % (str(self.imsize), str(imsize)))
806 phasecenter = self.phasecenter
807 if self.phasecenter == "" or \
808 len(imsize) == 0 or imsize[0] < 1:
809 map_param = self._get_pointing_extent()
810 # imsize
811 if len(imsize) == 0 or imsize[0] < 1:
812 imsize = self._get_imsize(map_param['width'], map_param['height'], cellx, celly)
813 if self.phasecenter != "":
814 casalog.post(
815 "You defined phasecenter but not imsize. "
816 "The image will cover as wide area as pointing in MS extends, "
817 "but be centered at phasecenter. "
818 "This could result in a strange image if your phasecenter is "
819 "apart from the center of pointings",
820 priority='WARN')
821 if imsize[0] > 1024 or imsize[1] > 1024:
822 casalog.post(
823 "The calculated image pixel number is larger than 1024. "
824 "It could take time to generate the image depending on "
825 "your computer resource. Please wait...",
826 priority='WARN')
828 # phasecenter
829 # if empty, it should be determined here...
830 if self.phasecenter == "":
831 phasecenter = map_param['center']
833 # imsize
834 (nx, ny) = sdutil.get_nx_ny(imsize)
835 self.imager_param['nx'] = nx
836 self.imager_param['ny'] = ny
838 # phasecenter
839 self.imager_param['phasecenter'] = phasecenter
841 self.imager_param['movingsource'] = self.ephemsrcname
843 # channel map
844 sorted_vislist, _ = mslisthelper.sort_mslist(self.infiles)
845 self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist]
846 imhelper = cleanhelper_minimal(self.imager, self.infiles, sort_index=self.sorted_idx, casalog=casalog)
847 spwsel = str(',').join([str(spwid) for spwid in selection_ids['spw']])
848 srestf = self.imager_param['restfreq'] \
849 if is_string_type(self.imager_param['restfreq']) \
850 else "%fHz" % self.imager_param['restfreq']
851 (imnchan, imstart, imwidth) = imhelper.setChannelizeDefault(
852 self.mode, spwsel, self.field, self.nchan,
853 self.start, self.width, self.imager_param['outframe'],
854 self.veltype, self.imager_param['phasecenter'], srestf)
855 del imhelper
857 # start and width
858 if self.mode == 'velocity':
859 startval = [self.imager_param['outframe'], imstart]
860 widthval = imwidth
861 elif self.mode == 'frequency':
862 startval = [self.imager_param['outframe'], imstart]
863 widthval = imwidth
864 else: # self.mode==channel
865 startval = int(self.start)
866 widthval = int(self.width)
868 if self.nchan < 0:
869 self.nchan = self.allchannels
870 self.imager_param['start'] = startval
871 self.imager_param['step'] = widthval
872 self.imager_param['nchan'] = imnchan # self.nchan
874 self.imager_param['projection'] = self.projection
876 def execute(self):
877 # imaging
878 casalog.post("Start imaging...", "INFO")
879 if len(self.infiles) == 1:
880 self.open_imager(self.infiles[0])
881 selection_ids = self.get_selection_idx_for_ms(0)
882 spwsel = self.get_selection_param_for_ms(0, self.spw)
883 if spwsel.strip() in ['', '*']:
884 spwsel = selection_ids['spw']
885 # TODO: channel selection based on spw
886 ok = self.imager.selectvis(field=selection_ids['field'],
887 # spw=selection_ids['spw'],
888 spw=spwsel,
889 nchan=-1, start=0, step=1,
890 baseline=selection_ids['baseline'],
891 scan=selection_ids['scan'],
892 intent=selection_ids['intent'])
893 if not ok:
894 raise ValueError("Selection is empty: you may want to review this MS selection")
895 else:
896 self.close_imager()
897 # self.sorted_idx.reverse()
898 for idx in self.sorted_idx.__reversed__():
899 name = self.infiles[idx]
900 selection_ids = self.get_selection_idx_for_ms(idx)
901 spwsel = self.get_selection_param_for_ms(idx, self.spw)
902 if spwsel.strip() in ['', '*']:
903 spwsel = selection_ids['spw']
904 # TODO: channel selection based on spw
905 self.imager.selectvis(vis=name, field=selection_ids['field'],
906 # spw=selection_ids['spw'],
907 spw=spwsel,
908 nchan=-1, start=0, step=1,
909 baseline=selection_ids['baseline'],
910 scan=selection_ids['scan'],
911 intent=selection_ids['intent'])
912 # need to do this
913 self.is_imager_opened = True
915 # it should be called after infiles are registered to imager
916 self._configure_map_property()
918 casalog.post(f"Using phasecenter {self.imager_param['phasecenter']}", "INFO")
920 self.imager.defineimage(**self.imager_param) # self.__get_param())
921 self.imager.setoptions(ftmachine='sd', gridfunction=self.gridfunction, freqinterp=self.interpolation)
922 self.imager.setsdoptions(
923 pointingcolumntouse=self.pointingcolumn,
924 convsupport=self.convsupport,
925 truncate=self.truncate,
926 gwidth=self.gwidth,
927 jwidth=self.jwidth,
928 minweight=0.,
929 clipminmax=self.clipminmax,
930 enablecache=self.enablecache,
931 convertfirst=self.convertfirst
932 )
934 # Create images
935 imgtype_suffix = {'singledish': '', 'coverage': '.weight'}
936 for img_type, img_suffix in imgtype_suffix.items():
937 img_file = self.outfile + img_suffix
938 msg_fmt = string.Template(f"$state {img_type} image {img_file}")
939 casalog.post(msg_fmt.substitute(state="Generating"), "INFO")
940 self.imager.makeimage(type=img_type, image=img_file)
941 if not os.path.exists(img_file):
942 raise RuntimeError(msg_fmt.substitute(state="Failed to generate"))
943 casalog.post(msg_fmt.substitute(state="Generated"), "INFO")
944 # Close imager
945 self.close_imager()
947 # Convert output images to proper output frame and set brightness unit (if necessary)
948 my_ia = image()
949 my_ia.open(self.outfile)
950 csys = my_ia.coordsys()
951 csys.setconversiontype(spectral=csys.referencecode('spectra')[0])
952 my_ia.setcoordsys(csys.torecord())
953 if len(self.image_unit) > 0:
954 casalog.post("Setting brightness unit '%s' to image." % self.image_unit)
955 my_ia.setbrightnessunit(self.image_unit)
956 csys.done()
957 my_ia.close()
959 # CAS-12984 set brightness unit for weight image to ''
960 weightfile = self.outfile + imgtype_suffix['coverage']
961 my_ia.open(weightfile)
962 my_ia.setbrightnessunit('')
963 my_ia.close()
965 # Mask image pixels whose weight are smaller than minweight.
966 # Weight image should have 0 weight for pixels below < minweight
967 casalog.post("Start masking the map using minweight = %f" %
968 self.minweight, "INFO")
969 my_ia.open(weightfile)
970 try:
971 stat = my_ia.statistics(mask="'" + weightfile + "' > 0.0", robust=True)
972 valid_pixels = stat['npts']
973 except RuntimeError as e:
974 if str(e).find('No valid data found.') >= 0:
975 valid_pixels = [0]
976 else:
977 raise
978 if len(valid_pixels) == 0 or valid_pixels[0] == 0:
979 my_ia.close()
980 casalog.post(
981 "All pixels weight zero. This indicates no data in MS is in image area. "
982 "Mask will not be set. Please check your image parameters.",
983 priority="WARN")
984 return
985 median_weight = stat['median'][0]
986 weight_threshold = median_weight * self.minweight
987 casalog.post("Median of weight in the map is %f" % median_weight,
988 "INFO")
989 casalog.post("Pixels in map with weight <= median(weight)*minweight = %f will be masked." %
990 (weight_threshold), "INFO")
991 # Leaving the original logic to calculate the number of masked pixels via
992 # product of median of and min_weight (which i don't understand the logic)
993 # if one wanted to find how many pixel were masked one could easily count the
994 # number of pixels set to false
995 # e.g after masking self.outfile below one could just do this
996 # nmasked_pixels=tb.calc(
997 # '[select from "'+self.outfile+'"/mask0'+'" giving [nfalse(PagedArray)]]')
998 my_tb = table()
999 nmask_pixels = 0
1000 nchan = stat['trc'][3] + 1
1001 casalog.filter('ERROR') # hide the useless message of tb.calc
1003 # doing it by channel to make sure it does not go out of memory
1004 # tab.calc try to load the whole chunk in ram
1005 for k in range(nchan):
1006 nmask_pixels += my_tb.calc(
1007 '[select from "' + weightfile + '" giving [ntrue(map[,,,' +
1008 str(k) + '] <=' + str(median_weight * self.minweight) + ')]]')['0'][0]
1009 casalog.filter() # set logging back to normal
1011 casalog.filter() # set logging back to normal
1012 imsize = numpy.prod(my_ia.shape())
1013 my_ia.close()
1014 # Modify default mask
1015 my_ia.open(self.outfile)
1016 my_ia.calcmask("'%s'>%f" % (weightfile, weight_threshold), asdefault=True)
1017 my_ia.close()
1018 masked_fraction = 100. * (1. - (imsize - nmask_pixels) / float(valid_pixels[0]))
1019 casalog.post("This amounts to %5.1f %% of the area with nonzero weight." %
1020 (masked_fraction), "INFO")
1021 casalog.post(
1022 f"The weight image '{weightfile}' is returned by this task, "
1023 "if the user wishes to assess the results in detail.",
1024 priority="INFO")
1026 # Calculate theoretical beam size
1027 casalog.post("Calculating image beam size.")
1028 if self.gridfunction.upper() not in ['SF']:
1029 casalog.post(
1030 f"Beam size definition for '{self.gridfunction}' kernel is experimental.",
1031 priority='WARN')
1032 casalog.post(
1033 "You may want to take careful look at the restoring beam in the image.",
1034 priority='WARN')
1035 my_msmd = msmetadata()
1036 # antenna diameter and blockage
1037 ref_ms_idx = self.sorted_idx[0]
1038 ref_ms_name = self.infiles[ref_ms_idx]
1039 selection_ids = self.get_selection_idx_for_ms(ref_ms_idx)
1040 ant_idx = selection_ids['antenna1']
1041 diameter = self._get_average_antenna_diameter(ant_idx)
1042 my_msmd.open(ref_ms_name)
1043 ant_name = my_msmd.antennanames(ant_idx)
1044 my_msmd.close()
1045 is_alma = False
1046 for name in ant_name:
1047 if name[0:2] in ["PM", "DV", "DA", "CM"]:
1048 is_alma = True
1049 break
1050 blockage = "0.75m" if is_alma else "0.0m" # unknown blockage diameter
1051 # output reference code
1052 my_ia.open(self.outfile)
1053 csys = my_ia.coordsys()
1054 my_ia.close()
1055 outref = csys.referencecode('direction')[0]
1056 cell = list(csys.increment(type='direction', format='s')['string'])
1057 csys.done()
1058 # pointing sampling
1059 ref_ms_spw = self.get_selection_param_for_ms(ref_ms_idx, self.spw)
1060 ref_ms_field = self.get_selection_param_for_ms(ref_ms_idx, self.field)
1061 ref_ms_scan = self.get_selection_param_for_ms(ref_ms_idx, self.scanno)
1062 # xSampling, ySampling, angle = sdutil.get_ms_sampling_arcsec(ref_ms_name, spw=ref_ms_spw,
1063 # antenna=selection_ids['baseline'],
1064 # field=ref_ms_field,
1065 # scan=ref_ms_scan,#timerange='',
1066 # outref=outref)
1068 # obtain sampling interval for beam calculation.
1069 self.open_imager(ref_ms_name)
1070 ok = self.imager.selectvis(field=ref_ms_field,
1071 # spw=selection_ids['spw'],
1072 spw=ref_ms_spw,
1073 nchan=-1, start=0, step=1,
1074 baseline=selection_ids['baseline'],
1075 scan=ref_ms_scan,
1076 intent=selection_ids['intent'])
1077 if len(ant_idx) > 1:
1078 casalog.post("Using only antenna %s to calculate sampling interval" % ant_name[0])
1079 ptg_samp = self.imager.pointingsampling(pattern='raster', ref=outref,
1080 movingsource=self.ephemsrcname,
1081 pointingcolumntouse=self.pointingcolumn,
1082 antenna=('%s&&&' % ant_name[0]))
1083 self.close_imager()
1084 my_qa = quanta()
1085 xSampling, ySampling = my_qa.getvalue(my_qa.convert(ptg_samp['sampling'], 'arcsec'))
1086 angle = my_qa.getvalue(my_qa.convert(ptg_samp['angle'], "deg"))[0]
1088 casalog.post("Detected raster sampling = [{x:f}, {y:f}] arcsec".format(
1089 x=xSampling,
1090 y=ySampling
1091 ))
1093 # handling of failed sampling detection
1094 valid_sampling = True
1095 sampling = [xSampling, ySampling]
1096 if abs(xSampling) < 2.2e-3 or not numpy.isfinite(xSampling):
1097 casalog.post(
1098 f"Invalid sampling={xSampling} arcsec. "
1099 f"Using the value of orthogonal direction={ySampling} arcsec",
1100 priority="WARN")
1101 sampling = [ySampling]
1102 angle = 0.0
1103 valid_sampling = False
1104 if abs(ySampling) < 1.0e-3 or not numpy.isfinite(ySampling):
1105 if valid_sampling:
1106 casalog.post(
1107 f"Invalid sampling={ySampling} arcsec. "
1108 f"Using the value of orthogonal direction={xSampling} arcsec",
1109 priority="WARN")
1110 sampling = [xSampling]
1111 angle = 0.0
1112 valid_sampling = True
1113 # reduce sampling and cell if it's possible
1114 if len(sampling) > 1 and abs(sampling[0] - sampling[1]) <= 0.01 * abs(sampling[0]):
1115 sampling = [sampling[0]]
1116 angle = 0.0
1117 if cell[0] == cell[1]:
1118 cell = [cell[0]]
1119 if valid_sampling:
1120 # actual calculation of beam size
1121 bu = sdbeamutil.TheoreticalBeam()
1122 bu.set_antenna(diameter, blockage)
1123 bu.set_sampling(sampling, "%fdeg" % angle)
1124 bu.set_image_param(cell, self.restfreq, self.gridfunction,
1125 self.convsupport, self.truncate, self.gwidth,
1126 self.jwidth, is_alma)
1127 bu.summary()
1128 imbeam_dict = bu.get_beamsize_image()
1129 casalog.post("Setting image beam: major={major}, minor={minor}, pa={pa}".format(
1130 major=imbeam_dict['major'],
1131 minor=imbeam_dict['minor'],
1132 pa=imbeam_dict['pa']
1133 ))
1134 # set beam size to image
1135 my_ia.open(self.outfile)
1136 my_ia.setrestoringbeam(**imbeam_dict)
1137 my_ia.close()
1138 else:
1139 # BOTH sampling was invalid
1140 casalog.post(
1141 "Could not detect valid raster sampling. "
1142 "Exitting without setting beam size to image",
1143 priority='WARN')
1145 def _calc_PB(self, antenna):
1146 """Calculate the primary beam size of antenna.
1148 Calculate the primary beam size of antenna, using dish diamenter
1149 and rest frequency
1150 Average antenna diamter and reference frequency are adopted for
1151 calculation.
1152 The input argument should be a list of antenna IDs.
1153 """
1154 casalog.post("Calculating Primary beam size:")
1155 # CAS-5410 Use private tools inside task scripts
1156 my_qa = quanta()
1158 pb_factor = 1.175
1159 # Reference frequency
1160 ref_freq = self.restfreq
1161 if type(ref_freq) in [float, numpy.float64]:
1162 ref_freq = my_qa.tos(my_qa.quantity(ref_freq, 'Hz'))
1163 if not my_qa.compare(ref_freq, 'Hz'):
1164 msg = "Could not get the reference frequency. " + \
1165 "Your data does not seem to have valid one in selected field.\n" + \
1166 "PB is not calculated.\n" + \
1167 "Please set restreq or cell manually to generate an image."
1168 raise Exception(msg)
1169 # Antenna diameter
1170 antdiam_ave = self._get_average_antenna_diameter(antenna)
1171 # Calculate PB
1172 wave_length = 0.2997924 / my_qa.convert(my_qa.quantity(ref_freq), 'GHz')['value']
1173 D_m = my_qa.convert(antdiam_ave, 'm')['value']
1174 lambda_D = wave_length / D_m * 3600. * 180 / numpy.pi
1175 PB = my_qa.quantity(pb_factor * lambda_D, 'arcsec')
1176 # Summary
1177 casalog.post("- Antenna diameter: %s m" % D_m)
1178 casalog.post("- Reference Frequency: %s" % ref_freq)
1179 casalog.post("PB size = %5.3f * lambda/D = %s" % (pb_factor, my_qa.tos(PB)))
1180 return PB
1182 def _get_imsize(self, width, height, dx, dy):
1183 casalog.post("Calculating pixel size.")
1184 # CAS-5410 Use private tools inside task scripts
1185 my_qa = quanta()
1186 ny = numpy.ceil((my_qa.convert(height, my_qa.getunit(dy))['value'] /
1187 my_qa.getvalue(dy)))
1188 nx = numpy.ceil((my_qa.convert(width, my_qa.getunit(dx))['value'] /
1189 my_qa.getvalue(dx)))
1190 casalog.post("- Map extent: [%s, %s]" % (my_qa.tos(width), my_qa.tos(height)))
1191 casalog.post("- Cell size: [%s, %s]" % (my_qa.tos(dx), my_qa.tos(dy)))
1192 casalog.post("Image pixel numbers to cover the extent: [%d, %d] (projected)" %
1193 (nx + 1, ny + 1))
1194 return (int(nx + 1), int(ny + 1))
1196 def _get_pointing_extent(self):
1197 # MS selection is ignored. This is not quite right.
1198 casalog.post("Calculating map extent from pointings.")
1199 # CAS-5410 Use private tools inside task scripts
1200 my_qa = quanta()
1201 ret_dict = {}
1203 colname = self.pointingcolumn.upper()
1205 # MSs should be registered to imager
1206 if not self.is_imager_opened:
1207 raise RuntimeError('Internal error: imager should be opened here.')
1209 if self.phasecenter == "":
1210 # defaut is J2000
1211 base_mref = 'J2000'
1212 elif isinstance(self.phasecenter, int) or self.phasecenter.isdigit():
1213 # may be field id
1214 self.open_table(self.field_table)
1215 base_mref = self.table.getcolkeyword('PHASE_DIR', 'MEASINFO')['Ref']
1216 self.close_table()
1217 else:
1218 # may be phasecenter is explicitly specified
1219 # numeric value: 3.14, -.3e1, etc.
1220 numeric_pattern = r'[-+]?([0-9]+(.[0-9]*)?|\.[0-9]+)([eE]-?[0-9])?'
1221 # HMS string: 9:15:29, -9h15m29
1222 hms_pattern = r'[-+]?[0-9]+[:h][0-9]+[:m][0-9.]+s?'
1223 # DMS string: 9.15.29, -9d15m29s
1224 dms_pattern = r'[-+]?[0-9]+[.d][0-9]+[.m][0-9.]+s?'
1225 # composite pattern
1226 pattern = fr'^({numeric_pattern}|{hms_pattern}|{dms_pattern})$'
1227 items = self.phasecenter.split()
1228 base_mref = 'J2000'
1229 for i in items:
1230 s = i.strip()
1231 if re.match(pattern, s) is None:
1232 base_mref = s
1233 break
1235 mapextent = self.imager.mapextent(ref=base_mref, movingsource=self.ephemsrcname,
1236 pointingcolumntouse=colname)
1237 if mapextent['status']:
1238 qheight = my_qa.quantity(mapextent['extent'][1], 'rad')
1239 qwidth = my_qa.quantity(mapextent['extent'][0], 'rad')
1240 qcent0 = my_qa.quantity(mapextent['center'][0], 'rad')
1241 qcent1 = my_qa.quantity(mapextent['center'][1], 'rad')
1242 scenter = '%s %s %s' % (base_mref, my_qa.formxxx(qcent0, 'hms'),
1243 my_qa.formxxx(qcent1, 'dms'))
1245 casalog.post("- Pointing center: %s" % scenter)
1246 casalog.post("- Pointing extent: [%s, %s] (projected)" % (my_qa.tos(qwidth),
1247 my_qa.tos(qheight)))
1248 ret_dict['center'] = scenter
1249 ret_dict['width'] = qwidth
1250 ret_dict['height'] = qheight
1251 else:
1252 casalog.post(
1253 'Failed to derive map extent from the MSs registered to the imager probably '
1254 'due to missing valid data.',
1255 priority='SEVERE')
1256 ret_dict['center'] = ''
1257 ret_dict['width'] = my_qa.quantity(0.0, 'rad')
1258 ret_dict['height'] = my_qa.quantity(0.0, 'rad')
1259 return ret_dict
1261 def _get_x_minmax(self, x):
1262 # assumed the x is in unit of rad.
1263 pi2 = 2. * numpy.pi
1264 x = (x % pi2)
1265 npart = 4
1266 dlon = pi2 / float(npart)
1267 pos = [int(v / dlon) for v in x]
1268 voids = [False for dummy in range(npart)]
1269 for ipos in range(npart):
1270 try:
1271 pos.index(ipos)
1272 except Exception:
1273 voids[ipos] = True
1274 if not any(voids):
1275 raise Exception(
1276 "Failed to find global pointing gap. "
1277 f"The algorithm requires at least 2PI/{npart} of pointing gap")
1278 rot_pos = []
1279 if (not voids[0]) and (not voids[npart - 1]):
1280 gmax = -1
1281 for idx in range(npart - 2, 0, -1):
1282 if voids[idx]:
1283 gmax = idx
1284 break
1285 if gmax < 0:
1286 raise Exception("Failed to detect gap max")
1287 rot_pos = range(gmax + 1, npart)
1288 for idx in range(len(x)):
1289 x[idx] = (x[idx] - pi2) if pos[idx] in rot_pos else x[idx]
1291 return (x.min(), x.max())
1293 def __update_subtable_name(self, msname):
1294 self.open_table(msname)
1295 keys = self.table.getkeywords()
1296 self.close_table()
1297 self.field_table = sdutil.get_subtable_name(keys['FIELD'])
1298 self.spw_table = sdutil.get_subtable_name(keys['SPECTRAL_WINDOW'])
1299 self.source_table = sdutil.get_subtable_name(keys['SOURCE'])
1300 self.antenna_table = sdutil.get_subtable_name(keys['ANTENNA'])
1301 self.polarization_table = sdutil.get_subtable_name(keys['POLARIZATION'])
1302 self.observation_table = sdutil.get_subtable_name(keys['OBSERVATION'])
1303 self.pointing_table = sdutil.get_subtable_name(keys['POINTING'])
1304 self.data_desc_table = sdutil.get_subtable_name(keys['DATA_DESCRIPTION'])
1306 def _get_average_antenna_diameter(self, antenna):
1307 my_qa = quanta()
1308 self.open_table(self.antenna_table)
1309 try:
1310 antdiam_unit = self.table.getcolkeyword('DISH_DIAMETER', 'QuantumUnits')[0]
1311 diams = self.table.getcol('DISH_DIAMETER')
1312 finally:
1313 self.close_table()
1315 if len(antenna) == 0:
1316 antdiam_ave = my_qa.quantity(diams.mean(), antdiam_unit)
1317 else:
1318 d_ave = sum([diams[idx] for idx in antenna]) / float(len(antenna))
1319 antdiam_ave = my_qa.quantity(d_ave, antdiam_unit)
1320 return antdiam_ave