Pipelines - design, implimentation, possible publishing to VO

Please add your work plan, developments, criticisms, suggestions for enhancements etc.

Participants

Cormac Purcell, Franocois Viallefond, Ed Fomalont, Renee Hlozek, John Morgan, John Lightfoot, Willem Baan, Eleni Vardoulaki

Goals

We wanted to explore the use of ParselTongue^? for creating pipelines.

Step 1: generating an image in ParselTongue^?

We then repeated this process in an interactive ParselTongue^? session. One of us recorded the steps of this interactive session as a script which will form the basis of our pipeline:

Download this script as a file

#!/usr/local/bin/ParselTongue

################################################################################
# This is a script of an interactive session we did in parseltongue. This is 
# the kind of thing that is probably best done in classical aips. However it's a
# nice way for us to learn about ParselTongue, and may form the basis for some
# more general pipeline script
################################################################################

#You can use this file from ParselTongue using 
#execfile('[pathtothisfile]/makemap.py')
#In that case the following is optional
from AIPS import AIPS
from AIPSTV import AIPSTV
from AIPSTask import AIPSTask, AIPSList
from AIPSData import AIPSUVData, AIPSImage, AIPSCat
AIPS.userno = 100


################################################################################
# First load in the UV file using UVLOD
################################################################################
uvlod = AIPSTask('uvlod')
#                 ^Tasks are always lower case

uvlod.inf = '/aips/FITS/3C338Y-XT.FITS'
#To generalise this script we could try getting the file from command line.
#That's pretty easy in a scripting language like python.

uvdata = AIPSUVData ('3c338_uv', 'uvdata', 1, 1)
#Here we are connecting the python object uvdata with an aips file with 
#data, class, sequence and disk; '3c338_uv', 'uvdata', 1, 1 

if uvdata.exists() : uvdata.zap()
#Parseltongue will throw up an error if you try to overwrite an aips file which
#already exists. So we have to get rid of any data which already has these 
#attributes.

uvlod.outdata = uvdata
#in/outdata replaces file, class, disk and sequence.
#Similarly there is in2data, out3data etc.

#We can check this has worked:
print uvlod.outdata
print uvlod.outclass
print uvlod.outdisk
print uvlod.outseq

#Or even get the whole lot
print uvlod.inputs()

uvlod.go()
#check it's worked
if uvdata.exists():
   print "success!"

################################################################################
# Next get the TV up
################################################################################

tv = AIPSTV() #The constructor for the tv, which is an action, hence ()
tv.start()
#An AIPS TV window should appear

################################################################################
# Now let's make a UV plot
################################################################################
uvplot = AIPSTask('uvplt')
uvplot.indata = uvdata
uvplot.xinc = 50
#only plot every 50th point
uvplot.dotv = 1
uvplot.go()

uvplot2 = AIPSTask('uvplt')
uvplot.indata = uvdata
uvplot.xinc = 50

uvplot.bparm[1:]=[6,7,2]
#Another way to do this would be:
#uvplot.bparm[1]=6
#uvplot.bparm[1]=7
#uvplot.bparm[1]=2
#Or
#uvplot.bparm=AIPSList[6,7,2]
#the 'opposite' of AIPSList is PythonList

uvplot.go()

################################################################################
# Now let's make our image
################################################################################

imager=AIPSTask('imagr')
imager.indata = uvdata
imager.cellsize[1:] = [0.5,0.5]
imager.imsize[1:] = [512,512]
imager.niter = 5000

#we cant use our fancy outdata with imagr because imagr doesn't accept an outclass
imager.outname = '3C338_PT'
imager.outseq = 1
imager.outdisk = 1
#we can make a parseltongue object pointing at our output image afterwards (see below)

imager.go()

#lets have a look at what we've made
print AIPSCat()

#now we can get our data into a parseltongue object
image_338 = AIPSImage('3C338_PT', 'ICL001',1 ,1)
if image_338.exists():
   print "Success!! We have our image!!"

################################################################################
# Now let's make a contour plot
################################################################################
tv.clear()

kntr = AIPSTask('kntr')
kntr.indata = image_338
kntr.dotv = 1
kntr.docont = -1
kntr.dogrey = 1
kntr.dovect = -1

kntr.blc[1:] = [156,156]
kntr.trc[1:] = [356,356]
kntr.go()

print image_338.tables
cctable = image_338.table('CC',1)

#if you want _loads_ of output, uncomment below
#there are currently some issues with slicing table objects in parseltongue
#so we can't e.g. print every 10th row (there are workarounds but it's complicated)

#for row in cctable:
#   print "type control-c to get out"
#   print row

Step 2: refine the pipeline.

We then refined this script by adding the functionality described below.

• Add a log file to which all aips output is copied, plus some of our own information

• Allow the script to read any fits file, chosen by the user on the command line

• Choose the parameters for IMAGR based on the attributes of the image.

The resulting script is reproduced below, it should be able to read in any calibrated UV file

Download this file

#! /usr/local/bin/ParselTongue

