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

Revision 1120, 6.2 KB checked in by klaus, 4 years ago (diff)

simulate menu entries

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