Coverage for /usr/local/lib/python3.11/dist-packages/grond/analysers/noise_analyser/analyser.py: 19%
147 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-11-27 15:15 +0000
1# https://pyrocko.org/grond - GPLv3
2#
3# The Grond Developers, 21st Century
5from pyrocko.client import catalog
7import logging
8import numpy as num
9from pyrocko.guts import Int, Bool, Float, String, StringChoice
10from pyrocko.gf.meta import OutOfBounds
11from ..base import Analyser, AnalyserConfig, AnalyserResult
12from grond.dataset import NotFound
13from grond.meta import store_t
15logger = logging.getLogger('grond.analysers.NoiseAnalyser')
18guts_prefix = 'grond'
21def get_phase_arrival_time(engine, source, target, wavename):
22 """
23 Get arrival time from Green's Function store for respective
24 :class:`pyrocko.gf.seismosizer.Target`,
25 :class:`pyrocko.gf.meta.Location` pair.
27 Parameters
28 ----------
29 engine : :class:`pyrocko.gf.seismosizer.LocalEngine`
30 source : :class:`pyrocko.gf.meta.Location`
31 can be therefore :class:`pyrocko.gf.seismosizer.Source` or
32 :class:`pyrocko.model.Event`
33 target : :class:`pyrocko.gf.seismosizer.Target`
34 wavename : string
35 of the tabulated phase_def that determines the phase arrival
37 Returns
38 -------
39 scalar, float of the arrival time of the wave
40 """
41 store = engine.get_store(target.store_id)
42 return store_t(store, wavename, source, target) + source.time
45def seismic_noise_variance(traces, engine, source, targets,
46 nwindows, pre_event_noise_duration,
47 check_events, phase_def):
48 """
49 Calculate variance of noise in a given time before P-Phase onset.
51 Optionally check the gCMT earthquake catalogue for M>5 events interfering.
53 Parameters
54 ----------
55 data_traces : list
56 of :class:`pyrocko.trace.Trace` containing observed data
57 engine : :class:`pyrocko.gf.seismosizer.LocalEngine`
58 processing object for synthetics calculation
59 source : :class:`pyrocko.gf.Source`
60 reference source
61 targets : list
62 of :class:`pyrocko.gf.seismosizer.Targets`
63 nwindows : integer
64 number of windows in which the noise trace is split. If not 1, the
65 variance is calculated for each window separately and a mean
66 variance is returned. Else, the variance is calculated on the
67 entire pre-event noise window.
68 pre_event_noise_duration : Time before the first arrival to include in the
69 noise analysis
70 phase_def : :class:'pyrocko.gf.Timing'
71 arrivals : list
72 of :class'pyrocko.gf.Timing' arrivals of waveforms
73 at station
75 Returns
76 -------
77 :class:`numpy.ndarray`
78 """
80 var_ds = []
81 global_cmt_catalog = catalog.GlobalCMT()
82 var_ds = []
83 ev_ws = []
84 for tr, target in zip(traces, targets):
85 stat_w = 1.
87 if tr is None:
88 var_ds.append(num.nan)
89 ev_ws.append(num.nan)
90 else:
92 arrival_time = get_phase_arrival_time(
93 engine=engine, source=source,
94 target=target, wavename=phase_def)
95 if check_events:
96 events = global_cmt_catalog.get_events(
97 time_range=(
98 arrival_time-pre_event_noise_duration-50.*60.,
99 arrival_time),
100 magmin=5.,)
101 for ev in events:
102 try:
103 arrival_time_pre = get_phase_arrival_time(
104 engine=engine,
105 source=ev,
106 target=target,
107 wavename=phase_def)
109 if arrival_time_pre > arrival_time \
110 - pre_event_noise_duration \
111 and arrival_time_pre < arrival_time:
113 stat_w = 0.
114 logger.info(
115 'Noise analyser found event "%s" phase onset '
116 'of "%s" for target "%s".' % (
117 ev.name, phase_def, target.name))
119 if arrival_time_pre > arrival_time-30.*60.\
120 and arrival_time_pre < arrival_time - \
121 pre_event_noise_duration:
122 stat_w *= 1.
123 logger.info(
124 'Noise analyser found event "%s" possibly '
125 'contaminating the noise.' % ev.name)
127 # this should be magnitude dependent
128 except Exception:
129 pass
130 ev_ws.append(stat_w)
132 if nwindows == 1:
133 vtrace_var = num.nanvar(tr.ydata)
134 var_ds.append(vtrace_var)
135 else:
136 win = arrival_time - (arrival_time -
137 pre_event_noise_duration)
138 win_len = win/nwindows
139 v_traces_w = []
140 for i in range(0, nwindows):
141 vtrace_w = tr.chop(
142 tmin=win+win_len*i,
143 tmax=arrival_time+win_len*i+1,
144 inplace=False)
145 v_traces_w.append(vtrace_w.ydata)
146 v_traces_w = num.nanmean(v_traces_w)
147 var_ds.append(v_traces_w)
148 var_ds = num.array(var_ds, dtype=float)
149 ev_ws = num.array(ev_ws, dtype=float)
150 return var_ds, ev_ws
153class NoiseAnalyser(Analyser):
154 '''
155 From the pre-event station noise variance-based trace weights are formed.
157 By default, the trace weights are the inverse of the noise variance. The
158 correlation of the noise is neglected. Optionally, using a the gCMT global
159 earthquake catalogue, the station data are checked for theoretical phase
160 arrivals of M>5 earthquakes. In case of a very probable contamination the
161 trace weights are set to zero. In case global earthquake phase arrivals are
162 within a 30 min time window before the start of the set pre-event noise
163 window, only a warning is thrown.
165 It is further possible to disregard data with a noise level exceeding the
166 median by a given ``cutoff`` factor. These weights are set to 0. This can
167 be done exclusively (``mode='weeding'``) such that noise weights are either
168 1 or 0, or in combination with weighting below the median-times-cutoff
169 noise level (``mode='weighting'``).
170 '''
172 def __init__(self, nwindows, pre_event_noise_duration,
173 check_events, phase_def, statistic, mode, cutoff,
174 cutoff_exception_on_high_snr):
176 Analyser.__init__(self)
177 self.nwindows = nwindows
178 self.pre_event_noise_duration = pre_event_noise_duration
179 self.check_events = check_events
180 self.phase_def = phase_def
181 self.statistic = statistic
182 self.mode = mode
183 self.cutoff = cutoff
184 self.cutoff_exception_on_high_snr = cutoff_exception_on_high_snr
186 def analyse(self, problem, ds):
188 tdur = self.pre_event_noise_duration
190 if tdur == 0:
191 return
193 if not problem.has_waveforms:
194 return
196 engine = problem.get_engine()
197 source = problem.base_source
199 paths = sorted(set(t.path for t in problem.waveform_targets))
201 for path in paths:
202 targets = [t for t in problem.waveform_targets if t.path == path]
204 deltats = set()
205 for target in targets: # deltat diff check?
206 store = engine.get_store(target.store_id)
207 deltats.add(float(store.config.deltat))
209 if len(deltats) > 1:
210 logger.warning(
211 'Differing sampling rates in stores used. Using highest.')
213 deltat = min(deltats)
215 data = []
216 for target in targets:
217 try:
218 freqlimits = list(target.get_freqlimits())
219 freqlimits = tuple(freqlimits)
221 source = problem.base_source
223 arrival_time = get_phase_arrival_time(
224 engine=engine,
225 source=source,
226 target=target,
227 wavename=self.phase_def)
229 tmin_fit, tmax_fit, tfade, tfade_taper = \
230 target.get_taper_params(engine, source)
232 data.append([
233 tmin_fit,
234 tmax_fit,
235 tfade_taper,
236 ds.get_waveform(
237 target.codes,
238 tmin=arrival_time-tdur-tfade,
239 tmax=tmax_fit+tfade,
240 tfade=tfade,
241 freqlimits=freqlimits,
242 deltat=deltat,
243 backazimuth=target.get_backazimuth_for_waveform(),
244 tinc_cache=1./freqlimits[0],
245 debug=False)])
247 except (NotFound, OutOfBounds) as e:
248 logger.debug(str(e))
249 data.append([None, None, None, None])
251 traces_noise = []
252 traces_signal = []
253 for tmin_fit, tmax_fit, tfade_taper, tr in data:
254 if tr:
255 traces_noise.append(
256 tr.chop(tr.tmin, tr.tmin + tdur, inplace=False))
257 traces_signal.append(
258 tr.chop(
259 tmin_fit-tfade_taper,
260 tmax_fit+tfade_taper,
261 inplace=False))
262 else:
263 traces_noise.append(None)
264 traces_signal.append(None)
266 var_ds, ev_ws = seismic_noise_variance(
267 traces_noise, engine, source, targets,
268 self.nwindows, tdur,
269 self.check_events, self.phase_def)
271 amp_maxs = num.array([
272 (tr.absmax()[1] if tr else num.nan) for tr in traces_signal])
274 if self.statistic == 'var':
275 noise = var_ds
276 elif self.statistic == 'std':
277 noise = num.sqrt(var_ds)
278 else:
279 assert False, 'invalid statistic argument'
281 ok = num.isfinite(noise)
283 if num.sum(ok) == 0:
284 norm_noise = 0.0
285 else:
286 norm_noise = num.median(noise[ok])
288 if norm_noise == 0.0:
289 logger.info(
290 'Noise Analyser returned a weight of 0 for all stations.')
292 assert num.all(noise[ok] >= 0.0)
294 ce_factor = self.cutoff_exception_on_high_snr
295 high_snr = num.zeros(ok.size, dtype=bool)
296 if ce_factor is not None:
297 high_snr[ok] = amp_maxs[ok] > ce_factor * num.sqrt(var_ds)[ok]
299 weights = num.zeros(noise.size)
300 if self.mode == 'weighting':
301 weights[ok] = norm_noise / noise[ok]
302 elif self.mode == 'weeding':
303 weights[ok] = 1.0
304 else:
305 assert False, 'invalid mode argument'
307 if self.cutoff is not None:
308 weights[ok] = num.where(
309 num.logical_or(
310 noise[ok] <= norm_noise * self.cutoff,
311 high_snr[ok]),
312 weights[ok], 0.0)
314 if self.check_events:
315 weights = weights*ev_ws
317 for itarget, target in enumerate(targets):
318 logger.info((
319 'Noise analysis for target "%s":\n'
320 ' var: %g\n'
321 ' std: %g\n'
322 ' max/std: %g\n'
323 ' %s/median(%s): %g\n'
324 ' contamination_weight: %g\n'
325 ' weight: %g') % (
326 target.string_id(),
327 var_ds[itarget],
328 num.sqrt(var_ds[itarget]),
329 amp_maxs[itarget] / num.sqrt(var_ds[itarget]),
330 self.statistic, self.statistic,
331 noise[itarget] / norm_noise,
332 ev_ws[itarget],
333 weights[itarget]))
335 for weight, target in zip(weights, targets):
336 target.analyser_results['noise'] = \
337 NoiseAnalyserResult(weight=float(weight))
340class NoiseAnalyserResult(AnalyserResult):
341 weight = Float.T(
342 help='The inverse of the pre-event data variance or standard '
343 'deviation. If traces were checked for other event phase '
344 'arrivals, the weight can be zero for contaminated traces.')
347class NoiseAnalyserConfig(AnalyserConfig):
348 """Configuration parameters for the pre-event noise analysis."""
350 nwindows = Int.T(
351 default=1,
352 help='number of windows for trace splitting')
354 pre_event_noise_duration = Float.T(
355 default=0.,
356 help='Total length of noise trace in the analysis')
358 phase_def = String.T(
359 default='P',
360 help='Onset of phase_def used for upper limit of window')
362 check_events = Bool.T(
363 default=False,
364 help='check the GlobalCMT for M>5 earthquakes'
365 ' that produce phase arrivals'
366 ' contaminating and affecting the noise analysis')
368 statistic = StringChoice.T(
369 choices=('var', 'std'),
370 default='var',
371 help='Set weight to inverse of noise variance (var) or standard '
372 'deviation (std).')
374 mode = StringChoice.T(
375 choices=('weighting', 'weeding'),
376 default='weighting',
377 help='Generate weights based on inverse of noise measure (weighting), '
378 'or discrete on/off style in combination with cutoff value '
379 '(weeding).')
381 cutoff = Float.T(
382 optional=True,
383 help='Set weight to zero, when noise level exceeds median by the '
384 'given cutoff factor.')
386 cutoff_exception_on_high_snr = Float.T(
387 optional=True,
388 help='Exclude from cutoff when max amplitude exceeds standard '
389 'deviation times this factor.')
391 def get_analyser(self):
392 return NoiseAnalyser(
393 nwindows=self.nwindows,
394 pre_event_noise_duration=self.pre_event_noise_duration,
395 check_events=self.check_events, phase_def=self.phase_def,
396 statistic=self.statistic, mode=self.mode, cutoff=self.cutoff,
397 cutoff_exception_on_high_snr=self.cutoff_exception_on_high_snr)
400__all__ = '''
401 NoiseAnalyser
402 NoiseAnalyserConfig
403'''.split()