V2 tt from velocity model

sourcespec2/input/augment_traces.py

@@ -19,6 +19,7 @@
 import contextlib
 from obspy.core.util import AttribDict
 from .instrument_type import get_instrument_type
+from ..setup import config
 logger = logging.getLogger(__name__.rsplit('.', maxsplit=1)[-1])
 
 
@@ -66,41 +67,68 @@ def _check_instrtype(trace):
     :param trace: ObsPy Trace object
     :type trace: :class:`obspy.core.trace.Trace`
     """
-    inv = trace.stats.inventory
-    if not inv:
-        raise RuntimeError(
-            f'{trace.id}: cannot get instrtype from inventory: '
-            'inventory is empty: skipping trace')
+    # If traces are already corrected, instrument response may not be present
+    # In this case, we just write the correct units to the trace headers
     instrtype = trace.stats.instrtype
-    new_instrtype = None
-    try:
-        units = inv.get_response(trace.id, trace.stats.starttime).\
-            instrument_sensitivity.input_units
-    except Exception as e:
-        # inventory attached to trace has only one channel
-        chan = inv[0][0][0]
-        start_date = chan.start_date
-        end_date = chan.end_date
-        raise RuntimeError(
-            f'{trace.id}: cannot get units from inventory.\n'
-            f'> {e.__class__.__name__}: {e}\n'
-            f'> Channel start/end date: {start_date} {end_date}\n'
-            f'> Trace start time: {trace.stats.starttime}\n'
-            '> Skipping trace'
-        ) from e
-    trace.stats.units = units
-    if units.lower() == 'm' and trace.stats.instrtype != 'disp':
-        new_instrtype = 'disp'
-    if units.lower() == 'm/s' and instrtype not in ['shortp', 'broadb']:
-        new_instrtype = 'broadb'
-    if units.lower() == 'm/s**2' and instrtype != 'acc':
-        new_instrtype = 'acc'
-    if new_instrtype is not None:
-        logger.warning(
-            f'{trace.id}: instrument response units are "{units}" but '
-            f'instrument type is "{instrtype}". Changing instrument type '
-            f'to "{new_instrtype}"')
-        trace.stats.instrtype = new_instrtype
+    if config.correct_instrumental_response is False:
+        # Override instrtype from config.trace_units, which may be
+        # different from native instrtype inferred from channel code
+        if config.trace_units:
+            config.trace_units = config.trace_units.lower()
+        if config.trace_units == 'vel':
+            if instrtype not in ('shortp', 'broadb'):
+                instrtype = 'broadb'
+        elif config.trace_units in ('acc', 'disp'):
+            instrtype = config.trace_units
+        trace.stats.instrtype = instrtype
+        if instrtype == 'disp':
+            trace.stats.units = 'm'
+        elif instrtype in ('shortp', 'broadb'):
+            trace.stats.units = 'm/s'
+        elif instrtype == 'acc':
+            trace.stats.units = 'm/s**2'
+        else:
+            raise RuntimeError(
+                f'{trace.id}: cannot get units for instrument type ' \
+                    '{instrtype}.\n'
+                f'> Trace start time: {trace.stats.starttime}\n'
+                '> Skipping trace'
+            )
+    else:
+        inv = trace.stats.inventory
+        if not inv:
+            raise RuntimeError(
+                f'{trace.id}: cannot get instrtype from inventory: '
+                'inventory is empty: skipping trace')
+        new_instrtype = None
+        try:
+            units = inv.get_response(trace.id, trace.stats.starttime).\
+                instrument_sensitivity.input_units
+        except Exception as e:
+            # inventory attached to trace has only one channel
+            chan = inv[0][0][0]
+            start_date = chan.start_date
+            end_date = chan.end_date
+            raise RuntimeError(
+                f'{trace.id}: cannot get units from inventory.\n'
+                f'> {e.__class__.__name__}: {e}\n'
+                f'> Channel start/end date: {start_date} {end_date}\n'
+                f'> Trace start time: {trace.stats.starttime}\n'
+                '> Skipping trace'
+            ) from e
+        trace.stats.units = units
+        if units.lower() == 'm' and trace.stats.instrtype != 'disp':
+            new_instrtype = 'disp'
+        if units.lower() == 'm/s' and instrtype not in ['shortp', 'broadb']:
+            new_instrtype = 'broadb'
+        if units.lower() == 'm/s**2' and instrtype != 'acc':
+            new_instrtype = 'acc'
+        if new_instrtype is not None:
+            logger.warning(
+                f'{trace.id}: instrument response units are "{units}" but '
+                f'instrument type is "{instrtype}". Changing instrument type '
+                f'to "{new_instrtype}"')
+            trace.stats.instrtype = new_instrtype
 
 
 def _add_coords(trace):

sourcespec2/source_spec.py

@@ -133,6 +133,10 @@ def ssp_run(st, inventory, ssp_event, picks, allow_exit=False):
     from .ssp_summary_statistics import compute_summary_statistics
     compute_summary_statistics(sspec_output)
 
+    # Add run_info to output object
+    from .ssp_output import _update_run_info
+    _update_run_info(sspec_output)
+
     return (proc_st, spec_st, specnoise_st, weight_st, sspec_output)
 
 

sourcespec2/ssp_output.py

@@ -437,14 +437,13 @@ def _make_symlinks():
         os.symlink(filename, linkname)
 
 
-def write_output(sspec_output):
+def _update_run_info(sspec_output):
     """
