Merge pull request #88 from megies/speed_up_magnitude_calculation

speed up updating of network magnitude

Author: megies

CHANGELOG.txt

@@ -19,6 +19,9 @@ master
    temporary directory for easy reuse (see #79)
  - taper settings (applied before any filtering) can now be controlled via
    config file (see #80)
+ - speed up magnitude calculation by reusing unchanged station magnitudes when
+   updating the network magnitude while changing magnitude picks on one station
+   (see #88)
 
 0.5.1
  - fix getting metadata via arclink (see #65)

obspyck/event_helper.py

@@ -74,6 +74,11 @@ def __init__(self, *args, **kwargs):
 
 
 class StationMagnitude(obspy.core.event.StationMagnitude, CommonEventHelper):
+    # we use _amplitude_ids to store amplitude IDS of used Amplitudes, since
+    # QuakeML only has room to store a single Amplitude for a StationMagnitude
+    # unfortunately
+    do_not_warn_on = ['extra', '_amplitude_ids']
+
     def __init__(self, *args, **kwargs):
         super(StationMagnitude, self).__init__()
         self.used = True
@@ -168,6 +173,11 @@ def __init__(self, origin=None, pick=None, *args, **kwargs):
 
 
 class Amplitude(obspy.core.event.Amplitude, CommonEventHelper):
+    # we use _station_magnitude_id to store ID of derived StationMagnitude,
+    # since QuakeML does not link back and this makes it easier to find the
+    # correct one
+    do_not_warn_on = ['extra', '_station_magnitude_id']
+
     def __init__(self, seed_string=None, *args, **kwargs):
         super(Amplitude, self).__init__()
         if seed_string:

obspyck/obspyck.py

@@ -2998,7 +2998,7 @@ def hypoDist(self, coords):
     def calculateStationMagnitudes(self):
         event = self.catalog[0]
         origin = event.origins[0]
-        event.station_magnitudes = []
+        station_magnitudes = []
 
         netstaloc = set([(amp.waveform_id.network_code,
                           amp.waveform_id.station_code,
@@ -3011,7 +3011,22 @@ def calculateStationMagnitudes(self):
             pazs = []
             channels = []
             p2ps = []
-            for amplitude in self.getAmplitudes(net, sta, loc):
+            amplitudes_ = self.getAmplitudes(net, sta, loc)
+
+            # check if we can reuse an existing StationMagnitude
+            amp_ids = set(str(amp.resource_id) for amp in amplitudes_)
+            for sm in event.station_magnitudes:
+                if sm.get('_amplitude_ids') == amp_ids:
+                    break
+            else:
+                sm = None
+            if sm is not None:
+                station_magnitudes.append(sm)
+                self.critical('reusing existing magnitude for %s: %0.2f' % (
+                    sta, sm.mag))
+                continue
+
+            for amplitude in amplitudes_:
                 self.debug(str(amplitude))
                 timedelta = amplitude.get_timedelta()
                 self.debug("Timedelta: " + str(timedelta))
@@ -3046,7 +3061,10 @@ def calculateStationMagnitudes(self):
             dist = self.hypoDist(tr.stats.coordinates)
             mag = estimate_magnitude(pazs, p2ps, timedeltas, dist)
             sm = StationMagnitude()
-            event.station_magnitudes.append(sm)
+            # store ids of used amplitudes, so we're able to tell when a
+            # station magnitude needs to be updated
+            sm._amplitude_ids = amp_ids
+            station_magnitudes.append(sm)
             sm.origin_id = origin.resource_id
             sm.method_id = "/".join(
                 [ID_ROOT, "station_magnitude_method", "obspyck", "2"])
@@ -3062,8 +3080,16 @@ def calculateStationMagnitudes(self):
             extra.amplitudeIDs = {'value': ",".join([str(a.resource_id)
                                                      for a in amplitudes]),
                                   'namespace': NAMESPACE}
+            # need to do this at the very end when StationMagnitude is not
+            # changed anymore, and we need to avoid setting a new random ID
+            # automatically (in event_helper.py)
+            for amplitude in amplitudes:
+                super(Amplitude, amplitude).__setattr__(
+                    '_station_magnitude_id', str(sm.resource_id))
             self.critical('calculated new magnitude for %s: %0.2f (channels: %s)' % (
                 sta, mag, ",".join(channels)))
+        # finally switch out station magnitudes
+        event.station_magnitudes = station_magnitudes
 
     #see http://www.scipy.org/Cookbook/LinearRegression for alternative routine
     #XXX replace with drawWadati()