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")

 Query Gaia tables in VizieR

 Query table I/337/gaia

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])


  _RAJ2000    _DEJ2000     _r    ...      ELON           ELAT     
    deg         deg       deg    ...      deg            deg      
----------- ----------- -------- ... -------------- --------------
130.1865314  19.3486212 0.328363 ... 127.5117262730   0.9888678171
130.1734120  19.3550861 0.319197 ... 127.4981014125   0.9919388024
130.1809087  19.3575670 0.318394 ... 127.5043015311   0.9961520570
130.1453588  19.3393683 0.330106 ... 127.4765479898   0.9699560692
130.1530759  19.3466008 0.323979 ... 127.4817313114   0.9788143258
130.1305209  19.3558631 0.312159 ... 127.4587835142   0.9823087783
130.1676459  19.3578603 0.315346 ... 127.4921305987   0.9932240029
130.1648495  19.3597605 0.312960 ... 127.4890925233   0.9943835979
130.1760724  19.3824866 0.293110 ... 127.4934939615   1.0190656351
130.1632544  19.3889285 0.284090 ... 127.4801523613   1.0221888630
        ...         ...      ... ...            ...            ...
130.1590402  19.3887521 0.283456 ... 127.4763549095   1.0209984334
130.1321700  19.3841594 0.284156 ... 127.4530305729   1.0100579748
130.1431663  19.4099695 0.259927 ... 127.4564370771   1.0376648584
130.1041954  19.3862586 0.280463 ... 127.4269815665   1.0053221432
130.1275363  19.4080501 0.259942 ... 127.4426783635   1.0320291265
130.1312447  19.4068515 0.261505 ... 127.4463670167   1.0317674990
130.1322694  19.4155339 0.252994 ... 127.4450744922   1.0404074369
130.0343128  19.3635230 0.309430 ... 127.3690687801   0.9664618369
130.0412745  19.3670033 0.304760 ... 127.3745279590   0.9715071170
130.0569750  19.3901419 0.279510 ... 127.3829236323   0.9976661910
Length = 50 rows

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


   RA_ICRS        DE_ICRS           Source       Plx
     deg            deg                          mas
-------------- -------------- ------------------ ---
130.1865313919  19.3486212439 661215094304448640  --
130.1734120176  19.3550860525 661215163021682304  --
130.1809086744  19.3575670451 661215163021683584  --
130.1453588017  19.3393683410 661215231741149824  --
130.1530758561  19.3466007875 661215266100770176  --
130.1305209099  19.3558630743 661215300460507904  --
130.1676459265  19.3578602741 661215369180114176  --
130.1648494742  19.3597605175 661215369180116096  --
130.1760723589  19.3824865661 661215472259347584  --
130.1632544126  19.3889284843 661215472259355008  --
           ...            ...                ... ...
130.3072612277  19.9280473548 661401362739275648  --
130.2779935036  19.9267150742 661401392804351616  --
130.2936865771  19.9351334294 661401427163456896  --
130.2921118352  19.9366005215 661401427164102144  --
130.2791825435  19.9420699604 661401598962153472  --
130.2639921413  19.9441509879 661401667682246656  --
130.2641519444  19.9601174381 661413453071905664  --
129.8526759744  19.9027306987 664318500231186816  --
129.8529399648  19.9034999794 664318500233554432  --
129.8750213384  19.9239969697 664318912548085120  --
Length = 1684 rows

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


   RA_ICRS        DE_ICRS           Source        Plx  
     deg            deg                           mas  
-------------- -------------- ------------------ ------
130.2344037741  19.5802836152 661246018066598656   5.14
130.4260703292  19.6604830268 661252993093483392   6.06
129.9905508004  19.5414242443 661268248817325312   5.58
129.7877457088  19.5923409653 661284019938140032   6.12
130.0752439503  19.5319393680 661290754445957248   5.65
130.1123817124  19.5447617762 661291132403077760   5.38
130.0918665674  19.6698786868 661297076637809024   5.37
130.1798638364  19.7192651397 661300546971382016   5.22
130.2943301698  19.8294836433 661305838371084288   5.80
129.9275637670  19.7784005964 661311439008441600   5.19
130.0976146683  19.8349117737 661316902206838528   6.69
130.2370701390  19.9347999163 661319547906689024   5.60
130.0637997945  19.9942491067 661323636715553152   5.13
130.3072283373  19.9219407004 661401358443742720   5.93

Send result to VO tools through SAMP


In [5]:
# TODO

 Query RRLyrae table and retrieve a light curve

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


       Source            P1        RA_ICRS        DE_ICRS    
                         d           deg            deg      
