How to use the gwpy.timeseries.TimeSeries.fetch_open_data function in gwpy

The LIGO Laboratory has publicly released the strain data around the time of
the GW150914 gravitational-wave detection; we can use these to calculate
and display the spectral sensitivity of each of the detectors at that time.
"""

__author__ = "Duncan Macleod "
__currentmodule__ = 'gwpy.frequencyseries'

# In order to generate a `FrequencySeries` we need to import the
# `~gwpy.timeseries.TimeSeries` and use
# :meth:`~gwpy.timeseries.TimeSeries.fetch_open_data` to download the strain
# records:

from gwpy.timeseries import TimeSeries
lho = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
llo = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478)

# We can then call the :meth:`~gwpy.timeseries.TimeSeries.asd` method to
# calculated the amplitude spectral density for each
# `~gwpy.timeseries.TimeSeries`:
lhoasd = lho.asd(4, 2)
lloasd = llo.asd(4, 2)

# We can then :meth:`~FrequencySeries.plot` the spectra using the 'standard'
# colour scheme:

plot = lhoasd.plot(label='LIGO-Hanford', color='gwpy:ligo-hanford')
ax = plot.gca()
ax.plot(lloasd, label='LIGO-Livingston', color='gwpy:ligo-livingston')
ax.set_xlim(10, 2000)
ax.set_ylim(5e-24, 1e-21)
ax.legend(frameon=False, bbox_to_anchor=(1., 1.), loc='lower right', ncol=2)

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.

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`:

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

power spectral density (PSD) of a given set of data.
It is used to measure the 'Gaussianity' of those data, where a value of 1
indicates Gaussian behaviour, less than 1 indicates coherent variations,
and greater than 1 indicates incoherent variation.
It is a useful measure of the quality of the strain data being generated
and recorded at a LIGO site.
"""

__author__ = "Duncan Macleod "
__currentmodule__ = 'gwpy.frequencyseries'

# To demonstate this, we can load some data from the LIGO Livingston
# intereferometer around the time of the GW151226 gravitational wave detection:

from gwpy.timeseries import TimeSeries
gwdata = TimeSeries.fetch_open_data('L1', 'Dec 26 2015 03:37',
                                    'Dec 26 2015 03:47', verbose=True)

# Next, we can calculate a Rayleigh statistic `FrequencySeries` using the
# :meth:`~gwpy.timeseries.TimeSeries.rayleigh_spectrum` method of the
# `~gwpy.timeseries.TimeSeries` with a 2-second FFT and 1-second overlap (50%):

rayleigh = gwdata.rayleigh_spectrum(2, 1)

# For easy comparison, we can calculate the spectral sensitivity ASD of the
# strain data and plot both on the same figure:

from gwpy.plot import Plot
plot = Plot(gwdata.asd(2, 1), rayleigh, geometry=(2, 1), sharex=True,
            xscale='log', xlim=(30, 1500))
asdax, rayax = plot.axes

is the distance to which a canonical binary neutron star (BNS) inspiral
(with two 1.4 solar mass components) would be detected with a
signal-to-noise ratio (SNR) of 8.

We can estimate using :func:`gwpy.astro.inspiral_range` after calculating
the power-spectral density (PSD) of the strain readout for a detector, and
can plot the variation over time by looping over a power spectral density
:class:`~gwpy.spectrogram.Spectrogram`.
"""

# First, we need to load some data, for this we can use the
# `public data `__
# around the GW150914 event:

from gwpy.timeseries import TimeSeries
h1 = TimeSeries.fetch_open_data('H1', 1126257414, 1126261510)
l1 = TimeSeries.fetch_open_data('L1', 1126257414, 1126261510)

# and then calculating the PSD spectrogram:

h1spec = h1.spectrogram(30, fftlength=4)
l1spec = l1.spectrogram(30, fftlength=4)

# To calculate the inspiral range variation, we need to create a
# :class:`~gwpy.timeseries.TimeSeries` in which to store the values, then
# loop over each PSD bin in the spectrogram, calculating the
# :func:`gwpy.astro.inspiral_range` for each one:

import numpy
from gwpy.astro import inspiral_range
h1range = TimeSeries(numpy.zeros(len(h1spec)),
                     dt=h1spec.dt, t0=h1spec.t0, unit='Mpc')

Below we use this algorithm to visualise GW170817, a gravitational-wave
signal from two merging neutron stars. In the LIGO-Livingston (L1) detector,
the end of this signal coincides with a very loud glitch (`Phys. Rev. Lett.
vol. 119, p. 161101 `_).
"""

__author__ = "Alex Urban "
__currentmodule__ = 'gwpy.timeseries'

# First, we need to download the `TimeSeries` record of L1 strain measurements
# from |GWOSC|_:

from gwosc import datasets
from gwpy.timeseries import TimeSeries
gps = datasets.event_gps('GW170817')
data = TimeSeries.fetch_open_data('L1', gps-34, gps+34, tag='C00')

# We can Q-transform these data and scan over time-frequency planes to
# find the one with the most significant tile near the time of merger:

from gwpy.segments import Segment
search = Segment(gps-0.25, gps+0.25)
qgram = data.q_gram(qrange=(4, 150), search=search, mismatch=0.35)

