PPSD: add option for cumulative plot (i.e. non-exceedence) + test

CHANGELOG.txt

@@ -61,6 +61,9 @@ master:
        (matplotlibs new viridis colormap). The old PQLX colormap is provided by
        `obspy.imaging.cm.pqlx` and can be used with
        `PPSD.plot(..., cmap=...)`.
+     - new option `PPSD.plot(..., cumulative=True)` for a cumulative plot of
+       the histogram, i.e. a non-exceedence percentage visualization, similar
+       to the `percentile` option.
 
 0.10.x:
   - obspy.fdsn:

misc/docs/source/tutorial/code_snippets/probabilistic_power_spectral_density.rst

@@ -85,6 +85,13 @@ window..
 
 .. plot:: tutorial/code_snippets/probabilistic_power_spectral_density.py
 
+A (for each frequency bin) cumulative version of the histogram can also be
+visualized:
+
+    >>> ppsd.plot(cumulative=True)
+
+.. plot:: tutorial/code_snippets/probabilistic_power_spectral_density3.py
+
 To use the colormap used by PQLX / [McNamara2004]_ you can import and use that
 colormap from :mod:`obspy.imaging.cm`:
 

misc/docs/source/tutorial/code_snippets/probabilistic_power_spectral_density3.py

@@ -0,0 +1,14 @@
+from obspy import read
+from obspy.signal import PPSD
+from obspy.io.xseed import Parser
+
+
+st = read("http://examples.obspy.org/BW.KW1..EHZ.D.2011.037")
+parser = Parser("http://examples.obspy.org/dataless.seed.BW_KW1")
+ppsd = PPSD(st[0].stats, metadata=parser)
+ppsd.add(st)
+
+st = read("http://examples.obspy.org/BW.KW1..EHZ.D.2011.038")
+ppsd.add(st)
+
+ppsd.plot(cumulative=True)

obspy/imaging/cm.py

@@ -26,14 +26,17 @@
 from matplotlib.colors import LinearSegmentedColormap
 
 
-def _get_cmap(name, reverse=False):
+def _get_cmap(name, lut=None, reverse=False):
     """
     Load a :class:`~matplotlib.colors.LinearSegmentedColormap` from
     `segmentdata` dictionary saved as numpy compressed binary data.
 
     :type name: str
     :param name: Name of colormap to load, same as filename in
         `obspy/imaging/data` without `.npz` file suffix.
+    :type lut: int
+    :param lut: Specifies the number of discrete color values in the colormap.
+        `None` to use matplotlib default value (continuous colormap).
     :type reverse: bool
     :param reverse: Whether to return the specified colormap reverted.
     :rtype: :class:`~matplotlib.colors.LinearSegmentedColormap`
@@ -51,7 +54,8 @@ def _get_cmap(name, reverse=False):
             # copied from matplotlib source, cm.py@f7a578656abc2b2c13 line 47
             data_r[key] = [(1.0 - x, y1, y0) for x, y0, y1 in reversed(val)]
         data = data_r
-    cmap = LinearSegmentedColormap(name=name, segmentdata=data)
+    kwargs = lut and {"N": lut} or {}
+    cmap = LinearSegmentedColormap(name=name, segmentdata=data, **kwargs)
     return cmap
 
 viridis = _get_cmap("viridis")

obspy/imaging/tests/test_ppsd.py

@@ -8,28 +8,25 @@
 
 import os
 import unittest
-from copy import deepcopy
 
 import matplotlib.pyplot as plt
 
 from obspy.core.util.testing import ImageComparison
 from obspy.signal.tests.test_spectral_estimation import _get_ppsd
 
 
-ppsd = _get_ppsd()
-
-
 class PPSDTestCase(unittest.TestCase):
     """
     Test cases for PPSD plotting.
     """
     def setUp(self):
         # directory where the test files are located
         self.path = os.path.join(os.path.dirname(__file__), 'images')
-        self.ppsd = deepcopy(ppsd)
+        self.ppsd = _get_ppsd()
 
     def test_ppsd_plot(self):
         """
