# WARNING: This document is obsolete.

## Please refer to the LIGO-Virgo Public Alerts User Guide, https://emfollow.docs.ligo.org/userguide/.

# LIGO-Virgo EM Follow-Up Tutorial

by Leo P. Singer (NASA/GSFC) leo.singer@ligo.org

This document is LIGO-G1500442-v11.

## Abstract

This document explains how to receive, filter, and process gravitational-wave (GW) detection candidate alerts from Advanced LIGO and Virgo. We provide sample code in Python and document alternatives for users of other programming environments. You can download this document and run the code samples in IPython Notebook.

## Table of Contents

Introduction

1. Sign up for GCN/TAN network
2. Sign up for a GraceDb robot password
3. Install some dependencies
4. Write GCN handler script
5. Working with probability sky maps
6. Basic observability calculations with Astropy
7. Submitting observation coordinates to GraceDB

Appendix: Full example code

## Introduction

The LIGO-Virgo data are analyzed in near real-time to search for GW transients due to compact binary coalescence (CBC) events or unmodeled "burst" sources. For each detection candidate, a series of alerts are produced and distributed via the Gamma-ray Coordinates Network/Transient Astronomy Network (GCN/TAN).

The GCN/TAN system may already be familiar to some, as it is has been in use since the early 1990s to transmit times and coordinates of gamma-ray bursts (GRBs) detected by gamma-ray space missions to ground-based observers. However, there are three peculiariaties of GW observations that veteran GCN users should be aware of:

1. Unlike GRBs, GW sources cannot (always) be localized to a unique sky location; GW position reconstructions are (sometimes) multimodal and non-Gaussian. As a consequence, LIGO-Virgo GCN notices do not contain a RA, Dec, or error radius; instead they contain a pointer URL to a FITS file containing a probability sky map in the HEALPix all-sky projection.
2. GCN notices will produce LIGO-Virgo candidates that span a range of significances. The significance of a candidate is reported as the false alarm rate (FAR). The FAR of a detection candidate measures approximately how frequently an event of similar strenght occurs due to chance noise fluctuations or instrumental glitches. A FAR of 1/century (∼ 3 × 10−10 Hz) is generally regarded as sufficient for a "confident" detection. However, alerts will produced for all events with FAR≥1/month (∼ 4 × 10−7 Hz), such that there should be at least one alert per month on average.
3. The first few GW alerts will be distributed privately to LIGO MOU partners, over a purpose-built GCN connection.

## 1. Sign up for GCN/TAN network

The first step is to sign up for the GCN network. Signing up for LIGO/Virgo GCN notices is slightly different from the standard signup process. However, if you are already receiving GCN notices (for example, for GRB follow-up), then you can reuse your existing GCN configuration and add the LIGO/Virgo notices.

There are several distribution methods for GCN notices. For the purpose of this tutorial, we will focus on VOEvent over VOEvent Transport Protocol, which is among the more convenient methods for autonomous operations. However, you can use any other distribution method of your choice.

To receive VOEvents, you will need a computer with a static IP address and you will need to register the IP address from which you will connect to the GCN network. Do the following two steps to submit a new GCN site configuration:

1. Go to http://gcn.gsfc.nasa.gov/lvc_config_builder.html, fill out the form, and then click the "Send LVC Notice request to GCN Admin" button.
2. If you are registering as a new GCN recipient (as opposed to adding LIGO/Virgo notices to an existing site configuration), then click the "New Site or New User" button and fill out the resulting form. Select the `VOEVENT2.0` option for the `Distribution Method` field.

## 2. Sign up for a GraceDb robot password

You will need a robot password to download files that are linked from the GCN notices. This is an automatically generated password that you can use to download files from GraceDb in a script. It is different from your log in for the GraceDb web site, but you use the GraceDb web site to manage your robot password. Here is how to obtain a robot password.

