Minor changes to refelect updates to correlator routines

Author: geojunky

seismic/xcorqc/xcorr_station_clock_analysis.py

@@ -29,6 +29,7 @@
 from sklearn.cluster import dbscan
 from scipy.interpolate import LSQUnivariateSpline
 from pandas.plotting import register_matplotlib_converters, deregister_matplotlib_converters
+from seismic.misc import print_exception
 
 # import obspy
 from netCDF4 import Dataset as NCDataset
@@ -37,7 +38,6 @@
 from seismic.ASDFdatabase import FederatedASDFDataSet
 from seismic.xcorqc.analytic_plot_utils import distance, timestamps_to_plottable_datetimes
 
-
 class XcorrClockAnalyzer:
     """
     Helper class for bundling of preprocessed cross-correlation data before plotting
@@ -189,8 +189,8 @@ def _preprocess(self):
         xc_start_times = xcdata.variables['IntervalStartTimes'][:]
         # xc_end_times = xcdata.variables['IntervalEndTimes'][:]
         xc_lag = xcdata.variables['lag'][:]
-        xc_xcorr = xcdata.variables['xcorr'][:, :]
-        xc_num_stacked_windows = xcdata.variables['NumStackedWindows'][:]
+        xc_xcorr = xcdata.variables['X'][:, :]
+        xc_num_stacked_windows = xcdata.variables['nX'][:]
         xcdata.close()
         xcdata = None
 
@@ -389,12 +389,12 @@ def plot_resampled_clusters(self, ax, ids, resampled_corrections, stn_code=''):
 # end class
 
 
-def station_codes(filename):
+def nsl_codes(filename):
     """
     Convert a netCDF4 .nc filename generated by correlator to the corresponding
     station codes in the format ``NETWORK.STATION``
 
-    Assumed format: ``NET1.STAT1.NET2.STA2.*.nc``
+    Assumed format: ``NET1.STAT1.LOC1.CHA1.NET2.STA2.LOC2.CHA2.nc``
 
     :param filename: The ``.nc`` filename from which to extract the station and network codes
     :type filename: str
@@ -403,34 +403,35 @@ def station_codes(filename):
     """
     _, fname = os.path.split(filename)
     parts = fname.split('.')
-    sta1 = '.'.join(parts[0:2])
-    sta2 = '.'.join(parts[2:4])
-    return (sta1, sta2)
+    nsl1 = '.'.join(parts[0:3])
+    nsl2 = '.'.join(parts[4:7])
+    return (nsl1, nsl2)
 
 
-def station_distance(federated_ds, code1, code2):
+def station_distance(federated_ds, nsl1, nsl2):
     """
     Determine the distance in km between a pair of stations at a given time.
 
     :param federated_ds: Federated dataset to query for station coordinates
     :type federated_ds: seismic.ASDFdatabase.FederatedASDFDataSet
