General updates and bugfixes (#97)

* added kwarg parameters  and  to let user determine if they want to show figures during pysep processing. defaults to false
bugfix: pysep load was not retaining keyword arguments for command line usage of pysep

* added version number to __init__ so that its accessible. bumped version number of 0.3.2 so that it is uniquely identifiable from the master version 0.3.1

* related to #95 display pysep version in debug log

* playing around with aspect ratio on maps

* changing default value of 'remove_clipped' parameter to True to be more restrictive by default, related to #80

* map plotting at local scales now auto scales based on the figure size, using obspy internal logic, which fixes the problem of having very narrow map plots or not showing the event, or not showing all the stations

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,