Credits:
Data: ESA Gaia DR2: http://cdn.gea.esac.esa.int/Gaia/gdr2/
# function to transform from barycentric to galactic rest frame
v_sun = coord.Galactocentric.galcen_v_sun.to_cartesian().get_xyz().to_value()
def rv_to_gsr(ra,dec,rv):
icrs = coord.ICRS(ra=ra*u.deg, dec=dec*u.deg,
radial_velocity=rv*u.km/u.s)
gal = icrs.transform_to(coord.Galactic)
cart_data = gal.data.to_cartesian().get_xyz().to_value()
unit_vector = cart_data / np.linalg.norm(cart_data)
return rv+np.dot(v_sun,unit_vector)
print(v_sun)
coord.Galactocentric.galcen_v_sun.norm()
That is the sun's velocity around the galaxy, effectively what we have to correct for when considering objects at galactic scale.
I downloaded the full RV CSV files since the ADQL interface was lagging, you can find these files here: http://cdn.gea.esac.esa.int/Gaia/gdr2/gaia_source_with_rv/csv/.
This is the code I used to combine everything, and then take a random subsample
# %%time
# rv_files = glob.glob('/scratch/rag394/gaia/cdn.gea.esac.esa.int/Gaia/gdr2/gaia_source_with_rv/csv/*.csv')
# file_list = []
# for file_ in rv_files:
# print('concatenating file: ',file_)
# frame = pd.read_csv(file_,index_col=None)
# file_list.append(frame)
# df = pd.concat(file_list)
# df.to_csv('/scratch/rag394/DR2_rv_complete.csv')
# df.sample(1000000).to_csv('/scratch/rag394/DR2_rv_subsample.csv')
# data = df.sample(500000)
# data['rv_corrected'] = data2[['ra','dec','radial_velocity']].apply(lambda x: rv_to_gsr(x.ra,x.dec,x.radial_velocity), axis=1)
# data.to_csv('/scratch/rag394/gaia/rv_corrected_subsample.csv')
# df = pd.read_csv('/scratch/rag394/DR2_rv_complete.csv')
# df = pd.read_csv('/scratch/rag394/DR2_rv_subsample.csv',index_col=None)
df = pd.read_csv('/scratch/rag394/gaia/rv_corrected_subsample.csv',index_col=None)
# remove bad solutions to parallax fit by filtering out negative parallaxes:
df = df[df.parallax>0.]
# only take high snr parallaxes
df = df[(df.parallax/df.parallax_error)>5]
# calc parallax based distance (careful with small (< 30pc) objects!)
df['distance'] = 1000./df['parallax']
data = df # I filtered positive parallax and high snr before saving previous csv file for corrected sample
print(data.shape[0])
ra = data.ra.values
dec = data.dec.values
ra = coord.Angle(ra*u.degree)
ra = ra.wrap_at(180*u.degree)
dec = coord.Angle(dec*u.degree)
fig = plt.figure(figsize=(18,14))
ax = fig.add_subplot(111,projection='mollweide')
s = ax.scatter(ra.radian,dec.radian,c=data.rv_corrected,cmap='viridis',
s=.007, alpha=1., vmin=-150,vmax=150)
ax.grid(alpha=0.2)
#plt.colorbar(s,ax=ax);
ax.set_xticklabels([]);