from __future__ import absolute_import
from glob import glob
import os
import re
import sys
import shutil
# the solution here to findCalModels should be reworked
# CASA5 uses the casa dict to set the default roots parameter values
# and that dict is not available in CASA6, where ctsys.resolve is
# available
# Here: defaultRoots is defined differently for CASA6 and CASA5 and
# that values is used as the default value for the roots parameter in
# findCalModels.
# Inside findCalModels there is a new section at the top to use
# ctsys.resolve to make the initial attempt to find the models.
# If that fails or is unavailable it falls back to the CASA5 code.
# get is_python3 and is_CASA6
from casatasks.private.casa_transition import *
if is_CASA6:
from .setjy_helper import *
from .parallel.parallel_data_helper import ParallelDataHelper
from .parallel.parallel_task_helper import ParallelTaskHelper
from .mstools import write_history
from casatools import ctsys, ms, imager, calibrater
from casatasks import casalog
# used in one type comparison
strType = str
# default roots argument for findCalModels
defaultRoots = ['.']
else:
from setjy_helper import *
from taskinit import *
from taskinit import mstool as ms
from taskinit import imtool as imager
from taskinit import cbtool as calibrater
from mstools import write_history
from parallel.parallel_data_helper import ParallelDataHelper
from parallel.parallel_task_helper import ParallelTaskHelper
# used in one type comparison
strType = string
# default roots argument for findCalModels
defaultRoots = ['.', casa['dirs']['data']]
# Helper class for Multi-MS processing (by SC)
class SetjyHelper():
def __init__(self, msfile=None):
self.__msfile = msfile
def resetModelCol(self):
rstatus = True
# Hide the log info
casalog.post("Resetting the log filter to WARN", "DEBUG")
casalog.post("Initialize the MODEL columns of sub-MSs to default 1", "DEBUG")
casalog.filter('WARN')
myms = ms()
try:
try:
myms.open(self.__msfile)
submslist = myms.getreferencedtables()
parentms = myms.name()
myms.close()
except:
if (os.path.exists(os.path.join(self.__msfile,'SUBMSS'))):
casalog.post("MS may be corrupted. Try to initialize sub-MSs anyway...","DEBUG")
submslist = [os.path.abspath(os.path.join(self.__msfile,'SUBMSS',fn))
for fn in os.listdir(os.path.join(self.__msfile,'SUBMSS')) if fn.endswith('.ms') ]
else:
rstatus = False
mycb = calibrater()
if parentms!= submslist[0]:
for subms in submslist:
mycb.open(subms, addcorr=False, addmodel=True)
mycb.close()
except:
rstatus = False
casalog.filter('INFO')
return rstatus
def setjy(vis=None, field=None, spw=None,
selectdata=None, timerange=None, scan=None, intent=None, observation=None,
scalebychan=None, standard=None, model=None,
listmodels=None, fluxdensity=None, spix=None, reffreq=None, polindex=None,
polangle=None, rotmeas=None, fluxdict=None,
useephemdir=None, interpolation=None, usescratch=None, ismms=None):
"""Fills the model column for flux density calibrators."""
casalog.origin('setjy')
casalog.post("standard="+standard,'DEBUG1')
mylocals = locals()
if not listmodels: # listmmodels=T does not require vis
sh = SetjyHelper(vis)
rstat = sh.resetModelCol()
# Take care of the trivial parallelization
if ( not listmodels and ParallelDataHelper.isMMSAndNotServer(vis) and usescratch):
# jagonzal: We actually operate in parallel when usescratch=True because only
# in this case there is a good trade-off between the parallelization overhead
# and speed up due to the load involved with MODEL_DATA column creation
# Create the default MODEL columns in all sub-MSs to avoid corruption of the MMS
# when there are NULL MS selections
#
# TT: Use ismms is used to change behavior of some of the execption handling
# for MMS case. It is a hidden task parameter only modified when input vis
# is identified as MMS via SetjyHelper.resetModel().
#sh = SetjyHelper(vis)
#rstat = sh.resetModelCol()
if not rstat:
raise RuntimeError("Could not initialize MODEL columns in sub-MSs", 'SEVERE')
ismms=rstat
mylocals['ismms']=ismms
helper = ParallelTaskHelper('setjy', mylocals)
helper._consolidateOutput = False
retval = helper.go()
# Remove the subMS names from the returned dictionary. Pick the first dictionary from
# the list of subMS dictionaries
if (any(isinstance(v, dict) for v in retval.values())):
for subMS in retval:
dict_i = retval[subMS]
# Take the first *non-empty dict* (CAS-13645)
if isinstance(dict_i,dict) and dict_i:
retval = dict_i
break
else:
raise RuntimeError("Error in parallel processing of MMS, retval: {}".
format(retval),'SEVERE')
else:
retval = setjy_core(vis, field, spw, selectdata, timerange,
scan, intent, observation, scalebychan, standard, model,
listmodels, fluxdensity, spix, reffreq,
polindex, polangle, rotmeas, fluxdict,
useephemdir, interpolation, usescratch, ismms)
return retval
def setjy_core(vis=None, field=None, spw=None,
selectdata=None, timerange=None, scan=None, intent=None, observation=None,
scalebychan=None, standard=None, model=None, listmodels=None,
fluxdensity=None, spix=None, reffreq=None,
polindex=None, polangle=None, rotmeas=None, fluxdict=None,
useephemdir=None, interpolation=None, usescratch=None, ismms=None):
"""Fills the model column for flux density calibrators."""
retval = {}
clnamelist=[]
# remove componentlist generated
deletecomp = True
#deletecomp = False
try:
# Here we only list the models available, but don't perform any operation
if listmodels:
casalog.post("Listing model candidates (listmodels == True).")
if vis:
casalog.post("%s is NOT being modified." % vis)
if standard=='Butler-JPL-Horizons 2012':
tpm_supported_objects = ['Ceres','Lutetia','Pallas','Vesta']
ssmoddirs = findCalModels(target='SolarSystemModels',
exts=['.im','.ms','tab'])
sstpmdirs = findCalModels(target='SolarSystemModels',
exts=['.im','.ms','tab'])
if ssmoddirs==set([]):
casalog.post("No models were found. Missing SolarSystemModels in the CASA data directory","WARN")
for d in ssmoddirs:
lsmodims(d,modpat='*Tb*.dat', header='Tb models of solar system objects available for %s' % standard)
for d in sstpmdirs:
lsmodims(d,modpat='*fd_time.dat', header='Time variable models of asteroids available for %s [only applicable for the observation date 2014.01.01 0UT and beyond]' % standard)
elif standard=='Butler-JPL-Horizons 2010':
availmodellist=['Venus', 'Mars', 'Jupiter', 'Uranus', 'Neptune', 'Pluto',
'Io', 'Europa', 'Ganymede', 'Callisto', 'Titan','Triton',
'Ceres', 'Pallas', 'Vesta', 'Juno', 'Victoria', 'Davida']
casalog.post("Solar system objects recognized by %s:" % standard)
casalog.posta(vailmodellist)
else:
lsmodims('.', modpat='*.im* *.mod*')
calmoddirs = findCalModels()
ssmoddirs=None
for d in calmoddirs:
lsmodims(d)
# Actual operation, when either the MODEL_DATA column or visibility model header are set
else:
if not os.path.isdir(vis):
#casalog.post(vis + " must be a valid MS unless listmodels is True.",
# "SEVERE")
raise Exception("%s is not a valid MS" % vis)
myms = ms()
myim = imager()
if ismms==None: ismms=False
if type(vis) == str and os.path.isdir(vis):
n_selected_rows = nselrows(vis, field, spw, observation, timerange, scan, intent, usescratch, ismms)
# jagonzal: When usescratch=True, creating the MODEL column only on a sub-set of
# Sub-MSs causes problems because ms::open requires all the tables in ConCatTable
# to have the same description (MODEL data column must exist in all Sub-MSs)
#
# This is a bit of an over-doing but it is necessary for the sake of transparency
#
# Besides, this is also the same behavior as when running setjy vs a normal MS
#
# Finally, This does not affect the normal MS case because nselrows throws an
# exception when the user enters an invalid data selection, but it works for the
# MMS case because every sub-Ms contains a copy of the entire MMS sub-tables
##if ((not n_selected_rows) and ((not usescratch) or (standard=="Butler-JPL-Horizons 2012"))) :
#mylocals = locals()
#if ((not n_selected_rows) and ((usescratch) or (standard=="Butler-JPL-Horizons 2012"))) :
if ((not n_selected_rows) and ((ismms) or (standard=="Butler-JPL-Horizons 2012"))) :
# jagonzal: Turn this SEVERE into WARNING, as explained above
casalog.post("No rows were selected.", "WARNING")
return {}
else:
if (not n_selected_rows):
raise Exception("No rows were selected. Please check your data selection")
myim.open(vis, usescratch=usescratch)
else:
raise Exception('Visibility data set not found - please verify the name')
# remove modimage handling to use 'model' parameter only (use modimage internally??)
#if modimage==None: # defined as 'hidden' with default '' in the xml
# # but the default value does not seem to set so deal
# # with it here...
# modimage=''
#if model:
# modimage=model
#elif not model and modimage:
# casalog.post("The modimage parameter is deprecated please use model instead", "WARNING")
# If modimage is not an absolute path, see if we can find exactly 1 match in the likely places.
#if modimage and modimage[0] != '/':
if model and model[0] != '/':
cwd = os.path.abspath('.')
calmoddirs = [cwd]
# casa dict unavailable in CASA6
if is_CASA6:
calmoddirs += findCalModels()
else:
calmoddirs += findCalModels(roots=[cwd,
casa['dirs']['data']])
candidates = []
for calmoddir in calmoddirs:
cand = os.path.join(calmoddir,model)
if os.path.isdir(cand):
candidates.append(cand)
if not candidates:
raise RuntimeError("%s was not found for model in %s." %(model,
', '.join(calmoddirs)))
elif len(candidates) > 1:
casalog.post("More than 1 candidate for model was found:",
'SEVERE')
for c in candidates:
casalog.post("\t" + c, 'SEVERE')
raise RuntimeError("Please pick 1 and use the absolute path (starting with /).")
else:
model = candidates[0]
casalog.post("Using %s for model." % model, 'INFO')
# Write the parameters to HISTORY before the tool writes anything.
try:
param_names = setjy.__code__.co_varnames[:setjy.__code__.co_argcount]
if is_python3:
vars = locals()
param_vals = [vars[p] for p in param_names]
else:
param_vals = [eval(p) for p in param_names]
write_history(myms, vis, 'setjy', param_names,
param_vals, casalog)
except Exception as instance:
casalog.post("*** Error \'%s\' updating HISTORY" % (instance),
'WARN')
# Split the process for solar system objects
# To maintain the behavior of SetJy such that if a flux density is specified
# we use it rather than the model, we need to check the state of the
# input fluxdensity JSK 10/25/2012
# TT comment 07/01/2013: some what redundant as fluxdensity is moved to
# a subparameter but leave fluxdensity=-1 for backward compatibility
userFluxDensity = fluxdensity is not None
if userFluxDensity and isinstance(fluxdensity, list):
userFluxDensity = fluxdensity[0] > 0.0
elif userFluxDensity:
userFluxDensity = fluxdensity > 0.0
#in mpirun somehow subparameter defaulting does not work properly,
#it seems to pick up default defined in "constraint" clause in the task xml instead
#of the one defined in "param", so userFluxDensity is not reliable to use in here
# (and somewhat redundant in this case).
#if standard=="Butler-JPL-Horizons 2012" and not userFluxDensity:
if standard=="Butler-JPL-Horizons 2012":
casalog.post("Using Butler-JPL-Horizons 2012")
ssmoddirs = findCalModels(target='SolarSystemModels',
exts=['.im','.ms','tab'])
if ssmoddirs==set([]):
raise Exception("Missing Tb or fd models in the data directory")
setjyutil=ss_setjy_helper(myim,vis,casalog)
retval=setjyutil.setSolarObjectJy(field=field,spw=spw,scalebychan=scalebychan,
timerange=timerange,observation=str(observation), scan=scan,
intent=intent, useephemdir=useephemdir,usescratch=usescratch)
clnamelist=setjyutil.getclnamelist()
else:
# Need to branch out the process for fluxscale since the input dictionary may
# contains multiple fields. Since fluxdensity parameter is just a vector contains
# IQUV flux densities and so need to run im.setjy per field
if standard=="fluxscale":
instandard="Perley-Butler 2010"
# function to return selected field, spw, etc
fieldidused=parse_fluxdict(fluxdict, vis, field, spw, observation, timerange,
scan, intent, usescratch)
if len(fieldidused):
retval={}
for selfld in fieldidused:
#selspix=fluxdict[selfld]["spidx"][1] # setjy only support alpha for now
selspix=fluxdict[selfld]["spidx"][1:] # omit c0 (=log(So))
# set all (even if fluxdensity = -1
if spw=='':
selspw = []
invalidspw = []
for ispw in fluxdict["spwID"].tolist():
# skip spw if fluxd=-1
if fluxdict[selfld][str(ispw)]["fluxd"][0]>-1:
selspw.append(ispw)
else:
invalidspw.append(ispw)
if len(invalidspw):
casalog.post("Fluxdict does not contains valid fluxdensity for spw="+
str(invalidspw)+
". Model will not set for those spws.", "WARN")
else: # use spw selection
selspw = ms.msseltoindex(vis,spw=spw)["spw"].tolist()
if fluxdict[selfld]["fitFluxd"]:
selfluxd = fluxdict[selfld]["fitFluxd"]
selreffreq = fluxdict[selfld]["fitRefFreq"]
else:
# use first selected spw's flux density
selfluxd = fluxdict[selfld][str(selspw[0])]['fluxd']
# assuming the fluxscale reporting the center freq of a given spw
selreffreq=fluxdict["freq"][selspw[0]]
casalog.post("Use fluxdensity=%s, reffreq=%s, spix=%s" %
(selfluxd,selreffreq,selspix))
curretval=myim.setjy(field=selfld,spw=selspw,modimage=model,
# enable spix in list
fluxdensity=selfluxd, spix=selspix, reffreq=selreffreq,
#fluxdensity=selfluxd, spix=[selspix], reffreq=selreffreq,
standard=instandard, scalebychan=scalebychan,
polindex=polindex, polangle=polangle, rotmeas=rotmeas,
time=timerange, observation=str(observation), scan=scan,
intent=intent, interpolation=interpolation)
retval.update(curretval)
else:
raise Exception("No field is selected. Check fluxdict and field selection.")
else:
influxdensity=fluxdensity
if standard=="manual":
instandard="Perley-Butler 2010" # default standard when fluxdensity=-1 or 1
else:
instandard=standard
# From the task, fluxdensity=-1 always for standard!="manual" (but possible to
# for im.setjy to run with both the stanadard and fluxdensity specified.
# Until CAS-6463 is fixed, need to catch and override inconsistent parameters.
if userFluxDensity and instandard!='manual':
influxdensity=-1
#raise Exception("Use standard=\"manual\" to set the specified fluxdensity.")
casalog.post("*** fluxdensity > 0 but standard != 'manual' (possibly interaction with globals), override to set fluxdensity=-1.", 'WARN')
if spix==[]: # handle the default
spix=0.0
# need to modify imager to accept double array for spix
retval=myim.setjy(field=field, spw=spw, modimage=model, fluxdensity=influxdensity,
spix=spix, reffreq=reffreq, standard=instandard, scalebychan=scalebychan,
polindex=polindex, polangle=polangle, rotmeas=rotmeas,
time=timerange, observation=str(observation), scan=scan, intent=intent,
interpolation=interpolation)
myim.close()
finally:
if standard=='Butler-JPL-Horizons 2012':
for cln in clnamelist:
if deletecomp and os.path.exists(cln) and os.path.isdir(cln):
shutil.rmtree(cln,True)
if os.path.exists(os.path.basename(vis)+"_tmp_setjyCLfile"):
shutil.rmtree(os.path.basename(vis)+"_tmp_setjyCLfile")
return retval
def better_glob(pats):
"""
Unlike ls, glob.glob('pat1 pat2') does not return
the union of matches to pat1 and pat2. This does.
"""
retset = set([])
patlist = pats.split()
for p in patlist:
retset.update(glob(p))
return retset
def lsmodims(path, modpat='*', header='Candidate models'):
"""
Does an ls -d of files or directories in path matching modpat.
Header describes what is being listed.
"""
if os.path.isdir(path):
if better_glob(os.path.join(path,modpat)):
casalog.post("\n%s (%s) in %s:" % (header, modpat, path))
sys.stdout.flush()
os.system('cd ' + path + ';ls -d ' + modpat)
else:
casalog.post("\nNo %s matching '%s' found in %s" % (header.lower(),
modpat, path))
def findCalModels(target='CalModels',
roots=defaultRoots,
#permexcludes = ['.svn', 'regression', 'ephemerides',
permexcludes = ['regression', 'ephemerides',
'geodetic', 'gui'],
exts=['.ms', '.im', '.tab']):
"""
Returns a set of directories ending in target that are in the trees of roots.
Because casa['dirs']['data'] can contain a lot, and CASA tables are
directories, branches matching permexcludes or exts are excluded for speed.
"""
retset = set([])
##
## first attempt to resolve using data path - only available in CASA6
##
if is_CASA6:
standard_locations = { 'CalModels': [ 'nrao/VLA/CalModels' ],
'SolarSystemModels': [ 'alma/SolarSystemModels' ] }
if target in standard_locations:
for p in standard_locations[target]:
candidate = ctsys.resolve(p)
if os.path.isdir(candidate):
retset.add(candidate)
if len(retset) > 0:
return retset
# exclulde all hidden(.xxx) directories
regex = re.compile('^\.\S+')
# extra skip for OSX
if sys.platform == "darwin":
permexcludes.append('Frameworks')
permexcludes.append('Library')
permexcludes.append('lib')
for root in roots:
# Do a walk to find target directories in root.
# 7/5/2011: glob('/export/data_1/casa/gnuactive/data/*/CalModels/*') doesn't work.
for path, dirs, fnames in os.walk(root, followlinks=True):
realpath = os.path.realpath(path)
if not realpath in retset:
excludes = permexcludes[:]
for ext in exts:
excludes += [d for d in dirs if ext in d]
excludes += [m.group(0) for d in dirs for m in [regex.search(d)] if m]
for d in excludes:
if d in dirs:
dirs.remove(d)
if path.split('/')[-1] == target:
# make path realpath to resolve the issue with duplicated path resulted from symlinks
retset.add(realpath)
return retset
def nselrows(vis, field='', spw='', obs='', timerange='', scan='', intent='', usescratch=None, ismms=False):
# modified to use ms.msselect. If no row is selected ms.msselect will
# raise an exception - TT 2013.12.13
retval = 0
myms = ms()
#msselargs = {'vis': vis}
msselargs = {}
if field:
msselargs['field'] = field
if spw:
msselargs['spw'] = spw
if intent:
msselargs['scanintent'] = intent
# only applicable for usescratch=T
if usescratch:
if obs:
if not type(obs) == strType:
sobs = str(obs)
msselargs['observation'] = sobs
if timerange:
msselargs['time'] = timerange
if scan:
msselargs['scan'] = scan
else:
warnstr='Data selection by '
datasels=[]
if timerange: datasels.append('timerange')
if scan: datasels.append('scan')
if obs: datasels.append('observation')
if len(datasels):
warnstr+=str(datasels)+' will be ignored for usescartch=False'
casalog.post(warnstr,'WARN')
# === skip this use ms.msselect instead ====
# ms.msseltoindex only goes by the subtables - it does NOT check
# whether the main table has any rows matching the selection.
# try:
# selindices = myms.msseltoindex(**msselargs)
# casalog.post("msselargs=",msselargs," selindices=",selindices)
# except Exception, instance:
# casalog.post('nselrowscore exception: %s' % instance,'SEVERE')
# raise instance
# query = []
# if field:
# query.append("FIELD_ID in " + str(selindices['field'].tolist()))
# if spw:
# query.append("DATA_DESC_ID in " + str(selindices['spw'].tolist()))
# if obs and usescratch:
# query.append("OBSERVATION_ID in " + str(selindices['obsids'].tolist()))
# I don't know why ms.msseltoindex takes a time argument
# - It doesn't seem to appear in the output.
# if scan and usescratch:
# query.append("SCAN_NUMBER in " + str(selindices['scan'].tolist()))
#if timerange and usescratch:
# query.append("TIME in [select TIME where
# for intent (OBS_MODE in STATE subtable), need a subquery part...
# if intent:
# query.append("STATE_ID in [select from ::STATE where OBS_MODE==pattern(\""+\
# intent+"\") giving [ROWID()]]")
# casalog.post("query=",query)
# mytb = tbtool()
# mytb.open(vis)
myms = ms()
casalog.post(str(msselargs))
myms.open(vis)
if (len(msselargs)==0):
retval = myms.nrow()
myms.close()
else:
try:
# st = mytb.query(' and '.join(query),style='python') # Does style matter here?
# retval = st.nrows()
# st.close() # needed to clear tablecache?
# mytb.close()
myms.msselect(msselargs)
retval = myms.nrow(True)
finally:
#if ismms:
# casalog.post('nselrows: %s' % instance,'WARN')
#else:
# this is probably redundant as the exception will be handle in setjy()
#casalog.post('nselrowscore exception: %s' % instance,'SEVERE')
# pass
myms.close()
#if not ismms:
# raise Exception, instance
#else:
# casalog.post('Proceed as it appears to be dealing with a MMS...','DEBUG')
return retval
def parse_fluxdict(fluxdict, vis, field='', spw='', observation='', timerange='', scan='', intent='', usescratch=None):
"""
Parser function for fluxdict (dictionary output from fluxscale) to set
fluxdensity, spix, and reffreq parameters (as in 'manual' mode)
"""
# set spix and reffreq if there
myms = ms()
# dictionary for storing modified (taking AND between the data selection
# and fluxdict content)
modmsselargs={}
msselargs={}
msselargs['field']=field
msselargs['spw']=spw
msselargs['scanintent']=intent
# only applicable for usescratch=T
if usescratch:
if observation:
msselargs['observation'] = observation
if timerange:
msselargs['time'] = timerange
if scan:
msselargs['scan'] = scan
try:
myms.open(vis)
myms.msselect(msselargs)
#selindices = myms.msselectedindices()
msselfldids=myms.range(["FIELD_ID"])['field_id']
except Exception as instance:
casalog.post('parse_fluxdict exception: %s' % instance,'SEVERE')
raise
finally:
myms.close()
# check if fluxdict is valid
if fluxdict=={}:
raise Exception("fluxdict is empty")
else:
msg=""
if not "freq" in fluxdict:
msg+="freq "
if not "spwID" in fluxdict:
msg+="spwID "
if len(msg):
raise Exception("Input fluxdict is missing keywords:"+msg)
# select fields only common to the dictionary and field selection
fieldids=[]
for ky in fluxdict.keys():
try:
int(ky)
fieldids.append(ky) # list in string field ids
except:
pass
if not len(fieldids):
casalog.post('No field ids was found in the dictionary given for fluxdict. Please check the input.', 'SEVERE')
#if len(selindices['field']):
if len(msselfldids):
# also check there is common field id, note field ids in selindices are in numpy.array
#selfieldids = [fd for fd in fieldids if int(fd) in selindices['field'].tolist()]
selfieldids = [fd for fd in fieldids if int(fd) in msselfldids]
if not len(selfieldids):
raise Exception("No field was found in fluxdict for the given data selection")
else:
selfieldids = fieldids
return selfieldids