Browse Source

initial import

tutorial
Sebastian Heimann 5 years ago
parent
commit
b24dc1caec
  1. 283
      apps/lassie
  2. 148
      scratch/lassie.svg
  3. BIN
      scratch/lassie1.png
  4. 144
      scratch/lassie1.svg
  5. 10
      setup.py
  6. 8
      src/__init__.py
  7. 19
      src/common.py
  8. 138
      src/config.py
  9. 220
      src/core.py
  10. 113
      src/geo.py
  11. 168
      src/grid.py
  12. 182
      src/ifc.py
  13. 305
      src/plot.py
  14. 42
      src/receiver.py
  15. 24
      src/shifter.py

283
apps/lassie

@ -0,0 +1,283 @@
#!/usr/bin/env python
import sys
import logging
from optparse import OptionParser
from pyrocko import util
import lassie
logger = logging.getLogger('main')
def d2u(d):
if isinstance(d, dict):
return dict((k.replace('-', '_'), v) for (k, v) in d.iteritems())
else:
return d.replace('-', '_')
subcommand_descriptions = {
'init': 'create initial configuration file',
'scan': 'detect seismic events',
'map-geometry': 'make station map'
}
subcommand_usages = {
'init': 'init',
'scan': 'scan [options]',
'map-geometry': 'map-geometry <configfile> [options]',
}
subcommands = subcommand_descriptions.keys()
program_name = 'lassie'
usage_tdata = d2u(subcommand_descriptions)
usage_tdata['program_name'] = program_name
usage = program_name + ''' <subcommand> [options] [--] <arguments> ...
Subcommands:
init %(init)s
scan %(scan)s
map-geometry %(map_geometry)s
To get further help and a list of available options for any subcommand run:
%(program_name)s <subcommand> --help
''' % usage_tdata
def add_common_options(parser):
parser.add_option(
'--loglevel',
action='store',
dest='loglevel',
type='choice',
choices=('critical', 'error', 'warning', 'info', 'debug'),
default='info',
help ='set logger level to '
'"critical", "error", "warning", "info", or "debug". '
'Default is "%default".')
def process_common_options(options):
util.setup_logging(program_name, options.loglevel)
def cl_parse(command, args, setup=None):
usage = subcommand_usages[command]
descr = subcommand_descriptions[command]
if isinstance(usage, basestring):
usage = [usage]
susage = '%s %s' % (program_name, usage[0])
for s in usage[1:]:
susage += '\n%s%s %s' % (' '*7, program_name, s)
parser = OptionParser(
usage=susage,
description=descr[0].upper() + descr[1:] + '.')
if setup:
setup(parser)
add_common_options(parser)
(options, args) = parser.parse_args(args)
process_common_options(options)
return parser, options, args
def die(message, err=''):
if err:
sys.exit('%s: error: %s \n %s' % (program_name, message, err))
else:
sys.exit('%s: error: %s' % (program_name, message))
def help_and_die(parser, message):
parser.print_help(sys.stderr)
sys.stderr.write('\n')
die(message)
def escape(s):
return s.replace("'", "\\'")
def command_init(args):
def setup(parser):
parser.add_option(
'--stations', dest='stations_path',
help='stations file')
parser, options, args = cl_parse('init', args, setup=setup)
if args:
data_path = '\n'.join("- '%s'" % escape(x) for x in args)
else:
data_path = "- 'DATA'"
if options.stations_path:
stations_path = options.stations_path
else:
stations_path = 'STATIONS_PATH'
print '''%%YAML 1.1
--- !lassie.Config
## Configuration file for Lassie, your friendly earthquake detector
##
## Receiver coordinates can be read from a stations file in Pyrocko format:
stations_path: '%(stations_path)s'
## Receivers can also be listed in the config file, lat/lon and carthesian
## (x/y/z) = (North/East/Down) coordinates are supported and may be combined
## (interpreted as reference + offset). Omitted values are treated as zero.
# receivers:
# - !lassie.Receiver
# codes: ['', 'ACC13', '']
# lat: 10.
# lon: 12.
# x: 2397.56
# y: 7331.94
# z: -404.1
## List of data directories. Lassie will recurse into subdirectories to find
## all contained waveform files.
data_paths:
%(data_path)s
## Processing time interval (default: use time interval of available data)
# tmin: '2016-02-05 16:00:00'
# tmax: '2016-02-05 17:40:00'
## Search grid, if not given here, a default grid will be chosen
# grid:
# !coarse.Carthesian3DGrid
# lat: 10.
# lon: 12.
# xmin: 2360.0
# xmax: 2450.0
# ymin: 7270.0
# ymax: 7360.0
# zmin: -415.0
# zmax: -400.0
# dx: 2.0
# dy: 2.0
# dz: 2.0
## Size factor to use when automatically choosing a grid size
autogrid_radius_factor: 1.5
## Grid density factor used when automatically choosing a grid
autogrid_density_factor: 10.0
## Image function contributions
ifcs:
- !lassie.WavePacketIFC
name: 'P'
weight: 1.0
fmin: 1.0
fmax: 50.0
fsmooth_factor: 0.1
shifter: !lassie.VelocityShifter
velocity: 3600.
- !lassie.WavePacketIFC
name: 'S'
weight: 1.0
fmin: 1.0
fmax: 10.0
fsmooth_factor: 0.1
shifter: !lassie.VelocityShifter
velocity: 2500.
## Whether to divide image function frames by their mean value
sharpness_normalization: true
## Threshold on detector function
detector_threshold: 200.0
## Output filename for detections
detections_path: 'detections.txt'
''' % dict(
stations_path=stations_path,
data_path=data_path)
def command_scan(args):
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser.add_option(
'--show-detections', dest='show_detections', action='store_true',
help='show plot for every detection found')
parser.add_option(
'--show-movie', dest='show_movie', action='store_true',
help='show movie when showing detections')
parser.add_option(
'--stop-after-first', dest='stop_after_first', action='store_true',
help='show plot for every detection found')
parser, options, args = cl_parse('scan', args, setup=setup)
if len(args) != 1:
help_and_die(parser, 'missing argument')
config_path = args[0]
config = lassie.read_config(config_path)
print config
try:
lassie.scan(
config,
force=options.force,
show_detections=options.show_detections,
show_movie=options.show_movie,
stop_after_first=options.stop_after_first)
except lassie.LassieError, e:
die(str(e))
def command_map_geometry(args):
parser, options, args = cl_parse('map-geometry', args)
if len(args) != 2:
help_and_die(parser, 'missing arguments')
config_path = args[0]
output_path = args[1]
config = lassie.read_config(config_path)
lassie.map_geometry(config, output_path)
if __name__ == '__main__':
usage_sub = 'fomosto %s [options]'
if len(sys.argv) < 2:
sys.exit('Usage: %s' % usage)
args = list(sys.argv)
args.pop(0)
command = args.pop(0)
if command in subcommands:
globals()['command_' + d2u(command)](args)
elif command in ('--help', '-h', 'help'):
if command == 'help' and args:
acommand = args[0]
if acommand in subcommands:
globals()['command_' + acommand](['--help'])
sys.exit('Usage: %s' % usage)
else:
die('no such subcommand: %s' % command)

