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