from astropy.table import Table,Column
from astropy.time import Time
from astropy import units
from astropy.coordinates import SkyCoord
from astroquery.vizier import Vizier
from astropy import units as u
#Create Vizier object, turn off default row limit
v = Vizier()
v.ROW_LIMIT = -1
from pywwt.mods import *
#Connect to WWT
wwt = WWTClient(host="127.0.0.1") #Can pass a IP address here if WWT is running on a remote machine
For out data catalog we'll choose The second Fermi/GBM GRB catalog (4yr) (von Kienlin+, 2014) Vizier catalog: J/ApJS/211/13/GBM Which contains Fermi events from July 2007 to July 2012
Cats = v.get_catalogs('J/ApJS/211/13/GBM')
grbCat=Cats[0]
grbCat.keep_columns(["GRB","RAJ2000","DEJ2000","Time","ObsTime","Fl.w","Fl.n"])
grbCat.rename_column('RAJ2000', 'RA')
grbCat.rename_column('DEJ2000', 'dec')
C:\Anaconda\lib\site-packages\astroquery\vizier\core.py:556: UserWarning: VOTABLE parsing raised exception: warnings.warn("VOTABLE parsing raised exception: {0}".format(ex))
#Set up matplotlib
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn
#Plot the Catalog
coordsCol=SkyCoord(grbCat['RA'],grbCat['dec'],unit=(u.degree, u.degree), frame= 'icrs')
fig = plt.figure (figsize=(12,6))
ax = fig.add_subplot(111,projection="mollweide")
ax.grid(True)
ax.get_xaxis().tick_bottom()
ax.scatter(coordsCol.galactic.l.wrap_at(180.*units.degree).radian,\
coordsCol.galactic.b.radian,s=6,lw=0,c='purple')
<matplotlib.collections.PathCollection at 0x18847eb8>
coordsCol=SkyCoord(grbCat['RA'],grbCat['dec'],unit=(u.degree, u.degree), frame= 'icrs')
del coordsCol
Extracting the event time from this table is tricky. The time of day is in the ObsTime column, but the date is embedded in the GRB name. We'll extract the date from the GRB name and combine that with ObsTime to make a astropy time object.
timeList=[]
grbList=grbCat['GRB']
for i in range(len(grbList)):
timeString= grbList[i][2:4]+'/'+grbList[i][4:6]+'/'+'20'+grbList[i][0:2]+' '+grbCat['ObsTime'][i]
timeList.append(timeString)
WWT contains its own time format, unfortunately one that astropy cannot write, so we'll have to create our own custom string.
grbList=grbCat['GRB']
timeList=[]
for i in range(len(grbList)):
timeString= grbList[i][2:4]+'/'+grbList[i][4:6]+'/'+'20'+grbList[i][0:2]+' '+grbCat['ObsTime'][i]
timeList.append(timeString)
grbCat.add_column(Column(timeList,name='TimeAndDate'))
grbCat
GRB | ObsTime | RA | dec | Time | Fl.w | Fl.n | TimeAndDate |
---|---|---|---|---|---|---|---|
"h:m:s" | deg | deg | ms | mJ / m2 | mJ / m2 | ||
080714B | 02:04:12.0534 | 41.9 | 8.5 | 512 | 6.8e-07 | 3.5e-07 | 07/14/2008 02:04:12.0534 |
080714C | 10:12:01.8376 | 187.5 | -74.0 | 4096 | 1.8e-06 | 9.8e-07 | 07/14/2008 10:12:01.8376 |
080714A | 17:52:54.0234 | 188.1 | -60.2 | 1024 | 6.3e-06 | 3.3e-06 | 07/14/2008 17:52:54.0234 |
080715A | 22:48:40.1634 | 214.7 | 9.9 | 256 | 5e-06 | 2.5e-06 | 07/15/2008 22:48:40.1634 |
080717A | 13:02:35.2207 | 147.3 | -70.0 | 4096 | 4.5e-06 | 2.4e-06 | 07/17/2008 13:02:35.2207 |
080719A | 12:41:40.9578 | 153.2 | -61.3 | 4096 | 7.7e-07 | 3.9e-07 | 07/19/2008 12:41:40.9578 |
080720A | 07:35:35.5476 | 98.5 | -43.9 | 8192 | -- | -- | 07/20/2008 07:35:35.5476 |
080723B | 13:22:21.3751 | 176.8 | -60.2 | 256 | 7.2e-05 | 3.9e-05 | 07/23/2008 13:22:21.3751 |
080723C | 21:55:23.0583 | 113.3 | -19.7 | 64 | 1.3e-07 | 7.5e-08 | 07/23/2008 21:55:23.0583 |
080723D | 23:37:42.7083 | 105.3 | 71.1 | 512 | 3.1e-05 | 1.6e-05 | 07/23/2008 23:37:42.7083 |
080724A | 09:37:40.6034 | 358.3 | 32.9 | 256 | 1.6e-05 | 8.7e-06 | 07/24/2008 09:37:40.6034 |
... | ... | ... | ... | ... | ... | ... | ... |
120701B | 15:41:48.3152 | 182.7 | -45.7 | 128 | 8.4e-08 | 4.9e-08 | 07/01/2012 15:41:48.3152 |
120702A | 21:23:19.1712 | 227.8 | 36.8 | 512 | 1.6e-06 | 1e-06 | 07/02/2012 21:23:19.1712 |
120703B | 10:01:11.6882 | 69.5 | 34.7 | 1024 | 1.1e-05 | 5.5e-06 | 07/03/2012 10:01:11.6882 |
120703C | 11:56:56.8702 | 210.5 | 46.3 | 1024 | 2.6e-06 | 1.5e-06 | 07/03/2012 11:56:56.8702 |
120703A | 17:25:17.0323 | 339.4 | -29.7 | 2048 | 8.3e-06 | 4.3e-06 | 07/03/2012 17:25:17.0323 |
120707A | 19:12:17.4295 | 291.1 | -34.4 | 4096 | 9.4e-05 | 5.2e-05 | 07/07/2012 19:12:17.4295 |
120709A | 21:11:40.3666 | 318.4 | -50.1 | 64 | 1.4e-05 | 6.3e-06 | 07/09/2012 21:11:40.3666 |
120710A | 02:23:17.0507 | 120.4 | -31.1 | 256 | 5.3e-06 | 2.7e-06 | 07/10/2012 02:23:17.0507 |
120711A | 02:44:53.2943 | 94.7 | -71.0 | 256 | 0.00019 | 6.6e-05 | 07/11/2012 02:44:53.2943 |
120711C | 10:42:54.5709 | 127.9 | -31.8 | 1024 | 1.9e-06 | 1e-06 | 07/11/2012 10:42:54.5709 |
#Set up WWT layer
grb_layer = wwt.new_layer("Dynamic Universe", "Gamma Ray Bursts", grbCat.colnames)
#Set visualization parameters in WWT
props_dict = {"CoordinatesType":"Spherical",\
"MarkerScale":"Screen",\
"PointScaleType":"Constant",\
"ScaleFactor":"64",\
"ShowFarSide":"True",\
"RaUnits":"Degrees",\
"PlotType":"Gaussian",\
"ColorValue":"ARGBColor:255:255:255:255",\
"TimeSeries":"False"}
grb_layer.set_properties(props_dict)
#Send data to WWT client
grb_layer.update(data=grbCat, purge_all=True, no_purge=False, show=True)
Now inside WWT we can choose how we visualize the data, we can show all the data at once or playback the events as they happen watching the GRB’s go off like popcorn across the sky.