------------------- ----------- -------------- --------------
4660711867541110784  0.21271385  84.0752731841 -65.7055413718
4660748774181753856  0.22183517  82.6928520423 -65.6441506678
5284152449583643648  0.22921618  90.3248461932 -66.8247183205
4663058229629271168  0.23316532  73.6300760013 -65.5173344610
5283973538437208320  0.23487892  91.7268314090 -66.8423299446
4662090178378689408  0.23565586  74.8057550936 -66.8600788745
4662010395066656384  0.23665484  77.0341004548 -66.4998070287
4663521502000875136  0.23766071  78.9541272098 -65.8416381584
4663686256931717376  0.23839866  78.2212926200 -64.6919818517
5276543829112460800  0.23944377 127.4767181695 -64.5756518278
                ...         ...            ...            ...
4675028024310945920  0.81744992  65.5170999930 -65.4086251949
4663900764773575680  0.82574879  76.6177392905 -64.8764986300
4663498201802520320  0.83220840  78.1237232923 -66.0479176389
4661920200755800064  0.83261467  75.9248629073 -66.6736981173
4659771819441025408  0.83377605  87.2597355588 -66.3909858124
4669806340510840320  0.86631745  60.9074757601 -65.1939895537
5283804144909441792  0.88004899  98.7339638672 -65.5468544963
4660458537484993920  0.88319490  79.5959922195 -66.5696616837
4662753626256671360  0.88504964  68.6887038774 -66.0138459756
4674099211863344256  0.89576731  52.7199112308 -62.5219921307
Length = 2595 rows

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


ObsTime        FG      
   d          e-/s     
------- ---------------
1687.99       480.50937
1688.24       491.64636
1688.49       437.73251
1688.74       395.67110
1688.99       404.05778
1689.24       404.28456
1689.49       495.58290
1689.74       475.91021
1689.99       418.01522
1690.24       394.27661
    ...             ...
1694.42       468.96957
1694.49       372.29648
1694.67       416.21899
1694.99       428.27182
1737.04       395.92665
1737.11       458.84048
1859.97       401.80493
1894.44       394.53520
1957.23       377.17774
1985.22       431.99535
Length = 39 rows

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)


Recreate HR diagram from Gaia DR1 paper

Get the data: query TAP VizieR


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)


    ra [1]       dec [1]        source_id [1]    hip [1]  g_mag_abs [1]  B-V [1]
     deg           deg                                         mag         mag  
------------- -------------- ------------------- ------- --------------- -------
274.642708474 -34.4457791552 4044586886081427072   89717   3.34145583118   0.632
274.662089117 -34.1855022022 4044599358666578560   89723   4.38266679476   0.636
275.991461847     -33.248192 4044723740920214400   90166   3.50045800088   0.493
275.393388974 -33.5303012248 4044738859203193856   89984  -1.91175539311   1.507
275.232281712 -33.3335833936 4044756588829994880   89922   0.35290495509   0.167
277.794151994 -32.2903031563 4045071427110885632   90772    2.1402805086   0.328
 277.77443737 -32.2121014617 4045075241041842560   90765   3.64113053028    0.83
276.555949538 -33.1860260588 4045088984941128960   90356  0.907499196734   1.028
276.645951999 -32.7752725491 4045132862326991488   90386  0.800896219733    1.05
276.153582954 -32.2911423871 4045208694270287872   90229  -2.46840008672   1.389
          ...            ...                 ...     ...             ...     ...
268.659441731 -33.9563424381 4043113952821300352   87686  0.541340221751   0.078
268.599795185 -33.6193706139 4043139997502975744   87664   4.34491649668    0.63
269.005238901 -33.6282051387 4043145976097452288   87791   3.17013545645   0.381
269.593437843 -32.9096635448 4043290252633204736   87985 -0.892133668191   1.606
269.208471055  -32.688904035 4043357254123008640   87853 -0.710303613051    0.07
268.468650088 -32.4780552313 4043549668660759936   87613  0.297844520171    0.96
270.604912638 -31.4898393177 4043843169543028096   88336   2.25088140888   0.491
269.532853394 -31.5399912992 4043992119012755328   87966   3.45570480341   0.727
269.814708083 -31.4420404466 4043996723217687424   88072   3.49207773779   0.505
275.357669156 -34.7438942742 4044333173774920064   89970   3.09010822652   0.674
Length = 74817 rows

Plot the HR diagram


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)


Find Gaia couterparts for a list of AllWISEsources and draw a color-color diagram

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))


313342

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))


209576

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 Gaia sources within HST and SDSS coverage

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')


0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

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')


Common coverage represents 0.00049% of the sky
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

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


        ra                 dec        
       deg                 deg        
------------------ -------------------
     45.1784271971  2.5477382024000002