148
scratch/lassie.svg

File diff suppressed because one or more lines are too long

After

Width:  |  Height:  |  Size: 227 KiB

BIN
scratch/lassie1.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 12 KiB

144
scratch/lassie1.svg

@ -0,0 +1,144 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:cc="http://creativecommons.org/ns#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns="http://www.w3.org/2000/svg"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
width="744.09448819"
height="1052.3622047"
id="svg3805"
version="1.1"
inkscape:version="0.48.3.1 r9886"
sodipodi:docname="lassie1.svg">
<defs
id="defs3807" />
<sodipodi:namedview
id="base"
pagecolor="#ffffff"
bordercolor="#666666"
borderopacity="1.0"
inkscape:pageopacity="0.0"
inkscape:pageshadow="2"
inkscape:zoom="1.4"
inkscape:cx="454.50577"
inkscape:cy="578.42402"
inkscape:document-units="px"
inkscape:current-layer="layer1"
showgrid="false"
inkscape:window-width="1074"
inkscape:window-height="1371"
inkscape:window-x="0"
inkscape:window-y="501"
inkscape:window-maximized="0" />
<metadata
id="metadata3810">
<rdf:RDF>
<cc:Work
rdf:about="">
<dc:format>image/svg+xml</dc:format>
<dc:type
rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
<dc:title></dc:title>
</cc:Work>
</rdf:RDF>
</metadata>
<g
inkscape:label="Layer 1"
inkscape:groupmode="layer"
id="layer1">
<path
transform="translate(70.535714,-71.964286)"
d="m 380.71429,498.61218 a 56.964287,56.964287 0 1 1 -113.92858,0 56.964287,56.964287 0 1 1 113.92858,0 z"
sodipodi:ry="56.964287"
sodipodi:rx="56.964287"
sodipodi:cy="498.61218"
sodipodi:cx="323.75"
id="path3844"
style="color:#000000;fill:#d5dded;stroke:none;stroke-width:2;marker:none;visibility:visible;display:inline;overflow:visible;enable-background:accumulate"
sodipodi:type="arc" />
<path
inkscape:connector-curvature="0"
style="fill:#f4ddcf;fill-opacity:1;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 362.74553,385.08929 c -1.38896,0.3788 -0.88007,7.35073 -0.375,19.09375 0.25783,5.99457 5.05059,13.11432 -0.40625,20.96875 -3.94019,5.67142 -2.26729,12.30631 -3.65625,16.03125 -1.67259,4.48558 -3.00126,16.94793 -3.8125,26.21875 10.26782,10.02588 24.29678,16.21875 39.78125,16.21875 11.61041,0 22.40089,-3.50422 31.40625,-9.46875 0.5671,-7.17552 1.47518,-19.50779 1.5625,-26.40625 0.12627,-9.97525 -5.78125,-18.4375 -5.78125,-18.4375 0,0 0.23731,0.125 1.5,0.125 1.26269,0 2.99357,-40.23926 -2.65625,-43.1875 -5.13725,-2.68077 -12.625,15.65625 -12.625,15.65625 -9.13153,-1.48998 -15.18322,-0.0731 -17.28125,-0.125 -2.09803,-0.0519 -2.03029,2.25603 -8.34375,-4.5625 -7.52055,-8.12219 -17.92354,-12.50381 -19.3125,-12.125 z"
id="path2996" />
<path
style="fill:#d89c66;stroke:none"
d="m 419.05257,477.02812 5.3033,-2.77792 c 0,0 2.2501,-19.83692 0.63134,-28.78935 -0.77887,-4.30748 -7.37219,-9.97898 -7.37219,-9.97898 0,0 0.66177,10.12979 -4.04061,13.92768 -2.53894,2.05059 0.36265,5.53936 0.30113,7.66806 -0.18641,6.45008 5.17703,19.95051 5.17703,19.95051 z"
id="path3779"
inkscape:connector-curvature="0"
sodipodi:nodetypes="ccscssc" />
<path
style="fill:#d89c66;fill-opacity:1;stroke:none"
d="m 399.25371,406.88816 c -3.91328,1.40839 3.99341,16.97482 0.89068,19.74437 -4.02382,3.59173 -3.33685,13.35914 -3.33685,13.35914 l 2.42078,5.72172 -8.83883,-4.04061 c 0,0 -10.07652,-10.03707 -14.89975,-7.32361 -4.12485,2.32056 -2.39778,9.54747 -1.26269,14.14214 1.15361,4.66964 7.91132,7.31252 7.82868,12.12183 -0.12156,7.07453 -12.12183,17.42513 -12.12183,17.42513 l -6.06092,-27.77919 7.57615,-18.94036 8.33376,-7.57615 -5.80838,-3.03046 -7.57614,-8.58629 -3.283,-16.92006 9.57851,-3.96664 15.6753,13.31055 18.68783,-2.0203 9.34391,-10.35406 5.05076,3.78807 2.77792,29.29442 -3.03046,4.54569 c 0,0 -5.98553,-2.65221 -10.18818,-4.46647 -2.47434,-1.06815 -9.22143,-19.36151 -11.75725,-18.44886 z"
id="path3008"
inkscape:connector-curvature="0"
sodipodi:nodetypes="sscccsssccccccccccccccss" />
<path
style="fill:#1f1a17;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 383.06588,421.09092 c 0,0 0.37881,-3.91434 4.41942,-5.17703 4.04061,-1.26269 4.41942,5.93465 4.41942,5.93465 0,0 -0.37881,1.38896 -4.04061,0.88388 -3.66181,-0.50508 -4.79823,-1.6415 -4.79823,-1.6415 z"
id="path2998"
inkscape:connector-curvature="0" />
<path
style="fill:#1f1a17;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 410.46627,417.17658 c 0,0 4.04061,0.25254 4.04061,1.89404 0,1.64149 -1.89403,3.53553 -2.90419,2.65165 -1.01015,-0.88389 -1.13642,-4.54569 -1.13642,-4.54569 z"
id="path3000"
inkscape:connector-curvature="0" />
<path
style="fill:#1f1a17;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 407.68835,432.07633 c 6.43972,-0.37881 8.20749,3.53553 7.44988,6.69226 -0.75762,3.15673 -1.76777,6.9448 -3.66181,7.44988 -1.89403,0.50507 -8.46002,-4.29315 -8.33375,-7.57615 0.12626,-3.28299 -0.12627,-5.55584 4.54568,-6.56599 z"
id="path3002"
inkscape:connector-curvature="0"
sodipodi:nodetypes="csssc" />
<path
style="fill:#1f1a17;stroke:none"
d="m 383.69722,436.18008 5.99779,4.54568 c 0,0 1.57836,8.01809 1.07329,8.7757 -0.50508,0.75762 3.03045,5.17703 3.28299,7.95495 0.25254,2.77792 3.66181,9.84899 7.95495,9.84899 4.29315,0 8.08122,-0.37881 8.08122,-0.37881 0,0 -0.75761,3.78808 -3.53553,4.04061 -2.77792,0.25254 -9.72272,-1.51523 -13.25825,-8.83883 -3.53554,-7.32361 -4.79823,-4.04061 -6.56599,-8.96511 -1.76777,-4.92449 -2.02031,-9.21764 -2.02031,-9.21764 z"
id="path3004"
inkscape:connector-curvature="0"
sodipodi:nodetypes="ccssscssscc" />
<path
style="fill:#f4b38d;fill-opacity:1;stroke:none"
d="m 397.06818,454.16273 3.0625,4.23224 4.41418,5.71669 6.96402,-0.77486 2.62942,3.33863 c 2.60054,3.30196 2.52951,4.29901 4.04473,3.28886 1.21367,-0.80911 1.57037,-4.30186 -2.91973,-7.04728 l -5.48366,-3.10532 c -0.60135,-2.34694 0.0621,-6.75868 1.32389,-7.38959 2.14286,-1.07143 -3.79553,-4.52031 -3.79553,-4.52031 l -7.125,1.4375 -1.72256,2.00147 -0.7124,-3.61752 -6.41772,-3.75002 4.39732,8.93303 z"
id="path3006"
inkscape:connector-curvature="0"
sodipodi:nodetypes="ccccssccsccccccc" />
<path
style="fill:#1f1a17;stroke:none"
d="m 362.09208,389.13303 c 0,0 8.54622,7.49128 13.00308,11.01231 3.39574,2.68272 10.04226,3.38871 10.45542,7.69654 0.52039,5.42583 -11.63751,6.59752 -11.0863,12.02031 0.23721,2.3337 5.5751,1.76555 5.71429,4.10714 0.22999,3.86908 -6.23915,5.04907 -8.44516,8.23595 -3.90677,5.64385 -5.45379,12.24469 -6.01912,19.08547 -0.72796,8.80861 4.55134,25.8512 4.30737,26.54587 -0.0962,0.27396 -12.81971,-7.26526 -13.05738,-10.11729 0,0 1.00828,-14.24657 2.32143,-21.25 1.10161,-5.87524 2.75572,-11.64951 4.64286,-17.32143 0.91077,-2.73737 2.72908,-5.1919 3.21428,-8.03571 0.47095,-2.76026 -0.0234,-5.61267 -0.35714,-8.39286 -0.9555,-7.95911 -4.69363,-23.5863 -4.69363,-23.5863 z"
id="path2995"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cssssssscssssc" />
<path
style="fill:#1f1a17;stroke:none"
d="m 372.22578,391.76519 -11.43753,-4.25507 2.6569,7.40437 z"
id="path3765"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cccc" />
<path
style="fill:#1f1a17;stroke:none"
d="m 416.56381,392.63639 7.10324,-3.51015 -2.34309,7.46038 z"
id="path3767"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cccc" />
<path
style="font-size:medium;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-indent:0;text-align:start;text-decoration:none;line-height:normal;letter-spacing:normal;word-spacing:normal;text-transform:none;direction:ltr;block-progression:tb;writing-mode:lr-tb;text-anchor:start;baseline-shift:baseline;color:#000000;fill:#1f1a17;fill-opacity:1;fill-rule:nonzero;stroke:none;stroke-width:1px;marker:none;visibility:visible;display:inline;overflow:visible;enable-background:accumulate;font-family:Sans;-inkscape-font-specification:Sans"
d="m 391.03566,444.79773 0.71875,0.77427 3.52237,7.54854 0.0312,0.0937 0.0625,0.0625 1.28125,1.1875 0.0625,0.0625 3,4.15625 0.0312,0 0,0.0312 4.40625,5.6875 0.15625,0.25 3.18544,-0.15747 3.78331,-0.62373 2.46875,3.125 c 1.28377,1.63002 1.88922,2.65618 2.46875,3.28125 0.28977,0.31253 0.63941,0.55734 1.0625,0.5625 0.42309,0.005 0.77861,-0.19615 1.1875,-0.46875 0.86431,-0.57621 1.2323,-1.84317 0.875,-3.3125 -0.3573,-1.46933 -1.48191,-3.13749 -3.8125,-4.5625 l -0.0312,0 -5.21875,-2.96875 c -0.24073,-1.0601 -0.26081,-2.60554 -0.0312,-3.96875 0.1193,-0.70844 0.2923,-1.35753 0.5,-1.84375 1.26488,-1.56914 1.05586,-3.15473 0.16196,-4.72256 -11.67466,-4.48039 -13.26586,-1.36973 -19.31816,-5.46494 l -1.05869,0.13452 z m 1.45987,0.44781 c 9.95602,5.05169 5.6192,5.88154 15.79204,3.81885 0,0 -1.48944,2.57667 -0.11977,2.98144 0.33145,0.0979 2.55236,-0.98393 2.79648,-0.55029 0.0947,0.17582 0.10063,0.28078 0.0937,0.3125 -0.007,0.0317 -0.009,0.067 -0.1875,0.15625 -0.50825,0.25413 -0.80971,0.75198 -1.0625,1.34375 -0.25279,0.59177 -1.69607,3.34733 -1.82519,4.11406 -1.06128,4.82716 3.35017,3.49229 7.01264,5.91714 l 0.0312,0 c 2.14205,1.31642 3.06137,2.77625 3.34375,3.9375 0.28346,1.16567 -0.11939,2.04835 -0.46875,2.28125 -0.34872,0.23248 -0.55398,0.28212 -0.625,0.28125 -0.071,-8.7e-4 -0.1368,-0.0268 -0.34375,-0.25 -0.41389,-0.44641 -1.08948,-1.54681 -2.40625,-3.21875 l -2.625,-3.34375 -0.1875,-0.21875 -2.14404,-1.10522 -4.79346,-0.13389 -3.49239,-5.24746 -3.82011,-2.45098 -0.0312,-0.0312 -0.0312,-0.0312 0.26523,-3.46034 -5.17148,-5.10216 z"
id="path3803"
inkscape:connector-curvature="0"
sodipodi:nodetypes="cccccccccccccccccsscccccccccccsccsscccscsccccccccccccc" />
<path
sodipodi:type="arc"
style="color:#000000;fill:none;stroke:#1f1a17;stroke-width:2;marker:none;visibility:visible;display:inline;overflow:visible;enable-background:accumulate"
id="path3780"
sodipodi:cx="323.75"
sodipodi:cy="498.61218"
sodipodi:rx="56.964287"
sodipodi:ry="56.964287"
d="m 380.71429,498.61218 a 56.964287,56.964287 0 1 1 -113.92858,0 56.964287,56.964287 0 1 1 113.92858,0 z"
transform="translate(70.535714,-71.964286)" />
</g>
</svg>

