# 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.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].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')