Basic Observability Calculations¶

Now we are going to teach our notice handler how to determine whether a gravitational-wave event is observable. We are going to use the astropy.coordinates module. (See also the Astropy example on observation planning in Python.) First, we will need to import a few extra Python modules:

import astropy.coordinates
import astropy.time
import astropy.units as u
import astropy_healpix as ah

The LIGO/Virgo/KAGRA probability sky maps are always in equatorial coordinates. Once we have looked up the coordinates of the HEALPix pixels, we will use Astropy to transform those coordinates to a horizontal (altitude–azimuth) frame for a particular site on the Earth at a particular time. Then we can quickly determine which pixels are visible from that site at that time, and integrate (sum) the probability contained in those pixels.

Note

You may want to do something more sophisticated like determine how much of the probability is visible for at least a certain length of time. This example will illustrate one key function of HEALPix (looking up coordinates of the grid with ah.healpix_to_lonlat) and some of the key positional astronomy functions with Astropy. For more advanced functionality, we recommend the astroplan package.

def prob_observable(skymap, time):
    """
    Determine the integrated probability contained in a gravitational-wave
    sky map that is observable from a particular ground-based site at a
    particular time.

    Bonus: make a plot of probability versus UTC time!
    """

    # Determine the sky area and probability contained in each pixel of the
    # multi-resolution HEALPix sky map.
    level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
    nside = ah.level_to_nside(level)
    pixel_area = ah.nside_to_pixel_area(nside)
    prob = skymap['PROBDENSITY'] * pixel_area

    # Look up the (celestial) coordinates of the center of each pixel.
    ra, dec = ah.healpix_to_lonlat(ipix, nside, order='nested')
    radecs = astropy.coordinates.SkyCoord(ra=ra, dec=dec)

    # Geodetic coordinates of observatory (example here: Mount Wilson)
    observatory = astropy.coordinates.EarthLocation(
        lat=34.2247*u.deg, lon=-118.0572*u.deg, height=1742*u.m)

    # Alt/az reference frame at observatory, at the given time
    frame = astropy.coordinates.AltAz(obstime=time, location=observatory)

    # Transform pixel coordinates to alt/az at the observatory
    altaz = radecs.transform_to(frame)

    # Where is the sun, at the given time?
    sun_altaz = astropy.coordinates.get_sun(time).transform_to(frame)

    # How likely is it that the (true, unknown) location of the source
    # is within the area that is visible? Demand that the sun is at
    # least 18 degrees below the horizon and that the airmass
    # (secant of zenith angle approximation) is at most 2.5.
    return prob[(sun_altaz.alt <= -18*u.deg) & (altaz.secz <= 2.5)].sum()

We will build on the JSON notice handler from the GCN sample code section. We need a few extra imports to decode the sky map and call our new function:

from base64 import b64decode
from io import BytesIO
import json

from astropy.table import Table
from gcn_kafka import Consumer

Now we update the handler to decode the sky map and call prob_observable:

def parse_notice(record):
    record = json.loads(record)

    # Only respond to mock events. Real events have GraceDB IDs like
    # S1234567, mock events have GraceDB IDs like M1234567.
    # NOTE NOTE NOTE replace the conditional below with this commented out
    # conditional to only parse real events.
    # if record['superevent_id'][0] != 'S':
    #    return
    if record['superevent_id'][0] != 'M':
        return

    if record['alert_type'] == 'RETRACTION':
        print(record['superevent_id'], 'was retracted')
        return

    # Respond only to 'CBC' events. Change 'CBC' to 'Burst' to respond to
    # only unmodeled burst events.
    if record['event']['group'] != 'CBC':
        return

    # Read the multi-resolution HEALPix sky map.
    skymap_str = record.get('event', {}).get('skymap')
    if skymap_str:
        skymap_bytes = b64decode(skymap_str)
        skymap = Table.read(BytesIO(skymap_bytes))

        prob = prob_observable(skymap, astropy.time.Time.now())
        print('Source has a {:d}% chance of being observable now'.format(
            int(round(100 * prob))))
        if prob > 0.5:
            pass # FIXME: perform some action

Let’s run the new handler now…

consumer = Consumer(client_id='fill me in', client_secret='fill me in')
consumer.subscribe(['igwn.gwalert'])

while True:
    for message in consumer.consume():
        parse_notice(message.value())

When you run this script, each time you receive a sample LIGO/Virgo/KAGRA notice, it will print something like the following (note that probability will change as a function of time):

Source has a 76% chance of being observable now