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

Revision 911, 5.7 KB checked in by marcus, 6 years ago (diff)
  • minor fixes
  • 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            newtraces.append(self.simulate(t, instrument))
77            log_message("debug.filter", msg % t.index(True))
78        add_traces(newtraces)
79
80    @staticmethod
81    def simulate(trace, instrument):
82        """
83        Simulate historical instrument. Removes simultaneously instrument
84        response and apply response of simulated instrument.
85
86        The data is automatically pre-filtered to match the frequency range
87        of the simulation instrument:
88
89        Wood-Anderson: band-pass between 0.1 and 10 Hz
90        WWSSN-SP:      band-pass between 0.01 and 100 Hz
91        WWSSN-LP:      band-pass between 0.001 and 1 Hz
92        KIRNOS:        band-pass between 0.001 and 100 Hz
93        SRO-LP:        band-pass between 0.001 and 1 Hz
94        LRSM-SP:       band-pass between 0.01 and 1000 Hz
95        LRSM-LP:       band-pass between 0.0005 and 10 Hz
96        """
97
98        instrument = instrument.lower()
99        tf = trace.copy()
100        tf.shx.original = trace
101
102        meta = Stations()
103        instruments = Instruments()
104
105        try:
106            remove = meta[(tf.id, tf.stats.starttime)]
107        except:
108            msg = "no meta data found for time starting at %s"
109            raise ShxError(msg % tf.stats.starttime)
110
111        try:
112            simul = instruments[instrument]
113        except Exceptions as E:
114            raise ShxError(str(E))
115
116        # demean data
117        tf.detrend("demean")
118
119        # pre-filter data
120        if instrument == "wood-anderson":
121            tf.filter('bandpass', freqmin=0.1, freqmax=10, corners=4,
122                                                                zerophase=True)
123        elif instrument == "wwssn-sp":
124            tf.filter('bandpass', freqmin=0.01, freqmax=100, corners=4,
125                                                                zerophase=True)
126        elif instrument == "wwssn-lp":
127            tf.filter('bandpass', freqmin=0.001, freqmax=1, corners=4,
128                                                                zerophase=True)
129        elif instrument == "kirnos":
130            tf.filter('bandpass', freqmin=0.001, freqmax=100, corners=4,
131                                                                zerophase=True)
132        elif instrument == "sro-lp":
133            tf.filter('bandpass', freqmin=0.001, freqmax=1, corners=4,
134                                                                zerophase=True)
135        elif instrument == "lrsm-lp":
136            tf.filter('bandpass', freqmin=0.0005, freqmax=10, corners=4,
137                                                                zerophase=True)
138        elif instrument == "lrsm-sp":
139            tf.filter('bandpass', freqmin=0.01, freqmax=1000, corners=4,
140                                                                zerophase=True)
141
142        simul = simul.parse()[0]
143
144        simfilter = {
145            'poles': simul.output_coeff,
146            'zeros': simul.input_coeff,
147            'gain': 1,
148            'sensitivity': simul.normalization,
149        }
150
151        remfilter = {
152            'poles': remove._poles,
153            'zeros': remove._zeros,
154            'gain': 1,
155            'sensitivity': 1,
156        }
157
158        tf.simulate(
159            paz_remove=remfilter,
160            paz_simulate=simfilter,
161            taper=False,
162            shsim=True
163        )
164
165        msg = "%s filtered by %s" % (tf.id, instrument)
166        log_message("debug.filter", msg)
167
168        # retain SH data standard
169        tf.data = np.require(tf.data, "float32")
170
171        return tf
172
173
174if __name__ == "__main__":
175    import doctest
176    doctest.testmod(exclude_empty=True)
Note: See TracBrowser for help on using the repository browser.