Commit 805cf696 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

restructured

- scan is now called search
- config paths are now relative to config file location
- use resample if downsample_to fails
parent 3f608e71
......@@ -27,14 +27,14 @@ def str_to_time(s):
subcommand_descriptions = {
'init': 'create initial configuration file',
'scan': 'detect seismic events',
'search': 'detect seismic events',
'map-geometry': 'make station map',
'snuffle': 'snuffle'
}
subcommand_usages = {
'init': 'init',
'scan': 'scan <configfile> [options]',
'search': 'search <configfile> [options]',
'map-geometry': 'map-geometry <configfile> [options] <output.(png|pdf)',
'snuffle': 'snuffle <configfile>',
}
......@@ -51,7 +51,7 @@ usage = program_name + ''' <subcommand> [options] [--] <arguments> ...
Subcommands:
init %(init)s
scan %(scan)s
search %(search)s
map-geometry %(map_geometry)s
snuffle %(snuffle)s
......@@ -231,7 +231,7 @@ stackmax_path: 'stackmax'
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',
......@@ -268,7 +268,7 @@ def command_scan(args):
'--nworkers', dest='nworkers', metavar="N",
help='use N cpus in parallel')
parser, options, args = cl_parse('scan', args, setup=setup)
parser, options, args = cl_parse('search', args, setup=setup)
if len(args) != 1:
help_and_die(parser, 'missing argument')
......@@ -288,7 +288,7 @@ def command_scan(args):
else:
nparallel = None
lassie.scan(
lassie.search(
config,
override_tmin=tmin,
override_tmax=tmax,
......
import os.path as op
from collections import defaultdict
from pyrocko.guts import Object, String
from pyrocko.gf import Earthmodel1D
......@@ -25,8 +26,87 @@ 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]
__all__ = [
'LassieError',
'Earthmodel',
'CakeEarthmodel',
'HasPaths',
'Path',
]
import logging
import os.path as op
from pyrocko.guts import Object, String, Float, Timestamp, List, Bool
from pyrocko.guts import load
from pyrocko import model
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
from lassie.common import Earthmodel, HasPaths, Path, LassieError
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 = String.T(
stations_stationxml_path = Path.T(
optional=True,
help='stations file in StationXML format')
......@@ -29,10 +29,10 @@ class Config(Object):
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 = String.T(
events_path = Path.T(
optional=True,
help='limit processing to time windows around given events')
......@@ -60,13 +60,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 detections')
help='output is saved to this directory')
stackmax_path = String.T(
optional=True,
help='output path name for stacked traces')
save_figures = Bool.T(
default=False,
help='flag to activate saving of detection figures')
grid = grid.Grid.T(
optional=True,
......@@ -107,7 +107,7 @@ class Config(Object):
TPDef.T(),
help='list of tabulated phase definitions usable shifters')
cache_path = String.T(
cache_path = Path.T(
default='lassie_cache',
help='directory where lassie stores tabulated phases etc.')
......@@ -117,13 +117,36 @@ class Config(Object):
self._grid = None
self._events = None
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(),
......@@ -132,7 +155,7 @@ class Config(Object):
z=station.depth))
if self.stations_stationxml_path:
sx = fs.load_xml(filename=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(
......@@ -148,7 +171,8 @@ class Config(Object):
return None
if self._events is None:
self._events = model.load_events(self.events_path)
self._events = model.load_events(self.expand_path(
self.events_path))
return self._events
......@@ -196,9 +220,21 @@ 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 '.')
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__ = [
......
import os
import math
import logging
import os.path as op
......@@ -8,11 +7,12 @@ from collections import defaultdict
import numpy as num
from pyrocko import pile, trace, util, io
from pyrocko import pile, trace, util, io, model
from pyrocko.parstack import parstack
from pyrocko.guts import Object, Timestamp, String, Float
from lassie import common, plot, grid as gridmod, geo
from lassie.config import write_config
logger = logging.getLogger('lassie.core')
......@@ -23,6 +23,18 @@ class Detection(Object):
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()
......@@ -101,7 +113,7 @@ def zero_fill(trs, tmin, tmax):
return trs_out
def scan(
def search(
config,
override_tmin=None,
override_tmax=None,
......@@ -112,32 +124,26 @@ def scan(
stop_after_first=False,
nparallel=None):
if config.detections_path:
if op.exists(config.detections_path):
if force:
os.unlink(config.detections_path)
if config.stackmax_path:
if op.exists(config.stackmax_path):
shutil.rmtree(config.stackmax_path)
os.mkdir(config.stackmax_path)
else:
raise common.LassieError(
'detections file already exists: %s'
% config.detections_path)
util.ensuredirs(config.detections_path)
fp = config.expand_path
else:
if config.stackmax_path:
if op.exists(config.stackmax_path):
if force:
shutil.rmtree(config.stackmax_path)
else:
raise common.LassieError(
'stackmax directory already exists: %s' %
config.stackmax_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.ensuredir(run_path)
write_config(config, op.join(run_path, 'config.yaml'))
util.ensuredirs(config.stackmax_path)
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:
......@@ -148,12 +154,13 @@ def scan(
norm_map = gridmod.geometrical_normalization(grid, receivers)
for data_path in config.data_paths:
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(config.data_paths, fileformat='detect')
p = pile.make_pile(data_paths, fileformat='detect')
if p.is_empty():
raise common.LassieError('no usable waveforms found')
......@@ -203,15 +210,6 @@ def scan(
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)
p = pile.make_pile(config.data_paths, fileformat='detect')
if p.is_empty():
raise common.LassieError('no usable waveforms found')
check_data_consistency(p, config)
deltat_cf = max(p.deltats.keys())
......@@ -418,7 +416,7 @@ def scan(
tr_stackmax_indx = tr_stackmax.copy(data=False)
imaxs = num.zeros(len(frames))
imaxs = num.zeros(frames.shape[1], dtype=num.int32)
for iframe, frame in enumerate(frames.T):
imaxs[iframe] = num.argmax(frame)
......@@ -449,21 +447,32 @@ def scan(
logger.info('detection: %s' % str(detection))
if config.detections_path:
f = open(config.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'),
apeak,
latpeak, lonpeak, xpeak, ypeak, zpeak))
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'),
apeak,
latpeak, lonpeak, xpeak, ypeak, zpeak))
f.close()
f.close()
if show_detections:
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,
......@@ -474,7 +483,9 @@ def scan(
wmin, wmax,
pdata, trs, fmin, fmax, idetection,
grid_station_shift_max=shift_max,
movie=show_movie)
movie=show_movie,
show=show_detections,
save_filename=fn)
except AttributeError as e:
logger.warn(e)
......@@ -484,11 +495,7 @@ def scan(
tr_stackmax.chop(wmin, wmax)
tr_stackmax_indx.chop(wmin, wmax)
if config.stackmax_path:
io.save(
[tr_stackmax, tr_stackmax_indx],
os.path.join(
config.stackmax_path, 'trace_%(tmin_ms)s.mseed'))
io.save([tr_stackmax, tr_stackmax_indx], ifm_path_template)
logger.info('end processing time window group: %s - %s' % (
util.time_to_str(tmin_win),
......@@ -496,5 +503,5 @@ def scan(
__all__ = [
'scan',
'search',
]
......@@ -37,6 +37,11 @@ def points_coords(points, system=None):
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):
return num.broadcast_arrays(*[
num.asarray(x, dtype=num.float) for x in args])
......
import sys
import logging
from collections import defaultdict
import numpy as num
......@@ -13,6 +12,17 @@ logger = logging.getLogger('lassie.ifc')
guts_prefix = 'lassie'
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 using NSLC-id lists and/or lists of regular
......@@ -40,7 +50,7 @@ class TraceSelector(Object):
self.__class__.__name__, self.white_list, self.white_list_regex)
class IFC(Object):
class IFC(common.HasPaths):
'''Image function contribution.'''
name = String.T()
......@@ -152,14 +162,7 @@ class WavePacketIFC(IFC):
tr.set_ydata(fftconvolve(tr.get_ydata(), taper))
tr.set_ydata(num.maximum(tr.ydata, 0.0))
tr.shift(-(n/2.*tr.deltat))
try:
tr.downsample_to(deltat_cf, snap=True, demean=False)
except util.UnavailableDecimation as e:
logger.fatal('%s, in_deltat = %g, out_deltat = %g' % (
e, tr.deltat, deltat_cf))
sys.exit(1)
downsample(tr, deltat_cf)
trs_filt.append(tr)
trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])
......@@ -293,13 +296,7 @@ class OnsetIFC(IFC):
sumtr.set_ydata(sumtr.get_ydata() / normtr.get_ydata())
try:
sumtr.downsample_to(deltat_cf, snap=True, demean=False)
except util.UnavailableDecimation as e:
logger.error('%s, in_deltat = %g, out_deltat = %g' % (
e, tr.deltat, deltat_cf))
continue
downsample(sumtr, deltat_cf)
sumtr.chop(wmin - tpad_new, wmax + tpad_new)
......@@ -336,7 +333,10 @@ class TemplateMatchingIFC(IFC):
return tmaster
def get_template_origin(self):
event = model.load_one_event(self.template_event_path)
event = model.load_one_event(self.expand_path(
self.template_event_path))
origin = geo.Point(
lat=event.lat,
lon=event.lon,
......@@ -350,7 +350,8 @@ class TemplateMatchingIFC(IFC):
def extract_template(self, p):
markers = gui_util.load_markers(self.template_markers_path)
markers = gui_util.load_markers(self.expand_path(
self.template_markers_path))
def trace_selector_global(tr):
return True
......@@ -377,10 +378,7 @@ class TemplateMatchingIFC(IFC):
masters = {}
for xtr in master_traces:
tr = xtr.copy()
if self.downsample_rate is not None:
tr.downsample_to(
1./self.downsample_rate, demean=False, snap=True,
allow_upsample_max=4)
downsample(tr, 1./self.downsample_rate)
tr.highpass(4, self.fmin, demean=False)
tr.lowpass(4, self.fmax, demean=False)
......@@ -416,10 +414,7 @@ class TemplateMatchingIFC(IFC):
nslc = b.nslc_id
a = self.masters[nslc]
if self.downsample_rate is not None:
b.downsample_to(
1./self.downsample_rate, demean=False, snap=True,
allow_upsample_max=4)
downsample(b, 1./self.downsample_rate)
b.highpass(4, self.fmin, demean=False)
b.lowpass(4, self.fmax, demean=False)
......@@ -467,7 +462,7 @@ class ManualPickIFC(IFC):
def __init__(self, *args, **kwargs):
IFC.__init__(self, *args, **kwargs)
markers = pmarker.load_markers(self.picks_path)
markers = pmarker.load_markers(self.expand_path(self.picks_path))
nsl_to_index = {}
picked_index = []
picked_times = []
......@@ -551,8 +546,8 @@ class ManualPickIFC(IFC):
sumtr.set_ydata(
num.convolve(sumtr.get_ydata()/len(trs_by_nsl), taper))