-    Write results into different formats.
+    Add run info to output object
 
     :param sspec_output: Output of the spectral inversion.
     :type sspec_output: :class:`~sourcespec.ssp_data_types.SourceSpecOutput`
     """
-    # Add run info to output object
     run_info = sspec_output.run_info
     run_info.SourceSpec_version = get_versions()['version']
     config.end_of_run = datetime.now()
@@ -458,6 +457,15 @@ def write_output(sspec_output):
     run_info.agency_full_name = config.agency_full_name
     run_info.agency_short_name = config.agency_short_name
     run_info.agency_url = config.agency_url
+
+
+def write_output(sspec_output):
+    """
+    Write results into different formats.
+
+    :param sspec_output: Output of the spectral inversion.
+    :type sspec_output: :class:`~sourcespec.ssp_data_types.SourceSpecOutput`
+    """
     # Symlink input files into output directory
     _make_symlinks()
     # Write to parfile (deprecated)

sourcespec2/ssp_plot_stations.py

@@ -57,7 +57,7 @@
         'ignore', category=RuntimeWarning, module='shapely')
 # Make pyproj calculations faster by using the global context for
 # pyproj < 3.7
-if pyproj.__version__ < '3.7':
+if '3.0' <= pyproj.__version__ < '3.7':
     pyproj.set_use_global_context(active=True)
 TILER = {
     'hillshade': EsriHillshade,
@@ -976,16 +976,16 @@ def _savefig(fig, vname):
     :param vname: name of the value.
     :type vname: str
     """
-    evid = config.event.event_id
-    figfile_base = os.path.join(config.options.outdir, evid)
-    figfile_base += f'.map_{vname}.'
-    fmt = config.plot_save_format
-    if fmt == 'pdf_multipage':
-        fmt = 'pdf'
-    figfile = figfile_base + fmt
     if config.plot_show:
         plt.show()
     if config.plot_save:
+        evid = config.event.event_id
+        figfile_base = os.path.join(config.options.outdir, evid)
+        figfile_base += f'.map_{vname}.'
+        fmt = config.plot_save_format
+        if fmt == 'pdf_multipage':
+            fmt = 'pdf'
+        figfile = figfile_base + fmt
         savefig(fig, figfile, fmt, quantize_colors=False, bbox_inches='tight')
         if vname == 'mag':
             logger.info(f'Station-magnitude map saved to: {figfile}')
@@ -1089,7 +1089,18 @@ def plot_stations(sspec_output):
     summary_mag_err = summary_uncertainties['Mw']
     basemap_data = _make_basemap(maxdist)
     # make a copy of the basemap data for the second plot, before modifying it
-    basemap_data2 = deepcopy(basemap_data)
+    try:
+        basemap_data2 = deepcopy(basemap_data)
+    except NotImplementedError:
+        # Workaround for older matplotlib versions
+        import io, pickle
+        basemap_data2 = {}
+        for key, ax in basemap_data.items():
+            buf = io.BytesIO()
+            pickle.dump(ax, buf)
+            buf.seek(0)
+            ax2 = pickle.load(buf)
+            basemap_data2[key] = ax2
     _make_station_map(
         basemap_data, lonlat_dist, st_ids,
         mag, mag_outliers, summary_mag, summary_mag_err, 'mag')