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 if show_progress:
448 pbar = util.progressbar('redeploying', dest.config.nrecords)
450 for i, args in enumerate(dest.config.iter_nodes()):
451 try:
452 tr = source.get(args, interpolation='off')
453 dest.put(args, tr)
455 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e:
456 logger.debug('skipping %s, (%s)' % (sindex(args), e))
458 except gf.store.StoreError as e:
459 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
461 if show_progress:
462 pbar.update(i+1)
464 if show_progress:
465 pbar.finish()
468def command_view(args):
469 def setup(parser):
470 parser.add_option(
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 trace.snuffle(traces, markers=markers, opengl=options.opengl)
569def command_extract(args):
570 def setup(parser):
571 parser.add_option(
572 '--format', dest='format', default='mseed',
573 choices=['mseed', 'sac', 'text', 'yaff'],
574 help='export to format "mseed", "sac", "text", or "yaff". '
575 'Default is "mseed".')
577 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
578 parser.add_option(
579 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
580 help='output path template [default: "%s"]' % fndfl)
582 parser, options, args = cl_parse('extract', args, setup=setup)
583 try:
584 sdef = args.pop()
585 except Exception:
586 parser.error('cannot get <selection> argument')
588 try:
589 gdef = gf.meta.parse_grid_spec(sdef)
590 except gf.meta.GridSpecError as e:
591 die(e)
593 store_dir = get_store_dir(args)
595 extensions = {
596 'mseed': 'mseed',
597 'sac': 'sac',
598 'text': 'txt',
599 'yaff': 'yaff'}
601 try:
602 store = gf.Store(store_dir)
603 for args in store.config.iter_extraction(gdef):
604 gtr = store.get(args)
605 if gtr:
606 tr = trace.Trace(
607 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
608 ydata=gtr.data,
609 deltat=gtr.deltat,
610 tmin=gtr.deltat * gtr.itmin)
612 additional = dict(
613 args='_'.join('%g' % x for x in args),
614 irecord=store.str_irecord(args),
615 extension=extensions[options.format])
617 io.save(
618 tr,
619 options.output_fn,
620 format=options.format,
621 additional=additional)
623 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
624 die(e)
627def command_import(args):
628 try:
629 from tunguska import gfdb
630 except ImportError:
631 die('the kiwi tools must be installed to use this feature')
633 parser, options, args = cl_parse('import', args)
635 show_progress = True
637 if not len(args) == 2:
638 sys.exit(parser.format_help())
640 source_path, dest_store_dir = args
642 if op.isdir(source_path):
643 source_path = op.join(source_path, 'db')
645 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
647 db = gfdb.Gfdb(source_path)
649 config = gf.meta.ConfigTypeA(
650 id='imported_gfs',
651 distance_min=db.firstx,
652 distance_max=db.firstx + (db.nx-1) * db.dx,
653 distance_delta=db.dx,
654 source_depth_min=db.firstz,
655 source_depth_max=db.firstz + (db.nz-1) * db.dz,
656 source_depth_delta=db.dz,
657 sample_rate=1.0/db.dt,
658 ncomponents=db.ng
659 )
661 try:
662 gf.store.Store.create(dest_store_dir, config=config)
663 dest = gf.Store(dest_store_dir, 'w')
664 if show_progress:
665 pbar = util.progressbar(
666 'importing', dest.config.nrecords/dest.config.ncomponents)
668 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
669 source_depth, distance = [float(x) for x in args]
670 traces = db.get_traces_pyrocko(distance, source_depth)
671 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
672 for ig in range(db.ng):
673 if ig in ig_to_trace:
674 tr = ig_to_trace[ig]
675 gf_tr = gf.store.GFTrace(
676 tr.get_ydata(),
677 int(round(tr.tmin / tr.deltat)),
678 tr.deltat)
680 else:
681 gf_tr = gf.store.Zero
683 dest.put((source_depth, distance, ig), gf_tr)
685 if show_progress:
686 pbar.update(i+1)
688 if show_progress:
689 pbar.finish()
691 dest.close()
693 except gf.StoreError as e:
694 die(e)
697def command_export(args):
698 from subprocess import Popen, PIPE
700 try:
701 from tunguska import gfdb
702 except ImportError as err:
703 die('the kiwi tools must be installed to use this feature', err)
705 def setup(parser):
706 parser.add_option(
707 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
708 help='split output gfdb into N chunks')
710 parser, options, args = cl_parse('export', args, setup=setup)
712 show_progress = True
714 if len(args) not in (1, 2):
715 sys.exit(parser.format_help())
717 target_path = args.pop()
718 if op.isdir(target_path):
719 target_path = op.join(target_path, 'kiwi_gfdb')
720 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
722 source_store_dir = get_store_dir(args)
724 source = gf.Store(source_store_dir, 'r')
725 config = source.config
727 if not isinstance(config, gf.meta.ConfigTypeA):
728 die('only stores of type A can be exported to Kiwi format')
730 if op.isfile(target_path + '.index'):
731 die('destation already exists')
733 cmd = [str(x) for x in [
734 'gfdb_build',
735 target_path,
736 options.nchunks,
737 config.ndistances,
738 config.nsource_depths,
739 config.ncomponents,
740 config.deltat,
741 config.distance_delta,
742 config.source_depth_delta,
743 config.distance_min,
744 config.source_depth_min]]
746 p = Popen(cmd, stdin=PIPE)
747 p.communicate()
749 out_db = gfdb.Gfdb(target_path)
751 if show_progress:
752 pbar = util.progressbar(
753 'exporting', config.nrecords/config.ncomponents)
755 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
757 data_out = []
758 for ig in range(config.ncomponents):
759 try:
760 tr = source.get((z, x, ig), interpolation='off')
761 data_out.append((tr.t, tr.data * config.factor))
763 except gf.store.StoreError as e:
764 logger.warning('cannot get %s, (%s)' % (sindex((z, x, ig)), e))
765 data_out.append(None)
767 # put a zero valued sample to no-data zero-traces at a compatible time
768 tmins = [
769 entry[0][0]
770 for entry in data_out
771 if entry is not None and entry[0].size != 0]
773 if tmins:
774 tmin = min(tmins)
775 for entry in data_out:
776 if entry is not None and entry[0].size == 0:
777 entry[0].resize(1)
778 entry[1].resize(1)
779 entry[0][0] = tmin
780 entry[1][0] = 0.0
782 out_db.put_traces_slow(x, z, data_out)
784 if show_progress:
785 pbar.update(i+1)
787 if show_progress:
788 pbar.finish()
790 source.close()
793def phasedef_or_horvel(x):
794 try:
795 return float(x)
796 except ValueError:
797 return cake.PhaseDef(x)
800def mkp(s):
801 return [phasedef_or_horvel(ps) for ps in s.split(',')]
804def stored_attribute_table_plots(phase_ids, options, args, attribute):
805 import numpy as num
806 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
808 plt = mpl_init()
810 np = 1
811 store_dir = get_store_dir(args)
812 for phase_id in phase_ids:
813 try:
814 store = gf.Store(store_dir)
816 if options.receiver_depth is not None:
817 receiver_depth = options.receiver_depth * 1000.0
818 else:
819 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
820 receiver_depth = store.config.receiver_depth
822 elif isinstance(store.config, gf.ConfigTypeB):
823 receiver_depth = store.config.receiver_depth_min
825 else:
826 receiver_depth = 0.0
828 phase = store.get_stored_phase(phase_id, attribute)
829 axes = plt.subplot(2, len(phase_ids), np)
830 labelspace(axes)
831 xscaled(1. / km, axes)
832 yscaled(1. / km, axes)
833 x = None
834 if isinstance(store.config, gf.ConfigTypeB):
835 x = (receiver_depth, None, None)
837 phase.plot_2d(axes, x=x)
838 axes.set_title(phase_id)
839 np += 1
840 except gf.StoreError as e:
841 die(e)
843 axes = plt.subplot(2, 1, 2)
844 num_d = 100
845 distances = num.linspace(store.config.distance_min,
846 store.config.distance_max,
847 num_d)
849 if options.source_depth is not None:
850 source_depth = options.source_depth * km
851 else:
852 source_depth = store.config.source_depth_min + (
853 store.config.source_depth_max - store.config.source_depth_min) / 2.
855 if isinstance(store.config, gf.ConfigTypeA):
856 attribute_vals = num.empty(num_d)
857 for phase_id in phase_ids:
858 attribute_vals[:] = num.NAN
859 for i, d in enumerate(distances):
860 if attribute == 'phase':
861 attribute_vals[i] = store.t(phase_id, (source_depth, d))
862 ylabel = 'TT [s]'
863 else:
864 attribute_vals[i] = store.get_stored_attribute(
865 phase_id, options.attribute, (source_depth, d))
866 ylabel = '%s [deg]' % options.attribute
868 axes.plot(distances / km, attribute_vals, label=phase_id)
870 axes.set_title('source source_depth %s km' % (source_depth / km))
871 axes.set_xlabel('distance [km]')
872 axes.set_ylabel(ylabel)
873 axes.legend()
875 plt.tight_layout()
876 mpl_show(plt)
879def command_ttt(args):
880 def setup(parser):
881 parser.add_option(
882 '--force', dest='force', action='store_true',
883 help='overwrite existing files')
885 parser, options, args = cl_parse('ttt', args, setup=setup)
887 store_dir = get_store_dir(args)
888 try:
889 store = gf.Store(store_dir)
890 store.make_travel_time_tables(force=options.force)
892 except gf.StoreError as e:
893 die(e)
896def command_tttview(args):
898 def setup(parser):
899 parser.add_option(
900 '--source-depth', dest='source_depth', type=float,
901 help='Source depth in km')
903 parser.add_option(
904 '--receiver-depth', dest='receiver_depth', type=float,
905 help='Receiver depth in km')
907 parser, options, args = cl_parse(
908 'tttview', args, setup=setup,
909 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
911 try:
912 phase_ids = args.pop().split(',')
913 except Exception:
914 parser.error('cannot get <phase-ids> argument')
916 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
919def command_sat(args):
920 def setup(parser):
921 parser.add_option(
922 '--force', dest='force', action='store_true',
923 help='overwrite existing files')
925 parser.add_option(
926 '--attribute',
927 action='store',
928 dest='attribute',
929 type='choice',
930 choices=gf.store.available_stored_tables[1::],
931 default='takeoff_angle',
932 help='calculate interpolation table for selected ray attributes.')
934 parser, options, args = cl_parse('sat', args, setup=setup)
936 store_dir = get_store_dir(args)
937 try:
938 store = gf.Store(store_dir)
939 store.make_stored_table(options.attribute, force=options.force)
941 except gf.StoreError as e:
942 die(e)
945def command_satview(args):
947 def setup(parser):
948 parser.add_option(
949 '--source-depth', dest='source_depth', type=float,
950 help='Source depth in km')
952 parser.add_option(
953 '--receiver-depth', dest='receiver_depth', type=float,
954 help='Receiver depth in km')
956 parser.add_option(
957 '--attribute',
958 action='store',
959 dest='attribute',
960 type='choice',
961 choices=gf.store.available_stored_tables[1::],
962 default='takeoff_angle',
963 help='view selected ray attribute.')
965 parser, options, args = cl_parse(
966 'satview', args, setup=setup,
967 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
969 try:
970 phase_ids = args.pop().split(',')
971 except Exception:
972 parser.error('cannot get <phase-ids> argument')
974 logger.info('Plotting stored attribute %s' % options.attribute)
976 stored_attribute_table_plots(
977 phase_ids, options, args, attribute=options.attribute)
980def command_tttextract(args):
981 def setup(parser):
982 parser.add_option(
983 '--output', dest='output_fn', metavar='TEMPLATE',
984 help='output to text files instead of stdout '
985 '(example TEMPLATE: "extracted/%(args)s.txt")')
987 parser, options, args = cl_parse('tttextract', args, setup=setup)
988 try:
989 sdef = args.pop()
990 except Exception:
991 parser.error('cannot get <selection> argument')
993 try:
994 sphase = args.pop()
995 except Exception:
996 parser.error('cannot get <phase> argument')
998 try:
999 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1000 except gf.meta.InvalidTimingSpecification:
1001 parser.error('invalid phase specification: "%s"' % sphase)
1003 try:
1004 gdef = gf.meta.parse_grid_spec(sdef)
1005 except gf.meta.GridSpecError as e:
1006 die(e)
1008 store_dir = get_store_dir(args)
1010 try:
1011 store = gf.Store(store_dir)
1012 for args in store.config.iter_extraction(gdef, level=-1):
1013 s = ['%e' % x for x in args]
1014 for phase in phases:
1015 t = store.t(phase, args)
1016 if t is not None:
1017 s.append('%e' % t)
1018 else:
1019 s.append('nan')
1021 if options.output_fn:
1022 d = dict(
1023 args='_'.join('%e' % x for x in args),
1024 extension='txt')
1026 fn = options.output_fn % d
1027 util.ensuredirs(fn)
1028 with open(fn, 'a') as f:
1029 f.write(' '.join(s))
1030 f.write('\n')
1031 else:
1032 print(' '.join(s))
1034 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1035 die(e)
1038def command_tttlsd(args):
1040 def setup(parser):
1041 pass
1043 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1045 try:
1046 sphase_ids = args.pop()
1047 except Exception:
1048 parser.error('cannot get <phase> argument')
1050 try:
1051 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1052 except gf.meta.InvalidTimingSpecification:
1053 parser.error('invalid phase specification: "%s"' % sphase_ids)
1055 store_dir = get_store_dir(args)
1057 try:
1058 store = gf.Store(store_dir)
1059 for phase_id in phase_ids:
1060 store.fix_ttt_holes(phase_id)
1062 except gf.StoreError as e:
1063 die(e)
1066def command_server(args):
1067 from pyrocko.gf import server
1069 def setup(parser):
1070 parser.add_option(
1071 '--port', dest='port', metavar='PORT', type='int', default=8080,
1072 help='serve on port PORT')
1074 parser.add_option(
1075 '--ip', dest='ip', metavar='IP', default='',
1076 help='serve on ip address IP')
1078 parser, options, args = cl_parse('server', args, setup=setup)
1080 engine = gf.LocalEngine(store_superdirs=args)
1081 server.run(options.ip, options.port, engine)
1084def command_download(args):
1085 from pyrocko.gf import ws
1087 details = '''
1089Browse pre-calculated Green's function stores online:
1091 https://greens-mill.pyrocko.org
1092'''
1094 def setup(parser):
1095 parser.add_option(
1096 '--force', dest='force', action='store_true',
1097 help='overwrite existing files')
1099 parser, options, args = cl_parse(
1100 'download', args, setup=setup, details=details)
1101 if not len(args) in (1, 2):
1102 sys.exit(parser.format_help())
1104 if len(args) == 2:
1105 site, store_id = args
1106 if not re.match(gf.meta.StringID.pattern, store_id):
1107 die('invalid store ID')
1108 else:
1109 site, store_id = args[0], None
1111 if site not in gf.ws.g_site_abbr:
1112 if -1 == site.find('://'):
1113 site = 'http://' + site
1115 try:
1116 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1117 except ws.DownloadError as e:
1118 die(str(e))
1121def command_modelview(args):
1123 import matplotlib.pyplot as plt
1124 import numpy as num
1125 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1126 mpl_init()
1128 neat_labels = {
1129 'vp': '$v_p$',
1130 'vs': '$v_s$',
1131 'qp': '$Q_p$',
1132 'qs': '$Q_s$',
1133 'rho': '$\\rho$'}
1135 def setup(parser):
1136 parser.add_option(
1137 '--parameters', dest='parameters',
1138 default='vp,vs', metavar='vp,vs,....',
1139 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1141 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1143 parser, options, args = cl_parse('modelview', args, setup=setup)
1145 store_dirs = get_store_dirs(args)
1147 parameters = options.parameters.split(',')
1149 fig, axes = plt.subplots(1,
1150 len(parameters),
1151 sharey=True,
1152 figsize=(len(parameters)*3, 5))
1154 if not isinstance(axes, num.ndarray):
1155 axes = [axes]
1157 axes = dict(zip(parameters, axes))
1159 for store_id in store_dirs:
1160 try:
1161 store = gf.Store(store_id)
1162 mod = store.config.earthmodel_1d
1163 z = mod.profile('z')
1165 for p in parameters:
1166 ax = axes[p]
1167 labelspace(ax)
1168 if '/' in p:
1169 p1, p2 = p.split('/')
1170 profile = mod.profile(p1)/mod.profile(p2)
1171 else:
1172 profile = mod.profile(p)
1174 ax.plot(profile, -z, label=store_id.split('/')[-1])
1175 if p in ['vs', 'vp', 'rho']:
1176 xscaled(1./1000., ax)
1178 yscaled(1./km, ax)
1180 except gf.StoreError as e:
1181 die(e)
1183 for p, ax in axes.items():
1184 ax.grid()
1185 if p in neat_labels:
1186 lab = neat_labels[p]
1187 elif all(x in neat_labels for x in p.split('/')):
1188 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1189 else:
1190 lab = p
1192 ax.set_title(lab, y=1.02)
1194 if p in units:
1195 lab += ' ' + units[p]
1197 ax.autoscale()
1198 ax.set_xlabel(lab)
1200 axes[parameters[0]].set_ylabel('Depth [km]')
1202 handles, labels = ax.get_legend_handles_labels()
1204 if len(store_dirs) > 1:
1205 ax.legend(
1206 handles,
1207 labels,
1208 bbox_to_anchor=(0.5, 0.12),
1209 bbox_transform=fig.transFigure,
1210 ncol=3,
1211 loc='upper center',
1212 fancybox=True)
1214 plt.subplots_adjust(bottom=0.22,
1215 wspace=0.05)
1217 mpl_show(plt)
1220def command_upgrade(args):
1221 parser, options, args = cl_parse('upgrade', args)
1222 store_dirs = get_store_dirs(args)
1223 try:
1224 for store_dir in store_dirs:
1225 store = gf.Store(store_dir)
1226 nup = store.upgrade()
1227 if nup == 0:
1228 print('%s: already up-to-date.' % store_dir)
1229 else:
1230 print('%s: %i file%s upgraded.' % (
1231 store_dir, nup, ['s', ''][nup == 1]))
1233 except gf.StoreError as e:
1234 die(e)
1237def command_addref(args):
1238 parser, options, args = cl_parse('addref', args)
1240 store_dirs = []
1241 filenames = []
1242 for arg in args:
1243 if op.isdir(arg):
1244 store_dirs.append(arg)
1245 elif op.isfile(arg):
1246 filenames.append(arg)
1247 else:
1248 die('invalid argument: %s' % arg)
1250 if not store_dirs:
1251 store_dirs.append('.')
1253 references = []
1254 try:
1255 for filename in args:
1256 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1257 except ImportError:
1258 die('pybtex module must be installed to use this function')
1260 if not references:
1261 die('no references found')
1263 for store_dir in store_dirs:
1264 try:
1265 store = gf.Store(store_dir)
1266 ids = [ref.id for ref in store.config.references]
1267 for ref in references:
1268 if ref.id in ids:
1269 die('duplicate reference id: %s' % ref.id)
1271 ids.append(ref.id)
1272 store.config.references.append(ref)
1274 store.save_config(make_backup=True)
1276 except gf.StoreError as e:
1277 die(e)
1280def command_qc(args):
1281 parser, options, args = cl_parse('qc', args)
1283 store_dir = get_store_dir(args)
1285 try:
1286 store = gf.Store(store_dir)
1287 s = store.stats()
1288 if s['empty'] != 0:
1289 print('has empty records')
1291 for aname in ['author', 'author_email', 'description', 'references']:
1293 if not getattr(store.config, aname):
1294 print('%s empty' % aname)
1296 except gf.StoreError as e:
1297 die(e)
1300def command_report(args):
1301 from pyrocko.fomosto.report import report_call
1302 report_call.run_program(args)
1305def main(args=None):
1306 if args is None:
1307 args = sys.argv[1:]
1309 if len(args) < 1:
1310 sys.exit('Usage: %s' % usage)
1312 command = args.pop(0)
1314 if command in subcommands:
1315 globals()['command_' + command](args)
1317 elif command in ('--help', '-h', 'help'):
1318 if command == 'help' and args:
1319 acommand = args[0]
1320 if acommand in subcommands:
1321 globals()['command_' + acommand](['--help'])
1323 sys.exit('Usage: %s' % usage)
1325 else:
1326 sys.exit('fomosto: error: no such subcommand: %s' % command)
1329if __name__ == '__main__':
1330 main(sys.argv[1:])