-    :param code1: Station and network code in the format ``NETWORK.STATION`` for first station
-    :type code1: str
-    :param code2: Station and network code in the format ``NETWORK.STATION`` for second station
-    :type code2: str
+    :param nsl1: ``NETWORK.STATION.LOCATION`` for first station
+    :type nsl1: str
+    :param nsl2: ``NETWORK.STATION.LOCATION`` for second station
+    :type nsl2: str
     :return: Distance between stations in kilometers
     :rtype: float
     """
-    coords1 = federated_ds.unique_coordinates[code1]
-    coords2 = federated_ds.unique_coordinates[code2]
+
+    ns1 = '.'.join(nsl1.split('.')[:2])
+    ns2 = '.'.join(nsl2.split('.')[:2])
+    coords1 = federated_ds.unique_coordinates[ns1]
+    coords2 = federated_ds.unique_coordinates[ns2]
     return distance(coords1, coords2)
 
 
 def plot_xcorr_time_series(ax, x_lag, y_times, xcorr_data, use_formatter=False):
-
     np_times = np.array([datetime.datetime.utcfromtimestamp(v) for v in y_times]).astype('datetime64[s]')
-    gx, gy = np.meshgrid(x_lag, np_times)
-    im = ax.pcolormesh(gx, gy, xcorr_data, cmap='RdYlBu_r', vmin=0, vmax=1, rasterized=True)
+    im = ax.pcolormesh(x_lag, np_times, xcorr_data, cmap='RdYlBu_r', vmin=0, vmax=1, rasterized=True)
 
     if use_formatter:
         date_formatter = matplotlib.dates.DateFormatter("%Y-%m-%d")
@@ -515,14 +516,13 @@ def plot_estimated_timeshift(ax, x_lag, y_times, correction, annotation=None, ro
     if row_rcf_crosscorr is not None:
         # Line plot laid over the top of RCF * CCF
         np_times = np.array([datetime.datetime.utcfromtimestamp(v) for v in y_times]).astype('datetime64[s]')
-        gx, gy = np.meshgrid(x_lag, np_times)
         plot_data = row_rcf_crosscorr
         crange_floor = 0.7
         plot_data[np.isnan(plot_data)] = crange_floor
         plot_data[(plot_data < crange_floor)] = crange_floor
         nan_mask = np.isnan(correction)
         plot_data[nan_mask, :] = np.nan
-        ax.pcolormesh(gx, gy, plot_data, cmap='RdYlBu_r', rasterized=True)
+        ax.pcolormesh(x_lag, np_times, plot_data, cmap='RdYlBu_r', rasterized=True)
         ax.set_ylim((min(np_times), max(np_times)))
         xrange = 1.2 * np.nanmax(np.abs(correction))
         ax.set_xlim((-xrange, xrange))
@@ -561,15 +561,15 @@ def plot_xcorr_file_clock_analysis(src_file, asdf_dataset, time_window, snr_thre
     rcf_corrected = xcorr_ca.rcf_corrected
     row_rcf_crosscorr = xcorr_ca.row_rcf_crosscorr
 
-    origin_code, dest_code = station_codes(src_file)
-    dist = station_distance(asdf_dataset, origin_code, dest_code)
+    nsl1, nsl2 = nsl_codes(src_file)
+    dist = station_distance(asdf_dataset, nsl1, nsl2)
 
     # -----------------------------------------------------------------------
     # Master layout and plotting code
 
     fig = plt.figure(figsize=(11.69, 16.53))
-    fig.suptitle("Station: {}, Dist. to {}: {:3.2f} km{}".format(
-        origin_code, dest_code, dist, '' if not title_tag else ' ' + title_tag), fontsize=16, y=1)
+    fig.suptitle("Station-pair: {}-{}, Separation: {:3.2f} km{}".format(
+        nsl1, nsl2, dist, '' if not title_tag else ' ' + title_tag), fontsize=16, y=1)
 
     ax1 = fig.add_axes([0.1, 0.075, 0.5, 0.725])
 
@@ -623,13 +623,9 @@ def setstr(name):
         pdf_out = PdfPages(pdf_file)
         pdf_out.savefig(plt.gcf(), dpi=600)
         pdf_out.close()
-        # Change permissions so that others in same group can read and write over this file.
-        os.chmod(pdf_file, (stat.S_IRGRP | stat.S_IWGRP))
 
     if png_file is not None:
         plt.savefig(png_file, dpi=150)
-        # Change permissions so that others in same group can read and write over this file.
-        os.chmod(png_file, (stat.S_IRGRP | stat.S_IWGRP))
 
     if show:
         plt.show()
@@ -663,8 +659,8 @@ def read_correlator_config(nc_file):
     folder, fname = os.path.split(nc_file)
     base, _ = os.path.splitext(fname)
     base_parts = base.split('.')
-    if len(base_parts) > 4:
-        timestamp = '.'.join(base_parts[4:])
+    if len(base_parts) > 8:
+        timestamp = '.'.join(base_parts[8:])
         config_filename = '.'.join(['correlator', timestamp, 'cfg'])
         config_file = os.path.join(folder, config_filename)
         if os.path.exists(config_file):
@@ -732,6 +728,7 @@ def batch_process_xcorr(src_files, dataset, time_window, snr_threshold, pearson_
             if save_plots:
                 basename, _ = os.path.splitext(src_file)
                 png_file = basename + ".png"
+                pdf_file = basename + ".pdf"
                 # If png file already exists and has later timestamp than src_file, then skip it.
                 if os.path.exists(png_file):
                     src_file_time = os.path.getmtime(src_file)
@@ -745,7 +742,8 @@ def batch_process_xcorr(src_files, dataset, time_window, snr_threshold, pearson_
                         pbar.update()
                         continue
                 plot_xcorr_file_clock_analysis(src_file, dataset, time_window, snr_threshold, pearson_cutoff_factor,
-                                               png_file=png_file, show=False, underlay_rcf_xcorr=underlay_rcf_xcorr,
+                                               pdf_file=pdf_file, show=False,
+                                               underlay_rcf_xcorr=underlay_rcf_xcorr,
                                                title_tag=title_tag, settings=settings)
             else:
                 plot_xcorr_file_clock_analysis(src_file, dataset, time_window, snr_threshold, pearson_cutoff_factor,
@@ -755,6 +753,7 @@ def batch_process_xcorr(src_files, dataset, time_window, snr_threshold, pearson_
             pbar.update()
 
         except Exception as e:
+            print_exception(e)
             tqdm.write("ERROR processing file {}".format(src_file))
             failed_files.append((src_file, str(e)))