General updates and bugfixes

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pysep-adjtomo"
-version = "0.3.1"
+version = "0.3.2"
 description = "Python Seismogram Extraction and Processing"
 readme = "README.md"
 requires-python = ">=3.7"

pysep/__init__.py

@@ -1,5 +1,7 @@
 import logging
 
+__version__ = "0.3.2"
+
 logger = logging.getLogger("pysep")
 logger.setLevel("INFO")
 logger.propagate = 0

pysep/pysep.py

@@ -29,7 +29,7 @@
 from obspy.core.event import Event, Origin, Magnitude
 from obspy.geodetics import kilometer2degrees
 
-from pysep import logger
+from pysep import logger, __version__
 from pysep.utils.cap_sac import (append_sac_headers, write_cap_weights_files,
                                  format_sac_header_w_taup_traveltimes,
                                  format_sac_headers_post_rotation)
@@ -55,7 +55,7 @@ def __init__(self, config_file=None, event_selection="default",
                  networks="*", stations="*", locations="*", channels="*",
                  event_latitude=None, event_longitude=None, event_depth_km=None,
                  event_magnitude=None, remove_response=True,
-                 remove_clipped=False, remove_insufficient_length=True,
+                 remove_clipped=True, remove_insufficient_length=True,
                  remove_masked_data=True,
                  water_level=60, detrend=True, demean=True, taper_percentage=0,
                  rotate=None, pre_filt="default", fill_data_gaps=False,
@@ -697,6 +697,8 @@ def load(self, config_file=None, overwrite_event=True):
                     if val != old_val:
                         logger.debug(f"{key}: {old_val} -> {val}")
                         setattr(self, key, val)
+                else:
+                    self.kwargs[key] = val
 
         # Reset log level based on the config file
         if self.log_level is not None:
@@ -1537,17 +1539,21 @@ def plot(self):
         """
         Plot map and record section if requested
         """
+        show_map = self.kwargs.get("show_map", False)
+        show_rs = self.kwargs.get("show_rs", False)
+
         if "map" in self.plot_files or "all" in self.plot_files:
             logger.info("plotting source receiver map")
             fid = os.path.join(self.output_dir, f"station_map.png")
-            plot_source_receiver_map(self.inv, self.event, save=fid)
+            plot_source_receiver_map(self.inv, self.event, save=fid, 
+                                     show=show_map)
 
         if "record_section" in self.plot_files or "all" in self.plot_files:
             fid = os.path.join(self.output_dir, f"record_section.png")
             # Default settings to create a general record section
             rs = RecordSection(st=self.st, sort_by="distance",
-                               scale_by="normalize", overwrite=True, show=False,
-                               save=fid)
+                               scale_by="normalize", overwrite=True, 
+                               show=show_rs, save=fid)
             rs.run()
 
     def _event_tag_and_output_dir(self):
@@ -1602,6 +1608,8 @@ def run(self, event=None, inv=None, st=None, **kwargs):
         :param st: optional user-provided strean object which will force a
             skip over waveform searching
         """
+        logger.debug(f"running PySEP version {__version__}")
+
         # Overload default parameters with event input file and check validity
         self.load(**kwargs)
         self.check()

pysep/utils/plot.py

@@ -200,21 +200,59 @@ def plot_source_receiver_map(inv, event, subset=None, save="./station_map.png",
     # Do some post-plot editing of the figure
     ax = plt.gca()
 
-    # Get the extent of stations in the Inventory to resize the figure in local
-    # since the original bounds will be just around the event
+    # Get the extent all points plotted to resize the figure in local
+    # since the original bounds will be confined to a small reigon around the
+    # event
     if projection == "local":
-        sta_lats, sta_lons = [], []
+        lats, lons = [], []
         for net in inv:
             for sta in net:
-                sta_lats.append(sta.latitude)
-                sta_lons.append(sta.longitude)
-
-        # Add a small buffer around stations so they're not directly on the edge
-        min_lat = min(sta_lats) * .95
-        max_lat = max(sta_lats) * 1.05
-        min_lon = min(sta_lons) * .95
-        max_lon = max(sta_lons) * 1.05
-        ax.set_extent([min_lon, max_lon, min_lat, max_lat])
+                lats.append(sta.latitude)
+                lons.append(sta.longitude)
+        lats.append(event.preferred_origin().latitude)
+        lons.append(event.preferred_origin().longitude)
+        lats = np.array(lats)
+        lons = np.array(lons)
+
+        # Copied verbatim from obspy.imaging.maps.plot_map(), find aspect ratio
+        # based on the points plotted. Need to re-do this operation since we
+        # have both event and stations
+        if min(lons) < -150 and max(lons) > 150:
+            max_lons = max(np.array(lons) % 360)
+            min_lons = min(np.array(lons) % 360)
+        else:
+            max_lons = max(lons)
+            min_lons = min(lons)
+        lat_0 = max(lats) / 2. + min(lats) / 2.
+        lon_0 = max_lons / 2. + min_lons / 2.
+        if lon_0 > 180:
+            lon_0 -= 360
+        deg2m_lat = 2 * np.pi * 6371 * 1000 / 360
+        deg2m_lon = deg2m_lat * np.cos(lat_0 / 180 * np.pi)
+        if len(lats) > 1:
+            height = (max(lats) - min(lats)) * deg2m_lat
+            width = (max_lons - min_lons) * deg2m_lon
+            margin = 0.2 * (width + height)
+            height += margin
+            width += margin
+        else:
+            height = 2.0 * deg2m_lat
+            width = 5.0 * deg2m_lon
+        # Do intelligent aspect calculation for local projection to adjust to
+        # figure dimensions
+        w, h = plt.gcf().get_size_inches()
+        aspect = w / h
+
+        if width / height < aspect:
+            width = height * aspect
+        else:
+            height = width / aspect
+
+        # We are ASSUMING that the event is located at the center of the figure
+        # (0, 0), which it should be since we plotted it first
+        x0, y0 = 0, 0  # modified from ObsPy function
+        ax.set_xlim(x0 - width / 2, x0 + width / 2)
+        ax.set_ylim(y0 - height / 2, y0 + height / 2)
 
     # Hijack the ObsPy plot and make some adjustments to make it look nicer
     gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False,