# sd task for imaging
from __future__ import absolute_import
import os
import re
import numpy
import shutil
import string
from casatasks.private.casa_transition import is_CASA6
if is_CASA6:
from casatools import quanta, ms, image, table, msmetadata, measures
from casatasks import casalog
from . import sdutil
from . import sdbeamutil
from . import mslisthelper
from .task_tsdimaging import conform_mslist
else:
from taskinit import casalog
from taskinit import msmdtool as msmetadata
from taskinit import tbtool as table
from taskinit import mstool as ms
from taskinit import iatool as image
from taskinit import qatool as quanta
from taskinit import metool as measures
import sdutil
import sdbeamutil
import recipes.mslisthelper as mslisthelper
from task_tsdimaging import conform_mslist
@sdutil.sdtask_decorator
def sdimaging(infiles, outfile, overwrite, field, spw, antenna, scan, intent,
mode, nchan, start, width, veltype, outframe,
gridfunction, convsupport, truncate, gwidth, jwidth,
imsize, cell, phasecenter, projection, ephemsrcname,
pointingcolumn, restfreq, stokes, minweight, brightnessunit, clipminmax):
with sdutil.sdtask_manager(sdimaging_worker, locals()) as worker:
worker.initialize()
worker.execute()
worker.finalize()
def is_string_type(val):
"""
Returns True if the argument is string type.
"""
return type(val) in [str, numpy.string_]
def smart_remove(path):
if os.path.isdir(path):
shutil.rmtree(path)
else:
os.remove(path)
def get_subtable_path(vis, name):
return os.path.join(vis, name)
# functions from cleanhelper
class cleanhelper_minimal(object):
def __init__(self, imtool='', vis='', sort_index=None, usescratch=False, casalog=casalog):
if((type(vis)==list) and (len(vis)==1)):
vis=vis[0]
####
if isinstance(vis, str):
vislist = [vis]
else:
vislist = vis
assert not isinstance(imtool, str)
self.im = imtool
self.vislist = vislist
assert sort_index is not None
self.sortedvisindx = sort_index
self.sortedvislist = [self.vislist[i] for i in sort_index]
self.dataspecframe='LSRK'
self.usespecframe=''
self.inframe=False
self._casalog = casalog
def convertvf(self,vf,frame,field,restf,veltype='radio'):
"""
returns doppler(velocity) or frequency in string
currently use first rest frequency
Assume input vf (velocity or fequency in a string) and
output are the same 'frame'.
"""
#pdb.set_trace()
docalcf=False
#if(frame==''): frame='LSRK'
#Use datasepcframe, it is cleanhelper initialized to set
#to LSRK
if(frame==''): frame=self.dataspecframe
qa = quanta()
if(qa.quantity(vf)['unit'].find('m/s') > -1):
docalcf=True
elif(qa.quantity(vf)['unit'].find('Hz') > -1):
docalcf=False
else:
if vf !=0:
raise TypeError("Unrecognized unit for the velocity or frequency parameter")
##fldinds=ms.msseltoindex(self.vis, field=field)['field'].tolist()
myms = ms()
fldinds=myms.msseltoindex(self.vislist[self.sortedvisindx[0]], field=field)['field'].tolist()
if(len(fldinds) == 0):
fldid0=0
else:
fldid0=fldinds[0]
if restf=='':
#tb.open(self.vis+'/FIELD')
fldtab=get_subtable_path(self.vislist[self.sortedvisindx[0]],'FIELD')
tb = table()
tb.open(fldtab)
nfld = tb.nrows()
if nfld >= fldid0:
srcid=tb.getcell('SOURCE_ID',fldid0)
else:
raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
"no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
tb.close()
# SOUECE_ID in FIELD table = -1 if no SOURCE table
if srcid==-1:
raise TypeError("Rest frequency info is not supplied")
#tb.open(self.vis+'/SOURCE')
sourcetab=get_subtable_path(self.vislist[self.sortedvisindx[0]], 'SOURCE')
tb.open(sourcetab)
tb2=tb.query('SOURCE_ID==%s' % srcid)
tb.close()
nsrc = tb2.nrows()
if nsrc > 0:
rfreq=tb2.getcell('REST_FREQUENCY',0)
else:
raise TypeError( "Cannot set REST_FREQUENCY from the data: "+
" no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 )
tb2.close()
if(rfreq<=0):
raise TypeError("Rest frequency does not seems to be properly set, check the data")
else:
if type(restf)==str: restf=[restf]
if(qa.quantity(restf[0])['unit'].find('Hz') > -1):
rfreq=[qa.convert(qa.quantity(restf[0]),'Hz')['value']]
#print("using user input rest freq=",rfreq)
else:
raise TypeError("Unrecognized unit or type for restfreq")
if(vf==0):
# assume just want to get a restfrequecy from the data
ret=str(rfreq[0])+'Hz'
else:
me = measures()
if(docalcf):
dop=me.doppler(veltype, qa.quantity(vf))
rvf=me.tofrequency(frame, dop, qa.quantity(rfreq[0],'Hz'))
else:
frq=me.frequency(frame, qa.quantity(vf))
rvf=me.todoppler(veltype, frq, qa.quantity(rfreq[0],'Hz'))
ret=str(rvf['m0']['value'])+rvf['m0']['unit']
return ret
def setChannelizeDefault(self,mode,spw,field,nchan,start,width,frame,veltype,phasec, restf,obstime=''):
"""
Determine appropriate values for channelization
parameters when default values are used
for mode='velocity' or 'frequency' or 'channel'
This makes use of ms.cvelfreqs.
"""
###############
# for debugging
###############
debug=False
###############
spectable=get_subtable_path(self.vislist[self.sortedvisindx[0]], "SPECTRAL_WINDOW")
tb = table()
tb.open(spectable)
chanfreqscol=tb.getvarcol('CHAN_FREQ')
chanwidcol=tb.getvarcol('CHAN_WIDTH')
spwframe=tb.getcol('MEAS_FREQ_REF');
tb.close()
# first parse spw parameter:
# use MSSelect if possible
if len(self.sortedvislist) > 0:
invis = self.sortedvislist[0]
inspw = self.vislist.index(self.sortedvislist[0])
else:
invis = self.vislist[0]
inspw = 0
myms = ms()
myms.open(invis)
if type(spw)==list:
spw=spw[inspw]
if spw in ('-1', '*', '', ' '):
spw='*'
if field=='':
field='*'
mssel=myms.msseltoindex(vis=invis, spw=spw, field=field)
selspw=mssel['spw']
selfield=mssel['field']
chaninds=mssel['channel'].tolist()
chanst0 = chaninds[0][1]
# frame
spw0=selspw[0]
chanfreqs=chanfreqscol['r'+str(spw0+1)].transpose()[0]
chanres = chanwidcol['r'+str(spw0+1)].transpose()[0]
# ascending or desending data frequencies?
# based on selected first spw's first CHANNEL WIDTH
# ==> some MS data may have positive chan width
# so changed to look at first two channels of chanfreq (TT)
#if chanres[0] < 0:
descending = False
if len(chanfreqs) > 1 :
if chanfreqs[1]-chanfreqs[0] < 0:
descending = True
else:
if chanres[0] < 0:
descending = True
# set dataspecframe:
elspecframe=["REST",
"LSRK",
"LSRD",
"BARY",
"GEO",
"TOPO",
"GALACTO",
"LGROUP",
"CMB"]
self.dataspecframe=elspecframe[spwframe[spw0]];
# set usespecframe: user's frame if set, otherwise data's frame
if(frame != ''):
self.usespecframe=frame
self.inframe=True
else:
self.usespecframe=self.dataspecframe
# some start and width default handling
if mode!='channel':
if width==1:
width=''
if start==0:
start=''
#get restfreq
if restf=='':
fldtab=get_subtable_path(invis,'FIELD')
tb.open(fldtab)
nfld=tb.nrows()
try:
if nfld >= selfield[0]:
srcid=tb.getcell('SOURCE_ID',selfield[0])
else:
if mode=='velocity':
raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
"no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
finally:
tb.close()
#SOUECE_ID in FIELD table = -1 if no SOURCE table
if srcid==-1:
if mode=='velocity':
raise TypeError("Rest frequency info is not supplied")
try:
srctab=get_subtable_path(invis, 'SOURCE')
tb.open(srctab)
tb2=tb.query('SOURCE_ID==%s' % srcid)
nsrc = tb2.nrows()
if nsrc > 0 and tb2.iscelldefined('REST_FREQUENCY',0):
rfqs = tb2.getcell('REST_FREQUENCY',0)
if len(rfqs)>0:
restf=str(rfqs[0])+'Hz'
else:
if mode=='velocity':
raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
"REST_FREQUENCY entry for ID %s in SOURCE table is empty, please supply restfreq" % srcid )
else:
if mode=='velocity':
raise TypeError( "Cannot set REST_FREQUENCY from the data: " +
"no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] )
finally:
tb.close()
tb2.close()
if type(phasec)==list:
inphasec=phasec[0]
else:
inphasec=phasec
if type(inphasec)==str and inphasec.isdigit():
inphasec=int(inphasec)
#if nchan==1:
# use data chan freqs
# newfreqs=chanfreqs
#else:
# obstime not included here
if debug: print("before ms.cvelfreqs (start,width,nchan)===>",start, width, nchan)
try:
newfreqs=myms.cvelfreqs(spwids=selspw,fieldids=selfield,mode=mode,nchan=nchan,
start=start,width=width,phasec=inphasec, restfreq=restf,
outframe=self.usespecframe,veltype=veltype).tolist()
except:
# ms must be closed here if ms.cvelfreqs failed with an exception
myms.close()
raise
myms.close()
#print(newfreqs)
descendingnewfreqs=False
if len(newfreqs)>1:
if newfreqs[1]-newfreqs[0] < 0:
descendingnewfreqs=True
try:
if((nchan in [-1, "-1", "", " "]) or (start in [-1, "-1", "", " "])):
frange=im.advisechansel(msname=invis, spwselection=spw, fieldid=selfield[0], getfreqrange=True)
startchan=0
endchan=len(newfreqs)-1
if(descendingnewfreqs):
startchan=numpy.min(numpy.where(frange['freqend'] < numpy.array(newfreqs)))
endchan=numpy.min(numpy.where(frange['freqstart'] < numpy.array(newfreqs)))
else:
startchan=numpy.max(numpy.where(frange['freqstart'] > numpy.array(newfreqs)))
endchan=numpy.max(numpy.where(frange['freqend'] > numpy.array(newfreqs)))
if((start not in [-1, "-1", "", " "]) and (mode=="channel")):
startchan=start
if(nchan not in [-1, "-1", "", " "]):
endchan=startchan+nchan-1
newfreqs=(numpy.array(newfreqs)[startchan:endchan]).tolist()
except:
pass
if debug: print("Mode, Start, width after cvelfreqs =",mode, start,width )
if type(newfreqs)==list and len(newfreqs) ==0:
raise TypeError( "Output frequency grid cannot be calculated: " +
" please check start and width parameters" )
if debug:
if len(newfreqs)>1:
print("FRAME=",self.usespecframe)
print("newfreqs[0]===>",newfreqs[0])
print("newfreqs[1]===>",newfreqs[1])
print("newfreqs[-1]===>",newfreqs[-1])
print("len(newfreqs)===>",len(newfreqs))
else:
print("newfreqs=",newfreqs)
# set output number of channels
if nchan ==1:
retnchan=1
else:
if len(newfreqs)>1:
retnchan=len(newfreqs)
else:
retnchan=nchan
newfreqs=chanfreqs.tolist()
# set start parameter
# first analyze data order etc
reverse=False
negativew=False
if descending:
# channel mode case (width always >0)
if width!="" and (type(width)==int or type(width)==float):
if descendingnewfreqs:
reverse=False
else:
reverse=True
elif width=="": #default width
if descendingnewfreqs and mode=="frequency":
reverse=False
else:
reverse=True
elif type(width)==str:
if width.lstrip().find('-')==0:
negativew=True
if descendingnewfreqs:
if negativew:
reverse=False
else:
reverse=True
else:
if negativew:
reverse=True
else:
reverse=False
else: #ascending data
# depends on sign of width only
# with CAS-3117 latest change(rev.15179), velocity start
# means lowest velocity for default width
if width=="" and mode=="velocity": #default width
# ms.cvelfreqs returns correct order so no reversing
reverse=False
elif type(width)==str:
if width.lstrip().find('-')==0:
reverse=True
else:
reverse=False
if reverse:
newfreqs.reverse()
#if (start!="" and mode=='channel') or \
# (start!="" and type(start)!=int and mode!='channel'):
# for now to avoid inconsistency later in imagecoordinates2 call
# user's start parameter is preserved for channel mode only.
# (i.e. the current code may adjust start parameter for other modes but
# this probably needs to be changed, especially for multiple ms handling.)
if ((start not in [-1, "", " "]) and mode=='channel'):
retstart=start
else:
# default cases
if mode=="frequency":
retstart=str(newfreqs[0])+'Hz'
elif mode=="velocity":
#startfreq=str(newfreqs[-1])+'Hz'
startfreq=(str(max(newfreqs))+'Hz') if(start=="") else (str(newfreqs[-1])+'Hz')
retstart=self.convertvf(startfreq,frame,field,restf,veltype)
elif mode=="channel":
# default start case, use channel selection from spw
retstart=chanst0
# set width parameter
if width!="":
retwidth=width
else:
if nchan==1:
finc = chanres[0]
else:
finc = newfreqs[1]-newfreqs[0]
if debug: print("finc(newfreqs1-newfreqs0)=",finc)
if mode=="frequency":
# It seems that this is no longer necessary... TT 2013-08-12
#if descendingnewfreqs:
# finc = -finc
retwidth=str(finc)+'Hz'
elif mode=="velocity":
# for default width assume it is vel<0 (incresing in freq)
if descendingnewfreqs:
ind1=-2
ind0=-1
else:
ind1=-1
ind0=-2
v1 = self.convertvf(str(newfreqs[ind1])+'Hz',frame,field,restf,veltype=veltype)
v0 = self.convertvf(str(newfreqs[ind0])+'Hz',frame,field,restf,veltype=veltype)
##v1 = self.convertvf(str(newfreqs[-1])+'Hz',frame,field,restf,veltype=veltype)
##v0 = self.convertvf(str(newfreqs[-2])+'Hz',frame,field,restf,veltype=veltype)
#v1 = self.convertvf(str(newfreqs[1])+'Hz',frame,field,restf,veltype=veltype)
#v0 = self.convertvf(str(newfreqs[0])+'Hz',frame,field,restf,veltype=veltype)
qa = quanta()
if(qa.lt(v0, v1) and start==""):
###user used "" as start make sure step is +ve in vel as start is min vel possible for freqs selected
retwidth=qa.tos(qa.sub(v1, v0))
else:
retwidth = qa.tos(qa.sub(v0, v1))
else:
retwidth=1
if debug: print("setChan retwidth=",retwidth)
return retnchan, retstart, retwidth
class sdimaging_worker(sdutil.sdtask_template_imaging):
def __init__(self, **kwargs):
super(sdimaging_worker,self).__init__(**kwargs)
self.imager_param = {}
self.sorted_idx = []
self.image_unit = ""
def parameter_check(self):
# outfile check
sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile,
'im',
self.overwrite)
sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile+'.weight',
'im',
self.overwrite)
# fix spw
if type(self.spw) == str:
self.spw = self.__format_spw_string(self.spw)
# check unit of start and width
# fix default
if self.mode == 'channel':
if self.start == '': self.start = 0
if self.width == '': self.width = 1
else:
if self.start == 0: self.start = ''
if self.width == 1: self.width = ''
# fix unit
if self.mode == 'frequency':
myunit = 'Hz'
elif self.mode == 'velocity':
myunit = 'km/s'
else: # channel
myunit = ''
for name in ['start', 'width']:
param = getattr(self, name)
new_param = self.__format_quantum_unit(param, myunit)
if new_param == None:
raise ValueError("Invalid unit for %s in mode %s: %s" % \
(name, self.mode, param))
setattr(self, name, new_param)
casalog.post("mode='%s': start=%s, width=%s, nchan=%d" % \
(self.mode, self.start, self.width, self.nchan))
# check length of selection parameters
if is_string_type(self.infiles):
nfile = 1
self.infiles = [ self.infiles ]
else:
nfile = len(self.infiles)
for name in ['field', 'spw', 'antenna', 'scanno']:
param = getattr(self, name)
if not self.__check_selection_length(param, nfile):
raise ValueError("Length of %s != infiles." % (name))
# set convsupport default
if self.convsupport >= 0 and self.gridfunction.upper() != 'SF':
casalog.post("user defined convsupport is ignored for %s kernel" % self.gridfunction, priority='WARN')
self.convsupport=-1
def __format_spw_string(self, spw):
"""
Returns formatted spw selection string which is accepted by imager.
"""
if type(spw) != str:
raise ValueError("The parameter should be string.")
if spw.strip() == '*': spw = ''
# WORKAROUND for CAS-6422, i.e., ":X~Y" fails while "*:X~Y" works.
if spw.startswith(":"): spw = '*' + spw
return spw
def __format_quantum_unit(self, data, unit):
"""
Returns False if data has an unit which in not a variation of
input unit.
Otherwise, returns input data as a quantum string. The input
unit is added to the return value if no unit is in data.
"""
my_qa = quanta()
if data == '' or my_qa.compare(data, unit):
return data
if my_qa.getunit(data) == '':
casalog.post("No unit specified. Using '%s'" % unit)
return '%f%s' % (data, unit)
return None
def __check_selection_length(self, data, nfile):
"""
Returns true if data is either a string, an array with length
1 or nfile
"""
if not is_string_type(data) and len(data) not in [1, nfile]:
return False
return True
def get_selection_param_for_ms(self, fileid, param):
"""
Returns valid selection string for a certain ms
Arguments
fileid : file idx in infiles list
param : string (array) selection value
"""
if is_string_type(param):
return param
elif len(param) == 1:
return param[0]
else:
return param[fileid]
def get_selection_idx_for_ms(self, file_idx):
"""
Returns a dictionary of selection indices for i-th MS in infiles
Argument: file idx in infiles list
"""
if file_idx < len(self.infiles) and file_idx > -1:
vis = self.infiles[file_idx]
field = self.get_selection_param_for_ms(file_idx, self.field)
spw = self.get_selection_param_for_ms(file_idx, self.spw)
spw = self.__format_spw_string(spw)
antenna = self.get_selection_param_for_ms(file_idx, self.antenna)
if antenna == -1: antenna = ''
scan = self.get_selection_param_for_ms(file_idx, self.scanno)
intent = self.get_selection_param_for_ms(file_idx, self.intent)
my_ms = ms()
sel_ids = my_ms.msseltoindex(vis=vis, spw=spw, field=field,
baseline=antenna, scan=scan)
fieldid = list(sel_ids['field']) if len(sel_ids['field']) > 0 else -1
baseline = self.format_ac_baseline(sel_ids['antenna1'])
scanid = list(sel_ids['scan']) if len(sel_ids['scan']) > 0 else ""
# SPW (need to get a list of valid spws instead of -1)
if len(sel_ids['channel']) > 0:
spwid = [ chanarr[0] for chanarr in sel_ids['channel'] ]
elif spw=="": # No spw selection
my_ms.open(vis)
try: spwinfo = my_ms.getspectralwindowinfo()
except: raise
finally: my_ms.close()
spwid = [int(idx) for idx in spwinfo.keys()]
else:
raise RuntimeError("Invalid spw selction, %s ,for MS %d" % (str(spw), file_idx))
return {'field': fieldid, 'spw': spwid, 'baseline': baseline, 'scan': scanid, 'intent': intent, 'antenna1': sel_ids['antenna1']}
else:
raise ValueError("Invalid file index, %d" % file_idx)
def format_ac_baseline(self, in_antenna):
""" format auto-correlation baseline string from antenna idx list """
# exact match string
if is_string_type(in_antenna):
if (len(in_antenna) != 0) and (in_antenna.find('&') == -1) \
and (in_antenna.find(';')==-1):
in_antenna =+ '&&&'
return in_antenna
# single integer -> list of int
if type(in_antenna)==int:
if in_antenna >=0:
in_antenna = [in_antenna]
else:
return -1
# format auto-corr string from antenna idices.
baseline = ''
for idx in in_antenna:
if len(baseline) > 0: baseline += ';'
if idx >= 0:
baseline += (str(idx) + '&&&')
return baseline
def compile(self):
# imaging mode
self.imager_param['mode'] = self.mode
# Work on selection of the first table in sorted list
# to get default restfreq and outframe
# chronological sort
sorted_vislist, sorted_timelist = mslisthelper.sort_mslist(self.infiles)
self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist]
mslisthelper.report_sort_result(sorted_vislist, sorted_timelist, self.sorted_idx, mycasalog=casalog)
# conform MS
conform_mslist(sorted_vislist, ignore_columns=[])
selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0])
self.__update_subtable_name(self.infiles[self.sorted_idx[0]])
# field
fieldid = selection_ids['field'][0] if type(selection_ids['field']) != int else selection_ids['field']
sourceid=-1
self.open_table(self.field_table)
source_ids = self.table.getcol('SOURCE_ID')
self.close_table()
if self.field == '' or fieldid ==-1:
sourceid = source_ids[0]
elif fieldid >= 0 and fieldid < len(source_ids):
sourceid = source_ids[fieldid]
else:
raise ValueError("No valid field in the first MS.")
# restfreq
if self.restfreq=='' and self.source_table != '':
self.open_table(self.source_table)
source_ids = self.table.getcol('SOURCE_ID')
for i in range(self.table.nrows()):
if sourceid == source_ids[i] \
and self.table.iscelldefined('REST_FREQUENCY',i) \
and (selection_ids['spw'] == -1 or \
self.table.getcell('SPECTRAL_WINDOW_ID', i) in selection_ids['spw']):
rf = self.table.getcell('REST_FREQUENCY',i)
if len(rf) > 0:
self.restfreq=self.table.getcell('REST_FREQUENCY',i)[0]
break
self.close_table()
casalog.post("restfreq set to %s" % self.restfreq, "INFO")
# REST_FREQUENCY column is optional (need retry if not exists)
self.imager_param['restfreq'] = self.restfreq
#
# spw (define representative spw id = spwid_ref)
spwid_ref = selection_ids['spw'][0] if type(selection_ids['spw']) != int else selection_ids['spw']
# Invalid spw selection should have handled at msselectiontoindex().
# -1 means all spw are selected.
self.open_table(self.spw_table)
if spwid_ref < 0:
for id in range(self.table.nrows()):
if self.table.getcell('NUM_CHAN',id) > 0:
spwid_ref = id
break
if spwid_ref < 0:
self.close_table()
msg='No valid spw id exists in the first table'
raise ValueError(msg)
self.allchannels = self.table.getcell('NUM_CHAN',spwid_ref)
freq_chan0 = self.table.getcell('CHAN_FREQ',spwid_ref)[0]
freq_inc0 = self.table.getcell('CHAN_WIDTH',spwid_ref)[0]
# in case rest frequency is not defined yet.
if self.restfreq=='':
self.restfreq = '%fHz' % self.table.getcell('CHAN_FREQ',spwid_ref).mean()
self.imager_param['restfreq'] = self.restfreq
casalog.post("Using mean freq of spw %d as restfreq: %s" %
(spwid_ref, self.restfreq), "INFO")
self.close_table()
self.imager_param['spw'] = -1 #spwid_ref
# outframe (force using the current frame)
self.imager_param['outframe'] = self.outframe
if self.outframe == '':
if len(self.infiles) > 1:
# The default will be 'LSRK'
casalog.post("Multiple MS inputs. The default outframe is set to 'LSRK'")
self.imager_param['outframe'] = 'LSRK'
else:
# get from MS
my_ms = ms()
my_ms.open(self.infiles[0])
spwinfo = my_ms.getspectralwindowinfo()
my_ms.close()
del my_ms
for key, spwval in spwinfo.items():
if spwval['SpectralWindowId'] == spwid_ref:
self.imager_param['outframe'] = spwval['Frame']
casalog.post("Using frequency frame of MS, '%s'" % self.imager_param['outframe'])
break
if self.imager_param['outframe'] == '':
raise Exception("Internal error of getting frequency frame of spw=%d." % spwid_ref)
else:
casalog.post("Using frequency frame defined by user, '%s'" % self.imager_param['outframe'])
# brightness unit
if len(self.brightnessunit) > 0:
if self.brightnessunit.lower() == 'k':
self.image_unit = 'K'
elif self.brightnessunit.lower() == 'jy/beam':
self.image_unit = 'Jy/beam'
else:
raise ValueError("Invalid brightness unit, %s" % self.brightnessunit)
# # antenna
# in_antenna = self.antenna # backup for future use
# if type(self.antenna)==int:
# if self.antenna >= 0:
# self.antenna=str(self.antenna)+'&&&'
# else:
# if (len(self.antenna) != 0) and (self.antenna.find('&') == -1) \
# and (self.antenna.find(';')==-1):
# self.antenna = self.antenna + '&&&'
def _configure_map_property(self):
selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0])
# stokes
if self.stokes == '':
self.stokes = 'I'
self.imager_param['stokes'] = self.stokes
# gridfunction
# outfile
if os.path.exists(self.outfile) and self.overwrite:
smart_remove(self.outfile)
if os.path.exists(self.outfile + '.weight') and self.overwrite:
smart_remove(self.outfile + '.weight')
# cell
cell = self.cell
if cell == '' or cell[0] == '':
# Calc PB
grid_factor = 3.
casalog.post("The cell size will be calculated using PB size of antennas in the first MS")
qPB = self._calc_PB(selection_ids['antenna1'])
cell = '%f%s' % (qPB['value']/grid_factor, qPB['unit'])
casalog.post("Using cell size = PB/%4.2F = %s" % (grid_factor, cell))
(cellx,celly) = sdutil.get_cellx_celly(cell, unit='arcmin')
self.imager_param['cellx'] = cellx
self.imager_param['celly'] = celly
# Calculate Pointing center and extent (if necessary)
# return a dictionary with keys 'center', 'width', 'height'
#imsize = self.imsize
imsize = sdutil._to_list(self.imsize, int) or \
sdutil._to_list(self.imsize, numpy.integer)
if imsize is None:
imsize = self.imsize if hasattr(self.imsize, '__iter__') else [ self.imsize ]
imsize = [ int(numpy.ceil(v)) for v in imsize ]
casalog.post("imsize is not integers. force converting to integer pixel numbers.", priority="WARN")
casalog.post("rounded-up imsize: %s --> %s" % (str(self.imsize), str(imsize)))
phasecenter = self.phasecenter
if self.phasecenter == "" or \
len(imsize) == 0 or imsize[0] < 1:
map_param = self._get_pointing_extent()
# imsize
if len(imsize) == 0 or imsize[0] < 1:
imsize = self._get_imsize(map_param['width'], map_param['height'], cellx, celly)
if self.phasecenter != "":
casalog.post("You defined phasecenter but not imsize. The image will cover as wide area as pointing in MS extends, but be centered at phasecenter. This could result in a strange image if your phasecenter is a part from the center of pointings", priority='WARN')
if imsize[0] > 1024 or imsize[1] > 1024:
casalog.post("The calculated image pixel number is larger than 1024. It could take time to generate the image depending on your computer resource. Please wait...", priority='WARN')
# phasecenter
# if empty, it should be determined here...
if self.phasecenter == "":
phasecenter = map_param['center']
# imsize
(nx,ny) = sdutil.get_nx_ny(imsize)
self.imager_param['nx'] = nx
self.imager_param['ny'] = ny
# phasecenter
self.imager_param['phasecenter'] = phasecenter
self.imager_param['movingsource'] = self.ephemsrcname
# channel map
sorted_vislist, _ = mslisthelper.sort_mslist(self.infiles)
self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist]
imhelper = cleanhelper_minimal(self.imager, self.infiles, sort_index=self.sorted_idx, casalog=casalog)
spwsel = str(',').join([str(spwid) for spwid in selection_ids['spw']])
srestf = self.imager_param['restfreq'] if is_string_type(self.imager_param['restfreq']) else "%fHz" % self.imager_param['restfreq']
(imnchan, imstart, imwidth) = imhelper.setChannelizeDefault(
self.mode,
spwsel,
self.field,
self.nchan,
self.start,
self.width,
self.imager_param['outframe'],
self.veltype,
self.imager_param['phasecenter'],
srestf
)
del imhelper
# start and width
if self.mode == 'velocity':
startval = [self.imager_param['outframe'], imstart]
widthval = imwidth
elif self.mode == 'frequency':
startval = [self.imager_param['outframe'], imstart]
widthval = imwidth
else: #self.mode==channel
startval = int(self.start)
widthval = int(self.width)
if self.nchan < 0: self.nchan = self.allchannels
self.imager_param['start'] = startval
self.imager_param['step'] = widthval
self.imager_param['nchan'] = imnchan #self.nchan
self.imager_param['projection'] = self.projection
def execute(self):
# imaging
casalog.post("Start imaging...", "INFO")
if len(self.infiles) == 1:
self.open_imager(self.infiles[0])
selection_ids = self.get_selection_idx_for_ms(0)
spwsel = self.get_selection_param_for_ms(0, self.spw)
if spwsel.strip() in ['', '*']: spwsel = selection_ids['spw']
### TODO: channel selection based on spw
ok = self.imager.selectvis(field=selection_ids['field'],
#spw=selection_ids['spw'],
spw=spwsel,
nchan=-1, start=0, step=1,
baseline=selection_ids['baseline'],
scan=selection_ids['scan'],
intent=selection_ids['intent'])
if not ok:
raise ValueError("Selection is empty: you may want to review this MS selection")
else:
self.close_imager()
#self.sorted_idx.reverse()
for idx in self.sorted_idx.__reversed__():
name = self.infiles[idx]
selection_ids = self.get_selection_idx_for_ms(idx)
spwsel = self.get_selection_param_for_ms(idx, self.spw)
if spwsel.strip() in ['', '*']: spwsel = selection_ids['spw']
### TODO: channel selection based on spw
self.imager.selectvis(vis=name, field=selection_ids['field'],
#spw=selection_ids['spw'],
spw = spwsel,
nchan=-1, start=0, step=1,
baseline=selection_ids['baseline'],
scan=selection_ids['scan'],
intent=selection_ids['intent'])
# need to do this
self.is_imager_opened = True
# it should be called after infiles are registered to imager
self._configure_map_property()
casalog.post("Using phasecenter \"%s\""%(self.imager_param['phasecenter']), "INFO")
self.imager.defineimage(**self.imager_param)#self.__get_param())
self.imager.setoptions(ftmachine='sd', gridfunction=self.gridfunction)
self.imager.setsdoptions(pointingcolumntouse=self.pointingcolumn, convsupport=self.convsupport, truncate=self.truncate, gwidth=self.gwidth, jwidth=self.jwidth, minweight = 0., clipminmax=self.clipminmax)
# Create images
imgtype_suffix = {'singledish': '', 'coverage' : '.weight'}
for img_type, img_suffix in imgtype_suffix.items():
img_file = self.outfile + img_suffix
msg_fmt = string.Template("$state {img_type} image {img_file}".format(
img_type=img_type,
img_file=img_file))
casalog.post(msg_fmt.substitute(state="Generating"), "INFO")
self.imager.makeimage(type=img_type, image=img_file)
if not os.path.exists(img_file):
raise RuntimeError(msg_fmt.substitute(state="Failed to generate"))
casalog.post(msg_fmt.substitute(state="Generated"), "INFO")
# Close imager
self.close_imager()
# Convert output images to proper output frame and set brightness unit (if necessary)
my_ia = image()
my_ia.open(self.outfile)
csys = my_ia.coordsys()
csys.setconversiontype(spectral=csys.referencecode('spectra')[0])
my_ia.setcoordsys(csys.torecord())
if len(self.image_unit) > 0:
casalog.post("Setting brightness unit '%s' to image." % self.image_unit)
my_ia.setbrightnessunit(self.image_unit)
csys.done()
my_ia.close()
# CAS-12984 set brightness unit for weight image to ''
weightfile = self.outfile + imgtype_suffix['coverage']
my_ia.open(weightfile)
my_ia.setbrightnessunit('')
my_ia.close()
# Mask image pixels whose weight are smaller than minweight.
# Weight image should have 0 weight for pixels below < minweight
casalog.post("Start masking the map using minweight = %f" % \
self.minweight, "INFO")
my_ia.open(weightfile)
try:
stat=my_ia.statistics(mask="'"+weightfile+"' > 0.0", robust=True)
valid_pixels=stat['npts']
except RuntimeError as e:
if str(e).find('No valid data found.') >= 0:
valid_pixels = [0]
else:
raise
if len(valid_pixels) == 0 or valid_pixels[0] == 0:
my_ia.close()
casalog.post("All pixels have zero weight. This means the imaged region contains no MS data. Mask will not be set. Please check your image parameters.","WARN")
return
median_weight = stat['median'][0]
weight_threshold = median_weight * self.minweight
casalog.post("Median of weight in the map is %f" % median_weight, \
"INFO")
casalog.post("Pixels in map with weight <= median(weight)*minweight = %f will be masked." % \
(weight_threshold),"INFO")
###Leaving the original logic to calculate the number of masked pixels via
###product of median of and min_weight (which i don't understand the logic)
### if one wanted to find how many pixel were masked one could easily count the
### number of pixels set to false
### e.g after masking self.outfile below one could just do this
### nmasked_pixels=tb.calc('[select from "'+self.outfile+'"/mask0'+'" giving [nfalse(PagedArray )]]')
my_tb = table()
nmask_pixels=0
nchan=stat['trc'][3]+1
casalog.filter('ERROR') ### hide the useless message of tb.calc
### doing it by channel to make sure it does not go out of memory
####tab.calc try to load the whole chunk in ram
for k in range(nchan):
nmask_pixels += my_tb.calc('[select from "'+weightfile+'" giving [ntrue(map[,,,'+str(k)+'] <='+str(median_weight*self.minweight)+')]]')['0'][0]
casalog.filter() ####set logging back to normal
casalog.filter() ####set logging back to normal
imsize=numpy.product(my_ia.shape())
my_ia.close()
# Modify default mask
my_ia.open(self.outfile)
my_ia.calcmask("'%s'>%f" % (weightfile,weight_threshold), asdefault=True)
my_ia.close()
masked_fraction = 100.*(1. - (imsize - nmask_pixels) / float(valid_pixels[0]) )
casalog.post("This amounts to %5.1f %% of the area with nonzero weight." % \
( masked_fraction ),"INFO")
casalog.post("The weight image '%s' is returned by this task, if the user wishes to assess the results in detail." \
% (weightfile), "INFO")
# Calculate theoretical beam size
casalog.post("Calculating image beam size.")
if self.gridfunction.upper() not in ['SF']:
casalog.post("Beam size definition for '%s' kernel is experimental." % self.gridfunction, priority='WARN')
casalog.post("You may want to take careful look at the restoring beam in the image.",priority='WARN')
my_msmd = msmetadata()
# antenna diameter and blockage
ref_ms_idx = self.sorted_idx[0]
ref_ms_name = self.infiles[ref_ms_idx]
selection_ids = self.get_selection_idx_for_ms(ref_ms_idx)
ant_idx = selection_ids['antenna1']
diameter = self._get_average_antenna_diameter(ant_idx)
my_msmd.open(ref_ms_name)
ant_name = my_msmd.antennanames(ant_idx)
my_msmd.close()
is_alma = False
for name in ant_name:
if name[0:2] in ["PM", "DV", "DA", "CM"]:
is_alma = True
break
blockage = "0.75m" if is_alma else "0.0m" # unknown blockage diameter
# output reference code
my_ia.open(self.outfile)
csys = my_ia.coordsys()
my_ia.close()
outref = csys.referencecode('direction')[0]
cell = list(csys.increment(type='direction',format='s')['string'])
csys.done()
# pointing sampling
ref_ms_spw = self.get_selection_param_for_ms(ref_ms_idx,self.spw)
ref_ms_field = self.get_selection_param_for_ms(ref_ms_idx,self.field)
ref_ms_scan = self.get_selection_param_for_ms(ref_ms_idx,self.scanno)
# xSampling, ySampling, angle = sdutil.get_ms_sampling_arcsec(ref_ms_name, spw=ref_ms_spw,
# antenna=selection_ids['baseline'],
# field=ref_ms_field,
# scan=ref_ms_scan,#timerange='',
# outref=outref)
# obtain sampling interval for beam calculation.
self.open_imager(ref_ms_name)
ok = self.imager.selectvis(field=ref_ms_field,
#spw=selection_ids['spw'],
spw=ref_ms_spw,
nchan=-1, start=0, step=1,
baseline=selection_ids['baseline'],
scan=ref_ms_scan,
intent=selection_ids['intent'])
if len(ant_idx) > 1:
casalog.post("Using only antenna %s to calculate sampling interval" % ant_name[0])
ptg_samp = self.imager.pointingsampling(pattern='raster', ref=outref,
movingsource=self.ephemsrcname,
pointingcolumntouse=self.pointingcolumn,
antenna=('%s&&&' % ant_name[0]))
self.close_imager()
my_qa = quanta()
xSampling, ySampling = my_qa.getvalue(my_qa.convert(ptg_samp['sampling'], 'arcsec'))
angle = my_qa.getvalue(my_qa.convert(ptg_samp['angle'], "deg"))[0]
casalog.post("Detected raster sampling = [{x:f}, {y:f}] arcsec".format(
x=xSampling,
y=ySampling
))
# handling of failed sampling detection
valid_sampling = True
sampling = [xSampling, ySampling]
if abs(xSampling) < 2.2e-3 or not numpy.isfinite(xSampling):
casalog.post("Invalid sampling=%s arcsec. Using the value of orthogonal direction=%s arcsec" % (xSampling, ySampling), priority="WARN")
sampling = [ ySampling ]
angle = 0.0
valid_sampling = False
if abs(ySampling) < 1.0e-3 or not numpy.isfinite(ySampling):
if valid_sampling:
casalog.post("Invalid sampling=%s arcsec. Using the value of orthogonal direction=%s arcsec" % (ySampling, xSampling), priority="WARN")
sampling = [ xSampling ]
angle = 0.0
valid_sampling = True
# reduce sampling and cell if it's possible
if len(sampling)>1 and abs(sampling[0]-sampling[1]) <= 0.01*abs(sampling[0]):
sampling = [sampling[0]]
angle = 0.0
if cell[0]==cell[1]: cell = [cell[0]]
if valid_sampling:
# actual calculation of beam size
bu = sdbeamutil.TheoreticalBeam()
bu.set_antenna(diameter,blockage)
bu.set_sampling(sampling, "%fdeg" % angle)
bu.set_image_param(cell, self.restfreq, self.gridfunction,
self.convsupport, self.truncate, self.gwidth,
self.jwidth,is_alma)
bu.summary()
imbeam_dict = bu.get_beamsize_image()
casalog.post("Setting image beam: major={major}, minor={minor}, pa={pa}".format(
major=imbeam_dict['major'],
minor=imbeam_dict['minor'],
pa=imbeam_dict['pa']
))
# set beam size to image
my_ia.open(self.outfile)
my_ia.setrestoringbeam(**imbeam_dict)
my_ia.close()
else:
#BOTH sampling was invalid
casalog.post("Could not detect valid raster sampling. Exitting without setting beam size to image", priority='WARN')
def _calc_PB(self, antenna):
"""
Calculate the primary beam size of antenna, using dish diamenter
and rest frequency
Average antenna diamter and reference frequency are adopted for
calculation.
The input argument should be a list of antenna IDs.
"""
casalog.post("Calculating Pirimary beam size:")
# CAS-5410 Use private tools inside task scripts
my_qa = quanta()
pb_factor = 1.175
# Reference frequency
ref_freq = self.restfreq
if type(ref_freq) in [float, numpy.float64]:
ref_freq = my_qa.tos(my_qa.quantity(ref_freq, 'Hz'))
if not my_qa.compare(ref_freq, 'Hz'):
msg = "Could not get the reference frequency. " + \
"Your data does not seem to have valid one in selected field.\n" + \
"PB is not calculated.\n" + \
"Please set restreq or cell manually to generate an image."
raise Exception(msg)
# Antenna diameter
antdiam_ave = self._get_average_antenna_diameter(antenna)
# Calculate PB
wave_length = 0.2997924 / my_qa.convert(my_qa.quantity(ref_freq),'GHz')['value']
D_m = my_qa.convert(antdiam_ave, 'm')['value']
lambda_D = wave_length / D_m * 3600. * 180 / numpy.pi
PB = my_qa.quantity(pb_factor*lambda_D, 'arcsec')
# Summary
casalog.post("- Antenna diameter: %s m" % D_m)
casalog.post("- Reference Frequency: %s" % ref_freq)
casalog.post("PB size = %5.3f * lambda/D = %s" % (pb_factor, my_qa.tos(PB)))
return PB
def _get_imsize(self, width, height, dx, dy):
casalog.post("Calculating pixel size.")
# CAS-5410 Use private tools inside task scripts
my_qa = quanta()
ny = numpy.ceil( ( my_qa.convert(height, my_qa.getunit(dy))['value'] / \
my_qa.getvalue(dy) ) )
nx = numpy.ceil( ( my_qa.convert(width, my_qa.getunit(dx))['value'] / \
my_qa.getvalue(dx) ) )
casalog.post("- Map extent: [%s, %s]" % (my_qa.tos(width), my_qa.tos(height)))
casalog.post("- Cell size: [%s, %s]" % (my_qa.tos(dx), my_qa.tos(dy)))
casalog.post("Image pixel numbers to cover the extent: [%d, %d] (projected)" % \
(nx+1, ny+1))
return (int(nx+1), int(ny+1))
def _get_pointing_extent(self):
### MS selection is ignored. This is not quite right.
casalog.post("Calculating map extent from pointings.")
# CAS-5410 Use private tools inside task scripts
my_qa = quanta()
ret_dict = {}
colname = self.pointingcolumn.upper()
# MSs should be registered to imager
if not self.is_imager_opened:
raise RuntimeError('Internal error: imager should be opened here.')
if self.phasecenter == "":
# defaut is J2000
base_mref = 'J2000'
elif isinstance(self.phasecenter, int) or self.phasecenter.isdigit():
# may be field id
self.open_table(self.field_table)
base_mref = self.table.getcolkeyword('PHASE_DIR', 'MEASINFO')['Ref']
self.close_table()
else:
# may be phasecenter is explicitly specified
pattern = '^([\-\+]?[0-9.]+([eE]?-?[0-9])?)|([\-\+]?[0-9][:h][0-9][:m][0-9.]s?)|([\-\+]?[0-9][.d][0-9][.d][0-9.]s?)$'
items = self.phasecenter.split()
base_mref = 'J2000'
for i in items:
s = i.strip()
if re.match(pattern, s) is None:
base_mref = s
break
mapextent = self.imager.mapextent(ref=base_mref, movingsource=self.ephemsrcname,
pointingcolumntouse=colname)
if mapextent['status'] is True:
qheight = my_qa.quantity(mapextent['extent'][1], 'rad')
qwidth = my_qa.quantity(mapextent['extent'][0], 'rad')
qcent0 = my_qa.quantity(mapextent['center'][0], 'rad')
qcent1 = my_qa.quantity(mapextent['center'][1], 'rad')
scenter = '{ref} {longitude} {latitude}'.format(
ref=base_mref,
longitude=my_qa.formxxx(qcent0, 'hms'),
latitude=my_qa.formxxx(qcent1, 'dms')
)
casalog.post("- Pointing center: {center}".format(center=scenter))
casalog.post("- Pointing extent: [{width}, {height}] (projected)".format(
width=my_qa.tos(qwidth),
height=my_qa.tos(qheight)
))
ret_dict['center'] = scenter
ret_dict['width'] = qwidth
ret_dict['height'] = qheight
else:
casalog.post('Failed to derive map extent from the MSs registered to the imager probably due to mising valid data.', priority='SEVERE')
ret_dict['center'] = ''
ret_dict['width'] = my_qa.quantity(0.0, 'rad')
ret_dict['height'] = my_qa.quantity(0.0, 'rad')
return ret_dict
def _get_x_minmax(self, x):
# assumed the x is in unit of rad.
pi2 = 2. * numpy.pi
x = (x % pi2)
npart = 4
dlon = pi2/float(npart)
pos = [int(v/dlon) for v in x]
voids = [False for dummy in range(npart)]
for ipos in range(npart):
try: dummy = pos.index(ipos)
except: voids[ipos] = True
if not any(voids):
raise Exception("Failed to find global pointing gap. The algorithm requires at least 2PI/%d of pointing gap" % npart)
rot_pos = []
if (not voids[0]) and (not voids[npart-1]):
gmax = -1
for idx in range(npart-2, 0, -1):
if voids[idx]:
gmax = idx
break
if gmax < 0:
raise Exception("Failed to detect gap max")
rot_pos = range(gmax+1, npart)
for idx in range(len(x)):
x[idx] = (x[idx] - pi2) if pos[idx] in rot_pos else x[idx]
return (x.min(), x.max())
def __update_subtable_name(self, msname):
self.open_table(msname)
keys = self.table.getkeywords()
self.close_table()
self.field_table = sdutil.get_subtable_name(keys['FIELD'])
self.spw_table = sdutil.get_subtable_name(keys['SPECTRAL_WINDOW'])
self.source_table = sdutil.get_subtable_name(keys['SOURCE'])
self.antenna_table = sdutil.get_subtable_name(keys['ANTENNA'])
self.polarization_table = sdutil.get_subtable_name(keys['POLARIZATION'])
self.observation_table = sdutil.get_subtable_name(keys['OBSERVATION'])
self.pointing_table = sdutil.get_subtable_name(keys['POINTING'])
self.data_desc_table = sdutil.get_subtable_name(keys['DATA_DESCRIPTION'])
def _get_average_antenna_diameter(self, antenna):
my_qa = quanta()
self.open_table(self.antenna_table)
try:
antdiam_unit = self.table.getcolkeyword('DISH_DIAMETER', 'QuantumUnits')[0]
diams = self.table.getcol('DISH_DIAMETER')
finally:
self.close_table()
if len(antenna) == 0:
antdiam_ave = my_qa.quantity(diams.mean(), antdiam_unit)
else:
d_ave = sum([diams[idx] for idx in antenna])/float(len(antenna))
antdiam_ave = my_qa.quantity(d_ave, antdiam_unit)
return antdiam_ave