Compare commits

...

11 Commits

  1. 3
      apps/lassie
  2. 14
      src/config.py
  3. 82
      src/core.py
  4. 83
      src/ifc.py
  5. 2
      src/plot.py
  6. 2
      src/snuffling.py

3
apps/lassie

@ -4,6 +4,9 @@ import sys
import logging
from optparse import OptionParser
import matplotlib
matplotlib.use('agg')
from pyrocko import util
import lassie

14
src/config.py

@ -6,7 +6,7 @@ from pyrocko import model, guts
from pyrocko.io import stationxml
from pyrocko.gf import TPDef
from lassie import receiver, ifc, grid, geo
from lassie import receiver, ifc, grid as lgrid, geo
from lassie.common import Earthmodel, HasPaths, Path, LassieError, \
expand_template
@ -69,7 +69,7 @@ class Config(HasPaths):
default=False,
help='flag to activate saving of detection figures')
grid = grid.Grid.T(
grid = lgrid.Grid.T(
optional=True,
help='definition of search grid, if not given, a default grid is '
'chosen')
@ -173,6 +173,10 @@ class Config(HasPaths):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'detections.list')
def get_markers_path_template(self):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'markers', '%(id)06i.markers')
def get_figures_path_template(self):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'figures', 'detection_%(id)06i.%(format)s')
@ -194,7 +198,9 @@ class Config(HasPaths):
z=station.depth))
if self.stations_stationxml_path:
sx = stationxml.load_xml(filename=fp(self.stations_stationxml_path))
sx = stationxml.load_xml(
filename=fp(self.stations_stationxml_path))
for station in sx.get_pyrocko_stations():
self._receivers.append(
receiver.Receiver(
@ -238,7 +244,7 @@ class Config(HasPaths):
*geo.points_coords(receivers),
scale=self.autogrid_radius_factor)
self._grid = grid.Carthesian3DGrid(
self._grid = lgrid.Carthesian3DGrid(
lat=lat0,
lon=lon0,
xmin=north[0],

82
src/core.py

@ -10,6 +10,7 @@ import numpy as num
from pyrocko import pile, trace, util, io, model
from pyrocko.parstack import parstack, argmax as pargmax
from pyrocko.gui import marker
from pyrocko.guts import Object, Timestamp, String, Float
@ -146,6 +147,7 @@ def search(
detections_path = config.get_detections_path()
events_path = config.get_events_path()
figures_path_template = config.get_figures_path_template()
markers_path_template = config.get_markers_path_template()
config.setup_image_function_contributions()
ifcs = config.image_function_contributions
@ -184,9 +186,10 @@ def search(
tpeaksearch = config.detector_tpeaksearch
else:
tpeaksearch = (tshift_max - tshift_min) + 1.0 / fsmooth_min
logger.info('using automatic detector_tpeaksearch: %g' % tpeaksearch)
tpad = max(ifc.get_tpad() for ifc in ifcs) + \
(tshift_max - tshift_min) + tpeaksearch
(tshift_max - tshift_min) + tpeaksearch / 2.0
tinc = (tshift_max - tshift_min) * 10. + 3.0 * tpad
tavail = p.tmax - p.tmin
@ -207,7 +210,7 @@ def search(
distance_min = num.min(distances)
distance_max = num.max(distances)
station_index = dict(
receiver_index = dict(
(rec.codes, i) for (i, rec) in enumerate(receivers)
if rec.codes not in blacklist and (
not whitelist or rec.codes in whitelist) and (
@ -273,9 +276,11 @@ def search(
iwin = -1
for trs in p.chopper(tmin=tmin_win, tmax=tmax_win, tinc=tinc, tpad=tpad,
for trs in p.chopper(
tmin=tmin_win, tmax=tmax_win, tinc=tinc, tpad=tpad,
want_incomplete=config.fill_incomplete_with_zeros,
trace_selector=lambda tr: tr.nslc_id[:3] in station_index):
trace_selector=lambda tr: tr.nslc_id[:3] in receiver_index):
iwin += 1
trs_ok = []
for tr in trs:
@ -311,8 +316,8 @@ def search(
trs = zero_fill(trs, wmin - tpad, wmax + tpad)
t0 = math.floor(wmin / deltat_cf) * deltat_cf
iwmin = int(round((wmin-tpeaksearch-t0) / deltat_cf))
iwmax = int(round((wmax+tpeaksearch-t0) / deltat_cf))
iwmin = int(round((wmin-tpeaksearch/2.-t0) / deltat_cf))
iwmax = int(round((wmax+tpeaksearch/2.-t0) / deltat_cf))
lengthout = iwmax - iwmin + 1
pdata = []
@ -320,13 +325,11 @@ def search(
parstack_params = []
for iifc, ifc in enumerate(ifcs):
dataset = ifc.preprocess(
trs, wmin-tpeaksearch, wmax+tpeaksearch,
trs, wmin-tpeaksearch/2., wmax+tpeaksearch/2.,
tshift_max - tshift_min, deltat_cf)
if not dataset:
continue
nstations_selected = len(dataset)
nsls_selected, trs_selected = zip(*dataset)
for tr in trs_selected:
@ -334,10 +337,13 @@ def search(
trs_debug.extend(trs + list(trs_selected))
istations_selected = num.array(
[station_index[nsl] for nsl in nsls_selected],
ireceivers_selected = num.array(
[receiver_index[nsl] for nsl in nsls_selected],
dtype=num.int)
receivers_selected = [
receivers[irec] for irec in ireceivers_selected]
arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
offsets = num.array(
@ -345,26 +351,22 @@ def search(
for tr in trs_selected],
dtype=num.int32)
w = ifc.get_weights(nsls_selected)
weights = num.ones((ngridpoints, nstations_selected))
weights *= w[num.newaxis, :]
weights *= ifc.weight
shift_table = shift_tables[iifc]
ok = num.isfinite(shift_table[:, istations_selected])
shifts_selected = shift_table[:, ireceivers_selected]
ok = num.isfinite(shifts_selected)
bad = num.logical_not(ok)
shifts = -num.round(
shift_table[:, istations_selected] /
shift_table[:, ireceivers_selected] /
deltat_cf).astype(num.int32)
weights[bad] = 0.0
shifts[bad] = num.max(shifts[ok])
pdata.append((list(trs_selected), shift_table, ifc))
weights = ifc.get_weights(
distances[:, ireceivers_selected],
receivers_selected,
bad)
pdata.append((list(trs_selected), shift_table, ifc))
parstack_params.append((arrays, offsets, shifts, weights))
if config.stacking_blocksize is not None:
@ -439,9 +441,11 @@ def search(
num.max(ydata_window),
num.median(ydata_window)))
tpeaks, apeaks = list(zip(*[(tpeak, apeak) for (tpeak, apeak) in zip(
*tr_stackmax.peaks(config.detector_threshold, tpeaksearch)) if
wmin <= tpeak and tpeak < wmax])) or ([], [])
tpeaks, apeaks = list(
zip(*[(tpeak, apeak) for (tpeak, apeak) in zip(
*tr_stackmax.peaks2(
config.detector_threshold, tpeaksearch)) if
wmin <= tpeak and tpeak < wmax])) or ([], [])
tr_stackmax_indx = tr_stackmax.copy(data=False)
tr_stackmax_indx.set_ydata(frame_argmaxs.astype(num.int32))
@ -489,6 +493,26 @@ def search(
model.dump_events([ev], stream=f)
f.close()
markers = [marker.EventMarker(ev)]
for _, shift_table, ifc in pdata:
for ireceiver, receiver in enumerate(receivers):
tshift = shift_table[imax, ireceiver]
if num.isfinite(tshift):
markers.append(
marker.PhaseMarker(
[receiver.codes + ('*',)],
tmin=tpeak + tshift,
tmax=tpeak + tshift,
event=ev,
phasename=ifc.name))
fn = markers_path_template % {
'id': idetection}
util.ensuredirs(fn)
marker.save_markers(markers, fn)
if show_detections or config.save_figures:
fmin = min(ifc.fmin for ifc in ifcs)
fmax = min(ifc.fmax for ifc in ifcs)
@ -507,10 +531,10 @@ def search(
else:
iframe_min = max(
0,
int(round(iframe - tpeaksearch/deltat_cf)))
int(round(iframe - tpeaksearch/2./deltat_cf)))
iframe_max = min(
lengthout-1,
int(round(iframe + tpeaksearch/deltat_cf)))
int(round(iframe + tpeaksearch/2./deltat_cf)))
ipsize = iframe_max - iframe_min + 1
frames_p = num.zeros((ngridpoints, ipsize))
@ -546,7 +570,7 @@ def search(
config.detector_threshold,
wmin, wmax,
pdata, trs, fmin, fmax, idetection,
tpeaksearch,
tpeaksearch/2.,
movie=show_movie,
show=show_detections,
save_filename=fn)

83
src/ifc.py

@ -3,7 +3,8 @@ import logging
import numpy as num
from collections import defaultdict
from scipy.signal import fftconvolve
from pyrocko.guts import Object, String, Float, Bool, StringChoice, List, Dict
from pyrocko.guts import Object, String, Float, Bool, StringChoice, List, \
Dict, Int
from pyrocko import trace, autopick, util, model
from pyrocko.gui import util as gui_util
from pyrocko import marker as pmarker
@ -25,6 +26,48 @@ def downsample(tr, deltat):
tr.resample(deltat)
class Weighting(Object):
def apply_weights(self, weights, distances, receivers):
raise NotImplemented('should be implemented in subclass')
class StationLimit(Weighting):
nstations = Int.T()
def apply_weights(self, weights, distances, receivers):
ii = num.argsort(distances, axis=1)
for ig in range(weights.shape[0]):
is_far = num.cumsum(weights[ig, ii[ig, :]] != 0) > self.nstations
weights[ig, ii[ig, is_far]] = 0.0
class DistanceDependent(Weighting):
distance_max = Float.T()
def apply_weights(self, weights, distances, receivers):
weights *= num.maximum(-distances/self.distance_max + 1.0, 0.0)
class StationWeighting(Weighting):
weights = Dict.T(
String.T(help='NSL regular expression identifying stations'),
Float.T(help='weighting factor'),
optional=True,
help='weight selected traces')
def apply_weights(self, weights, distances, receivers):
selectors = self.weights.keys()
for irec, rec in enumerate(receivers):
for selector in selectors:
if util.match_nslc(selector, rec.codes):
logger.debug('receiver weight for "%s" is: %s' % (
'.'.join(rec.codes), self.weights[selector]))
weights[:, irec] *= self.weights[selector]
break
class TraceSelector(Object):
'''
Filter traces used in an IFC by NSLC-id lists and/or lists of regular
@ -71,11 +114,7 @@ class IFC(common.HasPaths):
default=1.0,
help='global weight for this IFC')
weights = Dict.T(
String.T(help='NSL regular expression identifying stations'),
Float.T(help='weighting factor'),
optional=True,
help='weight selected traces')
weightings = List.T(Weighting.T())
fmin = Float.T()
fmax = Float.T()
@ -110,21 +149,13 @@ class IFC(common.HasPaths):
def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
pass
def get_weights(self, nsls):
if self.weights is None:
return num.ones(len(nsls), dtype=num.float)
else:
weights = num.empty(len(nsls))
selectors = self.weights.keys()
for insl, nsl in enumerate(nsls):
weights[insl] = 1.
for selector in selectors:
if util.match_nslc(selector, nsl):
weights[insl] = self.weights[selector]
break
def get_weights(self, distances, receivers, bad):
w = num.ones(distances.shape, dtype=num.float) * self.weight
w[bad] = 0.0
for weighting in self.weightings:
weighting.apply_weights(w, distances, receivers)
return weights
return w
class WavePacketIFC(IFC):
@ -139,6 +170,7 @@ class WavePacketIFC(IFC):
fsmooth = Float.T(optional=True)
fsmooth_factor = Float.T(default=0.1)
pre_downsample_deltat = Float.T(optional=True)
def get_tpad(self):
return 4. / self.get_fsmooth()
@ -164,6 +196,10 @@ class WavePacketIFC(IFC):
trs_filt = []
for orig_tr in trs:
tr = orig_tr.copy()
if self.pre_downsample_deltat is not None:
tr.downsample_to(self.pre_downsample_deltat)
tr.bandpass(4, self.fmin, self.fmax, demean=True)
tr.ydata = tr.ydata**2
@ -218,6 +254,7 @@ class OnsetIFC(IFC):
fsmooth_factor = Float.T(optional=True)
fnormalize = Float.T(optional=True)
fnormalize_factor = Float.T(optional=True)
pre_downsample_deltat = Float.T(optional=True)
def get_tpad(self):
return 2. / self.get_fnormalize()
@ -257,6 +294,9 @@ class OnsetIFC(IFC):
trs_filt = []
for orig_tr in trs:
tr = orig_tr.copy()
if self.pre_downsample_deltat is not None:
tr.downsample_to(self.pre_downsample_deltat)
tr.highpass(4, self.fmin, demean=True)
tr.lowpass(4, self.fmax, demean=False)
tr.ydata = tr.ydata**2
@ -592,6 +632,9 @@ class ManualPickIFC(IFC):
__all__ = [
'Weighting',
'StationLimit',
'DistanceDependent',
'IFC',
'WavePacketIFC',
'OnsetIFC',

2
src/plot.py

@ -389,6 +389,8 @@ def plot_detection(
if show:
plt.show()
else:
plt.close()
del ani

2
src/snuffling.py

@ -132,7 +132,7 @@ class LassieSnuffling(Snuffling):
tr_i = [None] * len(tr_smax)
for tr_i, tr_stackmax in zip(tr_i, tr_smax):
tpeaks, apeaks = tr_stackmax.peaks(
tpeaks, apeaks = tr_stackmax.peaks2(
self.detector_threshold, self.tsearch)
if self.level_trace:
ltrace = tr_stackmax.copy(data=False)

Loading…
Cancel
Save