Data-analysis and visualization examples

This notebook does not directly demonstrate voeventdb (see the tutorials!) but provides examples of what’s possible when data on astronomical transients is readily available:

Prelude - fetch 10 most-recent GRBs

from __future__ import print_function
import voeventdb.remote as vr
import voeventdb.remote.apiv1 as apiv1
from voeventdb.remote.apiv1 import FilterKeys
from voeventdb.remote.helpers import Synopsis
from datetime import datetime, timedelta
import pytz
import logging
logging.basicConfig(level=logging.DEBUG)

filters = { FilterKeys.ivorn_contains: 'BAT_GRB',
            FilterKeys.role: 'observation'}

recent_swift_grb_ivorns = apiv1.list_ivorn(
    filters,
    order=apiv1.OrderValues.author_datetime_desc,
    n_max=10,
)
recent_swift_grbs = [Synopsis(apiv1.packet_synopsis(i))
                     for i in recent_swift_grb_ivorns]

Plotting a timeline

In the last tutorial we retrieved and printed the timestamps of the most recent GRBs - but deciphering textual timestamps isn’t much fun. Let’s plot a timeline instead:

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

grb_dtimes=[s.author_datetime for s in recent_swift_grbs]

fonts = {'weight' : 'bold',
        'size'   : 14}

mpl.rc('font', **fonts)

now = pytz.UTC.localize((datetime.utcnow()))
week_markers = [now - timedelta(days=7)*w for w in range(0,5)]


markersize=95

plt.scatter(grb_dtimes, [1]*len(grb_dtimes),  marker='*',
            s=markersize, label='GRB')
plt.scatter(now, 1, marker='o', s=markersize, c='r')
for d in grb_dtimes:
    plt.axvline(d, ymax=0.5, ls='-')

first_label = True
for w in week_markers:
    plt.axvline(w, ymax=0.5, ls='--',c='r',
                label=('Week marker' if first_label else None))
    first_label=False

plt.xlim(min(grb_dtimes)-timedelta(days=2),
         max(max(grb_dtimes),now)+timedelta(days=2))
plt.gcf().autofmt_xdate()
ax = plt.gca()
ax.yaxis.set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.get_yaxis().set_ticklabels([])
plt.legend(loc='best')
plt.gcf().set_size_inches(12,4)
plt.gcf().suptitle("Recent Swift BAT alerts, as of {}".format(now),
                   fontsize=22)

Creating a sky-coordinate scatterplot

Let’s see where those GRBs lie on-sky (with credit to http://www.astropy.org/astropy-tutorials/plot-catalog.html):

from astropy.coordinates import Angle
import astropy.units as u
grb_ra_coords = Angle([grb.coords[0].ra for grb in recent_swift_grbs])
grb_dec_coords = Angle([grb.coords[0].dec for grb in recent_swift_grbs])
grb_ra_coords = grb_ra_coords.wrap_at(180*u.degree)

print(grb_ra_coords.deg)
print(grb_dec_coords.deg)

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection="mollweide")
ax.scatter(grb_ra_coords.radian, grb_dec_coords.radian,
           marker='*', s=13**2, label='GRB')
ax.grid(True)
plt.legend()
plt.gcf().suptitle('Locations of 10 most recent Swift BAT GRB alerts',
                   fontsize=24)

Planning an observation with astropy

Let’s suppose we want to observe a recently discovered GRB from Cambridge, UK. We have the location of the event, already converted to astropy co-ordinates. This means we can use astropy’s Alt-Az conversions to calculate the altitude of the target at any given time.

(This section borrows heavily from http://docs.astropy.org/en/stable/coordinates/observing-example.html, consult the original for more detail and background)

#Grab the latest IERS time-data first:
from astropy.utils.data import download_file
from astropy.utils import iers
iers.IERS.iers_table = iers.IERS_A.open(
    download_file(iers.IERS_A_URL, cache=True))

import numpy as np
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation, AltAz

Let’s pick a GRB which will actually be observable from the UK, i.e. one in the Northern hemisphere:

filters[FilterKeys.dec_greater_than]=0.0
filters

northern_grb_list = apiv1.list_ivorn(
    filters,
    order=apiv1.OrderValues.author_datetime_desc,
    n_max=1,
    )

grb_target = Synopsis(apiv1.packet_synopsis(northern_grb_list[0]))

print("Sample target GRB:")
print(grb_target)

grb_target_posn = grb_target.coords[0]
cambridge = EarthLocation(lat=52.205*u.deg, lon=0.118*u.deg, height=6*u.m)
utcoffset = -0*u.hour  # GMT
time = Time.now() - utcoffset
grb_altaz = grb_target_posn.transform_to(AltAz(obstime=time,location=cambridge))
print("GRB's current altitude = {0.alt:.1f}".format(grb_altaz)  )

from astropy.coordinates import get_sun
now = Time.now()
delta_24 = np.linspace(0, 24, 100)*u.hour
times = now + delta_24

altaz_frame = AltAz(obstime=times, location=cambridge)

grb_altaz_array = grb_altaz.transform_to(altaz_frame)
sun_altaz_array = get_sun(times).transform_to(altaz_frame)

fig = plt.figure(figsize=(12,8))
plt.plot(delta_24, grb_altaz_array.alt, lw=10, label='GRB')
plt.plot(delta_24, sun_altaz_array.alt, lw=10, c='y', label='Sun')
# plt.fill_between(delta_midnight, 0, 90, sun_altaz_array.alt < -0*u.deg, color='0.5', zorder=0)
# plt.fill_between(delta_midnight, 0, 90, sun_altaz_array.alt < -18*u.deg, color='k', zorder=0)

plt.axhline(0, ls='--', lw=8, c='k',
           label='Horizon',)
plt.xlabel('Hours from {}'.format(now))
plt.ylabel('Source altitude')
plt.legend(loc='best')
plt.gcf().suptitle("GRB target altitude over next 24 hours", fontsize=24)