Add function to compute distance between one event and a catalog

eqcorrscan/tests/clustering_test.py

@@ -18,7 +18,7 @@
 from eqcorrscan.utils.clustering import (
     cross_chan_correlation, distance_matrix, cluster, group_delays, svd,
     empirical_svd, svd_to_stream, corr_cluster, dist_mat_km, catalog_cluster,
-    space_time_cluster, remove_unclustered)
+    space_time_cluster, remove_unclustered, dist_array_km)
 from eqcorrscan.helpers.mock_logger import MockLoggingHandler
 
 
@@ -381,6 +381,12 @@ def test_dist_mat_km(self):
                         (slave_ori.latitude, slave_ori.longitude,
                          slave_ori.depth / 1000)), 6)
 
+    def test_dist_array_km(self):
+        """ Test that just computing one array works as expected """
+        dist_array = dist_array_km(master=self.cat[0], catalog=self.cat)
+        self.assertEqual(len(dist_array), len(self.cat))
+        self.assertEqual(dist_array[0], 0)
+
     def test_space_cluster(self):
         """Test the wrapper around dist_mat_km."""
         groups = catalog_cluster(

eqcorrscan/utils/catalog_to_dd.py

@@ -28,16 +28,60 @@
 
 # Some hypoDD specific event holders - classes were faster than named-tuples
 
+class SparseOrigin(object):
+    def __init__(self, resource_id, latitude, longitude, depth, time):
+        self.resource_id = resource_id
+        self.latitude = latitude
+        self.longitude = longitude
+        self.depth = depth
+        self.time = time
+
+    def __repr__(self):
+        return (f"SparseOrigin(resource_id={self.resource_id}, "
+                f"latitude={self.latitude}, longitude={self.longitude}, "
+                f"depth={self.depth}, time={self.time})")
+
+    @classmethod
+    def from_origin(cls, origin: Origin):
+        return cls(resource_id=origin.resource_id.id,
+                   latitude=origin.latitude,
+                   longitude=origin.longitude,
+                   depth=origin.depth,
+                   time=origin.time)
+
+
 class SparseEvent(object):
-    def __init__(self, resource_id, picks, origin_time):
+    def __init__(self, resource_id, picks, origin):
         self.resource_id = resource_id
         self.picks = picks
-        self.origin_time = origin_time
+        self.origin = origin
 
     def __repr__(self):
-        return ("SparseEvent(resource_id={0}, origin_time={1}, picks=[{2} "
+        return ("SparseEvent(resource_id={0}, origin={1}, picks=[{2} "
                 "picks])".format(
-                    self.resource_id, self.origin_time, len(self.picks)))
+                    self.resource_id, self.origin, len(self.picks)))
+
+    @property
+    def origin_time(self):
+        return self.origin.time
+
+    def preferred_origin(self):
+        return self.origin
+
+    @classmethod
+    def from_event(cls, event: Event):
+        origin = SparseOrigin.from_origin(
+            event.preferred_origin() or event.origins[0])
+        time_weight_dict = {
+            arr.pick_id: arr.time_weight or 1.0 for origin in event.origins
+            for arr in origin.arrivals}
+        picks = [
+            SparsePick.from_pick(
+                pick=pick, origin_time=origin.time,
+                time_weight=time_weight_dict.get(pick.resource_id, 1.0)
+            ) for pick in event.picks]
+        return cls(resource_id=event.resource_id.id, origin=origin,
+                   picks=picks)
 
 
 class SparsePick(object):
@@ -63,6 +107,29 @@
     def channel(self):
         return self.seed_id.split('.')[-1]
 
+    @property
+    def network(self):
+        return self.seed_id.split('.')[0]
+
+    @property
+    def location(self):
+        return self.seed_id.split('.')[2]
+
+    @classmethod
+    def from_pick(
+        cls,
+        pick: Pick,
+        origin_time: UTCDateTime,
+        time_weight: float = 1.0
+    ):
+        return cls(
+            tt=pick.time - origin_time,
+            time=pick.time,
+            seed_id=pick.waveform_id.get_seed_string(),
+            phase_hint=pick.phase_hint[0],  # Note, only using P or S hints
+            time_weight=time_weight,
+            waveform_id=pick.waveform_id)
+
 
 # Generic helpers
 
@@ -97,7 +164,7 @@
 
 
 class _EventPair(object):
-    """ Holder for event paid observations. """
+    """ Holder for event pair observations. """
 
     def __init__(self, event_id_1, event_id_2, obs=None):
         self.event_id_1 = event_id_1
@@ -146,22 +213,7 @@
 
 def _make_sparse_event(event):
     """ Make a sparse event with just the info hypoDD needs. """
-    origin_time = (event.preferred_origin() or event.origins[0]).time
-    time_weight_dict = {
-        arr.pick_id: arr.time_weight or 1.0 for origin in event.origins
-        for arr in origin.arrivals}
-    sparse_event = SparseEvent(
-        resource_id=event.resource_id.id,
-        origin_time=origin_time,
-        picks=[SparsePick(
-            tt=pick.time - origin_time,
-            time=pick.time,
-            seed_id=pick.waveform_id.get_seed_string(),
-            phase_hint=pick.phase_hint[0],  # Only use P or S hints.
-            time_weight=time_weight_dict.get(pick.resource_id, 1.0),
-            waveform_id=pick.waveform_id)
-            for pick in event.picks])
-    return sparse_event
+    return SparseEvent.from_event(event=event)
 
 
 def _prepare_stream(stream, event, extract_len, pre_pick, seed_pick_ids=None):

