3. Plotting a spectrogram of all open data for 1 day¶

In order to study interferometer performance, it is common in LIGO to plot all of the data for a day, in order to determine trends, and see data-quality issues.

This is done for the LIGO-Virgo detector network, with up-to-date plots available from GWOSC.

This example demonstrates how to download data segments from GWOSC, then use those to build a day-timescale spectrogram plot of LIGO-Hanford strain data.

3.1. Getting the segments¶

First, we need to fetch the Open Data timeline segments from GWOSC. To do that we can call the DataQualityFlag.fetch_open_data() method using 'H1_DATA' as the flag (for an explanation of what this means, read up on The S6 Data Release).

from gwpy.segments import DataQualityFlag
h1segs = DataQualityFlag.fetch_open_data('H1_DATA',
                                         'Sep 16 2010', 'Sep 17 2010')

For sanity, lets plot these segments:

splot = h1segs.plot(figsize=[12, 3])
splot.show()

(png)

We see that the LIGO Hanford Observatory detector was operating for the majority of the day, with a few outages of ~30 minutes or so.

We can use the abs() function to display the total amount of time spent taking data:

print(abs(h1segs.active))

3.2. Working with strain data¶

Now, we can loop through the active segments of 'H1_DATA' and fetch the strain TimeSeries for each segment, calculating a Spectrogram for each segment.

from gwpy.timeseries import TimeSeries
spectrograms = []
for start, end in h1segs.active:
    h1strain = TimeSeries.fetch_open_data('H1', start, end, verbose=True)
    specgram = h1strain.spectrogram(30, fftlength=4) ** (1/2.)
    spectrograms.append(specgram)

Finally, we can build a plot():

from gwpy.plot import Plot
plot = Plot(figsize=(12, 6))
ax = plot.gca()
for specgram in spectrograms:
    ax.imshow(specgram)
ax.set_xscale('auto-gps', epoch='Sep 16 2010')
ax.set_xlim('Sep 16 2010', 'Sep 17 2010')
ax.set_ylim(40, 2000)
ax.set_yscale('log')
ax.set_ylabel('Frequency [Hz]')
ax.set_title('LIGO-Hanford strain data')
ax.colorbar(cmap='viridis', norm='log', clim=(1e-23, 1e-19),
            label=r'Strain noise [1/$\sqrt{\mathrm{Hz}}$]')
plot.add_segments_bar(h1segs)
plot.show()

(png)