.. _gwpy-example-miscellaneous-open-data-spectrogram:

.. sectionauthor:: Duncan Macleod <duncan.macleod@ligo.org>

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 <https://gwosc.org/detector_status/>`__
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.

.. currentmodule:: gwpy.segments

Getting the segments
--------------------

First, we need to fetch the Open Data timeline segments from GWOSC.
To do that we can call the :meth:`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 <https://gwosc.org/S6/>`__).


.. plot::
   :context: reset
   :nofigs:
   :include-source:

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

For sanity, lets plot these segments:


.. plot::
   :context:
   :include-source:

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

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 :func:`abs` function to display the total amount of time
spent taking data:


.. plot::
   :context: close-figs
   :nofigs:
   :include-source:

   print(abs(h1segs.active))

.. currentmodule:: gwpy.timeseries

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
:class:`~gwpy.spectrogram.Spectrogram` for each segment.


.. plot::
   :context:
   :nofigs:
   :include-source:

   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 :meth:`~gwpy.spectrogram.Spectrogram.plot`:

.. plot::
   :context:
   :include-source:

   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()