In [1]:
%matplotlib inline
#Lets you play a movie inline in the ipython notebook 
%matplotlib nbagg
from astropy import units as u
from astroquery.xmatch import XMatch
from astroquery.vizier import Vizier
import matplotlib.pyplot as plt
import numpy as np
from astropy.table import Table, hstack, join
from astropy.io import ascii
import pyfits
In [2]:
#No limit on rows returned
Vizier.ROW_LIMIT = -1
In [3]:
#main Hipparcos catalog
dat = Vizier.get_catalogs('I/311/hip2')
dat = dat[0]
print dat
 _RAJ2000   _DEJ2000   HIP   n_HIP  Sn  So ...  B-V   e_B-V   V-I   HIP1 Phot
   deg        deg                          ...  mag    mag    mag            
---------- ---------- ------ ----- --- --- ... ------ ------ ------ ---- ----
  0.000901   1.089010      1         5   0 ...  0.482  0.025  0.550 HIP1 Phot
  0.004269 -19.498841      2        75   4 ...  0.999  0.002  1.040 HIP1 Phot
  0.005021  38.859278      3         5   0 ... -0.019  0.004  0.000 HIP1 Phot
  0.008629 -51.893546      4         5   0 ...  0.370  0.009  0.430 HIP1 Phot
  0.009971 -40.591205      5         5   0 ...  0.902  0.013  0.900 HIP1 Phot
  0.018690   3.946453      6         5   0 ...  1.336  0.020  1.550 HIP1 Phot
  0.022011  20.036112      7         5   0 ...  0.740  0.020  0.790 HIP1 Phot
  0.027348  25.886459      8     P   5   0 ...  1.102  0.051  3.920 HIP1 Phot
  0.035321  36.585959      9         5   0 ...  1.067  0.023  1.030 HIP1 Phot
  0.036414 -50.866974     10         5   0 ...  0.489  0.011  0.560 HIP1 Phot
       ...        ...    ...   ... ... ... ...    ...    ...    ...  ...  ...
102.462666  66.358178 120248         5   0 ...  1.500  0.020  0.930 HIP1 Phot
317.507030  -1.864537 120250         5   0 ...  0.990  0.005  1.030 HIP1 Phot
150.519995  79.705974 120276         5   0 ...  0.000  0.000  0.000 HIP1 Phot
254.420736  35.269710 120290         5   0 ...  0.718  0.037  0.710 HIP1 Phot
346.687308 -66.075232 120306        15   3 ...  1.250  0.015  1.210 HIP1 Phot
206.398681  17.742433 120313         5   0 ...  1.430  0.020  1.400 HIP1 Phot
119.382239 -60.630880 120401         5   0 ...  0.104  0.007  0.120 HIP1 Phot
119.448626 -60.609716 120402         5   0 ... -0.022  0.010  0.000 HIP1 Phot
119.454885 -60.683622 120403         5   0 ...  0.061  0.007  0.070 HIP1 Phot
119.512126 -60.614778 120404         5   0 ... -0.062  0.018 -0.040 HIP1 Phot
Length = 117955 rows
/home/jakub/anaconda/lib/python2.7/site-packages/astroquery/vizier/core.py:567: UserWarning: VOTABLE parsing raised exception: 
  warnings.warn("VOTABLE parsing raised exception: {0}".format(ex))
In [4]:
#Late G and K stars with parallax errors less than 100% 
keep, = np.where((dat['B-V'] > 0.65) & (dat['B-V'] < 0.80 )& ( dat['e_Plx']/dat['Plx'] < 1.0) )# &
print keep.size
7201
In [5]:
#rename columns for querying in Xquery
cut = dat[keep]
cut.rename_column('_RAJ2000','ra')
cut.rename_column('_DEJ2000','dec')
#Now Match to 2Mass Catalog
finder = cut['ra','dec']
#Mask needs to be False or Xquery won't work
finder = Table(finder,masked=False)
In [6]:
#use Xmatch to find the 2MASS ID and IR photometry for these HIP stars
tab = XMatch.query(cat1=finder,cat2='vizier:II/246/out',max_distance=0.5*u.arcsec,colRA1='ra',colDec1='dec')
print 'Total matches = {0:4.0f}, Unique Matches {1:4.0f}'.format(tab['2MASS'].size,np.unique(tab['2MASS']).size)
Total matches = 6873, Unique Matches 6873
/home/jakub/anaconda/lib/python2.7/site-packages/astropy/table/table.py:1990: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  result = self.as_array() == other
In [7]:
#print tab
#remove objects without 2MASS phot. (which there are none in this iteration)
keepers, = np.where(tab['2MASS'].mask == False)
tab = tab[keepers]
In [8]:
#Unfortunately Xmatch cuts the rows it doesn't find 
from astropy.coordinates import match_coordinates_sky, SkyCoord
c1 = SkyCoord(ra=cut['ra'],dec=cut['dec'],unit='deg')
c2 = SkyCoord(ra=tab['ra'],dec=tab['dec'],unit='deg')
#Get the indexes, which we can match
idx1,idx2, idis1, _ = c1.search_around_sky(c2, 0.5*u.arcsec)
#Create empty string array
cut['2MASS'] = ' '*18
#Fill string with values match in the previous lines using the index match from c2 to c1
cut['2MASS'][idx2] = tab['2MASS']
#Join Tables on 2MASS ID
tab = join(cut,tab,keys=['2MASS'])
#Make sure RA and Dec line up after match (there is an ra and dec component in cut and tab,
# so joining appends _1 and _2, respectively)
ddec = (tab['dec_1']-tab['dec_2'])
dra = tab['ra_1']-tab['ra_2']
print (np.sqrt(ddec**2+(dra*np.cos(np.radians(ddec)))**2)).max()
0.0
In [9]:
#Now let us check if the stars have apogee spectra (Apogee targeted HIP stars for calibration early on the 1m)
apdir = 'http://data.sdss3.org/sas/dr12/apogee/spectro/redux/r5/stars/apo1m/hip/'
#This is the location of the Apogee HIP stars list on the Apogee d12 server 
apohip = ascii.read(apdir+'r5_stars_apo1m_hip.sha1sum',
                    names=['num','2MASSfits'],exclude_names=['num'])
