source: SHX/trunk/SeismicHandler/modules/traces.py @ 1191

Revision 1191, 15.2 KB checked in by klaus, 4 years ago (diff)

changes due to new obspy version

  • Property svn:eol-style set to native
Line 
1# -*- coding: utf-8 -*-
2
3#    This file is part of Seismic Handler eXtended (SHX). For terms of use and
4#    license information please see license.txt and visit
5#    http://www.seismic-handler.org/portal/wiki/Shx/LicenseTerms
6
7import os
8import numpy as np
9from obspy.core.util import AttribDict
10from SeismicHandler.basics import Singleton, timeit, timestore
11from SeismicHandler.basics.messages import log_message, subscribe_event
12# disabled: _set_runtime_var
13from obspy.core.stream import Stream, Trace
14from SeismicHandler.config import Settings, set_dsptrcs, set_tottrcs
15from SeismicHandler.modules.stations import Stations
16from obspy.geodetics import gps2dist_azimuth
17
18
19META_STATUS_EMPTY    = 0
20META_STATUS_LOCATION = 1
21META_STATUS_COMPLETE = 2
22
23DEG_TO_KM = 111.19
24
25class Traces(object):
26    """
27    Supply access for two streams of traces. One is used for plotting, the other
28    holds hidden traces. This will be enhanced later, but for compatibility
29    reasons to existing command procedures it's done that way first.
30    """
31    __metaclass__ = Singleton
32
33    traces = AttribDict()
34
35    def __init__(self, *args, **kwargs):
36        self.traces.visible = Stream()
37        self.traces.hidden = Stream()
38        self.traces.created = Stream()
39        self.traces.overlays = OverlayTraces()
40        self.netdict = {}
41        self.optchan = self.readOptChans()
42
43    def __getattr__(self, name):
44        try:
45            return self.__class__.__dict__["traces"][name]
46        except:
47            raise AttributeError(name)
48
49    #@timeit
50    def addreplace(self, stream, applyGain=False, uniquechannel=False):
51        """
52        Adds or replaces visible streams depending on runtime setting.
53        """
54        # this is always a global switch!
55        if not Settings.swKeeptraces:
56            log_message("debug.traces", "replacing visible traces")
57            self.traces.visible.clear()
58
59        # check if traces contains relevant information
60        # single trace
61        if isinstance(stream, Trace):
62            stream = Stream(traces=[stream])
63        if isinstance(stream, list):
64            stream = Stream(traces=stream)
65
66        for t in stream:
67            if (Settings.swUniquechannel or uniquechannel)\
68                and self.alreadyReadIn(t):
69                #print "dbg: uniquechannel, dropping trace", t.stats.station, t.stats.channel
70                stream.remove(t)
71                continue
72            t.establish()
73            t.invalidate_cache()
74            self.cacheNetcode(t)
75
76            if applyGain:
77                # get gain from inventory and apply it to data
78                gain = None
79                try:
80                    gain = float(Stations()[t.id, t.stats.starttime].gain)
81                except KeyError:
82                    msg = "no meta data stored for station '%s' at %s"
83                    log_message("info.traces", msg % (t.id, t.get_info("start")))
84                # data from trace itself, this will override stored information
85                try:
86                    gain = 1e9/t.stats.paz.sensitivity
87                    msg = "gain information attached to trace: %f"
88                    log_message("info.traces", msg % gain)
89                except:
90                    pass
91
92                if gain is None:
93                    msg = "No gain configured for %s at %s"
94                    log_message("warning.traces", msg % (t.id, t.get_info("start")))
95                    gain = 1.
96
97                t.data *= gain
98                msg = "Applied gain factor of %f to %s"
99                log_message("debug.traces", msg % (gain, t.id))
100           
101            t.set_info( 'METASTATUS', self.getMetaStatus(t) )
102       
103        self.traces.visible += stream
104        self.traces.created.clear()
105        self.traces.created += stream
106        self.updateCounter()
107
108    #@timeit
109    def updateCounter(self):
110        dsptrcs = len(self.traces.visible)
111        hidtrcs = len(self.traces.hidden)
112        log_message("info.traces", "Now %u visible and %u hidden traces" % (
113                                                             dsptrcs, hidtrcs))
114
115        #_set_runtime_var("dsptrcs", dsptrcs)
116        #_set_runtime_var("tottrcs", dsptrcs+hidtrcs)
117        set_dsptrcs( dsptrcs )
118        set_tottrcs( dsptrcs+hidtrcs )
119
120    #@timeit
121    def visible_index(self, trace):
122        """
123        Returns position of trace inside visible traces list.
124        """
125        return self.traces.visible.traces.index(trace)
126   
127    def cacheNetcode( self, trace ):
128        """Store netcodes for known stations. Overwrites possibly existing
129        old entries.
130        """
131        upstat = trace.stats.station.upper()
132        self.netdict[upstat] = trace.stats.network.upper()
133   
134    def getNetcode( self, station ):
135        return self.netdict.get( station.upper(), None )
136   
137    def getOptChan( self, station ):
138        return self.optchan.get( station, 'BH' )
139   
140    def getMetaStatus( self, trace, addgain=False ):
141        "Returns 0, 1 or 2 depending on metadata status."
142        stations = Stations()
143        sname = "%s.%s.%s.%s" % (trace.stats.network, trace.stats.station,
144            trace.stats.location, trace.stats.channel)
145        metatime = trace.stats.starttime \
146            + (trace.stats.endtime-trace.stats.starttime)/2
147        try:
148            r = stations[(sname, metatime)]
149        except:
150            if addgain:
151                return (META_STATUS_EMPTY,None)
152            else:
153                return META_STATUS_EMPTY
154        try:
155            slat = float( r.latitude )
156            slon = float( r.longitude )
157        except:
158            if addgain:
159                return (META_STATUS_EMPTY,None)
160            else:
161                return META_STATUS_EMPTY
162        try:
163            gain = float( r.gain )
164        except:
165            if addgain:
166                return (META_STATUS_LOCATION,None)
167            else:
168                return META_STATUS_LOCATION
169        if gain == 1.0:
170            if addgain:
171                return (META_STATUS_LOCATION,None)
172            else:
173                return META_STATUS_LOCATION
174        if r.poles and r.poles != '[]':
175            if addgain:
176                return (META_STATUS_COMPLETE,gain)
177            else:
178                return META_STATUS_COMPLETE
179        else:
180            if addgain:
181                return (META_STATUS_LOCATION,gain)
182            else:
183                return META_STATUS_LOCATION
184   
185    def alreadyReadIn( self, trc ):
186        """
187        Returns true if this channel is already read in (independent of
188        time window)
189        """
190        sname = "%s.%s.%s.%s" % (trc.stats.network,trc.stats.station,
191            trc.stats.location,trc.stats.channel)
192        for t in self.traces.visible:
193            tname = "%s.%s.%s.%s" % (t.stats.network,t.stats.station,
194                t.stats.location,t.stats.channel)
195            if tname == sname:
196                return True
197        return False
198   
199    def readOptChans( self ):
200        optcname = os.path.join(os.path.expanduser("~"), ".shx", "optchans.txt")
201        if not os.path.exists(optcname):
202            return {}
203        optdict = {}
204        for line in file(optcname):
205            if line.strip() == '' or line.startswith('#'):
206                continue
207            try:
208                station, optchan = line.split()
209            except:
210                continue
211            optdict[station] = chan
212        return optdict
213
214
215class OverlayTraces:
216    """
217    Store info about traces drawn above others.
218    Example:
219        overlaydict = {'A':[1,2],'B':[3,4]}
220        allidx = [1,2,3,4]
221    """
222   
223    def __init__( self ):
224        self.clearOverlays()
225       
226    def numberOfOverlays( self ):
227        return (len(self.allidx) - len(self.overlaydict))
228   
229    def addOverlay( self, idxlist ):
230        if len(idxlist) < 2:
231            print "should not happen: short overlay list"
232            return
233        # Check for already specified indices, no index may occur twice.
234        for i in idxlist:
235            if i in self.allidx:
236                print "trace '%s' already in overlay" % i
237                return
238        self.overlaydict[self.nextname] = idxlist
239        self.nextname = chr(ord(self.nextname)+1)
240        self.allidx += idxlist
241        #self.dump()
242   
243    def clearOverlays( self ):
244        self.overlaydict = {}
245        self.allidx = []
246        self.nextname = 'A'
247        self.initPlot(0,0)
248   
249    def initPlot( self, numtraces, traceorder ):
250        "Reset plot parameters before each redraw."
251        self.posdict = {}
252        self.posstore = []
253        self.posidx = -1
254        self.numtraces = numtraces
255        self.traceorder = traceorder
256   
257    def plotPos( self, trcidx, ypos ):
258        "Change plot positions of traces according to overlays."
259        if len(self.overlaydict.keys()) == 0:
260            return ypos
261        if self.traceorder == 0:
262            # traces always plotted top-down
263            trcidx = self.numtraces + 1 - trcidx
264        self.posstore.append( ypos )
265        for ovk in self.overlaydict.keys():
266            if trcidx in self.overlaydict[ovk]:
267                if ovk in self.posdict.keys():
268                    return self.posdict[ovk]
269                else:
270                    self.posidx += 1
271                    self.posdict[ovk] = self.posstore[self.posidx]
272                    return self.posstore[self.posidx]
273        self.posidx += 1
274        return self.posstore[self.posidx]
275   
276    def dump( self ):
277        print 'overlay dict', self.overlaydict
278        print 'overlay allidx', self.allidx
279        print 'overlay pos', self.posdict
280
281#@timeit
282def traces_from_list(selection):
283    """
284    Expand comma separated list of traces into Stream of internal traces.
285
286    Final Stream may contain visible and hidden traces!
287    """
288
289    if isinstance(selection, int):
290        selection = str(selection)
291
292    selection = selection.lower()
293    traces = Traces()
294
295    if selection.startswith("h:"):
296        parents = [traces.traces.hidden]
297        selection = selection[2:]
298    else:
299        parents = [traces.traces.visible]
300
301    selection = selection.split(",")
302
303    if "all" in selection:
304        return parents[0][:]
305
306    if selection[0].startswith("_"):
307        identifier = selection[0][1:].lower()
308
309        if identifier == "created":
310            return [t for t in traces.traces.created]
311
312        # pay attention to visible and invisible traces!
313        # this makes no sense in all cases but sometimes it does...
314        parents = [traces.traces.visible, traces.traces.hidden]
315
316        # get limitations
317        identifier = identifier.split("(")
318        if len(identifier) != 2:
319            return []
320
321        limits = identifier[1][:-1].split(":")
322        if len(limits) == 1:
323            limits.append(limits[0])
324
325        if identifier[0].endswith("~"):
326            ident = identifier[0][:-1]
327            negate = True
328        else:
329            ident = identifier[0]
330            negate = False
331
332        straces = []
333        for tl in parents:
334            for t in tl:
335                try:
336                    check = t.get_info(ident)
337                except:
338                    continue
339                if hasattr(check,'lower'):
340                    check = check.lower()
341               
342                found = False
343                try:
344                    if check >= type(check)(limits[0]) and \
345                                               check <= type(check)(limits[1]):
346                        found = True
347                except:
348                    continue
349
350                if found == negate:
351                    continue
352
353                straces.append(t)
354
355        return straces
356
357    numbers = expandTracelist(selection)
358
359    straces = []
360    for i in numbers:
361        try:
362            # get trace
363            straces.append(parents[0][i])
364        except IndexError:
365            log_message("warning", "trace number %s not found" % (i+1))
366
367    return Stream(straces)
368
369
370def number_of_visible_traces():
371    return len(Traces().traces.visible)
372
373
374#@timeit
375def expandTracelist(selection):
376    """
377    Returns number of trace selection, e.g. 1,3-5 results in [0,2,3,4]. The
378    special value "all" must be treated outside this method.
379
380    >>> expandTracelist("1,3")
381    [0, 2]
382    >>> expandTracelist("3-5,3")
383    [2, 3, 4]
384    >>> expandTracelist(["1", "8-11"])
385    [0, 7, 8, 9, 10]
386    """
387    if isinstance(selection, basestring):
388        selection = selection.split(",")
389
390    numbers = []
391    for s in selection:
392        # range
393        if "-" in s:
394            try:
395                start, stop = map(int, s.split("-"))
396            except:
397                print "  skipping '%s'" % s
398                continue
399
400            if stop < start:
401                start, stop = stop, start
402
403            numbers.extend(range(start - 1, stop))
404        else:
405            if s.isdigit():
406                numbers.append(int(s) - 1)
407            else:
408                print "  skipping '%s'" % s
409
410    # unique
411    numbers = list(set(numbers))
412    numbers.sort()
413
414    return numbers
415
416
417#@timeit
418def new_starttime(src, dst, offset):
419    """
420    Set new starttime with respect to offset, trace's T-ORIGIN and source
421    trace.
422    """
423    dst.stats.starttime = src.stats.starttime
424
425    t_origin = dst.get_info("t-origin")
426    corr = offset - t_origin
427    if corr > 0:
428        dst.stats.starttime += corr
429
430
431#@timeit
432def new_trace_data(length, zero=False):
433    """
434    Return numpy array of length "length". If "zero" is True, all values
435    will be set to 0.
436    """
437
438    if zero:
439        return np.zeros(length, dtype="float32")
440    else:
441        return np.empty(length, dtype="float32")
442
443
444def add_traces(stream, gain=False, uniquechannel=False):
445    """
446    Helper method for "Traces().addreplace(...)
447    """
448    Traces().addreplace(stream, applyGain=gain, uniquechannel=uniquechannel)
449
450def get_meta_status( trace, addgain=False ):
451    return Traces().getMetaStatus( trace, addgain=addgain )
452
453
454#@timeit
455def invalidate_all_traces(payload=None):
456    """
457    Voids all trace caches.
458    """
459    for t in Traces().visible:
460        t.invalidate_cache()
461subscribe_event(invalidate_all_traces, "recache")
462
463def addNetcode( station ):
464    "Adds netcode to station if not yet specified and found in cache."
465    traces = Traces()
466    upstat = station.upper()
467    if '.' in station:
468        return upstat
469    netcode = traces.getNetcode( upstat )
470    if netcode == None:
471        return upstat
472    return "%s.%s" % (netcode,upstat)
473
474def getDistance( trc, lat, lon ):
475    """
476    Compute distance in degrees from epicenter to station of trace.
477    'trc' is an obspy trace or a station name (at least NET.STATION).
478    """
479    if trc == None:
480        return None
481    if hasattr(trc,'lower'):
482        # trc is a string -> station name
483        sname = trc
484        # have to take time from first trace on display
485        try:
486            trc = Traces().traces.visible[0]
487        except:
488            return None
489    else:
490        # trc is a obspy trace, take station name from trace metadata
491        sname = "%s.%s" % (trc.stats.network,trc.stats.station)
492    metatime = trc.stats.starttime \
493        + (trc.stats.endtime-trc.stats.starttime)/2.
494    stations = Stations()
495    try:
496        r = stations[(sname,metatime)]
497        slat = float( r.latitude )
498        slon = float( r.longitude )
499    except:
500        return None
501    dist, az, baz = gps2dist_azimuth( lat, lon, slat, slon )
502    return (dist/1000./DEG_TO_KM)
503   
504
505if __name__ == "__main__":
506    import doctest
507    doctest.testmod(exclude_empty=True)
Note: See TracBrowser for help on using the repository browser.