source: SHX/trunk/SeismicHandler/commands/simulate.py @ 931

Revision 931, 6.4 KB checked in by marcus, 6 years ago (diff)
  • 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 SeismicHandler.basics.command import BaseCommand
10from SeismicHandler.modules.traces import traces_from_list, add_traces
11from SeismicHandler.modules.stations import Stations
12from SeismicHandler.modules.filters import Instruments
13from SeismicHandler.basics.messages import log_message
14from SeismicHandler.basics.error import ShxError
15
16
17provides = {"simulate": "simulate"}
18class simulate(BaseCommand):
19    """
20    URI:http://www.seismic-handler.org/portal/wiki/ShSimulate
21    """
22    numberOfParameters = [0, 1, 2]
23    parameterQueries = [
24        {
25            "text": "trace list",
26            "type": "str",
27            "question": False,
28        },
29        {
30            "text": "simulation instrument",
31            "type": "str",
32            "question": False,
33        },
34    ]
35    known_qualifiers = [
36        "LIST",
37        "RESTORE",
38    ]
39
40    def __init__(self, *args, **kwargs):
41        # unroll args & kwargs
42        BaseCommand.__init__(self, *args, **kwargs)
43
44    def run(self):
45        if self.qualifiers["LIST"]:
46            # list all present transfer functions
47            path = os.path.join(os.path.split(__file__)[0], "..", "data")
48            for _, _, files in os.walk(path):
49                break
50
51            files.sort()
52            for _i in files:
53                print _i.lower().split('.')[0]
54
55            return
56
57        traces = traces_from_list(self.parameters[0])
58
59        if self.qualifiers["RESTORE"]:
60            for t in traces:
61                try:
62                    t.data = t.shx.original.data
63                    del t.shx.original
64                    t.invalidate_cache()
65                except:
66                    msg = "no unfiltered data for trace %u!"
67                    log_message("warning.simulate", msg % t.index(True))
68
69            return
70
71        instrument = self.parameters[1]
72
73        newtraces = []
74        msg = "applied filter to trace no %u"
75        for t in traces:
76            try:
77                newtraces.append(self.simulate(t, instrument))
78            except Exception:
79                continue
80            log_message("debug.filter", msg % t.index(True))
81        add_traces(newtraces)
82
83    @staticmethod
84    def simulate(trace, instrument):
85        """
86        Simulate historical instrument. Removes simultaneously instrument
87        response and apply response of simulated instrument.
88
89        The data is automatically pre-filtered to match the frequency range
90        of the simulation instrument:
91
92        Wood-Anderson: band-pass between 0.1 and 10 Hz
93        WWSSN-SP:      band-pass between 0.01 and 100 Hz
94        WWSSN-LP:      band-pass between 0.001 and 1 Hz
95        KIRNOS:        band-pass between 0.001 and 100 Hz
96        SRO-LP:        band-pass between 0.001 and 1 Hz
97        LRSM-SP:       band-pass between 0.01 and 1000 Hz
98        LRSM-LP:       band-pass between 0.0005 and 10 Hz
99
100        In any case the Nyquist freqency of the recording is used if it's
101        lower than the upper corner frequency.
102        """
103
104        instrument = instrument.lower()
105        tf = trace.copy()
106        tf.shx.original = trace
107
108        meta = Stations()
109        instruments = Instruments()
110
111        try:
112            remove = meta[(tf.id, tf.stats.starttime)]
113        except:
114            msg = "no meta data found for %s at start time %s"
115            raise ShxError(msg % (tf.id, tf.stats.starttime))
116
117        try:
118            simul = instruments[instrument]
119        except Exceptions as E:
120            raise ShxError(str(E))
121
122        # demean data
123        tf.detrend("demean")
124
125        # pre-filter data
126        nyquist = 1. / tf.stats.delta * .5
127        if instrument == "wood-anderson":
128            freqmax = nyquist > 10 and 10. or nyquist
129            tf.filter('bandpass', freqmin=0.1, freqmax=freqmax, corners=4,
130                                                                zerophase=True)
131        elif instrument == "wwssn-sp":
132            freqmax = nyquist > 100 and 100. or nyquist
133            tf.filter('bandpass', freqmin=0.01, freqmax=freqmax, corners=4,
134                                                                zerophase=True)
135        elif instrument == "wwssn-lp":
136            freqmax = nyquist > 1 and 1. or nyquist
137            tf.filter('bandpass', freqmin=0.001, freqmax=freqmax, corners=4,
138                                                                zerophase=True)
139        elif instrument == "kirnos":
140            freqmax = nyquist > 100 and 100. or nyquist
141            tf.filter('bandpass', freqmin=0.001, freqmax=freqmax, corners=4,
142                                                                zerophase=True)
143        elif instrument == "sro-lp":
144            freqmax = nyquist > 1 and 1. or nyquist
145            tf.filter('bandpass', freqmin=0.001, freqmax=freqmax, corners=4,
146                                                                zerophase=True)
147        elif instrument == "lrsm-lp":
148            freqmax = nyquist > 10 and 10. or nyquist
149            tf.filter('bandpass', freqmin=0.0005, freqmax=freqmax, corners=4,
150                                                                zerophase=True)
151        elif instrument == "lrsm-sp":
152            freqmax = nyquist > 1000 and 1000. or nyquist
153            tf.filter('bandpass', freqmin=0.01, freqmax=freqmax, corners=4,
154                                                                zerophase=True)
155        else:
156            msg = "Requested filter not found!"
157            raise ShxError(msg)
158
159        simul = simul.parse()[0]
160
161        simfilter = {
162            'poles': simul.output_coeff,
163            'zeros': simul.input_coeff,
164            'gain': 1,
165            'sensitivity': simul.normalization,
166        }
167
168        remfilter = {
169            'poles': remove._poles,
170            'zeros': remove._zeros,
171            'gain': 1,
172            'sensitivity': 1,
173        }
174
175        tf.simulate(
176            paz_remove=remfilter,
177            paz_simulate=simfilter,
178            taper=False,
179            shsim=True
180        )
181
182        msg = "%s filtered by %s" % (tf.id, instrument)
183        log_message("debug.filter", msg)
184
185        # retain SH data standard
186        tf.data = np.require(tf.data, "float32")
187        tf.invalidate_cache()
188
189        return tf
190
191
192if __name__ == "__main__":
193    import doctest
194    doctest.testmod(exclude_empty=True)
Note: See TracBrowser for help on using the repository browser.