#!/usr/bin/env python
# coding: utf-8

# # Milky Way Datasets for WWT Visualization

# In this notebook we will connect to the VizieR catalog service to pull Milky Way datasets of Globular Clusters, Open Clusters, Planetary Nebulae, Pulsars, Magnetars, Supernovae Remnants and OB Associations. As newer catalogs are made available the code can easily be updated. Catalogs are pulled from VizieR using Astroquery.
# 
# Mark SubbaRao<br> 
# msubbarao@adlerplanetarium.org<br>

# ###### Set up Astropy and Astroquery to connect to VizieR Catalog Service

# In[1]:


#Set up astropy and astroquery 
from astropy.table import Table,Column
from astropy.coordinates import SkyCoord, Distance, CartesianRepresentation
from astropy.io.votable import from_table, writeto
from astropy import units as u
from astroquery.vizier import Vizier
v = Vizier()
v.ROW_LIMIT = -1


# ###### Set up matplotlib for plotting

# In[2]:


get_ipython().run_line_magic('config', 'InlineBackend.rc = {}')
import matplotlib
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import seaborn


# ###### Connect to Worldwide Telescope and set up layer group

# In[3]:


#Connect to WWT and set up Layer Group
from pywwt.mods import *
wwt = WWTClient(host="127.0.0.1") #Can pass a IP address here if WWT is running on a remote machine
wwt.new_layer_group("Sky","Galactic Datasets")


# ## Planetary Nebulae

# Distances of Galactic Planetary Nebulae <a hhref="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?2008ApJ...689..194S&db_key=AST&nosetcookie=1">(Stanghellini+, 2008)</a><br>
# VizieR Catalog: <a href="http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/ApJ/689/194/table1&-out.max=50&-out.form=HTML%20Table&-out.add=_r&-out.add=_RAJ,_DEJ&-sort=_r&-oc.form=sexa">J/ApJ/689/194/table1</a><br>
# 563 Objects

# ###### Get Planetary Nebulae Catalog

# In[4]:


Cats = v.get_catalogs('J/ApJ/689/194/table1')
PNeCat=Cats[0]
PNeCat.rename_column('_RAJ2000', 'RA')
PNeCat.rename_column('_DEJ2000', 'dec')
del PNeCat['_RA','_DE','SimbadName']
PNeCoords=SkyCoord(PNeCat['RA'],PNeCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(PNeCat['Dist'],u.pc),frame='icrs')


# ###### Print Planetary Nebulae Catalog

# In[5]:


PNeCat 


# ###### Plot the catalog

# In[6]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Planetary Nebulae",fontsize=18)
ax.scatter(PNeCoords.galactic.l.wrap_at(180.*u.degree).radian,PNeCoords.galactic.b.radian,c="green")


# ###### Send planetary nebulae to WWT

# In[7]:


#Set up WWT layer
PNe_layer = wwt.new_layer("Galactic Datasets", "Planetary Nebulae", PNeCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"6",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:0:255:0",\
              "TimeSeries":"False"}
PNe_layer.set_properties(props_dict)
#Send data to WWT client
PNe_layer.update(data=PNeCat, purge_all=True, no_purge=False, show=True)


# ### Output a Digistar Style VOtable of Planetary Nebulae

# In[6]:


#Add the Cartesian coordinates as the first three columns in the data table
PNeCat.add_column(Column(PNeCoords.represent_as(CartesianRepresentation).x.to(u.pc),name='x_coord',meta={'ucd': 'pos.cartesian.x'}),0)
PNeCat.add_column(Column(PNeCoords.represent_as(CartesianRepresentation).y.to(u.pc),name='y_coord',meta={'ucd': 'pos.cartesian.y'}),1)
PNeCat.add_column(Column(PNeCoords.represent_as(CartesianRepresentation).z.to(u.pc),name='z_coord',meta={'ucd': 'pos.cartesian.z'}),2)
#Write VOtable for Digistar visualization
from astropy.io.votable import from_table, writeto
votable = from_table(PNeCat)
writeto(votable, "data/PNe.du")


# In[10]:


PNeCat.remove_columns(['x_coord', 'y_coord', 'z_coord'])


# ### Output a .speck file for Digital Sky / Uniview / Partiview

# In[18]:


#Add the Cartesian, Galactic coordinates as the first three columns in the data table
PNeCat.add_column(Column(PNeCoords.galactic.represent_as(CartesianRepresentation).x.to(u.pc),name='x_gal_coord',meta={'ucd': 'pos.cartesian.x'}),0
PNeCat.add_column(Column(PNeCoords.galactic.represent_as(CartesianRepresentation).y.to(u.pc),name='y_gal_coord',meta={'ucd': 'pos.cartesian.y'}),1)
PNeCat.add_column(Column(PNeCoords.galactic.represent_as(CartesianRepresentation).z.to(u.pc),name='z_gal_coord',meta={'ucd': 'pos.cartesian.z'}),2)
#Remove name Column since we can parse labels ***TO DO*** capture this in a .label file
PNeCat.remove_column("Name")
#Write .speck file for DigitalSky and Uniview
writeFile = open('data/PNe.speck','a')
#Define and write the header
headerText="\
#Planetary Nebulae From  (Stanghellini+, 2008)\n\
dataVar 0 ra\n\
dataVar 1 Dec\n\
dataVar 2 tau\n\
dataVar 3 Rad\n\
dataVar 4 F5GHz\n\
dataVar 5 Dist\n"
writeFile.write(headerText)
#Write the data
PNeCat.write(writeFile,format='ascii.no_header')
writeFile.close()


# ## Globular Clusters

# Properties of Galactic Globular Clusters<a hhref="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?2014MNRAS.441.1105F&db_key=AST&nosetcookie=1">((Francis+, 2014))</a><br>
# VizieR Catalog: <a href="http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=VII/271/catalog&-out.max=50&-out.form=HTML%20Table&-out.add=_r&-out.add=_RAJ,_DEJ&-sort=_r&-oc.form=sexa">VII/271/catalog</a><br>
# 157 Objects

# ###### Get Globular Cluster Catalog

# In[8]:


Cats = v.get_catalogs('VII/271/catalog')
GCCat=Cats[0]
del GCCat['RAJ2000','DEJ2000','SimbadName','Name','mu0','r_Dist','DistH','X','Y','Z','Rgc','E_B-V_','RV','r_RV','q_RV','r__Fe_H_']
GCCat.rename_column('__Fe_H_', 'Fe/H')
GCCat.rename_column('_RAJ2000', 'RA')
GCCat.rename_column('_DEJ2000', 'dec')
#Create Column with dist in pc, not kpc (for WWT)
distCol=Column(1000*GCCat['Dist'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
GCCat.add_column(distCol)
del GCCat['Dist']
GCCoords=SkyCoord(GCCat['RA'],GCCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(GCCat['distance'],u.pc),frame='icrs')


# ###### Print Globular Cluster Catalog

# In[9]:


GCCat #just first 10


# ###### Plot Globular Cluster Catalog

# In[10]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Globular Clusters",fontsize=18)
ax.scatter(GCCoords.galactic.l.wrap_at(180.*u.degree).radian,GCCoords.galactic.b.radian,c='red')


# ###### Send Globular Clusters to WWT

# In[11]:


#Set up WWT layer
GC_layer = wwt.new_layer("Galactic Datasets", "Globular Clusters", GCCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"4",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:255:255:255",\
              "TimeSeries":"False"}
GC_layer.set_properties(props_dict)
#Send data to WWT client
GC_layer.update(data=GCCat, purge_all=True, no_purge=False, show=True)


# ## Open Clusters

# Optically visible open clusters and Candidates<a href="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?2002A%26A...389..871D&db_key=AST&nosetcookie=1"> (Dias+ 2002-2013)</a><br>
# VizieR Catalog: <a href="http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=B/ocl/clusters&-out.max=50&-out.form=HTML%20Table&-out.add=_r&-out.add=_RAJ,_DEJ&-sort=_r&-oc.form=sexa">B/ocl/clusters</a><br>
# 2174 Objects

# ###### Get Open Cluster Catalog

# In[8]:


Cats = v.get_catalogs('B/ocl/clusters')
ocCat=Cats[0]
ocCat.keep_columns(["Cluster","_RAJ2000","_DEJ2000","Dist","Diam","Age"])
ocCat.rename_column('_RAJ2000', 'RA')
ocCat.rename_column('_DEJ2000', 'dec')
#Create Column with dist in pc (not kpc for WWT)
distCol=Column(ocCat['Dist'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
ocCat.add_column(distCol)
del ocCat['Dist']
ocCoords=SkyCoord(ocCat['RA'],ocCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(ocCat['distance'],u.pc),frame='icrs')


# ###### Print Open Cluster Calalog

# In[41]:


ocCat[:10] #just first 10


# ###### Plot Open Cluster Catalog

# In[33]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Open Clusters",fontsize=18)
ax.scatter(ocCoords.galactic.l.wrap_at(180.*u.degree).radian,ocCoords.galactic.b.radian,c='green')


# ###### Send Open Cluster catalog to WWT

# In[13]:


#Set up WWT layer
OC_layer = wwt.new_layer("Galactic Datasets", "Open Clusters", ocCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"5",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:255:255:0",\
              "TimeSeries":"False"}
OC_layer.set_properties(props_dict)
#Send data to WWT client
OC_layer.update(data=ocCat, purge_all=True, no_purge=False, show=True)


# ## Supernova Remnants

# A revised Galactic supernova remnant catalogue<a href="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?2014yCat.7272....0G&db_key=AST&nosetcookie=1http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?2009BASI...37...45G&db_key=AST&nosetcookie=1"> (Green 2014)</a><br>
# VizieR Catalog: <a href="http://vizier.ast.cam.ac.uk/viz-bin/VizieR?-source=VII/272">VII/272/snrs</a><br>
# 294 Objects<br><br>
# This catalog does not provide distances, so we will only provide code to visualize in WWT.

# ###### Get Supernova Remnant Catalog

# In[10]:


Cats = v.get_catalogs('VII/272/snrs')
snrCat=Cats[0]
snrCat.keep_columns(["SNR","_RAJ2000","_DEJ2000","MajDiam","S_1GHz_"])
snrCat.rename_column('S_1GHz_', 'S_1GHz')
snrCat.rename_column('_RAJ2000', 'RA')
snrCat.rename_column('_DEJ2000', 'dec')
snrCoords=SkyCoord(snrCat['RA'],snrCat['dec'],unit=(u.degree, u.degree),frame='icrs')


# ###### Print Supernova Remnant Catalog

# In[49]:


print len(snrCat)
snrCat[:10] #Just first 10


# ###### Plot Supernova Remnant Catalog

# In[34]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Supernova Remnants",fontsize=18)
ax.scatter(snrCoords.galactic.l.wrap_at(180.*u.degree).radian,snrCoords.galactic.b.radian,c='magenta')


# ###### Send Supernova Remnants to WWT

# In[15]:


#Set up WWT layer
SNR_layer = wwt.new_layer("Galactic Datasets", "Supernova Remnants", snrCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:255:0:255",\
              "TimeSeries":"False"}
SNR_layer.set_properties(props_dict)
#Send data to WWT client
SNR_layer.update(data=snrCat, purge_all=True, no_purge=False, show=True)


# ## Pulsars

# Catalog of Pulsars<a href="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?1993ApJS...88..529T&db_key=AST&nosetcookie=1">(Taylor+ 1995)</a><br>
# VizieR Catalog: <a href="http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=VII/189/table1&-out.max=50&-out.form=HTML%20Table&-out.add=_r&-out.add=_RAJ,_DEJ&-sort=_r&-oc.form=sexa">VII/189/table1</a><br>
# 706 Objects

# ###### Get Pulsar Catalog

# In[23]:


Cats = v.get_catalogs('VII/189/table1')
pulsarCat=Cats[0]
pulsarCat.keep_columns(["Jname","RAJ2000","DEJ2000","dist","P","Age"])
#Create Column with dist in pc (not kpc) for WWT
distCol=Column(1000*pulsarCat['dist'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
pulsarCat.add_column(distCol)
del pulsarCat['dist']
pulsarCat.rename_column('RAJ2000', 'RA')
pulsarCat.rename_column('DEJ2000', 'dec')
pulsarCoords=SkyCoord(pulsarCat['RA'],pulsarCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(pulsarCat['distance'],u.pc),frame='icrs')


# ###### Print Pulsar Catalog

# In[18]:


print len(pulsarCat)
pulsarCat[:10] #just first 10


# ###### Plot Pulsar Catalog

# In[35]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Pulsars",fontsize=18)
ax.scatter(pulsarCoords.galactic.l.wrap_at(180.*u.degree).radian,pulsarCoords.galactic.b.radian,c='orange')


# ###### Send Pulsar Catalog to WWT

# In[86]:


#Set up WWT layer
pulsar_layer = wwt.new_layer("Galactic Datasets", "Pulsars", pulsarCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"5",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:0:255:255",\
              "TimeSeries":"False"}
pulsar_layer.set_properties(props_dict)
#Send data to WWT client
pulsar_layer.update(data=pulsarCat, purge_all=True, no_purge=False, show=True)


# ## Magnetars

# The McGill magnetar catalog <a href="http://cdsbib.u-strasbg.fr/cgi-bin/cdsbib?2014ApJS..212....6O">(Olausen+, 2014)</a><br>
# VizieR Catalog: <a href="http://cdsbib.u-strasbg.fr/cgi-bin/cdsbib?2014ApJS..212....6O">J/ApJS/212/6/table9</a><br>
# 26 Objects

# ###### Get Magnetar Catalog

# In[14]:


Cats = v.get_catalogs('J/ApJS/212/6/table9')
magCat=Cats[0]
magCat.keep_columns(["Name","_RAJ2000","_DEJ2000","Dist","B","Per"])
magCat.rename_column('_RAJ2000', 'RA')
magCat.rename_column('_DEJ2000', 'dec')
#Create Column with dist in pc, not kpc (for WWT)
distCol=Column(1000*magCat['Dist'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
magCat.add_column(distCol)
del magCat['Dist']
magCoords=SkyCoord(magCat['RA'],magCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(magCat['distance'],u.pc),frame='icrs')


# ###### Print Magnetar Catalog

# In[62]:


print len(magCat)
magCat[:10] #just first 10


# ###### Plot Magnetar Catalog

# In[37]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("Magnetars",fontsize=18)
ax.scatter(magCoords.galactic.l.wrap_at(180.*u.degree).radian,magCoords.galactic.b.radian,c='yellow')


# ###### Send Magnetar Data to WWT

# In[15]:


#Set up WWT layer
Magnetar_layer = wwt.new_layer("Galactic Datasets", "Supernova Remnants", magCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"5",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:127:0:255",\
              "TimeSeries":"False"}
SNR_layer.set_properties(props_dict)
#Send data to WWT client
SNR_layer.update(data=amgCat, purge_all=True, no_purge=False, show=True)


# ## OB Associations

# New list of OB associations<a href="http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?1995PAZh...21...13M&db_key=AST&nosetcookie=1">(Melnik+, 1995) </a><br>
# VizieR Catalog: <a href="http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/PAZh/21/13/table&-out.max=50&-out.form=HTML%20Table&-out.add=_r&-out.add=_RAJ,_DEJ&-sort=_r&-oc.form=sexa">J/PAZh/21/13/table</a><br>
# 88 Objects

# ###### Get OB Association Catalog

# In[16]:


Cats = v.get_catalogs('J/PAZh/21/13/table')
OBCat=Cats[0]
OBCat.keep_columns(["AssME","_RA.icrs","_DE.icrs","D","Ntot"])
#Create Column with dist in pc (not kpc)
distCol=Column(1000*OBCat['D'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
OBCat.add_column(distCol)
del OBCat['D']
OBCat.rename_column('_RA.icrs', 'RA')
OBCat.rename_column('_DE.icrs', 'dec')
OBCoords=SkyCoord(OBCat['RA'],OBCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(OBCat['distance'],u.pc),frame='icrs')


# ###### Print OB Associations Catalog

# In[66]:


print len(OBCat)
OBCat[:10] #just first 10


# ###### Plot OB Associations

# In[38]:


fig = plt.figure (figsize=(13,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
plt.title("OB Associations",fontsize=18)
ax.scatter(OBCoords.galactic.l.wrap_at(180.*u.degree).radian,OBCoords.galactic.b.radian,c='cyan')


# In[87]:


#Set up WWT layer
OB_layer = wwt.new_layer("Galactic Datasets", "OB Associations", OBCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
              "MarkerScale":"Screen",\
              "PointScaleType":"Constant",\
              "ScaleFactor":"8",\
              "ShowFarSide":"True",\
              "AltColumn":"4",\
              "AltType":"Distance",\
              "AltUnit":"Parsecs",\
              "RaUnits":"Degrees",\
              "PlotType":"Circle",\
              "ColorValue":"ARGBColor:255:64:64:255",\
              "TimeSeries":"False"}
OB_layer.set_properties(props_dict)
#Send data to WWT client
OB_layer.update(data=OBCat, purge_all=True, no_purge=False, show=True)


# In[ ]: