1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import sys
7import re
8import os.path as op
9import logging
10import copy
11import shutil
12from optparse import OptionParser
14from pyrocko import util, trace, gf, cake, io, fomosto
15from pyrocko.gui.snuffler import marker
16from pyrocko.util import mpl_show
18logger = logging.getLogger('pyrocko.apps.fomosto')
19km = 1e3
22def d2u(d):
23 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
26subcommand_descriptions = {
27 'init': 'create a new empty GF store',
28 'build': 'compute GFs and fill into store',
29 'stats': 'print information about a GF store',
30 'check': 'check for problems in GF store',
31 'decimate': 'build decimated variant of a GF store',
32 'redeploy': 'copy traces from one GF store into another',
33 'view': 'view selected traces',
34 'extract': 'extract selected traces',
35 'import': 'convert Kiwi GFDB to GF store format',
36 'export': 'convert GF store to Kiwi GFDB format',
37 'ttt': 'create travel time tables',
38 'tttview': 'plot travel time table',
39 'tttextract': 'extract selected travel times',
40 'tttlsd': 'fix holes in travel time tables',
41 'sat': 'create stored ray attribute table',
42 'satview': 'plot stored ray attribute table',
43 'server': 'run seismosizer server',
44 'download': 'download GF store from a server',
45 'modelview': 'plot earthmodels',
46 'upgrade': 'upgrade store format to latest version',
47 'addref': 'import citation references to GF store config',
48 'qc': 'quality check',
49 'report': "report for Green's Function databases",
50}
52subcommand_usages = {
53 'init': ['init <type> <store-dir> [options]',
54 'init redeploy <source> <destination> [options]'],
55 'build': 'build [store-dir] [options]',
56 'stats': 'stats [store-dir] [options]',
57 'check': 'check [store-dir] [options]',
58 'decimate': 'decimate [store-dir] <factor> [options]',
59 'redeploy': 'redeploy <source> <destination> [options]',
60 'view': 'view [store-dir] ... [options]',
61 'extract': 'extract [store-dir] <selection>',
62 'import': 'import <source> <destination> [options]',
63 'export': 'export [store-dir] <destination> [options]',
64 'ttt': 'ttt [store-dir] [options]',
65 'tttview': 'tttview [store-dir] <phase-ids> [options]',
66 'tttextract': 'tttextract [store-dir] <phase> <selection>',
67 'tttlsd': 'tttlsd [store-dir] <phase>',
68 'sat': 'sat [store-dir] [options]',
69 'satview': 'satview [store-dir] <phase-ids> [options]',
70 'server': 'server [options] <store-super-dir> ...',
71 'download': 'download [options] <site> <store-id>',
72 'modelview': 'modelview <selection>',
73 'upgrade': 'upgrade [store-dir] ...',
74 'addref': 'addref [store-dir] ... <filename> ...',
75 'qc': 'qc [store-dir]',
76 'report': 'report <subcommand> <arguments>... [options]'
77}
79subcommands = subcommand_descriptions.keys()
81program_name = 'fomosto'
83usage = program_name + ''' <subcommand> <arguments> ... [options]
85Subcommands:
87 init %(init)s
88 build %(build)s
89 stats %(stats)s
90 check %(check)s
91 decimate %(decimate)s
92 redeploy %(redeploy)s
93 view %(view)s
94 extract %(extract)s
95 import %(import)s
96 export %(export)s
97 ttt %(ttt)s
98 tttview %(tttview)s
99 tttextract %(tttextract)s
100 tttlsd %(tttlsd)s
101 sat %(sat)s
102 satview %(satview)s
103 server %(server)s
104 download %(download)s
105 modelview %(modelview)s
106 upgrade %(upgrade)s
107 addref %(addref)s
108 qc %(qc)s
109 report %(report)s
111To get further help and a list of available options for any subcommand run:
113 fomosto <subcommand> --help
115''' % d2u(subcommand_descriptions)
118def add_common_options(parser):
119 parser.add_option(
120 '--loglevel',
121 action='store',
122 dest='loglevel',
123 type='choice',
124 choices=('critical', 'error', 'warning', 'info', 'debug'),
125 default='info',
126 help='set logger level to '
127 '"critical", "error", "warning", "info", or "debug". '
128 'Default is "%default".')
131def process_common_options(options):
132 util.setup_logging(program_name, options.loglevel)
135def cl_parse(command, args, setup=None, details=None):
136 usage = subcommand_usages[command]
137 descr = subcommand_descriptions[command]
139 if isinstance(usage, str):
140 usage = [usage]
142 susage = '%s %s' % (program_name, usage[0])
143 for s in usage[1:]:
144 susage += '\n%s%s %s' % (' '*7, program_name, s)
146 description = descr[0].upper() + descr[1:] + '.'
148 if details:
149 description = description + ' %s' % details
151 parser = OptionParser(usage=susage, description=description)
152 parser.format_description = lambda formatter: description
154 if setup:
155 setup(parser)
157 add_common_options(parser)
158 (options, args) = parser.parse_args(args)
159 process_common_options(options)
160 return parser, options, args
163def die(message, err=''):
164 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
167def fomo_wrapper_module(name):
168 try:
169 if not re.match(gf.meta.StringID.pattern, name):
170 raise ValueError('invalid name')
172 words = name.split('.', 1)
173 if len(words) == 2:
174 name, variant = words
175 else:
176 name = words[0]
177 variant = None
179 name_clean = re.sub(r'[.-]', '_', name)
180 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
181 mod = __import__(modname, level=0)
182 return getattr(mod.fomosto, name_clean), variant
184 except ValueError:
185 die('invalid modelling code wrapper name')
187 except ImportError:
188 die('''modelling code wrapper "%s" not available or not installed
189 (module probed: "%s")''' % (name, modname))
192def command_init(args):
194 details = '''
196Available modelling backends:
197%s
199 More information at
200 https://pyrocko.org/docs/current/apps/fomosto/backends.html
201''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
203 parser, options, args = cl_parse(
204 'init', args,
205 details=details)
207 if len(args) == 0:
208 sys.exit(parser.format_help())
210 if args[0] == 'redeploy':
211 if len(args) != 3:
212 parser.error('incorrect number of arguments')
214 source_dir, dest_dir = args[1:]
216 try:
217 source = gf.Store(source_dir)
218 except gf.StoreError as e:
219 die(e)
221 config = copy.deepcopy(source.config)
222 config.derived_from_id = source.config.id
223 try:
224 config_filenames = gf.store.Store.create_editables(
225 dest_dir, config=config)
227 except gf.StoreError as e:
228 die(e)
230 try:
231 dest = gf.Store(dest_dir)
232 except gf.StoreError as e:
233 die(e)
235 for k in source.extra_keys():
236 source_fn = source.get_extra_path(k)
237 dest_fn = dest.get_extra_path(k)
238 shutil.copyfile(source_fn, dest_fn)
240 logger.info(
241 '(1) configure settings in files:\n %s'
242 % '\n '.join(config_filenames))
244 logger.info(
245 '(2) run "fomosto redeploy <source> <dest>", as needed')
247 else:
248 if len(args) != 2:
249 parser.error('incorrect number of arguments')
251 (modelling_code_id, store_dir) = args
253 module, variant = fomo_wrapper_module(modelling_code_id)
254 try:
255 config_filenames = module.init(store_dir, variant)
256 except gf.StoreError as e:
257 die(e)
259 logger.info('(1) configure settings in files:\n %s'
260 % '\n '.join(config_filenames))
261 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
262 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
265def get_store_dir(args):
266 if len(args) == 1:
267 store_dir = op.abspath(args.pop(0))
268 else:
269 store_dir = op.abspath(op.curdir)
271 if not op.isdir(store_dir):
272 die('not a directory: %s' % store_dir)
274 return store_dir
277def get_store_dirs(args):
278 if len(args) == 0:
279 store_dirs = [op.abspath(op.curdir)]
280 else:
281 store_dirs = [op.abspath(x) for x in args]
283 for store_dir in store_dirs:
284 if not op.isdir(store_dir):
285 die('not a directory: %s' % store_dir)
287 return store_dirs
290def command_build(args):
292 def setup(parser):
293 parser.add_option(
294 '--force', dest='force', action='store_true',
295 help='overwrite existing files')
297 parser.add_option(
298 '--nworkers', dest='nworkers', type='int', metavar='N',
299 help='run N worker processes in parallel')
301 parser.add_option(
302 '--continue', dest='continue_', action='store_true',
303 help='continue suspended build')
305 parser.add_option(
306 '--step', dest='step', type='int', metavar='I',
307 help='process block number IBLOCK')
309 parser.add_option(
310 '--block', dest='iblock', type='int', metavar='I',
311 help='process block number IBLOCK')
313 parser, options, args = cl_parse('build', args, setup=setup)
315 store_dir = get_store_dir(args)
316 try:
317 if options.step is not None:
318 step = options.step - 1
319 else:
320 step = None
322 if options.iblock is not None:
323 iblock = options.iblock - 1
324 else:
325 iblock = None
327 store = gf.Store(store_dir)
328 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
329 module.build(
330 store_dir,
331 force=options.force,
332 nworkers=options.nworkers, continue_=options.continue_,
333 step=step,
334 iblock=iblock)
336 except gf.StoreError as e:
337 die(e)
340def command_stats(args):
342 parser, options, args = cl_parse('stats', args)
343 store_dir = get_store_dir(args)
345 try:
346 store = gf.Store(store_dir)
347 s = store.stats()
349 except gf.StoreError as e:
350 die(e)
352 for k in store.stats_keys:
353 print('%s: %s' % (k, s[k]))
356def command_check(args):
358 parser, options, args = cl_parse('check', args)
359 store_dir = get_store_dir(args)
361 try:
362 store = gf.Store(store_dir)
363 problems = store.check(show_progress=True)
364 if problems:
365 die('problems detected with gf store: %s' % store_dir)
367 except gf.StoreError as e:
368 die(e)
371def load_config(fn):
372 try:
373 config = gf.meta.load(filename=fn)
374 assert isinstance(config, gf.Config)
376 except Exception:
377 die('cannot load gf config from file: %s' % fn)
379 return config
382def command_decimate(args):
384 def setup(parser):
385 parser.add_option(
386 '--config', dest='config_fn', metavar='FILE',
387 help='use modified spacial sampling given in FILE')
389 parser.add_option(
390 '--force', dest='force', action='store_true',
391 help='overwrite existing files')
393 parser, options, args = cl_parse('decimate', args, setup=setup)
394 try:
395 decimate = int(args.pop())
396 except Exception:
397 parser.error('cannot get <factor> argument')
399 store_dir = get_store_dir(args)
401 config = None
402 if options.config_fn:
403 config = load_config(options.config_fn)
405 try:
406 store = gf.Store(store_dir)
407 store.make_decimated(decimate, config=config, force=options.force,
408 show_progress=True)
410 except gf.StoreError as e:
411 die(e)
414def sindex(args):
415 return '(%s)' % ', '.join('%g' % x for x in args)
418def command_redeploy(args):
420 parser, options, args = cl_parse('redeploy', args)
422 if not len(args) == 2:
423 sys.exit(parser.format_help())
425 source_store_dir, dest_store_dir = args
427 try:
428 source = gf.Store(source_store_dir)
429 except gf.StoreError as e:
430 die(e)
432 try:
433 gf.store.Store.create_dependants(dest_store_dir)
434 except gf.StoreError:
435 pass
437 try:
438 dest = gf.Store(dest_store_dir, 'w')
440 except gf.StoreError as e:
441 die(e)
443 show_progress = True
445 try:
446 if show_progress:
447 pbar = util.progressbar('redeploying', dest.config.nrecords)
449 for i, args in enumerate(dest.config.iter_nodes()):
450 try:
451 tr = source.get(args, interpolation='off')
452 dest.put(args, tr)
454 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa
455 logger.debug('skipping %s, (%s)' % (sindex(args), e))
457 except gf.store.StoreError as e:
458 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
460 if show_progress:
461 pbar.update(i+1)
463 finally:
464 if show_progress:
465 pbar.finish()
468def command_view(args):
469 def setup(parser):
470 parser.add_option(
471 '--extract',
472 dest='extract',
473 metavar='start:stop[:step|@num],...',
474 help='specify which traces to show')
476 parser.add_option(
477 '--show-phases',
478 dest='showphases',
479 default=None,
480 metavar='[phase_id1,phase_id2,...|all]',
481 help='add phase markers from ttt')
483 parser.add_option(
484 '--opengl',
485 dest='opengl',
486 action='store_true',
487 default=None,
488 help='use OpenGL for drawing')
490 parser.add_option(
491 '--no-opengl',
492 dest='opengl',
493 action='store_false',
494 default=None,
495 help='do not use OpenGL for drawing')
497 parser, options, args = cl_parse('view', args, setup=setup)
499 gdef = None
500 if options.extract:
501 try:
502 gdef = gf.meta.parse_grid_spec(options.extract)
503 except gf.meta.GridSpecError as e:
504 die(e)
506 store_dirs = get_store_dirs(args)
508 alpha = 'abcdefghijklmnopqrstxyz'.upper()
510 markers = []
511 traces = []
513 try:
514 for istore, store_dir in enumerate(store_dirs):
515 store = gf.Store(store_dir)
516 if options.showphases == 'all':
517 phasenames = [pn.id for pn in store.config.tabulated_phases]
518 elif options.showphases is not None:
519 phasenames = options.showphases.split(',')
521 ii = 0
522 for args in store.config.iter_extraction(gdef):
523 gtr = store.get(args)
525 loc_code = ''
526 if len(store_dirs) > 1:
527 loc_code = alpha[istore % len(alpha)]
529 if gtr:
531 sta_code = '%04i (%s)' % (
532 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
533 tmin = gtr.deltat * gtr.itmin
535 tr = trace.Trace(
536 '',
537 sta_code,
538 loc_code,
539 '%02i' % args[-1],
540 ydata=gtr.data,
541 deltat=gtr.deltat,
542 tmin=tmin)
544 if options.showphases:
545 for phasename in phasenames:
546 phase_tmin = store.t(phasename, args[:-1])
547 if phase_tmin:
548 m = marker.PhaseMarker(
549 [('',
550 sta_code,
551 loc_code,
552 '%02i' % args[-1])],
553 phase_tmin,
554 phase_tmin,
555 0,
556 phasename=phasename)
557 markers.append(m)
559 traces.append(tr)
561 ii += 1
563 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
564 die(e)
566 def setup(win):
567 viewer = win.get_view()
568 viewer.menuitem_showboxes.setChecked(False)
569 viewer.menuitem_colortraces.setChecked(True)
570 viewer.menuitem_demean.setChecked(False)
571 viewer.menuitems_sorting[5][0].setChecked(True)
572 viewer.menuitems_scaling[2][0].setChecked(True)
573 viewer.sortingmode_change()
574 viewer.scalingmode_change()
576 trace.snuffle(
577 traces, markers=markers, opengl=options.opengl, launch_hook=setup)
580def command_extract(args):
581 def setup(parser):
582 parser.add_option(
583 '--format', dest='format', default='mseed',
584 choices=['mseed', 'sac', 'text', 'yaff'],
585 help='export to format "mseed", "sac", "text", or "yaff". '
586 'Default is "mseed".')
588 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
589 parser.add_option(
590 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
591 help='output path template [default: "%s"]' % fndfl)
593 parser, options, args = cl_parse('extract', args, setup=setup)
594 try:
595 sdef = args.pop()
596 except Exception:
597 parser.error('cannot get <selection> argument')
599 try:
600 gdef = gf.meta.parse_grid_spec(sdef)
601 except gf.meta.GridSpecError as e:
602 die(e)
604 store_dir = get_store_dir(args)
606 extensions = {
607 'mseed': 'mseed',
608 'sac': 'sac',
609 'text': 'txt',
610 'yaff': 'yaff'}
612 try:
613 store = gf.Store(store_dir)
614 for args in store.config.iter_extraction(gdef):
615 gtr = store.get(args)
616 if gtr:
617 tr = trace.Trace(
618 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
619 ydata=gtr.data,
620 deltat=gtr.deltat,
621 tmin=gtr.deltat * gtr.itmin)
623 additional = dict(
624 args='_'.join('%g' % x for x in args),
625 irecord=store.str_irecord(args),
626 extension=extensions[options.format])
628 io.save(
629 tr,
630 options.output_fn,
631 format=options.format,
632 additional=additional)
634 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
635 die(e)
638def command_import(args):
639 try:
640 from tunguska import gfdb
641 except ImportError:
642 die('the kiwi tools must be installed to use this feature')
644 parser, options, args = cl_parse('import', args)
646 show_progress = True
648 if not len(args) == 2:
649 sys.exit(parser.format_help())
651 source_path, dest_store_dir = args
653 if op.isdir(source_path):
654 source_path = op.join(source_path, 'db')
656 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
658 db = gfdb.Gfdb(source_path)
660 config = gf.meta.ConfigTypeA(
661 id='imported_gfs',
662 distance_min=db.firstx,
663 distance_max=db.firstx + (db.nx-1) * db.dx,
664 distance_delta=db.dx,
665 source_depth_min=db.firstz,
666 source_depth_max=db.firstz + (db.nz-1) * db.dz,
667 source_depth_delta=db.dz,
668 sample_rate=1.0/db.dt,
669 ncomponents=db.ng
670 )
672 try:
673 gf.store.Store.create(dest_store_dir, config=config)
674 dest = gf.Store(dest_store_dir, 'w')
675 try:
676 if show_progress:
677 pbar = util.progressbar(
678 'importing', dest.config.nrecords/dest.config.ncomponents)
680 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
681 source_depth, distance = [float(x) for x in args]
682 traces = db.get_traces_pyrocko(distance, source_depth)
683 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
684 for ig in range(db.ng):
685 if ig in ig_to_trace:
686 tr = ig_to_trace[ig]
687 gf_tr = gf.store.GFTrace(
688 tr.get_ydata(),
689 int(round(tr.tmin / tr.deltat)),
690 tr.deltat)
692 else:
693 gf_tr = gf.store.Zero
695 dest.put((source_depth, distance, ig), gf_tr)
697 if show_progress:
698 pbar.update(i+1)
700 finally:
701 if show_progress:
702 pbar.finish()
704 dest.close()
706 except gf.StoreError as e:
707 die(e)
710def command_export(args):
711 from subprocess import Popen, PIPE
713 try:
714 from tunguska import gfdb
715 except ImportError as err:
716 die('the kiwi tools must be installed to use this feature', err)
718 def setup(parser):
719 parser.add_option(
720 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
721 help='split output gfdb into N chunks')
723 parser, options, args = cl_parse('export', args, setup=setup)
725 show_progress = True
727 if len(args) not in (1, 2):
728 sys.exit(parser.format_help())
730 target_path = args.pop()
731 if op.isdir(target_path):
732 target_path = op.join(target_path, 'kiwi_gfdb')
733 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
735 source_store_dir = get_store_dir(args)
737 source = gf.Store(source_store_dir, 'r')
738 config = source.config
740 if not isinstance(config, gf.meta.ConfigTypeA):
741 die('only stores of type A can be exported to Kiwi format')
743 if op.isfile(target_path + '.index'):
744 die('destation already exists')
746 cmd = [str(x) for x in [
747 'gfdb_build',
748 target_path,
749 options.nchunks,
750 config.ndistances,
751 config.nsource_depths,
752 config.ncomponents,
753 config.deltat,
754 config.distance_delta,
755 config.source_depth_delta,
756 config.distance_min,
757 config.source_depth_min]]
759 p = Popen(cmd, stdin=PIPE)
760 p.communicate()
762 out_db = gfdb.Gfdb(target_path)
764 try:
765 if show_progress:
766 pbar = util.progressbar(
767 'exporting', config.nrecords/config.ncomponents)
769 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
771 data_out = []
772 for ig in range(config.ncomponents):
773 try:
774 tr = source.get((z, x, ig), interpolation='off')
775 data_out.append((tr.t, tr.data * config.factor))
777 except gf.store.StoreError as e:
778 logger.warning(
779 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
780 data_out.append(None)
782 # put a zero valued sample to no-data zero-traces at a compatible
783 # time
784 tmins = [
785 entry[0][0]
786 for entry in data_out
787 if entry is not None and entry[0].size != 0]
789 if tmins:
790 tmin = min(tmins)
791 for entry in data_out:
792 if entry is not None and entry[0].size == 0:
793 entry[0].resize(1)
794 entry[1].resize(1)
795 entry[0][0] = tmin
796 entry[1][0] = 0.0
798 out_db.put_traces_slow(x, z, data_out)
800 if show_progress:
801 pbar.update(i+1)
803 source.close()
805 finally:
806 if show_progress:
807 pbar.finish()
810def phasedef_or_horvel(x):
811 try:
812 return float(x)
813 except ValueError:
814 return cake.PhaseDef(x)
817def mkp(s):
818 return [phasedef_or_horvel(ps) for ps in s.split(',')]
821def stored_attribute_table_plots(phase_ids, options, args, attribute):
822 import numpy as num
823 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
825 plt = mpl_init()
827 np = 1
828 store_dir = get_store_dir(args)
829 for phase_id in phase_ids:
830 try:
831 store = gf.Store(store_dir)
833 if options.receiver_depth is not None:
834 receiver_depth = options.receiver_depth * 1000.0
835 else:
836 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
837 receiver_depth = store.config.receiver_depth
839 elif isinstance(store.config, gf.ConfigTypeB):
840 receiver_depth = store.config.receiver_depth_min
842 else:
843 receiver_depth = 0.0
845 phase = store.get_stored_phase(phase_id, attribute)
846 axes = plt.subplot(2, len(phase_ids), np)
847 labelspace(axes)
848 xscaled(1. / km, axes)
849 yscaled(1. / km, axes)
850 x = None
851 if isinstance(store.config, gf.ConfigTypeB):
852 x = (receiver_depth, None, None)
854 phase.plot_2d(axes, x=x)
855 axes.set_title(phase_id)
856 np += 1
857 except gf.StoreError as e:
858 die(e)
860 axes = plt.subplot(2, 1, 2)
861 num_d = 100
862 distances = num.linspace(store.config.distance_min,
863 store.config.distance_max,
864 num_d)
866 if options.source_depth is not None:
867 source_depth = options.source_depth * km
868 else:
869 source_depth = store.config.source_depth_min + (
870 store.config.source_depth_max - store.config.source_depth_min) / 2.
872 if isinstance(store.config, gf.ConfigTypeA):
873 attribute_vals = num.empty(num_d)
874 for phase_id in phase_ids:
875 attribute_vals[:] = num.NAN
876 for i, d in enumerate(distances):
877 if attribute == 'phase':
878 attribute_vals[i] = store.t(phase_id, (source_depth, d))
879 ylabel = 'TT [s]'
880 else:
881 attribute_vals[i] = store.get_stored_attribute(
882 phase_id, options.attribute, (source_depth, d))
883 ylabel = '%s [deg]' % options.attribute
885 axes.plot(distances / km, attribute_vals, label=phase_id)
887 axes.set_title('source source_depth %s km' % (source_depth / km))
888 axes.set_xlabel('distance [km]')
889 axes.set_ylabel(ylabel)
890 axes.legend()
892 plt.tight_layout()
893 mpl_show(plt)
896def command_ttt(args):
897 def setup(parser):
898 parser.add_option(
899 '--force', dest='force', action='store_true',
900 help='overwrite existing files')
902 parser, options, args = cl_parse('ttt', args, setup=setup)
904 store_dir = get_store_dir(args)
905 try:
906 store = gf.Store(store_dir)
907 store.make_travel_time_tables(force=options.force)
909 except gf.StoreError as e:
910 die(e)
913def command_tttview(args):
915 def setup(parser):
916 parser.add_option(
917 '--source-depth', dest='source_depth', type=float,
918 help='Source depth in km')
920 parser.add_option(
921 '--receiver-depth', dest='receiver_depth', type=float,
922 help='Receiver depth in km')
924 parser, options, args = cl_parse(
925 'tttview', args, setup=setup,
926 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
928 try:
929 phase_ids = args.pop().split(',')
930 except Exception:
931 parser.error('cannot get <phase-ids> argument')
933 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
936def command_sat(args):
937 def setup(parser):
938 parser.add_option(
939 '--force', dest='force', action='store_true',
940 help='overwrite existing files')
942 parser.add_option(
943 '--attribute',
944 action='store',
945 dest='attribute',
946 type='choice',
947 choices=gf.store.available_stored_tables[1::],
948 default='takeoff_angle',
949 help='calculate interpolation table for selected ray attributes.')
951 parser, options, args = cl_parse('sat', args, setup=setup)
953 store_dir = get_store_dir(args)
954 try:
955 store = gf.Store(store_dir)
956 store.make_stored_table(options.attribute, force=options.force)
958 except gf.StoreError as e:
959 die(e)
962def command_satview(args):
964 def setup(parser):
965 parser.add_option(
966 '--source-depth', dest='source_depth', type=float,
967 help='Source depth in km')
969 parser.add_option(
970 '--receiver-depth', dest='receiver_depth', type=float,
971 help='Receiver depth in km')
973 parser.add_option(
974 '--attribute',
975 action='store',
976 dest='attribute',
977 type='choice',
978 choices=gf.store.available_stored_tables[1::],
979 default='takeoff_angle',
980 help='view selected ray attribute.')
982 parser, options, args = cl_parse(
983 'satview', args, setup=setup,
984 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
986 try:
987 phase_ids = args.pop().split(',')
988 except Exception:
989 parser.error('cannot get <phase-ids> argument')
991 logger.info('Plotting stored attribute %s' % options.attribute)
993 stored_attribute_table_plots(
994 phase_ids, options, args, attribute=options.attribute)
997def command_tttextract(args):
998 def setup(parser):
999 parser.add_option(
1000 '--output', dest='output_fn', metavar='TEMPLATE',
1001 help='output to text files instead of stdout '
1002 '(example TEMPLATE: "extracted/%(args)s.txt")')
1004 parser, options, args = cl_parse('tttextract', args, setup=setup)
1005 try:
1006 sdef = args.pop()
1007 except Exception:
1008 parser.error('cannot get <selection> argument')
1010 try:
1011 sphase = args.pop()
1012 except Exception:
1013 parser.error('cannot get <phase> argument')
1015 try:
1016 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1017 except gf.meta.InvalidTimingSpecification:
1018 parser.error('invalid phase specification: "%s"' % sphase)
1020 try:
1021 gdef = gf.meta.parse_grid_spec(sdef)
1022 except gf.meta.GridSpecError as e:
1023 die(e)
1025 store_dir = get_store_dir(args)
1027 try:
1028 store = gf.Store(store_dir)
1029 for args in store.config.iter_extraction(gdef, level=-1):
1030 s = ['%e' % x for x in args]
1031 for phase in phases:
1032 t = store.t(phase, args)
1033 if t is not None:
1034 s.append('%e' % t)
1035 else:
1036 s.append('nan')
1038 if options.output_fn:
1039 d = dict(
1040 args='_'.join('%e' % x for x in args),
1041 extension='txt')
1043 fn = options.output_fn % d
1044 util.ensuredirs(fn)
1045 with open(fn, 'a') as f:
1046 f.write(' '.join(s))
1047 f.write('\n')
1048 else:
1049 print(' '.join(s))
1051 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1052 die(e)
1055def command_tttlsd(args):
1057 def setup(parser):
1058 pass
1060 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1062 try:
1063 sphase_ids = args.pop()
1064 except Exception:
1065 parser.error('cannot get <phase> argument')
1067 try:
1068 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1069 except gf.meta.InvalidTimingSpecification:
1070 parser.error('invalid phase specification: "%s"' % sphase_ids)
1072 store_dir = get_store_dir(args)
1074 try:
1075 store = gf.Store(store_dir)
1076 for phase_id in phase_ids:
1077 store.fix_ttt_holes(phase_id)
1079 except gf.StoreError as e:
1080 die(e)
1083def command_server(args):
1084 from pyrocko.gf import server
1086 def setup(parser):
1087 parser.add_option(
1088 '--port', dest='port', metavar='PORT', type='int', default=8080,
1089 help='serve on port PORT')
1091 parser.add_option(
1092 '--ip', dest='ip', metavar='IP', default='',
1093 help='serve on ip address IP')
1095 parser, options, args = cl_parse('server', args, setup=setup)
1097 engine = gf.LocalEngine(store_superdirs=args)
1098 server.run(options.ip, options.port, engine)
1101def command_download(args):
1102 from pyrocko.gf import ws
1104 details = '''
1106Browse pre-calculated Green's function stores online:
1108 https://greens-mill.pyrocko.org
1109'''
1111 def setup(parser):
1112 parser.add_option(
1113 '--force', dest='force', action='store_true',
1114 help='overwrite existing files')
1116 parser, options, args = cl_parse(
1117 'download', args, setup=setup, details=details)
1118 if len(args) not in (1, 2):
1119 sys.exit(parser.format_help())
1121 if len(args) == 2:
1122 site, store_id = args
1123 if not re.match(gf.meta.StringID.pattern, store_id):
1124 die('invalid store ID')
1125 else:
1126 site, store_id = args[0], None
1128 if site not in gf.ws.g_site_abbr:
1129 if -1 == site.find('://'):
1130 site = 'http://' + site
1132 try:
1133 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1134 except ws.DownloadError as e:
1135 die(str(e))
1138def command_modelview(args):
1140 import matplotlib.pyplot as plt
1141 import numpy as num
1142 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1143 mpl_init()
1145 neat_labels = {
1146 'vp': '$v_p$',
1147 'vs': '$v_s$',
1148 'qp': '$Q_p$',
1149 'qs': '$Q_s$',
1150 'rho': '$\\rho$'}
1152 def setup(parser):
1153 parser.add_option(
1154 '--parameters', dest='parameters',
1155 default='vp,vs', metavar='vp,vs,....',
1156 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1158 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1160 parser, options, args = cl_parse('modelview', args, setup=setup)
1162 store_dirs = get_store_dirs(args)
1164 parameters = options.parameters.split(',')
1166 fig, axes = plt.subplots(1,
1167 len(parameters),
1168 sharey=True,
1169 figsize=(len(parameters)*3, 5))
1171 if not isinstance(axes, num.ndarray):
1172 axes = [axes]
1174 axes = dict(zip(parameters, axes))
1176 for store_id in store_dirs:
1177 try:
1178 store = gf.Store(store_id)
1179 mod = store.config.earthmodel_1d
1180 z = mod.profile('z')
1182 for p in parameters:
1183 ax = axes[p]
1184 labelspace(ax)
1185 if '/' in p:
1186 p1, p2 = p.split('/')
1187 profile = mod.profile(p1)/mod.profile(p2)
1188 else:
1189 profile = mod.profile(p)
1191 ax.plot(profile, -z, label=store_id.split('/')[-1])
1192 if p in ['vs', 'vp', 'rho']:
1193 xscaled(1./1000., ax)
1195 yscaled(1./km, ax)
1197 except gf.StoreError as e:
1198 die(e)
1200 for p, ax in axes.items():
1201 ax.grid()
1202 if p in neat_labels:
1203 lab = neat_labels[p]
1204 elif all(x in neat_labels for x in p.split('/')):
1205 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1206 else:
1207 lab = p
1209 ax.set_title(lab, y=1.02)
1211 if p in units:
1212 lab += ' ' + units[p]
1214 ax.autoscale()
1215 ax.set_xlabel(lab)
1217 axes[parameters[0]].set_ylabel('Depth [km]')
1219 handles, labels = ax.get_legend_handles_labels()
1221 if len(store_dirs) > 1:
1222 ax.legend(
1223 handles,
1224 labels,
1225 bbox_to_anchor=(0.5, 0.12),
1226 bbox_transform=fig.transFigure,
1227 ncol=3,
1228 loc='upper center',
1229 fancybox=True)
1231 plt.subplots_adjust(bottom=0.22,
1232 wspace=0.05)
1234 mpl_show(plt)
1237def command_upgrade(args):
1238 parser, options, args = cl_parse('upgrade', args)
1239 store_dirs = get_store_dirs(args)
1240 try:
1241 for store_dir in store_dirs:
1242 store = gf.Store(store_dir)
1243 nup = store.upgrade()
1244 if nup == 0:
1245 print('%s: already up-to-date.' % store_dir)
1246 else:
1247 print('%s: %i file%s upgraded.' % (
1248 store_dir, nup, ['s', ''][nup == 1]))
1250 except gf.StoreError as e:
1251 die(e)
1254def command_addref(args):
1255 parser, options, args = cl_parse('addref', args)
1257 store_dirs = []
1258 filenames = []
1259 for arg in args:
1260 if op.isdir(arg):
1261 store_dirs.append(arg)
1262 elif op.isfile(arg):
1263 filenames.append(arg)
1264 else:
1265 die('invalid argument: %s' % arg)
1267 if not store_dirs:
1268 store_dirs.append('.')
1270 references = []
1271 try:
1272 for filename in args:
1273 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1274 except ImportError:
1275 die('pybtex module must be installed to use this function')
1277 if not references:
1278 die('no references found')
1280 for store_dir in store_dirs:
1281 try:
1282 store = gf.Store(store_dir)
1283 ids = [ref.id for ref in store.config.references]
1284 for ref in references:
1285 if ref.id in ids:
1286 die('duplicate reference id: %s' % ref.id)
1288 ids.append(ref.id)
1289 store.config.references.append(ref)
1291 store.save_config(make_backup=True)
1293 except gf.StoreError as e:
1294 die(e)
1297def command_qc(args):
1298 parser, options, args = cl_parse('qc', args)
1300 store_dir = get_store_dir(args)
1302 try:
1303 store = gf.Store(store_dir)
1304 s = store.stats()
1305 if s['empty'] != 0:
1306 print('has empty records')
1308 for aname in ['author', 'author_email', 'description', 'references']:
1310 if not getattr(store.config, aname):
1311 print('%s empty' % aname)
1313 except gf.StoreError as e:
1314 die(e)
1317def command_report(args):
1318 from pyrocko.fomosto.report import report_call
1319 report_call.run_program(args)
1322def main(args=None):
1323 if args is None:
1324 args = sys.argv[1:]
1326 if len(args) < 1:
1327 sys.exit('Usage: %s' % usage)
1329 command = args.pop(0)
1331 if command in subcommands:
1332 globals()['command_' + command](args)
1334 elif command in ('--help', '-h', 'help'):
1335 if command == 'help' and args:
1336 acommand = args[0]
1337 if acommand in subcommands:
1338 globals()['command_' + acommand](['--help'])
1340 sys.exit('Usage: %s' % usage)
1342 else:
1343 sys.exit('fomosto: error: no such subcommand: %s' % command)
1346if __name__ == '__main__':
1347 main()