diff --git a/pyproject.toml b/pyproject.toml index 66995ba..41e3269 100644 --- a/pyproject.toml +++ b/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" diff --git a/pysep/__init__.py b/pysep/__init__.py index b59a3e3..94df49a 100644 --- a/pysep/__init__.py +++ b/pysep/__init__.py @@ -1,5 +1,7 @@ import logging +__version__ = "0.3.2" + logger = logging.getLogger("pysep") logger.setLevel("INFO") logger.propagate = 0 diff --git a/pysep/pysep.py b/pysep/pysep.py index 86a1483..88df5d8 100644 --- a/pysep/pysep.py +++ b/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() diff --git a/pysep/utils/plot.py b/pysep/utils/plot.py index 7e272d4..08bcddc 100644 --- a/pysep/utils/plot.py +++ b/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,