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, 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 try:
448 if show_progress:
449 pbar = util.progressbar('redeploying', dest.config.nrecords)
451 for i, args in enumerate(dest.config.iter_nodes()):
452 try:
453 tr = source.get(args, interpolation='off')
454 dest.put(args, tr)
456 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa
457 logger.debug('skipping %s, (%s)' % (sindex(args), e))
459 except gf.store.StoreError as e:
460 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
462 if show_progress:
463 pbar.update(i+1)
465 finally:
466 if show_progress:
467 pbar.finish()
470def command_view(args):
471 def setup(parser):
472 parser.add_option(
473 '--extract',
474 dest='extract',
475 metavar='start:stop[:step|@num],...',
476 help='specify which traces to show')
478 parser.add_option(
479 '--show-phases',
480 dest='showphases',
481 default=None,
482 metavar='[phase_id1,phase_id2,...|all]',
483 help='add phase markers from ttt')
485 parser.add_option(
486 '--opengl',
487 dest='opengl',
488 action='store_true',
489 default=None,
490 help='use OpenGL for drawing')
492 parser.add_option(
493 '--no-opengl',
494 dest='opengl',
495 action='store_false',
496 default=None,
497 help='do not use OpenGL for drawing')
499 parser, options, args = cl_parse('view', args, setup=setup)
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 try:
667 if show_progress:
668 pbar = util.progressbar(
669 'importing', dest.config.nrecords/dest.config.ncomponents)
671 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
672 source_depth, distance = [float(x) for x in args]
673 traces = db.get_traces_pyrocko(distance, source_depth)
674 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
675 for ig in range(db.ng):
676 if ig in ig_to_trace:
677 tr = ig_to_trace[ig]
678 gf_tr = gf.store.GFTrace(
679 tr.get_ydata(),
680 int(round(tr.tmin / tr.deltat)),
681 tr.deltat)
683 else:
684 gf_tr = gf.store.Zero
686 dest.put((source_depth, distance, ig), gf_tr)
688 if show_progress:
689 pbar.update(i+1)
691 finally:
692 if show_progress:
693 pbar.finish()
695 dest.close()
697 except gf.StoreError as e:
698 die(e)
701def command_export(args):
702 from subprocess import Popen, PIPE
704 try:
705 from tunguska import gfdb
706 except ImportError as err:
707 die('the kiwi tools must be installed to use this feature', err)
709 def setup(parser):
710 parser.add_option(
711 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
712 help='split output gfdb into N chunks')
714 parser, options, args = cl_parse('export', args, setup=setup)
716 show_progress = True
718 if len(args) not in (1, 2):
719 sys.exit(parser.format_help())
721 target_path = args.pop()
722 if op.isdir(target_path):
723 target_path = op.join(target_path, 'kiwi_gfdb')
724 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
726 source_store_dir = get_store_dir(args)
728 source = gf.Store(source_store_dir, 'r')
729 config = source.config
731 if not isinstance(config, gf.meta.ConfigTypeA):
732 die('only stores of type A can be exported to Kiwi format')
734 if op.isfile(target_path + '.index'):
735 die('destation already exists')
737 cmd = [str(x) for x in [
738 'gfdb_build',
739 target_path,
740 options.nchunks,
741 config.ndistances,
742 config.nsource_depths,
743 config.ncomponents,
744 config.deltat,
745 config.distance_delta,
746 config.source_depth_delta,
747 config.distance_min,
748 config.source_depth_min]]
750 p = Popen(cmd, stdin=PIPE)
751 p.communicate()
753 out_db = gfdb.Gfdb(target_path)
755 try:
756 if show_progress:
757 pbar = util.progressbar(
758 'exporting', config.nrecords/config.ncomponents)
760 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
762 data_out = []
763 for ig in range(config.ncomponents):
764 try:
765 tr = source.get((z, x, ig), interpolation='off')
766 data_out.append((tr.t, tr.data * config.factor))
768 except gf.store.StoreError as e:
769 logger.warning(
770 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
771 data_out.append(None)
773 # put a zero valued sample to no-data zero-traces at a compatible
774 # time
775 tmins = [
776 entry[0][0]
777 for entry in data_out
778 if entry is not None and entry[0].size != 0]
780 if tmins:
781 tmin = min(tmins)
782 for entry in data_out:
783 if entry is not None and entry[0].size == 0:
784 entry[0].resize(1)
785 entry[1].resize(1)
786 entry[0][0] = tmin
787 entry[1][0] = 0.0
789 out_db.put_traces_slow(x, z, data_out)
791 if show_progress:
792 pbar.update(i+1)
794 source.close()
796 finally:
797 if show_progress:
798 pbar.finish()
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:])