# .. note::
#
#    To recover as much signal as possible, in practice we should suppress
#    background noise by whitening the data before the Q-transform. This
#    can be done with :meth:`TimeSeries.whiten`.

# Finally, we can plot the loudest time-frequency plane, focusing on a
# specific window around the merger time:

One of the most useful methods of visualising gravitational-wave data is to
use a spectrogram, highlighting the frequency-domain content of some data
over a number of time steps.

For this example we can use the public data around the GW150914 detection.
"""

__author__ = "Duncan Macleod "
__currentmodule__ = 'gwpy.timeseries'

# First, we import the `TimeSeries` and call
# :meth:`TimeSeries.fetch_open_data` the download the strain
# data for the LIGO-Hanford interferometer
from gwpy.timeseries import TimeSeries
data = TimeSeries.fetch_open_data(
    'H1', 'Sep 14 2015 09:45', 'Sep 14 2015 09:55')

# Next, we can calculate a `~gwpy.spectrogram.Spectrogram` using the
# :meth:`spectrogram` method of the `TimeSeries` over a 2-second stride
# with a 1-second FFT and # .5-second overlap (50%):
specgram = data.spectrogram(2, fftlength=1, overlap=.5) ** (1/2.)

# .. note::
#    :meth:`TimeSeries.spectrogram` returns a Power Spectral Density (PSD)
#    `~gwpy.spectrogram.Spectrogram` by default, so we use the ``** (1/2.)``
#    to convert this into a (more familiar) Amplitude Spectral Density.

# Finally, we can make a plot using the
# :meth:`~gwpy.spectrogram.Spectrogram.plot` method
plot = specgram.imshow(norm='log', vmin=5e-24, vmax=1e-19)
ax = plot.gca()

to calculate discrete PSDs for each stride. This is fine for long-duration
data, but give poor resolution when studying short-duration phenomena.

The `~TimeSeries.spectrogram2` method allows for highly-overlapping FFT
calculations to over-sample the frequency content of the input `TimeSeries`
to produce a much more feature-rich output.
"""

__author__ = "Duncan Macleod "
__currentmodule__ = 'gwpy.timeseries'

# To demonstrate this, we can download some data associated with the
# gravitational-wave event GW510914:

from gwpy.timeseries import TimeSeries
lho = TimeSeries.fetch_open_data('H1', 1126259458, 1126259467, verbose=True)

# and can :meth:`~TimeSeries.highpass` and :meth:`~TimeSeries.whiten`
# the remove low-frequency noise and try and enhance low-amplitude signals
# across the middle of the frequency band:

hp = lho.highpass(20)
white = hp.whiten(4, 2).crop(1126259460, 1126259465)

# .. note::
#
#    We chose to :meth:`~TimeSeries.crop` out the leading and trailing 2
#    seconds of the the whitened data series here to remove any filtering
#    artefacts that may have been introduced.

# Now we can call the `~TimeSeries.spectrogram2` method of `gwdata` to
# calculate our over-dense `~gwpy.spectrogram.Spectrogram`, using a

However, an actual astrophysical search algorithm detects signals by
calculating the signal-to-noise ratio (SNR) of data for each in a large bank
of signal models, known as templates.

Using |pycbc|_ (the actual search code), we can do that.
"""

__author__ = "Duncan Macleod "
__currentmodule__ = 'gwpy.timeseries'

# First, as always, we fetch some of the public data from the LIGO Open
# Science Center:

from gwpy.timeseries import TimeSeries
data = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)

# and condition it by applying a highpass filter at 15 Hz
high = data.highpass(15)

# This is important to remove noise at lower frequencies that isn't
# accurately calibrated, and swamps smaller noises at higher frequencies.

# For this example, we want to calculate the SNR over a 4 second segment, so
# we calculate a Power Spectral Density with a 4 second FFT length (using all
# of the data), then :meth:`~TimeSeries.crop` the data:

psd = high.psd(4, 2)
zoom = high.crop(1126259460, 1126259464)

# In order to calculate signal-to-noise ratio, we need a signal model
# against which to compare our data.

ax = plot.gca()
ax.set_title('LIGO-Hanford strain data around GW150914')
ax.set_ylabel('Amplitude [strain]')
ax.set_xlim(1126259462, 1126259462.6)
ax.set_xscale('seconds', epoch=1126259462)
plot.show()
plot.close()  # hide

# Congratulations, you have succesfully filtered LIGO data to uncover the
# first ever directly-detected gravitational wave signal, GW150914!
# But wait, what about LIGO-Livingston?
# We can easily add that to our figure by following the same procedure.
#
# First, we load the new data

ldata = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478)
lfilt = ldata.filter(zpk, filtfilt=True)

# The article announcing the detector told us that the signals were
# separated by 6.9 ms between detectors, and that the relative orientation
# of those detectors means that we need to invert the data from one
# before comparing them, so we apply those corrections:

lfilt.shift('6.9ms')
lfilt *= -1

# and finally make a new plot with both detectors:

plot = Plot(figsize=[12, 4])
ax = plot.gca()
ax.plot(hfilt, label='LIGO-Hanford', color='gwpy:ligo-hanford')
ax.plot(lfilt, label='LIGO-Livingston', color='gwpy:ligo-livingston')