Defining some utilitary functions and importing some modules
In [1]:
def query_TAP(tap_endpoint, adql_query, table_to_upload=None):
"""
Query a TAP service (designated by its tap_endpoint)
with a given ADQL query
Query is performed synchronously
Return an AstroPy Table object
"""
import requests
from astropy.table import Table
from astropy.io.votable import parse_single_table
import os
import tempfile
import warnings
r = requests.post(tap_endpoint + '/sync', data={'query': adql_query, 'request': 'doQuery', 'lang': 'adql', 'format': 'votable', 'phase': 'run'})
with warnings.catch_warnings():
warnings.simplefilter("ignore")
tmp_vot = tempfile.NamedTemporaryFile(delete = False)
with open(tmp_vot.name, 'w') as h:
for line in r.iter_lines():
if line:
h.write(line.decode(r.encoding)+'\n')
table = parse_single_table(tmp_vot.name).to_table()
# finally delete temp files
os.unlink(tmp_vot.name)
return table
from astroquery.xmatch import XMatch
import types
import six
from astropy.io import ascii
# monkey patching XMatch
def patched_is_table_available(self, table_id):
if isinstance(table_id, six.string_types) and (table_id[:7] == 'vizier:'):
table_id = table_id[7:]
if not isinstance(table_id, six.string_types):
return False
return table_id in self.get_available_tables()
XMatch.is_table_available = types.MethodType(patched_is_table_available, XMatch)
from astropy import units as u
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib
def plot_scatter_density(xdata, ydata, xlabel, ylabel, title, xlim=None, ylim=None, cmap='viridis', invert_yaxis = True, s=2, grid=False):
from scipy.stats import gaussian_kde
import numpy as np
x = np.reshape(np.array(xdata, copy=False).astype('float'), (len(xdata)))
y = np.reshape(np.array(ydata, copy=False).astype('float'), len(ydata))
xy = np.vstack([x,y])
z = gaussian_kde(xy)(xy)
# Sort the points by density, so that the densest points are plotted last
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
w = h = 8
fig, ax = plt.subplots(figsize = (w, h))
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if xlim:
ax.set_xlim(*xlim)
if ylim:
ax.set_ylim(*ylim)
ax.set_title(title)
ax.scatter(x, y, c=z, s=s, edgecolor='', cmap=cmap)
if invert_yaxis:
ax.invert_yaxis()
ax.grid(grid)
# disable warnings as not to clutter notebook
import warnings
warnings.simplefilter("ignore")
Let us query table I/337/gaia 20 arcmin around Messier 44
In [2]:
from astroquery.vizier import Vizier
result = Vizier.query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia')
print(result[0])
By default, output is limited to 50 rows. Let's change this and redo the query.
In [3]:
Vizier.ROW_LIMIT = 100000
result = Vizier.query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia')
print(result[0]['RA_ICRS', 'DE_ICRS', 'Source', 'Plx'])
We add now a constraint to keep only sources with a parallax > 5 mas
In [4]:
result = Vizier(column_filters={"Plx":">5"}).query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia', )
print(result[0]['RA_ICRS', 'DE_ICRS', 'Source', 'Plx'])
Send result to VO tools through SAMP
In [5]:
# TODO
Retrieve list of RRLyrae and sort it by period
In [6]:
rrlyrae = Vizier.get_catalogs(catalog='I/337/rrlyrae')[0]
rrlyrae.sort('P1')
print(rrlyrae['Source', 'P1', 'RA_ICRS', 'DE_ICRS'])
Retrieve light curve points for Source = 4660711867541110784
In [7]:
lc_points = Vizier(column_filters={"Source":"4660711867541110784"}).get_catalogs(catalog='I/337/fov')[0]
print(lc_points['ObsTime', 'FG'])
Plot the light curve
In [8]:
x = lc_points['ObsTime']
y = lc_points['FG']
y_err = lc_points['e_FG']
plt.scatter(x, y, color='r')
plt.errorbar(x, y, yerr=y_err, color='r', linestyle='None')
plt.grid(True)
Fold the light curve, folded by the period (P1 = 0.21271385 d)
In [9]:
period = 0.21271385
#(time+8)%16 - 8
x = (lc_points['ObsTime'] % period) / period
y = lc_points['FG']
y_err = lc_points['e_FG']
plt.scatter(x, y, color='r')
plt.errorbar(x, y, yerr=y_err, color='r', linestyle='None')
plt.grid(True)
In [10]:
adql_query = """SELECT gaia.ra, gaia.dec, gaia.source_id, gaia.hip, gaia.phot_g_mean_mag+5*log10(gaia.parallax)-10 as g_mag_abs,
hip."B-V"
FROM "I/337/tgas" as gaia
inner join "I/311/hip2" as hip
on gaia.hip= hip.HIP
where gaia.parallax/gaia.parallax_error >= 5 and hip."e_B-V" > 0.0 and hip."e_B-V" <= 0.05 and
2.5/log(10)*gaia.phot_g_mean_flux_error/gaia.phot_g_mean_flux <= 0.05"""
tap_vizier_endpoint = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap/'
hr_data = query_TAP(tap_vizier_endpoint, adql_query)
print(hr_data)
In [11]:
plot_scatter_density(hr_data['B-V'], hr_data['g_mag_abs'],
title='Gaia DR1 HR diagram', xlabel = 'B-V', ylabel=r'$M_G$',
cmap='viridis', invert_yaxis = True, s=6, grid=True)
Load sample of AllWISE sources in Small Magellanic Cloud
In [12]:
from astropy.io import votable
allwise_sources = votable.parse_single_table('allwise-SMC-sample.vot').to_table()
print(len(allwise_sources))
Cross-match with Gaia DR1
In [13]:
allwise_X_gaia = XMatch.query(cat1 = allwise_sources,
cat2 = 'vizier:I/337/gaia',
max_distance = 1 * u.arcsec,
colRA1 = 'RAJ2000',
colDec1= 'DEJ2000', cache=False)
print(len(allwise_X_gaia))
Create color-color plot
In [ ]:
plot_scatter_density(allwise_X_gaia['phot_g_mean_mag'] - allwise_X_gaia['W1mag'],
allwise_X_gaia['W1mag'] - allwise_X_gaia['W4mag'], title='Color-color plot',
xlabel = 'G-W1', ylabel='W1-W4',
cmap='viridis', invert_yaxis = True, s=6, grid=True)
Retrieve SDSS DR9 and HST-R MOCs (coverage)
In [14]:
from mocpy import MOC
moc_sdss = MOC.from_ivorn('ivo://CDS/P/SDSS9/u', nside=512)
moc_hst = MOC.from_ivorn('ivo://CDS/P/HST/R', nside=512)
moc_sdss.plot('SDSS coverage')
moc_hst.degrade_to_order(6).plot('HST coverage')
Compute intersection
In [15]:
common_coverage = moc_sdss.intersection(moc_hst)
print('Common coverage represents %.5f%% of the sky' % (common_coverage.sky_fraction))
common_coverage.degrade_to_order(6).plot('Common coverage between SDSS and HST')
Retrieve Gaia DR1 sources in the coverage of the computed MOC
In [16]:
gaia_sources = common_coverage.query_vizier_table('GAIA-DR1', max_rows=200000)
print(gaia_sources['ra', 'dec'])
In [17]:
stars_in_cluster = query_TAP('http://simbad.u-strasbg.fr/simbad/sim-tap',
"""SELECT ra, dec, pmra, pmdec, plx_value, plx_err FROM basic
WHERE otype_txt = '*iC'
and plx_value>1 and plx_value/plx_err > 5""")
print(stars_in_cluster)
We now query the xmatch service to find counterparts in TGAS (max search radius: 1 arcsec)
In [18]:
xmatch_result = XMatch.query(cat1 = stars_in_cluster,
cat2 = 'vizier:I/337/tgasptyc',
max_distance = 1 * u.arcsec,
colRA1 = 'ra',
colDec1= 'dec', cache=False)
print(xmatch_result['ra', 'dec', 'plx_value', 'Plx', 'plx_err', 'e_Plx'])
In [19]:
histvals, binvals, patches = plt.hist(xmatch_result['Plx'], bins=100, facecolor='g', alpha=0.75)
plt.xlabel('Plx')
plt.ylabel('count')
plt.title('Gaia parallax distribution for cluster members')
plt.grid(True)
In [20]:
histvals, binvals, patches = plt.hist(100 * abs(xmatch_result['Plx'] - xmatch_result['plx_value']) / xmatch_result['Plx'], bins=100, facecolor='g', alpha=0.75)
plt.xlabel('Percentage')
plt.ylabel('count')
plt.title('Relative difference between parallax in Simbad and parallax measures by Gaia')
plt.grid(True)
In [21]:
plot_scatter_density(xmatch_result['ra'] - xmatch_result['_RAJ2000'],
xmatch_result['dec'] - xmatch_result['_DEJ2000'],
r'$\Delta$RA', r'$\Delta$dec', 'Position differences between Simbad and Gaia', xlim=[-5e-5, 5e-5], ylim=[-5e-5, 5e-5],
cmap='viridis', invert_yaxis = False, s=6, grid=True)