1. Got to https://gracedb.ligo.org/ and log in.
2. Click the Options tab (see screenshot below).
3. Click the Manage Password link (see screenshot below).
4. Click Get me a password! (See screenshot below.) Warning: if you had an old robot password, this will reset it and generate a new one.
5. Record the robot user name (which is the same as your user name for the GraceDb web site) and the generated robot password (see screenshot below).
6. Optional, but recommended: For convenience, we suggset that you add the following line to your ~/.netrc file:

`machine gracedb.ligo.org login` albert.einstein@LIGO.ORG `password` ABCDEabcde0123456789

replacing albert.einstein@LIGO.ORG with the robot user name and ABCDEabcde0123456789 with the robot password. This will make it so that you can call `curl` or other HTTP downloader programs without explicitly providing the user name and password. To test your user name and password, you can run the following command:

    curl --netrc https://gracedb.ligo.org/apibasic/

## 3. Install some dependencies

You will need to install a few third-party Python packages to run the example code in this tutorial. These include:

- pygcn for connecting to GCN (alternatives: comet)
- requests for easy HTTP downloads in Python (many alternatives in Python standard library)
- healpy, for decoding HEALPix images (alternatives: DS9, Aladin, HEALPix bindings for C/C++/Fortran/Java/IDL, reprojection with the Python reproject package)
- astropy version 1.0.1 or newer (optional, for computing observability windows, etc.)
- numpy and matplotlib, popular math and plotting packages for Python

If you are on a Mac and use the MacPorts package manager, you can install all of the above with the following command:

    $ sudo port install py27-gcn py27-lxml py27-requests py27-healpy py27-astropy

Otherwise, the fastest way to install the dependencies is with `pip`, a package manager that comes with most Python distributions. To install these packages with `pip`, run the following command:

    $ pip install pygcn lxml requests healpy astropy

## 4. Write GCN handler script

Now we'll write a GCN handler script. First, some imports...

    # Python standard library imports
    import tempfile
    import shutil
    import sys
    import glob

    # Third-party imports
    import gcn
    import gcn.handlers
    import gcn.notice_types
    import requests
    import healpy as hp
    import numpy as np

IPython Notebook only: If you are following along and running commands in the IPython Notebook, you should also run the following command so that plots are rendered inline inside the notebook:

    %matplotlib inline

Next, we'll write a function that we want to get called every time a GCN notice is received. We will use the `@gcn.handlers.include_notice_types` function decorator to specify that we only want to process certain notice types. There are three notice types:

1. `LVC_PRELIMINARY`: Provides the time, significance, and basic parameters about a GW detection candidate. No localization information. Sent with a latency of a minute or so.
2. `LVC_INITIAL`: A rapid sky localization is available. Sent with a latency of a few minutes.
3. `LVC_UDPATE`: A refined sky localization is availaable. Sent with a latency of hours or more.

In the following example, we will process only the last two types (`LVC_INITIAL` and `LVC_UPDATE`), which contain links to sky map FITS files. The following handler function will parse out the URL of the FITS file, download it, and extract the probability sky map.

Events come in two very general flavors: 'CBC' or compact binary coalescence candidates detected by matched filtering, and generic 'Burst' candidates detected by model-independent methods. Most users will want to receive only 'CBC' or only 'Burst' events. In this example code, we are going to keep only 'CBC' events.

