Skip to content
Snippets Groups Projects
Commit 65b2e292 authored by Ger van Diepen's avatar Ger van Diepen
Browse files

bug 1634: Make cone searching possible

parent 9385fa33
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CEP/pyparmdb/src/gsm.py
+ 127
− 74
View file @ 65b2e292
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
# cat_id = 3 => NVSS # cat_id = 3 => NVSS
# cat_id = 4 => VLSS # cat_id = 4 => VLSS
# cat_id = 5 => WENSS # 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 os
import monetdb.sql as db import monetdb.sql as db
...@@ -21,30 +21,24 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids): ...@@ -21,30 +21,24 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
db_passwd = "msss" db_passwd = "msss"
#db_lang = "sql" #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 a single catalog, use =
if len(cat_ids) == 1: 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: else:
# Multiple catalogs. # Multiple catalogs.
# If a single minimum flux, use it for all catalogs. # If a single minimum flux, use it for all catalogs.
# Note that the IN clause needs (), not []. # Note that the IN clause needs (), not [].
if len(fluxi_mins) == 1: 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: else:
# We have a minimum flux per catalog. # There is a minimum flux per catalog.
cat_where = '(' catSelect = ''
for i in range(len(cat_ids)): for i in range(len(cat_ids)):
if i > 0: if i > 0:
cat_where += ' OR ' catSelect += ' OR '
cat_where += '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[i], fluxi_mins[i]) catSelect += '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[i], fluxi_mins[i])
cat_where += ')'
# Form the entire where clause. # 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 \ conn = db.connect(hostname=db_host \
,port=db_port \ ,port=db_port \
...@@ -83,36 +77,36 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids): ...@@ -83,36 +77,36 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
" END" + \ " END" + \
")))))))))))))))))))))))" + \ ")))))))))))))))))))))))" + \
" ) AS VARCHAR(300)) AS line " + \ " ) 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 " + \ " ,CASE WHEN c1.pa IS NULL " + \
" THEN CAST('POINT' AS VARCHAR(20)) " + \ " THEN CAST('POINT' AS VARCHAR(23)) " + \
" ELSE CAST('GAUSSIAN' AS VARCHAR(20)) " + \ " ELSE CAST('GAUSSIAN' AS VARCHAR(23)) " + \
" END AS src_type " + \ " END AS src_type " + \
" ,CAST(c1.ra AS VARCHAR(20)) AS ra " + \ " ,CAST(c1.ra AS VARCHAR(23)) AS ra " + \
" ,CAST(c1.decl AS VARCHAR(20)) AS decl " + \ " ,CAST(c1.decl AS VARCHAR(23)) AS decl " + \
" ,CAST(c1.i_int_avg AS VARCHAR(20)) AS i " + \ " ,CAST(c1.i_int_avg AS VARCHAR(23)) AS i " + \
" ,CAST(0 AS VARCHAR(20)) AS q " + \ " ,CAST(0 AS VARCHAR(23)) AS q " + \
" ,CAST(0 AS VARCHAR(20)) AS u " + \ " ,CAST(0 AS VARCHAR(23)) AS u " + \
" ,CAST(0 AS VARCHAR(20)) AS v " + \ " ,CAST(0 AS VARCHAR(23)) AS v " + \
" ,CASE WHEN c1.pa IS NULL " + \ " ,CASE WHEN c1.pa IS NULL " + \
" THEN CAST('' AS VARCHAR(20)) " + \ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CASE WHEN c1.major IS NULL " + \ " ELSE CASE WHEN c1.major IS NULL " + \
" THEN CAST('' AS VARCHAR(20)) " + \ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CAST(c1.major AS varchar(20)) " + \ " ELSE CAST(c1.major AS varchar(23)) " + \
" END " + \ " END " + \
" END AS MajorAxis " + \ " END AS MajorAxis " + \
" ,CASE WHEN c1.pa IS NULL " + \ " ,CASE WHEN c1.pa IS NULL " + \
" THEN CAST('' AS VARCHAR(20)) " + \ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CASE WHEN c1.minor IS NULL " + \ " ELSE CASE WHEN c1.minor IS NULL " + \
" THEN CAST('' AS VARCHAR(20)) " + \ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CAST(c1.minor AS varchar(20)) " + \ " ELSE CAST(c1.minor AS varchar(23)) " + \
" END " + \ " END " + \
" END AS MinorAxis " + \ " END AS MinorAxis " + \
" ,CASE WHEN c1.pa IS NULL " + \ " ,CASE WHEN c1.pa IS NULL " + \
" THEN CAST('' AS VARCHAR(20)) " + \ " THEN CAST('' AS VARCHAR(23)) " + \
" ELSE CAST(c1.pa AS varchar(20)) " + \ " ELSE CAST(c1.pa AS varchar(23)) " + \
" END AS Orientation " + \ " 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 " + \ " ,CASE WHEN si.spindx_degree IS NULL " + \
" THEN 0 " + \ " THEN 0 " + \
" ELSE si.spindx_degree " + \ " ELSE si.spindx_degree " + \
...@@ -149,55 +143,108 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids): ...@@ -149,55 +143,108 @@ def gsmSelect (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
return len(y) 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): 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 ''
print 'Insufficient arguments given; run as:' print 'Insufficient arguments given; run as:'
print '' 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 ''
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 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] minFlux = [0.5]
if len(argv) > 5 and len(argv[5]) > 0: if len(argv) > optarg and len(argv[optarg]) > 0:
minFlux = [float(x) for x in argv[5].split(',')] minFlux = [float(x) for x in argv[optarg].split(',')]
# cat_id = 3 => NVSS # cat_id = 3 => NVSS
# cat_id = 4 => VLSS # cat_id = 4 => VLSS
# cat_id = 5 => WENSS # cat_id = 5 => WENSS
cats = ['WENSS'] cats = ['WENSS']
if len(argv) > 6 and len(argv[6]) > 0: if len(argv) > optarg+1 and len(argv[optarg+1]) > 0:
cats = argv[6].split (',') cats = argv[optarg+1].split (',')
cat_ids = [] cat_ids = []
for cat in cats: for cat in cats:
cat = cat.upper() cat = cat.upper()
...@@ -214,13 +261,19 @@ def gsmMain (name, argv): ...@@ -214,13 +261,19 @@ def gsmMain (name, argv):
print 'Nr of minFlux values (%s) must be 1 or match nr of catalogs' % argv[5] print 'Nr of minFlux values (%s) must be 1 or match nr of catalogs' % argv[5]
return False; return False;
# Do the selection and create the makesourcedb file. # 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: if nr < 0:
return False return False
if nr == 0: elif nr == 0:
print 'Warning: no matching sources found in GSM for' print 'Warning: no matching sources found in GSM for'
print ' stRA=%f endRA=%f stDEC=%f endDEC=%f' % (stRA,endRA,stDEC,endDEC) if useBox:
print ' minFluxI=%f catalogs=%s' % (minFlux, cats) 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: else:
print '%d sources written into %s' % (nr, outfile) print '%d sources written into %s' % (nr, outfile)
return True return True
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment