1#!/usr/bin/env python
2# http://pyrocko.org - GPLv3
3#
4# The Pyrocko Developers, 21st Century
5# ---|P------/S----------~Lg----------
7import sys
8import re
9import os.path as op
10import logging
11import copy
12import shutil
13from optparse import OptionParser
15from pyrocko import util, trace, gf, cake, io, fomosto
16from pyrocko.gui import marker
17from pyrocko.util import mpl_show
19logger = logging.getLogger('pyrocko.apps.fomosto')
20km = 1e3
23def d2u(d):
24 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
27subcommand_descriptions = {
28 'init': 'create a new empty GF store',
29 'build': 'compute GFs and fill into store',
30 'stats': 'print information about a GF store',
31 'check': 'check for problems in GF store',
32 'decimate': 'build decimated variant of a GF store',
33 'redeploy': 'copy traces from one GF store into another',
34 'view': 'view selected traces',
35 'extract': 'extract selected traces',
36 'import': 'convert Kiwi GFDB to GF store format',
37 'export': 'convert GF store to Kiwi GFDB format',
38 'ttt': 'create travel time tables',
39 'tttview': 'plot travel time table',
40 'tttextract': 'extract selected travel times',
41 'tttlsd': 'fix holes in travel time tables',
42 'sat': 'create stored ray attribute table',
43 'satview': 'plot stored ray attribute table',
44 'server': 'run seismosizer server',
45 'download': 'download GF store from a server',
46 'modelview': 'plot earthmodels',
47 'upgrade': 'upgrade store format to latest version',
48 'addref': 'import citation references to GF store config',
49 'qc': 'quality check',
50 'report': 'report for Green\'s Function databases',
51}
53subcommand_usages = {
54 'init': ['init <type> <store-dir> [options]',
55 'init redeploy <source> <destination> [options]'],
56 'build': 'build [store-dir] [options]',
57 'stats': 'stats [store-dir] [options]',
58 'check': 'check [store-dir] [options]',
59 'decimate': 'decimate [store-dir] <factor> [options]',
60 'redeploy': 'redeploy <source> <destination> [options]',
61 'view': 'view [store-dir] ... [options]',
62 'extract': 'extract [store-dir] <selection>',
63 'import': 'import <source> <destination> [options]',
64 'export': 'export [store-dir] <destination> [options]',
65 'ttt': 'ttt [store-dir] [options]',
66 'tttview': 'tttview [store-dir] <phase-ids> [options]',
67 'tttextract': 'tttextract [store-dir] <phase> <selection>',
68 'tttlsd': 'tttlsd [store-dir] <phase>',
69 'sat': 'sat [store-dir] [options]',
70 'satview': 'satview [store-dir] <phase-ids> [options]',
71 'server': 'server [options] <store-super-dir> ...',
72 'download': 'download [options] <site> <store-id>',
73 'modelview': 'modelview <selection>',
74 'upgrade': 'upgrade [store-dir] ...',
75 'addref': 'addref [store-dir] ... <filename> ...',
76 'qc': 'qc [store-dir]',
77 'report': 'report <subcommand> <arguments>... [options]'
78}
80subcommands = subcommand_descriptions.keys()
82program_name = 'fomosto'
84usage = program_name + ''' <subcommand> <arguments> ... [options]
86Subcommands:
88 init %(init)s
89 build %(build)s
90 stats %(stats)s
91 check %(check)s
92 decimate %(decimate)s
93 redeploy %(redeploy)s
94 view %(view)s
95 extract %(extract)s
96 import %(import)s
97 export %(export)s
98 ttt %(ttt)s
99 tttview %(tttview)s
100 tttextract %(tttextract)s
101 tttlsd %(tttlsd)s
102 sat %(sat)s
103 satview %(satview)s
104 server %(server)s
105 download %(download)s
106 modelview %(modelview)s
107 upgrade %(upgrade)s
108 addref %(addref)s
109 qc %(qc)s
110 report %(report)s
112To get further help and a list of available options for any subcommand run:
114 fomosto <subcommand> --help
116''' % d2u(subcommand_descriptions)
119def add_common_options(parser):
120 parser.add_option(
121 '--loglevel',
122 action='store',
123 dest='loglevel',
124 type='choice',
125 choices=('critical', 'error', 'warning', 'info', 'debug'),
126 default='info',
127 help='set logger level to '
128 '"critical", "error", "warning", "info", or "debug". '
129 'Default is "%default".')
132def process_common_options(options):
133 util.setup_logging(program_name, options.loglevel)
136def cl_parse(command, args, setup=None, details=None):
137 usage = subcommand_usages[command]
138 descr = subcommand_descriptions[command]
140 if isinstance(usage, str):
141 usage = [usage]
143 susage = '%s %s' % (program_name, usage[0])
144 for s in usage[1:]:
145 susage += '\n%s%s %s' % (' '*7, program_name, s)
147 description = descr[0].upper() + descr[1:] + '.'
149 if details:
150 description = description + ' %s' % details
152 parser = OptionParser(usage=susage, description=description)
153 parser.format_description = lambda formatter: description
155 if setup:
156 setup(parser)
158 add_common_options(parser)
159 (options, args) = parser.parse_args(args)
160 process_common_options(options)
161 return parser, options, args
164def die(message, err=''):
165 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
168def fomo_wrapper_module(name):
169 try:
170 if not re.match(gf.meta.StringID.pattern, name):
171 raise ValueError('invalid name')
173 words = name.split('.', 1)
174 if len(words) == 2:
175 name, variant = words
176 else:
177 name = words[0]
178 variant = None
180 name_clean = re.sub(r'[.-]', '_', name)
181 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
182 mod = __import__(modname, level=0)
183 return getattr(mod.fomosto, name_clean), variant
185 except ValueError:
186 die('invalid modelling code wrapper name')
188 except ImportError:
189 die('''modelling code wrapper "%s" not available or not installed
190 (module probed: "%s")''' % (name, modname))
193def command_init(args):
195 details = '''
197Available modelling backends:
198%s
200 More information at
201 https://pyrocko.org/docs/current/apps/fomosto/backends.html
202''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
204 parser, options, args = cl_parse(
205 'init', args,
206 details=details)
208 if len(args) == 0:
209 sys.exit(parser.format_help())
211 if args[0] == 'redeploy':
212 if len(args) != 3:
213 parser.error('incorrect number of arguments')
215 source_dir, dest_dir = args[1:]
217 try:
218 source = gf.Store(source_dir)
219 except gf.StoreError as e:
220 die(e)
222 config = copy.deepcopy(source.config)
223 config.derived_from_id = source.config.id
224 try:
225 config_filenames = gf.store.Store.create_editables(
226 dest_dir, config=config)
228 except gf.StoreError as e:
229 die(e)
231 try:
232 dest = gf.Store(dest_dir)
233 except gf.StoreError as e:
234 die(e)
236 for k in source.extra_keys():
237 source_fn = source.get_extra_path(k)
238 dest_fn = dest.get_extra_path(k)
239 shutil.copyfile(source_fn, dest_fn)
241 logger.info(
242 '(1) configure settings in files:\n %s'
243 % '\n '.join(config_filenames))
245 logger.info(
246 '(2) run "fomosto redeploy <source> <dest>", as needed')
248 else:
249 if len(args) != 2:
250 parser.error('incorrect number of arguments')
252 (modelling_code_id, store_dir) = args
254 module, variant = fomo_wrapper_module(modelling_code_id)
255 try:
256 config_filenames = module.init(store_dir, variant)
257 except gf.StoreError as e:
258 die(e)
260 logger.info('(1) configure settings in files:\n %s'
261 % '\n '.join(config_filenames))
262 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
263 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
266def get_store_dir(args):
267 if len(args) == 1:
268 store_dir = op.abspath(args.pop(0))
269 else:
270 store_dir = op.abspath(op.curdir)
272 if not op.isdir(store_dir):
273 die('not a directory: %s' % store_dir)
275 return store_dir
278def get_store_dirs(args):
279 if len(args) == 0:
280 store_dirs = [op.abspath(op.curdir)]
281 else:
282 store_dirs = [op.abspath(x) for x in args]
284 for store_dir in store_dirs:
285 if not op.isdir(store_dir):
286 die('not a directory: %s' % store_dir)
288 return store_dirs
291def command_build(args):
293 def setup(parser):
294 parser.add_option(
295 '--force', dest='force', action='store_true',
296 help='overwrite existing files')
298 parser.add_option(
299 '--nworkers', dest='nworkers', type='int', metavar='N',
300 help='run N worker processes in parallel')
302 parser.add_option(
303 '--continue', dest='continue_', action='store_true',
304 help='continue suspended build')
306 parser.add_option(
307 '--step', dest='step', type='int', metavar='I',
308 help='process block number IBLOCK')
310 parser.add_option(
311 '--block', dest='iblock', type='int', metavar='I',
312 help='process block number IBLOCK')
314 parser, options, args = cl_parse('build', args, setup=setup)
316 store_dir = get_store_dir(args)
317 try:
318 if options.step is not None:
319 step = options.step - 1
320 else:
321 step = None
323 if options.iblock is not None:
324 iblock = options.iblock - 1
325 else:
326 iblock = None
328 store = gf.Store(store_dir)
329 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
330 module.build(
331 store_dir,
332 force=options.force,
333 nworkers=options.nworkers, continue_=options.continue_,
334 step=step,
335 iblock=iblock)
337 except gf.StoreError as e:
338 die(e)
341def command_stats(args):
343 parser, options, args = cl_parse('stats', args)
344 store_dir = get_store_dir(args)
346 try:
347 store = gf.Store(store_dir)
348 s = store.stats()
350 except gf.StoreError as e:
351 die(e)
353 for k in store.stats_keys:
354 print('%s: %s' % (k, s[k]))
357def command_check(args):
359 parser, options, args = cl_parse('check', args)
360 store_dir = get_store_dir(args)
362 try:
363 store = gf.Store(store_dir)
364 problems = store.check(show_progress=True)
365 if problems:
366 die('problems detected with gf store: %s' % store_dir)
368 except gf.StoreError as e:
369 die(e)
372def load_config(fn):
373 try:
374 config = gf.meta.load(filename=fn)
375 assert isinstance(config, gf.Config)
377 except Exception:
378 die('cannot load gf config from file: %s' % fn)
380 return config
383def command_decimate(args):
385 def setup(parser):
386 parser.add_option(
387 '--config', dest='config_fn', metavar='FILE',
388 help='use modified spacial sampling given in FILE')
390 parser.add_option(
391 '--force', dest='force', action='store_true',
392 help='overwrite existing files')
394 parser, options, args = cl_parse('decimate', args, setup=setup)
395 try:
396 decimate = int(args.pop())
397 except Exception:
398 parser.error('cannot get <factor> argument')
400 store_dir = get_store_dir(args)
402 config = None
403 if options.config_fn:
404 config = load_config(options.config_fn)
406 try:
407 store = gf.Store(store_dir)
408 store.make_decimated(decimate, config=config, force=options.force,
409 show_progress=True)
411 except gf.StoreError as e:
412 die(e)
415def sindex(args):
416 return '(%s)' % ', '.join('%g' % x for x in args)
419def command_redeploy(args):
421 parser, options, args = cl_parse('redeploy', args)
423 if not len(args) == 2:
424 sys.exit(parser.format_help())
426 source_store_dir, dest_store_dir = args
428 try:
429 source = gf.Store(source_store_dir)
430 except gf.StoreError as e:
431 die(e)
433 try:
434 gf.store.Store.create_dependants(dest_store_dir)
435 except gf.StoreError:
436 pass
438 try:
439 dest = gf.Store(dest_store_dir, 'w')
441 except gf.StoreError as e:
442 die(e)
444 show_progress = True
446 try:
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: # noqa
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 finally:
465 if show_progress:
466 pbar.finish()
469def command_view(args):
470 def setup(parser):
471 parser.add_option(
472 '--extract',
473 dest='extract',
474 metavar='start:stop[:step|@num],...',
475 help='specify which traces to show')
477 parser.add_option(
478 '--show-phases',
479 dest='showphases',
480 default=None,
481 metavar='[phase_id1,phase_id2,...|all]',
482 help='add phase markers from ttt')
484 parser.add_option(
485 '--opengl',
486 dest='opengl',
487 action='store_true',
488 default=None,
489 help='use OpenGL for drawing')
491 parser.add_option(
492 '--no-opengl',
493 dest='opengl',
494 action='store_false',
495 default=None,
496 help='do not use OpenGL for drawing')
498 parser, options, args = cl_parse('view', args, setup=setup)
500 gdef = None
501 if options.extract:
502 try:
503 gdef = gf.meta.parse_grid_spec(options.extract)
504 except gf.meta.GridSpecError as e:
505 die(e)
507 store_dirs = get_store_dirs(args)
509 alpha = 'abcdefghijklmnopqrstxyz'.upper()
511 markers = []
512 traces = []
514 try:
515 for istore, store_dir in enumerate(store_dirs):
516 store = gf.Store(store_dir)
517 if options.showphases == 'all':
518 phasenames = [pn.id for pn in store.config.tabulated_phases]
519 elif options.showphases is not None:
520 phasenames = options.showphases.split(',')
522 ii = 0
523 for args in store.config.iter_extraction(gdef):
524 gtr = store.get(args)
526 loc_code = ''
527 if len(store_dirs) > 1:
528 loc_code = alpha[istore % len(alpha)]
530 if gtr:
532 sta_code = '%04i (%s)' % (
533 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
534 tmin = gtr.deltat * gtr.itmin
536 tr = trace.Trace(
537 '',
538 sta_code,
539 loc_code,
540 '%02i' % args[-1],
541 ydata=gtr.data,
542 deltat=gtr.deltat,
543 tmin=tmin)
545 if options.showphases:
546 for phasename in phasenames:
547 phase_tmin = store.t(phasename, args[:-1])
548 if phase_tmin:
549 m = marker.PhaseMarker(
550 [('',
551 sta_code,
552 loc_code,
553 '%02i' % args[-1])],
554 phase_tmin,
555 phase_tmin,
556 0,
557 phasename=phasename)
558 markers.append(m)
560 traces.append(tr)
562 ii += 1
564 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
565 die(e)
567 def setup(win):
568 viewer = win.get_view()
569 viewer.menuitem_showboxes.setChecked(False)
570 viewer.menuitem_colortraces.setChecked(True)
571 viewer.menuitem_demean.setChecked(False)
572 viewer.menuitems_sorting[5][0].setChecked(True)
573 viewer.menuitems_scaling[2][0].setChecked(True)
574 viewer.sortingmode_change()
575 viewer.scalingmode_change()
577 trace.snuffle(
578 traces, markers=markers, opengl=options.opengl, launch_hook=setup)
581def command_extract(args):
582 def setup(parser):
583 parser.add_option(
584 '--format', dest='format', default='mseed',
585 choices=['mseed', 'sac', 'text', 'yaff'],
586 help='export to format "mseed", "sac", "text", or "yaff". '
587 'Default is "mseed".')
589 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
590 parser.add_option(
591 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
592 help='output path template [default: "%s"]' % fndfl)
594 parser, options, args = cl_parse('extract', args, setup=setup)
595 try:
596 sdef = args.pop()
597 except Exception:
598 parser.error('cannot get <selection> argument')
600 try:
601 gdef = gf.meta.parse_grid_spec(sdef)
602 except gf.meta.GridSpecError as e:
603 die(e)
605 store_dir = get_store_dir(args)
607 extensions = {
608 'mseed': 'mseed',
609 'sac': 'sac',
610 'text': 'txt',
611 'yaff': 'yaff'}
613 try:
614 store = gf.Store(store_dir)
615 for args in store.config.iter_extraction(gdef):
616 gtr = store.get(args)
617 if gtr:
618 tr = trace.Trace(
619 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
620 ydata=gtr.data,
621 deltat=gtr.deltat,
622 tmin=gtr.deltat * gtr.itmin)
624 additional = dict(
625 args='_'.join('%g' % x for x in args),
626 irecord=store.str_irecord(args),
627 extension=extensions[options.format])
629 io.save(
630 tr,
631 options.output_fn,
632 format=options.format,
633 additional=additional)
635 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
636 die(e)
639def command_import(args):
640 try:
641 from tunguska import gfdb
642 except ImportError:
643 die('the kiwi tools must be installed to use this feature')
645 parser, options, args = cl_parse('import', args)
647 show_progress = True
649 if not len(args) == 2:
650 sys.exit(parser.format_help())
652 source_path, dest_store_dir = args
654 if op.isdir(source_path):
655 source_path = op.join(source_path, 'db')
657 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
659 db = gfdb.Gfdb(source_path)
661 config = gf.meta.ConfigTypeA(
662 id='imported_gfs',
663 distance_min=db.firstx,
664 distance_max=db.firstx + (db.nx-1) * db.dx,
665 distance_delta=db.dx,
666 source_depth_min=db.firstz,
667 source_depth_max=db.firstz + (db.nz-1) * db.dz,
668 source_depth_delta=db.dz,
669 sample_rate=1.0/db.dt,
670 ncomponents=db.ng
671 )
673 try:
674 gf.store.Store.create(dest_store_dir, config=config)
675 dest = gf.Store(dest_store_dir, 'w')
676 try:
677 if show_progress:
678 pbar = util.progressbar(
679 'importing', dest.config.nrecords/dest.config.ncomponents)
681 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
682 source_depth, distance = [float(x) for x in args]
683 traces = db.get_traces_pyrocko(distance, source_depth)
684 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
685 for ig in range(db.ng):
686 if ig in ig_to_trace:
687 tr = ig_to_trace[ig]
688 gf_tr = gf.store.GFTrace(
689 tr.get_ydata(),
690 int(round(tr.tmin / tr.deltat)),
691 tr.deltat)
693 else:
694 gf_tr = gf.store.Zero
696 dest.put((source_depth, distance, ig), gf_tr)
698 if show_progress:
699 pbar.update(i+1)
701 finally:
702 if show_progress:
703 pbar.finish()
705 dest.close()
707 except gf.StoreError as e:
708 die(e)
711def command_export(args):
712 from subprocess import Popen, PIPE
714 try:
715 from tunguska import gfdb
716 except ImportError as err:
717 die('the kiwi tools must be installed to use this feature', err)
719 def setup(parser):
720 parser.add_option(
721 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
722 help='split output gfdb into N chunks')
724 parser, options, args = cl_parse('export', args, setup=setup)
726 show_progress = True
728 if len(args) not in (1, 2):
729 sys.exit(parser.format_help())
731 target_path = args.pop()
732 if op.isdir(target_path):
733 target_path = op.join(target_path, 'kiwi_gfdb')
734 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
736 source_store_dir = get_store_dir(args)
738 source = gf.Store(source_store_dir, 'r')
739 config = source.config
741 if not isinstance(config, gf.meta.ConfigTypeA):
742 die('only stores of type A can be exported to Kiwi format')
744 if op.isfile(target_path + '.index'):
745 die('destation already exists')
747 cmd = [str(x) for x in [
748 'gfdb_build',
749 target_path,
750 options.nchunks,
751 config.ndistances,
752 config.nsource_depths,
753 config.ncomponents,
754 config.deltat,
755 config.distance_delta,
756 config.source_depth_delta,
757 config.distance_min,
758 config.source_depth_min]]
760 p = Popen(cmd, stdin=PIPE)
761 p.communicate()
763 out_db = gfdb.Gfdb(target_path)
765 try:
766 if show_progress:
767 pbar = util.progressbar(
768 'exporting', config.nrecords/config.ncomponents)
770 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
772 data_out = []
773 for ig in range(config.ncomponents):
774 try:
775 tr = source.get((z, x, ig), interpolation='off')
776 data_out.append((tr.t, tr.data * config.factor))
778 except gf.store.StoreError as e:
779 logger.warning(
780 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
781 data_out.append(None)
783 # put a zero valued sample to no-data zero-traces at a compatible
784 # time
785 tmins = [
786 entry[0][0]
787 for entry in data_out
788 if entry is not None and entry[0].size != 0]
790 if tmins:
791 tmin = min(tmins)
792 for entry in data_out:
793 if entry is not None and entry[0].size == 0:
794 entry[0].resize(1)
795 entry[1].resize(1)
796 entry[0][0] = tmin
797 entry[1][0] = 0.0
799 out_db.put_traces_slow(x, z, data_out)
801 if show_progress:
802 pbar.update(i+1)
804 source.close()
806 finally:
807 if show_progress:
808 pbar.finish()
811def phasedef_or_horvel(x):
812 try:
813 return float(x)
814 except ValueError:
815 return cake.PhaseDef(x)
818def mkp(s):
819 return [phasedef_or_horvel(ps) for ps in s.split(',')]
822def stored_attribute_table_plots(phase_ids, options, args, attribute):
823 import numpy as num
824 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
826 plt = mpl_init()
828 np = 1
829 store_dir = get_store_dir(args)
830 for phase_id in phase_ids:
831 try:
832 store = gf.Store(store_dir)
834 if options.receiver_depth is not None:
835 receiver_depth = options.receiver_depth * 1000.0
836 else:
837 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
838 receiver_depth = store.config.receiver_depth
840 elif isinstance(store.config, gf.ConfigTypeB):
841 receiver_depth = store.config.receiver_depth_min
843 else:
844 receiver_depth = 0.0
846 phase = store.get_stored_phase(phase_id, attribute)
847 axes = plt.subplot(2, len(phase_ids), np)
848 labelspace(axes)
849 xscaled(1. / km, axes)
850 yscaled(1. / km, axes)
851 x = None
852 if isinstance(store.config, gf.ConfigTypeB):
853 x = (receiver_depth, None, None)
855 phase.plot_2d(axes, x=x)
856 axes.set_title(phase_id)
857 np += 1
858 except gf.StoreError as e:
859 die(e)
861 axes = plt.subplot(2, 1, 2)
862 num_d = 100
863 distances = num.linspace(store.config.distance_min,
864 store.config.distance_max,
865 num_d)
867 if options.source_depth is not None:
868 source_depth = options.source_depth * km
869 else:
870 source_depth = store.config.source_depth_min + (
871 store.config.source_depth_max - store.config.source_depth_min) / 2.
873 if isinstance(store.config, gf.ConfigTypeA):
874 attribute_vals = num.empty(num_d)
875 for phase_id in phase_ids:
876 attribute_vals[:] = num.NAN
877 for i, d in enumerate(distances):
878 if attribute == 'phase':
879 attribute_vals[i] = store.t(phase_id, (source_depth, d))
880 ylabel = 'TT [s]'
881 else:
882 attribute_vals[i] = store.get_stored_attribute(
883 phase_id, options.attribute, (source_depth, d))
884 ylabel = '%s [deg]' % options.attribute
886 axes.plot(distances / km, attribute_vals, label=phase_id)
888 axes.set_title('source source_depth %s km' % (source_depth / km))
889 axes.set_xlabel('distance [km]')
890 axes.set_ylabel(ylabel)
891 axes.legend()
893 plt.tight_layout()
894 mpl_show(plt)
897def command_ttt(args):
898 def setup(parser):
899 parser.add_option(
900 '--force', dest='force', action='store_true',
901 help='overwrite existing files')
903 parser, options, args = cl_parse('ttt', args, setup=setup)
905 store_dir = get_store_dir(args)
906 try:
907 store = gf.Store(store_dir)
908 store.make_travel_time_tables(force=options.force)
910 except gf.StoreError as e:
911 die(e)
914def command_tttview(args):
916 def setup(parser):
917 parser.add_option(
918 '--source-depth', dest='source_depth', type=float,
919 help='Source depth in km')
921 parser.add_option(
922 '--receiver-depth', dest='receiver_depth', type=float,
923 help='Receiver depth in km')
925 parser, options, args = cl_parse(
926 'tttview', args, setup=setup,
927 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
929 try:
930 phase_ids = args.pop().split(',')
931 except Exception:
932 parser.error('cannot get <phase-ids> argument')
934 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
937def command_sat(args):
938 def setup(parser):
939 parser.add_option(
940 '--force', dest='force', action='store_true',
941 help='overwrite existing files')
943 parser.add_option(
944 '--attribute',
945 action='store',
946 dest='attribute',
947 type='choice',
948 choices=gf.store.available_stored_tables[1::],
949 default='takeoff_angle',
950 help='calculate interpolation table for selected ray attributes.')
952 parser, options, args = cl_parse('sat', args, setup=setup)
954 store_dir = get_store_dir(args)
955 try:
956 store = gf.Store(store_dir)
957 store.make_stored_table(options.attribute, force=options.force)
959 except gf.StoreError as e:
960 die(e)
963def command_satview(args):
965 def setup(parser):
966 parser.add_option(
967 '--source-depth', dest='source_depth', type=float,
968 help='Source depth in km')
970 parser.add_option(
971 '--receiver-depth', dest='receiver_depth', type=float,
972 help='Receiver depth in km')
974 parser.add_option(
975 '--attribute',
976 action='store',
977 dest='attribute',
978 type='choice',
979 choices=gf.store.available_stored_tables[1::],
980 default='takeoff_angle',
981 help='view selected ray attribute.')
983 parser, options, args = cl_parse(
984 'satview', args, setup=setup,
985 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
987 try:
988 phase_ids = args.pop().split(',')
989 except Exception:
990 parser.error('cannot get <phase-ids> argument')
992 logger.info('Plotting stored attribute %s' % options.attribute)
994 stored_attribute_table_plots(
995 phase_ids, options, args, attribute=options.attribute)
998def command_tttextract(args):
999 def setup(parser):
1000 parser.add_option(
1001 '--output', dest='output_fn', metavar='TEMPLATE',
1002 help='output to text files instead of stdout '
1003 '(example TEMPLATE: "extracted/%(args)s.txt")')
1005 parser, options, args = cl_parse('tttextract', args, setup=setup)
1006 try:
1007 sdef = args.pop()
1008 except Exception:
1009 parser.error('cannot get <selection> argument')
1011 try:
1012 sphase = args.pop()
1013 except Exception:
1014 parser.error('cannot get <phase> argument')
1016 try:
1017 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1018 except gf.meta.InvalidTimingSpecification:
1019 parser.error('invalid phase specification: "%s"' % sphase)
1021 try:
1022 gdef = gf.meta.parse_grid_spec(sdef)
1023 except gf.meta.GridSpecError as e:
1024 die(e)
1026 store_dir = get_store_dir(args)
1028 try:
1029 store = gf.Store(store_dir)
1030 for args in store.config.iter_extraction(gdef, level=-1):
1031 s = ['%e' % x for x in args]
1032 for phase in phases:
1033 t = store.t(phase, args)
1034 if t is not None:
1035 s.append('%e' % t)
1036 else:
1037 s.append('nan')
1039 if options.output_fn:
1040 d = dict(
1041 args='_'.join('%e' % x for x in args),
1042 extension='txt')
1044 fn = options.output_fn % d
1045 util.ensuredirs(fn)
1046 with open(fn, 'a') as f:
1047 f.write(' '.join(s))
1048 f.write('\n')
1049 else:
1050 print(' '.join(s))
1052 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1053 die(e)
1056def command_tttlsd(args):
1058 def setup(parser):
1059 pass
1061 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1063 try:
1064 sphase_ids = args.pop()
1065 except Exception:
1066 parser.error('cannot get <phase> argument')
1068 try:
1069 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1070 except gf.meta.InvalidTimingSpecification:
1071 parser.error('invalid phase specification: "%s"' % sphase_ids)
1073 store_dir = get_store_dir(args)
1075 try:
1076 store = gf.Store(store_dir)
1077 for phase_id in phase_ids:
1078 store.fix_ttt_holes(phase_id)
1080 except gf.StoreError as e:
1081 die(e)
1084def command_server(args):
1085 from pyrocko.gf import server
1087 def setup(parser):
1088 parser.add_option(
1089 '--port', dest='port', metavar='PORT', type='int', default=8080,
1090 help='serve on port PORT')
1092 parser.add_option(
1093 '--ip', dest='ip', metavar='IP', default='',
1094 help='serve on ip address IP')
1096 parser, options, args = cl_parse('server', args, setup=setup)
1098 engine = gf.LocalEngine(store_superdirs=args)
1099 server.run(options.ip, options.port, engine)
1102def command_download(args):
1103 from pyrocko.gf import ws
1105 details = '''
1107Browse pre-calculated Green's function stores online:
1109 https://greens-mill.pyrocko.org
1110'''
1112 def setup(parser):
1113 parser.add_option(
1114 '--force', dest='force', action='store_true',
1115 help='overwrite existing files')
1117 parser, options, args = cl_parse(
1118 'download', args, setup=setup, details=details)
1119 if not len(args) in (1, 2):
1120 sys.exit(parser.format_help())
1122 if len(args) == 2:
1123 site, store_id = args
1124 if not re.match(gf.meta.StringID.pattern, store_id):
1125 die('invalid store ID')
1126 else:
1127 site, store_id = args[0], None
1129 if site not in gf.ws.g_site_abbr:
1130 if -1 == site.find('://'):
1131 site = 'http://' + site
1133 try:
1134 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1135 except ws.DownloadError as e:
1136 die(str(e))
1139def command_modelview(args):
1141 import matplotlib.pyplot as plt
1142 import numpy as num
1143 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1144 mpl_init()
1146 neat_labels = {
1147 'vp': '$v_p$',
1148 'vs': '$v_s$',
1149 'qp': '$Q_p$',
1150 'qs': '$Q_s$',
1151 'rho': '$\\rho$'}
1153 def setup(parser):
1154 parser.add_option(
1155 '--parameters', dest='parameters',
1156 default='vp,vs', metavar='vp,vs,....',
1157 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1159 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1161 parser, options, args = cl_parse('modelview', args, setup=setup)
1163 store_dirs = get_store_dirs(args)
1165 parameters = options.parameters.split(',')
1167 fig, axes = plt.subplots(1,
1168 len(parameters),
1169 sharey=True,
1170 figsize=(len(parameters)*3, 5))
1172 if not isinstance(axes, num.ndarray):
1173 axes = [axes]
1175 axes = dict(zip(parameters, axes))
1177 for store_id in store_dirs:
1178 try:
1179 store = gf.Store(store_id)
1180 mod = store.config.earthmodel_1d
1181 z = mod.profile('z')
1183 for p in parameters:
1184 ax = axes[p]
1185 labelspace(ax)
1186 if '/' in p:
1187 p1, p2 = p.split('/')
1188 profile = mod.profile(p1)/mod.profile(p2)
1189 else:
1190 profile = mod.profile(p)
1192 ax.plot(profile, -z, label=store_id.split('/')[-1])
1193 if p in ['vs', 'vp', 'rho']:
1194 xscaled(1./1000., ax)
1196 yscaled(1./km, ax)
1198 except gf.StoreError as e:
1199 die(e)
1201 for p, ax in axes.items():
1202 ax.grid()
1203 if p in neat_labels:
1204 lab = neat_labels[p]
1205 elif all(x in neat_labels for x in p.split('/')):
1206 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1207 else:
1208 lab = p
1210 ax.set_title(lab, y=1.02)
1212 if p in units:
1213 lab += ' ' + units[p]
1215 ax.autoscale()
1216 ax.set_xlabel(lab)
1218 axes[parameters[0]].set_ylabel('Depth [km]')
1220 handles, labels = ax.get_legend_handles_labels()
1222 if len(store_dirs) > 1:
1223 ax.legend(
1224 handles,
1225 labels,
1226 bbox_to_anchor=(0.5, 0.12),
1227 bbox_transform=fig.transFigure,
1228 ncol=3,
1229 loc='upper center',
1230 fancybox=True)
1232 plt.subplots_adjust(bottom=0.22,
1233 wspace=0.05)
1235 mpl_show(plt)
1238def command_upgrade(args):
1239 parser, options, args = cl_parse('upgrade', args)
1240 store_dirs = get_store_dirs(args)
1241 try:
1242 for store_dir in store_dirs:
1243 store = gf.Store(store_dir)
1244 nup = store.upgrade()
1245 if nup == 0:
1246 print('%s: already up-to-date.' % store_dir)
1247 else:
1248 print('%s: %i file%s upgraded.' % (
1249 store_dir, nup, ['s', ''][nup == 1]))
1251 except gf.StoreError as e:
1252 die(e)
1255def command_addref(args):
1256 parser, options, args = cl_parse('addref', args)
1258 store_dirs = []
1259 filenames = []
1260 for arg in args:
1261 if op.isdir(arg):
1262 store_dirs.append(arg)
1263 elif op.isfile(arg):
1264 filenames.append(arg)
1265 else:
1266 die('invalid argument: %s' % arg)
1268 if not store_dirs:
1269 store_dirs.append('.')
1271 references = []
1272 try:
1273 for filename in args:
1274 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1275 except ImportError:
1276 die('pybtex module must be installed to use this function')
1278 if not references:
1279 die('no references found')
1281 for store_dir in store_dirs:
1282 try:
1283 store = gf.Store(store_dir)
1284 ids = [ref.id for ref in store.config.references]
1285 for ref in references:
1286 if ref.id in ids:
1287 die('duplicate reference id: %s' % ref.id)
1289 ids.append(ref.id)
1290 store.config.references.append(ref)
1292 store.save_config(make_backup=True)
1294 except gf.StoreError as e:
1295 die(e)
1298def command_qc(args):
1299 parser, options, args = cl_parse('qc', args)
1301 store_dir = get_store_dir(args)
1303 try:
1304 store = gf.Store(store_dir)
1305 s = store.stats()
1306 if s['empty'] != 0:
1307 print('has empty records')
1309 for aname in ['author', 'author_email', 'description', 'references']:
1311 if not getattr(store.config, aname):
1312 print('%s empty' % aname)
1314 except gf.StoreError as e:
1315 die(e)
1318def command_report(args):
1319 from pyrocko.fomosto.report import report_call
1320 report_call.run_program(args)
1323def main(args=None):
1324 if args is None:
1325 args = sys.argv[1:]
1327 if len(args) < 1:
1328 sys.exit('Usage: %s' % usage)
1330 command = args.pop(0)
1332 if command in subcommands:
1333 globals()['command_' + command](args)
1335 elif command in ('--help', '-h', 'help'):
1336 if command == 'help' and args:
1337 acommand = args[0]
1338 if acommand in subcommands:
1339 globals()['command_' + acommand](['--help'])
1341 sys.exit('Usage: %s' % usage)
1343 else:
1344 sys.exit('fomosto: error: no such subcommand: %s' % command)
1347if __name__ == '__main__':
1348 main(sys.argv[1:])