After

Width:  |  Height:  |  Size: 10 KiB

10
setup.py

@ -0,0 +1,10 @@
#!/usr/bin/env python
from distutils.core import setup
setup(
name='lassie',
packages=['lassie'],
package_dir={'lassie': 'src'},
scripts=['apps/lassie'],
package_data={'lassie': []})

8
src/__init__.py

@ -0,0 +1,8 @@
from lassie.common import * # noqa
from lassie.receiver import * # noqa
from lassie.grid import * # noqa
from lassie.shifter import * # noqa
from lassie.ifc import * # noqa
from lassie.core import * # noqa
from lassie.plot import * # noqa
from lassie.config import * # noqa

19
src/common.py

@ -0,0 +1,19 @@
from collections import defaultdict
guts_prefix = 'lassie'
class LassieError(Exception):
pass
def grouped_by(l, key):
d = defaultdict(list)
for x in l:
d[key(x)].append(x)
return d
__all__ = [
'LassieError',
]

138
src/config.py

@ -0,0 +1,138 @@
import logging
from pyrocko.guts import Object, String, Float, Timestamp, List, Bool, load
from pyrocko import model
from lassie import receiver, ifc, grid, geo
guts_prefix = 'lassie'
logger = logging.getLogger('lassie.config')
class Config(Object):
stations_path = String.T(
optional=True,
help='stations file in Pyrocko format')
receivers = List.T(
receiver.Receiver.T(),
help='receiver coordinates if not read from file')
data_paths = List.T(
String.T(),
help='list of directories paths to search for data')
blacklist = List.T(
String.T(),
help='codes in the form NET.STA.LOC of receivers to be excluded')
tmin = Timestamp.T(
optional=True,
help='beginning of time interval to be processed')
tmax = Timestamp.T(
optional=True,
help='end of time interval to be processed')
detections_path = String.T(
optional=True,
help='output file name for for detections')
grid = grid.Grid.T(
optional=True,
help='definition of search grid, if not given, a default grid is '
'chosen')
autogrid_radius_factor = Float.T(
default=1.5,
help='size factor to use when automatically choosing a grid size')
autogrid_density_factor = Float.T(
default=10.0,
help='grid density factor used when automatically choosing a grid '
'spacing')
ifcs = List.T(
ifc.IFC.T(),
help='image function contributions')
sharpness_normalization = Bool.T(
default=True,
help='whether to divide image function frames by their mean value')
detector_threshold = Float.T(
default=200.,
help='threshold on detector function')
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._receivers = None
self._grid = None
def get_receivers(self):
'''Aggregate receivers from different sources.'''
if self._receivers is None:
self._receivers = list(self.receivers)
if self.stations_path:
for station in model.load_stations(self.stations_path):
self._receivers.append(
receiver.Receiver(
codes=station.nsl(),
lat=station.lat,
lon=station.lon))
return self._receivers
def get_grid(self):
'''Get grid or make default grid.'''
if self._grid is None:
if not self.grid:
receivers = self.get_receivers()
fsmooth_max = max(
ifc.get_fsmooth() for ifc in self.ifcs)
vmin = min(ifc.shifter.get_vmin() for ifc in self.ifcs)
spacing = vmin / fsmooth_max / self.autogrid_density_factor
lat0, lon0, north, east, depth = geo.bounding_box_square(
*receiver.receivers_coords(receivers),
scale=self.autogrid_radius_factor)
self._grid = grid.Carthesian3DGrid(
lat=lat0,
lon=lon0,
xmin=north[0],
xmax=north[1],
dx=spacing,
ymin=east[0],
ymax=east[1],
dy=spacing,
zmin=depth[0],
zmax=depth[1],
dz=spacing)
logger.info('automatic grid:\n%s' % self._grid)
else:
self._grid = self.grid
self._grid.update()
return self._grid
def read_config(file_path):
conf = load(filename=file_path)
return conf
__all__ = [
'Config',
'read_config',
]