Very important: Note that mock or 'test' observations are denoted by the `role="test"` VOEvent attribute. Alerts resulting from real LIGO/Virgo science data will always have `role="observation"`. The sample code below will respond only to 'test' events. When preparing for actual observations, you must remember to switch to 'observation' events.

    def get_skymap(root):
        """
        Look up URL of sky map in VOEvent XML document,
        download sky map, and parse FITS file.
        """
        # Read out URL of sky map.
        # This will be something like
        # https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gz
        skymap_url = root.find(
            "./What/Param[@name='SKYMAP_URL_FITS_BASIC']").attrib['value']

        # Send HTTP request for sky map
        response = requests.get(skymap_url, stream=True)

        # Uncomment to save VOEvent payload to file
        # open('example.xml', 'w').write(payload)

        # Raise an exception unless the download succeeded (HTTP 200 OK)
        response.raise_for_status()

        # Create a temporary file to store the downloaded FITS file
        with tempfile.NamedTemporaryFile() as tmpfile:
            # Save the FITS file to the temporary file
            shutil.copyfileobj(response.raw, tmpfile)
            tmpfile.flush()

            # Uncomment to save FITS payload to file
            # shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))

            # Read HEALPix data from the temporary file
            skymap, header = hp.read_map(tmpfile.name, h=True, verbose=False)
            header = dict(header)

        # Done!
        return skymap, header

    # Function to call every time a GCN is received.
    # Run only for notices of type LVC_INITIAL or LVC_UPDATE.
    @gcn.handlers.include_notice_types(
        gcn.notice_types.LVC_INITIAL,
        gcn.notice_types.LVC_UPDATE)
    def process_gcn(payload, root):
        # Print the alert
        print('Got VOEvent:')
        print(payload)

        # Respond only to 'test' events.
        # VERY IMPORTANT! Replce with the following line of code
        # to respond to only real 'observation' events.
        # if root.attrib['role'] != 'observation': return
        if root.attrib['role'] != 'test': return

        # Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only
        # unmodeled burst events.
        if root.find("./What/Param[@name='Group']").attrib['value'] != 'CBC': return

        # Read out integer notice type (note: not doing anythin with this right now)
        notice_type = int(root.find("./What/Param[@name='Packet_Type']").attrib['value'])

        # Read sky map
        skymap, header = get_skymap(root)

Finally, we will listen for GCNs. You need to tell the `gcn.listen` function on what port you want to connect to the GCN network; this will be port 8096. You also need to tell `gcn.listen` which function to call whenever it receives an GCN; this will be the `process_gcn` function that we just defined.

When you run the following code snippet, the `gcn.listen` function will continue until you interrupt the program (by pressing the stop button in IPython Notebook, typing Control-C in the terminal, or sending a kill signal to your Python script).

Note: `gcn.listen` will automatically reconnect to the GCN network if the network connection is broken.

    # Listen for GCNs until the program is interrupted
    # (killed or interrupted with control-C).
    gcn.listen(port=8096, handler=process_gcn)

That was fun! Now kill the listener: if you are running the example code in an IPython Notebook, press the stop button. If you copied the example code into a standalone Python script, then kill the script by typing control-C.

## 5. Working with probability sky maps

Let's take a look at what is inside one of the LIGO/Virgo probability sky maps. They are FITS image files and can be manipulated and viewed with many commonplace FITS tools. However, they are a little unusual in two regards. First, since they are all-sky images, they are stored in the HEALPix projection, a format that is used for Planck all-sky CMB maps and by Aladin for archival all-sky survey images. Second, the value stored at each pixel is the probability that the gravitational-wave source is within that pixel.

Let's download one of them with `curl`:

    $ curl --netrc -O https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz

or with `wget`:

    $ wget --auth-no-challenge https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz

