Browse Source

Merge branch 'master' of https://gitext.gfz-potsdam.de/heimann/lassie

template_parallel
miili 4 years ago
parent
commit
3af729e02e
16 changed files with 1968 additions and 356 deletions
  1. +15
    -0
      README.md
  2. +117
    -27
      apps/lassie
  3. BIN
      doc/gx/snuffling_example.png
  4. +1
    -1
      setup.py
  5. +1
    -0
      src/__init__.py
  6. +114
    -0
      src/common.py
  7. +130
    -13
      src/config.py
  8. +412
    -133
      src/core.py
  9. BIN
      src/data/bark.wav
  10. +42
    -2
      src/geo.py
  11. +53
    -4
      src/grid.py
  12. +438
    -53
      src/ifc.py
  13. +183
    -88
      src/plot.py
  14. +3
    -32
      src/receiver.py
  15. +166
    -3
      src/shifter.py
  16. +293
    -0
      src/snuffling.py

+ 15
- 0
README.md View File

@ -189,3 +189,18 @@ Now run the detector with
lassie scan config.yaml
```
#### Scrutinize detections
After the detection has finished have a look at the results using Pyrocko's
[Snuffler](http://emolch.github.io/pyrocko/current/snuffler.html):
```bash
lassie snuffle config.yaml
```
![example snuffling](doc/gx/snuffling_example.png "Test with synthetic events")
Snuffler opens loading waveform data together with station meta data and
detections as event markers. Load a reference catalog to compare the detections
to and scrutinize detection performance at different detection thresholds,
interactively.

+ 117
- 27
apps/lassie View File

@ -27,14 +27,16 @@ def str_to_time(s):
subcommand_descriptions = {
'init': 'create initial configuration file',
'scan': 'detect seismic events',
'map-geometry': 'make station map'
'search': 'detect seismic events',
'map-geometry': 'make station map',
'snuffle': 'snuffle'
}
subcommand_usages = {
'init': 'init',
'scan': 'scan <configfile> [options]',
'map-geometry': 'map-geometry <configfile> [options]',
'search': 'search <configfile> [options]',
'map-geometry': 'map-geometry <configfile> [options] <output.(png|pdf)',
'snuffle': 'snuffle <configfile>',
}
subcommands = subcommand_descriptions.keys()
@ -49,8 +51,9 @@ usage = program_name + ''' <subcommand> [options] [--] <arguments> ...
Subcommands:
init %(init)s
scan %(scan)s
search %(search)s
map-geometry %(map_geometry)s
snuffle %(snuffle)s
To get further help and a list of available options for any subcommand run:
@ -67,9 +70,9 @@ def add_common_options(parser):
type='choice',
choices=('critical', 'error', 'warning', 'info', 'debug'),
default='info',
help ='set logger level to '
'"critical", "error", "warning", "info", or "debug". '
'Default is "%default".')
help='set logger level to '
'"critical", "error", "warning", "info", or "debug". '
'Default is "%default".')
def process_common_options(options):
@ -169,6 +172,10 @@ stations_path: '%(stations_path)s'
data_paths:
%(s_data_paths)s
## name template for Lassie's output directory. The placeholder
## "${config_name}" will be replaced with the basename of the config file.
run_path: '${config_name}.turd'
## Processing time interval (default: use time interval of available data)
# tmin: '2012-02-06 04:20:00'
# tmax: '2012-02-06 04:30:00'
@ -193,39 +200,90 @@ autogrid_radius_factor: 1.5
## Grid density factor used when automatically choosing a grid
autogrid_density_factor: 10.0
## Image function composition
## Composition of image function
image_function_contributions:
- !lassie.WavePacketIFC
- !lassie.OnsetIFC
name: 'P'
weight: 1.0
weight: 30.0
fmin: 1.0
fmax: 15.0
fsmooth_factor: 0.1
shifter: !lassie.VelocityShifter
velocity: 3500.
short_window: 1.0
window_ratio: 8.0
fsmooth: 0.2
fnormalize: 0.02
#shifter: !lassie.VelocityShifter
# velocity: 6000.
shifter: !lassie.CakePhaseShifter
timing: '{stored:p}'
earthmodel_id: 'swiss'
- !lassie.WavePacketIFC
name: 'S'
weight: 1.0
fmin: 1.0
fmax: 10.0
fsmooth_factor: 0.1
shifter: !lassie.VelocityShifter
velocity: 2500.
fmax: 8.0
fsmooth: 0.05
#shifter: !lassie.VelocityShifter
# velocity: 3300.
shifter: !lassie.CakePhaseShifter
# factor: 1.0
# offset: 1.0
timing: '{stored:s}'
earthmodel_id: 'swiss'
## Whether to divide image function frames by their mean value
sharpness_normalization: true
sharpness_normalization: false
## Threshold on detector function
detector_threshold: 100.0
## Output filename for detections
detections_path: 'detections.txt'
detector_threshold: 150.
## Whether to create a figure for every detection and save it in the output
## directory
save_figures: true
## Mapping of phase ID to phase definition in cake syntax (used e.g. in the
## CakePhaseShifter config sections)
tabulated_phases:
- !pf.TPDef
id: 'p'
definition: 'P,p'
- !pf.TPDef
id: 's'
definition: 'S,s'
## Mapping of earthmodel ID to the actual earth model in nd format (used in
## the CakePhaseShifter config sections)
earthmodels:
- !lassie.CakeEarthmodel
id: 'swiss'
earthmodel_1d: |2
0.0 5.53 3.10 2.75
2.0 5.53 3.10 2.75
2.0 5.80 3.25 2.75
5.0 5.80 3.25 2.75
5.0 5.83 3.27 2.75
8.0 5.83 3.27 2.75
8.0 5.95 3.34 2.8
13.0 5.95 3.34 2.8
13.0 5.96 3.34 2.8
22.0 5.96 3.34 2.8
22.0 6.53 3.66 2.8
30.0 6.53 3.66 2.8
30.0 7.18 4.03 3.3
40.0 7.18 4.03 3.3
40.0 7.53 4.23 3.3
50.0 7.53 4.23 3.3
50.0 7.83 4.39 3.3
60.0 7.83 4.39 3.3
60.0 8.15 4.57 3.3
120.0 8.15 4.57 3.3
''' % dict(
stations_path=stations_path,
s_data_paths=s_data_paths)
def command_scan(args):
def command_search(args):
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
@ -239,6 +297,11 @@ def command_scan(args):
'--show-movie', dest='show_movie', action='store_true',
help='show movie when showing detections')
parser.add_option(
'--show-window-traces', dest='show_window_traces',
action='store_true',
help='show preprocessed traces for every processing time window')
parser.add_option(
'--stop-after-first', dest='stop_after_first', action='store_true',
help='show plot for every detection found')
@ -253,7 +316,15 @@ def command_scan(args):
help='end of processing time window '
'(overrides config file settings)')
parser, options, args = cl_parse('scan', args, setup=setup)
parser.add_option(
'--nworkers', dest='nworkers', metavar="N",
help='use N cpus in parallel')
parser.add_option(
'--speak', dest='bark', action='store_true',
help='alert on detection of events')
parser, options, args = cl_parse('search', args, setup=setup)
if len(args) != 1:
help_and_die(parser, 'missing argument')
@ -268,14 +339,22 @@ def command_scan(args):
if options.tmax:
tmax = str_to_time(options.tmax)
lassie.scan(
if options.nworkers:
nparallel = int(options.nworkers)
else:
nparallel = None
lassie.search(
config,
override_tmin=tmin,
override_tmax=tmax,
force=options.force,
show_detections=options.show_detections,
show_movie=options.show_movie,
stop_after_first=options.stop_after_first)
show_window_traces=options.show_window_traces,
stop_after_first=options.stop_after_first,
nparallel=nparallel,
bark=options.bark)
except lassie.LassieError, e:
die(str(e))
@ -292,6 +371,17 @@ def command_map_geometry(args):
lassie.map_geometry(config, output_path)
def command_snuffle(args):
parser, options, args = cl_parse('snuffle', args)
if len(args) != 1:
help_and_die(parser, 'missing arguments')
config_path = args[0]
config = lassie.read_config(config_path)
lassie.snuffle(config)
if __name__ == '__main__':
usage_sub = 'fomosto %s [options]'


BIN
doc/gx/snuffling_example.png View File

Before After
Width: 660  |  Height: 665  |  Size: 65 KiB

+ 1
- 1
setup.py View File

@ -7,4 +7,4 @@ setup(
packages=['lassie'],
package_dir={'lassie': 'src'},
scripts=['apps/lassie'],
package_data={'lassie': []})
package_data={'lassie': ['data/*.wav']})

+ 1
- 0
src/__init__.py View File

@ -6,6 +6,7 @@ from lassie.ifc import * # noqa
from lassie.core import * # noqa
from lassie.plot import * # noqa
from lassie.config import * # noqa
from lassie.snuffling import * # noqa
__version__ = '0.2'

+ 114
- 0
src/common.py View File

@ -1,12 +1,29 @@
import os.path as op
from string import Template
from collections import defaultdict
from pyrocko.guts import Object, String
from pyrocko.gf import Earthmodel1D
guts_prefix = 'lassie'
def data_file(fn):
return op.join(op.split(__file__)[0], 'data', fn)
class LassieError(Exception):
pass
class Earthmodel(Object):
id = String.T()
class CakeEarthmodel(Earthmodel):
earthmodel_1d = Earthmodel1D.T()
def grouped_by(l, key):
d = defaultdict(list)
for x in l:
@ -14,6 +31,103 @@ def grouped_by(l, key):
return d
def xjoin(basepath, path):
if path is None and basepath is not None:
return basepath
elif op.isabs(path) or basepath is None:
return path
else:
return op.join(basepath, path)
def xrelpath(path, start):
if op.isabs(path):
return path
else:
return op.relpath(path, start)
class Path(String):
pass
class HasPaths(Object):
path_prefix = Path.T(optional=True)
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._basepath = None
self._parent_path_prefix = None
def set_basepath(self, basepath, parent_path_prefix=None):
self._basepath = basepath
self._parent_path_prefix = parent_path_prefix
for (prop, val) in self.T.ipropvals(self):
if isinstance(val, HasPaths):
val.set_basepath(
basepath, self.path_prefix or self._parent_path_prefix)
def get_basepath(self):
assert self._basepath is not None
return self._basepath
def change_basepath(self, new_basepath, parent_path_prefix=None):
assert self._basepath is not None
self._parent_path_prefix = parent_path_prefix
if self.path_prefix or not self._parent_path_prefix:
self.path_prefix = op.normpath(xjoin(xrelpath(
self._basepath, new_basepath), self.path_prefix))
for val in self.T.ivals(self):
if isinstance(val, HasPaths):
val.change_basepath(
new_basepath, self.path_prefix or self._parent_path_prefix)
self._basepath = new_basepath
def expand_path(self, path, extra=None):
assert self._basepath is not None
if extra is None:
def extra(path):
return path
path_prefix = self.path_prefix or self._parent_path_prefix
if path is None:
return None
elif isinstance(path, basestring):
return extra(
op.normpath(xjoin(self._basepath, xjoin(path_prefix, path))))
else:
return [
extra(
op.normpath(xjoin(self._basepath, xjoin(path_prefix, p))))
for p in path]
def expand_template(template, d):
try:
return Template(template).substitute(d)
except KeyError as e:
raise LassieError(
'invalid placeholder "%s" in template: "%s"' % (str(e), template))
except ValueError:
raise LassieError(
'malformed placeholder in template: "%s"' % template)
def bark():
import subprocess
subprocess.call(['aplay', data_file('bark.wav')])
__all__ = [
'LassieError',
'Earthmodel',
'CakeEarthmodel',
'HasPaths',
'Path',
]

+ 130
- 13
src/config.py View File

@ -1,33 +1,58 @@
import logging
import os.path as op
from pyrocko.guts import Object, String, Float, Timestamp, List, Bool, load
from pyrocko import model
from pyrocko.guts import Object, String, Float, Timestamp, List, Bool
from pyrocko import model, guts
from pyrocko.fdsn import station as fs
from pyrocko.gf import TPDef
from lassie import receiver, ifc, grid, geo
from lassie.common import Earthmodel, HasPaths, Path, LassieError, \
expand_template
guts_prefix = 'lassie'
logger = logging.getLogger('lassie.config')
class Config(Object):
class Config(HasPaths):
stations_path = String.T(
stations_path = Path.T(
optional=True,
help='stations file in Pyrocko format')
stations_stationxml_path = Path.T(
optional=True,
help='stations file in StationXML format')
receivers = List.T(
receiver.Receiver.T(),
help='receiver coordinates if not read from file')
data_paths = List.T(
String.T(),
Path.T(),
help='list of directories paths to search for data')
events_path = Path.T(
optional=True,
help='limit processing to time windows around given events')
event_time_window_factor = Float.T(
default=2.,
help='controls length of time windows for event-wise processing')
blacklist = List.T(
String.T(),
help='codes in the form NET.STA.LOC of receivers to be excluded')
whitelist = List.T(
String.T(),
help='codes in the form NET.STA.LOC of receivers to be included')
distance_max = Float.T(
optional=True,
help='receiver maximum distance from grid')
tmin = Timestamp.T(
optional=True,
help='beginning of time interval to be processed')
@ -36,9 +61,13 @@ class Config(Object):
optional=True,
help='end of time interval to be processed')
detections_path = String.T(
run_path = Path.T(
optional=True,
help='output file name for for detections')
help='output is saved to this directory')
save_figures = Bool.T(
default=False,
help='flag to activate saving of detection figures')
grid = grid.Grid.T(
optional=True,
@ -66,26 +95,99 @@ class Config(Object):
default=200.,
help='threshold on detector function')
fill_incomplete_with_zeros = Bool.T(
default=False,
help='fill incomplete trace time windows with zeros '
'(and let edge effects ruin your day)')
earthmodels = List.T(
Earthmodel.T(),
help='list of earthmodels usable in shifters')
tabulated_phases = List.T(
TPDef.T(),
help='list of tabulated phase definitions usable shifters')
cache_path = Path.T(
default='lassie_phases.cache',
help='directory where lassie stores tabulated phases etc.')
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._receivers = None
self._grid = None
self._events = None
self._config_name = 'untitled'
def set_config_name(self, config_name):
self._config_name = config_name
def expand_path(self, path):
def extra(path):
return expand_template(path, dict(
config_name=self._config_name))
return HasPaths.expand_path(self, path, extra=extra)
def get_events_path(self):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'events.list')
def get_ifm_dir_path(self):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'ifm')
def get_ifm_path_template(self):
return op.join(
self.get_ifm_dir_path(),
'%(station)s_%(tmin_ms)s.mseed')
def get_detections_path(self):
run_path = self.expand_path(self.run_path)
return op.join(run_path, 'detections.list')
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')
def get_receivers(self):
'''Aggregate receivers from different sources.'''
fp = self.expand_path
if self._receivers is None:
self._receivers = list(self.receivers)
if self.stations_path:
for station in model.load_stations(self.stations_path):
for station in model.load_stations(fp(self.stations_path)):
self._receivers.append(
receiver.Receiver(
codes=station.nsl(),
lat=station.lat,
lon=station.lon,
z=station.depth))
if self.stations_stationxml_path:
sx = fs.load_xml(filename=fp(self.stations_stationxml_path))
for station in sx.get_pyrocko_stations():
self._receivers.append(
receiver.Receiver(
codes=station.nsl(),
lat=station.lat,
lon=station.lon))
lon=station.lon,
z=station.depth))
return self._receivers
def get_events(self):
if self.events_path is None:
return None
if self._events is None:
self._events = model.load_events(self.expand_path(
self.events_path))
return self._events
def get_grid(self):
'''Get grid or make default grid.'''
@ -104,7 +206,7 @@ class Config(Object):
spacing = vmin / fsmooth_max / self.autogrid_density_factor
lat0, lon0, north, east, depth = geo.bounding_box_square(
*receiver.receivers_coords(receivers),
*geo.points_coords(receivers),
scale=self.autogrid_radius_factor)
self._grid = grid.Carthesian3DGrid(
@ -130,9 +232,24 @@ class Config(Object):
return self._grid
def read_config(file_path):
conf = load(filename=file_path)
return conf
def read_config(path):
config = guts.load(filename=path)
if not isinstance(config, Config):
raise LassieError('invalid Lassie configuration in file "%s"' % path)
config.set_basepath(op.dirname(path) or '.')
config.set_config_name(op.splitext(op.basename(path))[0])
return config
def write_config(config, path):
basepath = config.get_basepath()
dirname = op.dirname(path) or '.'
config.change_basepath(dirname)
guts.dump(config, filename=path)
config.change_basepath(basepath)
__all__ = [
'Config',


+ 412
- 133
src/core.py View File

@ -1,231 +1,510 @@
import sys
import os
import math
import logging
import os.path as op
import shutil
from collections import defaultdict
import numpy as num
from pyrocko import pile, trace, util, io
from pyrocko.parstack import parstack
from pyrocko import pile, trace, util, io, model
from pyrocko.parstack import parstack, argmax as pargmax
from pyrocko.guts import Object, Timestamp, String, Float
from lassie import common, plot, grid as gridmod
from lassie import common, plot, grid as gridmod, geo
from lassie.config import write_config
logger = logging.getLogger('lassie.core')
def scan(
class Detection(Object):
id = String.T()
time = Timestamp.T()
location = geo.Point.T()
ifm = Float.T()
def get_event(self):
lat, lon = geo.point_coords(self.location, system='latlon')
event = model.Event(
lat=lat,
lon=lon,
depth=self.location.z,
time=self.time,
name='%s (%g)' % (self.id, self.ifm))
return event
def check_data_consistency(p, config):
receivers = config.get_receivers()
nslc_ids = p.nslc_ids.keys()
nsl_ids = [nslc_id[:3] for nslc_id in nslc_ids]
r_ids = [r.codes for r in receivers]
r_not_in_p = []
t_not_in_r = []
to_be_blacklisted = []
for r in receivers:
if r.codes[:3] not in nsl_ids:
r_not_in_p.append(r.codes)
if '.'.join(r.codes[:3]) in config.blacklist:
to_be_blacklisted.append('.'.join(r.codes))
for nsl_id in nsl_ids:
if nsl_id not in r_ids:
t_not_in_r.append(nsl_id)
if len(r_not_in_p) != 0.:
logger.warn('Following receivers have no traces in data set:')
for nsl_id in r_not_in_p:
logger.warn(' %s' % '.'.join(nsl_id))
logger.warn('-' * 40)
if len(t_not_in_r) != 0.:
logger.warn('Following traces have no associated receivers:')
for nsl_id in t_not_in_r:
logger.warn(' %s' % '.'.join(nsl_id))
logger.warn('-' * 40)
if len(to_be_blacklisted):
logger.info('Blacklisted receivers:')
for code in to_be_blacklisted:
logger.info(' %s' % code)
logger.info('-' * 40)
if len(config.blacklist) and\
len(to_be_blacklisted) != len(config.blacklist):
logger.warn('Blacklist NSL codes that did not match any receiver:')
for code in config.blacklist:
if code not in to_be_blacklisted:
logger.warn(' %s' % code)
logger.info('-' * 40)
def zero_fill(trs, tmin, tmax):
trs = trace.degapper(trs)
d = defaultdict(list)
for tr in trs:
d[tr.nslc_id].append(tr)
trs_out = []
for nslc, trs_group in d.iteritems():
if not all(tr.deltat == trs_group[0].deltat for tr in trs_group):
logger.warn('inconsistent sample rate, cannot merge traces')
continue
if not all(tr.ydata.dtype == trs_group[0].ydata.dtype
for tr in trs_group):
logger.warn('inconsistent data type, cannot merge traces')
continue
tr_combi = trs_group[0].copy()
tr_combi.extend(tmin, tmax, fillmethod='zeros')
for tr in trs_group[1:]:
tr_combi.add(tr)
trs_out.append(tr_combi)
return trs_out
def search(
config,
override_tmin=None,
override_tmax=None,
show_detections=False,
show_movie=False,
show_window_traces=False,
force=False,
stop_after_first=False):
stop_after_first=False,
nparallel=None,
bark=False):
fp = config.expand_path
if config.detections_path:
if op.exists(config.detections_path):
if force:
os.unlink(config.detections_path)
else:
raise common.LassieError(
'detections file already exists: %s' % config.detections_path)
run_path = fp(config.run_path)
if op.exists(run_path):
if force:
shutil.rmtree(run_path)
else:
raise common.LassieError(
'run directory already exists: %s' %
run_path)
util.ensuredirs(config.detections_path)
util.ensuredir(run_path)
write_config(config, op.join(run_path, 'config.yaml'))
ifm_path_template = config.get_ifm_path_template()
detections_path = config.get_detections_path()
events_path = config.get_events_path()
figures_path_template = config.get_figures_path_template()
ifcs = config.image_function_contributions
for ifc in ifcs:
ifc.setup(config)
grid = config.get_grid()
receivers = config.get_receivers()
norm_map = gridmod.geometrical_normalization(grid, receivers)
ifcs = config.image_function_contributions
data_paths = fp(config.data_paths)
for data_path in fp(data_paths):
if not op.exists(data_path):
raise common.LassieError(
'waveform data path does not exist: %s' % data_path)
p = pile.make_pile(data_paths, fileformat='detect')
if p.is_empty():
raise common.LassieError('no usable waveforms found')
for ifc in ifcs:
ifc.prescan(p)
shift_tables = []
tshift_maxs = []
tshift_minmaxs = []
for ifc in ifcs:
shift_tables.append(ifc.shifter.get_table(grid, receivers))
tshift_maxs.append(shift_tables[-1].max())
tshift_minmaxs.append(num.nanmin(shift_tables[-1]))
tshift_minmaxs.append(num.nanmax(shift_tables[-1]))
fsmooth_min = min(ifc.get_fsmooth() for ifc in ifcs)
tshift_max = max(tshift_maxs) * 1.0
tshift_min = min(tshift_minmaxs)
tshift_max = max(tshift_minmaxs)
tinc = tshift_max * 10
tpeaksearch = (tshift_max - tshift_min) + 1.0 / fsmooth_min
tpad = max(ifc.get_tpad() for ifc in ifcs) + tshift_max
tpad = max(ifc.get_tpad() for ifc in ifcs) + \
(tshift_max - tshift_min) + tpeaksearch
fsmooth_min = min(ifc.get_fsmooth() for ifc in ifcs)
tinc = (tshift_max - tshift_min) * 10. + 3.0 * tpad
tavail = p.tmax - p.tmin
tinc = min(tinc, tavail - 2.0 * tpad)
if tinc <= 0:
raise common.LassieError(
'available waveforms too short \n'
'required: %g s\n'
'available: %g s\n' % (2.*tpad, tavail))
blacklist = set(tuple(s.split('.')) for s in config.blacklist)
whitelist = set(tuple(s.split('.')) for s in config.whitelist)
distances = grid.distances(receivers)
distances_to_grid = num.min(distances, axis=0)
distance_min = num.min(distances)
distance_max = num.max(distances)
station_index = dict(
(rec.codes, i) for (i, rec) in enumerate(receivers)
if rec.codes not in blacklist)
if rec.codes not in blacklist and (
not whitelist or rec.codes in whitelist) and (
config.distance_max is None or
distances_to_grid[i] <= config.distance_max))
for data_path in config.data_paths:
if not op.exists(data_path):
raise common.LassieError(
'waveform data path does not exist: %s' % data_path)
check_data_consistency(p, config)
p = pile.make_pile(config.data_paths, fileformat='detect')
if p.is_empty():
raise common.LassieError('no usable waveforms found')
deltat_cf = max(p.deltats.keys())
assert deltat_cf > 0.0
while True:
if not all(ifc.deltat_cf_is_available(deltat_cf * 2) for ifc in ifcs):
break
deltat_cf *= 2
logger.info('CF sampling interval (rate): %g s (%g Hz)' % (
deltat_cf, 1.0/deltat_cf))
ngridpoints = grid.size()
logger.info('number of grid points: %i' % ngridpoints)
logger.info('minimum source-receiver distance: %g m' % distance_min)
logger.info('maximum source-receiver distance: %g m' % distance_max)
logger.info('minimum travel-time: %g s' % tshift_min)
logger.info('maximum travel-time: %g s' % tshift_max)
idetection = 0
station_weights = {}
tmin = override_tmin or config.tmin or p.tmin + tpad
tmax = override_tmax or config.tmax or p.tmax - tpad
tmin = override_tmin or config.tmin
tmax = override_tmax or config.tmax
for trs in p.chopper(
tmin=tmin, tmax=tmax, tinc=tinc, tpad=tpad,
want_incomplete=False,
trace_selector=lambda tr: tr.nslc_id[:3] in station_index):
events = config.get_events()
twindows = []
if events is not None:
for ev in events:
if tmin <= ev.time <= tmax:
twindows.append((
ev.time + tshift_min - (tshift_max - tshift_min) *
config.event_time_window_factor,
ev.time + tshift_min + (tshift_max - tshift_min) *
config.event_time_window_factor))
if not trs:
continue
else:
twindows.append((tmin, tmax))
for iwindow_group, (tmin_win, tmax_win) in enumerate(twindows):
nwin = int(math.ceil((tmax_win - tmin_win) / tinc))
logger.info('start processing time window group %i/%i: %s - %s' % (
iwindow_group + 1, len(twindows),
util.time_to_str(tmin_win),
util.time_to_str(tmax_win)))
logger.info('number of time windows: %i' % nwin)
logger.info('time window length: %g s' % (tinc + 2.0*tpad))
logger.info('time window payload: %g s' % tinc)
logger.info('time window padding: 2 x %g s' % tpad)
logger.info('time window overlap: %g%%' % (
100.0*2.0*tpad / (tinc+2.0*tpad)))
iwin = -1
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):
iwin += 1
trs_ok = []
for tr in trs:
if tr.ydata.size == 0:
logger.warn(
'skipping empty trace: %s.%s.%s.%s' % tr.nslc_id)
continue
if not num.all(num.isfinite(tr.ydata)):
logger.warn(
'skipping trace because of invalid values: '
'%s.%s.%s.%s' % tr.nslc_id)
continue
logger.info('processing time window %s - %s' % (
util.time_to_str(trs[0].wmin),
util.time_to_str(trs[0].wmax)))
trs_ok.append(tr)
wmin = trs[0].wmin
wmax = trs[0].wmax
trs = trs_ok
frames = None
pdata = []
trs_debug = []
shift_maxs = []
for iifc, ifc in enumerate(ifcs):
dataset = ifc.preprocess(trs, wmin, wmax, tshift_max)
if not dataset:
if not trs:
continue
nstations_selected = len(dataset)
logger.info('processing time window %i/%i: %s - %s' % (
iwin + 1, nwin,
util.time_to_str(trs[0].wmin),
util.time_to_str(trs[0].wmax)))
nsls_selected, trs_selected = zip(*dataset)
wmin = trs[0].wmin
wmax = trs[0].wmax
trs_debug.extend(trs + list(trs_selected))
if config.fill_incomplete_with_zeros:
trs = zero_fill(trs, wmin - tpad, wmax + tpad)
deltat_cf = trs_selected[0].deltat
frames = None
pdata = []
trs_debug = []
shift_maxs = []
for iifc, ifc in enumerate(ifcs):
dataset = ifc.preprocess(
trs, wmin-tpeaksearch, wmax+tpeaksearch,
tshift_max - tshift_min, deltat_cf)
t0 = (wmin / deltat_cf) * deltat_cf
if not dataset:
continue
istations_selected = num.array(
[station_index[nsl] for nsl in nsls_selected], dtype=num.int)
nstations_selected = len(dataset)
arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
nsls_selected, trs_selected = zip(*dataset)
offsets = num.array(
[int(round((tr.tmin-t0) / deltat_cf)) for tr in trs_selected],
dtype=num.int32)
for tr in trs_selected:
tr.meta = {'tabu': True}
w = num.array(
[station_weights.get(nsl, 1.0) for nsl in nsls_selected],
dtype=num.float)
trs_debug.extend(trs + list(trs_selected))
t0 = (wmin / deltat_cf) * deltat_cf
weights = num.ones((ngridpoints, nstations_selected))
weights *= w[num.newaxis, :]
weights *= ifc.weight
istations_selected = num.array(
[station_index[nsl] for nsl in nsls_selected],
dtype=num.int)
shift_table = shift_tables[iifc]
arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
shifts = -num.round(
shift_table[:, istations_selected] /
deltat_cf).astype(num.int32)
offsets = num.array(
[int(round((tr.tmin-t0) / deltat_cf))
for tr in trs_selected],
dtype=num.int32)
pdata.append((list(trs_selected), shift_table, ifc))
w = ifc.get_weights(nsls_selected)
iwmin = int(round((wmin-t0) / deltat_cf))
iwmax = int(round((wmax-t0) / deltat_cf))
weights = num.ones((ngridpoints, nstations_selected))
weights *= w[num.newaxis, :]
weights *= ifc.weight
lengthout = iwmax - iwmin + 1
shift_table = shift_tables[iifc]
shift_maxs.append(num.max(-shifts) * deltat_cf)
ok = num.isfinite(shift_table[:, istations_selected])
bad = num.logical_not(ok)
frames, ioff = parstack(
arrays, offsets, shifts, weights, 0,
offsetout=iwmin,
lengthout=lengthout,
result=frames,
impl='openmp')
shifts = -num.round(
shift_table[:, istations_selected] /
deltat_cf).astype(num.int32)
shift_max = max(shift_maxs)
weights[bad] = 0.0
shifts[bad] = num.max(shifts[ok])
pdata.append((list(trs_selected), shift_table, ifc))
iwmin = int(round((wmin-tpeaksearch-t0) / deltat_cf))
iwmax = int(round((wmax+tpeaksearch-t0) / deltat_cf))
lengthout = iwmax - iwmin + 1
shift_maxs.append(num.max(-shifts) * deltat_cf)
frames, ioff = parstack(
arrays, offsets, shifts, weights, 0,
offsetout=iwmin,
lengthout=lengthout,
result=frames,
nparallel=nparallel,
impl='openmp')
shift_max = max(shift_maxs)
if config.sharpness_normalization:
frame_maxs = frames.max(axis=0)
frame_means = num.abs(frames).mean(axis=0)
frames *= (frame_maxs / frame_means)[num.newaxis, :]
frames *= norm_map[:, num.newaxis]
if config.sharpness_normalization:
frame_maxs = frames.max(axis=0)
frame_means = num.abs(frames).mean(axis=0)
frames *= (frame_maxs / frame_means)[num.newaxis, :]
frames *= norm_map[:, num.newaxis]
frame_maxs = frames.max(axis=0)
tmin_frames = t0 + ioff * deltat_cf
tmin_frames = t0 + ioff * deltat_cf
tr_stackmax = trace.Trace(
'', 'SMAX', '', '',
tmin=tmin_frames,
deltat=deltat_cf,
ydata=frame_maxs)
tr_stackmax = trace.Trace(
'', 'SMAX', '', '',
tmin=tmin_frames,
deltat=deltat_cf,
ydata=frame_maxs)
tr_stackmax.meta = {'tabu': True}
trs_debug.append(tr_stackmax)
trs_debug.append(tr_stackmax)
# trace.snuffle(trs_debug)
if show_window_traces:
trace.snuffle(trs_debug)
ydata_window = tr_stackmax.chop(wmin, wmax, inplace=False).get_ydata()
ydata_window = tr_stackmax.chop(
wmin, wmax, inplace=False).get_ydata()
logger.info('CF stats: min %g, max %g, median %g' % (
num.min(ydata_window),
num.max(ydata_window),
num.median(ydata_window)))
logger.info('CF stats: min %g, max %g, median %g' % (
num.min(ydata_window),
num.max(ydata_window),
num.median(ydata_window)))
tpeaks, apeaks = tr_stackmax.peaks(
config.detector_threshold, shift_max + 1.0/fsmooth_min)
tpeaks, apeaks = zip(*[(tpeak, apeak) for (tpeak, apeak) in zip(
*tr_stackmax.peaks(config.detector_threshold, tpeaksearch)) if
wmin <= tpeak and tpeak < wmax]) or ([], [])
for (tpeak, apeak) in zip(tpeaks, apeaks):
if not (wmin <= tpeak and tpeak < wmax):
continue
tr_stackmax_indx = tr_stackmax.copy(data=False)
imaxs = pargmax(frames)
tr_stackmax_indx.set_ydata(imaxs.astype(num.int32))
tr_stackmax_indx.set_location('i')
logger.info('detection: %s %g' % (
util.time_to_str(tpeak),
apeak))
for (tpeak, apeak) in zip(tpeaks, apeaks):
iframe = int(round(((tpeak-t0) - ioff*deltat_cf) / deltat_cf))
frame = frames[:, iframe]
imax = num.argmax(frame)
iframe = int(round(((tpeak-t0) - ioff*deltat_cf) / deltat_cf))
frame = frames[:, iframe]
imax = num.argmax(frame)
latpeak, lonpeak, xpeak, ypeak, zpeak = grid.index_to_location(imax)
latpeak, lonpeak, xpeak, ypeak, zpeak = \
grid.index_to_location(imax)
idetection += 1
idetection += 1
if config.detections_path:
f = open(config.detections_path, 'a')
detection = Detection(
id='%06i' % idetection,
time=tpeak,
location=geo.Point(
lat=float(latpeak),
lon=float(lonpeak),
x=float(xpeak),
y=float(ypeak),
z=float(zpeak)),
ifm=float(apeak))
if bark:
common.bark()
logger.info('detection: %s' % str(detection))
f = open(detections_path, 'a')
f.write('%06i %s %g %g %g %g %g %g\n' % (
idetection,
util.time_to_str(tpeak, format='%Y-%m-%d %H:%M:%S.6FRAC'),
util.time_to_str(
tpeak,
format='%Y-%m-%d %H:%M:%S.6FRAC'),
apeak,
latpeak, lonpeak, xpeak, ypeak, zpeak))
f.close()
if show_detections:
fmin = min(ifc.fmin for ifc in ifcs)
fmax = min(ifc.fmax for ifc in ifcs)
plot.plot_detection(
grid, receivers, frames, tmin_frames,
deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak, zpeak,
tr_stackmax, tpeaks, apeaks, config.detector_threshold,
pdata, trs, fmin, fmax, idetection,
grid_station_shift_max=shift_max,
movie=show_movie)
ev = detection.get_event()
f = open(events_path, 'a')
model.dump_events([ev], stream=f)
f.close()
if show_detections or config.save_figures:
fmin = min(ifc.fmin for ifc in ifcs)
fmax = min(ifc.fmax for ifc in ifcs)
fn = figures_path_template % {
'id': idetection,
'format': 'png'}
util.ensuredirs(fn)
try:
plot.plot_detection(
grid, receivers, frames, tmin_frames,
deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak,
zpeak,
tr_stackmax, tpeaks, apeaks,
config.detector_threshold,
wmin, wmax,
pdata, trs, fmin, fmax, idetection,
grid_station_shift_max=shift_max,
movie=show_movie,
show=show_detections,
save_filename=fn)
except AttributeError as e:
logger.warn(e)
if stop_after_first:
return
tr_stackmax.chop(wmin, wmax)
tr_stackmax_indx.chop(wmin, wmax)
if stop_after_first:
return
io.save([tr_stackmax, tr_stackmax_indx], ifm_path_template)
tr_stackmax.chop(wmin, wmax)
io.save([tr_stackmax], 'stackmax/trace_%(tmin)s.mseed')
logger.info('end processing time window group: %s - %s' % (
util.time_to_str(tmin_win),
util.time_to_str(tmax_win)))
__all__ = [
'scan',
'search',
]

BIN
src/data/bark.wav View File


+ 42
- 2
src/geo.py View File

@ -1,5 +1,45 @@
import numpy as num
from pyrocko import orthodrome as od
from pyrocko.guts import Object, Float
class Point(Object):
lat = Float.T(default=0.0)
lon = Float.T(default=0.0)
x = Float.T(default=0.0)
y = Float.T(default=0.0)
z = Float.T(default=0.0)
@property
def vec(self):
return self.lat, self.lon, self.x, self.y, self.z
def points_coords(points, system=None):
a = num.zeros((len(points), 5))
for ir, r in enumerate(points):
a[ir, :] = r.vec
lats, lons, xs, ys, zs = a.T
if system is None:
return (lats, lons, xs, ys, zs)
elif system == 'latlon':
return od.ne_to_latlon(lats, lons, xs, ys)
elif system[0] == 'ne':
lat0, lon0 = system[1:3]
if num.all(lats == lat0) and num.all(lons == lon0):
return xs, ys
else:
elats, elons = od.ne_to_latlon(lats, lons, xs, ys)
return od.latlon_to_ne_numpy(lat0, lon0, elats, elons)
def point_coords(point, system=None):
coords = points_coords([point], system=system)
return [x[0] for x in coords]
def float_array_broadcast(*args):
@ -72,8 +112,8 @@ def bounding_box(lat, lon, north, east, depth, scale=1.0):
def bounding_box_square(lat, lon, north, east, depth, scale=1.0):
lat, lon, (north_min, north_max), (east_min, east_max), (depth_min, depth_max) = \
bounding_box(lat, lon, north, east, depth)
lat, lon, (north_min, north_max), (east_min, east_max), \
(depth_min, depth_max) = bounding_box(lat, lon, north, east, depth)
dnorth = north_max - north_min
mnorth = 0.5 * (north_min + north_max)


+ 53
- 4
src/grid.py View File

@ -3,7 +3,7 @@ import numpy as num
from pyrocko.guts import Object, Float
from pyrocko import orthodrome as od
from lassie import receiver
from lassie import geo
guts_prefix = 'lassie'
@ -58,17 +58,39 @@ class Carthesian3DGrid(Grid):
return self._coords
def depths(self):
nx, ny, _ = self._shape()
_, _, z = self._get_coords()
return num.repeat(z, nx*ny)
def index_to_location(self, i):
nx, ny, nz = self._shape()
iz, iy, ix = num.unravel_index(i, (nz, ny, nx))
x, y, z = self._get_coords()
return self.lat, self.lon, x[ix], y[iy], z[iz]
def lateral_distances(self, receivers):
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
rx, ry = geo.points_coords(
receivers, system=('ne', self.lat, self.lon))
na = num.newaxis
nr = len(receivers)
distances = num.sqrt(
(x[na, na, :, na] - rx[na, na, na, :])**2 +
(y[na, :, na, na] - ry[na, na, na, :])**2).reshape((nx*ny*nr))
return num.tile(distances, nz).reshape((nx*ny*nz, nr))
def distances(self, receivers):
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
rx, ry = receiver.receivers_coords(
rx, ry = geo.points_coords(
receivers, system=('ne', self.lat, self.lon))
rz = num.array([r.z for r in receivers], dtype=num.float)
@ -84,6 +106,12 @@ class Carthesian3DGrid(Grid):
return distances
def distance_max(self):
return math.sqrt(
(self.xmax - self.xmin)**2 +
(self.ymax - self.ymin)**2 +
(self.zmax - self.zmin)**2)
def surface_points(self, system='latlon'):
x, y, z = self._get_coords()
xs = num.tile(x, y.size)
@ -108,7 +136,10 @@ class Carthesian3DGrid(Grid):
amax=None,
z_slice=None,
cmap=None,
system='latlon'):
system='latlon',
shading='gouraud',
units=1.,
artists=[]):
if system == 'latlon':
assert False, 'not implemented yet'
@ -127,7 +158,21 @@ class Carthesian3DGrid(Grid):
else:
a2d = num.max(a3d, axis=0)
axes.pcolormesh(y, x, a2d.T, vmin=amin, vmax=amax, cmap=cmap)
if artists:
if shading == 'gouraud':
artists[0].set_array(a2d.T.ravel())
elif shading == 'flat':
artists[0].set_array(a2d.T[:-1, :-1].ravel())
else:
assert False, 'unknown shading option'
return artists
else:
return [
axes.pcolormesh(
y/units, x/units, a2d.T,
vmin=amin, vmax=amax, cmap=cmap, shading=shading)]
def geometrical_normalization(grid, receivers):
@ -154,11 +199,15 @@ def geometrical_normalization(grid, receivers):
nexpect = math.pi * ((dist + delta_ring)**2 - dist**2) /\
delta_grid**2
# nexpect = math.pi * ((dist + delta_ring)**3 - dist**3) /\
# delta_grid**3
norm_map[indices] += indices.size / nexpect
dist += delta_ring
norm_map /= nstations
return norm_map


+ 438
- 53
src/ifc.py View File

@ -1,30 +1,129 @@
import sys
import logging
from collections import defaultdict
import numpy as num
from pyrocko.guts import Object, String, Float
from pyrocko import trace, autopick, util
from lassie import shifter, common
from scipy.signal import fftconvolve
from pyrocko.guts import Object, String, Float, Bool, StringChoice, List, Dict
from pyrocko import trace, autopick, util, model, gui_util
from pyrocko import marker as pmarker
from lassie import shifter, common, geo
logger = logging.getLogger('lassie.ifc')
guts_prefix = 'lassie'
class IFC(Object):
def downsample(tr, deltat):
try:
tr.downsample_to(
deltat, demean=False, snap=True,
allow_upsample_max=4)
except util.UnavailableDecimation:
logger.warn('using resample instead of decimation')
tr.resample(deltat)
class TraceSelector(Object):
'''
Filter traces used in an IFC by NSLC-id lists and/or lists of regular
expressions.
'''
white_list = List.T(
optional=True,
default=[],
help='list of NSLC ids'
)
white_list_regex = List.T(
String.T(
default=[],
optional=True,
),
help='list of regular expressions'
)
def __call__(self, trs):
matched = []
for tr in trs:
nslc = '.'.join(tr.nslc_id)
if self.white_list and nslc in self.white_list:
matched.append(tr)
continue
if self.white_list_regex and util.match_nslc(
self.white_list_regex, nslc):
matched.append(tr)
return matched
def __str__(self):
return '%s:\n wl: %s\n wlre: %s\n' % (
self.__class__.__name__, self.white_list, self.white_list_regex)
class IFC(common.HasPaths):
'''Image function contribution.'''
name = String.T()
weight = Float.T(default=1.0)
weight = Float.T(
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')
fmin = Float.T()
fmax = Float.T()
shifter = shifter.Shifter.T()
shifter = shifter.Shifter.T(optional=True)
trace_selector = TraceSelector.T(
optional=True, help='select traces to be treated by this IFC')
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self.shifter.t_tolerance = 1./(self.fmax * 2.)
def setup(self, config):
if self.shifter:
self.shifter.setup(config)
def get_table(self, grid, receivers):
return self.shifter.get_table(grid, receivers)
def get_tpad(self):
return 4. / self.fmin
def preprocess(self, trs, wmin, wmax, tpad_new):
def get_fsmooth(self):
return self.fmin
def prescan(self, p):
pass
def deltat_cf_is_available(self, deltat_cf):
return False
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
return weights