Updated Clustering Routines (#249)

* Clustering is now being performed on observed travel-times, corrected
  for ellipticity and elevation
* Residuals are now recomputed in SSST_Results

seismic/pick_harvester/cluster_arrivals.py

@@ -81,15 +81,15 @@ def cluster_helper(c_imask):
 
         result = []
         for k, rows in tqdm(cdict.items()):
-            ott = sr.ott[rows[:, 2]]
+            ctt = sr.corrected_travel_time[rows[:, 2]]
 
-            med_idx = np.argwhere(ott == np.quantile(ott, 0.5, interpolation='nearest'))[0][0]
+            med_idx = np.argwhere(ctt == np.quantile(ctt, 0.5, interpolation='nearest'))[0][0]
 
             src_block, sta_block, g_med_idx = rows[med_idx, :]
             tup = (src_block, sta_block, sr.residual[g_med_idx], sr.eorigin_ts[g_med_idx],
                    sr.elons[g_med_idx], sr.elats[g_med_idx], sr.edepths_km[g_med_idx],
                    sr.slons[g_med_idx], sr.slats[g_med_idx], sr.selevs_km[g_med_idx],
-                   (sr.arrivals['arrival_ts'][g_med_idx] - sr.eorigin_ts[g_med_idx] - sr.tcorr[g_med_idx]),
+                   sr.corrected_travel_time[g_med_idx],
                    sr.ecdists[g_med_idx], sr.arrivals['phase'][g_med_idx],
                    1 if sr.is_P[g_med_idx] else 2)
 
@@ -126,7 +126,7 @@ def cluster_helper(c_imask):
     # 1. arrival is part of an event marked as being of
     #    good quality
     # 2. its residual lies within cutoff value
-    # 3. if an automatic arrival, it should meed the
+    # 3. if an automatic arrival, it should meet the
     #    slope-ratio criteria
     # 4. its phase must belong to the set of phases being
     #    used for clustering

seismic/pick_harvester/parametric_data.py

@@ -24,6 +24,9 @@
 from seismic.pick_harvester.utils import split_list
 import h5py
 
+EARTH_RADIUS_KM = 6371.
+DEG2KM = np.pi / 180 * EARTH_RADIUS_KM
+
 class ParametricData:
     def __init__(self, csv_catalog, auto_pick_files=[], auto_pick_phases=[],
                  events_only=False, phase_list='P Pg Pb Pn S Sg Sb Sn', temp_dir='./'):
@@ -41,8 +44,6 @@ def __init__(self, csv_catalog, auto_pick_files=[], auto_pick_phases=[],
                          Note that this temporary folder must be accessible by all MPI ranks, e.g.
                          within a project folder on the NCI.
         """
-        self.EARTH_RADIUS_KM = 6371.
-        self.DEG2KM = np.pi/180 * self.EARTH_RADIUS_KM
         self.STATION_DIST_M = 1e3 # distance in metres within which neighbouring stations are coalesced
         self.csv_catalog = csv_catalog
         self.auto_pick_files = auto_pick_files
@@ -192,7 +193,7 @@ def _coalesce_network_codes(self):
                         lon2, lat2 = coordsdict[net2 + b'.' + sta]
 
                         _, _, dist = self._geod.inv(lon1, lat1, lon2, lat2)
-                        dist *= self.DEG2KM * 1e3 #m
+                        dist *= DEG2KM * 1e3 #m
                         if(dist < self.STATION_DIST_M):
                             swap_map[(net1, sta)] = net2
                         else:
@@ -217,7 +218,7 @@ def _coalesce_network_codes(self):
             iter_count += 1
             if(len(swap_map) == 0): break
         # wend
-        if(self.rank==0): np.save('coalesced_net.npy', self.arrivals['net'])
+        #if(self.rank==0): np.save('coalesced_net.npy', self.arrivals['net'])
     # end func
 
     def _make_temp_dir(self):

seismic/pick_harvester/ssst_relocator.py

@@ -19,11 +19,14 @@
 from collections import defaultdict
 from seismic.pick_harvester.utils import split_list
 
-from seismic.pick_harvester.parametric_data import ParametricData
+from seismic.pick_harvester.parametric_data import ParametricData, EARTH_RADIUS_KM
 import os
 from tqdm import tqdm
 from seismic.pick_harvester.travel_time import TTInterpolator
 from seismic.pick_harvester.ellipticity import Ellipticity
+from seismic.pick_harvester.ssst_utils import compute_ellipticity_correction, \
+                                                compute_elevation_correction, \
+                                                VP_SURF, VS_SURF
 import h5py
 import json
 
@@ -88,13 +91,11 @@ def __init__(self, csv_catalog, auto_pick_files=[], auto_pick_phases=[], config_
 
         # initialize surface velocity array (only needs to be computed once, since
         # elevations are fixed and P/S arrivals are not interchanged)
-        vs_surf = 3.46        # Sg velocity km/s for elevation corrections
-        vp_surf = 5.8         # Pg velocity km/s for elevation corrections
         pidx = self.is_P(self.arrivals['phase'])
         sidx = self.is_S(self.arrivals['phase'])
         self.vsurf = np.zeros(len(self.arrivals), dtype='f4')
-        self.vsurf[pidx] = vp_surf
-        self.vsurf[sidx] = vs_surf
+        self.vsurf[pidx] = VP_SURF
+        self.vsurf[sidx] = VS_SURF
     # end func
 
     def is_P(self, phase):
@@ -158,8 +159,8 @@ def _compute_residual_helper(self, phase, imask=None):
 
         azim, _, ecdist = self._geod.inv(elon, elat, slon, slat)
 
-        elev_corr = self.compute_elevation_correction(phase, vsurf, elev_km, ecdist, edepth_km)
-        ellip_corr = self.compute_ellipticity_correction(phase, ecdist, edepth_km, elat, azim)
+        elev_corr = compute_elevation_correction(self._tti, phase, vsurf, elev_km, ecdist, edepth_km)
+        ellip_corr = compute_ellipticity_correction(self._ellipticity, phase, ecdist, edepth_km, elat, azim)
         ptt = self._tti.get_tt(phase, ecdist, edepth_km)
 
         # compute empirical tt
@@ -232,46 +233,6 @@ def redefine_phases(self, imask=None):
         self.residual = self._sync_var(local_new_residual)
     # end func
 
-    def compute_elevation_correction(self, phase, vsurf, elev_km, ecdist, edepth_km):
-        """
-        Operates on input parameters irrespective of local/global indices
-
-        @param phase: dim(any)
-        @param vsurf:  dim(any)
-        @param elev_km:  dim(any)
-        @param ecdist:  dim(any)
-        @param edepth_km:  dim(any)
-        @return:
-        """
-        dtdd = self._tti.get_dtdd(phase, ecdist, edepth_km)
-        elev_corr = vsurf * (dtdd / self.DEG2KM)
-        elev_corr = np.power(elev_corr, 2.)
-        elev_corr[elev_corr > 1.] = 1./elev_corr[elev_corr > 1.]
-        elev_corr = np.sqrt(1. - elev_corr)
-        elev_corr = elev_corr * elev_km / vsurf
-
-        elev_corr = np.array(elev_corr.astype('f4'))
-
-        return elev_corr
-    # end func
-
-    def compute_ellipticity_correction(self, phase, ecdist, edepth_km, elat, azim):
-        """
-        Operates on input parameters irrespective of local/global indices
-
-        @param phase:  dim(any)
-        @param ecdist:  dim(any)
-        @param edepth_km:  dim(any)
-        @param elat:  dim(any)
-        @param azim:  dim(any)
-        @return:
-        """
-        ellip_corr = self._ellipticity.get_correction(phase, ecdist, edepth_km, elat, azim)
-        ellip_corr = np.array(ellip_corr.astype('f4'))
-
-        return ellip_corr
-    # end func
-
     def compute_sst_correction(self, min_slope_ratio=5):
         """
         Operates on global data
@@ -466,9 +427,10 @@ def relocate_events(self, imask=None, reloc_niter=10, reloc_dlon=0.25,
                             edepth_km = edepth_km0 + idep * ddep
                             azim, _, ecdist = self._geod.inv(elon * ones, elat * ones, slon, slat)
 
-                            elev_corr = self.compute_elevation_correction(phase, vsurf, elev_km, ecdist, edepth_km * ones)
-                            ellip_corr = self.compute_ellipticity_correction(phase, ecdist, edepth_km * ones, elat * ones,
-                                                                           azim)
+                            elev_corr = compute_elevation_correction(self._tti, phase, vsurf,
+                                                                     elev_km, ecdist, edepth_km * ones)
+                            ellip_corr = compute_ellipticity_correction(self._ellipticity, phase, ecdist,
+                                                                        edepth_km * ones, elat * ones, azim)
                             tt = self._tti.get_tt(phase, ecdist, edepth_km * ones)
                             tt = np.ma.masked_array(tt, mask=tt == self._tti.fill_value)
 
@@ -618,7 +580,7 @@ def dump_h5(fn, iter):
     def _initialize_spatial_functors(self, ellipsoidal_distance=False):
         if(self.rank == 0): print('Initializing spatial functors..')
 
-        ER = self.EARTH_RADIUS_KM * 1e3 #m
+        ER = EARTH_RADIUS_KM * 1e3 #m
 
         if(ellipsoidal_distance):
             transformer = pyproj.Transformer.from_crs(

seismic/pick_harvester/ssst_utils.py

@@ -3,8 +3,56 @@
 from collections import defaultdict
 import h5py
 from seismic.pick_harvester.ellipticity import Ellipticity
+from seismic.pick_harvester.travel_time import TTInterpolator
+from seismic.pick_harvester.parametric_data import EARTH_RADIUS_KM, DEG2KM
 from pyproj import Geod
 
+VS_SURF = 3.46  # Sg velocity km/s for elevation corrections
+VP_SURF = 5.8  # Pg velocity km/s for elevation corrections
+
+def compute_elevation_correction(tti:TTInterpolator, phase, vsurf, elev_km, ecdist, edepth_km):
+    """
+    Computes elevation corrections
+
+    @param tti: travel-time interpolator
+    @param phase: dim(any)
+    @param vsurf:  dim(any)
+    @param elev_km:  dim(any)
+    @param ecdist:  dim(any)
+    @param edepth_km:  dim(any)
+    @return:
+    """
+
+    dtdd = tti.get_dtdd(phase, ecdist, edepth_km)
+    elev_corr = vsurf * (dtdd / DEG2KM)
+    elev_corr = np.power(elev_corr, 2.)
+    elev_corr[elev_corr > 1.] = 1./elev_corr[elev_corr > 1.]
+    elev_corr = np.sqrt(1. - elev_corr)
+    elev_corr = elev_corr * elev_km / vsurf
+
+    elev_corr = np.array(elev_corr.astype('f4'))
+
+    return elev_corr
+# end func
+
+def compute_ellipticity_correction(ellipticity:Ellipticity, phase, ecdist, edepth_km, elat, azim):
+    """
+    Computes ellipticity corrections
+
+    @param ellipticity
+    @param phase:  dim(any)
+    @param ecdist:  dim(any)
+    @param edepth_km:  dim(any)
+    @param elat:  dim(any)
+    @param azim:  dim(any)
+    @return:
+    """
+    ellip_corr = ellipticity.get_correction(phase, ecdist, edepth_km, elat, azim)
+    ellip_corr = np.array(ellip_corr.astype('f4'))
+
+    return ellip_corr
+# end func
+
 def h5_to_named_array(hfn, key):
     h = h5py.File(hfn, 'r')
 
@@ -55,8 +103,6 @@ def __init__(self, h5_fn, iteration=-1):
         self.arrivals = h5_to_named_array(self.h5_fn, '{}/arrivals'.format(self.iter))
         self.events = h5_to_named_array(self.h5_fn, '{}/events'.format(self.iter))
         self.event_id_to_idx = h5_to_named_array(self.h5_fn, '0/events/event_id_to_idx')
-        self.residual = h5_to_named_array(self.h5_fn, '{}/arrivals/residual'.format(self.iter))
-        self.tcorr = h5_to_named_array(self.h5_fn, '{}/arrivals/tcorr'.format(self.iter))
         self.is_AUTO_arrival = h5_to_named_array(self.h5_fn, '0/arrivals/is_AUTO_arrival')
 
         self.elons = self.events['lon'][self.event_id_to_idx[self.arrivals['event_id']]]
@@ -78,18 +124,27 @@ def __init__(self, h5_fn, iteration=-1):
         self.is_S = self.phase.astype('S1') == b'S'
         self.slope_ratio = self.arrivals['quality_measure_slope']
 
-        # ===========================================================
+        self._ellipticity = Ellipticity()
+        self._tti = TTInterpolator()
+
         # Compute ellipticity-correction
-        # ===========================================================
-        self.ellipticity = Ellipticity()
+        self.ellip_corr = compute_ellipticity_correction(self._ellipticity, self.phase, self.ecdists,
+                                                         self.edepths_km, self.elats, self.azims)
+
+        # Compute elevation-correction
+        self.vsurf = np.zeros(len(self.arrivals), dtype='f4')
+        self.vsurf[self.is_P] = VP_SURF
+        self.vsurf[self.is_S] = VS_SURF
 
-        self.ellip_corr = self.ellipticity.get_correction(self.phase, self.ecdists,
-                                                           self.edepths_km, self.elats, self.azims)
+        self.elev_corr = compute_elevation_correction(self._tti, self.phase, self.vsurf,
+                                                      self.selevs_km, self.ecdists, self.edepths_km)
 
-        # ===========================================================
         # Compute observed travel-times
-        # ===========================================================
-        self.ott = self.arrival_ts - self.eorigin_ts - self.tcorr + self.ellip_corr
+        self.corrected_travel_time = self.arrival_ts - self.eorigin_ts - self.ellip_corr - self.elev_corr
+
+        # Compute residual
+        ptt = self._tti.get_tt(self.phase, self.ecdists, self.edepths_km)
+        self.residual = self.corrected_travel_time - ptt
     # end func
 
     def find_paired_S(self):