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 'sat': 'create stored ray attribute table',
44 'satview': 'plot stored ray attribute table',
45 'server': 'run seismosizer server',
46 'download': 'download GF store from a server',
47 'modelview': 'plot earthmodels',
48 'upgrade': 'upgrade store format to latest version',
49 'addref': 'import citation references to GF store config',
50 'qc': 'quality check',
51 'report': 'report for Green\'s Function databases',
52}
54subcommand_usages = {
55 'init': ['init <type> <store-dir> [options]',
56 'init redeploy <source> <destination> [options]'],
57 'build': 'build [store-dir] [options]',
58 'stats': 'stats [store-dir] [options]',
59 'check': 'check [store-dir] [options]',
60 'decimate': 'decimate [store-dir] <factor> [options]',
61 'redeploy': 'redeploy <source> <destination> [options]',
62 'view': 'view [store-dir] ... [options]',
63 'extract': 'extract [store-dir] <selection>',
64 'import': 'import <source> <destination> [options]',
65 'export': 'export [store-dir] <destination> [options]',
66 'ttt': 'ttt [store-dir] [options]',
67 'tttview': 'tttview [store-dir] <phase-ids> [options]',
68 'tttextract': 'tttextract [store-dir] <phase> <selection>',
69 'tttlsd': 'tttlsd [store-dir] <phase>',
70 'sat': 'sat [store-dir] [options]',
71 'satview': 'satview [store-dir] <phase-ids> [options]',
72 'server': 'server [options] <store-super-dir> ...',
73 'download': 'download [options] <site> <store-id>',
74 'modelview': 'modelview <selection>',
75 'upgrade': 'upgrade [store-dir] ...',
76 'addref': 'addref [store-dir] ... <filename> ...',
77 'qc': 'qc [store-dir]',
78 'report': 'report <subcommand> <arguments>... [options]'
79}
81subcommands = subcommand_descriptions.keys()
83program_name = 'fomosto'
85usage = program_name + ''' <subcommand> <arguments> ... [options]
87Subcommands:
89 init %(init)s
90 build %(build)s
91 stats %(stats)s
92 check %(check)s
93 decimate %(decimate)s
94 redeploy %(redeploy)s
95 view %(view)s
96 extract %(extract)s
97 import %(import)s
98 export %(export)s
99 ttt %(ttt)s
100 tttview %(tttview)s
101 tttextract %(tttextract)s
102 tttlsd %(tttlsd)s
103 sat %(sat)s
104 satview %(satview)s
105 server %(server)s
106 download %(download)s
107 modelview %(modelview)s
108 upgrade %(upgrade)s
109 addref %(addref)s
110 qc %(qc)s
111 report %(report)s
113To get further help and a list of available options for any subcommand run:
115 fomosto <subcommand> --help
117''' % d2u(subcommand_descriptions)
120def add_common_options(parser):
121 parser.add_option(
122 '--loglevel',
123 action='store',
124 dest='loglevel',
125 type='choice',
126 choices=('critical', 'error', 'warning', 'info', 'debug'),
127 default='info',
128 help='set logger level to '
129 '"critical", "error", "warning", "info", or "debug". '
130 'Default is "%default".')
133def process_common_options(options):
134 util.setup_logging(program_name, options.loglevel)
137def cl_parse(command, args, setup=None, details=None):
138 usage = subcommand_usages[command]
139 descr = subcommand_descriptions[command]
141 if isinstance(usage, str):
142 usage = [usage]
144 susage = '%s %s' % (program_name, usage[0])
145 for s in usage[1:]:
146 susage += '\n%s%s %s' % (' '*7, program_name, s)
148 description = descr[0].upper() + descr[1:] + '.'
150 if details:
151 description = description + ' %s' % details
153 parser = OptionParser(usage=susage, description=description)
154 parser.format_description = lambda formatter: description
156 if setup:
157 setup(parser)
159 add_common_options(parser)
160 (options, args) = parser.parse_args(args)
161 process_common_options(options)
162 return parser, options, args
165def die(message, err=''):
166 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
169def fomo_wrapper_module(name):
170 try:
171 if not re.match(gf.meta.StringID.pattern, name):
172 raise ValueError('invalid name')
174 words = name.split('.', 1)
175 if len(words) == 2:
176 name, variant = words
177 else:
178 name = words[0]
179 variant = None
181 name_clean = re.sub(r'[.-]', '_', name)
182 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
183 mod = __import__(modname, level=0)
184 return getattr(mod.fomosto, name_clean), variant
186 except ValueError:
187 die('invalid modelling code wrapper name')
189 except ImportError:
190 die('''modelling code wrapper "%s" not available or not installed
191 (module probed: "%s")''' % (name, modname))
194def command_init(args):
196 details = '''
198Available modelling backends:
199%s
201 More information at
202 https://pyrocko.org/docs/current/apps/fomosto/backends.html
203''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
205 parser, options, args = cl_parse(
206 'init', args,
207 details=details)
209 if len(args) == 0:
210 sys.exit(parser.format_help())
212 if args[0] == 'redeploy':
213 if len(args) != 3:
214 parser.error('incorrect number of arguments')
216 source_dir, dest_dir = args[1:]
218 try:
219 source = gf.Store(source_dir)
220 except gf.StoreError as e:
221 die(e)
223 config = copy.deepcopy(source.config)
224 config.derived_from_id = source.config.id
225 try:
226 config_filenames = gf.store.Store.create_editables(
227 dest_dir, config=config)
229 except gf.StoreError as e:
230 die(e)
232 try:
233 dest = gf.Store(dest_dir)
234 except gf.StoreError as e:
235 die(e)
237 for k in source.extra_keys():
238 source_fn = source.get_extra_path(k)
239 dest_fn = dest.get_extra_path(k)
240 shutil.copyfile(source_fn, dest_fn)
242 logger.info(
243 '(1) configure settings in files:\n %s'
244 % '\n '.join(config_filenames))
246 logger.info(
247 '(2) run "fomosto redeploy <source> <dest>", as needed')
249 else:
250 if len(args) != 2:
251 parser.error('incorrect number of arguments')
253 (modelling_code_id, store_dir) = args
255 module, variant = fomo_wrapper_module(modelling_code_id)
256 try:
257 config_filenames = module.init(store_dir, variant)
258 except gf.StoreError as e:
259 die(e)
261 logger.info('(1) configure settings in files:\n %s'
262 % '\n '.join(config_filenames))
263 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
264 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
267def get_store_dir(args):
268 if len(args) == 1:
269 store_dir = op.abspath(args.pop(0))
270 else:
271 store_dir = op.abspath(op.curdir)
273 if not op.isdir(store_dir):
274 die('not a directory: %s' % store_dir)
276 return store_dir
279def get_store_dirs(args):
280 if len(args) == 0:
281 store_dirs = [op.abspath(op.curdir)]
282 else:
283 store_dirs = [op.abspath(x) for x in args]
285 for store_dir in store_dirs:
286 if not op.isdir(store_dir):
287 die('not a directory: %s' % store_dir)
289 return store_dirs
292def command_build(args):
294 def setup(parser):
295 parser.add_option(
296 '--force', dest='force', action='store_true',
297 help='overwrite existing files')
299 parser.add_option(
300 '--nworkers', dest='nworkers', type='int', metavar='N',
301 help='run N worker processes in parallel')
303 parser.add_option(
304 '--continue', dest='continue_', action='store_true',
305 help='continue suspended build')
307 parser.add_option(
308 '--step', dest='step', type='int', metavar='I',
309 help='process block number IBLOCK')
311 parser.add_option(
312 '--block', dest='iblock', type='int', metavar='I',
313 help='process block number IBLOCK')
315 parser, options, args = cl_parse('build', args, setup=setup)
317 store_dir = get_store_dir(args)
318 try:
319 if options.step is not None:
320 step = options.step - 1
321 else:
322 step = None
324 if options.iblock is not None:
325 iblock = options.iblock - 1
326 else:
327 iblock = None
329 store = gf.Store(store_dir)
330 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
331 module.build(
332 store_dir,
333 force=options.force,
334 nworkers=options.nworkers, continue_=options.continue_,
335 step=step,
336 iblock=iblock)
338 except gf.StoreError as e:
339 die(e)
342def command_stats(args):
344 parser, options, args = cl_parse('stats', args)
345 store_dir = get_store_dir(args)
347 try:
348 store = gf.Store(store_dir)
349 s = store.stats()
351 except gf.StoreError as e:
352 die(e)
354 for k in store.stats_keys:
355 print('%s: %s' % (k, s[k]))
358def command_check(args):
360 parser, options, args = cl_parse('check', args)
361 store_dir = get_store_dir(args)
363 try:
364 store = gf.Store(store_dir)
365 problems = store.check(show_progress=True)
366 if problems:
367 die('problems detected with gf store: %s' % store_dir)
369 except gf.StoreError as e:
370 die(e)
373def load_config(fn):
374 try:
375 config = gf.meta.load(filename=fn)
376 assert isinstance(config, gf.Config)
378 except Exception:
379 die('cannot load gf config from file: %s' % fn)
381 return config
384def command_decimate(args):
386 def setup(parser):
387 parser.add_option(
388 '--config', dest='config_fn', metavar='FILE',
389 help='use modified spacial sampling given in FILE')
391 parser.add_option(
392 '--force', dest='force', action='store_true',
393 help='overwrite existing files')
395 parser, options, args = cl_parse('decimate', args, setup=setup)
396 try:
397 decimate = int(args.pop())
398 except Exception:
399 parser.error('cannot get <factor> argument')
401 store_dir = get_store_dir(args)
403 config = None
404 if options.config_fn:
405 config = load_config(options.config_fn)
407 try:
408 store = gf.Store(store_dir)
409 store.make_decimated(decimate, config=config, force=options.force,
410 show_progress=True)
412 except gf.StoreError as e:
413 die(e)
416def sindex(args):
417 return '(%s)' % ', '.join('%g' % x for x in args)
420def command_redeploy(args):
422 parser, options, args = cl_parse('redeploy', args)
424 if not len(args) == 2:
425 sys.exit(parser.format_help())
427 source_store_dir, dest_store_dir = args
429 try:
430 source = gf.Store(source_store_dir)
431 except gf.StoreError as e:
432 die(e)
434 try:
435 gf.store.Store.create_dependants(dest_store_dir)
436 except gf.StoreError:
437 pass
439 try:
440 dest = gf.Store(dest_store_dir, 'w')
442 except gf.StoreError as e:
443 die(e)
445 show_progress = True
447 if show_progress:
448 pbar = util.progressbar('redeploying', dest.config.nrecords)
450 for i, args in enumerate(dest.config.iter_nodes()):
451 try:
452 tr = source.get(args, interpolation='off')
453 dest.put(args, tr)
455 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e:
456 logger.debug('skipping %s, (%s)' % (sindex(args), e))
458 except gf.store.StoreError as e:
459 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
461 if show_progress:
462 pbar.update(i+1)
464 if show_progress:
465 pbar.finish()
468def command_view(args):
469 def setup(parser):
470 parser.add_option('--extract',
471 dest='extract',
472 metavar='start:stop[:step|@num],...',
473 help='specify which traces to show')
475 parser.add_option('--show-phases',
476 dest='showphases',
477 default=None,
478 metavar='[phase_id1,phase_id2,...|all]',
479 help='add phase markers from ttt')
481 parser.add_option('--qt5',
482 dest='gui_toolkit_qt5',
483 action='store_true',
484 default=False,
485 help='use Qt5 for the GUI')
487 parser.add_option('--qt4',
488 dest='gui_toolkit_qt4',
489 action='store_true',
490 default=False,
491 help='use Qt4 for the GUI')
493 parser.add_option('--opengl',
494 dest='opengl',
495 action='store_true',
496 default=False,
497 help='use OpenGL for drawing')
499 parser, options, args = cl_parse('view', args, setup=setup)
501 if options.gui_toolkit_qt4:
502 config.override_gui_toolkit = 'qt4'
504 if options.gui_toolkit_qt5:
505 config.override_gui_toolkit = 'qt5'
507 gdef = None
508 if options.extract:
509 try:
510 gdef = gf.meta.parse_grid_spec(options.extract)
511 except gf.meta.GridSpecError as e:
512 die(e)
514 store_dirs = get_store_dirs(args)
516 alpha = 'abcdefghijklmnopqrstxyz'.upper()
518 markers = []
519 traces = []
521 try:
522 for istore, store_dir in enumerate(store_dirs):
523 store = gf.Store(store_dir)
524 if options.showphases == 'all':
525 phasenames = [pn.id for pn in store.config.tabulated_phases]
526 elif options.showphases is not None:
527 phasenames = options.showphases.split(',')
529 ii = 0
530 for args in store.config.iter_extraction(gdef):
531 gtr = store.get(args)
533 loc_code = ''
534 if len(store_dirs) > 1:
535 loc_code = alpha[istore % len(alpha)]
537 if gtr:
539 sta_code = '%04i (%s)' % (
540 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
541 tmin = gtr.deltat * gtr.itmin
543 tr = trace.Trace(
544 '',
545 sta_code,
546 loc_code,
547 '%02i' % args[-1],
548 ydata=gtr.data,
549 deltat=gtr.deltat,
550 tmin=tmin)
552 if options.showphases:
553 for phasename in phasenames:
554 phase_tmin = store.t(phasename, args[:-1])
555 if phase_tmin:
556 m = marker.PhaseMarker(
557 [('',
558 sta_code,
559 loc_code,
560 '%02i' % args[-1])],
561 phase_tmin,
562 phase_tmin,
563 0,
564 phasename=phasename)
565 markers.append(m)
567 traces.append(tr)
569 ii += 1
571 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
572 die(e)
574 trace.snuffle(traces, markers=markers, opengl=options.opengl)
577def command_extract(args):
578 def setup(parser):
579 parser.add_option(
580 '--format', dest='format', default='mseed',
581 choices=['mseed', 'sac', 'text', 'yaff'],
582 help='export to format "mseed", "sac", "text", or "yaff". '
583 'Default is "mseed".')
585 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
586 parser.add_option(
587 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
588 help='output path template [default: "%s"]' % fndfl)
590 parser, options, args = cl_parse('extract', args, setup=setup)
591 try:
592 sdef = args.pop()
593 except Exception:
594 parser.error('cannot get <selection> argument')
596 try:
597 gdef = gf.meta.parse_grid_spec(sdef)
598 except gf.meta.GridSpecError as e:
599 die(e)
601 store_dir = get_store_dir(args)
603 extensions = {
604 'mseed': 'mseed',
605 'sac': 'sac',
606 'text': 'txt',
607 'yaff': 'yaff'}
609 try:
610 store = gf.Store(store_dir)
611 for args in store.config.iter_extraction(gdef):
612 gtr = store.get(args)
613 if gtr:
614 tr = trace.Trace(
615 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
616 ydata=gtr.data,
617 deltat=gtr.deltat,
618 tmin=gtr.deltat * gtr.itmin)
620 additional = dict(
621 args='_'.join('%g' % x for x in args),
622 irecord=store.str_irecord(args),
623 extension=extensions[options.format])
625 io.save(
626 tr,
627 options.output_fn,
628 format=options.format,
629 additional=additional)
631 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
632 die(e)
635def command_import(args):
636 try:
637 from tunguska import gfdb
638 except ImportError:
639 die('the kiwi tools must be installed to use this feature')
641 parser, options, args = cl_parse('import', args)
643 show_progress = True
645 if not len(args) == 2:
646 sys.exit(parser.format_help())
648 source_path, dest_store_dir = args
650 if op.isdir(source_path):
651 source_path = op.join(source_path, 'db')
653 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
655 db = gfdb.Gfdb(source_path)
657 config = gf.meta.ConfigTypeA(
658 id='imported_gfs',
659 distance_min=db.firstx,
660 distance_max=db.firstx + (db.nx-1) * db.dx,
661 distance_delta=db.dx,
662 source_depth_min=db.firstz,
663 source_depth_max=db.firstz + (db.nz-1) * db.dz,
664 source_depth_delta=db.dz,
665 sample_rate=1.0/db.dt,
666 ncomponents=db.ng
667 )
669 try:
670 gf.store.Store.create(dest_store_dir, config=config)
671 dest = gf.Store(dest_store_dir, 'w')
672 if show_progress:
673 pbar = util.progressbar(
674 'importing', dest.config.nrecords/dest.config.ncomponents)
676 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
677 source_depth, distance = [float(x) for x in args]
678 traces = db.get_traces_pyrocko(distance, source_depth)
679 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
680 for ig in range(db.ng):
681 if ig in ig_to_trace:
682 tr = ig_to_trace[ig]
683 gf_tr = gf.store.GFTrace(
684 tr.get_ydata(),
685 int(round(tr.tmin / tr.deltat)),
686 tr.deltat)
688 else:
689 gf_tr = gf.store.Zero
691 dest.put((source_depth, distance, ig), gf_tr)
693 if show_progress:
694 pbar.update(i+1)
696 if show_progress:
697 pbar.finish()
699 dest.close()
701 except gf.StoreError as e:
702 die(e)
705def command_export(args):
706 from subprocess import Popen, PIPE
708 try:
709 from tunguska import gfdb
710 except ImportError as err:
711 die('the kiwi tools must be installed to use this feature', err)
713 def setup(parser):
714 parser.add_option(
715 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
716 help='split output gfdb into N chunks')
718 parser, options, args = cl_parse('export', args, setup=setup)
720 show_progress = True
722 if len(args) not in (1, 2):
723 sys.exit(parser.format_help())
725 target_path = args.pop()
726 if op.isdir(target_path):
727 target_path = op.join(target_path, 'kiwi_gfdb')
728 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
730 source_store_dir = get_store_dir(args)
732 source = gf.Store(source_store_dir, 'r')
733 config = source.config
735 if not isinstance(config, gf.meta.ConfigTypeA):
736 die('only stores of type A can be exported to Kiwi format')
738 if op.isfile(target_path + '.index'):
739 die('destation already exists')
741 cmd = [str(x) for x in [
742 'gfdb_build',
743 target_path,
744 options.nchunks,
745 config.ndistances,
746 config.nsource_depths,
747 config.ncomponents,
748 config.deltat,
749 config.distance_delta,
750 config.source_depth_delta,
751 config.distance_min,
752 config.source_depth_min]]
754 p = Popen(cmd, stdin=PIPE)
755 p.communicate()
757 out_db = gfdb.Gfdb(target_path)
759 if show_progress:
760 pbar = util.progressbar(
761 'exporting', config.nrecords/config.ncomponents)
763 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
765 data_out = []
766 for ig in range(config.ncomponents):
767 try:
768 tr = source.get((z, x, ig), interpolation='off')
769 data_out.append((tr.t, tr.data * config.factor))
771 except gf.store.StoreError as e:
772 logger.warning('cannot get %s, (%s)' % (sindex((z, x, ig)), e))
773 data_out.append(None)
775 # put a zero valued sample to no-data zero-traces at a compatible time
776 tmins = [
777 entry[0][0]
778 for entry in data_out
779 if entry is not None and entry[0].size != 0]
781 if tmins:
782 tmin = min(tmins)
783 for entry in data_out:
784 if entry is not None and entry[0].size == 0:
785 entry[0].resize(1)
786 entry[1].resize(1)
787 entry[0][0] = tmin
788 entry[1][0] = 0.0
790 out_db.put_traces_slow(x, z, data_out)
792 if show_progress:
793 pbar.update(i+1)
795 if show_progress:
796 pbar.finish()
798 source.close()
801def phasedef_or_horvel(x):
802 try:
803 return float(x)
804 except ValueError:
805 return cake.PhaseDef(x)
808def mkp(s):
809 return [phasedef_or_horvel(ps) for ps in s.split(',')]
812def stored_attribute_table_plots(phase_ids, options, args, attribute):
813 import numpy as num
814 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
816 plt = mpl_init()
818 np = 1
819 store_dir = get_store_dir(args)
820 for phase_id in phase_ids:
821 try:
822 store = gf.Store(store_dir)
824 if options.receiver_depth is not None:
825 receiver_depth = options.receiver_depth * 1000.0
826 else:
827 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
828 receiver_depth = store.config.receiver_depth
830 elif isinstance(store.config, gf.ConfigTypeB):
831 receiver_depth = store.config.receiver_depth_min
833 else:
834 receiver_depth = 0.0
836 phase = store.get_stored_phase(phase_id, attribute)
837 axes = plt.subplot(2, len(phase_ids), np)
838 labelspace(axes)
839 xscaled(1. / km, axes)
840 yscaled(1. / km, axes)
841 x = None
842 if isinstance(store.config, gf.ConfigTypeB):
843 x = (receiver_depth, None, None)
845 phase.plot_2d(axes, x=x)
846 axes.set_title(phase_id)
847 np += 1
848 except gf.StoreError as e:
849 die(e)
851 axes = plt.subplot(2, 1, 2)
852 num_d = 100
853 distances = num.linspace(store.config.distance_min,
854 store.config.distance_max,
855 num_d)
857 if options.source_depth is not None:
858 source_depth = options.source_depth * km
859 else:
860 source_depth = store.config.source_depth_min + (
861 store.config.source_depth_max - store.config.source_depth_min) / 2.
863 if isinstance(store.config, gf.ConfigTypeA):
864 attribute_vals = num.empty(num_d)
865 for phase_id in phase_ids:
866 attribute_vals[:] = num.NAN
867 for i, d in enumerate(distances):
868 if attribute == 'phase':
869 attribute_vals[i] = store.t(phase_id, (source_depth, d))
870 ylabel = 'TT [s]'
871 else:
872 attribute_vals[i] = store.get_stored_attribute(
873 phase_id, options.attribute, (source_depth, d))
874 ylabel = '%s [deg]' % options.attribute
876 axes.plot(distances / km, attribute_vals, label=phase_id)
878 axes.set_title('source source_depth %s km' % (source_depth / km))
879 axes.set_xlabel('distance [km]')
880 axes.set_ylabel(ylabel)
881 axes.legend()
883 plt.tight_layout()
884 mpl_show(plt)
887def command_ttt(args):
888 def setup(parser):
889 parser.add_option(
890 '--force', dest='force', action='store_true',
891 help='overwrite existing files')
893 parser, options, args = cl_parse('ttt', args, setup=setup)
895 store_dir = get_store_dir(args)
896 try:
897 store = gf.Store(store_dir)
898 store.make_travel_time_tables(force=options.force)
900 except gf.StoreError as e:
901 die(e)
904def command_tttview(args):
906 def setup(parser):
907 parser.add_option(
908 '--source-depth', dest='source_depth', type=float,
909 help='Source depth in km')
911 parser.add_option(
912 '--receiver-depth', dest='receiver_depth', type=float,
913 help='Receiver depth in km')
915 parser, options, args = cl_parse(
916 'tttview', args, setup=setup,
917 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
919 try:
920 phase_ids = args.pop().split(',')
921 except Exception:
922 parser.error('cannot get <phase-ids> argument')
924 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
927def command_sat(args):
928 def setup(parser):
929 parser.add_option(
930 '--force', dest='force', action='store_true',
931 help='overwrite existing files')
933 parser.add_option(
934 '--attribute',
935 action='store',
936 dest='attribute',
937 type='choice',
938 choices=gf.store.available_stored_tables[1::],
939 default='takeoff_angle',
940 help='calculate interpolation table for selected ray attributes.')
942 parser, options, args = cl_parse('sat', args, setup=setup)
944 store_dir = get_store_dir(args)
945 try:
946 store = gf.Store(store_dir)
947 store.make_stored_table(options.attribute, force=options.force)
949 except gf.StoreError as e:
950 die(e)
953def command_satview(args):
955 def setup(parser):
956 parser.add_option(
957 '--source-depth', dest='source_depth', type=float,
958 help='Source depth in km')
960 parser.add_option(
961 '--receiver-depth', dest='receiver_depth', type=float,
962 help='Receiver depth in km')
964 parser.add_option(
965 '--attribute',
966 action='store',
967 dest='attribute',
968 type='choice',
969 choices=gf.store.available_stored_tables[1::],
970 default='takeoff_angle',
971 help='view selected ray attribute.')
973 parser, options, args = cl_parse(
974 'satview', args, setup=setup,
975 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
977 try:
978 phase_ids = args.pop().split(',')
979 except Exception:
980 parser.error('cannot get <phase-ids> argument')
982 logger.info('Plotting stored attribute %s' % options.attribute)
984 stored_attribute_table_plots(
985 phase_ids, options, args, attribute=options.attribute)
988def command_tttextract(args):
989 def setup(parser):
990 parser.add_option(
991 '--output', dest='output_fn', metavar='TEMPLATE',
992 help='output to text files instead of stdout '
993 '(example TEMPLATE: "extracted/%(args)s.txt")')
995 parser, options, args = cl_parse('tttextract', args, setup=setup)
996 try:
997 sdef = args.pop()
998 except Exception:
999 parser.error('cannot get <selection> argument')
1001 try:
1002 sphase = args.pop()
1003 except Exception:
1004 parser.error('cannot get <phase> argument')
1006 try:
1007 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1008 except gf.meta.InvalidTimingSpecification:
1009 parser.error('invalid phase specification: "%s"' % sphase)
1011 try:
1012 gdef = gf.meta.parse_grid_spec(sdef)
1013 except gf.meta.GridSpecError as e:
1014 die(e)
1016 store_dir = get_store_dir(args)
1018 try:
1019 store = gf.Store(store_dir)
1020 for args in store.config.iter_extraction(gdef, level=-1):
1021 s = ['%e' % x for x in args]
1022 for phase in phases:
1023 t = store.t(phase, args)
1024 if t is not None:
1025 s.append('%e' % t)
1026 else:
1027 s.append('nan')
1029 if options.output_fn:
1030 d = dict(
1031 args='_'.join('%e' % x for x in args),
1032 extension='txt')
1034 fn = options.output_fn % d
1035 util.ensuredirs(fn)
1036 with open(fn, 'a') as f:
1037 f.write(' '.join(s))
1038 f.write('\n')
1039 else:
1040 print(' '.join(s))
1042 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1043 die(e)
1046def command_tttlsd(args):
1048 def setup(parser):
1049 pass
1051 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1053 try:
1054 sphase_ids = args.pop()
1055 except Exception:
1056 parser.error('cannot get <phase> argument')
1058 try:
1059 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1060 except gf.meta.InvalidTimingSpecification:
1061 parser.error('invalid phase specification: "%s"' % sphase_ids)
1063 store_dir = get_store_dir(args)
1065 try:
1066 store = gf.Store(store_dir)
1067 for phase_id in phase_ids:
1068 store.fix_ttt_holes(phase_id)
1070 except gf.StoreError as e:
1071 die(e)
1074def command_server(args):
1075 from pyrocko.gf import server
1077 def setup(parser):
1078 parser.add_option(
1079 '--port', dest='port', metavar='PORT', type='int', default=8080,
1080 help='serve on port PORT')
1082 parser.add_option(
1083 '--ip', dest='ip', metavar='IP', default='',
1084 help='serve on ip address IP')
1086 parser, options, args = cl_parse('server', args, setup=setup)
1088 engine = gf.LocalEngine(store_superdirs=args)
1089 server.run(options.ip, options.port, engine)
1092def command_download(args):
1093 from pyrocko.gf import ws
1095 details = '''
1097Browse pre-calculated Green's function stores online:
1099 https://greens-mill.pyrocko.org
1100'''
1102 def setup(parser):
1103 parser.add_option(
1104 '--force', dest='force', action='store_true',
1105 help='overwrite existing files')
1107 parser, options, args = cl_parse(
1108 'download', args, setup=setup, details=details)
1109 if not len(args) in (1, 2):
1110 sys.exit(parser.format_help())
1112 if len(args) == 2:
1113 site, store_id = args
1114 if not re.match(gf.meta.StringID.pattern, store_id):
1115 die('invalid store ID')
1116 else:
1117 site, store_id = args[0], None
1119 if site not in gf.ws.g_site_abbr:
1120 if -1 == site.find('://'):
1121 site = 'http://' + site
1123 try:
1124 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1125 except ws.DownloadError as e:
1126 die(str(e))
1129def command_modelview(args):
1131 import matplotlib.pyplot as plt
1132 import numpy as num
1133 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1134 mpl_init()
1136 neat_labels = {
1137 'vp': '$v_p$',
1138 'vs': '$v_s$',
1139 'qp': '$Q_p$',
1140 'qs': '$Q_s$',
1141 'rho': '$\\rho$'}
1143 def setup(parser):
1144 parser.add_option(
1145 '--parameters', dest='parameters',
1146 default='vp,vs', metavar='vp,vs,....',
1147 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1149 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1151 parser, options, args = cl_parse('modelview', args, setup=setup)
1153 store_dirs = get_store_dirs(args)
1155 parameters = options.parameters.split(',')
1157 fig, axes = plt.subplots(1,
1158 len(parameters),
1159 sharey=True,
1160 figsize=(len(parameters)*3, 5))
1162 if not isinstance(axes, num.ndarray):
1163 axes = [axes]
1165 axes = dict(zip(parameters, axes))
1167 for store_id in store_dirs:
1168 try:
1169 store = gf.Store(store_id)
1170 mod = store.config.earthmodel_1d
1171 z = mod.profile('z')
1173 for p in parameters:
1174 ax = axes[p]
1175 labelspace(ax)
1176 if '/' in p:
1177 p1, p2 = p.split('/')
1178 profile = mod.profile(p1)/mod.profile(p2)
1179 else:
1180 profile = mod.profile(p)
1182 ax.plot(profile, -z, label=store_id.split('/')[-1])
1183 if p in ['vs', 'vp', 'rho']:
1184 xscaled(1./1000., ax)
1186 yscaled(1./km, ax)
1188 except gf.StoreError as e:
1189 die(e)
1191 for p, ax in axes.items():
1192 ax.grid()
1193 if p in neat_labels:
1194 lab = neat_labels[p]
1195 elif all(x in neat_labels for x in p.split('/')):
1196 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1197 else:
1198 lab = p
1200 ax.set_title(lab, y=1.02)
1202 if p in units:
1203 lab += ' ' + units[p]
1205 ax.autoscale()
1206 ax.set_xlabel(lab)
1208 axes[parameters[0]].set_ylabel('Depth [km]')
1210 handles, labels = ax.get_legend_handles_labels()
1212 if len(store_dirs) > 1:
1213 ax.legend(
1214 handles,
1215 labels,
1216 bbox_to_anchor=(0.5, 0.12),
1217 bbox_transform=fig.transFigure,
1218 ncol=3,
1219 loc='upper center',
1220 fancybox=True)
1222 plt.subplots_adjust(bottom=0.22,
1223 wspace=0.05)
1225 mpl_show(plt)
1228def command_upgrade(args):
1229 parser, options, args = cl_parse('upgrade', args)
1230 store_dirs = get_store_dirs(args)
1231 try:
1232 for store_dir in store_dirs:
1233 store = gf.Store(store_dir)
1234 nup = store.upgrade()
1235 if nup == 0:
1236 print('%s: already up-to-date.' % store_dir)
1237 else:
1238 print('%s: %i file%s upgraded.' % (
1239 store_dir, nup, ['s', ''][nup == 1]))
1241 except gf.StoreError as e:
1242 die(e)
1245def command_addref(args):
1246 parser, options, args = cl_parse('addref', args)
1248 store_dirs = []
1249 filenames = []
1250 for arg in args:
1251 if op.isdir(arg):
1252 store_dirs.append(arg)
1253 elif op.isfile(arg):
1254 filenames.append(arg)
1255 else:
1256 die('invalid argument: %s' % arg)
1258 if not store_dirs:
1259 store_dirs.append('.')
1261 references = []
1262 try:
1263 for filename in args:
1264 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1265 except ImportError:
1266 die('pybtex module must be installed to use this function')
1268 if not references:
1269 die('no references found')
1271 for store_dir in store_dirs:
1272 try:
1273 store = gf.Store(store_dir)
1274 ids = [ref.id for ref in store.config.references]
1275 for ref in references:
1276 if ref.id in ids:
1277 die('duplicate reference id: %s' % ref.id)
1279 ids.append(ref.id)
1280 store.config.references.append(ref)
1282 store.save_config(make_backup=True)
1284 except gf.StoreError as e:
1285 die(e)
1288def command_qc(args):
1289 parser, options, args = cl_parse('qc', args)
1291 store_dir = get_store_dir(args)
1293 try:
1294 store = gf.Store(store_dir)
1295 s = store.stats()
1296 if s['empty'] != 0:
1297 print('has empty records')
1299 for aname in ['author', 'author_email', 'description', 'references']:
1301 if not getattr(store.config, aname):
1302 print('%s empty' % aname)
1304 except gf.StoreError as e:
1305 die(e)
1308def command_report(args):
1309 from pyrocko.fomosto.report import report_call
1310 report_call.run_program(args)
1313def main(args=None):
1314 if args is None:
1315 args = sys.argv[1:]
1317 if len(args) < 1:
1318 sys.exit('Usage: %s' % usage)
1320 command = args.pop(0)
1322 if command in subcommands:
1323 globals()['command_' + command](args)
1325 elif command in ('--help', '-h', 'help'):
1326 if command == 'help' and args:
1327 acommand = args[0]
1328 if acommand in subcommands:
1329 globals()['command_' + acommand](['--help'])
1331 sys.exit('Usage: %s' % usage)
1333 else:
1334 sys.exit('fomosto: error: no such subcommand: %s' % command)
1337if __name__ == '__main__':
1338 main(sys.argv[1:])