220
src/core.py

@ -0,0 +1,220 @@
import sys
import os
import logging
import os.path as op
import numpy as num
from pyrocko import pile, trace, util, io
from pyrocko.parstack import parstack
from lassie import common, plot, grid as gridmod
logger = logging.getLogger('lassie.core')
def scan(
config,
show_detections=False,
show_movie=False,
force=False,
stop_after_first=False):
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)
util.ensuredirs(config.detections_path)
grid = config.get_grid()
receivers = config.get_receivers()
norm_map = gridmod.geometrical_normalization(grid, receivers)
shift_tables = []
tshift_maxs = []
for ifc in config.ifcs:
shift_tables.append(ifc.shifter.get_table(grid, receivers))
tshift_maxs.append(shift_tables[-1].max())
tshift_max = max(tshift_maxs) * 1.0
tinc = tshift_max * 10
tpad = max(ifc.get_tpad() for ifc in config.ifcs) + tshift_max
fsmooth_min = min(ifc.get_fsmooth() for ifc in config.ifcs)
blacklist = set(tuple(s.split('.')) for s in config.blacklist)
station_index = dict(
(rec.codes, i) for (i, rec) in enumerate(receivers)
if rec.codes not in blacklist)
p = pile.make_pile(config.data_paths, fileformat='detect')
ngridpoints = grid.size()
idetection = 0
station_weights = {}
tmin = config.tmin
tmax = 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):
if not trs:
continue
logger.info('processing time window %s - %s' % (
util.time_to_str(trs[0].wmin),
util.time_to_str(trs[0].wmax)))
wmin = trs[0].wmin
wmax = trs[0].wmax
frames = None
pdata = []
trs_debug = []
shift_maxs = []
for iifc, ifc in enumerate(config.ifcs):
dataset = ifc.preprocess(trs, wmin, wmax, tshift_max)
if not dataset:
continue
nstations_selected = len(dataset)
nsls_selected, trs_selected = zip(*dataset)
trs_debug.extend(trs + list(trs_selected))
deltat_cf = trs_selected[0].deltat
t0 = (wmin / deltat_cf) * deltat_cf
istations_selected = num.array(
[station_index[nsl] for nsl in nsls_selected], dtype=num.int)
arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
offsets = num.array(
[int(round((tr.tmin-t0) / deltat_cf)) for tr in trs_selected],
dtype=num.int32)
w = num.array(
[station_weights.get(nsl, 1.0) for nsl in nsls_selected],
dtype=num.float)
weights = num.ones((ngridpoints, nstations_selected))
weights *= w[num.newaxis, :]
weights *= ifc.weight
shift_table = shift_tables[iifc]
shifts = -num.round(
shift_table[:, istations_selected] /
deltat_cf).astype(num.int32)
pdata.append((list(trs_selected), shift_table, ifc))
iwmin = int(round((wmin-t0) / deltat_cf))
iwmax = int(round((wmax-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,
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]
frame_maxs = frames.max(axis=0)
tmin_frames = t0 + ioff * deltat_cf
tr_stackmax = trace.Trace(
'', 'SMAX', '', '',
tmin=tmin_frames,
deltat=deltat_cf,
ydata=frame_maxs)
trs_debug.append(tr_stackmax)
# trace.snuffle(trs_debug)
ydata_window = tr_stackmax.chop(wmin, wmax, inplace=False).get_ydata()
logger.info('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)
for (tpeak, apeak) in zip(tpeaks, apeaks):
if not (wmin <= tpeak and tpeak < wmax):
continue
logger.info('detection: %s %g' % (
util.time_to_str(tpeak),
apeak))
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)
idetection += 1
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.close()
if show_detections:
fmin = min(ifc.fmin for ifc in config.ifcs)
fmax = min(ifc.fmax for ifc in config.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)
if stop_after_first:
return
tr_stackmax.chop(wmin, wmax)
io.save([tr_stackmax], 'stackmax/trace_%(tmin)s.mseed')
__all__ = [
'scan',
]

113
src/geo.py

@ -0,0 +1,113 @@
import numpy as num
from pyrocko import orthodrome as od
def float_array_broadcast(*args):
return num.broadcast_arrays(*[
num.asarray(x, dtype=num.float) for x in args])
def surface_distance(alat, alon, anorth, aeast, blat, blon, bnorth, beast):
args = float_array_broadcast(
alat, alon, anorth, aeast, blat, blon, bnorth, beast)
want_scalar = False
if len(args[0].shape) == 0:
want_scalar = True
args = [num.atleast_1d(x) for x in args]
(alat, alon, anorth, aeast, blat, blon, bnorth, beast) = args
eqref = num.logical_and(alat == blat, alon == blon)
neqref = num.logical_not(eqref)
dist = num.empty_like(alat)
dist[eqref] = num.sqrt(
(bnorth[eqref] - anorth[eqref])**2 + (beast[eqref] - aeast[eqref])**2)
aalat, aalon = od.ne_to_latlon(
alat[neqref], alon[neqref], anorth[neqref], aeast[neqref])
bblat, bblon = od.ne_to_latlon(
blat[neqref], blon[neqref], bnorth[neqref], beast[neqref])
dist[neqref] = od.distance_accurate50m_numpy(aalat, aalon, bblat, bblon)
if want_scalar:
return dist[0]
else:
return dist
def bounding_box(lat, lon, north, east, depth, scale=1.0):
lat, lon, north, east, depth = float_array_broadcast(
lat, lon, north, east, depth)
if num.all(lat[0] == lat) and num.all(lon[0] == lon):
return _scaled_bb(
lat[0], lon[0],
(num.min(north), num.max(north)),
(num.min(east), num.max(east)),
(num.min(depth), num.max(depth)), scale)
else:
elat, elon = od.ne_to_latlon(lat, lon, north, east)
enorth, eeast = od.latlon_to_ne_numpy(elat[0], elon[0], elat, elon)
enorth_min = num.min(enorth)
enorth_max = num.max(enorth)
eeast_min = num.min(eeast)
eeast_max = num.max(eeast)
mnorth = 0.5*(enorth_min + enorth_max)
meast = 0.5*(eeast_min + eeast_max)
lat0, lon0 = od.ne_to_latlon(elat[0], elon[0], mnorth, meast)
return _scaled_bb(
lat0, lon0,
(enorth_min - mnorth, enorth_max - mnorth),
(eeast_min - meast, eeast_max - meast),
(num.min(depth), num.max(depth)),
scale)
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)
dnorth = north_max - north_min
mnorth = 0.5 * (north_min + north_max)
deast = east_max - east_min
meast = 0.5 * (east_min + east_max)
dmax = max(dnorth, deast)
if dnorth < dmax:
north_min = mnorth - 0.5*dmax
north_max = mnorth + 0.5*dmax
if deast < dmax:
east_min = meast - 0.5*dmax
east_max = meast + 0.5*dmax
return _scaled_bb(
lat, lon,
(north_min, north_max),
(east_min, east_max),
(depth_min, depth_max),
scale)
def _scaled_range(ra, scale):
mi, ma = ra
d = ma - mi
m = 0.5 * (ma + mi)
return m - 0.5*scale*d, m + 0.5*scale*d
def _scaled_bb(lat, lon, north, east, depth, scale):
return (
lat, lon,
_scaled_range(north, scale),
_scaled_range(east, scale),
_scaled_range(depth, scale))

168
src/grid.py

@ -0,0 +1,168 @@
import math
import numpy as num
from pyrocko.guts import Object, Float
from pyrocko import orthodrome as od
from lassie import receiver
guts_prefix = 'lassie'
class Grid(Object):
pass
class Carthesian3DGrid(Grid):
lat = Float.T(default=0.0)
lon = Float.T(default=0.0)
xmin = Float.T()
xmax = Float.T()
ymin = Float.T()
ymax = Float.T()
zmin = Float.T()
zmax = Float.T()
dx = Float.T()
dy = Float.T()
dz = Float.T()
def __init__(self, **kwargs):
for k in kwargs.keys():
kwargs[k] = float(kwargs[k])
Grid.__init__(self, **kwargs)
self.update()
def size(self):
nx, ny, nz = self._shape()
return nx * ny * nz
def max_delta(self):
return max(self.dx, self.dy, self.dz)
def update(self):
self._coords = None
def _shape(self):
nx = int(round(((self.xmax - self.xmin) / self.dx))) + 1
ny = int(round(((self.ymax - self.ymin) / self.dy))) + 1
nz = int(round(((self.zmax - self.zmin) / self.dz))) + 1
return nx, ny, nz
def _get_coords(self):
if self._coords is None:
nx, ny, nz = self._shape()
x = num.linspace(self.xmin, self.xmax, nx)
y = num.linspace(self.ymin, self.ymax, ny)
z = num.linspace(self.zmin, self.zmax, nz)
self._coords = x, y, z
return self._coords
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 distances(self, receivers):
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
rx, ry = receiver.receivers_coords(
receivers, system=('ne', self.lat, self.lon))
rz = num.array([r.z for r in receivers], dtype=num.float)
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 +
(z[:, na, na, na] - rz[na, na, na, :])**2).reshape((nx*ny*nz, nr))
return distances
def surface_points(self, system='latlon'):
x, y, z = self._get_coords()
xs = num.tile(x, y.size)
ys = num.repeat(y, x.size)
if system == 'latlon':
return od.ne_to_latlon(self.lat, self.lon, xs, ys)
elif system[0] == 'ne':
lat0, lon0 = system[1:]
if lat0 == self.lat and lon0 == self.lon:
return xs, ys
else:
elats, elons = od.ne_to_latlon(self.lat, self.lon, xs, ys)
return od.latlon_to_ne_numpy(lat0, lon0, elats, elons)
def plot_points(self, axes, system='latlon'):
x, y = self.surface_points(system=system)
axes.plot(y, x, '.', color='black', ms=1.0)
def plot(
self, axes, a,
amin=None,
amax=None,
z_slice=None,
cmap=None,
system='latlon'):
if system == 'latlon':
assert False, 'not implemented yet'
if not (system[1] == self.lat and system[2] == self.lon):
assert False, 'not implemented yet'
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
a3d = a.reshape((nz, ny, nx))
if z_slice is not None:
iz = num.argmin(num.abs(z-z_slice))
a2d = a3d[iz, :, :]
else:
a2d = num.max(a3d, axis=0)
axes.pcolormesh(y, x, a2d.T, vmin=amin, vmax=amax, cmap=cmap)
def geometrical_normalization(grid, receivers):
distances = grid.distances(receivers)
delta_grid = grid.max_delta()
delta_ring = delta_grid * 3.0
ngridpoints, nstations = distances.shape
norm_map = num.zeros(ngridpoints)
for istation in xrange(nstations):
dists_station = distances[:, istation]
dist_min = num.floor(num.min(dists_station) / delta_ring) * delta_ring
dist_max = num.ceil(num.max(dists_station) / delta_ring) * delta_ring
dist = dist_min
while dist < dist_max:
indices = num.where(num.logical_and(
dist <= dists_station,
dists_station < dist + delta_ring))[0]
nexpect = math.pi * ((dist + delta_ring)**2 - dist**2) /\
delta_grid**2
norm_map[indices] += indices.size / nexpect
dist += delta_ring
norm_map /= nstations
return norm_map
__all__ = [
'Grid',
'Carthesian3DGrid',
]

182
src/ifc.py

@ -0,0 +1,182 @@
import numpy as num
from pyrocko.guts import Object, String, Float
from pyrocko import trace, autopick
from lassie import shifter, common
guts_prefix = 'lassie'
class IFC(Object):
'''Image function contribution.'''
name = String.T()
weight = Float.T(default=1.0)
fmin = Float.T()
fmax = Float.T()
shifter = shifter.Shifter.T()
def get_tpad(self):
return 4. / self.fmin
def preprocess(self, trs, wmin, wmax, tpad_new):
pass
class WavePacketIFC(IFC):
'''
An image function useful as a simple fast general purpose event detector.
This detector type image function is sensible to transient wave packets
traveling with a given velocity through the network, e.g. P phase together
with P coda or S together with S coda or surface waves packets.
'''
fsmooth_factor = Float.T(default=0.1)
def get_tpad(self):
return 4. / self.get_fsmooth()
def get_fsmooth(self):
return self.fmin * self.fsmooth_factor
def preprocess(self, trs, wmin, wmax, tpad_new):
fsmooth = self.fmin * self.fsmooth_factor
if not trs:
return []
deltat_cf = trs[0].deltat
while deltat_cf * 2 < 0.25 / self.fmin:
deltat_cf *= 2
trs_filt = []
for orig_tr in trs:
tr = orig_tr.copy()
tr.highpass(4, self.fmin, demean=True)
tr.lowpass(4, self.fmax, demean=False)
tr.ydata = tr.ydata**2
n = int(num.round(1./fsmooth / tr.deltat))
taper = num.hanning(n)
tr.set_ydata(num.convolve(tr.get_ydata(), taper))
tr.set_ydata(num.maximum(tr.ydata, 0.0))
tr.shift(-(n/2.*tr.deltat))
tr.downsample_to(deltat_cf, snap=True, demean=False)
trs_filt.append(tr)
trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])
dataset = []
for nsl in sorted(trs_by_nsl.keys()):
sumtr = None
for tr in trs_by_nsl[nsl]:
if sumtr is None:
sumtr = tr.copy()
sumtr.set_codes(channel='CF')
else:
sumtr.add(tr)
navg = 5./fsmooth / deltat_cf
yavg = trace.moving_avg(sumtr.ydata, navg)
if num.any(yavg == 0.0):
continue
sumtr.ydata /= yavg
# sumtr.ydata -= 1.0
sumtr.chop(wmin - tpad_new, wmax + tpad_new)
print sumtr.max()
dataset.append((nsl, sumtr))
return dataset
class OnsetIFC(IFC):
'''
An image function sensible to sharp onsets in the signal.
This image function is based on an STA/LTA.
'''
short_window_factor = Float.T(default=1.0)
window_ratio = Float.T(default=5.0)
fsmooth_factor = Float.T(default=0.1)
def get_fsmooth(self):
return self.fmin * self.fsmooth_factor
def preprocess(self, trs, wmin, wmax, tpad_new):
fsmooth = self.fmin * self.fsmooth_factor
if not trs:
return []
deltat_cf = trs[0].deltat
while deltat_cf * 2 < 0.25 / self.fmin:
deltat_cf *= 2
trs_filt = []
for orig_tr in trs:
tr = orig_tr.copy()
tr.highpass(4, self.fmin, demean=True)
tr.lowpass(4, self.fmax, demean=False)
tr.ydata = tr.ydata**2
trs_filt.append(tr)
trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])
swin = self.short_window_factor / self.fmin
lwin = self.window_ratio * swin
duration = 1.0/self.fmin
n = int(num.round(duration / trs_filt[0].deltat))
dataset = []
for nsl in sorted(trs_by_nsl.keys()):
sumtr = None
for tr in trs_by_nsl[nsl]:
if sumtr is None:
sumtr = tr.copy()
sumtr.set_codes(channel='CF')
else:
sumtr.add(tr)
sumtr.set_ydata(sumtr.get_ydata().astype(num.float32))
autopick.recursive_stalta(swin, lwin, 1.0, 4.0, 3.0, sumtr)
sumtr.shift(-(swin + n/2.*sumtr.deltat))
ntap = int(num.round(1./fsmooth / tr.deltat))
taper = num.hanning(ntap)
sumtr.shift(-(ntap/2.*tr.deltat))
sumtr.set_ydata(
num.convolve(sumtr.get_ydata()/len(trs_by_nsl), taper))
# sumtr.lowpass(4, fsmooth*4, demean=False)
# sumtr.set_ydata(num.maximum(tr.ydata, 0.0))
# sumtr.downsample_to(deltat_cf, snap=True, demean=False)
sumtr.downsample_to(deltat_cf, snap=True, demean=False)
sumtr.chop(wmin - tpad_new, wmax + tpad_new)
print sumtr.max()
dataset.append((nsl, sumtr))
return dataset
__all__ = [
'IFC',
'WavePacketIFC',
'OnsetIFC',
]

305
src/plot.py

@ -0,0 +1,305 @@
import time
import numpy as num
from pyrocko import automap
from lassie import receiver, grid as gridmod, geo
def map_gmt(
receivers, grid, center_lat, center_lon, radius, output_path,
width=25., height=25.,
show_station_labels=False):
station_lats = num.array([r.lat for r in receivers])
station_lons = num.array([r.lon for r in receivers])
map = automap.Map(
width=width,
height=height,
lat=center_lat,
lon=center_lon,
radius=radius,
show_rivers=False,
show_topo=False,
illuminate_factor_land=0.35,
color_dry=(240, 240, 235),
topo_cpt_wet='white_sea_land',
topo_cpt_dry='white_sea_land')
map.gmt.psxy(
in_columns=(station_lons, station_lats),
S='t8p',
G='black',
*map.jxyr)
if grid:
map.gmt.psxy(
in_columns=(grid.lons, grid.lats),
S='c1p',
G='black',
*map.jxyr)
if show_station_labels:
for r in receivers:
map.add_label(r.lat, r.lon, '%s' % r.station)
map.save(output_path)
def map_geometry(config, output_path):
receivers = config.get_receivers()
grid = config.get_grid()
lat0, lon0, north, east, depth = geo.bounding_box_square(
*receiver.receivers_coords(receivers),
scale=config.autogrid_radius_factor)
radius = max((north[1] - north[0]), (east[1] - east[0]))
radius *= config.autogrid_radius_factor * 1.5
map_gmt(receivers, grid, lat0, lon0, radius, output_path,
show_station_labels=False)
def plot_receivers(axes, receivers, system='latlon'):
x, y = receiver.receivers_coords(receivers, system=system)
axes.plot(y, x, '^', color='black', ms=10.)
def plot_geometry_carthesian(grid, receivers):
from matplotlib import pyplot as plt
from pyrocko.cake_plot import mpl_init, labelspace
mpl_init()
plt.figure(figsize=(9, 9))
axes = plt.subplot2grid((1, 1), (0, 0), aspect=1.0)
labelspace(axes)
grid.plot_points(axes, system=('ne', grid.lat, grid.lon))
plot_receivers(axes, receivers, system=('ne', grid.lat, grid.lon))
distances = grid.distances(receivers)
delta_grid = max(grid.dx, grid.dy, grid.dz)
norm_map = gridmod.geometrical_normalization(distances, delta_grid)
grid.plot(axes, norm_map, system=('ne', grid.lat, grid.lon))
plt.show()
def plot_detection(
grid, receivers, frames, tmin_frames, deltat_cf, imax, iframe,
fsmooth_min, xpeak, ypeak, zpeak, tr_stackmax, tpeaks, apeaks,
detector_threshold, pdata, trs_raw, fmin, fmax, idetection,
grid_station_shift_max,
movie=False):
from matplotlib import pyplot as plt
from matplotlib import cm
from pyrocko.cake_plot import mpl_init, labelspace, colors, \
str_to_mpl_color as scolor
distances = grid.distances(receivers)
mpl_init()
fig = plt.figure(figsize=(16, 9))
axes = plt.subplot2grid((2, 3), (0, 2), aspect=1.0)
labelspace(axes)
axes2 = plt.subplot2grid((2, 3), (1, 2))
labelspace(axes2)
axes3 = plt.subplot2grid((2, 3), (0, 1), rowspan=2)
axes4 = plt.subplot2grid((2, 3), (0, 0), rowspan=2)
axes.set_xlabel('X [m]')
axes.set_ylabel('Y [m]')
axes.locator_params(nbins=6, tight=True)
axes2.set_xlabel('Time [s]')
axes2.set_ylabel('Detector level')
axes3.set_xlabel('Time [s]')
for el in axes3.get_yticklabels():
el.set_visible(False)
axes4.set_xlabel('Time [s]')
axes4.set_ylabel('Distance [m]')
tpeak_current = tmin_frames + deltat_cf * iframe
t0 = tpeak_current
tduration = 2.0*grid_station_shift_max + 1./fsmooth_min
for tpeak, apeak in zip(tpeaks, apeaks):
axes2.axvline(
tpeak-t0,
color=scolor('aluminium3'),
lw=2.)
axes2.axvspan(
tpeak_current-0.5*tduration - t0,
tpeak_current+0.5*tduration - t0,
color=scolor('aluminium2'),
lw=2.)
axes2.axhline(
detector_threshold,
color=scolor('aluminium6'),
lw=2.)
t = tr_stackmax.get_xdata()
amp = tr_stackmax.get_ydata()
axes2.plot(t - t0, amp, color=scolor('scarletred2'), lw=2.)
for tpeak, apeak in zip(tpeaks, apeaks):
axes2.plot(
tpeak-t0, apeak, '*',
ms=20.,
mfc='white',
mec='black')
station_index = dict(
(rec.codes, i) for (i, rec) in enumerate(receivers))
dists = []
amps = []
shifts = []
pdata2 = []
for trs, shift_table, shifter in pdata:
trs = [tr.copy() for tr in trs]
for tr in trs:
istation = station_index[tr.nslc_id[:3]]
shift = shift_table[imax, istation]
tr2 = tr.chop(
tpeak_current - 0.5*tduration + shift,
tpeak_current + 0.5*tduration + shift,
inplace=False)
dists.append(distances[imax, istation])
amp = tr2.get_ydata() * shifter.weight
amps.append(num.max(num.abs(amp)))
shifts.append(shift)
pdata2.append((trs, shift_table, shifter))
dist_min = min(dists)
dist_max = max(dists)
shift_min = min(shifts)
shift_max = max(shifts)
amp_max = max(amps)
scalefactor = (dist_max - dist_min) / len(trs) * 0.5
axes3.set_xlim(-0.5*tduration + shift_min, 0.5*tduration + shift_max)
axes3.set_ylim(dist_min - scalefactor, dist_max + scalefactor)
axes4.set_xlim(-0.5*tduration + shift_min, 0.5*tduration + shift_max)
axes4.set_ylim(dist_min - scalefactor, dist_max + scalefactor)
axes3.axvline(
0.,
color=scolor('aluminium3'),
lw=2.)
for ishifter, (trs, shift_table, shifter) in enumerate(pdata2):
color = colors[ishifter % len(colors)]
for tr, dist in zip(trs, dists):
tr = tr.chop(
tpeak_current - 0.5*tduration + shift_min,
tpeak_current + 0.5*tduration + shift_max, inplace=False)
istation = station_index[tr.nslc_id[:3]]
shift = shift_table[imax, istation]
print ishifter, shift
axes3.plot(shift, dist, '|', mew=2, mec=color, ms=10)
t = tr.get_xdata()
amp = tr.get_ydata() * shifter.weight
amp /= amp_max
axes3.plot(
t-t0,
dist + scalefactor*amp + ishifter*scalefactor*0.1,
color=color)
for tr in trs_raw:
istation = station_index[tr.nslc_id[:3]]
dist = distances[imax, istation]
shift = shift_table[imax, istation]
tr = tr.copy()
tr.highpass(4, fmin, demean=True)
tr.lowpass(4, fmax, demean=False)
tr.chop(
tpeak_current - 0.5*tduration + shift_min,
tpeak_current + 0.5*tduration + shift_max)
t = tr.get_xdata()
amp = tr.get_ydata().astype(num.float)
amp = amp / num.max(num.abs(amp))
axes4.plot(t-t0, dist + scalefactor*amp, color='black')
axes4.plot(shift, dist, '|', mew=2, mec=color, ms=10)
nframes = frames.shape[1]
iframe_min = max(0, int(round(iframe - 0.5*tduration/deltat_cf)))
iframe_max = min(nframes-1, int(round(iframe + 0.5*tduration/deltat_cf)))
amax = frames[imax, iframe]
cmap = cm.YlOrBr
system = ('ne', grid.lat, grid.lon)
if movie:
plt.ion()
plt.show()
for iframe in xrange(iframe_min, iframe_max+1):
frame = frames[:, iframe]
axes.cla()
grid.plot(
axes, frame,
amin=0.0,
amax=amax,
cmap=cmap,
system=system)
plot_receivers(axes, receivers, system=system)
axes.plot(
ypeak, xpeak, '*',
ms=20.,
mfc='white',
mec='black')
plt.draw()
time.sleep(0.05)
plt.ioff()
else:
frame = num.max(frames[:, iframe_min:iframe_max+1], axis=1)
grid.plot(axes, frame, amin=0.0, amax=amax, cmap=cmap, system=system)
plot_receivers(axes, receivers, system=system)
axes.plot(
ypeak, xpeak, '*',
ms=20.,
mec='black',
mfc='white')
fig.tight_layout()
plt.show()
__all__ = [
'map_geometry',
'map_gmt',
'plot_detection',
'plot_geometry_carthesian',
'plot_receivers',
]

42
src/receiver.py

@ -0,0 +1,42 @@
import numpy as num
from pyrocko.guts import Object, Tuple, Float, String
import pyrocko.orthodrome as od
guts_prefix = 'lassie'
class Receiver(Object):
codes = Tuple.T(3, String.T())