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 | """ |
---|
23 | This modules defines extensions to obspy. They are patched into the |
---|
24 | original obspy modules at a very basic level in SeismicHandler/__init__.py |
---|
25 | |
---|
26 | This module is named ObsPy.py to prevent import loops. |
---|
27 | """ |
---|
28 | |
---|
29 | import numpy as np |
---|
30 | import math |
---|
31 | from obspy.core import Stream, Trace |
---|
32 | from obspy.core.util import AttribDict |
---|
33 | from SeismicHandler.utils.pubsub import pub as msgs |
---|
34 | #from obspy.core import UTCDateTime |
---|
35 | |
---|
36 | __all__ = [ |
---|
37 | "patched", |
---|
38 | ] |
---|
39 | |
---|
40 | def _streamAutoMerge(self, *args, **kwargs): |
---|
41 | """ |
---|
42 | Call obspy's Stream merge method to fill up gaps using the latest pre-gap |
---|
43 | sample value. This destroys the spectra of the data but it's common. |
---|
44 | """ |
---|
45 | if self.traces: |
---|
46 | self.merge(method=1, fill_value='latest') |
---|
47 | |
---|
48 | def streamAutoMergeInit(self, *args, **kwargs): |
---|
49 | self.__shx_init__(*args, **kwargs) |
---|
50 | _streamAutoMerge(self) |
---|
51 | |
---|
52 | def streamAutoMergeAppend(self, *args, **kwargs): |
---|
53 | self.__shx_append(*args, **kwargs) |
---|
54 | _streamAutoMerge(self, *args, **kwargs) |
---|
55 | |
---|
56 | def streamAutoMergeExtend(self, *args, **kwargs): |
---|
57 | self.__shx_extend(*args, **kwargs) |
---|
58 | _streamAutoMerge(self, *args, **kwargs) |
---|
59 | |
---|
60 | def streamAutoMergeInsert(self, *args, **kwargs): |
---|
61 | self.__shx_insert(*args, **kwargs) |
---|
62 | _streamAutoMerge(self, *args, **kwargs) |
---|
63 | |
---|
64 | def streamAutoMergeAdd(self, *args, **kwargs): |
---|
65 | self.__shx_add__(*args, **kwargs) |
---|
66 | _streamAutoMerge(self, *args, **kwargs) |
---|
67 | |
---|
68 | def tracePrepareDataForImage(self, width, height, timescale=None, zoom=1, norm="window"): |
---|
69 | """ |
---|
70 | Preparing trace data for fast plotting. Using numpy's minmax feature - this |
---|
71 | will not just downsample the data! |
---|
72 | |
---|
73 | width : total width of plotting window in pixels |
---|
74 | height : hight of plotting window in pixels |
---|
75 | zoom : global zoom factor (y extension) |
---|
76 | timescale : extension of current timescale (tuple of min and max) |
---|
77 | norm : norm method: "window", "total" or value (effects y extension) |
---|
78 | """ |
---|
79 | msgs.sendMessage('log.debug.patches', |
---|
80 | message=["Preparing plotting data...", self.id]) |
---|
81 | |
---|
82 | if not timescale: |
---|
83 | # use trace's total dimensions |
---|
84 | windowstart, windowend = (self.stats.starttime, self.stats.endtime) |
---|
85 | else: |
---|
86 | windowstart, windowend = timescale |
---|
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 | |
---|
97 | # remember values for caching |
---|
98 | CacheTraceID = \ |
---|
99 | (windowstart, windowend, "%.5f" % pixel_width, height, zoom, norm) |
---|
100 | try: |
---|
101 | cachedX = CacheTraceID[:3] == self.shx.CacheTraceID[:3] |
---|
102 | cachedY = cachedX and CacheTraceID[3:] == self.shx.CacheTraceID[3:] |
---|
103 | except: |
---|
104 | cachedX = False |
---|
105 | cachedY = False |
---|
106 | |
---|
107 | if not cachedX: |
---|
108 | # save new ID |
---|
109 | self.shx.CacheTraceID = CacheTraceID |
---|
110 | |
---|
111 | # data points per pixel |
---|
112 | npts = pt.stats.npts |
---|
113 | dpp = int(npts / pixel_width) |
---|
114 | |
---|
115 | # get copy of data |
---|
116 | self.shx.PlotData = pt.data.copy() |
---|
117 | else: |
---|
118 | msgs.sendMessage('log.devel.patches', message=["using cachedX data for", |
---|
119 | CacheTraceID[:3]]) |
---|
120 | |
---|
121 | # use minmax approach if more than 4 data points per pixel |
---|
122 | # self.shx_plotdata will be changed |
---|
123 | if not cachedX and dpp >= 4: |
---|
124 | # do integer division since array shape values cannot be floats |
---|
125 | dimx = npts // dpp |
---|
126 | covered = dpp * dimx |
---|
127 | |
---|
128 | # get data for reshape operation |
---|
129 | data = self.shx.PlotData[:covered] |
---|
130 | |
---|
131 | # extra treatment for remaining values |
---|
132 | remaining = self.shx.PlotData[covered:] |
---|
133 | |
---|
134 | # reshape data and get min/max values |
---|
135 | data.shape = dimx, dpp |
---|
136 | _min = data.min(axis=1) |
---|
137 | _max = data.max(axis=1) |
---|
138 | # combine data |
---|
139 | joined = np.append([_min], [_max], axis=0).T.flatten() |
---|
140 | # handle remaining |
---|
141 | if len(remaining): |
---|
142 | joined = np.append(joined, [remaining.min(), remaining.max()]) |
---|
143 | |
---|
144 | msgs.sendMessage( |
---|
145 | 'log.devel.patches', |
---|
146 | message=[ |
---|
147 | "pixel_width", pixel_width, "minmax", "npts", npts, "dpp", |
---|
148 | dpp, "dimx", dimx, "covered", covered, "len_rest", |
---|
149 | len(remaining), "data_len", len(data), "joined", len(joined), |
---|
150 | "width", width |
---|
151 | ] |
---|
152 | ) |
---|
153 | |
---|
154 | # At this stage the x-transformed data can be cached! If pixel_width |
---|
155 | # doesn't change, the data can be reused to save time. |
---|
156 | self.shx.PlotData = joined |
---|
157 | |
---|
158 | msgs.sendMessage('log.devel.patches', |
---|
159 | message=[width, height, norm, type(self.shx.PlotData)]) |
---|
160 | |
---|
161 | if cachedY: |
---|
162 | msgs.sendMessage('log.devel.patches', |
---|
163 | message=["using cachedY data for", CacheTraceID[3:]]) |
---|
164 | return |
---|
165 | |
---|
166 | # renew cache id |
---|
167 | self.shx.CacheTraceID = CacheTraceID |
---|
168 | |
---|
169 | # get basis for normation |
---|
170 | if hasattr(norm, "lower"): |
---|
171 | norm = norm.lower() |
---|
172 | if norm == "window": |
---|
173 | norm = abs(pt.max()) |
---|
174 | elif norm == "total": |
---|
175 | norm = self.max() |
---|
176 | else: |
---|
177 | raise ValueError("Invalid input for normation!") |
---|
178 | |
---|
179 | # Calculate y extension, leads to raw "image" data. |
---|
180 | # Use a copy to save time if only y changes occur (e.g. zooming). |
---|
181 | # Apply normation, calibration and total height factor |
---|
182 | y = self.shx.PlotData.copy() |
---|
183 | y *= self.stats.calib / norm * height / 2 * zoom |
---|
184 | |
---|
185 | # print self.shxPlotData.min(), self.shxPlotData.max(), y.min(), y.max() |
---|
186 | |
---|
187 | # factor for x stretching (trade-off between pixel_width and array shape) |
---|
188 | stretch = pixel_width / len(y) |
---|
189 | width = int(math.ceil(pixel_width)) |
---|
190 | |
---|
191 | # calculate line data |
---|
192 | # do not rely on framework axis' management -> shift/flip y-axis |
---|
193 | offset = height // 2 |
---|
194 | lines = [] |
---|
195 | oldx, oldy = 0, -y[0] + offset |
---|
196 | |
---|
197 | for i in xrange(1, len(y)): |
---|
198 | newx = i * stretch |
---|
199 | newy = -y[i] + offset |
---|
200 | lines.append([oldx, oldy, newx, newy]) |
---|
201 | oldx, oldy = newx, newy |
---|
202 | |
---|
203 | self.shx.ImageData = lines |
---|
204 | self.shx.PlotPixels = pixel_width |
---|
205 | |
---|
206 | def traceAddShxInit(self, *args, **kwargs): |
---|
207 | self.__shx_init__(*args, **kwargs) |
---|
208 | |
---|
209 | # container for SHX releated data |
---|
210 | self.shx = AttribDict() |
---|
211 | |
---|
212 | # monkey patching obspy stream class: traces are merge automatically if stream |
---|
213 | # get's altered. |
---|
214 | Stream.__shx_init__, Stream.__init__ = Stream.__init__, streamAutoMergeInit |
---|
215 | Stream.__shx_append, Stream.append = Stream.append, streamAutoMergeAppend |
---|
216 | Stream.__shx_extend, Stream.extend = Stream.extend, streamAutoMergeExtend |
---|
217 | Stream.__shx_insert, Stream.insert = Stream.insert, streamAutoMergeInsert |
---|
218 | Stream.__shx_add__, Stream.__add__ = Stream.__add__, streamAutoMergeAdd |
---|
219 | |
---|
220 | # monkey patching obspy trace class: add shx container and method for fast |
---|
221 | # plotting |
---|
222 | Trace.__shx_init__, Trace.__init__ = Trace.__init__, traceAddShxInit |
---|
223 | Trace._shxPrepareImageData = tracePrepareDataForImage |
---|
224 | |
---|
225 | # just an indicator |
---|
226 | patched = True |
---|