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