source: SHX/trunk/SeismicHandler/patches/ObsPy.py @ 529

Revision 529, 8.6 KB checked in by marcus, 11 years ago (diff)
  • accept tailing comment within command line
  • access to trace info entries
  • Property svn:eol-style set to native
Line 
1# -*- coding: utf-8 -*-
2#
3# Copyright (C) 2008-2011 Marcus Walther (walther@szgrf.bgr.de)
4#
5# This file is part of Seismic Handler eXtended (SHX)
6# Full details can be found at project website http://www.seismic-handler.org/
7#
8# SHX is free software; you can redistribute it and/or modify
9# it under the terms of the GNU LESSER GENERAL PUBLIC LICENSE as published by
10# the Free Software Foundation; either version 3 of the License, or
11# (at your option) any later version.
12#
13# SHX is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16# GNU Lesser General Public License for more details.
17#
18# You should have received a copy of the GNU Lesser General Public License
19# along with SHX (see license.txt); if not, write to the Free Software
20# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
21
22"""
23This modules defines extensions to obspy. They are patched into the
24original obspy modules at a very basic level in SeismicHandler/__init__.py
25
26This module is named ObsPy.py to prevent import loops.
27"""
28
29import numpy as np
30import math
31from obspy.core import Stream, Trace
32from obspy.core.util import AttribDict
33from obspy.sh.core import fromUTCDateTime
34from SeismicHandler.basics.messages import logMessage, setStatus
35
36__all__ = []
37
38def _streamAutoMerge(self, *args, **kwargs):
39    """
40    Call obspy's Stream merge method to fill up gaps using the latest pre-gap
41    sample value. This destroys the spectra of the data but it's common.
42    """
43
44    # Change value data type to float32 as in former versions of SH.
45    # This will be changed in later versions.
46    for t in self.traces:
47        if t.data.dtype != np.dtype('float32'):
48            t.data = t.data.astype('float32')
49
50    if self.traces:
51        self.merge(method=1, fill_value='latest')
52
53def streamAutoMergeInit(self, *args, **kwargs):
54    self.__shx_init__(*args, **kwargs)
55    _streamAutoMerge(self)
56
57def streamAutoMergeAppend(self, *args, **kwargs):
58    self.__shx_append(*args, **kwargs)
59    _streamAutoMerge(self, *args, **kwargs)
60
61def streamAutoMergeExtend(self, *args, **kwargs):
62    self.__shx_extend(*args, **kwargs)
63    _streamAutoMerge(self, *args, **kwargs)
64
65def streamAutoMergeInsert(self, *args, **kwargs):
66    self.__shx_insert(*args, **kwargs)
67    _streamAutoMerge(self, *args, **kwargs)
68
69def streamAutoMergeAdd(self, *args, **kwargs):
70    self.__shx_add__(*args, **kwargs)
71    _streamAutoMerge(self, *args, **kwargs)
72
73def tracePrepareDataForImage(self, width, height, timescale=None, zoom=1, norm="window"):
74    """
75    Preparing trace data for fast plotting. Using numpy's minmax feature - this
76    will not just downsample the data!
77
78    width : total width of plotting window in pixels
79    height : hight of plotting window in pixels
80    zoom : global zoom factor (y extension)
81    timescale : extension of current timescale (tuple of min and max)
82    norm : norm method: "window", "total" or value (effects y extension)
83    """
84    logMessage('debug.patches', "Preparing plotting data...", self.id)
85
86    if not timescale:
87        # use trace's total dimensions
88        windowstart, windowend = (self.stats.starttime, self.stats.endtime)
89    else:
90        windowstart, windowend = timescale
91
92    # get slice for time window
93    pt = self.slice(windowstart, windowend)
94
95    duration_total = windowend - windowstart
96    duration_trace = pt.stats.endtime - pt.stats.starttime
97
98    # width of the total screen that is covered by trace data (pixels)
99    pixel_width = duration_trace * width / duration_total
100
101    # remember values for caching
102    CacheTraceID = \
103              (windowstart, windowend, "%.5f" % pixel_width, height, zoom, norm)
104    try:
105        cachedX = CacheTraceID[:3] == self.shx.CacheTraceID[:3]
106        cachedY = cachedX and CacheTraceID[3:] == self.shx.CacheTraceID[3:]
107    except:
108        cachedX = False
109        cachedY = False
110
111    if not cachedX:
112        # save new ID
113        self.shx.CacheTraceID = CacheTraceID
114
115        # data points per pixel
116        npts = pt.stats.npts
117        dpp = int(npts / pixel_width)
118
119        # get copy of data
120        self.shx.PlotData = pt.data.copy()
121    else:
122        logMessage('debug.patches.devel',
123                                    "using cachedX data for", CacheTraceID[:3])
124
125    # use minmax approach if more than 4 data points per pixel
126    # self.shx_plotdata will be changed
127    if not cachedX and dpp >= 4:
128        # do integer division since array shape values cannot be floats
129        dimx = npts // dpp
130        covered = dpp * dimx
131
132        # get data for reshape operation
133        data = self.shx.PlotData[:covered]
134
135        # extra treatment for remaining values
136        remaining = self.shx.PlotData[covered:]
137
138        # reshape data and get min/max values
139        data.shape = dimx, dpp
140        _min = data.min(axis=1)
141        _max = data.max(axis=1)
142        # combine data
143        joined = np.append([_min], [_max], axis=0).T.flatten()
144        # handle remaining
145        if len(remaining):
146            joined = np.append(joined, [remaining.min(), remaining.max()])
147
148        logMessage(
149            'debug.patches.devel',
150             "pixel_width", pixel_width, "minmax", "npts", npts, "dpp",
151             dpp, "dimx", dimx, "covered", covered, "len_rest",
152             len(remaining), "data_len", len(data), "joined", len(joined),
153             "width", width
154        )
155
156        # At this stage the x-transformed data can be cached! If pixel_width
157        # doesn't change, the data can be reused to save time.
158        self.shx.PlotData = joined
159
160    logMessage('debug.patches.devel',
161                                  width, height, norm, type(self.shx.PlotData))
162
163    if cachedY:
164        logMessage('debug.patches.devel',
165                                    "using cachedY data for", CacheTraceID[3:])
166        return
167   
168    # renew cache id
169    self.shx.CacheTraceID = CacheTraceID
170
171    # get basis for normation
172    if hasattr(norm, "lower"):
173        norm = norm.lower()
174        if norm == "window":
175            norm = abs(pt.max())
176        elif norm == "total":
177            norm = self.max()
178        else:
179            raise ValueError("Invalid input for normation!")
180
181    # Calculate y extension, leads to raw "image" data.
182    # Use a copy to save time if only y changes occur (e.g. zooming).
183    # Apply normation, calibration and total height factor
184    y = self.shx.PlotData.copy()
185    y *= self.stats.calib / norm * height / 2 * zoom
186
187#    print self.shxPlotData.min(), self.shxPlotData.max(), y.min(), y.max()
188
189    # factor for x stretching (trade-off between pixel_width and array shape)
190    stretch = pixel_width / len(y)
191    width = int(math.ceil(pixel_width))
192
193    # calculate line data
194    # do not rely on framework axis' management -> shift/flip y-axis
195    offset = height // 2
196    lines = []
197    oldx, oldy = 0, -y[0] + offset
198
199    for i in xrange(1, len(y)):
200        newx = i * stretch
201        newy = -y[i] + offset
202        lines.append([oldx, oldy, newx, newy])
203        oldx, oldy = newx, newy
204
205    self.shx.ImageData = lines
206    self.shx.PlotPixels = pixel_width
207
208def traceAddShxInit(self, *args, **kwargs):
209    self.__shx_init__(*args, **kwargs)
210
211    # container for SHX releated data
212    self.shx = AttribDict()
213
214def traceGetStatusInfo(self, name):
215    """
216    Mapping of Seismic Handler trace information entries to ObsPy status info.
217    """
218    name = name.lower()
219
220    # identical name
221    if name in self.stats:
222        return self.stats[name]
223
224    # mapping
225    if name == "length":
226        return self.stats.npts
227
228    if name == "chan1":
229        return self.stats.channel[0]
230
231    if name == "chan2":
232        return self.stats.channel[1]
233
234    if name == "comp":
235        return self.stats.channel[2]
236
237    if name == "start":
238        return fromUTCDateTime(self.stats.starttime)
239
240    # finally look into self defined entries
241    if name in self.shx:
242        return self.shx[name]
243
244    raise NameError
245
246# monkey patching obspy stream class: traces are merge automatically if stream
247# get's altered.
248Stream.__shx_init__, Stream.__init__ = Stream.__init__, streamAutoMergeInit
249#Stream.__shx_append, Stream.append = Stream.append, streamAutoMergeAppend
250#Stream.__shx_extend, Stream.extend = Stream.extend, streamAutoMergeExtend
251#Stream.__shx_insert, Stream.insert = Stream.insert, streamAutoMergeInsert
252#Stream.__shx_add__, Stream.__add__ = Stream.__add__, streamAutoMergeAdd
253
254# monkey patching obspy trace class: add shx container and method for fast
255# plotting
256Trace.__shx_init__, Trace.__init__ = Trace.__init__, traceAddShxInit
257Trace._shxPrepareImageData = tracePrepareDataForImage
258Trace._shxInfo = traceGetStatusInfo
Note: See TracBrowser for help on using the repository browser.