Downloading http://data.sdss3.org/sas/dr12/apogee/spectro/redux/r5/stars/apo1m/hip/r5_stars_apo1m_hip.sha1sum [Done]
In [10]:
#remove all repeat spectra because we are only looking for certain reduced apogee spectra
reducer = np.core.defchararray.find(apohip['2MASSfits'],'apStar-r5')
#Only keep file names which star with 'apStar-r5'
reduced, = np.where(reducer == 0)
apohipfile = apohip[reduced]
In [11]:
#Create 2MASS ID list from Apogee file names
apohipfile['2MASS'] = apohipfile['2MASSfits']
for j,i in enumerate(apohipfile['2MASSfits']):
    apohipfile['2MASS'][j] = str(i[12:-11])
#Fix formating so we can match to tab['2MASS'], which is a string column
apohipfile = Table(apohipfile,names=['2MASSfits','2MASS',],dtype=('S48','S16'))
In [12]:
#Combine the Apogee and Hipparcos list
fhipap = join(tab,apohipfile,keys=['2MASS'])
print 'Final number of HIP stars with Apogee Spectra {0:4.0f}'.format(fhipap['2MASS'].size)
#You can now also double check 2MASS and HIP IDs in simbad and they are the same (i.e. the Xmatch query worked)
print fhipap['HIP','2MASS']
Final number of HIP stars with Apogee Spectra    7
 HIP        2MASS      
------ ----------------
 16226 03290380+0314553
 19740 04135637+0915497
 84934 17213124+2845290
 99170 20080116+4223061
101916 20390779+1005104
108090 21535776+0651530
110991 22291029+5824549
In [13]:
#Now let's grab the spectra from the SDSS archive
spec = []
for j,i in enumerate(fhipap['2MASSfits']):
    spec.append(pyfits.open(apdir+i)[1])
In [14]:
#add tick marks to plots
def fancy_plot(ax):
    #Turn minor ticks on
    ax.minorticks_on()
    #set the width of the ticks
    ax.tick_params(which='both',width=1)
    #set the length of the major ticks
    ax.tick_params(which='major',length=7)
    #set length of the minor ticks
    ax.tick_params(which='minor',length=3)
In [15]:
#And for the final scene, a movie
import matplotlib.animation as an
#First create a plot with nothing on it
fig, ax = plt.subplots()
l, = ax.plot(10**(spec[0].header['CRVAL1']+spec[0].header['CDELT1']*np.arange(spec[0].data.size)),
             spec[0].data,color=None)
ax.set_xlabel('Wavelength [$\AA$]')
ax.set_ylabel('Flux  [10$^{-17}$ erg/s/cm$^2$/$\AA$]')
fancy_plot(ax)

#create a function, which updates the plot information
def update_movie(num):
    spec0 = spec[num].data
    wave0 = 10.**(spec[num].header['CRVAL1'] + spec[num].header['CDELT1']*np.arange(spec0.size))
    l.set_data(wave0,spec0)
    l.set_color('black')
    ax.set_ylim([0,spec0.max()])

ani = an.FuncAnimation(fig,update_movie,frames=7)
#must have ffmpeg or mencoder installed to save (this step can be nontrivial) 
ani.save('test.mp4',writer='ffmpeg',fps=.5)
#press the X at the top right hand corner to stop video playing
In [ ]: