Commit afc9db0c authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

improved time window determination and reporting

parent 0c84b7ae
import os
import math
import logging
import os.path as op
import shutil
......@@ -139,20 +140,23 @@ def scan(
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(num.nanmax(shift_tables[-1]))
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)
tpeaksearch = tshift_max + 1.0 / fsmooth_min
tpeaksearch = (tshift_max - tshift_min) + 1.0 / fsmooth_min
tpad = max(ifc.get_tpad() for ifc in ifcs) + tshift_max + tpeaksearch
tpad = max(ifc.get_tpad() for ifc in ifcs) + \
(tshift_max - tshift_min) + tpeaksearch
tinc = tshift_max * 10. + 3.0 * tpad
tinc = (tshift_max - tshift_min) * 10. + 3.0 * tpad
tavail = p.tmax - p.tmin
tinc = min(tinc, tavail - 2.0 * tpad)
......@@ -168,6 +172,9 @@ def scan(
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 and (
......@@ -195,16 +202,23 @@ def scan(
deltat_cf *= 2
logger.info('sampling interval: %g s' % deltat_cf)
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
tmax = override_tmax or config.tmax
tmin = override_tmin or config.tmin or p.tmin + tpad
tmax = override_tmax or config.tmax or p.tmax - tpad
events = config.get_events()
twindows = []
......@@ -212,18 +226,38 @@ def scan(
for ev in events:
if tmin <= ev.time <= tmax:
twindows.append((
ev.time - tshift_max * config.event_time_window_factor,
ev.time + tshift_max * config.event_time_window_factor))
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))
else:
twindows.append((tmin, tmax))
for tmin_win, tmax_win in twindows:
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:
......@@ -246,7 +280,8 @@ def scan(
if not trs:
continue
logger.info('processing time window %s - %s' % (
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)))
......@@ -263,7 +298,7 @@ def scan(
for iifc, ifc in enumerate(ifcs):
dataset = ifc.preprocess(
trs, wmin-tpeaksearch, wmax+tpeaksearch,
tshift_max, deltat_cf)
tshift_max - tshift_min, deltat_cf)
if not dataset:
continue
......@@ -425,6 +460,10 @@ def scan(
os.path.join(
config.stackmax_path, 'trace_%(tmin_ms)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',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment