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 trace.snuffle(traces, markers=markers, opengl=options.opengl)
570def command_extract(args):
571 def setup(parser):
572 parser.add_option(
573 '--format', dest='format', default='mseed',
574 choices=['mseed', 'sac', 'text', 'yaff'],
575 help='export to format "mseed", "sac", "text", or "yaff". '
576 'Default is "mseed".')
578 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
579 parser.add_option(
580 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
581 help='output path template [default: "%s"]' % fndfl)
583 parser, options, args = cl_parse('extract', args, setup=setup)
584 try:
585 sdef = args.pop()
586 except Exception:
587 parser.error('cannot get <selection> argument')
589 try:
590 gdef = gf.meta.parse_grid_spec(sdef)
591 except gf.meta.GridSpecError as e:
592 die(e)
594 store_dir = get_store_dir(args)
596 extensions = {
597 'mseed': 'mseed',
598 'sac': 'sac',
599 'text': 'txt',
600 'yaff': 'yaff'}
602 try:
603 store = gf.Store(store_dir)
604 for args in store.config.iter_extraction(gdef):
605 gtr = store.get(args)
606 if gtr:
607 tr = trace.Trace(
608 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
609 ydata=gtr.data,
610 deltat=gtr.deltat,
611 tmin=gtr.deltat * gtr.itmin)
613 additional = dict(
614 args='_'.join('%g' % x for x in args),
615 irecord=store.str_irecord(args),
616 extension=extensions[options.format])
618 io.save(
619 tr,
620 options.output_fn,
621 format=options.format,
622 additional=additional)
624 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
625 die(e)
628def command_import(args):
629 try:
630 from tunguska import gfdb
631 except ImportError:
632 die('the kiwi tools must be installed to use this feature')
634 parser, options, args = cl_parse('import', args)
636 show_progress = True
638 if not len(args) == 2:
639 sys.exit(parser.format_help())
641 source_path, dest_store_dir = args
643 if op.isdir(source_path):
644 source_path = op.join(source_path, 'db')
646 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
648 db = gfdb.Gfdb(source_path)
650 config = gf.meta.ConfigTypeA(
651 id='imported_gfs',
652 distance_min=db.firstx,
653 distance_max=db.firstx + (db.nx-1) * db.dx,
654 distance_delta=db.dx,
655 source_depth_min=db.firstz,
656 source_depth_max=db.firstz + (db.nz-1) * db.dz,
657 source_depth_delta=db.dz,
658 sample_rate=1.0/db.dt,
659 ncomponents=db.ng
660 )
662 try:
663 gf.store.Store.create(dest_store_dir, config=config)
664 dest = gf.Store(dest_store_dir, 'w')
665 try:
666 if show_progress:
667 pbar = util.progressbar(
668 'importing', dest.config.nrecords/dest.config.ncomponents)
670 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
671 source_depth, distance = [float(x) for x in args]
672 traces = db.get_traces_pyrocko(distance, source_depth)
673 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
674 for ig in range(db.ng):
675 if ig in ig_to_trace:
676 tr = ig_to_trace[ig]
677 gf_tr = gf.store.GFTrace(
678 tr.get_ydata(),
679 int(round(tr.tmin / tr.deltat)),
680 tr.deltat)
682 else:
683 gf_tr = gf.store.Zero
685 dest.put((source_depth, distance, ig), gf_tr)
687 if show_progress:
688 pbar.update(i+1)
690 finally:
691 if show_progress:
692 pbar.finish()
694 dest.close()
696 except gf.StoreError as e:
697 die(e)
700def command_export(args):
701 from subprocess import Popen, PIPE
703 try:
704 from tunguska import gfdb
705 except ImportError as err:
706 die('the kiwi tools must be installed to use this feature', err)
708 def setup(parser):
709 parser.add_option(
710 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
711 help='split output gfdb into N chunks')
713 parser, options, args = cl_parse('export', args, setup=setup)
715 show_progress = True
717 if len(args) not in (1, 2):
718 sys.exit(parser.format_help())
720 target_path = args.pop()
721 if op.isdir(target_path):
722 target_path = op.join(target_path, 'kiwi_gfdb')
723 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
725 source_store_dir = get_store_dir(args)
727 source = gf.Store(source_store_dir, 'r')
728 config = source.config
730 if not isinstance(config, gf.meta.ConfigTypeA):
731 die('only stores of type A can be exported to Kiwi format')
733 if op.isfile(target_path + '.index'):
734 die('destation already exists')
736 cmd = [str(x) for x in [
737 'gfdb_build',
738 target_path,
739 options.nchunks,
740 config.ndistances,
741 config.nsource_depths,
742 config.ncomponents,
743 config.deltat,
744 config.distance_delta,
745 config.source_depth_delta,
746 config.distance_min,
747 config.source_depth_min]]
749 p = Popen(cmd, stdin=PIPE)
750 p.communicate()
752 out_db = gfdb.Gfdb(target_path)
754 try:
755 if show_progress:
756 pbar = util.progressbar(
757 'exporting', config.nrecords/config.ncomponents)
759 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
761 data_out = []
762 for ig in range(config.ncomponents):
763 try:
764 tr = source.get((z, x, ig), interpolation='off')
765 data_out.append((tr.t, tr.data * config.factor))
767 except gf.store.StoreError as e:
768 logger.warning(
769 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
770 data_out.append(None)
772 # put a zero valued sample to no-data zero-traces at a compatible
773 # time
774 tmins = [
775 entry[0][0]
776 for entry in data_out
777 if entry is not None and entry[0].size != 0]
779 if tmins:
780 tmin = min(tmins)
781 for entry in data_out:
782 if entry is not None and entry[0].size == 0:
783 entry[0].resize(1)
784 entry[1].resize(1)
785 entry[0][0] = tmin
786 entry[1][0] = 0.0
788 out_db.put_traces_slow(x, z, data_out)
790 if show_progress:
791 pbar.update(i+1)
793 source.close()
795 finally:
796 if show_progress:
797 pbar.finish()
800def phasedef_or_horvel(x):
801 try:
802 return float(x)
803 except ValueError:
804 return cake.PhaseDef(x)
807def mkp(s):
808 return [phasedef_or_horvel(ps) for ps in s.split(',')]
811def stored_attribute_table_plots(phase_ids, options, args, attribute):
812 import numpy as num
813 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
815 plt = mpl_init()
817 np = 1
818 store_dir = get_store_dir(args)
819 for phase_id in phase_ids:
820 try:
821 store = gf.Store(store_dir)
823 if options.receiver_depth is not None:
824 receiver_depth = options.receiver_depth * 1000.0
825 else:
826 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
827 receiver_depth = store.config.receiver_depth
829 elif isinstance(store.config, gf.ConfigTypeB):
830 receiver_depth = store.config.receiver_depth_min
832 else:
833 receiver_depth = 0.0
835 phase = store.get_stored_phase(phase_id, attribute)
836 axes = plt.subplot(2, len(phase_ids), np)
837 labelspace(axes)
838 xscaled(1. / km, axes)
839 yscaled(1. / km, axes)
840 x = None
841 if isinstance(store.config, gf.ConfigTypeB):
842 x = (receiver_depth, None, None)
844 phase.plot_2d(axes, x=x)
845 axes.set_title(phase_id)
846 np += 1
847 except gf.StoreError as e:
848 die(e)
850 axes = plt.subplot(2, 1, 2)
851 num_d = 100
852 distances = num.linspace(store.config.distance_min,
853 store.config.distance_max,
854 num_d)
856 if options.source_depth is not None:
857 source_depth = options.source_depth * km
858 else:
859 source_depth = store.config.source_depth_min + (
860 store.config.source_depth_max - store.config.source_depth_min) / 2.
862 if isinstance(store.config, gf.ConfigTypeA):
863 attribute_vals = num.empty(num_d)
864 for phase_id in phase_ids:
865 attribute_vals[:] = num.NAN
866 for i, d in enumerate(distances):
867 if attribute == 'phase':
868 attribute_vals[i] = store.t(phase_id, (source_depth, d))
869 ylabel = 'TT [s]'
870 else:
871 attribute_vals[i] = store.get_stored_attribute(
872 phase_id, options.attribute, (source_depth, d))
873 ylabel = '%s [deg]' % options.attribute
875 axes.plot(distances / km, attribute_vals, label=phase_id)
877 axes.set_title('source source_depth %s km' % (source_depth / km))
878 axes.set_xlabel('distance [km]')
879 axes.set_ylabel(ylabel)
880 axes.legend()
882 plt.tight_layout()
883 mpl_show(plt)
886def command_ttt(args):
887 def setup(parser):
888 parser.add_option(
889 '--force', dest='force', action='store_true',
890 help='overwrite existing files')
892 parser, options, args = cl_parse('ttt', args, setup=setup)
894 store_dir = get_store_dir(args)
895 try:
896 store = gf.Store(store_dir)
897 store.make_travel_time_tables(force=options.force)
899 except gf.StoreError as e:
900 die(e)
903def command_tttview(args):
905 def setup(parser):
906 parser.add_option(
907 '--source-depth', dest='source_depth', type=float,
908 help='Source depth in km')
910 parser.add_option(
911 '--receiver-depth', dest='receiver_depth', type=float,
912 help='Receiver depth in km')
914 parser, options, args = cl_parse(
915 'tttview', args, setup=setup,
916 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
918 try:
919 phase_ids = args.pop().split(',')
920 except Exception:
921 parser.error('cannot get <phase-ids> argument')
923 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
926def command_sat(args):
927 def setup(parser):
928 parser.add_option(
929 '--force', dest='force', action='store_true',
930 help='overwrite existing files')
932 parser.add_option(
933 '--attribute',
934 action='store',
935 dest='attribute',
936 type='choice',
937 choices=gf.store.available_stored_tables[1::],
938 default='takeoff_angle',
939 help='calculate interpolation table for selected ray attributes.')
941 parser, options, args = cl_parse('sat', args, setup=setup)
943 store_dir = get_store_dir(args)
944 try:
945 store = gf.Store(store_dir)
946 store.make_stored_table(options.attribute, force=options.force)
948 except gf.StoreError as e:
949 die(e)
952def command_satview(args):
954 def setup(parser):
955 parser.add_option(
956 '--source-depth', dest='source_depth', type=float,
957 help='Source depth in km')
959 parser.add_option(
960 '--receiver-depth', dest='receiver_depth', type=float,
961 help='Receiver depth in km')
963 parser.add_option(
964 '--attribute',
965 action='store',
966 dest='attribute',
967 type='choice',
968 choices=gf.store.available_stored_tables[1::],
969 default='takeoff_angle',
970 help='view selected ray attribute.')
972 parser, options, args = cl_parse(
973 'satview', args, setup=setup,
974 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
976 try:
977 phase_ids = args.pop().split(',')
978 except Exception:
979 parser.error('cannot get <phase-ids> argument')
981 logger.info('Plotting stored attribute %s' % options.attribute)
983 stored_attribute_table_plots(
984 phase_ids, options, args, attribute=options.attribute)
987def command_tttextract(args):
988 def setup(parser):
989 parser.add_option(
990 '--output', dest='output_fn', metavar='TEMPLATE',
991 help='output to text files instead of stdout '
992 '(example TEMPLATE: "extracted/%(args)s.txt")')
994 parser, options, args = cl_parse('tttextract', args, setup=setup)
995 try:
996 sdef = args.pop()
997 except Exception:
998 parser.error('cannot get <selection> argument')
1000 try:
1001 sphase = args.pop()
1002 except Exception:
1003 parser.error('cannot get <phase> argument')
1005 try:
1006 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1007 except gf.meta.InvalidTimingSpecification:
1008 parser.error('invalid phase specification: "%s"' % sphase)
1010 try:
1011 gdef = gf.meta.parse_grid_spec(sdef)
1012 except gf.meta.GridSpecError as e:
1013 die(e)
1015 store_dir = get_store_dir(args)
1017 try:
1018 store = gf.Store(store_dir)
1019 for args in store.config.iter_extraction(gdef, level=-1):
1020 s = ['%e' % x for x in args]
1021 for phase in phases:
1022 t = store.t(phase, args)
1023 if t is not None:
1024 s.append('%e' % t)
1025 else:
1026 s.append('nan')
1028 if options.output_fn:
1029 d = dict(
1030 args='_'.join('%e' % x for x in args),
1031 extension='txt')
1033 fn = options.output_fn % d
1034 util.ensuredirs(fn)
1035 with open(fn, 'a') as f:
1036 f.write(' '.join(s))
1037 f.write('\n')
1038 else:
1039 print(' '.join(s))
1041 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1042 die(e)
1045def command_tttlsd(args):
1047 def setup(parser):
1048 pass
1050 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1052 try:
1053 sphase_ids = args.pop()
1054 except Exception:
1055 parser.error('cannot get <phase> argument')
1057 try:
1058 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1059 except gf.meta.InvalidTimingSpecification:
1060 parser.error('invalid phase specification: "%s"' % sphase_ids)
1062 store_dir = get_store_dir(args)
1064 try:
1065 store = gf.Store(store_dir)
1066 for phase_id in phase_ids:
1067 store.fix_ttt_holes(phase_id)
1069 except gf.StoreError as e:
1070 die(e)
1073def command_server(args):
1074 from pyrocko.gf import server
1076 def setup(parser):
1077 parser.add_option(
1078 '--port', dest='port', metavar='PORT', type='int', default=8080,
1079 help='serve on port PORT')
1081 parser.add_option(
1082 '--ip', dest='ip', metavar='IP', default='',
1083 help='serve on ip address IP')
1085 parser, options, args = cl_parse('server', args, setup=setup)
1087 engine = gf.LocalEngine(store_superdirs=args)
1088 server.run(options.ip, options.port, engine)
1091def command_download(args):
1092 from pyrocko.gf import ws
1094 details = '''
1096Browse pre-calculated Green's function stores online:
1098 https://greens-mill.pyrocko.org
1099'''
1101 def setup(parser):
1102 parser.add_option(
1103 '--force', dest='force', action='store_true',
1104 help='overwrite existing files')
1106 parser, options, args = cl_parse(
1107 'download', args, setup=setup, details=details)
1108 if not len(args) in (1, 2):
1109 sys.exit(parser.format_help())
1111 if len(args) == 2:
1112 site, store_id = args
1113 if not re.match(gf.meta.StringID.pattern, store_id):
1114 die('invalid store ID')
1115 else:
1116 site, store_id = args[0], None
1118 if site not in gf.ws.g_site_abbr:
1119 if -1 == site.find('://'):
1120 site = 'http://' + site
1122 try:
1123 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1124 except ws.DownloadError as e:
1125 die(str(e))
1128def command_modelview(args):
1130 import matplotlib.pyplot as plt
1131 import numpy as num
1132 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1133 mpl_init()
1135 neat_labels = {
1136 'vp': '$v_p$',
1137 'vs': '$v_s$',
1138 'qp': '$Q_p$',
1139 'qs': '$Q_s$',
1140 'rho': '$\\rho$'}
1142 def setup(parser):
1143 parser.add_option(
1144 '--parameters', dest='parameters',
1145 default='vp,vs', metavar='vp,vs,....',
1146 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1148 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1150 parser, options, args = cl_parse('modelview', args, setup=setup)
1152 store_dirs = get_store_dirs(args)
1154 parameters = options.parameters.split(',')
1156 fig, axes = plt.subplots(1,
1157 len(parameters),
1158 sharey=True,
1159 figsize=(len(parameters)*3, 5))
1161 if not isinstance(axes, num.ndarray):
1162 axes = [axes]
1164 axes = dict(zip(parameters, axes))
1166 for store_id in store_dirs:
1167 try:
1168 store = gf.Store(store_id)
1169 mod = store.config.earthmodel_1d
1170 z = mod.profile('z')
1172 for p in parameters:
1173 ax = axes[p]
1174 labelspace(ax)
1175 if '/' in p:
1176 p1, p2 = p.split('/')
1177 profile = mod.profile(p1)/mod.profile(p2)
1178 else:
1179 profile = mod.profile(p)
1181 ax.plot(profile, -z, label=store_id.split('/')[-1])
1182 if p in ['vs', 'vp', 'rho']:
1183 xscaled(1./1000., ax)
1185 yscaled(1./km, ax)
1187 except gf.StoreError as e:
1188 die(e)
1190 for p, ax in axes.items():
1191 ax.grid()
1192 if p in neat_labels:
1193 lab = neat_labels[p]
1194 elif all(x in neat_labels for x in p.split('/')):
1195 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1196 else:
1197 lab = p
1199 ax.set_title(lab, y=1.02)
1201 if p in units:
1202 lab += ' ' + units[p]
1204 ax.autoscale()
1205 ax.set_xlabel(lab)
1207 axes[parameters[0]].set_ylabel('Depth [km]')
1209 handles, labels = ax.get_legend_handles_labels()
1211 if len(store_dirs) > 1:
1212 ax.legend(
1213 handles,
1214 labels,
1215 bbox_to_anchor=(0.5, 0.12),
1216 bbox_transform=fig.transFigure,
1217 ncol=3,
1218 loc='upper center',
1219 fancybox=True)
1221 plt.subplots_adjust(bottom=0.22,
1222 wspace=0.05)
1224 mpl_show(plt)
1227def command_upgrade(args):
1228 parser, options, args = cl_parse('upgrade', args)
1229 store_dirs = get_store_dirs(args)
1230 try:
1231 for store_dir in store_dirs:
1232 store = gf.Store(store_dir)
1233 nup = store.upgrade()
1234 if nup == 0:
1235 print('%s: already up-to-date.' % store_dir)
1236 else:
1237 print('%s: %i file%s upgraded.' % (
1238 store_dir, nup, ['s', ''][nup == 1]))
1240 except gf.StoreError as e:
1241 die(e)
1244def command_addref(args):
1245 parser, options, args = cl_parse('addref', args)
1247 store_dirs = []
1248 filenames = []
1249 for arg in args:
1250 if op.isdir(arg):
1251 store_dirs.append(arg)
1252 elif op.isfile(arg):
1253 filenames.append(arg)
1254 else:
1255 die('invalid argument: %s' % arg)
1257 if not store_dirs:
1258 store_dirs.append('.')
1260 references = []
1261 try:
1262 for filename in args:
1263 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1264 except ImportError:
1265 die('pybtex module must be installed to use this function')
1267 if not references:
1268 die('no references found')
1270 for store_dir in store_dirs:
1271 try:
1272 store = gf.Store(store_dir)
1273 ids = [ref.id for ref in store.config.references]
1274 for ref in references:
1275 if ref.id in ids:
1276 die('duplicate reference id: %s' % ref.id)
1278 ids.append(ref.id)
1279 store.config.references.append(ref)
1281 store.save_config(make_backup=True)
1283 except gf.StoreError as e:
1284 die(e)
1287def command_qc(args):
1288 parser, options, args = cl_parse('qc', args)
1290 store_dir = get_store_dir(args)
1292 try:
1293 store = gf.Store(store_dir)
1294 s = store.stats()
1295 if s['empty'] != 0:
1296 print('has empty records')
1298 for aname in ['author', 'author_email', 'description', 'references']:
1300 if not getattr(store.config, aname):
1301 print('%s empty' % aname)
1303 except gf.StoreError as e:
1304 die(e)
1307def command_report(args):
1308 from pyrocko.fomosto.report import report_call
1309 report_call.run_program(args)
1312def main(args=None):
1313 if args is None:
1314 args = sys.argv[1:]
1316 if len(args) < 1:
1317 sys.exit('Usage: %s' % usage)
1319 command = args.pop(0)
1321 if command in subcommands:
1322 globals()['command_' + command](args)
1324 elif command in ('--help', '-h', 'help'):
1325 if command == 'help' and args:
1326 acommand = args[0]
1327 if acommand in subcommands:
1328 globals()['command_' + acommand](['--help'])
1330 sys.exit('Usage: %s' % usage)
1332 else:
1333 sys.exit('fomosto: error: no such subcommand: %s' % command)
1336if __name__ == '__main__':
1337 main(sys.argv[1:])