"""
Catalog table tools
"""
import os
import inspect
from collections import OrderedDict
import glob
import traceback
import numpy as np
import astropy.io.fits as pyfits
import astropy.units as u
import astropy.coordinates as coord
from astropy.table import Table
from . import utils
__all__ = ["table_to_radec",
"table_to_regions",
"randomize_segmentation_labels",
"get_ukidss_catalog",
"get_sdss_catalog",
"get_twomass_catalog",
"get_irsa_catalog",
"get_gaia_radec_at_time",
"get_gaia_DR2_vizier_columns",
"get_gaia_DR2_vizier",
"gaia_dr2_conesearch_query",
"get_gaia_DR2_catalog",
"gen_tap_box_query",
"query_tap_catalog",
"get_hubble_source_catalog",
"get_nsc_catalog",
"get_desdr1_catalog",
"get_skymapper_catalog",
"get_vexas_catalog",
"get_vhs_catalog",
"get_panstarrs_catalog",
"get_radec_catalog"]
[docs]def table_to_radec(table, output='coords.radec'):
"""Make a "radec" ascii file with ra, dec columns from a table object
"""
if 'X_WORLD' in table.colnames:
rc, dc = 'X_WORLD', 'Y_WORLD'
elif 'x_world' in table.colnames:
rc, dc = 'x_world', 'y_world'
elif 'ra' in table.colnames:
rc, dc = 'ra', 'dec'
else:
raise ValueError("Couldn't identify sky coordinates (x_world, ra)")
table[rc].format = '.8f'
table[dc].format = '.8f'
table[rc, dc].write(output, format='ascii.commented_header',
overwrite=True)
[docs]def table_to_regions(table, output='ds9.reg', comment=None, header='global color=green width=1', size=0.5, use_ellipse=False, use_world=True, verbose=True, unit_offset_colums=['x'], scale_major=1, extra=None):
"""Make a DS9 region file from a table object
"""
fp = open(output, 'w')
fp.write(header+'\n')
if ('X_WORLD' in table.colnames) & (use_world):
xc, yc = 'X_WORLD', 'Y_WORLD'
is_world = True
elif ('ra' in table.colnames) & (use_world):
xc, yc = 'ra', 'dec'
is_world = True
elif ('x_world' in table.colnames) & (use_world):
xc, yc = 'x_world', 'y_world'
is_world = True
elif 'x_image' in table.colnames:
xc, yc = 'x_image','y_image'
is_world = False
elif 'x' in table.colnames:
xc, yc = 'x','y'
is_world = False
else:
raise ValueError("Couldn't identify ra/dec or x/y columns")
xdata = table[xc]*1
ydata = table[yc]*1
if xc in unit_offset_colums:
if verbose:
print(f'Unit offset for columns {xc}, {yc}')
xdata += 1
ydata += 1
if is_world:
fp.write('icrs\n')
sec='"'
else:
fp.write('image\n')
sec=''
# GAIA
if 'solution_id' in table.colnames:
e = np.sqrt(table['ra_error']**2+table['dec_error']**2)/1000.
e = np.maximum(e, 0.1)
else:
e = np.ones(len(table))*size
if use_ellipse:
if is_world:
if 'a_world' in table.colnames:
amaj = table['a_world']*3600*scale_major
amin = table['b_world']*3600*scale_major
etheta = table['theta_world']/np.pi*180#+90
else:
use_ellipse = False
else:
if 'a_image' in table.colnames:
amaj = table['a_image']*scale_major
amin = table['b_image']*scale_major
etheta = table['theta_image']/np.pi*180#+90
else:
use_ellipse = False
if verbose:
print(f'{output}: x = {xc}, y={yc}, ellipse={use_ellipse}')
if use_ellipse:
regstr = 'ellipse({x:.7f}, {y:.7f}, {a:.3f}{sec}, {b:.3f}{sec}, {theta:.1f})'
lines = [regstr.format(x=xdata[i], y=ydata[i],
a=amaj[i], b=amin[i],
sec=sec, theta=etheta[i])
for i in range(len(table))]
else:
regstr = 'circle({0:.7f}, {1:.7f}, {2:.3f}{3})'
lines = [regstr.format(xdata[i], ydata[i], e[i], sec)
for i in range(len(table))]
if comment is not None:
for i in range(len(table)):
lines[i] += ' # text={{{0}}}'.format(comment[i])
if extra is not None:
for i, li in enumerate(lines):
if '#' not in li:
lines[i] += ' #'
lines[i] += ' '+extra[i]
# newline
lines = [l+'\n' for l in lines]
fp.writelines(lines)
fp.close()
[docs]def randomize_segmentation_labels(seg, random_seed=1, fill_value=np.nan):
"""
Randomize labels on a segmentation image for easier visualization
Parameters
----------
seg : str or array
Segmentation filename or integer array
random_seed : int
Random number seed for `numpy.random.seed`
fill_value : scalar
Value to insert where ``seg == 0``
Returns
-------
rand_ids : array
1D array of the randomized labels with size `max(seg)+1`
rand_seg : array
2D rray with randomized labels
"""
np.random.seed(random_seed)
if isinstance(seg, str):
seg_im = pyfits.open(seg)
seg_data = seg_im[0].data
else:
seg_data = seg
imax = seg_data.max()
rand_ids = np.append(fill_value, np.argsort(np.random.rand(imax))+1)
return rand_ids, rand_ids[seg_data]
[docs]def get_ukidss_catalog(ra=165., dec=34.8, radius=3, database='UKIDSSDR9PLUS',
programme_id='LAS'):
"""Query for objects in the UKIDSS catalogs
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
Returns
-------
table : `~astropy.table.Table`
Result of the query
"""
from astroquery.ukidss import Ukidss
coo = coord.SkyCoord(ra*u.deg, dec*u.deg)
table = Ukidss.query_region(coo, radius=radius*u.arcmin,
database=database, programme_id=programme_id)
return table
[docs]def get_sdss_catalog(ra=165.86, dec=34.829694, radius=3):
"""Query for objects in the SDSS photometric catalog
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
Returns
-------
table : `~astropy.table.Table`
Result of the query
"""
from astroquery.sdss import SDSS
coo = coord.SkyCoord(ra*u.deg, dec*u.deg)
fields = ['ra', 'dec', 'raErr', 'decErr', 'petroMag_r', 'petroMagErr_r']
# print fields
fields = None
table = SDSS.query_region(coo, radius=radius*u.arcmin, spectro=False,
photoobj_fields=fields)
return table
[docs]def get_twomass_catalog(ra=165.86, dec=34.829694, radius=3, catalog='allwise_p3as_psd'):
return get_irsa_catalog(ra=ra, dec=dec, radius=radius, catalog='fp_psc', wise=False, twomass=True)
[docs]def get_irsa_catalog(ra=165.86, dec=34.829694, tab=None, radius=3, catalog='allwise_p3as_psd', wise=False, twomass=False, ROW_LIMIT=500000, TIMEOUT=3600):
"""Query for objects in the `AllWISE <http://wise2.ipac.caltech.edu/docs/release/allwise/>`_ source catalog
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
Returns
-------
table : `~astropy.table.Table`
Result of the query
"""
from astroquery.irsa import Irsa
Irsa.ROW_LIMIT = ROW_LIMIT
Irsa.TIMEOUT = TIMEOUT
#all_wise = 'wise_allwise_p3as_psd'
#all_wise = 'allwise_p3as_psd'
if wise:
catalog = 'allwise_p3as_psd'
elif twomass:
catalog = 'fp_psc'
coo = coord.SkyCoord(ra*u.deg, dec*u.deg)
table = Irsa.query_region(coo, catalog=catalog, spatial="Cone",
radius=radius*u.arcmin, get_query_payload=False)
return table
#
[docs]def get_gaia_radec_at_time(gaia_tbl, date=2015.5, format='decimalyear'):
"""
Use `~astropy.coordinates.SkyCoord.apply_space_motion` to compute GAIA positions at a specific observation date
Parameters
----------
gaia_tbl : `~astropy.table.Table`
GAIA table query, e.g., provided by `get_gaia_DR2_catalog`.
date : e.g., float
Observation date that can be parsed with `~astropy.time.Time`
format : str
Date format, see `~astropy.time.Time.FORMATS`.
Returns
-------
coords : `~astropy.coordinates.SkyCoord`
Projected sky coordinates.
"""
from astropy.time import Time
from astropy.coordinates import SkyCoord
# Distance and radial_velocity are dummy numbers needed
# to get the space motion correct
try:
# Try to use pyia
import pyia
g = pyia.GaiaData(gaia_tbl)
coord = g.get_skycoord(distance=1*u.kpc, frame='icrs',
radial_velocity=0.*u.km/u.second)
except:
# From table itself
if 'ref_epoch' in gaia_tbl.colnames:
ref_epoch = Time(gaia_tbl['ref_epoch'].data,
format='decimalyear')
else:
ref_epoch = Time(2015.5, format='decimalyear')
coord = SkyCoord(ra=gaia_tbl['ra'], dec=gaia_tbl['dec'],
pm_ra_cosdec=gaia_tbl['pmra'],
pm_dec=gaia_tbl['pmdec'],
obstime=ref_epoch,
frame='icrs',
distance=1*u.kpc, radial_velocity=0.*u.km/u.second)
new_obstime = Time(date, format=format)
coord_at_time = coord.apply_space_motion(new_obstime=new_obstime)
return(coord_at_time)
#GAIA_DR2_COLUMNS = None
GAIA_DR2_COLUMNS = ['RA_ICRS','DE_ICRS','Epoch','e_RA_ICRS','e_DE_ICRS',
'pmRA','e_pmRA','pmDE','e_pmDE',
'NAL','NAC','Solved','APF','WAL']
[docs]def get_gaia_DR2_vizier_columns():
"""
Get translation of Vizier GAIA DR2 columns. The subset of desired
output columns can be specified in the GAIA_DR2_COLUMNS global variable,
which, if `None`, will be all available columns.
"""
file = os.path.join(os.path.dirname(__file__), 'data/gaia_dr2_vizier_columns.txt')
lines = open(file).readlines()[1:]
gdict = OrderedDict()
for line in lines:
viz_col = line.split()[0]
gaia_col = line.split()[1][1:-1]
if GAIA_DR2_COLUMNS is not None:
if viz_col in GAIA_DR2_COLUMNS:
gdict[viz_col] = gaia_col
else:
gdict[viz_col] = gaia_col
return gdict
[docs]def get_gaia_DR2_vizier(ra=165.86, dec=34.829694, radius=3., max=100000,
catalog="I/355/gaiadr3", server='vizier.u-strasbg.fr',
use_mirror=False, keys=None, mjd=None, clean_mjd=True):
"""
Query GAIA catalog from Vizier
DR2: ``catalog="I/345/gaia2"``
EDR3: ``catalog="I/350/gaiaedr3"``
DR3: ``catalog="I/355/gaiadr3"``
"""
from astroquery.vizier import Vizier
import astropy.units as u
import astropy.coordinates as coord
coo = coord.SkyCoord(ra=ra, dec=dec, unit=(u.deg, u.deg),
frame='icrs')
gdict = get_gaia_DR2_vizier_columns()
if keys is None:
keys = list(gdict.keys())
try:
# Hack, Vizier object doesn't seem to allow getting all keys
# simultaneously (astroquery v0.3.7)
N = 9
for i in range(len(keys)//N+1):
v = Vizier(catalog=catalog, columns=['+_r']+keys[i*N:(i+1)*N])
v.VIZIER_SERVER = server
v.ROW_LIMIT = max
tab = v.query_region(coo, radius="{0}m".format(radius),
catalog=catalog)[0]
if i == 0:
result = tab
else:
for k in tab.colnames:
#print(i, k)
result[k] = tab[k]
for k in gdict:
if k in result.colnames:
result.rename_column(k, gdict[k])
except:
utils.log_exception(utils.LOGFILE, traceback)
return False
if mjd is not None:
rd = get_gaia_radec_at_time(result, date=mjd, format='mjd')
result['ra_time'] = rd.ra.deg
result['dec_time'] = rd.dec.deg
result.meta['cootime'] = mjd, 'Specified MJD for ra/dec_time'
if clean_mjd:
ok = np.isfinite(rd.ra.deg + rd.dec.deg)
result = result[ok]
return result
def gaia_catalog_for_assoc(assoc_name='j191132p1652_hen-2-427-f335m_00007'):
"""
Get GAIA catalog
"""
from grizli.aws import db
from grizli import utils
res = db.SQL(f"""select t_min, t_max, filter, footprint
from assoc_table where assoc_name = '{assoc_name}'
""")
sr = utils.SRegion(res['footprint'][0])
for fp in res['footprint'][1:]:
sri = utils.SRegion(fp)
sr.xy.extend(sri.xy)
un = sr.union()
radius = np.sqrt(un.sky_area()[0]).value
gaia = get_gaia_DR2_vizier(*un.centroid[0], radius=radius)
rd = get_gaia_radec_at_time(gaia, np.mean(res['t_min']), format='mjd')
gaia['ra_orig'] = gaia['ra']
gaia['ra_dec'] = gaia['dec']
gaia['ra'] = rd.ra.deg
gaia['dec'] = rd.dec.deg
coo = np.array([gaia['ra'], gaia['dec']])
in_footprint = un.path[0].contains_points(coo.T)
gaia['in_assoc'] = in_footprint
gaia.meta['assocnam'] = assoc_name
gaia['valid'] = np.isfinite(gaia['ra'] + gaia['dec'])
return gaia
[docs]def gaia_dr2_conesearch_query(ra=165.86, dec=34.829694, radius=3., max=100000):
"""
Generate a query string for the TAP servers
TBD
Parameters
----------
ra, dec : float
RA, Dec in decimal degrees
radius : float
Search radius, in arc-minutes.
Returns
-------
query : str
Query string
"""
query = "SELECT TOP {3} * FROM gaiadr2.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',{0},{1},{2:.2f}))=1".format(ra, dec, radius/60., max)
return query
[docs]def get_gaia_DR2_catalog(ra=165.86, dec=34.829694, radius=3.,
use_mirror=True, max_wait=20,
max=100000, output_file='gaia.fits'):
"""Query GAIA DR2 astrometric catalog
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
use_mirror : bool
If True, use the mirror at `gaia.ari.uni-heidelberg.de`. Otherwise
use `gea.esac.esa.int`.
Returns
-------
table : `~astropy.table.Table`
Result of the query
"""
try:
import httplib
from urllib import urlencode
except:
# python 3
import http.client as httplib
from urllib.parse import urlencode
# import http.client in Python 3
# import urllib.parse in Python 3
import time
from xml.dom.minidom import parseString
host = "gea.esac.esa.int"
port = 80
pathinfo = "/tap-server/tap/async"
if use_mirror:
host = "gaia.ari.uni-heidelberg.de"
pathinfo = "/tap/async"
# -------------------------------------
# Create job
query = gaia_dr2_conesearch_query(ra=ra, dec=dec, radius=radius, max=max) # "SELECT TOP 100000 * FROM gaiadr2.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',{0},{1},{2:.2f}))=1".format(ra, dec, radius/60.)
print(query)
params = urlencode({
"REQUEST": "doQuery",
"LANG": "ADQL",
"FORMAT": "fits",
"PHASE": "RUN",
"QUERY": query
})
headers = {
"Content-type": "application/x-www-form-urlencoded",
"Accept": "text/plain"
}
connection = httplib.HTTPConnection(host, port)
connection.request("POST", pathinfo, params, headers)
# Status
response = connection.getresponse()
print("Status: " + str(response.status), "Reason: " + str(response.reason))
# Server job location (URL)
location = response.getheader("location")
print("Location: " + location)
# Jobid
jobid = location[location.rfind('/')+1:]
print("Job id: " + jobid)
connection.close()
# -------------------------------------
# Check job status, wait until finished
tcount = 0
while True:
connection = httplib.HTTPConnection(host, port)
connection.request("GET", pathinfo+"/"+jobid)
response = connection.getresponse()
data = response.read()
# XML response: parse it to obtain the current status
dom = parseString(data)
if use_mirror:
phaseElement = dom.getElementsByTagName('phase')[0]
else:
phaseElement = dom.getElementsByTagName('uws:phase')[0]
phaseValueElement = phaseElement.firstChild
phase = phaseValueElement.toxml()
print("Status: " + phase)
# Check finished
if phase == 'COMPLETED':
break
#wait and repeat
time.sleep(0.2)
tcount += 0.2
if (phase == 'ERROR') | (tcount > max_wait):
return False
# print "Data:"
# print data
connection.close()
# -------------------------------------
# Get results
connection = httplib.HTTPConnection(host, port)
connection.request("GET", pathinfo+"/"+jobid+"/results/result")
response = connection.getresponse()
data = response.read()
outputFileName = output_file + (not use_mirror)*".gz"
try:
outputFile = open(outputFileName, "w")
outputFile.write(data)
except:
# Python 3
outputFile = open(outputFileName, "wb")
outputFile.write(data)
outputFile.close()
connection.close()
print("Data saved in: " + outputFileName)
if not use_mirror:
# ESA archive returns gzipped
try:
os.remove(output_file)
except:
pass
os.system('gunzip {output_file}.gz'.format(output_file=output_file))
table = Table.read(output_file, format='fits')
return table
[docs]def gen_tap_box_query(ra=165.86, dec=34.829694, radius=3., corners=None, max=100000, db='ls_dr7.tractor_primary', columns=['*'], rd_colnames=['ra', 'dec'], wcs_pad=0.5):
"""
Generate a query string for the NOAO Legacy Survey TAP server
Parameters
----------
ra, dec : float
RA, Dec in decimal degrees
radius : float
Search radius, in arc-minutes.
corners : 4-tuple, `~astropy.wcs.WCS` or None
ra_min, ra_max, dec_min, dec_max of a query box to use instead of
`radius`. Or if a `~astropy.wcs.WCS` object, get limits from the
`~astropy.wcs.WCS.calc_footprint` method
Returns
-------
query : str
Query string
"""
rmi = radius/60/2
cosd = np.cos(dec/180*np.pi)
if max is not None:
maxsel = 'TOP {0}'.format(max)
else:
maxsel = ''
if corners is not None:
if hasattr(corners, 'calc_footprint'):
foot = corners.calc_footprint()
left = foot[:,0].min()
right = foot[:,0].max()
bottom = foot[:,1].min()
top = foot[:,1].max()
dx = (right-left)
dy = (top-bottom)
left -= wcs_pad*dx
right += wcs_pad*dx
bottom -= wcs_pad*dy
top += wcs_pad*dy
elif len(corners) != 4:
msg = 'corners needs 4 values (ra_min, ra_max, dec_min, dec_max)'
raise ValueError(msg)
else:
left, right, bottom, top = corners
else:
left = ra - rmi / cosd
right = ra + rmi / cosd
bottom = dec - rmi
top = dec + rmi
fmt = dict(rc=rd_colnames[0], dc=rd_colnames[1],
left=left, right=right,
top=top, bottom=bottom,
maxsel=maxsel,
db=db,
output_columns=', '.join(columns))
if not np.isfinite(ra+dec):
query = "SELECT {maxsel} {output_columns} FROM {db} "
else:
query = ("SELECT {maxsel} {output_columns} FROM {db} WHERE " +
"{rc} > {left} AND {rc} < {right} AND " +
"{dc} > {bottom} AND {dc} < {top} ")
return query.format(**fmt)
[docs]def query_tap_catalog(ra=165.86, dec=34.829694, radius=3., corners=None,
max_wait=20,
db='ls_dr9.tractor', columns=['*'], extra='',
rd_colnames=['ra', 'dec'],
tap_url='https://datalab.noirlab.edu/tap',
max=1000000, clean_xml=True, verbose=True,
des=False, gaia=False, nsc=False, vizier=False,
skymapper=False,
hubble_source_catalog=False, tap_kwargs={},
**kwargs):
"""Query NOAO Catalog holdings
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
corners : 4-tuple, `~astropy.wcs.WCS` or None
ra_min, ra_max, dec_min, dec_max of a query box to use instead of
`radius`. Or if a `WCS` object, get limits from the
`~astropy.wcs.WCS.calc_footprint` method
db : str
Parent database (https://datalab.noirlab.edu/query.php).
columns : list of str
List of columns to output. Default ['*'] returns all columns.
extra : str
String to add to the end of the positional box query, e.g.,
'AND mag_auto_i > 16 AND mag_auto_i < 16.5'.
rd_colnames : str, str
Column names in `db` corresponding to ra/dec (degrees).
tap_url : str
TAP hostname
des : bool
Query `des_dr1.main` from NOAO.
gaia : bool
Query `gaiadr2.gaia_source` from http://gea.esac.esa.int.
nsc : bool
Query the NOAO Source Catalog (Nidever et al. 2018), `nsc_dr1.object`.
vizier : bool
Use the VizieR TAP server at http://tapvizier.u-strasbg.fr/TAPVizieR/tap, see http://tapvizier.u-strasbg.fr/adql/about.html.
hubble_source_catalog : bool
Query the Hubble Source Catalog (v3). If no 'NumImages' criteria is
found in `extra`, then add an additional requirement:
>>> extra += 'AND NumImages > 1'
Returns
-------
table : `~astropy.table.Table`
Result of the query
"""
from astroquery.utils.tap.core import TapPlus
# DES DR1
if des:
if verbose:
print('Query DES DR1 from NOAO')
db = 'des_dr1.main'
tap_url = 'https://datalab.noirlab.edu/tap'
# NOAO source catalog, seems to have some junk
if nsc:
if verbose:
print('Query NOAO source catalog')
db = 'nsc_dr1.object'
tap_url = 'https://datalab.noirlab.edu/tap'
extra += ' AND nsc_dr1.object.flags = 0'
# GAIA DR2
if gaia:
if verbose:
print('Query GAIA DR2 from ESA')
db = 'gaiadr2.gaia_source'
tap_url = 'http://gea.esac.esa.int/tap-server/tap'
# VizieR TAP server
if vizier:
if verbose:
print('Query {0} from VizieR TAP server'.format(db))
tap_url = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap'
rd_colnames = ['RAJ2000', 'DEJ2000']
if skymapper:
if verbose:
print('Query {0} from VizieR TAP server'.format(db))
tap_url = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap'
rd_colnames = ['RAICRS', 'DEICRS']
if hubble_source_catalog:
if db is None:
db = 'dbo.SumPropMagAutoCat'
elif 'dbo' not in db:
db = 'dbo.SumPropMagAutoCat'
tap_url = 'http://vao.stsci.edu/HSCTAP/tapservice.aspx'
rd_colnames = ['MatchRA', 'MatchDec']
if 'NumImages' not in extra:
extra += 'AND NumImages > 1'
tap = TapPlus(url=tap_url, **tap_kwargs)
query = gen_tap_box_query(ra=ra, dec=dec, radius=radius, max=max,
db=db, columns=columns,
rd_colnames=rd_colnames,
corners=corners)
job = tap.launch_job(query+extra, dump_to_file=True, verbose=verbose)
try:
table = job.get_results()
if clean_xml:
if hasattr(job, 'outputFile'):
jobFile = job.outputFile
else:
jobFile = job.get_output_file()
os.remove(jobFile)
# Provide ra/dec columns
for c, cc in zip(rd_colnames, ['ra', 'dec']):
if (c in table.colnames) & (cc not in table.colnames):
table[cc] = table[c]
table.meta['TAPURL'] = tap_url, 'TAP URL'
table.meta['TAPDB'] = db, 'TAP database name'
table.meta['TAPQUERY'] = query+extra, 'TAP query'
table.meta['RAQUERY'] = ra, 'Query central RA'
table.meta['DECQUERY'] = dec, 'Query central Dec'
table.meta['RQUERY'] = radius, 'Query radius, arcmin'
if hubble_source_catalog:
for col in table.colnames:
if table[col].dtype == 'object':
print('Reformat column: {0}'.format(col))
strcol = list(table[col])
table[col] = strcol
except:
if hasattr(job, 'outputFile'):
jobFile = job.outputFile
else:
jobFile = job.get_output_file()
print('Query failed, check {0} for error messages'.format(jobFile))
table = None
return table
# Limit Hubble Source Catalog query to brighter sources in limited bands
HSCv3_FILTER_LIMITS = {'W3_F160W': 23.5,
'W3_F140W': 23.5,
'W3_F125W': 23.5,
'W3_F110W': 23.5,
'W3_F098M': 23.5,
'W3_F105W': 23.5,
'A_F814W': 23.5,
'W3_F814W': 23.5,
'A_F606W': 23.5,
'W3_F606W': 23.5,
'A_F850LP': 23.5,
'W3_F850LP': 23.5,
'A_F775W': 23.5,
'W3_F775W': 23.5}
HSCv3_COLUMNS = ['MatchRA', 'MatchDec', 'CI', 'CI_Sigma',
'KronRadius', 'KronRadius_Sigma', 'Extinction',
'TargetName', 'NumImages', 'NumFilters', 'NumVisits',
'DSigma']
[docs]def get_hubble_source_catalog(ra=0., dec=0., radius=3, corners=None, max=int(1e7), extra=' AND NumImages > 0', kron_max=0.45, dsigma_max=100, clip_singles=10*u.arcsec, verbose=True, columns=HSCv3_COLUMNS, filter_limits=HSCv3_FILTER_LIMITS):
"""
Query NOAO Source Catalog, which is aligned to GAIA DR1.
The default `extra` query returns well-detected sources in red bands.
filter_limits : query on individual HSC filter magnitudes
"""
import astropy.table
msg = 'Query NOAO Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')'
print(msg.format(ra=ra, dec=dec, radius=radius))
if kron_max is not None:
extra += ' AND KronRadius < {0}'.format(kron_max)
if dsigma_max is not None:
extra += ' AND DSigma < {0}'.format(dsigma_max)
if filter_limits is not None:
limit_list = ['{0} < {1}'.format(f, filter_limits[f])
for f in filter_limits]
filter_selection = ' AND ({0})'.format(' OR '.join(limit_list))
extra += filter_selection
columns += [f for f in filter_limits]
db = 'dbo.SumPropMagAutoCat p join dbo.SumMagAutoCat m on p.MatchID = m.MatchID'
else:
db = 'dbo.SumPropMagAutoCat'
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius,
corners=corners, max=max,
extra=extra, hubble_source_catalog=True,
verbose=verbose, db=db, columns=columns)
if clip_singles not in [None, False]:
rr = tab['NumImages'] > 1
if (rr.sum() > 0) & ((~rr).sum() > 0):
r0, r1 = tab[rr], tab[~rr]
idx, dr = utils.GTable(r0).match_to_catalog_sky(r1)
new = dr > clip_singles
xtab = astropy.table.vstack([r0, r1[new]])
if verbose:
msg = ('HSCv3: Remove {0} NumImages == 1 sources ' +
'with tolerance {1}')
print(msg.format((~new).sum(), clip_singles))
return tab
[docs]def get_nsc_catalog(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (rerr < 0.08 OR ierr < 0.08 OR zerr < 0.08) AND raerr < 0.2 AND decerr < 0.2', verbose=True):
"""
Query NOAO Source Catalog, which is aligned to GAIA DR1.
The default `extra` query returns well-detected sources in red bands.
"""
msg = 'Query NOAO Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')'
print(msg.format(ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners,
extra=extra, nsc=True, verbose=verbose, max=max)
return tab
[docs]def get_desdr1_catalog(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (e_rmag < 0.15 OR e_imag < 0.15)', verbose=True):
"""
Query DES DR1 Catalog from Vizier
The default `extra` query returns well-detected sources in one or more
red bands.
"""
msg = 'Query DES Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')'
print(msg.format(ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners,
extra=extra, db='"II/357/des_dr1"',
vizier=True,
verbose=verbose, max=max)
return tab
def get_desdr1_catalog_old(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (magerr_auto_r < 0.15 OR magerr_auto_i < 0.15)', verbose=True):
"""
Query DES DR1 Catalog.
The default `extra` query returns well-detected sources in one or more
red bands.
"""
msg = 'Query DES Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')'
print(msg.format(ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners,
extra=extra, des=True, verbose=verbose, max=max)
return tab
[docs]def get_vhs_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='', table='vhs_dr4'):
"""
VEXAS DR2 from vizier
https://vizier.cds.unistra.fr/viz-bin/VizieR?-source=II/369
table: 'vexasds', 'vexasps', 'vexassm'
"""
msg = 'Query VHS DR4 {table} catalog ({ra},{dec},{radius})'
print(msg.format(table=table, ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2,
corners=corners, extra=extra,
vizier=True,
db=f'"II/359/{table}"',
verbose=verbose, max=max_records)
return tab
[docs]def get_vexas_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='', table='vexasdes'):
"""
VEXAS DR2 from vizier
https://vizier.cds.unistra.fr/viz-bin/VizieR?-source=II/369
table: 'vexasds', 'vexasps', 'vexassm'
"""
msg = 'Query VEXAS DR2 {table} catalog ({ra},{dec},{radius})'
print(msg.format(table=table, ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2,
corners=corners, extra=extra,
vizier=True,
db=f'"II/369/{table}"',
verbose=verbose, max=max_records)
return tab
[docs]def get_skymapper_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra=''):
"""
Get Skymapper DR1 from Vizier
"""
msg = 'Query Skymapper DR1 catalog ({ra},{dec},{radius})'
print(msg.format(ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2,
corners=corners, extra=extra,
skymapper=True, db='"II/358/smss"',
verbose=verbose, max=max_records)
return tab
def get_legacysurveys_catalog(ra=0., dec=0., radius=3., verbose=True, db='ls_dr9.tractor', sn_lim=('r',10), **kwargs):
"""
Query LegacySurveys TAP catalog
"""
if verbose:
msg = 'Query LegacySurveys ({db}) catalog ({ra},{dec},{radius:.2f})'
print(msg.format(ra=ra, dec=dec, radius=radius, db=db))
if (sn_lim is not None):
if len(sn_lim) >= 2:
b = sn_lim[0]
extra = f' AND flux_{b}*SQRT(flux_ivar_{b}) > {sn_lim[1]}'
if 'extra' in kwargs:
kwargs['extra'] += extra
else:
kwargs['extra'] = extra
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, db=db, **kwargs)
# if (sn_lim is not None) & (len(tab) > 1):
# if len(sn_lim) >= 2:
# sn = tab[f'flux_{sn_lim[0]}']
# sn *= np.sqrt(tab[f'flux_ivar_{sn_lim[0]}'])
#
# valid = sn > sn_lim[1]
# return tab[valid]
return tab
[docs]def get_panstarrs_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='AND "II/349/ps1".e_imag < 0.2 AND "II/349/ps1".e_RAJ2000 < 0.15 AND "II/349/ps1".e_DEJ2000 < 0.15'):
"""
Get PS1 from Vizier
"""
msg = 'Query PanSTARRS catalog ({ra},{dec},{radius})'
print(msg.format(ra=ra, dec=dec, radius=radius))
tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2,
corners=corners, extra=extra,
vizier=True, db='"II/349/ps1"', verbose=verbose,
max=max_records)
return tab
[docs]def get_radec_catalog(ra=0., dec=0., radius=3., product='cat', verbose=True, reference_catalogs=['GAIA', 'LS_DR9', 'PS1', 'Hubble', 'NSC', 'SDSS', 'WISE', 'DES'], use_self_catalog=False, **kwargs):
"""Decide what reference astrometric catalog to use
First search SDSS, then WISE looking for nearby matches.
Parameters
----------
ra, dec : float
Center of the query region, decimal degrees
radius : float
Radius of the query, in arcmin
product : str
Basename of the drizzled product. If a locally-created catalog with
filename that startswith `product` is found, use that one instead of
the external (low precision) catalogs so that you're matching
HST-to-HST astrometry.
reference_catalogs : list
Order in which to query reference catalogs. Options are 'GAIA',
'PS1' (STScI PanSTARRS), 'SDSS', 'WISE', 'NSC' (NOAO Source Catalog),
'DES' (Dark Energy Survey DR1), 'Hubble' (Hubble Source Catalog v3),
'LS_DR9' (LegacySurveys DR9).
Returns
-------
radec : str
Filename of the RA/Dec list derived from the parent catalog
ref_catalog : str, {'SDSS', 'WISE', 'VISIT'}
Provenance of the `radec` list.
"""
query_functions = {'SDSS': get_sdss_catalog,
'GAIA_TAP': get_gaia_DR2_catalog,
'PS1': get_panstarrs_catalog,
'WISE': get_irsa_catalog,
'2MASS': get_twomass_catalog,
'GAIA_Vizier': get_gaia_DR2_vizier,
'GAIA': get_gaia_DR2_vizier,
'NSC': get_nsc_catalog,
'DES': get_desdr1_catalog,
'Hubble': get_hubble_source_catalog,
'Skymapper': get_skymapper_catalog,
'VEXAS': get_vexas_catalog,
'LS_DR9': get_legacysurveys_catalog}
# Try queries
has_catalog = False
ref_catalog = 'None'
ref_cat = []
for ref_src in reference_catalogs:
try:
if ref_src == 'GAIA':
try:
ref_cat = query_functions[ref_src](ra=ra, dec=dec,
radius=radius, use_mirror=False)
except:
try:
ref_cat = query_functions[ref_src](ra=ra, dec=dec,
radius=radius)
except:
ref_cat = False
# Try GAIA mirror at Heidelberg
if ref_cat is False:
ref_cat = query_functions[ref_src](ra=ra, dec=dec,
radius=radius, use_mirror=True)
else:
ref_cat = query_functions[ref_src](ra=ra, dec=dec,
radius=radius)
# #
# ref_cat = query_functions[ref_src](ra=ra, dec=dec,
# radius=radius)
valid = np.isfinite(ref_cat['ra']+ref_cat['dec'])
ref_cat = ref_cat[valid]
if len(ref_cat) < 2:
raise ValueError
table_to_regions(ref_cat, output='{0}_{1}.reg'.format(product,
ref_src.lower()))
ref_cat['ra', 'dec'].write('{0}_{1}.radec'.format(product,
ref_src.lower()),
format='ascii.commented_header',
overwrite=True)
radec = '{0}_{1}.radec'.format(product, ref_src.lower())
ref_catalog = ref_src
has_catalog = True
if len(ref_cat) > 0:
break
except:
print('{0} query failed'.format(ref_src))
has_catalog = False
if (ref_src.startswith('GAIA')) & ('date' in kwargs) & has_catalog:
if kwargs['date'] is not None:
if 'date_format' in kwargs:
date_format = kwargs['date_format']
else:
date_format = 'mjd'
gaia_tbl = ref_cat # utils.GTable.gread('gaia.fits')
coo = get_gaia_radec_at_time(gaia_tbl, date=kwargs['date'],
format=date_format)
coo_tbl = utils.GTable()
coo_tbl['ra'] = coo.ra
coo_tbl['dec'] = coo.dec
ok = np.isfinite(coo_tbl['ra']) & np.isfinite(coo_tbl['dec'])
coo_tbl.meta['date'] = kwargs['date']
coo_tbl.meta['datefmt'] = date_format
msg = 'Apply observation ({0},{1}) to GAIA catalog'
print(msg.format(kwargs['date'], date_format))
table_to_regions(coo_tbl[ok], output='{0}_{1}.reg'.format(product,
ref_src.lower()))
coo_tbl['ra', 'dec'][ok].write('{0}_{1}.radec'.format(product,
ref_src.lower()),
format='ascii.commented_header',
overwrite=True)
if not has_catalog:
return False
# WISP, check if a catalog already exists for a given rootname and use
# that if so.
cat_files = glob.glob('-f1'.join(product.split('-f1')[:-1]) + '-f*.cat*')
if (len(cat_files) > 0) & (use_self_catalog):
ref_cat = utils.GTable.gread(cat_files[0])
root = cat_files[0].split('.cat')[0]
ref_cat['X_WORLD', 'Y_WORLD'].write('{0}.radec'.format(root),
format='ascii.commented_header',
overwrite=True)
radec = '{0}.radec'.format(root)
ref_catalog = 'VISIT'
if verbose:
msg = '{0} - Reference RADEC: {1} [{2}] N={3}'
print(msg.format(product, radec, ref_catalog, len(ref_cat)))
return radec, ref_catalog