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 The Gravitational-Wave Open Science Centre (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)
