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

Revision 394, 6.5 KB checked in by marcus, 13 years ago (diff)

more plotting tests

  • 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
31import wx
32from obspy.core import Stream, Trace
33#from obspy.core import UTCDateTime
34
35__all__ = [
36    "patched",
37]
38
39def _streamAutoMerge(self, *args, **kwargs):
40    """
41    Call obspy's Stream merge method to fill up gaps using the latest pre-gap
42    sample value. This destroys the spectra of the data but it's common.
43    """
44    if self.traces:
45        self.merge(method=1, fill_value='latest')
46
47def streamAutoMergeInit(self, *args, **kwargs):
48    self.__shx_init__(*args, **kwargs)
49    _streamAutoMerge(self)
50
51def streamAutoMergeAppend(self, *args, **kwargs):
52    self.__shx_append(*args, **kwargs)
53    _streamAutoMerge(self, *args, **kwargs)
54
55def streamAutoMergeExtend(self, *args, **kwargs):
56    self.__shx_extend(*args, **kwargs)
57    _streamAutoMerge(self, *args, **kwargs)
58
59def streamAutoMergeInsert(self, *args, **kwargs):
60    self.__shx_insert(*args, **kwargs)
61    _streamAutoMerge(self, *args, **kwargs)
62
63def streamAutoMergeAdd(self, *args, **kwargs):
64    self.__shx_add__(*args, **kwargs)
65    _streamAutoMerge(self, *args, **kwargs)
66
67def tracePrepareDataForImage(self, width, height=None, timescale=None, zoom=1, norm="window"):
68    """
69    Preparing trace data for fast plotting. Using numpy's minmax feature - this
70    will not just downsample the data!
71
72    width : total width of plotting window in pixels
73    height : hight of plotting window in pixels
74    zoom : global zoom factor (y extension)
75    timescale : extension of current timescale (tuple of min and max)
76    norm : norm method: "window", "total" or value (effects y extension)
77    """
78    if not timescale:
79        # use trace's total dimensions
80        windowstart, windowend = (self.stats.starttime, self.stats.endtime)
81    else:
82        windowstart, windowend = timescale
83
84    # XXX
85    if not height:
86        height = 1000
87
88    # get slice for time window
89    pt = self.slice(windowstart, windowend)
90
91    duration_total = windowend - windowstart
92    duration_trace = pt.stats.endtime - pt.stats.starttime
93
94    # width of the total screen that is covered by trace data (pixels)
95    pixel_width = duration_trace * width / duration_total
96    self.shxPixelWidth = pixel_width
97   
98    # data points per pixel
99    npts = pt.stats.npts
100    dpp = int(npts / pixel_width)
101
102    # get copy of data
103    self.shxPlotData = pt.data.copy()
104
105    # use minmax approach if more than 4 data points per pixel
106    # self.shx_plotdata will be changed
107    if dpp >= 4:
108        # do integer division since array shape values cannot be floats
109        dimx = npts // dpp
110        covered = dpp * dimx
111
112        # get data for reshape operation
113        data = self.shxPlotData[:covered]
114
115        # extra treatment for remaining values
116        remaining = self.shxPlotData[covered:]
117
118        # reshape data and get min/max values
119        data.shape = dimx, dpp
120        _min = data.min(axis=1)
121        _max = data.max(axis=1)
122        # combine data
123        joined = np.append([_min], [_max], axis=0).T.flatten()
124        # handle remaining
125        if len(remaining):
126            joined = np.append(joined, [remaining.min(), remaining.max()])
127
128        print "pixel_width", pixel_width, "minmax", "npts", npts, "dpp", dpp, "dimx", dimx, \
129              "covered", covered, "len_rest", len(remaining), "data_len", \
130              len(data), "joined", len(joined)
131
132        self.shxPlotData = joined
133
134    # At this stage the x-transformed data can be cached! If pixel_width doesn't
135    # change, the data can be reused to save time.
136
137    print width, height, norm, duration_trace, duration_total, pixel_width, dpp, type(self.shxPlotData)
138
139    # get basis for normation
140    if hasattr(norm, "lower"):
141        norm = norm.lower()
142        if norm == "window":
143            norm = abs(pt.max())
144        elif norm == "total":
145            norm = self.max()
146        else:
147            raise ValueError("Invalid input for normation!")
148
149    # Calculate y extension, leads to raw "image" data.
150    # Use a copy to save time if only y changes occur (e.g. zooming).
151    # Apply normation, calibration and total height factor
152    y = self.shxPlotData.copy()
153    y *= self.stats.calib / norm * height / 2 * zoom
154
155    print self.shxPlotData.min(), self.shxPlotData.max(), y.min(), y.max()
156
157    # factor for x stretching (trade-off between pixel_width and array shape)
158    stretch = pixel_width / len(y)
159    width = int(math.ceil(pixel_width))
160
161    # calculate line data
162    # do not rely on framework axis' management -> shift/flip y-axis
163    offset = height // 2
164    lines = []
165    oldx, oldy = 0, -y[0] + offset
166
167    for i in xrange(1, len(y)):
168        newx = i * stretch
169        newy = -y[i] + offset
170        lines.append([oldx, oldy, newx, newy])
171        oldx, oldy = newx, newy
172
173    self.shxImageData = lines
174
175def __traceX(parameters):
176    """
177    Only the x axis is adopted here. Zoom, norm etc. regarding the y extension
178    is done ...
179    """
180    pass
181
182# monkey patching obspy stream class: traces are merge automatically if stream
183# get's altered.
184Stream.__shx_init__, Stream.__init__ = Stream.__init__, streamAutoMergeInit
185Stream.__shx_append, Stream.append = Stream.append, streamAutoMergeAppend
186Stream.__shx_extend, Stream.extend = Stream.extend, streamAutoMergeExtend
187Stream.__shx_insert, Stream.insert = Stream.insert, streamAutoMergeInsert
188Stream.__shx_add__, Stream.__add__ = Stream.__add__, streamAutoMergeAdd
189
190# add method for faster plotting
191Trace.shxPrepareImageData = tracePrepareDataForImage
192
193# just an indicator
194patched = True
Note: See TracBrowser for help on using the repository browser.