speed up updating of network magnitude

..by finding out which already existing station magnitudes can be
reused. It was taking ages with many stations because even
unchanged stations were recalculated when picking magnitudes on one
station.

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()