1#!/usr/bin/env python
2# http://pyrocko.org - GPLv3
3#
4# The Pyrocko Developers, 21st Century
5# ---|P------/S----------~Lg----------
6from __future__ import print_function
8import sys
9import re
10import os.path as op
11import logging
12import copy
13import shutil
14from optparse import OptionParser
16from pyrocko import util, trace, gf, cake, io, config, fomosto
17from pyrocko.gui import marker
18from pyrocko.util import mpl_show
20logger = logging.getLogger('pyrocko.apps.fomosto')
21km = 1e3
24def d2u(d):
25 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
28subcommand_descriptions = {
29 'init': 'create a new empty GF store',
30 'build': 'compute GFs and fill into store',
31 'stats': 'print information about a GF store',
32 'check': 'check for problems in GF store',
33 'decimate': 'build decimated variant of a GF store',
34 'redeploy': 'copy traces from one GF store into another',
35 'view': 'view selected traces',
36 'extract': 'extract selected traces',
37 'import': 'convert Kiwi GFDB to GF store format',
38 'export': 'convert GF store to Kiwi GFDB format',
39 'ttt': 'create travel time tables',
40 'tttview': 'plot travel time table',
41 'tttextract': 'extract selected travel times',
42 'tttlsd': 'fix holes in travel time tables',
43 'server': 'run seismosizer server',
44 'download': 'download GF store from a server',
45 'modelview': 'plot earthmodels',
46 'upgrade': 'upgrade store format to latest version',
47 'addref': 'import citation references to GF store config',
48 'qc': 'quality check',
49 'report': 'report for Green\'s Function databases',
50}
52subcommand_usages = {
53 'init': ['init <type> <store-dir> [options]',
54 'init redeploy <source> <destination> [options]'],
55 'build': 'build [store-dir] [options]',
56 'stats': 'stats [store-dir] [options]',
57 'check': 'check [store-dir] [options]',
58 'decimate': 'decimate [store-dir] <factor> [options]',
59 'redeploy': 'redeploy <source> <destination> [options]',
60 'view': 'view [store-dir] ... [options]',
61 'extract': 'extract [store-dir] <selection>',
62 'import': 'import <source> <destination> [options]',
63 'export': 'export [store-dir] <destination> [options]',
64 'ttt': 'ttt [store-dir] [options]',
65 'tttview': 'tttview [store-dir] <phase-ids> [options]',
66 'tttextract': 'tttextract [store-dir] <phase> <selection>',
67 'tttlsd': 'tttlsd [store-dir] <phase>',
68 'server': 'server [options] <store-super-dir> ...',
69 'download': 'download [options] <site> <store-id>',
70 'modelview': 'modelview <selection>',
71 'upgrade': 'upgrade [store-dir] ...',
72 'addref': 'addref [store-dir] ... <filename> ...',
73 'qc': 'qc [store-dir]',
74 'report': 'report <subcommand> <arguments>... [options]'
75}
77subcommands = subcommand_descriptions.keys()
79program_name = 'fomosto'
81usage = program_name + ''' <subcommand> <arguments> ... [options]
83Subcommands:
85 init %(init)s
86 build %(build)s
87 stats %(stats)s
88 check %(check)s
89 decimate %(decimate)s
90 redeploy %(redeploy)s
91 view %(view)s
92 extract %(extract)s
93 import %(import)s
94 export %(export)s
95 ttt %(ttt)s
96 tttview %(tttview)s
97 tttextract %(tttextract)s
98 tttlsd %(tttlsd)s
99 server %(server)s
100 download %(download)s
101 modelview %(modelview)s
102 upgrade %(upgrade)s
103 addref %(addref)s
104 qc %(qc)s
105 report %(report)s
107To get further help and a list of available options for any subcommand run:
109 fomosto <subcommand> --help
111''' % d2u(subcommand_descriptions)
114def add_common_options(parser):
115 parser.add_option(
116 '--loglevel',
117 action='store',
118 dest='loglevel',
119 type='choice',
120 choices=('critical', 'error', 'warning', 'info', 'debug'),
121 default='info',
122 help='set logger level to '
123 '"critical", "error", "warning", "info", or "debug". '
124 'Default is "%default".')
127def process_common_options(options):
128 util.setup_logging(program_name, options.loglevel)
131def cl_parse(command, args, setup=None, details=None):
132 usage = subcommand_usages[command]
133 descr = subcommand_descriptions[command]
135 if isinstance(usage, str):
136 usage = [usage]
138 susage = '%s %s' % (program_name, usage[0])
139 for s in usage[1:]:
140 susage += '\n%s%s %s' % (' '*7, program_name, s)
142 description = descr[0].upper() + descr[1:] + '.'
144 if details:
145 description = description + ' %s' % details
147 parser = OptionParser(usage=susage, description=description)
148 parser.format_description = lambda formatter: description
150 if setup:
151 setup(parser)
153 add_common_options(parser)
154 (options, args) = parser.parse_args(args)
155 process_common_options(options)
156 return parser, options, args
159def die(message, err=''):
160 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
163def fomo_wrapper_module(name):
164 try:
165 if not re.match(gf.meta.StringID.pattern, name):
166 raise ValueError('invalid name')
168 words = name.split('.', 1)
169 if len(words) == 2:
170 name, variant = words
171 else:
172 name = words[0]
173 variant = None
175 name_clean = re.sub(r'[.-]', '_', name)
176 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
177 mod = __import__(modname, level=0)
178 return getattr(mod.fomosto, name_clean), variant
180 except ValueError:
181 die('invalid modelling code wrapper name')
183 except ImportError:
184 die('''modelling code wrapper "%s" not available or not installed
185 (module probed: "%s")''' % (name, modname))
188def command_init(args):
190 details = '''
192Available modelling backends:
193%s
195 More information at
196 https://pyrocko.org/docs/current/apps/fomosto/backends.html
197''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
199 parser, options, args = cl_parse(
200 'init', args,
201 details=details)
203 if len(args) == 0:
204 sys.exit(parser.format_help())
206 if args[0] == 'redeploy':
207 if len(args) != 3:
208 parser.error('incorrect number of arguments')
210 source_dir, dest_dir = args[1:]
212 try:
213 source = gf.Store(source_dir)
214 except gf.StoreError as e:
215 die(e)
217 config = copy.deepcopy(source.config)
218 config.derived_from_id = source.config.id
219 try:
220 config_filenames = gf.store.Store.create_editables(
221 dest_dir, config=config)
223 except gf.StoreError as e:
224 die(e)
226 try:
227 dest = gf.Store(dest_dir)
228 except gf.StoreError as e:
229 die(e)
231 for k in source.extra_keys():
232 source_fn = source.get_extra_path(k)
233 dest_fn = dest.get_extra_path(k)
234 shutil.copyfile(source_fn, dest_fn)
236 logger.info(
237 '(1) configure settings in files:\n %s'
238 % '\n '.join(config_filenames))
240 logger.info(
241 '(2) run "fomosto redeploy <source> <dest>", as needed')
243 else:
244 if len(args) != 2:
245 parser.error('incorrect number of arguments')
247 (modelling_code_id, store_dir) = args
249 module, variant = fomo_wrapper_module(modelling_code_id)
250 try:
251 config_filenames = module.init(store_dir, variant)
252 except gf.StoreError as e:
253 die(e)
255 logger.info('(1) configure settings in files:\n %s'
256 % '\n '.join(config_filenames))
257 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
258 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
261def get_store_dir(args):
262 if len(args) == 1:
263 store_dir = op.abspath(args.pop(0))
264 else:
265 store_dir = op.abspath(op.curdir)
267 if not op.isdir(store_dir):
268 die('not a directory: %s' % store_dir)
270 return store_dir
273def get_store_dirs(args):
274 if len(args) == 0:
275 store_dirs = [op.abspath(op.curdir)]
276 else:
277 store_dirs = [op.abspath(x) for x in args]
279 for store_dir in store_dirs:
280 if not op.isdir(store_dir):
281 die('not a directory: %s' % store_dir)
283 return store_dirs
286def command_build(args):
288 def setup(parser):
289 parser.add_option(
290 '--force', dest='force', action='store_true',
291 help='overwrite existing files')
293 parser.add_option(
294 '--nworkers', dest='nworkers', type='int', metavar='N',
295 help='run N worker processes in parallel')
297 parser.add_option(
298 '--continue', dest='continue_', action='store_true',
299 help='continue suspended build')
301 parser.add_option(
302 '--step', dest='step', type='int', metavar='I',
303 help='process block number IBLOCK')
305 parser.add_option(
306 '--block', dest='iblock', type='int', metavar='I',
307 help='process block number IBLOCK')
309 parser, options, args = cl_parse('build', args, setup=setup)
311 store_dir = get_store_dir(args)
312 try:
313 if options.step is not None:
314 step = options.step - 1
315 else:
316 step = None
318 if options.iblock is not None:
319 iblock = options.iblock - 1
320 else:
321 iblock = None
323 store = gf.Store(store_dir)
324 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
325 module.build(
326 store_dir,
327 force=options.force,
328 nworkers=options.nworkers, continue_=options.continue_,
329 step=step,
330 iblock=iblock)
332 except gf.StoreError as e:
333 die(e)
336def command_stats(args):
338 parser, options, args = cl_parse('stats', args)
339 store_dir = get_store_dir(args)
341 try:
342 store = gf.Store(store_dir)
343 s = store.stats()
345 except gf.StoreError as e:
346 die(e)
348 for k in store.stats_keys:
349 print('%s: %s' % (k, s[k]))
352def command_check(args):
354 parser, options, args = cl_parse('check', args)
355 store_dir = get_store_dir(args)
357 try:
358 store = gf.Store(store_dir)
359 problems = store.check(show_progress=True)
360 if problems:
361 die('problems detected with gf store: %s' % store_dir)
363 except gf.StoreError as e:
364 die(e)
367def load_config(fn):
368 try:
369 config = gf.meta.load(filename=fn)
370 assert isinstance(config, gf.Config)
372 except Exception:
373 die('cannot load gf config from file: %s' % fn)
375 return config
378def command_decimate(args):
380 def setup(parser):
381 parser.add_option(
382 '--config', dest='config_fn', metavar='FILE',
383 help='use modified spacial sampling given in FILE')
385 parser.add_option(
386 '--force', dest='force', action='store_true',
387 help='overwrite existing files')
389 parser, options, args = cl_parse('decimate', args, setup=setup)
390 try:
391 decimate = int(args.pop())
392 except Exception:
393 parser.error('cannot get <factor> argument')
395 store_dir = get_store_dir(args)
397 config = None
398 if options.config_fn:
399 config = load_config(options.config_fn)
401 try:
402 store = gf.Store(store_dir)
403 store.make_decimated(decimate, config=config, force=options.force,
404 show_progress=True)
406 except gf.StoreError as e:
407 die(e)
410def sindex(args):
411 return '(%s)' % ', '.join('%g' % x for x in args)
414def command_redeploy(args):
416 parser, options, args = cl_parse('redeploy', args)
418 if not len(args) == 2:
419 sys.exit(parser.format_help())
421 source_store_dir, dest_store_dir = args
423 try:
424 source = gf.Store(source_store_dir)
425 except gf.StoreError as e:
426 die(e)
428 try:
429 gf.store.Store.create_dependants(dest_store_dir)
430 except gf.StoreError:
431 pass
433 try:
434 dest = gf.Store(dest_store_dir, 'w')
436 except gf.StoreError as e:
437 die(e)
439 show_progress = True
441 if show_progress:
442 pbar = util.progressbar('redeploying', dest.config.nrecords)
444 for i, args in enumerate(dest.config.iter_nodes()):
445 try:
446 tr = source.get(args, interpolation='off')
447 dest.put(args, tr)
449 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e:
450 logger.debug('skipping %s, (%s)' % (sindex(args), e))
452 except gf.store.StoreError as e:
453 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
455 if show_progress:
456 pbar.update(i+1)
458 if show_progress:
459 pbar.finish()
462def command_view(args):
463 def setup(parser):
464 parser.add_option('--extract',
465 dest='extract',
466 metavar='start:stop[:step|@num],...',
467 help='specify which traces to show')
469 parser.add_option('--show-phases',
470 dest='showphases',
471 default=None,
472 metavar='[phase_id1,phase_id2,...|all]',
473 help='add phase markers from ttt')
475 parser.add_option('--qt5',
476 dest='gui_toolkit_qt5',
477 action='store_true',
478 default=False,
479 help='use Qt5 for the GUI')
481 parser.add_option('--qt4',
482 dest='gui_toolkit_qt4',
483 action='store_true',
484 default=False,
485 help='use Qt4 for the GUI')
487 parser.add_option('--opengl',
488 dest='opengl',
489 action='store_true',
490 default=False,
491 help='use OpenGL for drawing')
493 parser, options, args = cl_parse('view', args, setup=setup)
495 if options.gui_toolkit_qt4:
496 config.override_gui_toolkit = 'qt4'
498 if options.gui_toolkit_qt5:
499 config.override_gui_toolkit = 'qt5'
501 gdef = None
502 if options.extract:
503 try:
504 gdef = gf.meta.parse_grid_spec(options.extract)
505 except gf.meta.GridSpecError as e:
506 die(e)
508 store_dirs = get_store_dirs(args)
510 alpha = 'abcdefghijklmnopqrstxyz'.upper()
512 markers = []
513 traces = []
515 try:
516 for istore, store_dir in enumerate(store_dirs):
517 store = gf.Store(store_dir)
518 if options.showphases == 'all':
519 phasenames = [pn.id for pn in store.config.tabulated_phases]
520 elif options.showphases is not None:
521 phasenames = options.showphases.split(',')
523 ii = 0
524 for args in store.config.iter_extraction(gdef):
525 gtr = store.get(args)
527 loc_code = ''
528 if len(store_dirs) > 1:
529 loc_code = alpha[istore % len(alpha)]
531 if gtr:
533 sta_code = '%04i (%s)' % (
534 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
535 tmin = gtr.deltat * gtr.itmin
537 tr = trace.Trace(
538 '',
539 sta_code,
540 loc_code,
541 '%02i' % args[-1],
542 ydata=gtr.data,
543 deltat=gtr.deltat,
544 tmin=tmin)
546 if options.showphases:
547 for phasename in phasenames:
548 phase_tmin = store.t(phasename, args[:-1])
549 if phase_tmin:
550 m = marker.PhaseMarker(
551 [('',
552 sta_code,
553 loc_code,
554 '%02i' % args[-1])],
555 phase_tmin,
556 phase_tmin,
557 0,
558 phasename=phasename)
559 markers.append(m)
561 traces.append(tr)
563 ii += 1
565 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
566 die(e)
568 trace.snuffle(traces, markers=markers, opengl=options.opengl)
571def command_extract(args):
572 def setup(parser):
573 parser.add_option(
574 '--format', dest='format', default='mseed',
575 choices=['mseed', 'sac', 'text', 'yaff'],
576 help='export to format "mseed", "sac", "text", or "yaff". '
577 'Default is "mseed".')
579 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
580 parser.add_option(
581 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
582 help='output path template [default: "%s"]' % fndfl)
584 parser, options, args = cl_parse('extract', args, setup=setup)
585 try:
586 sdef = args.pop()
587 except Exception:
588 parser.error('cannot get <selection> argument')
590 try:
591 gdef = gf.meta.parse_grid_spec(sdef)
592 except gf.meta.GridSpecError as e:
593 die(e)
595 store_dir = get_store_dir(args)
597 extensions = {
598 'mseed': 'mseed',
599 'sac': 'sac',
600 'text': 'txt',
601 'yaff': 'yaff'}
603 try:
604 store = gf.Store(store_dir)
605 for args in store.config.iter_extraction(gdef):
606 gtr = store.get(args)
607 if gtr:
608 tr = trace.Trace(
609 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
610 ydata=gtr.data,
611 deltat=gtr.deltat,
612 tmin=gtr.deltat * gtr.itmin)
614 additional = dict(
615 args='_'.join('%g' % x for x in args),
616 irecord=store.str_irecord(args),
617 extension=extensions[options.format])
619 io.save(
620 tr,
621 options.output_fn,
622 format=options.format,
623 additional=additional)
625 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
626 die(e)
629def command_import(args):
630 try:
631 from tunguska import gfdb
632 except ImportError:
633 die('the kiwi tools must be installed to use this feature')
635 parser, options, args = cl_parse('import', args)
637 show_progress = True
639 if not len(args) == 2:
640 sys.exit(parser.format_help())
642 source_path, dest_store_dir = args
644 if op.isdir(source_path):
645 source_path = op.join(source_path, 'db')
647 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
649 db = gfdb.Gfdb(source_path)
651 config = gf.meta.ConfigTypeA(
652 id='imported_gfs',
653 distance_min=db.firstx,
654 distance_max=db.firstx + (db.nx-1) * db.dx,
655 distance_delta=db.dx,
656 source_depth_min=db.firstz,
657 source_depth_max=db.firstz + (db.nz-1) * db.dz,
658 source_depth_delta=db.dz,
659 sample_rate=1.0/db.dt,
660 ncomponents=db.ng
661 )
663 try:
664 gf.store.Store.create(dest_store_dir, config=config)
665 dest = gf.Store(dest_store_dir, 'w')
666 if show_progress:
667 pbar = util.progressbar(
668 'importing', dest.config.nrecords/dest.config.ncomponents)
670 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
671 source_depth, distance = [float(x) for x in args]
672 traces = db.get_traces_pyrocko(distance, source_depth)
673 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
674 for ig in range(db.ng):
675 if ig in ig_to_trace:
676 tr = ig_to_trace[ig]
677 gf_tr = gf.store.GFTrace(
678 tr.get_ydata(),
679 int(round(tr.tmin / tr.deltat)),
680 tr.deltat)
682 else:
683 gf_tr = gf.store.Zero
685 dest.put((source_depth, distance, ig), gf_tr)
687 if show_progress:
688 pbar.update(i+1)
690 if show_progress:
691 pbar.finish()
693 dest.close()
695 except gf.StoreError as e:
696 die(e)
699def command_export(args):
700 from subprocess import Popen, PIPE
702 try:
703 from tunguska import gfdb
704 except ImportError as err:
705 die('the kiwi tools must be installed to use this feature', err)
707 def setup(parser):
708 parser.add_option(
709 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
710 help='split output gfdb into N chunks')
712 parser, options, args = cl_parse('export', args, setup=setup)
714 show_progress = True
716 if len(args) not in (1, 2):
717 sys.exit(parser.format_help())
719 target_path = args.pop()
720 if op.isdir(target_path):
721 target_path = op.join(target_path, 'kiwi_gfdb')
722 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
724 source_store_dir = get_store_dir(args)
726 source = gf.Store(source_store_dir, 'r')
727 config = source.config
729 if not isinstance(config, gf.meta.ConfigTypeA):
730 die('only stores of type A can be exported to Kiwi format')
732 if op.isfile(target_path + '.index'):
733 die('destation already exists')
735 cmd = [str(x) for x in [
736 'gfdb_build',
737 target_path,
738 options.nchunks,
739 config.ndistances,
740 config.nsource_depths,
741 config.ncomponents,
742 config.deltat,
743 config.distance_delta,
744 config.source_depth_delta,
745 config.distance_min,
746 config.source_depth_min]]
748 p = Popen(cmd, stdin=PIPE)
749 p.communicate()
751 out_db = gfdb.Gfdb(target_path)
753 if show_progress:
754 pbar = util.progressbar(
755 'exporting', config.nrecords/config.ncomponents)
757 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
759 data_out = []
760 for ig in range(config.ncomponents):
761 try:
762 tr = source.get((z, x, ig), interpolation='off')
763 data_out.append((tr.t, tr.data * config.factor))
765 except gf.store.StoreError as e:
766 logger.warning('cannot get %s, (%s)' % (sindex((z, x, ig)), e))
767 data_out.append(None)
769 # put a zero valued sample to no-data zero-traces at a compatible time
770 tmins = [
771 entry[0][0]
772 for entry in data_out
773 if entry is not None and entry[0].size != 0]
775 if tmins:
776 tmin = min(tmins)
777 for entry in data_out:
778 if entry is not None and entry[0].size == 0:
779 entry[0].resize(1)
780 entry[1].resize(1)
781 entry[0][0] = tmin
782 entry[1][0] = 0.0
784 out_db.put_traces_slow(x, z, data_out)
786 if show_progress:
787 pbar.update(i+1)
789 if show_progress:
790 pbar.finish()
792 source.close()
795def phasedef_or_horvel(x):
796 try:
797 return float(x)
798 except ValueError:
799 return cake.PhaseDef(x)
802def mkp(s):
803 return [phasedef_or_horvel(ps) for ps in s.split(',')]
806def command_ttt(args):
807 def setup(parser):
808 parser.add_option(
809 '--force', dest='force', action='store_true',
810 help='overwrite existing files')
812 parser, options, args = cl_parse('ttt', args, setup=setup)
814 store_dir = get_store_dir(args)
815 try:
816 store = gf.Store(store_dir)
817 store.make_ttt(force=options.force)
819 except gf.StoreError as e:
820 die(e)
823def command_tttview(args):
824 import numpy as num
825 import matplotlib.pyplot as plt
826 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
827 mpl_init()
829 def setup(parser):
830 parser.add_option(
831 '--source-depth', dest='source_depth', type=float,
832 help='Source depth in km')
834 parser.add_option(
835 '--receiver-depth', dest='receiver_depth', type=float,
836 help='Receiver depth in km')
838 parser, options, args = cl_parse(
839 'tttview', args, setup=setup,
840 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
842 try:
843 phase_ids = args.pop().split(',')
844 except Exception:
845 parser.error('cannot get <phase-ids> argument')
847 np = 1
848 store_dir = get_store_dir(args)
849 for phase_id in phase_ids:
850 try:
851 store = gf.Store(store_dir)
853 if options.receiver_depth is not None:
854 receiver_depth = options.receiver_depth * 1000.0
855 else:
856 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
857 receiver_depth = store.config.receiver_depth
859 elif isinstance(store.config, gf.ConfigTypeB):
860 receiver_depth = store.config.receiver_depth_min
862 else:
863 receiver_depth = 0.0
865 phase = store.get_stored_phase(phase_id)
866 axes = plt.subplot(2, len(phase_ids), np)
867 labelspace(axes)
868 xscaled(1./km, axes)
869 yscaled(1./km, axes)
870 x = None
871 if isinstance(store.config, gf.ConfigTypeB):
872 x = (receiver_depth, None, None)
874 phase.plot_2d(axes, x=x)
875 axes.set_title(phase_id)
876 np += 1
877 except gf.StoreError as e:
878 die(e)
880 axes = plt.subplot(2, 1, 2)
881 num_d = 100
882 distances = num.linspace(store.config.distance_min,
883 store.config.distance_max,
884 num_d)
886 if options.source_depth is not None:
887 source_depth = options.source_depth * 1000.0
888 else:
889 source_depth = store.config.source_depth_min + (
890 store.config.source_depth_max - store.config.source_depth_min)/2.
892 if isinstance(store.config, gf.ConfigTypeA):
893 arrivals = num.empty(num_d)
894 for phase_id in phase_ids:
895 arrivals[:] = num.NAN
896 for i, d in enumerate(distances):
897 arrivals[i] = store.t(phase_id, (source_depth, d))
898 axes.plot(distances/1000.0, arrivals, label=phase_id)
899 axes.set_title('source source_depth %s km' % (source_depth/1000.0))
900 axes.set_xlabel('distance [km]')
901 axes.set_ylabel('TT [s]')
902 axes.legend()
904 plt.tight_layout()
905 mpl_show(plt)
908def command_tttextract(args):
909 def setup(parser):
910 parser.add_option(
911 '--output', dest='output_fn', metavar='TEMPLATE',
912 help='output to text files instead of stdout '
913 '(example TEMPLATE: "extracted/%(args)s.txt")')
915 parser, options, args = cl_parse('tttextract', args, setup=setup)
916 try:
917 sdef = args.pop()
918 except Exception:
919 parser.error('cannot get <selection> argument')
921 try:
922 sphase = args.pop()
923 except Exception:
924 parser.error('cannot get <phase> argument')
926 try:
927 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
928 except gf.meta.InvalidTimingSpecification:
929 parser.error('invalid phase specification: "%s"' % sphase)
931 try:
932 gdef = gf.meta.parse_grid_spec(sdef)
933 except gf.meta.GridSpecError as e:
934 die(e)
936 store_dir = get_store_dir(args)
938 try:
939 store = gf.Store(store_dir)
940 for args in store.config.iter_extraction(gdef, level=-1):
941 s = ['%e' % x for x in args]
942 for phase in phases:
943 t = store.t(phase, args)
944 if t is not None:
945 s.append('%e' % t)
946 else:
947 s.append('nan')
949 if options.output_fn:
950 d = dict(
951 args='_'.join('%e' % x for x in args),
952 extension='txt')
954 fn = options.output_fn % d
955 util.ensuredirs(fn)
956 with open(fn, 'a') as f:
957 f.write(' '.join(s))
958 f.write('\n')
959 else:
960 print(' '.join(s))
962 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
963 die(e)
966def command_tttlsd(args):
968 def setup(parser):
969 pass
971 parser, options, args = cl_parse('tttlsd', args, setup=setup)
973 try:
974 sphase_ids = args.pop()
975 except Exception:
976 parser.error('cannot get <phase> argument')
978 try:
979 phase_ids = [x.strip() for x in sphase_ids.split(',')]
980 except gf.meta.InvalidTimingSpecification:
981 parser.error('invalid phase specification: "%s"' % sphase_ids)
983 store_dir = get_store_dir(args)
985 try:
986 store = gf.Store(store_dir)
987 for phase_id in phase_ids:
988 store.fix_ttt_holes(phase_id)
990 except gf.StoreError as e:
991 die(e)
994def command_server(args):
995 from pyrocko.gf import server
997 def setup(parser):
998 parser.add_option(
999 '--port', dest='port', metavar='PORT', type='int', default=8080,
1000 help='serve on port PORT')
1002 parser.add_option(
1003 '--ip', dest='ip', metavar='IP', default='',
1004 help='serve on ip address IP')
1006 parser, options, args = cl_parse('server', args, setup=setup)
1008 engine = gf.LocalEngine(store_superdirs=args)
1009 server.run(options.ip, options.port, engine)
1012def command_download(args):
1013 from pyrocko.gf import ws
1015 details = '''
1017Browse pre-calculated Green's function stores online:
1019 https://greens-mill.pyrocko.org
1020'''
1022 def setup(parser):
1023 parser.add_option(
1024 '--force', dest='force', action='store_true',
1025 help='overwrite existing files')
1027 parser, options, args = cl_parse(
1028 'download', args, setup=setup, details=details)
1029 if not len(args) in (1, 2):
1030 sys.exit(parser.format_help())
1032 if len(args) == 2:
1033 site, store_id = args
1034 if not re.match(gf.meta.StringID.pattern, store_id):
1035 die('invalid store ID')
1036 else:
1037 site, store_id = args[0], None
1039 if site not in gf.ws.g_site_abbr:
1040 if -1 == site.find('://'):
1041 site = 'http://' + site
1043 try:
1044 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1045 except ws.DownloadError as e:
1046 die(str(e))
1049def command_modelview(args):
1051 import matplotlib.pyplot as plt
1052 import numpy as num
1053 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1054 mpl_init()
1056 neat_labels = {
1057 'vp': '$v_p$',
1058 'vs': '$v_s$',
1059 'qp': '$Q_p$',
1060 'qs': '$Q_s$',
1061 'rho': '$\\rho$'}
1063 def setup(parser):
1064 parser.add_option(
1065 '--parameters', dest='parameters',
1066 default='vp,vs', metavar='vp,vs,....',
1067 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1069 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1071 parser, options, args = cl_parse('modelview', args, setup=setup)
1073 store_dirs = get_store_dirs(args)
1075 parameters = options.parameters.split(',')
1077 fig, axes = plt.subplots(1,
1078 len(parameters),
1079 sharey=True,
1080 figsize=(len(parameters)*3, 5))
1082 if not isinstance(axes, num.ndarray):
1083 axes = [axes]
1085 axes = dict(zip(parameters, axes))
1087 for store_id in store_dirs:
1088 try:
1089 store = gf.Store(store_id)
1090 mod = store.config.earthmodel_1d
1091 z = mod.profile('z')
1093 for p in parameters:
1094 ax = axes[p]
1095 labelspace(ax)
1096 if '/' in p:
1097 p1, p2 = p.split('/')
1098 profile = mod.profile(p1)/mod.profile(p2)
1099 else:
1100 profile = mod.profile(p)
1102 ax.plot(profile, -z, label=store_id.split('/')[-1])
1103 if p in ['vs', 'vp', 'rho']:
1104 xscaled(1./1000., ax)
1106 yscaled(1./km, ax)
1108 except gf.StoreError as e:
1109 die(e)
1111 for p, ax in axes.items():
1112 ax.grid()
1113 if p in neat_labels:
1114 lab = neat_labels[p]
1115 elif all(x in neat_labels for x in p.split('/')):
1116 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1117 else:
1118 lab = p
1120 ax.set_title(lab, y=1.02)
1122 if p in units:
1123 lab += ' ' + units[p]
1125 ax.autoscale()
1126 ax.set_xlabel(lab)
1128 axes[parameters[0]].set_ylabel('Depth [km]')
1130 handles, labels = ax.get_legend_handles_labels()
1132 if len(store_dirs) > 1:
1133 ax.legend(
1134 handles,
1135 labels,
1136 bbox_to_anchor=(0.5, 0.12),
1137 bbox_transform=fig.transFigure,
1138 ncol=3,
1139 loc='upper center',
1140 fancybox=True)
1142 plt.subplots_adjust(bottom=0.22,
1143 wspace=0.05)
1145 mpl_show(plt)
1148def command_upgrade(args):
1149 parser, options, args = cl_parse('upgrade', args)
1150 store_dirs = get_store_dirs(args)
1151 try:
1152 for store_dir in store_dirs:
1153 store = gf.Store(store_dir)
1154 nup = store.upgrade()
1155 if nup == 0:
1156 print('%s: already up-to-date.' % store_dir)
1157 else:
1158 print('%s: %i file%s upgraded.' % (
1159 store_dir, nup, ['s', ''][nup == 1]))
1161 except gf.StoreError as e:
1162 die(e)
1165def command_addref(args):
1166 parser, options, args = cl_parse('addref', args)
1168 store_dirs = []
1169 filenames = []
1170 for arg in args:
1171 if op.isdir(arg):
1172 store_dirs.append(arg)
1173 elif op.isfile(arg):
1174 filenames.append(arg)
1175 else:
1176 die('invalid argument: %s' % arg)
1178 if not store_dirs:
1179 store_dirs.append('.')
1181 references = []
1182 try:
1183 for filename in args:
1184 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1185 except ImportError:
1186 die('pybtex module must be installed to use this function')
1188 if not references:
1189 die('no references found')
1191 for store_dir in store_dirs:
1192 try:
1193 store = gf.Store(store_dir)
1194 ids = [ref.id for ref in store.config.references]
1195 for ref in references:
1196 if ref.id in ids:
1197 die('duplicate reference id: %s' % ref.id)
1199 ids.append(ref.id)
1200 store.config.references.append(ref)
1202 store.save_config(make_backup=True)
1204 except gf.StoreError as e:
1205 die(e)
1208def command_qc(args):
1209 parser, options, args = cl_parse('qc', args)
1211 store_dir = get_store_dir(args)
1213 try:
1214 store = gf.Store(store_dir)
1215 s = store.stats()
1216 if s['empty'] != 0:
1217 print('has empty records')
1219 for aname in ['author', 'author_email', 'description', 'references']:
1221 if not getattr(store.config, aname):
1222 print('%s empty' % aname)
1224 except gf.StoreError as e:
1225 die(e)
1228def command_report(args):
1229 from pyrocko.fomosto.report import report_call
1230 report_call.run_program(args)
1233def main(args=None):
1234 if args is None:
1235 args = sys.argv[1:]
1237 if len(args) < 1:
1238 sys.exit('Usage: %s' % usage)
1240 command = args.pop(0)
1242 if command in subcommands:
1243 globals()['command_' + command](args)
1245 elif command in ('--help', '-h', 'help'):
1246 if command == 'help' and args:
1247 acommand = args[0]
1248 if acommand in subcommands:
1249 globals()['command_' + acommand](['--help'])
1251 sys.exit('Usage: %s' % usage)
1253 else:
1254 sys.exit('fomosto: error: no such subcommand: %s' % command)
1257if __name__ == '__main__':
1258 main(sys.argv[1:])