################################################################################
# This is a script of an interactive session we did in ParselTongue. It reads in
# a fits file and attempts to image it using parameters which it calculates from
# the image header.
#
# 2.0 - Now with Added functionality. 
#
#Additional comments can be found in makemap.py
################################################################################

from AIPS import AIPS
from AIPSTV import AIPSTV
from AIPSTask import AIPSTask, AIPSList
from AIPSData import AIPSUVData, AIPSImage, AIPSCat
AIPS.userno = 100

#These are other python modules which we use for reading from the command line
#and extracting the filename from the path. 
import sys, os


################################################################################
# Open up a log file
################################################################################
AIPS.log = open('ParselTongue.log', 'a')
#All aips messages now go to the file.
#However we will also direct some of our own messages there too. 

################################################################################
# read the path from the command line and extract the filename
################################################################################
   path = sys.argv[1]
   print >> AIPS.log, "pipln :path= " + path 

   name = os.path.basename(path).split(".")
   name = name[0]
   print >> AIPS.log"pipln :name = " + name

################################################################################
# First load in the UV file using UVLOD
################################################################################
uvdata = AIPSUVData (name, 'UVDATA', 1, 1)
if uvdata.exists() : uvdata.zap()

uvlod = AIPSTask('uvlod')
uvlod.inf = path
uvlod.outdata = uvdata
uvlod.go()

if uvdata.exists():print >> AIPS.log, "pipln :UV data read from fits file"

#Next print some interesting information about the file we've just loaded
print >> AIPS.log, "pipln : RA  =         " + repr(uvdata.header.crval[4])
print >> AIPS.log, "pipln : Dec =         " + repr(uvdata.header.crval[5])
print >> AIPS.log, "pipln : Visibilities= " + repr(len(uvdata))

################################################################################
# Next get the TV up
################################################################################
tv = AIPSTV()

try:
   tv.clear()
except:
   tv.start()
#check if it's possible to communicate with the tv by trying to clear it. If 
#this doesn't work then make a new tv

################################################################################
# Now let's make a UV plot
################################################################################
uvplot = AIPSTask('uvplt')
uvplot.indata = uvdata
uvplot.xinc = int(len(uvdata)/2000)
#only plot approximately 1000 points
uvplot.dotv = 1
uvplot.go()

uvplot.indata = uvdata
uvplot.xinc = int(len(uvdata)/2000)
uvplot.bparm[1:]=[6,7,2]
uvplot.go()


################################################################################
# Now calculate some values to use as inputs to imager
################################################################################

#Next calculate a sensible value for the cellsize based on the maximum baseline.
#We can find this out using the task uvprm.
uvprm = AIPSTask('uvprm')
uvprm.indata = uvdata
uvprm.go()
MaxBaseline = uvdata.keywords['UVPRAMAX']
res =( 206000 /(MaxBaseline * 4))

#We can also calculate a sensible threshold using the RMS value.
flux = (4 * uvdata.keywords['UVPRMSF0'])

print >> AIPS.log, "pipln : MaxBaseline= " + repr(MaxBaseline)
print >> AIPS.log, "pipln : cellsize   = " + repr(res)
print >> AIPS.log, "pipln : 4*rms      = " + repr(flux) 

################################################################################
# Now let's make our image
################################################################################

image = AIPSImage(name, 'ICL001',1 ,1)
beam = AIPSImage(name, 'IBM001',1 ,1)
if image.exists(): image.zap()
if beam.exists(): beam.zap()

imager=AIPSTask('imagr')
imager.indata = uvdata
imager.cellsize[1:] = [res,res]
imager.imsize[1:] = [512,512]
imager.niter = 10000
#this is a maximum number of iterations, hopefully the clean will stop earlier 
#when the flux reaches our calculated value:
imager.flux =  flux

imager.outname = name
imager.outseq = 1
imager.outdisk = 1

imager.go()

image = AIPSImage(name, 'ICL001',1 ,1)
if image.exists():
   print "Success!! We have our image!!"

################################################################################
# Now let's make a contour plot
################################################################################
tv.clear()

kntr = AIPSTask('kntr')
kntr.indata = image
kntr.dotv = 1
kntr.docont = -1
kntr.dogrey = 1
kntr.dovect = -1

kntr.blc[1:] = [156,156]
kntr.trc[1:] = [356,356]
kntr.go()

print image.tables
cctable = image.table('CC',1)

AIPS.log.flush()

Future

The next obvious step is to calibrate the data as well.

Another possibility we were considering was to move away from AIPSTV. The UV data can be plotted from ParselTongue^? using matplotlib, a python package which provides similar plotting capabilities to matlab.

It would also be possible to read the images out to fits files, then open these using a python fits viewer. The images could then be saved as gifs.

Python could even be used to directly make a web page from these images.

Other Files

• aips_wrappers.py.txt: Alternative functions to call AIPS tasks

• CalibrateCORNISH.py.txt: Seed of a parseltongue AIPS pipeline

• sample.py.txt: Updated config file reading

• ParselMark.py.txt: Functions for running y2k-like tests from parseltongue