From 65b2e292649ad839471331a8cbf50003a738cec7 Mon Sep 17 00:00:00 2001
From: Ger van Diepen <diepen@astron.nl>
Date: Thu, 27 Jan 2011 14:21:10 +0000
Subject: [PATCH] bug 1634: Make cone searching possible
---
CEP/pyparmdb/src/gsm.py | 201 +++++++++++++++++++++++++---------------
1 file changed, 127 insertions(+), 74 deletions(-)
diff --git a/CEP/pyparmdb/src/gsm.py b/CEP/pyparmdb/src/gsm.py
index fc54364905c..b30365c8d34 100755
--- a/CEP/pyparmdb/src/gsm.py
+++ b/CEP/pyparmdb/src/gsm.py
@@ -7,7 +7,7 @@
# cat_id = 3 => NVSS
# cat_id = 4 => VLSS
# cat_id = 5 => WENSS
-def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
+def doSelect (outfile, dirSelect, fluxi_mins, cat_ids):
import os
import monetdb.sql as db
@@ -21,30 +21,24 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
db_passwd = "msss"
#db_lang = "sql"
- # If RA crosses 360 degrees, split in two parts.
- if ra_st <= ra_end:
- ra_where = '(c1.ra BETWEEN %f AND %f)' % (ra_st, ra_end)
- else:
- ra_where = '(c1.ra BETWEEN %f AND 360 OR c1.ra BETWEEN 0 AND %f)' % (ra_st, ra_end)
# If a single catalog, use =
if len(cat_ids) == 1:
- cat_where = '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[0], fluxi_mins[0])
+ catSelect = 'c1.cat_id = %d AND c1.i_int_avg >= %f' % (cat_ids[0], fluxi_mins[0])
else:
# Multiple catalogs.
# If a single minimum flux, use it for all catalogs.
# Note that the IN clause needs (), not [].
if len(fluxi_mins) == 1:
- cat_where = '(c1.cat_id IN (%s) AND c1.i_int_avg >= %f)' % (str(cat_ids)[1:-1], fluxi_mins[0])
+ catSelect = 'c1.cat_id IN (%s) AND c1.i_int_avg >= %f' % (str(cat_ids)[1:-1], fluxi_mins[0])
else:
- # We have a minimum flux per catalog.
- cat_where = '('
+ # There is a minimum flux per catalog.
+ catSelect = ''
for i in range(len(cat_ids)):
if i > 0:
- cat_where += ' OR '
- cat_where += '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[i], fluxi_mins[i])
- cat_where += ')'
+ catSelect += ' OR '
+ catSelect += '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[i], fluxi_mins[i])
# Form the entire where clause.
- where = 'WHERE %s AND %s AND c1.decl BETWEEN %f AND %f' % (cat_where, ra_where, dec_st, dec_end)
+ where = 'WHERE (%s) AND (%s)' % (catSelect, dirSelect)
conn = db.connect(hostname=db_host \
,port=db_port \
@@ -83,36 +77,36 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
" END" + \
")))))))))))))))))))))))" + \
" ) AS VARCHAR(300)) AS line " + \
- " FROM (SELECT CAST(TRIM(c1.catsrcname) AS VARCHAR(20)) AS catsrcname " + \
+ " FROM (SELECT CAST(TRIM(c1.catsrcname) AS VARCHAR(23)) AS catsrcname " + \
" ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('POINT' AS VARCHAR(20)) " + \
- " ELSE CAST('GAUSSIAN' AS VARCHAR(20)) " + \
+ " THEN CAST('POINT' AS VARCHAR(23)) " + \
+ " ELSE CAST('GAUSSIAN' AS VARCHAR(23)) " + \
" END AS src_type " + \
- " ,CAST(c1.ra AS VARCHAR(20)) AS ra " + \
- " ,CAST(c1.decl AS VARCHAR(20)) AS decl " + \
- " ,CAST(c1.i_int_avg AS VARCHAR(20)) AS i " + \
- " ,CAST(0 AS VARCHAR(20)) AS q " + \
- " ,CAST(0 AS VARCHAR(20)) AS u " + \
- " ,CAST(0 AS VARCHAR(20)) AS v " + \
+ " ,CAST(c1.ra AS VARCHAR(23)) AS ra " + \
+ " ,CAST(c1.decl AS VARCHAR(23)) AS decl " + \
+ " ,CAST(c1.i_int_avg AS VARCHAR(23)) AS i " + \
+ " ,CAST(0 AS VARCHAR(23)) AS q " + \
+ " ,CAST(0 AS VARCHAR(23)) AS u " + \
+ " ,CAST(0 AS VARCHAR(23)) AS v " + \
" ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(20)) " + \
+ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CASE WHEN c1.major IS NULL " + \
- " THEN CAST('' AS VARCHAR(20)) " + \
- " ELSE CAST(c1.major AS varchar(20)) " + \
+ " THEN CAST('' AS VARCHAR(23)) " + \
+ " ELSE CAST(c1.major AS varchar(23)) " + \
" END " + \
" END AS MajorAxis " + \
" ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(20)) " + \
+ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CASE WHEN c1.minor IS NULL " + \
- " THEN CAST('' AS VARCHAR(20)) " + \
- " ELSE CAST(c1.minor AS varchar(20)) " + \
+ " THEN CAST('' AS VARCHAR(23)) " + \
+ " ELSE CAST(c1.minor AS varchar(23)) " + \
" END " + \
" END AS MinorAxis " + \
" ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(20)) " + \
- " ELSE CAST(c1.pa AS varchar(20)) " + \
+ " THEN CAST('' AS VARCHAR(23)) " + \
+ " ELSE CAST(c1.pa AS varchar(23)) " + \
" END AS Orientation " + \
- " ,CAST(c1.freq_eff AS VARCHAR(20)) AS ReferenceFrequency " + \
+ " ,CAST(c1.freq_eff AS VARCHAR(23)) AS ReferenceFrequency " + \
" ,CASE WHEN si.spindx_degree IS NULL " + \
" THEN 0 " + \
" ELSE si.spindx_degree " + \
@@ -149,55 +143,108 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
return len(y)
+def gsmSelectCone (outfile, ra, dec, radius, fluxi_mins, cat_ids):
+ if ra < 0 or ra > 360:
+ print 'RA', ra, 'is invalid: <0 or >360 degrees'
+ return -1
+ if dec < -90 or dec > 90:
+ print 'DEC', dec, 'is invalid: <-90 or >90 degrees'
+ return -1
+ if radius < 0 or radius > 90:
+ print 'radius', radius, 'is invalid: <0 or >90 degrees'
+ return -1
+ # Determine DEC boundaries (for hopefully faster query).
+ # Convert degrees to radians.
+ import math
+ mindec = dec - radius
+ maxdec = dec + radius
+ d2r = math.pi / 180.
+ ra *= d2r
+ dec *= d2r
+ radius *= d2r
+ # Source (srcRA,srcDEC) is inside the cone if:
+ # sin(dec)*sin(srcDEC) + cos(dec)*cos(srcDEC)*cos(ra-srcRA) >= cos(radius)
+ select = '(c1.decl BETWEEN %16.12f AND %16.12f) AND (%16.14f*sin(c1.decl*%16.14f) + %16.14f*cos(c1.decl*%16.14f)*cos(%16.14f-c1.ra*%16.14f) > %16.14f)' % (mindec, maxdec, math.sin(dec), d2r, math.cos(dec), d2r, ra, d2r, math.cos(radius))
+ return doSelect (outfile, select, fluxi_mins, cat_ids)
+
+def gsmSelectBox (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
+ for ra in [ra_st,ra_end]:
+ if ra < 0 or ra > 360:
+ print 'RA', ra, 'is invalid: <0 or >360 degrees'
+ return -1
+ for dec in [dec_st,dec_end]:
+ if dec < -90 or dec > 90:
+ print 'DEC', dec, 'is invalid: <-90 or >90 degrees'
+ return -1
+ if dec_st > dec_end:
+ print 'invalid DEC: start', dec_st, '> end', dec_end
+ return -1
+ # Form the selection for the RA,DEC.
+ # If RA crosses 360 degrees, split in two parts.
+ if ra_st <= ra_end:
+ select = '(c1.ra BETWEEN %16.12f AND %16.12f)' % (ra_st, ra_end)
+ else:
+ select = '(c1.ra BETWEEN %16.12f AND 360 OR c1.ra BETWEEN 0 AND %16.12f)' % (ra_st, ra_end)
+ select += ' AND (c1.decl BETWEEN %16.12f AND %16.12f)' % (dec_st, dec_end)
+ return doSelect (outfile, select, fluxi_mins, cat_ids)
+
+# Interpret the arguments and do the selection.
def gsmMain (name, argv):
- if len(argv) < 5:
+ useBox = False
+ starg = 0
+ optarg = 4
+ if len(argv) > 0 and argv[0] == '-b':
+ # select using a box
+ useBox = True
+ starg = 1
+ optarg = 6
+ if len(argv) < optarg:
print ''
print 'Insufficient arguments given; run as:'
print ''
- print ' %s outfile stRA endRA stDEC endDEC [minFluxI [catalogs]]' % name
+ print ' %s outfile RA DEC radius [minFluxI [catalogs]]' % name
+ print 'to select using a cone'
+ print ''
+ print ' %s -b outfile stRA endRA stDEC endDEC [minFluxI [catalogs]]' % name
+ print 'to select using a box'
+ print ''
+ print ' outfile path-name of the output file'
+ print ' It will be overwritten if already existing'
+ print ' RA cone center Right Ascension (J2000, degrees)'
+ print ' DEC cone center Declination (J2000, degrees)'
+ print ' radius cone radius (degrees)'
+ print ' stRA box start Right Ascension (J2000, degrees)'
+ print ' endRA box end Right Ascension (J2000, degrees)'
+ print ' stRA can be > endRA indicating crossing 360 deg'
+ print ' minDEC box start Declination (J2000, degrees)'
+ print ' maxDEC box end Declination (J2000, degrees)'
+ print ' minFluxI minimum flux (integrated Stokes I flux in Jy)'
+ print ' If a single value is given, it is used for all catalogs.'
+ print ' If multiple values given (separated by commas), each applies'
+ print ' to the corresponding catalog.'
+ print ' Default is 0.5'
+ print ' catalogs names of catalogs to search (case-insensitive)'
+ print ' (NVSS, VLSS, and/or WENSS)'
+ print ' If multiple, separate by commas'
+ print ' Default is WENSS'
print ''
- print ' outfile path-name of the output file'
- print ' It will be overwritten if already existing'
- print ' stRA start Right Ascension (J2000, degrees)'
- print ' endRA end Right Ascension (J2000, degrees)'
- print ' stRA can be > endRA indicating crossing 360 deg'
- print ' minDEC start Declination (J2000, degrees)'
- print ' maxDEC end Declination (J2000, degrees)'
- print ' minFluxI minimum flux (integrated Stokes I in Jy)'
- print ' default = 0.5'
- print ' If a single value is given, it is used for all catalogs.'
- print ' If multiple values (separated by commas), each applies'
- print ' to the corresponding catalog.'
- print ' catalogs names of catalogs to search (case-insensitive)'
- print ' (NVSS, VLSS, and/or WENSS)'
- print ' If multiple, separate by commas'
- print ' Default is WENSS'
- return False
- outfile = argv[0]
- stRA = float(argv[1])
- endRA = float(argv[2])
- stDEC = float(argv[3])
- endDEC = float(argv[4])
- for ra in [stRA,endRA]:
- if ra < 0 or ra > 360:
- print 'RA', ra, 'is invalid: <0 or >360 degrees'
- return False
- for dec in [stDEC,endDEC]:
- if dec < -90 or dec > 90:
- print 'DEC', dec, 'is invalid: <-90 or >90 degrees'
- return False
- if stDEC > endDEC:
- print 'invalid DEC: start', stDEC, '> end', endDEC
return False
+ # Get the arguments.
+ outfile = argv[starg]
+ a1 = float(argv[starg+1])
+ a2 = float(argv[starg+2])
+ a3 = float(argv[starg+3])
+ if useBox:
+ a4 = float(argv[starg+4])
minFlux = [0.5]
- if len(argv) > 5 and len(argv[5]) > 0:
- minFlux = [float(x) for x in argv[5].split(',')]
+ if len(argv) > optarg and len(argv[optarg]) > 0:
+ minFlux = [float(x) for x in argv[optarg].split(',')]
# cat_id = 3 => NVSS
# cat_id = 4 => VLSS
# cat_id = 5 => WENSS
cats = ['WENSS']
- if len(argv) > 6 and len(argv[6]) > 0:
- cats = argv[6].split (',')
+ if len(argv) > optarg+1 and len(argv[optarg+1]) > 0:
+ cats = argv[optarg+1].split (',')
cat_ids = []
for cat in cats:
cat = cat.upper()
@@ -214,13 +261,19 @@ def gsmMain (name, argv):
print 'Nr of minFlux values (%s) must be 1 or match nr of catalogs' % argv[5]
return False;
# Do the selection and create the makesourcedb file.
- nr = gsmSelect (outfile, stRA, endRA, stDEC, endDEC, minFlux, cat_ids)
+ if useBox:
+ nr = gsmSelectBox (outfile, a1, a2, a3, a4, minFlux, cat_ids)
+ else:
+ nr = gsmSelectCone (outfile, a1, a2, a3, minFlux, cat_ids)
if nr < 0:
return False
- if nr == 0:
+ elif nr == 0:
print 'Warning: no matching sources found in GSM for'
- print ' stRA=%f endRA=%f stDEC=%f endDEC=%f' % (stRA,endRA,stDEC,endDEC)
- print ' minFluxI=%f catalogs=%s' % (minFlux, cats)
+ if useBox:
+ print ' stRA=%f endRA=%f stDEC=%f endDEC=%f' % (a1,a2,a3,a4)
+ else:
+ print ' RA=%f DEC=%f radius=%f' % (a1,a2,a3)
+ print ' minFluxI=%s catalogs=%s' % (str(minFlux), cats)
else:
print '%d sources written into %s' % (nr, outfile)
return True
--
GitLab