add option to plot uniform color in sds plot, do not sort events by e…

…vent id in all_sds plot

CHANGELOG

@@ -3,6 +3,7 @@ v4.0-dev:
   * add source subcommand to calculate sorce spectra and derive source parameters including moment magnitude with the help of previous results
   * add command line arguments to dump opt and fit pkl files for later plotting
   * add command line arguments to turn all ploting on or off
+  * add option to plot uniform color in sds plot, do not sort events by event id in all_sds plot
 v3.0:
   * prefix option is more stringent now and adds prefix also to log file
   * add option to filter event catalog, add option to resolve full SEED id, add corresponding command line options

qopen/core.py

@@ -1,4 +1,4 @@
-# Copyright 2015-2017 Tom Eulenfeld, MIT license
+# Copyright 2015-2020 Tom Eulenfeld, MIT license
 """
 Qopen command line script and routines
 

qopen/imaging.py

@@ -372,13 +372,14 @@ def get_Emod(G, t):
 def plot_sds(freq, result, ax=None,
              annotate=False, va='bottom',
              seismic_moment_method=None, seismic_moment_options={},
-             cmap = 'viridis_r', vmin=None, vmax=None, max_nobs=None,
+             cmap='viridis_r', vmin=None, vmax=None, max_nobs=None,
+             color=None,
              **kwargs):
     """Plot source displacement spectrum and fitted source model"""
 
     kw = {'s': 4 * MS ** 2, 'marker': 'o', 'zorder': 10, #'linewidth': 0.5,
           'c': 'k'}
-    if 'nstations' in result and cmap is not None:
+    if color is None and 'nstations' in result and cmap is not None:
 #        Rvals = np.array(list(result['R'].values()), dtype='float')
 #        nobs = np.sum(~np.isnan(Rvals), axis=0)
         nobs = result['nstations']
@@ -391,6 +392,9 @@ def plot_sds(freq, result, ax=None,
             vmin = 0.5
         norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
         kw.update({'c': nobs, 'norm': norm, 'cmap': cmap})
+    elif color is not None:
+        kw['c'] = None
+        kw['color'] = color
     freq = np.array(freq)
     omM = np.array(result['sds'], dtype=np.float)
     if all(np.isnan(omM)):
@@ -662,8 +666,10 @@ def plot_all_sds(result, seismic_moment_method=None,
                  seismic_moment_options=None,
                  xlim=None, ylim=None, nx=None,
                  annotate=None, va='top',
+                 annotate_evid=True,
                  plot_only_ids=None,
-                 cmap = 'viridis_r', vmin=None, vmax=None,
+                 cmap='viridis_r', vmin=None, vmax=None,
+                 colors=None,
                  **kwargs):
     """Plot all source displacement spectra with fitted source models"""
     freq = np.array(result['freq'])
@@ -674,8 +680,9 @@ def plot_all_sds(result, seismic_moment_method=None,
     result = result['events']
     if plot_only_ids:
         result = {id_: r for id_, r in result.items() if id_ in plot_only_ids}
-    if 'nstations' not in list(result.values())[0] or cmap is None:
-        max_nobs = 1  # single inversion
+    if ('nstations' not in list(result.values())[0] or cmap is None or
+            colors is not None):
+        max_nobs = 1  # single inversion or uniform color
     else:
         nobs = [evres['nstations'] for evres in result.values()]
         max_nobs = np.max(nobs)
@@ -687,13 +694,22 @@ def plot_all_sds(result, seismic_moment_method=None,
     share = None
     if annotate is None:
         annotate = nx < 7
-    for i, evid in enumerate(sorted(result)):
+    color = None
+    for i, evid in enumerate(result):
+        if colors is not None:
+            try:
+                color = colors[evid]
+            except:
+                color = colors
         ax = plt.subplot(gs[i // nx, i % nx], sharex=share, sharey=share)
         sc = plot_sds(freq, result[evid], seismic_moment_method=smm, va=va,
                       seismic_moment_options=smo, ax=ax, annotate=annotate,
-                      cmap=cmap, vmin=vmin, vmax=vmax, max_nobs=max_nobs)
-        ax.annotate(evid, (0, 0), (5, 5), 'axes fraction',
-                    'offset points', ha='left', va='bottom', size='x-small')
+                      cmap=cmap, vmin=vmin, vmax=vmax, max_nobs=max_nobs,
+                      color=color)
+        if annotate_evid:
+            ax.annotate(evid, (0, 0), (5, 5), 'axes fraction',
+                        'offset points', ha='left', va='bottom',
+                        size='x-small')
         _set_gridlabels(ax, i, nx, ny, N, ylabel=r'$\omega$M (Nm)')
         if share is None:
             share = ax

qopen/util.py

@@ -4,7 +4,7 @@
 import functools
 
 import numpy as np
-from statsmodels.regression.linear_model import OLS
+from statsmodels.regression.linear_model import OLS, WLS
 from statsmodels.robust.robust_linear_model import RLM
 
 
@@ -216,22 +216,22 @@ def smooth(x, window_len=None, window='flat', method='zeros'):
     return np.convolve(w / w.sum(), s, mode='valid')
 
 
-def linear_fit(y, x, m=None, method='robust'):
+def linear_fit(y, x, m=None, method='robust', **kw):
     """Linear fit between x and y
 
     :param y,x: data
     :param m: fix slope at specific value
-    :param method: one of ('least_squares', 'robust')
+    :param method: one of ('least_squares', 'weighted', 'robust')
     :return: slope a and intercept b of y = ax + b
     """
-    Model = RLM if method == 'robust' else OLS
+    Model = RLM if method == 'robust' else WLS if method == 'weighted' else OLS
     if m is None:
         X = np.empty((len(y), 2))
         X[:, 0] = x
         X[:, 1] = 1
-        res = Model(y, X).fit()
+        res = Model(y, X, **kw).fit()
         return res.params
     else:
         X = np.ones(len(y))
-        res = Model(np.array(y) - m * np.array(x), X).fit()
+        res = Model(np.array(y) - m * np.array(x), X, **kw).fit()
         return m, res.params[0]