+        Test plot of ppsd example data, normal (non-cumulative) style.
         """
         with ImageComparison(self.path, 'ppsd.png') as ic:
             self.ppsd.plot(
@@ -43,6 +40,21 @@ def test_ppsd_plot(self):
             plt.draw()
             fig.savefig(ic.name)
 
+    def test_ppsd_plot_cumulative(self):
+        """
+        Test plot of ppsd example data, cumulative style.
+        """
+        with ImageComparison(self.path, 'ppsd_cumulative.png') as ic:
+            self.ppsd.plot(
+                show=False, show_coverage=True, show_histogram=True,
+                show_noise_models=True, grid=True, period_lim=(0.02, 100),
+                cumulative=True)
+            fig = plt.gcf()
+            ax = fig.axes[0]
+            ax.set_ylim(-160, -130)
+            plt.draw()
+            fig.savefig(ic.name)
+
 
 def suite():
     return unittest.makeSuite(PPSDTestCase, 'test')

obspy/signal/spectral_estimation.py

@@ -32,6 +32,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib import mlab
+from matplotlib.colors import LinearSegmentedColormap
 from matplotlib.dates import date2num
 from matplotlib.mlab import detrend_none, window_hanning
 from matplotlib.ticker import FormatStrFormatter
@@ -888,7 +889,8 @@ def plot(self, filename=None, show_coverage=True, show_histogram=True,
              show_percentiles=False, percentiles=[0, 25, 50, 75, 100],
              show_noise_models=True, grid=True, show=True,
              max_percentage=30, period_lim=(0.01, 179), show_mode=False,
-             show_mean=False, cmap=obspy_sequential):
+             show_mean=False, cmap=obspy_sequential, cumulative=False,
+             cumulative_number_of_colors=20):
         """
         Plot the 2D histogram of the current PPSD.
         If a filename is specified the plot is saved to this file, otherwise
@@ -927,6 +929,16 @@ def plot(self, filename=None, show_coverage=True, show_histogram=True,
         :type cmap: :class:`matplotlib.colors.Colormap`
         :param cmap: Colormap to use for the plot. To use the color map like in
             PQLX, [McNamara2004]_ use :const:`obspy.imaging.cm.pqlx`.
+        :type cumulative: bool
+        :param cumulative: Can be set to `True` to show a cumulative
+            representation of the histogram, i.e. showing color coded for each
+            frequency/amplitude bin at what percentage in time the value is
+            not exceeded by the data (similar to the `percentile` option but
+            continuously and color coded over the whole area). `max_percentage`
+            is ignored when this option is specified.
+        :type cumulative_number_of_colors: int
+        :param cumulative_number_of_colors: Number of discrete color shades to
+            use, `None` for a continuous colormap.
         """
         # check if any data has been added yet
         if self.hist_stack is None:
@@ -945,9 +957,24 @@ def plot(self, filename=None, show_coverage=True, show_histogram=True,
             ax = fig.add_subplot(111)
 
         if show_histogram:
-            ppsd = ax.pcolormesh(X, Y, hist_stack.T, cmap=cmap)
+            label = "[%]"
+            data = hist_stack
+            if cumulative:
+                data = data.cumsum(axis=1)
+                data = np.multiply(data.T, 100.0/data.max(axis=1)).T
+                if max_percentage is not None:
+                    msg = ("Parameter 'max_percentage' is ignored when "
+                           "'cumulative=True'.")
+                    warnings.warn(msg)
+                max_percentage = 100
+                label = "non-exceedance (cumulative) [%]"
+                if cumulative_number_of_colors is not None:
+                    cmap = LinearSegmentedColormap(
+                        name=cmap.name, segmentdata=cmap._segmentdata,
+                        N=cumulative_number_of_colors)
+            ppsd = ax.pcolormesh(X, Y, data.T, cmap=cmap)
             cb = plt.colorbar(ppsd, ax=ax)
-            cb.set_label("[%]")
+            cb.set_label(label)
             if max_percentage is not None:
                 color_limits = (0, max_percentage)
                 ppsd.set_clim(*color_limits)