##################################################################
# v 1.01 - 2011 Nov 06
##################################################################
# Task to process wideband continuum EVLA datasets
#
# v 1.01: added parameter "doplot" to control whether plots of data
# and calibration products are produced. This is in response to the
# current plotcal bug (now ~15 x slower than before) as well as
# display-back issues for plotms
#
# Order of operations:
#
# -> Import data
# -> Check for missing scans
# -> Check for zeros; plot
# -> Apply on-line flags, zero flagging, shadow flagging
# -> Plots of raw data for flux calibrator
# -> Calibrate & plot
# -> Apply calibration & plot
# -> Image
#
# to-do:
# - somehow, sniff for configuration -- if D, need to image
# more primary beams, and perhaps time-averaging needs to take this
# into account for the calculation of time-average smearing
# - fix ":" in HTML files
# - calculate dynamic range using peak / rms instead of fit at phase
# center
# - check time-binning of data to decide on solution lengths
# - get flux at peak of image as well as fitting to phase center
# - add option of using nterms > 1
# - channel-average data
# - import calibrators separately and merge MSs after calibration (?)
# - add refant choosing code: be sure it's config-dependent
# - make padding of online flags integration-time dependent
# - make calibrated channels observation-specific (not hard-coded);
# perhaps dependent on detected RMS noise/channel
# - add histogram / rms plots
import os, time
import numpy as np
import pylab as pl
from taskinit import *
from listsdm_cli import listsdm_cli as listsdm
from listobs_cli import listobs_cli as listobs
from importevla_cli import importevla_cli as importevla
from flagmanager_cli import flagmanager_cli as flagmanager
from flagdata_cli import flagdata_cli as flagdata
from flagcmd_cli import flagcmd_cli as flagcmd
from plotms_cli import plotms_cli as plotms
from plotcal_cli import plotcal_cli as plotcal
from bandpass_cli import bandpass_cli as bandpass
from gaincal_cli import gaincal_cli as gaincal
from fluxscale_cli import fluxscale_cli as fluxscale
from plotants_cli import plotants_cli as plotants
from applycal_cli import applycal_cli as applycal
from clean_cli import clean_cli as clean
import casac
from tasks import *
from math import pi,floor
#=====================================================================
# Task inputs
def mkpipeline(mode,
msname,
sdmname,
rootname,
band,
dummy,
checkzeros,
timeave,
flaglist,
doplot,
calimage,
email,
webdir,
http):
version = "v 1.01, 6 November 2011"
timeStart = time.time()
#==> Check that needed parameters are set:
if mode == 'initproc':
if rootname == '':
casalog.post("ERROR: you must specify a root name for output files", \
priority='ERROR')
return
if band == '':
casalog.post("ERROR: you must specify a band for processing", \
priority='ERROR')
return
if mode == 'reproc':
if msname == '':
casalog.post("ERROR: you must specify an MS for reprocessing", \
priority='ERROR')
return
#==> Remove trailing '/' which will cause problems later on:
if mode == 'initproc':
sdmname = sdmname.strip('/')
if mode == 'reproc':
msname = msname.strip('/')
#==> Make a time string for later use:
timeArr = time.localtime()
timeStr = '%4.0f-%02.0f-%02.0f.%02.0f:%02.0f' % (timeArr[0], timeArr[1], timeArr[2], timeArr[3], timeArr[4])
#==> Definitions to use later:
if mode == 'initproc':
plotDir = rootname + '.mkpipeline'
msName = rootname + '.ms'
repro = False
else:
plotDir = msname + '.mkpipeline.' + timeStr
msName = msname
rootname = msname
repro = True
#==> Log file handling:
initLog = casalog.logfile()
if repro:
logName = rootname + '.mkpipeline.' + timeStr + '.log'
casalog.post("Running mkpipeline " + version)
casalog.post("See %s for ongoing log output" % logName)
casalog.setlogfile(logName)
else:
logName = rootname + '.mkpipeline.log'
casalog.post("Running mkpipeline " + version)
casalog.post("See %s for ongoing log output" % logName)
casalog.setlogfile(logName)
casalog.origin('mkpipeline')
casalog.post("Task inputs:")
casalog.post("mode = " + mode)
if mode == 'initproc':
casalog.post("rootname = " + rootname)
casalog.post("sdmname = " + sdmname)
casalog.post("band = " + band)
casalog.post("dummy = " + str(dummy))
casalog.post("checkzeros = " + str(checkzeros))
casalog.post("timeave = " + str(timeave))
if mode == 'reproc':
casalog.post("msname = " + msname)
casalog.post("doplot = " + str(doplot))
casalog.post("calimage = " + str(calimage))
casalog.post("email = " + email)
casalog.post("webdir = " + webdir)
casalog.post("http = " + http)
casalog.post(" ")
# Reference antenna to use
refAnt='ea21'
#==> Given band, retrieve other needed information from SDM:
if mode == 'initproc':
infoDict = getSDMInfo(sdmName=sdmname, band=band, dummy=dummy, \
calImage=calimage)
else:
infoDict = getMSInfo(msName=msName, calImage=calimage)
band = infoDict['band']
if not infoDict:
return
#==> Make output plot directory & check
os.system('mkdir '+plotDir)
if not os.path.isdir(plotDir):
casalog.post("mkpipeline: Unable to create output plot directory -- ", priority='ERROR')
casalog.post("mkpipeline: check permissions.", priority='ERROR')
return
#==> Import the data
if mode == 'initproc':
importData(sdmName=sdmname, msName=msName, timeStr=timeStr,
baseBand=infoDict['baseBand'], plotDir=plotDir,
scanStr=infoDict['scanStr'],
checkZeros=checkzeros)
bBandOrig = infoDict['baseBand']
# **** Perhaps should put this after the time-average step?
#==> Plot the raw data
if ((mode == 'initproc') and (doplot)):
plotRawData(msName=msName, baseBand=infoDict['baseBand'],
plotScan=infoDict['plotScan'], refAnt=refAnt,
plotDir=plotDir)
#==> Time-average the data, if requested
if ((mode == 'initproc') and (timeave)):
timeAveDict = timeAverage(msName=msName, band=band,
baseBand=infoDict['baseBand'],
msRoot=rootname)
# Time-averaging was not successful (zeros were produced):
if timeAveDict['success'] == 0:
return
msName = timeAveDict['msAve']
infoDict['baseBand'] = timeAveDict['baseBand']
infoDict['SPWs'] = timeAveDict['SPWs']
## #== *** This is a hack to be removed later! For now, need to flag
## # subband 0 lower edges or bandpass fails completely.
## if mode == 'initproc':
## for bBand in infoDict['baseBand']:
## testBand = bBand.find('~')
## if (testBand > -1):
## (bLow, bHigh) = str.split(bBand, '~')
## spwStr = bLow+':0~5'
## flagdata(vis=msName, spw=spwStr, flagbackup=False)
#==> Apply additional flags, if requested:
# *** put in its own subroutine & add timing
if flaglist:
flagSave = 'before_mf_' + timeStr
flagmanager(vis=msName, mode='save', versionname=flagSave)
flagcmd(vis=msName, flagmode='cmd', command=flaglist, flagbackup=False)
#==> If this is a reprocessing run, clear calibration first:
if mode == 'reproc':
clearcal(msName)
#==> Calibrate the data
calTables = calibData(msName=msName, rootName=rootname,
fluxFields=infoDict['fluxSrc'],
baseBand=infoDict['baseBand'],
spwCalib=infoDict['SPWs'], doPlot=doplot,
plotDir=plotDir, useGC=infoDict['useGC'],
calFields=infoDict['calSrc'], band=band,
refAnt=refAnt, repro=repro, timeStr=timeStr,
nChan=infoDict['nChan'])
#==> Apply the calibration
applyCalib(msName=msName, clearPrev=False, rootName=rootname,
# calChan='0~63', calTables=calTables, useGC=infoDict['useGC'],
calTables=calTables, useGC=infoDict['useGC'],
spwCalib=infoDict['SPWs'], plotField=infoDict['calSrc'],
timeBin='10s', chanBin='5', doCal=1, timeStr=timeStr,
doPlot=1, plotDir=plotDir, targetField=infoDict['targetSrc'],
fluxFields=infoDict['fluxSrc'])
#==> Image the data
mkImages(msName=msName, rootName=rootname, repro=repro,
imageFields=infoDict['imageSrc'], timeStr=timeStr,
plotDir=plotDir, baseBand=infoDict['baseBand'])
timeEnd = time.time()
timeTotal = timeEnd - timeStart
casalog.origin('mkpipeline')
casalog.post("Complete mkpipeline run took %0.2f seconds (%0.2f minutes)" % (timeTotal, timeTotal/60))
casalog.post(" ")
#==> Copy the CASA log information to appropriately-named files
os.system("grep 'mkpipeline::' " + logName + " >> " + plotDir + "/mkpipeline.log")
os.system('mv %s %s/' % (logName, plotDir))
casalog.setlogfile(initLog)
#==> Write the summary web page:
mkWebPage(msName=msName, plotDir=plotDir, logName=logName,
plotScan=infoDict['plotScan'], calFields=infoDict['calSrc'],
baseBand=infoDict['baseBand'], checkZeros=checkzeros,
bBandOrig=bBandOrig, fluxFields=infoDict['fluxSrc'],
repro=repro, doPlot=doplot)
#==> Copy plots to web space, if requested:
if webdir:
os.system('cp -rp ' + plotDir + ' ' + webdir)
#==> Send an email notification, if requested:
if email:
subject = 'mkpipeline complete, plots made: %s' % msName
if http:
os.system('echo "Finished processing %s to %s; see %s/%s" > tempmail' % (sdmname, msName, http, plotDir+'/index.html'))
else:
os.system('echo "Finished processing %s to %s; see %s" > tempmail' % (sdmname, msName, plotDir+'/index.html'))
os.system('cat tempmail | mail -s "%s" %s' % (subject, email))
os.system('rm tempmail')
########################################################################
# Get information needed for processing from the SDM files
def getSDMInfo(sdmName, band, dummy, calImage):
# EVLA bands - ultimately, the band names should be in the SDM.
# This is not currently the case, so we must make some
# assumptions about the frequency ranges, which may overlap.
#
# Last boolean controls whether or not gain curves should be
# used during calibration.
bandDict = { 'L': [0.9, 2.0, False],
'S': [2.0, 4.0, False],
'C': [4.0, 8.0, True],
'X': [8.0, 12.0, True],
'Ku': [12.0, 18.0, True],
'K': [18.0, 26.5, True],
'Ka': [26.5, 40.0, True],
'Q': [40.0, 50.0, True] }
scanList = listsdm(sdmName)
if (dummy):
scan1 = scanList.pop(1)
scan2 = scanList.pop(2)
myFreqRange = bandDict[band]
freqMin = myFreqRange[0]*1e9
freqMax = myFreqRange[1]*1e9
scanStr = ''
bpassSrc = []
gainSrc = []
targetSrc = []
fluxSrc = []
spwObs = []
nChan = []
for scan in scanList:
# Use the midpoint of the AC frequency range to avoid band
# edges:
testEl = int(len(scanList[scan]['reffreq'][0])/4)
if ((scanList[scan]['reffreq'][0][testEl] > freqMin) and
(scanList[scan]['reffreq'][0][testEl] < freqMax) and
(string.find(scanList[scan]['intent'], 'POINTING') == -1)):
scanStr = scanStr + str(scan) + ','
spwObs.extend(scanList[scan]['spws'][0])
nChan.extend(scanList[scan]['nchan'][0])
if (string.find(scanList[scan]['intent'], 'BANDPASS') > -1):
# Will plot the last bandpass scan
plotScan = str(scan)
bpassSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'PHASE') > -1):
gainSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'TARGET') > -1):
targetSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'AMPLI') > -1):
fluxSrc.append(scanList[scan]['field'])
if scanStr:
# Check that needed sources are present, make unique lists:
if not gainSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no gain calibration sources present. This could be due to lack of proper SCAN_INTENT definition.", priority='ERROR')
return
else:
gainSrc = list(np.unique(gainSrc))
if not bpassSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no bandpass calibration sources present. This could be due to lack of proper SCAN_INTENT definition.", priority='ERROR')
return
else:
bpassSrc = list(np.unique(bpassSrc))
if not targetSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no science target scans are present. This could be due to lack of proper SCAN_INTENT definition. Proceeding with only calibration sources", priority='WARNING')
else:
targetSrc = list(np.unique(targetSrc))
if not fluxSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no flux calibration sources present. Proceeding with calibration; will not perform flux scaling.", priority='WARNING')
else:
fluxSrc = list(np.unique(fluxSrc))
# Need to get rid of trailing ','
scanStr = str.rstrip(scanStr, ',')
spwObs = list(np.unique(spwObs))
# If there are spectral windows with varying numbers of channels, fail:
nChan = np.unique(nChan)
if nChan.shape[0] > 1:
casalog.origin('mkpipeline')
casalog.post("ERROR: Spectral windows have varying numbers of channels; cannot proceed.", priority='ERROR')
return
else:
nChan = nChan[0]
# Determine baseband SPW sets
spwObs.sort()
baseband1 = spwObs[0:len(spwObs)/2]
baseband2 = spwObs[len(spwObs)/2:len(spwObs)]
bb1str = str(min(baseband1))+'~'+str(max(baseband1))
bb2str = str(min(baseband2))+'~'+str(max(baseband2))
# Image calibrators if requested; otherwise, just image the
# science target:
if calImage:
tmp1 = np.append(bpassSrc,gainSrc)
tmp2 = np.append(tmp1,targetSrc)
imageSrc = list(np.unique(np.append(tmp2,fluxSrc)))
else:
imageSrc = list(np.unique(targetSrc))
# Make a list of all the calibrators:
tmp1 = np.append(bpassSrc,gainSrc)
if fluxSrc:
calSrc = list(np.unique(np.append(tmp1,fluxSrc)))
else:
calSrc = list(np.unique(tmp1))
# Post the sources to the logger:
casalog.origin('mkpipeline')
casalog.post("Bandpass calibrator fields: " + str(bpassSrc))
casalog.post("Gain calibrator fields: " + str(gainSrc))
casalog.post("Flux calibrator fields: " + str(fluxSrc))
casalog.post("All calibrator fields: " + str(calSrc))
casalog.post("Science target fields: " + str(targetSrc))
casalog.post("Fields to image: " + str(imageSrc))
casalog.post(" ")
infoDict = { 'SPWs': spwObs,
'nChan': nChan,
'baseBand': [bb1str, bb2str],
'useGC': bandDict[band][2],
'scanStr': scanStr,
'plotScan': plotScan,
'bpassSrc': bpassSrc,
'gainSrc': gainSrc,
'fluxSrc': fluxSrc,
'calSrc': calSrc,
'targetSrc': targetSrc,
'imageSrc': imageSrc }
return(infoDict)
else:
casalog.origin('mkpipeline')
casalog.post("There are no scans present for %s-band" % band, priority='ERROR')
return
########################################################################
# Get information needed for processing from the MS files. Note that
# the MS must have the verbatim xml files from the SDM directory.
def getMSInfo(msName, calImage):
# EVLA bands - ultimately, the band names should be in the SDM.
# This is not currently the case, so we must make some
# assumptions about the frequency ranges, which may overlap.
#
# Last boolean controls whether or not gain curves should be
# used during calibration.
bandDict = { 'L': [0.9, 2.0, False],
'S': [2.0, 4.0, False],
'C': [4.0, 8.0, True],
'X': [8.0, 12.0, True],
'Ku': [12.0, 18.0, True],
'K': [18.0, 26.5, True],
'Ka': [26.5, 40.0, True],
'Q': [40.0, 50.0, True] }
# Return the SPW range:
tb.open(msName + '/SPECTRAL_WINDOW')
refFreq = tb.getcol('REF_FREQUENCY')
nChan = tb.getcol('NUM_CHAN')
spwObs = range(0, len(refFreq))
baseband1 = spwObs[0:len(spwObs)/2]
baseband2 = spwObs[len(spwObs)/2:len(spwObs)]
bb1str = str(min(baseband1))+'~'+str(max(baseband1))
bb2str = str(min(baseband2))+'~'+str(max(baseband2))
# If there are spectral windows with varying numbers of channels, fail:
nChan = np.unique(nChan)
if nChan.shape[0] > 1:
casalog.origin('mkpipeline')
casalog.post("ERROR: Spectral windows have varying numbers of channels; cannot proceed.", priority='ERROR')
return
else:
nChan = nChan[0]
# Determine the band:
freqMinGHz = min(refFreq)/1e9
freqMaxGHz = max(refFreq)/1e9
freqMid = freqMinGHz + (freqMaxGHz - freqMinGHz) / 2
for bandName in bandDict:
if (freqMid > bandDict[bandName][0]) and (freqMid < bandDict[bandName][1]):
band = bandName
break
freqMin = min(refFreq)
freqMax = max(refFreq)
bpassSrc = []
gainSrc = []
targetSrc = []
fluxSrc = []
scanList = listsdm(msName)
for scan in scanList:
# Use the midpoint of the frequency range to avoid band
# edges:
testEl = int(len(scanList[scan]['reffreq'][0])/4)
if ((scanList[scan]['reffreq'][0][testEl] > freqMin) and
(scanList[scan]['reffreq'][0][testEl] < freqMax) and
(string.find(scanList[scan]['intent'], 'POINTING') == -1)):
if (string.find(scanList[scan]['intent'], 'BANDPASS') > -1):
bpassSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'PHASE') > -1):
gainSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'TARGET') > -1):
targetSrc.append(scanList[scan]['field'])
if (string.find(scanList[scan]['intent'], 'AMPLI') > -1):
fluxSrc.append(scanList[scan]['field'])
# Check that needed sources are present, make unique lists:
if not gainSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no gain calibration sources present. This could be due to lack of proper SCAN_INTENT definition.", priority='ERROR')
return
else:
gainSrc = list(np.unique(gainSrc))
if not bpassSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no bandpass calibration sources present. This could be due to lack of proper SCAN_INTENT definition.", priority='ERROR')
return
else:
bpassSrc = list(np.unique(bpassSrc))
if not targetSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no science target scans are present. This could be due to lack of proper SCAN_INTENT definition. Proceeding with only calibration sources", priority='WARNING')
else:
targetSrc = list(np.unique(targetSrc))
if not fluxSrc:
casalog.origin('mkpipeline')
casalog.post("ERROR: no flux calibration sources present. Proceeding with calibration; will not perform flux scaling.", priority='WARNING')
else:
fluxSrc = list(np.unique(fluxSrc))
# Image calibrators if requested; otherwise, just image the
# science target:
if calImage:
tmp1 = np.append(bpassSrc,gainSrc)
tmp2 = np.append(tmp1,targetSrc)
imageSrc = list(np.unique(np.append(tmp2,fluxSrc)))
else:
imageSrc = list(np.unique(targetSrc))
# Make a list of all the calibrators:
tmp1 = np.append(bpassSrc,gainSrc)
if fluxSrc:
calSrc = list(np.unique(np.append(tmp1,fluxSrc)))
else:
calSrc = list(np.unique(tmp1))
# Post the sources to the logger:
casalog.origin('mkpipeline')
casalog.post("Bandpass calibrator fields: " + str(bpassSrc))
casalog.post("Gain calibrator fields: " + str(gainSrc))
casalog.post("Flux calibrator fields: " + str(fluxSrc))
casalog.post("All calibrator fields: " + str(calSrc))
casalog.post("Science target fields: " + str(targetSrc))
casalog.post("Fields to image: " + str(imageSrc))
casalog.post(" ")
infoDict = { 'SPWs': spwObs,
'nChan': nChan,
'baseBand': [bb1str, bb2str],
'useGC': bandDict[band][2],
'band': band,
'plotScan': '',
'bpassSrc': bpassSrc,
'gainSrc': gainSrc,
'fluxSrc': fluxSrc,
'calSrc': calSrc,
'targetSrc': targetSrc,
'imageSrc': imageSrc }
return(infoDict)
########################################################################
# Import the data
def importData(sdmName, timeStr, msName,
baseBand, scanStr, plotDir, checkZeros):
casalog.origin('mkpipeline')
casalog.post("Importing %s to MS format..." % (sdmName))
timeImportStart = time.time()
importevla(asdm=sdmName, vis=msName, applyflags=False, \
switchedpower=True, scans=scanStr, flagzero=False, \
shadow=False)
# Copy over the XML tables needed for listsdm
xmlTab = ['Scan.xml', 'Main.xml', 'ConfigDescription.xml', 'DataDescription.xml', 'SpectralWindow.xml', 'Field.xml', 'Antenna.xml', 'Station.xml', 'ExecBlock.xml']
for xmlFile in xmlTab:
os.system('cp ' + sdmName + '/' + xmlFile + ' ' + msName)
timeImportEnd = time.time()
timeImportTotal = timeImportEnd - timeImportStart
casalog.origin('mkpipeline')
casalog.post("Importing data took %0.2f seconds (%0.2f minutes)" % (timeImportTotal, timeImportTotal/60))
casalog.post(" ")
# Plot antenna distribution and choose a reference antenna
## casalog.origin('mkpipeline')
## casalog.post("Selecting a reference antenna...")
## timeRefAntStart = time.time()
pl.clf()
plotants(vis=msName, figfile=plotDir+'/antpos.png')
## refAnt = chooseRefAnt(msName)
## timeRefAntEnd = time.time()
## timeRefAntTotal = timeRefAntEnd - timeRefAntStart
## casalog.origin('mkpipeline')
## casalog.post("Selecting a reference antenna took %0.2f seconds (%0.2f minutes)" % (timeRefAntTotal, timeRefAntTotal/60))
## casalog.post(" ")
## return(refAnt)
# Check for and plot zeros, if requested:
if checkZeros:
casalog.origin('mkpipeline')
casalog.post("Checking for zeros...")
timeZeroCheckStart = time.time()
findZeros(sdmName=sdmName, timeStr=timeStr, msName=msName, plotDir=plotDir)
timeZeroCheckEnd = time.time()
timeZeroCheckTotal = timeZeroCheckEnd - timeZeroCheckStart
casalog.origin('mkpipeline')
casalog.post("Checking for and plotting zeros took %0.2f seconds (%0.2f minutes)" % (timeZeroCheckTotal, timeZeroCheckTotal/60))
casalog.post(" ")
## # Flag antenna(s) with known problems
## if antenna:
## flagdata(vis=msName, mode='manualflag', selectdata=True, \
## antenna=antenna, flagbackup=False)
# Check for missing scans
# Note that this check is not meaningful if only a subset of the
# SDM file was imported
casalog.origin('mkpipeline')
casalog.post("Checking for missing scans...")
checkScans(msName=msName)
# Import the online flags, make plots, apply flags, and save a flag
# version:
casalog.origin('mkpipeline')
casalog.post("Importing and applying online flags...")
timeOnlineFlagStart = time.time()
plotName = plotDir + '/onlineFlags.png'
# plot flags
flagcmd(vis=msName, flagmode='xml', optype='plot', outfile=plotName)
# apply online flags, 0.5 s padding (should really be tint / 2).
flagcmd(vis=msName, flagmode='xml', tbuff=0.5, optype='apply', \
flagbackup=False)
timeOnlineFlagEnd = time.time()
timeOnlineFlagTotal = timeOnlineFlagEnd - timeOnlineFlagStart
casalog.origin('mkpipeline')
casalog.post("Applying online flags took %0.2f seconds (%0.2f minutes)" % (timeOnlineFlagTotal, timeOnlineFlagTotal/60))
casalog.post(" ")
# Apply zero and shadow flags; save flags:
casalog.origin('mkpipeline')
casalog.post("Applying zero and shadow flags...")
timeFlagStart = time.time()
flagdata(vis=msName, mode='shadow', flagbackup=False)
flagdata(vis=msName, mode='manualflag', clipminmax=[0.0, 1E-10], \
clipoutside=False, clipcolumn='DATA', \
clipexpr=['ABS RR','ABS LL','ABS RL','ABS LR'], \
flagbackup=False)
flagmanager(vis=msName, mode='save', \
versionname='flags.shadow.zeros_'+timeStr, \
comment='Flag save after online, shadow, and zero flags applied')
timeFlagEnd = time.time()
timeFlagTotal = timeFlagEnd - timeFlagStart
casalog.origin('mkpipeline')
casalog.post("Applying zero and shadow flags took %0.2f seconds (%0.2f minutes)" % (timeFlagTotal, timeFlagTotal/60))
casalog.post(" ")
########################################################################
# Time-average the data
def timeAverage(msName, band, baseBand, msRoot):
casalog.origin('mkpipeline')
casalog.post("Time-averaging data...")
timeTimeAveStart = time.time()
tb.open(msName+'/ANTENNA')
antPos = tb.getcol('POSITION')
tb.close()
antDist = []
for antenna1 in range(0,antPos.shape[1]):
xRef = antPos[0][antenna1]
yRef = antPos[1][antenna1]
zRef = antPos[2][antenna1]
for antenna2 in range(antenna1+1,antPos.shape[1]):
x = antPos[0][antenna2]
y = antPos[1][antenna2]
z = antPos[2][antenna2]
antDist.append(np.sqrt((xRef-x)**2.0+(yRef-y)**2.0+(zRef-z)**2.0))
# Maximum distance between antennas (m)
maxDist = max(antDist)
# Tolerance for time-average smearing at the edge of the image
# field (assume this is 1 primary beam), fraction of ideal response:
Rtol = 0.01
# Ratio of synthesized beam to primary beam, assuming
#
# thetaSB(rad) = 1.2 * lambda / Dmax(m)
# thetaPB(arcmin) = 45 / nu(GHz)
#
# where Dmax(m) is the max. distance between antennas calculated
# gives thetaSB / thetaPB = 27.5 / Dmax(m):
thetaRat = 27.5 / maxDist
# From Chapter 18 of Synthesis Imaging in Radio Astronomy II, find
# an acceptable average time tau; round down the the nearest
# integer number of seconds:
tau = int(np.sqrt(1e9 * Rtol) * thetaRat)
casalog.post("Assuming a 1% tolerance for time-average smearing,")
casalog.post(" can time-average to %i seconds." % tau)
casalog.post(" ")
# If the observation is high-frequency, don't want to average too
# much (consider atmospheric effects); for low-frequency don't
# want to average beyond the timescale of possible instrumental
# variations. (Could improve this in the future, perhaps, by
# using information in the weather table.)
bandAve = { 'L': 60,
'S': 60,
'C': 60,
'X': 30,
'Ku': 20,
'K': 10,
'Ka': 5,
'Q': 5 }
tauBand = bandAve[band];
casalog.post("Observations were performed in %s-band, so recommended" % band)
casalog.post(" maximum time bin is %i seconds." % tauBand)
casalog.post(" ")
# Use the shorter of these two:
tauApply = min(tau, tauBand)
tauStr = '%is' % tauApply
casalog.post("Averaging data using %i-second bins." % tauApply)
# Average the data:
msAve = msRoot + '.' + str(tauApply) + 's.ms'
# Explicitly set the SPW range to avoid a bug related to
# gain curve selection (CAS-2705):
spwStr = ','.join(baseBand)
split(vis=msName, outputvis=msAve, datacolumn='data', timebin=tauStr, \
keepflags=False, spw=spwStr)
# Check that split didn't write all zeros to a column:
ms.open(msAve)
testRange = ms.range(['u','v','w','sigma','weight'])
# Takes a lot of time to check whether the data are zeros.
# Just test the uvw, sigma, and weight columns for now.
## dataRange = ms.range(['real'])
ms.close
dataOK = 1
if (((max(testRange['u']) == 0) and (max(testRange['v']) == 0) and (max(testRange['w']) == 0))):
casalog.post("All UVW values in " + msAve + " are zero.", \
priority='ERROR')
dataOK = 0
if (max(testRange['sigma']) == 0):
casalog.post("All SIGMA values in " + msAve + " are zero.", \
priority='ERROR')
dataOK = 0
if (max(testRange['weight']) == 0):
casalog.post("All WEIGHT values in " + msAve + " are zero.", \
priority='ERROR')
dataOK = 0
if (dataOK == 0):
aveDict = { 'success': 0 }
return aveDict
# Copy over the XML tables needed for listsdm
xmlTab = ['Scan.xml', 'Main.xml', 'ConfigDescription.xml', 'DataDescription.xml', 'SpectralWindow.xml', 'Field.xml', 'Antenna.xml', 'Station.xml', 'ExecBlock.xml']
for xmlFile in xmlTab:
os.system('cp ' + msName + '/' + xmlFile + ' ' + msAve)
# Return a corrected version of the SPW range:
tb.open(msAve + '/SPECTRAL_WINDOW')
newSPW = tb.getcol('REF_FREQUENCY')
spwObs = range(0, len(newSPW))
baseband1 = spwObs[0:len(spwObs)/2]
baseband2 = spwObs[len(spwObs)/2:len(spwObs)]
bb1str = str(min(baseband1))+'~'+str(max(baseband1))
bb2str = str(min(baseband2))+'~'+str(max(baseband2))
timeTimeAveEnd = time.time()
timeTimeAveTotal = timeTimeAveEnd - timeTimeAveStart
casalog.origin('mkpipeline')
casalog.post("Time-averaging data took %0.2f seconds (%0.2f minutes)" % (timeTimeAveTotal, timeTimeAveTotal/60))
casalog.post(" ")
aveDict = { 'msAve': msAve,
'tauApply': tauApply,
'baseBand': [bb1str, bb2str],
'SPWs': spwObs,
'success': 1 }
return aveDict
########################################################################
# Plot the raw data
def plotRawData(msName, baseBand, plotScan, refAnt, plotDir):
casalog.origin('mkpipeline')
casalog.post("Plotting uncalibrated phase and amplitude data...")
timePlotStart = time.time()
# Only plot antennas for which there is data:
tb.open(msName + '/ANTENNA')
antArr = tb.getcol('NAME')
tb.close()
for i in range(0, len(baseBand)):
spwStr = baseBand[i]
html = open(plotDir+'/scan_%s.amp-phase_v_freq.spw%s.html' % (plotScan, spwStr), 'w')
htmlHead = '''
Phase/amplitude vs. freq for '''+msName+''', scan '''+plotScan+''', spw '''+spwStr+'''
| Amplitude vs. Frequency: |
Phase vs. Frequency: |
'''
html.write(htmlHead)
for j in range(0, len(antArr)):
plotMade = False
if antArr[j] == refAnt:
continue
antStr = '%s&%s' % (antArr[j], refAnt)
plotNameAmp = plotDir + '/' + antStr + '.spw' + spwStr + '.freq.amp.png'
plotNamePh = plotDir + '/' + antStr + '.spw' + spwStr + '.freq.phase.png'
htmlNameAmp = antStr + '.spw' + spwStr + '.freq.amp.png'
htmlNamePh = antStr + '.spw' + spwStr + '.freq.phase.png'
htmlStr = '''
%s:
 |
|
''' % (antStr, htmlNameAmp, htmlNameAmp, htmlNamePh, htmlNamePh)
html.write(htmlStr)
# Make amplitude vs. frequency plots:
plotms(vis=msName, xaxis='freq', yaxis='amp', selectdata=True, \
## spw=spwStr, antenna=antStr, avgtime='5s', \
spw=spwStr, antenna=antStr, \
plotfile=plotNameAmp, correlation='RR,LL', scan=plotScan, \
coloraxis='corr')
while (plotMade == False):
time.sleep(1)
print "Testing for plot %s..." % plotNameAmp
plotMade = os.path.isfile(plotNameAmp)
print "Plot not found."
plotMade = False
# Make phase vs. frequency plots:
plotms(vis=msName, xaxis='freq', yaxis='phase', selectdata=True, \
## spw=spwStr, antenna=antStr, avgtime='5s', \
spw=spwStr, antenna=antStr, \
plotfile=plotNamePh, correlation='RR,LL', \
scan=plotScan, plotrange=[0,0,-180,180], \
coloraxis='corr')
while (plotMade == False):
time.sleep(1)
print "Testing for plot %s..." % plotNamePh
plotMade = os.path.isfile(plotNamePh)
print "Plot not found."
htmlTail = '''
'''
html.write(htmlTail)
html.close()
timePlotEnd = time.time()
timePlotTotal = timePlotEnd - timePlotStart
casalog.origin('mkpipeline')
casalog.post("Plotting uncalibrated data took %0.2f seconds (%0.2f minutes)" % (timePlotTotal, timePlotTotal/60))
casalog.post(" ")
########################################################################
# Calibrate data
#
# msName = name of ms to be calibrated
# fluxFields = list of field IDs for flux calibrator, e.g. [2,3,4]
# refAnt = desired reference antenna
# baseBand = array of spw strings to plot together (oft. basebands),
# e.g. ['2~9', '10~17']
# plotDir = output plotting directory
# rootName = root for output filenames
# band = observing band (e.g., 'C' or 'Ku')
def calibData(msName, fluxFields, calFields, nChan,
refAnt, spwCalib, baseBand, useGC, doPlot,
plotDir, rootName, band, repro, timeStr):
#=====================================================================
# Set up useful bits for later
# Look for the model image directory
casaPath = os.environ.get('CASAPATH').split()[0]
modImagePath = casaPath + '/data/nrao/VLA/CalModels/'
if not os.path.isdir(modImagePath):
casalog.origin('mkpipeline')
casalog.post("ERROR: The calibrator model image directory not present at "+modImagePath, priority='ERROR')
return
# Get antenna information
tb.open(msName + '/ANTENNA')
antArr = tb.getcol('NAME')
tb.close()
#=====================================================================
# Determine calibration intervals based on observation information
# Set up dictionary of phase and amplitude solution intervals,
# dependent on observing band. Order is:
# bpassPhInt, gainPhInt, gainAmpInt, gainPhaseSmo, gainCurve
solIntDict = { 'L': ['5s', '20s', 'inf', 40],
'S': ['5s', '20s', 'inf', 40],
'C': ['5s', '20s', 'inf', 40],
'X': ['5s', '15s', 'inf', 30],
'Ku': ['5s', '10s', 'inf', 20],
'K': ['5s', '10s', 'inf', 20],
'Ka': ['5s', '10s', 'inf', 20],
'Q': ['5s', '10s', 'inf', 20] }
# Solution interval and channel range for initial phase solutions:
bpassPhInt = solIntDict[band][0]
# Use middle 10% of channels for the initial phase solution:
phChanNum = int(nChan * 0.1)
phChanMin = int((nChan / 2) - (phChanNum / 2))
phChanMax = int((nChan / 2) + (phChanNum / 2) - 1)
bpassPhChan = str(phChanMin) + '~' + str(phChanMax)
# bpassPhChan = '29~34'
# Solution interval and channel range for gain phase solutions,
# as well as smoothing time interval:
gainPhInt = solIntDict[band][1]
gainPhSmo = solIntDict[band][3]
# gainChan = '0~63'
# Solution interval for gain amplitude solutions:
gainAmpInt = solIntDict[band][2]
#=====================================================================
# Determine needed opacity corrections, per spw, and plot the weather
# table:
casalog.origin('mkpipeline')
casalog.post("Determining opacity values and plotting weather table...")
tauMS = plotWX(vis=msName, seasonal_weight=0.5, doPlot=True, plotDir=plotDir)
#=====================================================================
# setjy for any calibrator fields that are standard flux calibrators:
casalog.origin('mkpipeline')
casalog.post("Setting flux scale for any standard flux calibrators...")
timeSetjyStart = time.time()
# Calibrate all spectral windows
tb.open(msName+'/SPECTRAL_WINDOW')
freqList = tb.getcol('REF_FREQUENCY')/1e9
tb.close()
# Get information about the fields
tb.open(msName+'/FIELD')
srcNamesArr = tb.getcol('NAME')
srcIDArrField = tb.getcol('SOURCE_ID')
srcPosArr = tb.getcol('PHASE_DIR')
tb.close()
# Get information about observed spectral windows
tb.open(msName+'/SOURCE')
srcIDArrSource = tb.getcol('SOURCE_ID')
srcSPWArr = tb.getcol('SPECTRAL_WINDOW_ID')
tb.close()
# EVLA bands - ultimately, these should be in the MS and sniffed out
# from there. This is not currently the case, so we must make some
# assumptions about the frequency ranges, which may overlap.
#
# Since there aren't currently Ka images, use K or Q images for this
# band instead.
#
# Since there aren't currently S images, use L or C images for this
# band instead.
#
# Note that this is a rather fragile way of doing things. Need to
# change once band indications are available in the MS.
bandDict = { 'L': [0.9, 3.0],
## 'S': [2.0, 4.0],
'C': [3.0, 8.0],
'X': [8.0, 12.0],
'U': [12.0, 18.0],
'K': [18.0, 34.0],
## 'Ka': [26.5, 40.0],
'Q': [34.0, 50.0] }
# List of positions for flux calibrators:
fluxcalDict = { '3C138': [np.array([80.2911858, 16.6394103])*np.pi/180.0, \
np.array(['L','C','X','U','K','Q'])],
'3C147': [np.array([85.6511025, 49.8519675])*np.pi/180.0, \
np.array(['L','C','X','U','K','Q'])],
'3C286': [np.array([202.7845246, 30.5091444])*np.pi/180.0, \
np.array(['L','C','X','U','K','Q'])],
'3C48': [np.array([24.4220712, 33.1597550])*np.pi/180.0, \
np.array(['L','C','X','U','K','Q'])] }
# Open the MS for writing
im.open(msName, usescratch=True)
# Comparing the positions of flux calibrators in the observation to find
# model images & run setjy, once for each source:
for field in calFields:
srcID = srcIDArrField[field]
spwArr = srcSPWArr[np.where(srcIDArrSource == srcID)]
fluxBandArr = []
# Make a band list that matches the spw list:
for spw in spwArr:
for keys in bandDict:
if (freqList[spw] >= bandDict[keys][0]) and \
(freqList[spw] <= bandDict[keys][1]):
fluxBandArr.append(keys)
break
# RA, Dec of flux calibrator in radians:
RA = np.array(srcPosArr[0][0])[field]
Dec = np.array(srcPosArr[1][0])[field]
srcX = np.cos(RA) * np.cos(Dec)
srcY = np.sin(RA) * np.cos(Dec)
srcZ = np.sin(Dec)
# Check distance from standard calibrators to match source:
for keys in fluxcalDict:
calRA = fluxcalDict[keys][0][0]
calDec = fluxcalDict[keys][0][1]
calX = np.cos(calRA) * np.cos(calDec)
calY = np.sin(calRA) * np.cos(calDec)
calZ = np.sin(calDec)
testDist = np.arccos(srcX*calX + srcY*calY + srcZ*calZ)
# If the source is less than 60 arcsec away, assume it is
# a match and loop over spectral windows to set the flux:
if (testDist < 2.9e-4):
for i in range(0,len(spwArr)):
if ((fluxcalDict[keys][1] == fluxBandArr[i]).any()):
refImgBand = fluxcalDict[keys][1][np.where(fluxcalDict[keys][1] == fluxBandArr[i])]
modImg = modImagePath + keys + '_' + refImgBand[0] + '.im'
im.setjy(field=str(field), spw=str(spwArr[i]),
modimage=modImg, standard='Perley-Butler 2010',
fluxdensity=-1, scalebychan=True)
else:
im.setjy(field=str(field), spw=str(spwArr[i]),
standard='Perley-Butler 2010', fluxdensity=-1,
scalebychan=True)
im.close()
timeSetjyEnd = time.time()
timeSetjyTotal = timeSetjyEnd - timeSetjyStart
casalog.origin('mkpipeline')
casalog.post("Setting flux scale took %0.2f seconds (%0.2f minutes)" % (timeSetjyTotal, timeSetjyTotal/60))
casalog.post(" ")
#=====================================================================
# Bandpass calibration:
casalog.origin('mkpipeline')
casalog.post("Calibrating antenna bandpasses...")
timeBandpassStart = time.time()
if repro:
bpassPh = rootName + '.' + timeStr + '.' + bpassPhInt + '.bcal.phase'
bpassTab = rootName + '.' + timeStr + '.bcal'
else:
bpassPh = rootName + '.' + bpassPhInt + '.bcal.phase'
bpassTab = rootName + '.bcal'
# Loop through spectral windows to create initial phase solutions for
# all bandpass calibration scans:
for n in spwCalib:
zVal = tauMS[n]
bpassSpwStr = '%i:%s' % (n, bpassPhChan)
gaincal(vis=msName, caltable=bpassPh, intent='*BANDPASS*', \
selectdata=True, spw=bpassSpwStr, gaintype='G', \
solint=bpassPhInt, calmode='p', refant=refAnt, \
append=True, minsnr=3.0)
## gaincal(vis=msName, caltable=bpassPh, intent='*BANDPASS*', \
## selectdata=True, spw=bpassSpwStr, gaintype='G', \
## solint=bpassPhInt, calmode='p', refant=refAnt, \
## append=True, \
## minsnr=5.0)
# Loop through spectral windows to create bandpass solutions:
for i in spwCalib:
zVal = tauMS[i]
bpassSpwStr = '%i' % i
bandpass(vis=msName, caltable=bpassTab, gaintable=[bpassPh], \
spw=bpassSpwStr, solint='inf', refant=refAnt, \
combine='scan', solnorm=False, selectdata=True, \
append=True, gaincurve=useGC, opacity=zVal, \
intent='*BANDPASS*', minsnr=10.0)
## bandpass(vis=msName, caltable=bpassTab, gaintable=[bpassPh], \
## spw=bpassSpwStr, solint='inf', refant=refAnt, \
## combine='scan', solnorm=True, selectdata=True, \
## append=True, gaincurve=True, opacity=zVal, \
## intent='*BANDPASS*', minsnr=10.0, \
## bandtype='BPOLY')
timeBandpassEnd = time.time()
timeBandpassTotal = timeBandpassEnd - timeBandpassStart
casalog.origin('mkpipeline')
casalog.post("Calibrating antenna bandpasses took %0.2f seconds (%0.2f minutes)" % (timeBandpassTotal, timeBandpassTotal/60))
casalog.post(" ")
#=====================================================================
# Gain calibration:
casalog.origin('mkpipeline')
casalog.post("Performing gain calibration...")
timeGainStart = time.time()
if repro:
gainPh = rootName + '.' + timeStr + '.' + gainPhInt + '.gain.phase'
gainPhSmooth = rootName + '.' + timeStr + '.' + gainPhInt + '.gain.phase.sm'
gainPhAmp = rootName + '.' + timeStr + '.' + gainAmpInt + '.gain.phase.amp'
else:
gainPh = rootName + '.' + gainPhInt + '.gain.phase'
gainPhSmooth = rootName + '.' + gainPhInt + '.gain.phase.sm'
gainPhAmp = rootName + '.' + gainAmpInt + '.gain.phase.amp'
# Loop through spectral windows to create phase-only solutions:
for n in spwCalib:
zVal = tauMS[n]
gainSpwStr = '%i' % n
# gainSpwStr = '%i:%s' % (n, gainChan)
## gaincal(vis=msName, caltable=gainPh, \
## intent='*BANDPASS*,*PHASE*,*AMPLI*', \
## selectdata=False, spw=gainSpwStr, interp=['nearest'], \
## solint=gainPhInt, calmode='p', refant=refAnt, \
## gaintable=[bpassTab], minsnr=10.0, minblperant=3, \
## ## append=True, gaincurve=True, opacity=zVal)
## append=True)
if fluxFields:
gaincal(vis=msName, caltable=gainPh, \
intent='*BANDPASS*,*PHASE*,*AMPLI*', \
selectdata=False, spw=gainSpwStr, interp=['nearest'], \
solint=gainPhInt, calmode='p', refant=refAnt, \
gaintable=[bpassTab], minsnr=10.0, \
## append=True, gaincurve=True, opacity=zVal)
append=True)
else:
gaincal(vis=msName, caltable=gainPh, \
intent='*BANDPASS*,*PHASE*', \
selectdata=False, spw=gainSpwStr, interp=['nearest'], \
solint=gainPhInt, calmode='p', refant=refAnt, \
gaintable=[bpassTab], minsnr=10.0, \
## append=True, gaincurve=True, opacity=zVal)
append=True)
# Smooth this solution:
smoothcal(vis=msName, tablein=gainPh, caltable=gainPhSmooth, \
smoothtype='median', smoothtime=gainPhSmo)
# Loop through spectral windows to create phase+amplitude solutions:
for n in spwCalib:
zVal = tauMS[n]
gainSpwStr = '%i' % n
# gainSpwStr = '%i:%s' % (n, gainChan)
## gaincal(vis=msName, caltable=gainPhAmp, \
## intent='*BANDPASS*,*PHASE*,*AMPLI*', \
## selectdata=False, spw=gainSpwStr, interp=['nearest','nearest'], \
## solint=gainAmpInt, calmode='ap', refant=refAnt, \
## gaintable=[bpassTab, gainPh], minsnr=10.0, \
## minblperant=3, append=True, gaincurve=True, opacity=zVal)
if fluxFields:
gaincal(vis=msName, caltable=gainPhAmp, \
intent='*BANDPASS*,*PHASE*,*AMPLI*', \
selectdata=False, spw=gainSpwStr, interp=['nearest','nearest'], \
solint=gainAmpInt, calmode='ap', refant=refAnt, \
gaintable=[bpassTab, gainPhSmooth], minsnr=10.0, \
append=True, gaincurve=useGC, opacity=zVal)
else:
gaincal(vis=msName, caltable=gainPhAmp, \
intent='*BANDPASS*,*PHASE*', \
selectdata=False, spw=gainSpwStr, interp=['nearest','nearest'], \
solint=gainAmpInt, calmode='ap', refant=refAnt, \
gaintable=[bpassTab, gainPhSmooth], minsnr=10.0, \
append=True, gaincurve=useGC, opacity=zVal)
timeGainEnd = time.time()
timeGainTotal = timeGainEnd - timeGainStart
casalog.origin('mkpipeline')
casalog.post("Performing gain calibration took %0.2f seconds (%0.2f minutes)" % (timeGainTotal, timeGainTotal/60))
casalog.post(" ")
#=====================================================================
# Flux scaling:
if fluxFields:
casalog.origin('mkpipeline')
casalog.post("Applying flux scale to gain solutions...")
timeFluxscaleStart = time.time()
fluxStr = ''
for s in range(0, len(fluxFields)):
fluxStr = str(fluxFields[s]) + ',' + fluxStr
fluxStr = str.rstrip(fluxStr,',') # remove trailing comma
if repro:
fluxTab = rootName + '.' + timeStr + '.' + gainAmpInt + '.fscale.phase.amp'
else:
fluxTab = rootName + '.' + gainAmpInt + '.fscale.phase.amp'
cb.open(msName)
# Though we don't currently do anything with it, scaledFlux
# is a list containing one entry per source ID / SPW combination
# of the returned flux values.
scaledFlux = cb.fluxscale(tablein=gainPhAmp, tableout=fluxTab,
reference=fluxStr);
cb.close();
timeFluxscaleEnd = time.time()
timeFluxscaleTotal = timeFluxscaleEnd - timeFluxscaleStart
casalog.origin('mkpipeline')
casalog.post("Applying flux scale took %0.2f seconds (%0.2f minutes)" % (timeFluxscaleTotal, timeFluxscaleTotal/60))
casalog.post(" ")
#=====================================================================
# Plot solutions and make web pages
if doPlot:
casalog.origin('mkpipeline')
casalog.post("Making calibration plots...")
timeCalplotStart = time.time()
# Phase solutions
html = open(plotDir+'/bpphase.html', 'w')
htmlHead = '''
Initial phase solutions vs. time for '''+msName+'''
'''
html.write(htmlHead)
for j in range(0, len(antArr)):
antStr = '%s' % antArr[j]
htmlNamePh = antStr + '.time.phase.png'
plotNamePh = plotDir + '/' + htmlNamePh
htmlStr = '''
%s:
|
''' % (antStr, htmlNamePh, htmlNamePh)
html.write(htmlStr)
# Make phase vs. time plots:
plotcal(caltable=bpassPh, xaxis='time', yaxis='phase', \
antenna=antStr, \
showgui=False, figfile=plotNamePh, \
plotrange=[-1,-1,-180,180])
htmlTail = '''
'''
html.write(htmlTail)
html.close()
# Bandpass solutions, one per baseband per antenna:
for i in range(0, len(baseBand)):
html = open(plotDir+'/bandpass.spw%s.html' % (baseBand[i]), 'w')
htmlHead = '''
Bandpass phase/amplitude vs. freq for '''+msName+''', band '''+baseBand[i]+'''
'''
html.write(htmlHead)
for j in range(0, len(antArr)):
antStr = '%s' % antArr[j]
htmlNameAmp = antStr + '.spw' + baseBand[i] + '.freq.amp.png'
htmlNamePh = antStr + '.spw' + baseBand[i] + '.freq.phase.png'
plotNameAmp = plotDir + '/' + htmlNameAmp
plotNamePh = plotDir + '/' + htmlNamePh
htmlStr = '''
%s:
|
|
''' % (antStr, htmlNameAmp, htmlNameAmp, htmlNamePh, htmlNamePh)
html.write(htmlStr)
# Make amplitude vs. frequency plots:
plotcal(caltable=bpassTab, xaxis='freq', yaxis='amp', \
antenna=antStr, spw=baseBand[i], \
showgui=False, figfile=plotNameAmp)
# Make phase vs. frequency plots:
plotcal(caltable=bpassTab, xaxis='freq', yaxis='phase', \
antenna=antStr, spw=baseBand[i], \
showgui=False, figfile=plotNamePh, \
plotrange=[-1,-1,-180,180])
htmlTail = '''
'''
html.write(htmlTail)
html.close()
# Make plots of gain solutions, one per baseband per antenna:
for i in range(0, len(baseBand)):
html = open(plotDir+'/gain.spw%s.html' % (baseBand[i]), 'w')
htmlHead = '''
Gain phase/amplitude vs. freq for '''+msName+''', band '''+baseBand[i]+'''
'''
html.write(htmlHead)
for j in range(0, len(antArr)):
antStr = '%s' % antArr[j]
htmlNameAmp = antStr + '.spw' + baseBand[i] + '.time.amp.png'
htmlNamePh = antStr + '.spw' + baseBand[i] + '.time.phase.png'
plotNameAmp = plotDir + '/' + htmlNameAmp
plotNamePh = plotDir + '/' + htmlNamePh
htmlStr = '''
%s:
|
|
''' % (antStr, htmlNameAmp, htmlNameAmp, htmlNamePh, htmlNamePh)
html.write(htmlStr)
# Make amplitude vs. frequency plots:
plotcal(caltable=gainPhAmp, xaxis='time', yaxis='amp', \
antenna=antStr, spw=baseBand[i], \
showgui=False, figfile=plotNameAmp)
# Make phase vs. frequency plots:
plotcal(caltable=gainPhSmooth, xaxis='time', yaxis='phase', \
antenna=antStr, spw=baseBand[i], \
showgui=False, figfile=plotNamePh, \
plotrange=[-1,-1,-180,180])
htmlTail = '''
'''
html.write(htmlTail)
html.close()
# Make flux-scaled amplitude vs. frequency plots, if flux-scaled calibration exists:
if fluxFields:
for i in range(0, len(baseBand)):
html = open(plotDir+'/fscale.gain.spw%s.html' % (baseBand[i]), 'w')
htmlHead = '''
Gain phase/flux-scaled amplitude vs. freq for '''+msName+''', band '''+baseBand[i]+'''
'''
html.write(htmlHead)
for j in range(0, len(antArr)):
antStr = '%s' % antArr[j]
htmlNameAmp = antStr + '.spw' + baseBand[i] + '.time.amp.fscale.png'
htmlNamePh = antStr + '.spw' + baseBand[i] + '.time.phase.png'
plotNameAmp = plotDir + '/' + htmlNameAmp
htmlStr = '''
%s:
|
|
''' % (antStr, htmlNameAmp, htmlNameAmp, htmlNamePh, htmlNamePh)
html.write(htmlStr)
plotcal(caltable=fluxTab, xaxis='time', yaxis='amp', \
antenna=antStr, spw=baseBand[i], \
showgui=False, figfile=plotNameAmp)
htmlTail = '''
'''
html.write(htmlTail)
html.close()
timeCalplotEnd = time.time()
timeCalplotTotal = timeCalplotEnd - timeCalplotStart
casalog.origin('mkpipeline')
casalog.post("Calibration plotting took %0.2f seconds (%0.2f minutes)" % (timeCalplotTotal, timeCalplotTotal/60))
casalog.post(" ")
if not fluxFields:
fluxTab = gainPhAmp
else:
fluxTab = fluxTab
calTables = [bpassTab, gainPhSmooth, fluxTab]
return calTables
########################################################################
# Apply calibration
#
# msName = name of file to be calibrated
# clearPrev = run clearcal to clear prior calibration?
# calChan = channels to calibrate, e.g. '4~58'
# plotDir = output/plotting directory
# rootName = root name for output filenames
# calInterp = list of interpolation methods to use
# plotField = list of field(s) to plot, e.g. [1,2]
# timeBin = time binning to use for plots, e.g. '5'
# chanBin = channel binning to use for plots, e.g. '5'
# doCal = apply calibration?
# doPlot = make plots?
def applyCalib(msName, clearPrev, plotDir, rootName,
calTables, plotField, timeBin, chanBin, timeStr,
doCal, doPlot, useGC, spwCalib, targetField,
fluxFields):
#=====================================================================
# Clear prior calibration, if requested:
if doCal:
casalog.origin('mkpipeline')
casalog.post("Applying calibration...")
timeApplycalStart = time.time()
if clearPrev == True:
clearcal(vis=msName)
#=====================================================================
# Determine needed opacity corrections, per spw:
tauMS = plotWX(vis=msName, seasonal_weight=0.5, doPlot=False, plotDir='')
# Save a version of the flags before applying calibration:
flagStr = 'preCalib_'+timeStr
flagComm = 'Flag save before calibration on '+timeStr
flagmanager(vis=msName, mode='save', versionname=flagStr, \
comment=flagComm)
tb.open(msName+'/SPECTRAL_WINDOW')
freqList = tb.getcol('REF_FREQUENCY')/1e9
tb.close()
# Apply calibration, per spw:
for n in spwCalib:
# spwStr = str(n)+':'+calChan
spwStr = str(n)
zVal = tauMS[n]
calTableLen = len(calTables)
interpArr = ['nearest']*calTableLen
if fluxFields:
applycal(vis=msName, spw=spwStr, selectdata=False, \
intent='*BANDPASS*,*PHASE*,*AMPLI*', \
interp=interpArr, gaintable=calTables, gaincurve=useGC, \
opacity=zVal, calwt=False, flagbackup=False)
else:
applycal(vis=msName, spw=spwStr, selectdata=False, \
intent='*BANDPASS*,*PHASE*', \
interp=interpArr, gaintable=calTables, gaincurve=useGC, \
opacity=zVal, calwt=False, flagbackup=False)
if targetField:
applycal(vis=msName, spw=spwStr, selectdata=False, \
intent='*TARGET*', \
gaintable=calTables, gaincurve=useGC, \
opacity=zVal, calwt=False, flagbackup=False)
timeApplycalEnd = time.time()
timeApplycalTotal = timeApplycalEnd - timeApplycalStart
casalog.origin('mkpipeline')
casalog.post("Applying calibration took %0.2f seconds (%0.2f minutes)" % (timeApplycalTotal, timeApplycalTotal/60))
casalog.post(" ")
#=====================================================================
# Make plots of calibrated data, one per spw of each phase vs. amp and
# amp vs. baseline:
if doPlot:
casalog.origin('mkpipeline')
casalog.post("Plotting calibrated data...")
timePlotCaldataStart = time.time()
# Loop over fields and spws:
for i in range(0, len(plotField)):
plotStr = str(plotField[i])
html = open(plotDir+'/field_%s.html' % plotStr, 'w')
htmlHead = '''
Calibrated phase vs. amp and amp vs. baseline for '''+msName+''', field '''+plotStr+'''
| Phase vs. Amplitude: |
Amplitude vs. Baseline: |
'''
html.write(htmlHead)
for s in spwCalib:
# spwStr = str(s)+':'+calChan
spwStr = str(s)
plotMade = False
htmlAmpPh = 'field_'+plotStr+'.spw'+spwStr+'.amp.phase.png'
htmlAmpBl = 'field_'+plotStr+'.spw'+spwStr+'.amp.bl.png'
plotAmpPh = plotDir+'/'+htmlAmpPh
plotAmpBl = plotDir+'/'+htmlAmpBl
htmlStr = '''
spw%s:
|
|
''' % (spwStr, htmlAmpPh, htmlAmpPh, htmlAmpBl, htmlAmpBl)
html.write(htmlStr)
# Make phase vs. amplitude plots:
plotms(vis=msName, xaxis='phase', yaxis='amp', \
selectdata=True, \
field=plotStr, spw=spwStr, avgtime=timeBin, \
avgchannel=chanBin, \
plotfile=plotAmpPh, correlation='RR,LL', \
xdatacolumn='corrected', ydatacolumn='corrected', \
coloraxis='antenna1')
while (plotMade == False):
time.sleep(1)
print "Testing for plot %s..." % plotAmpPh
plotMade = os.path.isfile(plotAmpPh)
print "Plot not found."
plotMade = False
# Make amplitude vs. baseline plots:
plotms(vis=msName, xaxis='baseline', \
yaxis='amp', selectdata=True, \
field=plotStr, spw=spwStr, avgtime=timeBin, \
avgchannel=chanBin, \
plotfile=plotAmpBl, correlation='RR,LL', \
ydatacolumn='corrected', \
coloraxis='antenna1')
while (plotMade == False):
time.sleep(1)
print "Testing for plot %s..." % plotAmpBl
plotMade = os.path.isfile(plotAmpBl)
print "Plot not found."
htmlTail = '''
'''
html.write(htmlTail)
html.close()
timePlotCaldataEnd = time.time()
timePlotCaldataTotal = timePlotCaldataEnd - timePlotCaldataStart
casalog.origin('mkpipeline')
casalog.post("Plotting calibrated data took %0.2f seconds (%0.2f minutes)" % (timePlotCaldataTotal, timePlotCaldataTotal/60))
casalog.post(" ")
########################################################################
# Imaging and measurement of dynamic range -- one image per source name
# and spw for which there are valid calibration or target data.
def mkImages(msName, imageFields, rootName, plotDir, baseBand,
repro, timeStr):
casalog.origin('mkpipeline')
casalog.post("Imaging calibrated data...")
timeImgStart = time.time()
# Speed of light, m/s
C = 2.998e8
tb.open(msName+'/ANTENNA')
antPos = tb.getcol('POSITION')
tb.close()
tb.open(msName+'/SPECTRAL_WINDOW')
freqList = tb.getcol('REF_FREQUENCY')/1e9
tb.close()
tb.open(msName+'/FIELD')
srcNamesArr = tb.getcol('NAME')
tb.close()
antDist = []
for antenna1 in range(0,antPos.shape[1]):
xRef = antPos[0][antenna1]
yRef = antPos[1][antenna1]
zRef = antPos[2][antenna1]
for antenna2 in range(antenna1+1,antPos.shape[1]):
x = antPos[0][antenna2]
y = antPos[1][antenna2]
z = antPos[2][antenna2]
antDist.append(np.sqrt((xRef-x)**2.0+(yRef-y)**2.0+(zRef-z)**2.0))
# Maximum distance between antennas (m)
maxDist = max(antDist)
# Set up arrays to be used for image statistics
allFlux = np.zeros((len(imageFields),len(baseBand)),float)
## allFluxErr = np.zeros((len(imageFields),len(baseBand)),float)
allMedianRMS = np.zeros((len(imageFields),len(baseBand)),float)
allFreq = np.zeros((len(imageFields),len(baseBand)),float)
for j in range(0,len(imageFields)):
for s in range(0, len(baseBand)):
imgField = imageFields[j]
imgSPWstr = baseBand[s]
srcName = str(srcNamesArr[imgField])
# chanMax = chanNum[imgSPW] - 4
# srcSPWstr = '%s:3~%d' % (imgSPWstr,chanMax)
# Use the center SPW's frequency:
testBand = baseBand[s].find('~')
if (testBand > -1):
(bLow, bHigh) = str.split(baseBand[s], '~')
bLow = int(bLow)
bHigh = int(bHigh)
cntrSPW = int(bLow + (bHigh - bLow + 1)/2)
else:
cntrSPW = int(baseBand[s])
srcFreq = freqList[cntrSPW]
try:
# Get source name, spw frequency:
if repro:
imName = rootName+'.'+timeStr+'.field'+str(imageFields[j])+'.spw'+imgSPWstr
else:
imName = rootName+'.field'+str(imageFields[j])+'.spw'+imgSPWstr
# Determine the cell and image size (note using ref.
# frequency that is the lower bound of the spw).
# Want image to be 6xprimary beam and pixels to be
# 1/5 of the synthesized beam.
wavelength = C/(srcFreq*1e9)
maxDistL = maxDist/wavelength
# Synthesized beam in arcseconds
synthBeam=(1.2/maxDistL)*(180/np.pi)*60*60
# Pixel size (arcsec)
cellSize = synthBeam/5
# Primary beam in arcseconds
primaryBeam = 60.0*45.0/srcFreq
# Image the larger of one full primary beam or 120 synthesized beams.
# Note that the number of synthesized beams this will include depends on the
# configuration:
# A-config -> theta(PB) ~ 1200 * theta(SB)
# D-config -> theta(PB) ~ 40 * theta(SB) (so we will image 3 primary beams)
imSize = np.max([np.int(primaryBeam / cellSize), 640])
# Don't use a mask for cleaning for now
# imMask = [ctrPix-(imSize/8),ctrPix-(imSize/8), \
# ctrPix+(imSize/8),ctrPix+(imSize/8)]
casalog.origin('mkpipeline')
casalog.post("Calculated a synthesized beam size of %.3f arcsec and a primary beam size of %.2f arcmin" % (synthBeam, primaryBeam/60))
casalog.post("Imaging %s in spw%s (%s GHz) using a cell size of %.5f arcsec and" % (srcName,imgSPWstr,str(srcFreq),cellSize))
casalog.post("an image size of %i pixels, corresponding to %.1f primary beam(s)." % (imSize, (imSize / primaryBeam / cellSize)))
clean(vis=msName, imagename=imName, field=str(imgField), \
spw=imgSPWstr, imsize=[imSize], imagermode='csclean', \
cell=[str(cellSize)+'arcsec'], calready=True, \
niter=100, psfmode='clark', pbcor=True, minpb=0.05)
# Print out the images for summary page
outImg = plotDir + '/' + imName + '.jpg'
viewer(infile=imName+'.image', zoom=4, gui=False, outfile=outImg)
# Sample RMS at eight points around phase center, each of which is
# 1/4 PB from center of image and 1/10 PB on a side
ctrPix = imSize/2
rmsPix = imSize/4
lenPix = imSize/20
# Offset in X, Y:
eightPtX = [-rmsPix, 0, rmsPix, rmsPix, rmsPix, 0, -rmsPix, -rmsPix]
eightPtY = [rmsPix, rmsPix, rmsPix, 0, -rmsPix, -rmsPix, -rmsPix, 0]
rmsList = []
for i in range(0,8):
Xmin = ctrPix + eightPtX[i] - lenPix
Xmax = ctrPix + eightPtX[i] + lenPix
Ymin = ctrPix + eightPtY[i] - lenPix
Ymax = ctrPix + eightPtY[i] + lenPix
boxStr = str(Xmin)+','+str(Ymin)+','+str(Xmax)+','+str(Ymax)
statArr = imstat(imagename=imName+'.image', \
box=boxStr)
rmsList.append(statArr['rms'][0])
allMedianRMS[j,s] = np.median(rmsList)
# Find maximum in image, as well as frequency:
imMax = imstat(imagename=imName+'.image')['max'][0]
srcFreq = imhead(imName+'.image', mode='get', hdkey='crval4')['value'] / 1e9
# Image statistics: peak in image as well as noise
# statistics from several different locations. Use
# the median value of the latter to determine the
# dynamic range of the image. Use a box that is 12 x
# synthesized beam on a side:
## boxStr = str(ctrPix - 30) + ',' + str(ctrPix - 30) + ',' + \
## str(ctrPix + 30) + ',' + str(ctrPix + 30)
## casalog.origin('mkpipeline')
## casalog.post("Fitting gaussian in region %s" % boxStr)
## gaussFit = imfit(imagename=imName+'.image', \
## box=boxStr)
## srcFlux = gaussFit['results']['component0']['flux']['value'][0]
## srcFluxErr = gaussFit['results']['component0']['flux']['error'][0]
## srcFreq = gaussFit['results']['component0']['spectrum']['frequency']['m0']['value']
# Estimate dynamic range:
casalog.origin('mkpipeline')
casalog.post("Evaluating image "+imName+".image:")
## casalog.post(" source flux is "+str(srcFlux)+"("+str(srcFluxErr)+") Jy")
casalog.post(" maximum value is "+str(imMax)+" Jy / beam")
casalog.post(" median RMS is "+str(np.median(rmsList))+" Jy / beam:")
## casalog.post(" measured dynamic range is "+str(srcFlux/np.median(rmsList)))
casalog.post(" estimated dynamic range is "+str(imMax/np.median(rmsList)))
casalog.post(" ")
## allFlux[j,s] = srcFlux
## allFluxErr[j,s] = srcFluxErr
allFlux[j,s] = imMax
allFreq[j,s] = srcFreq
except Exception, instance:
casalog.origin('mkpipeline')
casalog.post("*** Error *** "+str(instance),priority='ERROR')
timeImgEnd = time.time()
timeImgTotal = timeImgEnd - timeImgStart
casalog.origin('mkpipeline')
casalog.post("Imaging calibrated data took %0.2f seconds (%0.2f minutes)" % (timeImgTotal, timeImgTotal/60))
casalog.post(" ")
# Make web page with plots
linkHTML = open(plotDir+'/imagePlots.html', 'w')
htmlText = '''
Cleaned images for '''+msName+'''
Image median RMS values and dynamic ranges:
| Source (field ID): |
SPW (frequency): |
Median image RMS: |
Source flux: |
Dynamic range: |
'''
linkHTML.write(htmlText)
for j in range(0,len(imageFields)):
srcName = str(srcNamesArr[imageFields[j]])
for b in range(0,len(baseBand)):
srcSPW = baseBand[b]
srcFreq = allFreq[j,b]
if (allMedianRMS[j,b] == 0):
htmlText = '''
| %s (%i) |
spw %s, %.2f GHz |
INVALID |
INVALID |
INVALID |
''' % (srcName,imageFields[j],srcSPW,srcFreq)
else:
htmlText = '''
| %s (%i) |
spw %s, %.2f GHz |
%d uJy / beam |
%.3f mJy / beam |
%d |
''' % (srcName,imageFields[j],srcSPW,srcFreq,allMedianRMS[j,b]*1e6,allFlux[j,b]*1e3,allFlux[j,b]/allMedianRMS[j,b])
linkHTML.write(htmlText)
htmlText = '''
Images:
'''
linkHTML.write(htmlText)
for j in range(0,len(imageFields)):
srcName = str(srcNamesArr[imageFields[j]])
srcID = str(imageFields[j])
for b in range(0,len(baseBand)):
srcSPW = baseBand[b]
if repro:
imName = rootName + '.' + timeStr + '.field' + srcID + '.spw' + srcSPW + '.jpg'
else:
imName = rootName + '.field' + srcID + '.spw' + srcSPW + '.jpg'
htmlText = ' | spw%s, %s (field ID %s): |
|---|
|
\n' % (srcSPW, srcName, srcID, imName)
linkHTML.write(htmlText)
htmlText = '''
'''
linkHTML.write(htmlText)
linkHTML.close()
########################################################################
# Make a summary web page:
def mkWebPage(msName, plotDir, plotScan, logName, fluxFields, doPlot,
repro, checkZeros, baseBand, bBandOrig, calFields):
linkHTML = open(plotDir+'/index.html', 'w')
htmlText = '''
mkpipeline: '''+msName+'''
Plots and log files for '''+msName+''':
'''
linkHTML.write(htmlText)
if not repro:
htmlText = '''
Observation metadata:
Raw data:
'''
linkHTML.write(htmlText)
htmlText = '''
Calibration:
Imaging:
Log files:
'''
linkHTML.write(htmlText)
linkHTML.close()
os.system("cat " + plotDir + "/mkpipeline.log >> " + plotDir + "/index.html")
linkHTML = open(plotDir+'/index.html', 'a')
htmlText = '''
'''
linkHTML.write(htmlText)
linkHTML.close()
########################################################################
# Check for missing scans.
def checkScans(msName):
tb.open(msName)
scanNums = sorted(np.unique(tb.getcol('SCAN_NUMBER')))
tb.close()
missingScans = 0
missingScanStr = ''
for i in range(max(scanNums)):
if scanNums.count(i+1) == 1: pass
else:
casalog.origin('mkpipeline')
casalog.post("WARNING: Scan "+str(i+1)+" is not present",priority='WARNING')
missingScans += 1
missingScanStr = missingScanStr+str(i+1)+', '
if (missingScans > 0):
casalog.origin('mkpipeline')
casalog.post("WARNING: There were "+str(missingScans)+" missing scans in this MS",priority='WARNING')
casalog.post(" ")
else:
casalog.origin('mkpipeline')
casalog.post("No missing scans found.")
casalog.post(" ")
########################################################################
# Find zeros. For now, use a clip expression that clips *outside* a
# particular range. Back up the flags first, then restore when done.
def findZeros(sdmName, timeStr, msName, plotDir):
zerosave = sdmName + '.' + timeStr + '.flagzero.backup'
zerocomment = 'widartest - saved flags before zero flagging'
fg.done()
fg.open(msName)
fg.saveflagversion(versionname=zerosave,comment=zerocomment,merge='replace')
fg.done()
flagdata(vis=msName, mode='manualflag', \
clipexpr=['ABS RR','ABS LL','ABS RL','ABS LR'], \
clipminmax=[0.0,1e-8], clipoutside=True, flagbackup=False)
zeroStats = flagdata(vis=msName, mode='summary', flagbackup=False)
# Want to make plots of zeros,
# First by spectral window:
SPWs = []
zerosSPW = []
nonZerosSPW = []
for spwStr in zeroStats['spw'].keys():
SPWs.append(int(spwStr))
zerosSPW.append(zeroStats['spw'][spwStr]['total'] - zeroStats['spw'][spwStr]['flagged'])
nonZerosSPW.append(zeroStats['spw'][spwStr]['flagged'])
spwInd = np.argsort(SPWs)
zerosSPW = np.array(zerosSPW)[spwInd]
nonZerosSPW = np.array(nonZerosSPW)[spwInd]
SPWs = np.array(SPWs)[spwInd]
# Next by baseline:
Baselines = []
zerosBL = []
nonZerosBL = []
sortedKeys = sorted(zeroStats['baseline'].keys())
for baseline in sortedKeys:
Baselines.append(baseline)
zerosBL.append(zeroStats['baseline'][baseline]['total'] - zeroStats['baseline'][baseline]['flagged'])
nonZerosBL.append(zeroStats['baseline'][baseline]['flagged'])
# And by antenna:
Antennas = []
zerosAnt = []
nonZerosAnt = []
sortedKeys = sorted(zeroStats['antenna'].keys())
for antenna in sortedKeys:
Antennas.append(antenna)
zerosAnt.append(zeroStats['antenna'][antenna]['total'] - zeroStats['antenna'][antenna]['flagged'])
nonZerosAnt.append(zeroStats['antenna'][antenna]['flagged'])
# Finally, by scan:
Scans = []
zerosScan = []
nonZerosScan = []
for scanStr in zeroStats['scan'].keys():
Scans.append(int(scanStr))
zerosScan.append(zeroStats['scan'][scanStr]['total'] - zeroStats['scan'][scanStr]['flagged'])
nonZerosScan.append(zeroStats['scan'][scanStr]['flagged'])
scanInd = np.argsort(Scans)
zerosScan = np.array(zerosScan)[scanInd]
nonZerosScan = np.array(nonZerosScan)[scanInd]
Scans = np.array(Scans)[scanInd]
# Check that there are zeros present; if not, print a message that
# asks if perhaps they were already flagged:
if ((np.sum(zerosSPW) == 0) and (np.sum(zerosBL) == 0) and (np.sum(zerosAnt) == 0) and (np.sum(zerosScan) == 0)):
casalog.origin('mkpipeline')
casalog.post("Did not find any zeros in the data. This could be because they were already flagged.")
# Plot everything:
pl.clf()
pl.subplot(221)
pl.subplots_adjust(left=0.1,right=0.97,bottom=0.1,top=0.95,wspace=0.3,hspace=0.2)
pl.title('Spectral Window')
pl.ylabel('% Zeros')
pl.xticks(pl.arange(len(SPWs)),SPWs)
zerosPctSPW = 100*np.array(zerosSPW,float)/(np.array(nonZerosSPW,float)+np.array(zerosSPW,float))
pl.plot(SPWs, zerosPctSPW, 'ro')
pl.axis([SPWs[0]-0.5, SPWs[len(SPWs)-1]+0.5, -5, 105])
pl.subplot(222)
pl.title('Scan')
pl.ylabel('% Zeros')
zerosPctScan = 100*np.array(zerosScan,float)/(np.array(nonZerosScan,float)+np.array(zerosScan,float))
pl.plot(Scans, zerosPctScan, 'r.')
pl.axis([Scans[0]-0.5, Scans[len(Scans)-1]+0.5, -5, 105])
pl.subplot(223)
pl.title('Baseline')
pl.ylabel('% Zeros')
zerosPctBL = 100*np.array(zerosBL,float)/(np.array(nonZerosBL,float)+np.array(zerosBL,float))
pl.plot(zerosPctBL, 'r.')
pl.axis([0.5, len(zerosBL)+0.5, -5, 105])
pl.subplot(224)
pl.title('Antenna')
pl.ylabel('% Zeros')
pl.xticks(pl.arange(len(Antennas)),Antennas,rotation=90)
zerosPctAnt = 100*np.array(zerosAnt,float)/(np.array(nonZerosAnt,float)+np.array(zerosAnt,float))
pl.plot(zerosPctAnt, 'r.')
pl.axis([-0.5, len(zerosAnt), -5, 105])
pl.savefig(plotDir+'/zeros.png', format='png')
# Revert to the old flags
fg.done()
fg.open(msName)
fg.restoreflagversion(versionname=zerosave, merge='replace')
fg.done()
########################################################################
# Choose a reference antenna
def chooseRefAnt(msName):
# Get antenna information
tb.open(msName+'/ANTENNA')
antStation = tb.getcol('STATION')
antName = tb.getcol('NAME')
antPos = tb.getcol('POSITION')
tb.close()
# A list of pads, in order of preference (note: may wish to tie this
# to the array configuration in the future):
bestStationList = ['E01','W02','E02','W03','E03','N02','W01','N01', \
'N03','E04','W04','N04','E05','W05','N05','E06', \
'W06','N06','E07','W07','N07','E08','W08','N08', \
'E09','W09','N09']
bestAntList = []
for i in range(0,len(bestStationList)):
if antStation.__contains__(bestStationList[i]):
name = antName[np.where(antStation == bestStationList[i])][0]
bestAntList.append(name)
# Get flagging statistics to make sure antenna isn't "overly" flagged:
fg.done()
fg.clearflagselection(0)
fg.open(msName)
fg.setdata()
fg.setflagsummary()
msFlagStats = fg.run()
fg.done()
fracFlagAnt = []
fracAntList = []
for i in range(0,len(bestAntList)):
flagged = float(msFlagStats['antenna'][bestAntList[i].upper()]['flagged'])
total = float(msFlagStats['antenna'][bestAntList[i].upper()]['total'])
fracFlagged = float(flagged/total)
if (fracFlagged < 1.0): # Antennas which aren't completely flagged
fracFlagAnt.append(fracFlagged)
fracAntList.append(bestAntList[i])
casalog.origin('mkpipeline')
casalog.post("For antenna %s there are %d flagged points, %d total points, and %.2f fraction" % (bestAntList[i],flagged,total,fracFlagged))
# Choose the first antenna with fewer flags than the mean number:
flagAntMean = np.mean(fracFlagAnt)
refAnt = ''
for i in range(0,len(fracAntList)):
if (fracFlagAnt[i] < flagAntMean):
refAnt = fracAntList[i]
break
if (refAnt == ''):
casalog.origin('mkpipeline')
casalog.post("ERROR: Unable to find a reference antenna",priority='ERROR')
return flagAntMean, refAnt, fracFlagAnt
else:
casalog.origin('mkpipeline')
casalog.post("Using reference antenna %s" % refAnt)
return refant
########################################################################
# Weather table plotting & opacity determination - J. Marvil 3.21.11
# gets and plots data from the weather table of the given MS
def plotWX(vis, seasonal_weight, doPlot, plotDir):
myMS=vis
##retrieve center frequency for each sub-band
tb.open(myMS+'/SPECTRAL_WINDOW')
spwFreqs=tb.getcol('REF_FREQUENCY') * 1e-9
tb.close()
##retrieve stuff from weather table
tb.open(myMS+'/WEATHER')
mytime=tb.getcol('TIME')
mytemp=tb.getcol('TEMPERATURE') - 273.15
mydew=tb.getcol('DEW_POINT') - 273.15
mywinds=tb.getcol('WIND_SPEED')
mywindd=tb.getcol('WIND_DIRECTION')*(180.0/pi)
mypres=tb.getcol('PRESSURE')
myhum=tb.getcol('REL_HUMIDITY')
tb.close()
##calculate the elevation of the sun
sunEL=[]
for time in mytime:
t1= qa.quantity(time,'s')
myday=qa.convert(t1,'d')
sunEL1=jm_sunEL(myday)
sunEL.append(sunEL1)
##convert time to a string of date/time
myTimestr = []
myTimestr2=[]
for time in mytime:
q1=qa.quantity(time,'s')
time1=qa.time(q1,form='ymd')
time2=qa.time(q1,form='local')
myTimestr.append(time1)
myTimestr2.append(time2)
##convert time to a decimal
numtime=pl.datestr2num(myTimestr)
##calculate opacity as in EVLA memo 143
thisday= 30*(float(myTimestr[0][5:7])-1)+float(myTimestr[0][8:10])
thisday=thisday + 5 * (thisday / 365.)
myTauZ=Tau_K_Calc(mydew,mytemp,thisday)
myTauZ1=Tau_K_Calc(mydew,mytemp,thisday, weights=(0,1))
myTauZ2=Tau_K_Calc(mydew,mytemp,thisday, weights=(1,0))
myTauZ_custom=Tau_K_Calc(mydew,mytemp,thisday,weights=(1-seasonal_weight,seasonal_weight))
myPWV = -1.71 + 1.3647*myTauZ
myPWV1 = -1.71 + 1.3647*myTauZ1
myPWV2 = -1.71 + 1.3647*myTauZ2
myPWV_custom=-1.71+1.3647*myTauZ_custom
tmp = casac.Quantity(270.0,'K')
pre = casac.Quantity(790.0,'mbar')
alt = casac.Quantity(2125,'m')
h0 = casac.Quantity(2.0,'km')
wvl = casac.Quantity(-5.6, 'K/km')
mxA = casac.Quantity(48,'km')
dpr = casac.Quantity(10.0,'mbar')
dpm = 1.2
att = 1
nb = 1
fC=casac.Quantity(25.0,'GHz')
fW=casac.Quantity(50.,'GHz')
fR=casac.Quantity(0.25,'GHz')
attool=casac.homefinder.find_home_by_name('atmosphereHome')
at=attool.create()
hum=20.0
myatm=at.initAtmProfile(alt,tmp,pre,mxA,hum,wvl,dpr,dpm,h0,att)
at.initSpectralWindow(nb,fC,fW,fR)
sg=at.getSpectralWindow()
mysg = sg.value
nstep = 20
pwv = []
opac = pl.zeros((len(mysg),nstep))
for i in range(nstep):
hum = 20.0*(i+1)
myatm = at.initAtmProfile(alt,tmp,pre,mxA,hum,wvl,dpr,dpm,h0,att)
w=at.getGroundWH2O()
pwv.append(w.value)
at.initSpectralWindow(nb,fC,fW,fR)
at.setUserWH2O(w)
sg=at.getSpectralWindow()
mysg = sg.value
sdry=at.getDryOpacitySpec()
swet=at.getWetOpacitySpec()
sd=sdry['dryOpacity']
sw=swet['wetOpacity'].value
stot = pl.array(sd)+pl.array(sw)
opac[:,i]=stot
pwv_coef=pl.zeros((len(mysg),2))
for i in range(len(mysg)):
myfit=pl.polyfit(pwv,opac[i,:],1)
pwv_coef[i,:]=myfit
freqs=pl.array(mysg)/1e9
coef0=pwv_coef[:,1]/1e-3
coef1=pwv_coef[:,0]/1e-3
#interpolate between nearest table entries for each spw center frequency
meanTau=[]
meanTau1=[]
meanTau2=[]
meanTau_custom=[]
for i in range(len(spwFreqs)):
mysearch=(pl.array(freqs)-spwFreqs[i])**2
hits=pl.find(mysearch == min(mysearch))
if len(hits) > 1: hits=hits[0]
tau_interp = (pl.array(coef0[hits-2:hits+2])+pl.array(coef1[hits-2:hits+2])*pl.mean(myPWV)) * 1e-1 #taupercent
tau_F = pl.interp(spwFreqs[i],freqs[hits-2:hits+2],tau_interp)
meanTau.append(pl.mean(tau_F))
tau_interp1 = (pl.array(coef0[hits-2:hits+2])+pl.array(coef1[hits-2:hits+2])*pl.mean(myPWV1)) * 1e-1 #taupercent
tau_F1 = pl.interp(spwFreqs[i],freqs[hits-2:hits+2],tau_interp1)
meanTau1.append(pl.mean(tau_F1))
tau_interp2 = (pl.array(coef0[hits-2:hits+2])+pl.array(coef1[hits-2:hits+2])*pl.mean(myPWV2)) * 1e-1 #taupercent
tau_F2 = pl.interp(spwFreqs[i],freqs[hits-2:hits+2],tau_interp2)
meanTau2.append(pl.mean(tau_F2))
tau_interp_c = (pl.array(coef0[hits-2:hits+2])+pl.array(coef1[hits-2:hits+2])*pl.mean(myPWV_custom)) * 1e-1 #taupercent
tau_F_c = pl.interp(spwFreqs[i],freqs[hits-2:hits+2],tau_interp_c)
meanTau_custom.append(pl.mean(tau_F_c)/100.)
tau_allF = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV)) * 1e-1
tau_allF1 = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV1)) *1e-1
tau_allF2 = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV2)) *1e-1
tau_allF_c = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV_custom)) *1e-1
##make the plot
if doPlot==False:
casalog.origin('mkpipeline')
casalog.post('SPW : Frequency (GHz) : Zenith opacity (nepers)')
for i in range(len(meanTau_custom)):
myStr = str(i).rjust(3) + ' : '
myStr2 = '%.3f'%(spwFreqs[i])
myStr += myStr2.rjust(7) + ' : ' +str(round(meanTau_custom[i], 3))
casalog.post(myStr)
casalog.post(" ")
return meanTau_custom
myColor1='#A6A6A6'
myBlueGreyColor='#92B5F2'
myColor2='#4D4DFF'
myOrangeColor='#FF6600'
myYellowColor='#FFCC00'
myWeirdColor='#006666'
myLightBrown='#996633'
myDarkGreay='#333333'
pl.rcdefaults()
pl.rc('axes',labelsize=12)
thisfig=pl.figure(1)
thisfig.clf()
thisfig.set_size_inches(8.5,10)
Xrange=numtime[-1]-numtime[0]
Yrange=max(mywinds)-min(mywinds)
Xtextoffset=-Xrange*.01
Ytextoffset=-Yrange*.08
Xplotpad=Xrange*.03
Yplotpad=Yrange*.03
sp1=thisfig.add_axes([.13,.8,.8,.15])
pl.ylabel('solar el')
nsuns=30
myj=pl.array(pl.linspace(0,len(sunEL)-1,nsuns),dtype='int')
for i in myj:
if sunEL[i]<0: pl.plot([numtime[i],numtime[i]],[(180/pi)*sunEL[i],(180/pi)*sunEL[i]],'kH')
else: pl.plot([numtime[i],numtime[i]],[(180/pi)*sunEL[i],(180/pi)*sunEL[i]],'H',color=myYellowColor)
pl.plot([numtime[0],numtime[-1]],[0,0],'-',color='brown')
xa=pl.gca(); xa.set_xticklabels('')
jm_set_Ylim_ticks(myMin=-90,myMax=90)
jm_set_Ylabel_pos(pos=(0,.5))
pl.title('Weather Summary for '+myMS)
pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad)
xa.set_xticks(pl.linspace(min(numtime),max(numtime),3))
sp2=thisfig.add_axes([.13,.65,.8,.15])
pl.ylabel('wind (m/s)')
nwind=60
myj=pl.array(pl.linspace(0,len(mywinds)-1,nwind),dtype='int')
for i in myj:
pl.text(numtime[i]+Xtextoffset,Ytextoffset+mywinds[i],'-->',rotation=mywindd[i], alpha=1,color='purple',fontsize=12)
pl.plot(numtime, .3+mywinds,'.', color='black', ms=2, alpha=0)
jm_set_Ylabel_pos(pos=(0,.5))
jm_set_Yvar_ticks(5)
xa=pl.gca(); xa.set_xticklabels('')
pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad)
pl.ylim(min(mywinds)-Yplotpad,max(mywinds)+Yplotpad)
xa.set_xticks(pl.linspace(min(numtime),max(numtime),3))
sp4=thisfig.add_axes([.13,.5,.8,.15])
pl.plot(numtime, mytemp,'-', color=myOrangeColor,lw=2)
pl.plot(numtime, mydew,'-', color=myWeirdColor,lw=2)
pl.ylabel('temp,dew')
jm_set_Ylabel_pos(pos=(0, .5))
xa=pl.gca(); xa.set_xticklabels('')
jm_set_Yvar_ticks(5)
pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad)
xa.set_xticks(pl.linspace(min(numtime),max(numtime),3))
sp7=thisfig.add_axes([.13,.35,.8,.15])
pl.ylabel('PWV (mm)')
pl.plot(numtime, myPWV, color=myBlueGreyColor,lw=2, label='average')
pl.plot(numtime, myPWV1, color=myColor1, lw=2, label='weather station')
pl.plot(numtime, myPWV2, color=myColor2, lw=2, label='seasonal')
thismin=min([min(myPWV),min(myPWV1),min(myPWV2)])
thismax=max([max(myPWV),max(myPWV1),max(myPWV2)])
pl.ylim(.8*thismin,1.2*thismax)
jm_set_Ylabel_pos(pos=(0,.5))
jm_set_Yvar_ticks(5)
xa=pl.gca(); xa.set_xticklabels('')
pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad)
middletimei=int(floor(len(myTimestr)/2.))
middletimes=str(myTimestr[middletimei])[11:]
endtimes=myTimestr[-1][11:]
ax=pl.gca()
axt=ax.get_xticks()
ax.set_xticks(pl.linspace(min(numtime),max(numtime),3))
ax.set_xticklabels([myTimestr[0],middletimes,endtimes ])
sp8=thisfig.add_axes([.13,.1,.8,.2])
pl.plot(freqs,.01*tau_allF2,'-', color=myColor2, lw=2, label='weather station')
pl.plot(freqs,.01*tau_allF1,'-', color=myColor1, lw=2, label='seasonal model')
pl.plot(freqs,.01*tau_allF,'-', color=myBlueGreyColor, lw=2,label='average')
sp8.legend(loc=2, borderaxespad=0)
pl.ylabel('Tau_Z (nepers)')
pl.xlabel('Frequency (GHz)')
pl.ylim(0,.25)
jm_set_Yvar_ticks(6)
jm_set_Ylabel_pos(pos=(0,.5))
pl.savefig(plotDir+'/weatherInfo.png',dpi=150)
pl.close()
casalog.post('wrote weather figure: weatherInfo.png')
casalog.post('SPW : Frequency (GHz) : Zenith opacity (nepers)')
for i in range(len(meanTau_custom)):
myStr = str(i).rjust(3) + ' : '
myStr2 = '%.3f'%(spwFreqs[i])
myStr += myStr2.rjust(7) + ' : ' +str(round(meanTau_custom[i], 3))
casalog.post(myStr)
casalog.post(" ")
return meanTau_custom
###############
## hides the extreme Y-axis ticks, helps stack plots close together without labels overlaping
def jm_clip_Yticks():
xa=pl.gca()
nlabels=0
for label in xa.yaxis.get_ticklabels():
nlabels+=1
thislabel=0
if nlabels>3:
for label in xa.yaxis.get_ticklabels():
if thislabel==0: label.set_alpha(0)
if thislabel==nlabels-1: label.set_alpha(0)
thislabel+=1
##############
## sets the position of the y-axis label to the right side of the plot, can also move up/down
def jm_set_Ylabel_pos(pos=(0.5,0.5)):
ax=pl.gca();
ax.yaxis.set_label_position('right')
ax.yaxis.label.set_rotation(270)
ax.yaxis.label.set_position(pos)
###############
## fixed y-ticks for plotting phase, in this case, wind direction
def jm_set_Yphase_ticks(myScale=180):
myYlocs=pl.linspace(-myScale,myScale,5)
myLocator = pl.FixedLocator(myYlocs)
ax=pl.gca()
ax.yaxis.set_major_locator( myLocator )
pl.ylim([-1.1*myScale,1.1*myScale])
jm_clip_Yticks()
###############
## fixed y-ticks, from zero to myScale, used for data that can't be negative
def jm_set_Ygain_ticks(myScale=.2):
myYlocs=pl.linspace(0,myScale,5)
myLocator = pl.FixedLocator(myYlocs)
ax=pl.gca()
ax.yaxis.set_major_locator( myLocator )
pl.ylim([-.02*myScale,1.1*myScale])
jm_clip_Yticks()
###############
## fixed y-ticks, from myMin to myMax
def jm_set_Ylim_ticks(myMin=-1,myMax=1):
myYlocs=pl.linspace(round(myMin,1),round(myMax,1),5)
myLocator = pl.FixedLocator(myYlocs)
ax=pl.gca()
ax.yaxis.set_major_locator( myLocator )
pl.ylim(myMin,myMax)
jm_clip_Yticks()
###############
## variable y-ticks, but not more than 1+ this argument
def jm_set_Yvar_ticks(myScale=4):
xa=pl.gca()
xa.yaxis.set_major_locator(pl.MaxNLocator(myScale))
jm_clip_Yticks()
###############
## calculates K-band zenith opacity from temperature and dewpoint
def Tau_K_Calc(D,T,day, weights=(.5,.5)):
P = pl.exp(1.81+(17.27*D)/(D+237.3)) # water vapor partial pressure
h = 324.7*P/(T+273.15) # PWV in mm
tau_w = 3.8 + 0.23*h + 0.065*h**2 # tau from weather, in %, at 22GHz
if day > 199: day = day - 365.
m = day + 165. # modified day of the year
tau_d = 22.1 - 0.178*m + 0.00044*m**2 # tau from seaonal model, in %
tau_k = weights[0]*tau_w + weights[1]*tau_d # the average, with equal weights (as in the AIPS default)
return tau_k
################
## calculates elevation of the sun
def jm_sunEL(mytime):
me.doframe(me.observatory('VLA'))
me.doframe(me.epoch('utc',mytime))
mysun=me.measure(me.direction('SUN'),'AZELGEO')
return mysun['m1']['value']