45.173226863099998  2.5691562219000001
45.178255202099997  2.5583101777000001
45.170533466599998  2.5612363282000001
45.200298225099999        2.5629864296
45.206767306800003  2.5837332433000002
45.152746384899999        2.5660274183
45.162382379599997  2.5777730894999999
45.154739107600001  2.5785405207999998
45.177905306299998        2.5910490701
               ...                 ...
323.37174070409998 -7.1343296060999997
323.32617711850003 -7.1445170831000002
323.31559257459998       -7.1490559908
323.31195597990001 -7.1426892760999996
    323.3653256618 -7.1236426100000001
323.36278926199998       -7.1185015879
323.33427399589999 -7.1260084710999996
323.34350891129998 -7.1260929279000003
    323.3364825365 -7.1285626640000004
323.35063311760001 -7.1154539418000002
Length = 119286 rows

 Query SIMBAD to retrieve stars member of a cluster and cross-match with TGAS

Using astroquery, we will retrieve SIMBAD stars members of a cluster and having a parallax measurement greater than 1 mas and an error smaller than 1/5th of the parallax value

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)


        ra                 dec         ...     plx_value       plx_err  
       deg                 deg         ...        mas            mas    
------------------ ------------------- ... ------------------ ----------
161.56899012567726   -64.5145598145577 ... 6.8499999999999996 0.20999999
51.912372155824507   49.59991488038947 ... 6.2199999999999998 0.46000001
 48.34947130671857  48.176952442309371 ... 8.5800000000000001       0.37
 54.77802151852503  24.702873017617172 ... 8.3200000000000003 0.79000002
66.060726124423113  21.736244737907192 ...              22.59 0.76999998
184.81643444292581  25.060243297651738 ... 6.4199999999999999 0.66000003
75.203695332014647  4.7331138426158796 ... 18.379999999999999  1.9400001
51.636065921855888  47.266422576498691 ... 4.1600000000000001       0.62
51.545221196674149  48.384040625007195 ... 6.2199999999999998       0.75
 186.6002684910917  27.268237616650165 ...              11.82 0.23999999
               ...                 ... ...                ...        ...
56.197494121107724 -71.703478734107406 ... 9.1300000000000008 0.56999999
128.79217841185607 -54.206188431387794 ... 4.1200000000000001 0.33000001
129.59964017006868 -53.721777193711247 ... 7.0599999999999996 0.43000001
138.59472358043055 -65.053948093652892 ... 6.3099999999999996       0.41
265.36513432472361  -53.80298816604045 ...              15.59       1.79
305.67183091100219  -78.96091174141533 ... 9.2300000000000004 0.36000001
306.70067777432962 -78.194736900703418 ... 6.7800000000000002 0.76999998
307.50825968817594 -78.906364714305454 ...              11.74 0.60000002
63.606877838721537  14.625034704847046 ... 21.620000000000001        1.1
188.11431864818221    28.0846215372977 ... 20.870000000000001       1.16
Length = 498 rows

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


      ra           dec       plx_value  Plx     plx_err     e_Plx
------------- -------------- --------- ----- -------------- -----
51.9123721558  49.5999148804      6.22  5.68 0.460000008345   0.4
48.3494713067  48.1769524423      8.58  9.97 0.370000004768  0.46
66.0607261244  21.7362447379     22.59 23.36 0.769999980927  0.32
184.816434443  25.0602432977      6.42  6.24 0.660000026226  0.45
 75.203695332  4.73311384262     18.38 19.25  1.94000005722  0.24
51.6360659219  47.2664225765      4.16  6.19 0.620000004768  0.35
51.5452211967   48.384040625      6.22  5.87           0.75  0.27
122.950714152 -12.6787998036      5.48  3.21  1.01999998093  0.26
52.7773740974  50.4814932869      7.58  6.64 0.550000011921   0.3
128.831023404  19.5900653065      4.93  5.24 0.550000011921  0.41
          ...            ...       ...   ...            ...   ...
147.795155141 -53.1829827912      7.15  7.82 0.709999978542  0.26
6.90935584206 -71.8826491879      11.3 11.22 0.699999988079  0.29
56.1974941211 -71.7034787341      9.13  9.52 0.569999992847  0.24
128.792178412 -54.2061884314      4.12  4.34 0.330000013113  0.29
 129.59964017 -53.7217771937      7.06   7.0 0.430000007153  0.23
 138.59472358 -65.0539480937      6.31  5.65 0.409999996424  0.22
305.671830911 -78.9609117414      9.23   9.1 0.360000014305  0.36
306.700677774 -78.1947369007      6.78  6.76 0.769999980927  0.26
307.508259688 -78.9063647143     11.74 10.98 0.600000023842  0.23
63.6068778387  14.6250347048     21.62 20.56  1.10000002384  0.28
188.114318648  28.0846215373     20.87 21.44  1.15999996662  0.37
Length = 387 rows

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)