# Tested against CASA 3.3 r14258, 3/10/2011
# This task queries Vizier to find nearby sources to the input
# field and writes an outlierfile that can be fed into CASA
# imaging tools. Similar to the AIPS task SETFC.
#
# Requires that you have the 'vizquery' tool installed and in
# the path; this may be obtained at:
# http://cdsarc.u-strasbg.fr/doc/cdsclient.html
#
# Note: currently the 'flux' cutoff won't work for magnitudes,
# since in that case, one would want to select on values < the
# cutoff value
#
# v0.3 - mkrauss, 3/10/2011: fixed formatting issue w/returned values
# v0.2 - mkrauss, 10/6/2010
# v0.1 - mkrauss, 11/25/2009
import os, sys, re, time
from taskinit import *
def findsources(outlierfile = None,
querycenter = None,
vis = None,
field = None,
ra = None,
dec = None,
catalog = None,
fluxparam = None,
fluxlim = None,
minradius = None,
maxradius = None,
imsize = None,
clobber = None):
casalog.origin('findsources')
# Checks that needed parameters are set
if outlierfile=='':
casalog.post("ERROR: you must specify a filename for the outlierfile", \
priority='ERROR')
return
if catalog=='':
casalog.post("ERROR: you must specify a VizieR catalog to query", \
priority='ERROR')
return
if (fluxparam == '') and (fluxlim > 0):
casalog.post("ERROR: fluxparam is not set, but there is a flux limit.", \
priority='ERROR')
casalog.post(" Please set fluxparam or set fluxlim = 0.", \
priority='ERROR')
return
casalog.post("Input parameters:")
if vis: casalog.post("vis = "+vis)
if field: casalog.post("field = %s" % field)
if ra: casalog.post("ra = %.6f" % ra)
if dec: casalog.post("dec = %.6f" % dec)
casalog.post("outlierfile = %s" % outlierfile)
casalog.post("catalog = %s" % catalog)
casalog.post("fluxparam = %s " % fluxparam)
casalog.post("fluxlim = %.3f" % fluxlim)
casalog.post("minradius = %.2f arcmin" % minradius)
casalog.post("maxradius = %.2f arcmin" % maxradius)
casalog.post("imsize = [ %i, %i ]" % (imsize[0], imsize[1]))
casalog.post("clobber = %s" % clobber)
# Check status of output file existence; clobber if requested or exit with
# error:
outFileExists = os.path.isfile(outlierfile)
if not outFileExists:
pass
elif outFileExists & clobber:
casalog.post("Found a file called %s; deleting before proceeding..." % \
outlierfile, priority = 'WARN')
os.system('rm %s' % outlierfile)
elif not clobber:
casalog.post("A file called %s already exists and clobber = False." % \
outlierfile, priority = 'ERROR')
return
# If coordinates are requested by field:
if (querycenter == 'byfield'):
tb.open(vis + '/FIELD')
srccoords = tb.getcol('REFERENCE_DIR')
srcnames = tb.getcol('NAME')
tb.close()
# Choose the appropriate SRC_ID based on field index or name:
fieldIndex = ms.msseltoindex(vis, field = field)['field'][0]
RAInp = {'unit': 'rad', 'value': srccoords[0][0][fieldIndex]}
DecInp = {'unit': 'rad', 'value': srccoords[1][0][fieldIndex]}
casalog.post("Will use phase center of %s (field %s) for VizieR query" % \
(srcnames[fieldIndex], fieldIndex), priority = 'NORMAL')
# If coordinates are requested by RA, Dec (decimal degrees):
if (querycenter == 'bycoord'):
RAInp = {'unit': 'deg', 'value': ra}
DecInp = {'unit': 'deg', 'value': dec}
# Query VizieR online catalogs given input coordinates
RAHMS = re.sub(':', ' ', qa.formxxx(RAInp, format='hms'))
DecDMS = re.sub('\.', ' ', qa.formxxx(DecInp, format='dms'), count=2)
casalog.post("Querying VizieR catalog %s around position RA = %s, Dec = %s" % \
(catalog, RAHMS, DecDMS), priority = 'NORMAL')
# Get current time for naming query and return files from VizieR:
timeArr = time.localtime()
vizQueryFile = 'VizieR_query-%i.%i.%i-%ih%im%is.txt' % \
(timeArr[0], timeArr[1], timeArr[2], timeArr[3], \
timeArr[4], timeArr[5])
vizOutFile = 'VizieR_return-%i.%i.%i-%ih%im%is.txt' % \
(timeArr[0], timeArr[1], timeArr[2], timeArr[3], \
timeArr[4], timeArr[5])
casalog.post("Writing VizieR query file to %s and " % vizQueryFile, priority = 'NORMAL')
casalog.post(" VizieR output to %s" % vizOutFile, priority = 'NORMAL')
# Write the VizieR query file
f = open(vizQueryFile,'w')
print >>f,'-source='+catalog
print >>f,'-c='+RAHMS+' '+DecDMS
print >>f,'-c.rm='+str(maxradius)
print >>f,'-out= _1 _r _RA _DE '+fluxparam
print >>f,'-sort=_r'
f.close()
# Query VizieR and parse output
os.system('vizquery -mime=csv %s > %s ' % (vizQueryFile, vizOutFile))
f = open(vizOutFile, 'r')
lines = f.readlines()
f.close()
# Write the outlierfile
f = open(outlierfile, 'w')
# Test to see if the flux column was found; ignore sources closer than
# minradius from query center:
srcIndex = 0
for i in range(0,len(lines)):
if lines[i].lstrip().startswith('1'):
arr = lines[i].rstrip().split(';')
arrlen = len(arr)
# Case: fluxlim is > 0, but there is no flux column in the output.
if (arrlen < 5) and (fluxlim > 0):
#casalog.post("case 1")
casalog.post("The supplied flux density key (%s) does not seem to be correct for %s:" % (fluxparam, catalog), priority = 'ERROR')
casalog.post(" please check the catalog and try again.", priority = 'ERROR')
f.close()
return
# Case: fluxlim is > 0, and there are exactly 5 columns, as
# expected (the fifth contains the flux information).
elif (arrlen == 5) and (fluxlim > 0):
#casalog.post("case 2")
if (float(arr[1]) < minradius) and (float(arr[4]) > fluxlim):
casalog.post("There is a source above the flux limit at RA=%s, Dec=%s; %s arcmin from the query center." % \
(arr[2], arr[3], arr[1]), priority = 'NORMAL')
casalog.post(" This is less than the requested minimum radius: not including in the outlier file", priority = 'NORMAL')
elif (float(arr[1]) > minradius) and (float(arr[4]) > fluxlim):
print >>f,'C %s %s %s %s %s' % \
(srcIndex, imsize[0], str(imsize[1]), arr[2], arr[3])
srcIndex += 1
# Case: fluxlim = 0 and there is no flux column.
elif (arrlen < 5) and (fluxlim == 0):
#casalog.post("case 3")
if float(arr[1]) < minradius:
casalog.post("There is a source at RA=%s, Dec=%s; %s arcmin from the query center." % \
(arr[2], arr[3], arr[1]), priority = 'NORMAL')
casalog.post(" This is less than the requested minimum radius: not including in the outlier file", priority = 'NORMAL')
else:
print >>f,'C %s %s %s %s %s' % \
(srcIndex, imsize[0], str(imsize[1]), arr[2], arr[3])
srcIndex += 1
# Case: fluxlim = 0 but there is a flux column.
elif (arrlen == 5) and (fluxlim == 0):
#casalog.post("case 3")
if float(arr[1]) < minradius:
casalog.post("There is a source at RA=%s, Dec=%s; %s arcmin from the query center." % \
(arr[2], arr[3], arr[1]), priority = 'NORMAL')
casalog.post(" This is less than the requested minimum radius: not including in the outlier file", priority = 'NORMAL')
else:
print >>f,'C %s %s %s %s %s' % \
(srcIndex, imsize[0], str(imsize[1]), arr[2], arr[3])
srcIndex += 1
# Case: something is not as expected. Return with an error.
else:
#casalog.post("case 4")
casalog.post("ERROR: findsources is unable to parse the VizieR output.", priority = 'ERROR')
casalog.post(" Please check the VizieR input and output files", priority = 'ERROR')
casalog.post(" %s and %s. " % (vizQueryFile, vizOutFile), priority = 'ERROR')
f.close()
return
f.close()
# Give a warning if no sources were found for the outlierfile. Else, print
# how many sources found:
if srcIndex == 0:
casalog.post("WARNING: no sources found that satisfy the input parameters.", priority = 'WARN')
else:
casalog.post("%i sources written to outlierfile %s" % (srcIndex, outlierfile), priority = 'NORMAL')