eqcorrscan/utils/clustering.py

@@ -11,6 +11,7 @@
 import os
 import logging
 import ctypes
+from math import radians
 from multiprocessing import cpu_count
 
 import matplotlib.pyplot as plt
@@ -963,6 +964,65 @@ def remove_unclustered(catalog, distance_cutoff, num_threads=None):
     return catalog
 
 
+def dist_array_km(master, catalog, num_threads=None):
+    """
+    Compute an array of distances relative to one core event.
+
+    Will give phyiscal hypocentral distances in km.
+
+    :type master: obspy.core.event.Event
+    :param master: Master event to compute distances relative to
+    :type catalog: obspy.core.event.Catalog
+    :param catalog: Catalog for which to compute the distance matrix
+
+    :returns: distance array
+    :rtype: :class:`numpy.ndarray`
+    """
+    utilslib = _load_cdll('libutils')
+
+    utilslib.distance_array.argtypes = [
+        ctypes.c_float, ctypes.c_float, ctypes.c_float,
+        np.ctypeslib.ndpointer(dtype=np.float32,
+                               flags='C_CONTIGUOUS'),
+        np.ctypeslib.ndpointer(dtype=np.float32,
+                               flags='C_CONTIGUOUS'),
+        np.ctypeslib.ndpointer(dtype=np.float32,
+                               flags='C_CONTIGUOUS'),
+        ctypes.c_long,
+        np.ctypeslib.ndpointer(dtype=np.float32,
+                               flags='C_CONTIGUOUS'),
+        ctypes.c_int]
+    utilslib.distance_array.restype = ctypes.c_int
+
+    # Initialize empty array
+    dist_array = np.zeros(len(catalog), dtype=np.float32)
+    latitudes, longitudes, depths = (
+        np.empty(len(catalog)), np.empty(len(catalog)), np.empty(len(catalog)))
+    for i, event in enumerate(catalog):
+        origin = event.preferred_origin() or event.origins[0]
+        latitudes[i] = origin.latitude
+        longitudes[i] = origin.longitude
+        depths[i] = origin.depth
+    depths /= 1000.0
+    depths = np.ascontiguousarray(depths, dtype=np.float32)
+    latitudes = np.ascontiguousarray(np.radians(latitudes), dtype=np.float32)
+    longitudes = np.ascontiguousarray(np.radians(longitudes), dtype=np.float32)
+
+    # Generally this is better single threaded
+    num_threads = num_threads or 1
+
+    master_origin = master.preferred_origin() or master.origins[0]
+
+    ret = utilslib.distance_array(
+        radians(master_origin.latitude), radians(master_origin.longitude),
+        master_origin.depth / 1000., latitudes, longitudes, depths,
+        len(catalog), dist_array, num_threads)
+
+    if ret != 0:  # pragma: no cover
+        raise Exception("Internal error while computing distance matrix")
+    return dist_array
+
+
 def dist_mat_km(catalog, num_threads=None):
     """
     Compute the distance matrix for a catalog using hypocentral separation.

eqcorrscan/utils/src/distance_cluster.c

@@ -32,6 +32,8 @@
 
 float dist_calc(float, float, float, float, float, float);
 
+int distance_array(float, float, float, float*, float*, float*, long, float*, int);
+
 int distance_matrix(float*, float*, float*, long, float*, int);
 
 int remove_unclustered(float*, float*, float*, long, unsigned char*, float, int);
@@ -67,6 +69,32 @@ float dist_calc(float lat1, float lon1, float depth1, float lat2, float lon2, fl
     return distance;
 }
 
+int distance_array(float lat1, float lon1, float depth1, float *latitudes,
+                   float *longitudes, float *depths, long n_locs,
+                   float *dist_array, int n_threads){
+    /* Calculate the distances of all locations relative to one location.
+    *
+    *  :type lat1: Latitude of core location in radians
+    *  :type lon1: Longitude of core location in radians
+    *  :type depth1: Depth of core location in km (positive down)
+    *  :type latitudes: Array of floats of latitudes in radians
+    *  :type longitudes: Array of floats of longitudes in radians
+    *  :type depths: Array of floats of depths in km (positive down)
+    *  :type n_locs: Number of locations
+    *  :type dist_array: Array of floats of for output - should be initialized as zeros, and should be n_locs long
+    */
+    int out = 0;
+    long n;
+
+    #pragma omp parallel for num_threads(n_threads)
+    for (n = 0; n < n_locs; ++n){
+        dist_array[n] = dist_calc(
+            lat1, lon1, depth1, latitudes[n], longitudes[n], depths[n]);
+    }
+    return out;
+}
+
+
 int distance_matrix(float *latitudes, float *longitudes, float *depths, long n_locs,
                     float *dist_mat, int n_threads){
     /* Calculate the distance matrix for a set of locations.

eqcorrscan/utils/src/libutils.def

@@ -13,5 +13,6 @@ EXPORTS
     multi_normxcorr_time
     multi_normxcorr_time_threaded
     dist_calc
+    distance_array
     distance_matrix
     remove_unclustered