Merge pull request #716 from duncanmmacleod/timeseries-shift

Added new TimeSeries.shift method

examples/signal/gw150914.py

@@ -34,7 +34,7 @@
 # First we download the raw strain data from the LOSC public archive:
 
 from gwpy.timeseries import TimeSeries
-data = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
+hdata = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
 
 # Next we can design a zero-pole-gain filter to remove the extranious noise.
 
@@ -43,13 +43,13 @@
 # high frequency content
 
 from gwpy.signal import filter_design
-bp = filter_design.bandpass(50, 250, data.sample_rate)
+bp = filter_design.bandpass(50, 250, hdata.sample_rate)
 
 # Now we want to combine the bandpass with a series of
 # :meth:`~gwpy.signal.filter_design.notch` filters, so we create those
 # for the first three harmonics of the 60 Hz AC mains power:
 
-notches = [filter_design.notch(line, data.sample_rate) for
+notches = [filter_design.notch(line, hdata.sample_rate) for
            line in (60, 120, 180)]
 
 # and concatenate each of our filters together to create a single ZPK:
@@ -60,7 +60,7 @@
 # to filter both backwards and forwards to preserve the correct phase
 # at all frequencies
 
-b = data.filter(zpk, filtfilt=True)
+hfilt = hdata.filter(zpk, filtfilt=True)
 
 # .. note::
 #
@@ -69,37 +69,73 @@
 #    analogue filters (perhaps by passing `analog=True` to the filter design
 #    method), the easiest application would be `TimeSeries.zpk`
 
+# The :mod:`~gwpy.signal.filter_design` methods return infinite impulse
+# response filters by default, which, when applied, corrupt a small amount of
+# data at the beginning and the end of our original `TimeSeries`.
+# We can discard those data using the :meth:`~TimeSeries.crop` method
+# (for consistency we apply this to both data series):
+
+hdata = hdata.crop(*hdata.span.contract(1))
+hfilt = hfilt.crop(*hfilt.span.contract(1))
+
 # Finally, we can :meth:`~TimeSeries.plot` the original and filtered data,
 # adding some code to prettify the figure:
 
 from gwpy.plotter import TimeSeriesPlot
-plot = TimeSeriesPlot(
-    data.crop(*data.span.contract(1)),
-    b.crop(*b.span.contract(1)),
-    figsize=[12, 8], sep=True, sharex=True, color='gwpy:ligo-hanford')
-plot.axes[0].set_title('LIGO-Hanford strain data around GW150914')
-plot.axes[0].text(
-    1.0, 1.0, 'Unfiltered data',
-    transform=plot.axes[0].transAxes, ha='right')
-plot.axes[0].set_ylabel('Amplitude [strain]', y=-0.2)
-plot.axes[1].text(
-    1.0, 1.0, '50-250\,Hz bandpass, notches at 60, 120, 180 Hz',
-    transform=plot.axes[1].transAxes, ha='right')
+plot = TimeSeriesPlot(hdata, hfilt, figsize=[12, 8], sep=True, sharex=True,
+                      color='gwpy:ligo-hanford')
+ax1, ax2 = plot.axes
+ax1.set_title('LIGO-Hanford strain data around GW150914')
+ax1.text(1.0, 1.0, 'Unfiltered data', transform=ax1.transAxes, ha='right')
+ax1.set_ylabel('Amplitude [strain]', y=-0.2)
+ax2.set_ylabel('')
+ax2.text(1.0, 1.0, '50-250\,Hz bandpass, notches at 60, 120, 180 Hz',
+         transform=ax2.transAxes, ha='right')
 plot.show()
 
 # We see now a spike around 16 seconds into the data, so let's zoom into
-# that time by using :meth:`~TimeSeries.crop` and :meth:`~TimeSeries.plot`:
+# that time (and prettify):
 
-plot = b.crop(1126259462, 1126259462.6).plot(
-    figsize=[12, 4], color='gwpy:ligo-hanford')
+plot = hfilt.plot(figsize=[12, 4], color='gwpy:ligo-hanford')
 ax = plot.gca()
 ax.set_title('LIGO-Hanford strain data around GW150914')
 ax.set_ylabel('Amplitude [strain]')
-ax.set_epoch(1126259462.427)
+ax.set_xlim(1126259462, 1126259462.6)
+ax.set_xscale('seconds', epoch=1126259462)
 plot.show()
 
 # 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 = TimeSeriesPlot(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')
+ax.set_title('LIGO strain data around GW150914')
+ax.set_xlim(1126259462, 1126259462.6)
+ax.set_xscale('seconds', epoch=1126259462)
+ax.set_ylabel('Amplitude [strain]')
+ax.set_ylim(-1e-21, 1e-21)
+ax.legend()
+plot.show()
+
 # The above filtering is (almost) the same as what was applied to LIGO data
 # to produce Figure 1 in
 # `Abbott et al. (2016) <https://doi.org/10.1103/PhysRevLett.116.061102>`_

gwpy/plotter/timeseries.py

@@ -32,6 +32,7 @@
 from ..time import (Time, LIGOTimeGPS, to_gps)
 from ..segments import SegmentList
 from .decorators import auto_refresh
+from .gps import GPSScale
 from .series import (SeriesPlot, SeriesAxes)
 
 __author__ = "Duncan Macleod <duncan.macleod@ligo.org>"
@@ -62,8 +63,10 @@ def draw(self, *args, **kwargs):
 
         # dynamically set x-axis label
         nolabel = self.get_xlabel() == '_auto'
-        if nolabel:
+        if nolabel and isinstance(self.xaxis._scale, GPSScale):
             self.auto_gps_label()
+        elif nolabel:
+            self.set_xlabel('')
 
         # draw
         super(TimeSeriesAxes, self).draw(*args, **kwargs)
@@ -78,8 +81,7 @@ def draw(self, *args, **kwargs):
 
     def set_xscale(self, scale, *args, **kwargs):
         super(TimeSeriesAxes, self).set_xscale(scale, *args, **kwargs)
-        if scale != 'auto-gps' and self.get_xlabel() == '_auto':
-            self.set_xlabel('')
+
     set_xscale.__doc__ = SeriesAxes.set_xscale.__doc__
 
     def auto_gps_label(self):

gwpy/tests/test_timeseries.py

@@ -200,6 +200,21 @@ def test_duration(self, array):
 
     # -- test methods ---------------------------
 
+    def test_shift(self):
+        a = self.create()
+        t0 = a.t0.copy()
+        a.shift(5)
+        assert a.t0 == t0 + 5 * t0.unit
+
+        a.shift('1 hour')
+        assert a.t0 == t0 + 3605 * t0.unit
+
+        a.shift(-0.007)
+        assert a.t0 == t0 + (3604.993) * t0.unit
+
+        with pytest.raises(ValueError):
+            a.shift('1 Hz')
+
     def test_plot(self, array):
         with rc_context(rc={'text.usetex': False}):
             plot = array.plot()

gwpy/timeseries/core.py

@@ -593,6 +593,32 @@ def get(cls, channel, start, end, pad=None, dtype=None, verbose=False,
 
     # -- utilities ------------------------------
 
+    def shift(self, delta):
+        """Shift this `TimeSeries` forward in time by ``delta``
+
+        This modifies the series in-place.
+
+        Parameters
+        ----------
+        delta : `float`, `~astropy.units.Quantity`, `str`
+            The amount by which to shift (in seconds if `float`), give
+            a negative value to shift backwards in time
+
+        Examples
+        --------
+        >>> from gwpy.timeseries import TimeSeries
+        >>> a = TimeSeries([1, 2, 3, 4, 5], t0=0, dt=1)
+        >>> print(a.t0)
+        0.0 s
+        >>> a.shift(5)
+        >>> print(a.t0)
+        5.0 s
+        >>> a.shift('-1 hour')
+        -3595.0 s
+        """
+        delta = units.Quantity(delta, 's')
+        self.t0 += delta
+
     def plot(self, **kwargs):
         """Plot the data for this timeseries