We can look at the metadata inside the FITS file by printing its header with tools like funhead from Funtools, imhead from WCSTools, or fitsheader from Astropy:

    $ funhead -a bayestar.fits.gz
    SIMPLE  =                    T / conforms to FITS standard
    BITPIX  =                    8 / array data type
    NAXIS   =                    0 / number of array dimensions
    EXTEND  =                    T
    END
    		Extension: XTENSION

    XTENSION= 'BINTABLE'           / binary table extension
    BITPIX  =                    8 / array data type
    NAXIS   =                    2 / number of array dimensions
    NAXIS1  =                 4096 / length of dimension 1
    NAXIS2  =                 3072 / length of dimension 2
    PCOUNT  =                    0 / number of group parameters
    GCOUNT  =                    1 / number of groups
    TFIELDS =                    1 / number of table fields
    TTYPE1  = 'PROB    '
    TFORM1  = '1024E   '
    TUNIT1  = 'pix-1   '
    PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
    ORDERING= 'NESTED  '           / Pixel ordering scheme, either RING or NESTED
    COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
    EXTNAME = 'XTENSION'           / name of this binary table extension
    NSIDE   =                  512 / Resolution parameter of HEALPIX
    FIRSTPIX=                    0 / First pixel # (0 based)
    LASTPIX =              3145727 / Last pixel # (0 based)
    INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
    OBJECT  = 'T125706 '           / Unique identifier for this event
    REFERENC= 'https://gracedb.ligo.org/events/T125706' / URL of this event
    INSTRUME= 'H1,L1   '           / Instruments that triggered this event
    DATE-OBS= '2010-08-30T10:11:34.322149' / UTC date of the observation
    MJD-OBS =     55438.4247028027 / modified Julian date of the observation
    DATE    = '2014-08-13T19:21:12' / UTC date of file creation
    CREATOR = 'bayestar_localize_lvalert' / Program that created this file
    ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
    COMMENT
    COMMENT This simulated detection candidate is a copy of event 11178
    COMMENT from the 2015 scenario of the "First Two Years" paper
    COMMENT (http://www.ligo.org/scientists/first2years/).
    END

There are several useful pieces of information here:

- `COORDSYS=C`, telling you that the HEALPix projection is in the Celestial (equatorial) frame (as all LIGO/Virgo probability sky maps will be).
- `OBJECT`, the unique LIGO/Virgo identifier for the event.
- `REFERENC`, a link to the candidate page in the GraceDb gravitational-wave candidate event database.
- `INSTRUME`, a list of gravitational-wave sites that triggered on the event (`H1` represents LIGO Hanford, `L1` is LIGO Livingston, and `V1` is Virgo).
- `DATE-OBS`, the UTC time of the event. In the case of a compact binary coalescence candidate, this is the time that the signal from the merger passed through the geocenter.
- `MJD-OBS`, same as `DATE-OBS`, but given as a modified Julian day.

You can view the sky map in many common FITS image viewers such as Aladin:

or DS9 (although DS9 shows HEALPix sky maps in an unusual orientation; see Figure 4 of Calabretta & Roukema 2007 for information):

Now, let's go through some examples of manipulating HEALPix sky maps programmatically. The HEALPix project provides official libraries for many languages, including C, C++, Fortran, IDL, Java, and Python. However, since this is a Python tutorial, we are going to demonstrate how to manipulate HEALPix maps with the official Python library, Healpy.

First, if you have not already downloaded an example sky map, you can do so now by having Python call `curl` on the command line:

    # Download sky map
    import subprocess
    subprocess.check_call([
        'curl', '-O', '--netrc',
        'https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz'])

    0

Next, we need to read in the file with Healpy:

    hpx = hp.read_map('bayestar.fits.gz')

    NSIDE = 512
    ORDERING = NESTED in fits file
    Ordering converted to RING

You can suppress printing informational messages while loading the file by passing the keyword argument `verbose=False`. You can read both the HEALPix image data and the FITS header by passing the `h=True` keyword argument:

    hpx, header = hp.read_map('bayestar.fits.gz', h=True, verbose=False)

The image data is a 1D array of values:

    hpx

    array([  3.08905283e-22,   3.04303258e-22,   1.21389563e-22, ...,
             9.39631031e-22,   3.19591130e-22,   2.66302276e-22])

Healpy has several useful plotting routines including `mollview` for plotting a Mollweide-projection all-sky map:

    hp.mollview(hpx)

Each entry in the array represents the probability contained within a quadrilateral pixel whose position on the sky is uniquely specified by the index in the array and the array's length. Because HEALPix pixels are equal area, we can find the number of pixels per square degree just from the length of the HEALPix array:

    npix = len(hpx)
    sky_area = 4 * 180**2 / np.pi
    sky_area / npix

    0.013113963206424481

The `pix2ang` function converts from pixel index to spherical polar coordinates; the function `ang2pix` does the reverse.

Both `pix2ang` and `ang2pix` take, as their first argument, `nside`, the lateral resolution fo the HEALPix map. You can find `nside` from the length of the image array by calling `npix2nside`:

    nside = hp.npix2nside(npix)
    nside

    512

Let's look up the right ascension and declination of pixel number 123. We'll call `pix2ang` to get the spherical polar coordinates (θ,ϕ) in radians, and then use Numpy's `rad2deg` function to convert these to right ascension and declination in degrees.

    ipix = 123
    theta, phi = hp.pix2ang(nside, ipix)
    ra = np.rad2deg(phi)
    dec = np.rad2deg(0.5 * np.pi - theta)
    ra, dec

    (129.375, 89.269029291573901)

Let's find which pixel contains the point RA=194.95, Dec=27.98.

    ra = 194.95
    dec = 27.98
    theta = 0.5 * np.pi - np.deg2rad(dec)
    phi = np.deg2rad(ra)
    ipix = hp.ang2pix(nside, theta, phi)
    ipix

    833621

Let's find the highest probability pixel. What is the probability inside it?

    ipix_max = np.argmax(hpx)
    hpx[ipix_max]

    5.7702549383975565e-05

Where is the highest probability pixel on the sky? Use `pix2ang`.

    theta, phi = hp.pix2ang(nside, ipix_max)
    ra = np.rad2deg(phi)
    dec = np.rad2deg(0.5 * np.pi - theta)
    ra, dec

    (213.22265625, -37.450292350169015)

How do we find the probability that the source is contained within a circle on the sky? First we find the pixels that are contained within the circle using `query_disc`. Note that `query_disc` takes as its arguments the Cartesian coordinates of the center of the circle, and its radius in radians.

Then, we sum the values of the HEALPix image array contained at those pixels.

    # RA, Dec, and radius of circle in degrees
    ra = 213.22
    dec = -37.45
    radius = 3.1

    # Spherical polar coordinates and radius of circle in radians
    theta = 0.5 * np.pi - np.deg2rad(dec)
    phi = np.deg2rad(ra)
    radius = np.deg2rad(radius)

    # Cartesian coordinates of center of circle
    xyz = hp.ang2vec(theta, phi)

    # Array of indices of pixels inside circle
    ipix_disc = hp.query_disc(nside, xyz, radius)

    # Probability that source is within circle
    hpx[ipix_disc].sum()

    0.056844148994230181

Similarly, we can use the `query_polygon` function to look up the indices of the pixels within a polygon (defined by the Cartesian coordinates of its vertices), and then compute the probability that the source is inside that polygon by summing the values of the pixels.

    # Vertices of polygon
    xyz = [[-0.69601758, -0.41315628, -0.58724902],
           [-0.68590811, -0.40679797, -0.60336181],
           [-0.69106913, -0.39820114, -0.60320752],
           [-0.7011786 , -0.40455945, -0.58709473]]

    # Array of indices of pixels inside polygon
    ipix_poly = hp.query_polygon(nside, xyz)

    # Probability that source is within polygon
    hpx[ipix_poly].sum()

    0.0011913816551896161

These are all of the HEALPix functions from Healpy that we will need for the remainder of the this tutorial.

Other useful Healpy functions include `ud_grade` for upsampling or downsampling a sky map, and `get_interp_val` for performing bilinear interpolation between pixels. See the Healpy tutorial for other useful operations.

## 6. Basic observability calculations with Astropy

Now we are going to teach our GCN handler how to determine whether a gravitational-wave event is observable. We are goint to use the astropy.coordinates module of the Astropy package to do 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

The LIGO/Virgo probability sky maps are always in equatorial coordinates. Once we have looked up the coordinates of the HEALPix pixels, we will use the positional astronomy features of Astropy to transform those coordinates to an alt/az 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: users 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 `hp.pix2ang`) and some of the key positional astronomy functions with Astropy.

    def prob_observable(m, header):
        """
        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 resolution of sky map
        npix = len(m)
        nside = hp.npix2nside(npix)

        # Get time now
        time = astropy.time.Time.now()
        # Or at the time of the gravitational-wave event...
        # time = astropy.time.Time(header['MJD-OBS'], format='mjd')
        # Or at a particular time...
        # time = astropy.time.Time('2015-03-01 13:55:27')

        # 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, now
        frame = astropy.coordinates.AltAz(obstime=time, location=observatory)

        # Look up (celestial) spherical polar coordinates of HEALPix grid.
        theta, phi = hp.pix2ang(nside, np.arange(npix))
        # Convert to RA, Dec.
        radecs = astropy.coordinates.SkyCoord(
            ra=phi*u.rad, dec=(0.5*np.pi - theta)*u.rad)

        # Transform grid to alt/az coordinates at observatory, now
        altaz = radecs.transform_to(frame)

        # Where is the sun, now?
        sun_altaz = astropy.coordinates.get_sun(time).transform_to(altaz)

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

        # Done!
        return prob

Finally, we need to update our GCN handler to call this function.

    # Function to call every time a GCN is received.
    # Run only for notices of type LVC_INITIAL or LVC_UPDATE.
    @gcn.handlers.include_notice_types(
        gcn.notice_types.LVC_INITIAL,
        gcn.notice_types.LVC_UPDATE)
    def process_gcn(payload, root):
        # Print the alert
        print('Got VOEvent:')
        print(payload)

        # Respond only to 'test' events.
        # VERY IMPORTANT! Replce with the following line of code
        # to respond to only real 'observation' events.
        # if root.attrib['role'] != 'observation': return
        if root.attrib['role'] != 'test': return

        # Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only
        # unmodeled burst events.
        if root.find("./What/Param[@name='Group']").attrib['value'] != 'CBC': return

        skymap, header = get_skymap(root)
        prob = prob_observable(skymap, header)
        print('Source has a %d%% chance of being observable now' % round(100 * prob))
        if prob > 0.5:
            pass # FIXME: perform some action

Let's run the new GCN handler now...

    # Listen for GCNs until the program is interrupted
    # (killed or interrupted with control-C).
    gcn.listen(port=8096, handler=process_gcn)

## 7. Submitting observation coordinates to GraceDB

Suppose you have performed some EM observations to follow up on a candidate GW event, and you now want to supply the coordinates of those observations to the LV-EM group. The first step is to obtain a robotic access password from GraceDB and add it to a netrc file (see Section 2). One can then use the Python GraceDB client in order to submit a list of coordinates to GraceDB. To install the GraceDB client package:

    $ pip install ligo-gracedb

Note: It is highly recommended to install the GraceDB client (as well as other packages described used in this tutorial) inside a virtual environment. This isolates the packages required for interacting with GraceDB from other packages on the system. Now that the GraceDB client is installed, one can use a script such as this to submit observation records to GraceDB:

    from ligo.gracedb.rest import GraceDbBasic, HTTPError

    service = 'https://gracedb.ligo.org/apibasic/'
    g = GraceDbBasic(service)

    graceid = 'T125706'
    raList = [45.0, 47.0, 49.0]
    raWidthList = 2.0
    decList = [45.0, 47.0, 49.0]
    decWidthList = 2.0
    startTimeList = ['2015-05-01T12:30:10.95', '2015-05-01T12:31:10.95', '2015-05-01T12:32:10.95']
    durationList = 100.0
    comment = 'Some text comment goes here. This is optional.'

    g.writeEMObservation(graceid, 'Test', raList, raWidthList,
        decList, decWidthList, startTimeList, durationList, comment)

In the above example, the lists of ra and dec values were chosen arbitrarily, as well as the start times and durations (or exposure times). Since the comma-separated lists have three elements, there will be three footprints associated with this observation. And since only one value was provided for the raWidth and decWidth, these widths will be assumed to apply to all three footprints.

It is possible to create the same record using curl instead:

    $ curl -X POST --netrc --data "group=Test&raList=45,47,49&raWidthList=2&decList=45,47,49&decWidthList=2&startTimeList=2015-05-01T12:30:10.95,2015-05-01T12:31:10.95,2015-05-01T12:32:10.95&durationList=100,100,100" https://gracedb.ligo.org/apibasic/events/T125706/emobservation/

However, the data argument in the above curl command is rather unwieldy--and this is only with three simple footprints. So using Python (as in the above code example) or your own favorite language and HTTP library is much preferred. We are also planning an email submission mechanism which may be eaiser for some users.

# Appendix: Full example code

Here is a complete, working GCN processing script. Copy it into a `.py` file and customize as needed.

    # Python standard library imports
    import tempfile
    import shutil
    import sys
    import glob

    # Third-party imports
    import gcn
    import gcn.handlers
    import gcn.notice_types
    import requests
    import healpy as hp
    import numpy as np
    import astropy.coordinates
    import astropy.time
    import astropy.units as u

    def get_skymap(root):
        """
        Look up URL of sky map in VOEvent XML document,
        download sky map, and parse FITS file.
        """
        # Read out URL of sky map.
        # This will be something like
        # https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gz
        skymap_url = root.find(
            "./What/Param[@name='SKYMAP_URL_FITS_BASIC']").attrib['value']

        # Send HTTP request for sky map
        response = requests.get(skymap_url, stream=True)

        # Uncomment to save VOEvent payload to file
        # open('example.xml', 'w').write(payload)

        # Raise an exception unless the download succeeded (HTTP 200 OK)
        response.raise_for_status()

        # Create a temporary file to store the downloaded FITS file
        with tempfile.NamedTemporaryFile() as tmpfile:
            # Save the FITS file to the temporary file
            shutil.copyfileobj(response.raw, tmpfile)
            tmpfile.flush()

            # Uncomment to save FITS payload to file
            # shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))

            # Read HEALPix data from the temporary file
            skymap, header = hp.read_map(tmpfile.name, h=True, verbose=False)
            header = dict(header)

        # Done!
        return skymap, header

    def prob_observable(m, header):
        """
        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 resolution of sky map
        npix = len(m)
        nside = hp.npix2nside(npix)

        # Get time now
        time = astropy.time.Time.now()
        # Or at the time of the gravitational-wave event...
        # time = astropy.time.Time(header['MJD-OBS'], format='mjd')
        # Or at a particular time...
        # time = astropy.time.Time('2015-03-01 13:55:27')

        # 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, now
        frame = astropy.coordinates.AltAz(obstime=time, location=observatory)

        # Look up (celestial) spherical polar coordinates of HEALPix grid.
        theta, phi = hp.pix2ang(nside, np.arange(npix))
        # Convert to RA, Dec.
        radecs = astropy.coordinates.SkyCoord(
            ra=phi*u.rad, dec=(0.5*np.pi - theta)*u.rad)

        # Transform grid to alt/az coordinates at observatory, now
        altaz = radecs.transform_to(frame)

        # Where is the sun, now?
        sun_altaz = astropy.coordinates.get_sun(time).transform_to(altaz)

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

        # Done!
        return prob

    # Function to call every time a GCN is received.
    # Run only for notices of type LVC_INITIAL or LVC_UPDATE.
    @gcn.handlers.include_notice_types(
        gcn.notice_types.LVC_INITIAL,
        gcn.notice_types.LVC_UPDATE)
    def process_gcn(payload, root):
        # Print the alert
        print('Got VOEvent:')
        print(payload)

        # Respond only to 'test' events.
        # VERY IMPORTANT! Replce with the following line of code
        # to respond to only real 'observation' events.
        # if root.attrib['role'] != 'observation': return
        if root.attrib['role'] != 'test': return

        # Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only
        # unmodeled burst events.
        if root.find("./What/Param[@name='Group']").attrib['value'] != 'CBC': return

        skymap, header = get_skymap(root)
        prob = prob_observable(skymap, header)
        print('Source has a %d%% chance of being observable now' % round(100 * prob))
        if prob > 0.5:
            pass # FIXME: perform some action

    # Listen for GCNs until the program is interrupted
    # (killed or interrupted with control-C).
    gcn.listen(port=8096, handler=process_gcn)