Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/apps/fomosto.py: 60%
705 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-12 12:00 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Fomosto - Green's function database management.
8'''
10import sys
11import re
12import os.path as op
13import logging
14import copy
15import shutil
16from optparse import OptionParser
18from pyrocko import util, trace, gf, cake, io, fomosto
19from pyrocko.gui.snuffler import marker
20from pyrocko.util import mpl_show
22logger = logging.getLogger('pyrocko.apps.fomosto')
23km = 1e3
26def d2u(d):
27 return dict((k.replace('-', '_'), v) for (k, v) in d.items())
30subcommand_descriptions = {
31 'init': 'create a new empty GF store',
32 'build': 'compute GFs and fill into store',
33 'stats': 'print information about a GF store',
34 'check': 'check for problems in GF store',
35 'decimate': 'build decimated variant of a GF store',
36 'redeploy': 'copy traces from one GF store into another',
37 'view': 'view selected traces',
38 'extract': 'extract selected traces',
39 'import': 'convert Kiwi GFDB to GF store format',
40 'export': 'convert GF store to Kiwi GFDB format',
41 'ttt': 'create travel time tables',
42 'tttview': 'plot travel time table',
43 'tttextract': 'extract selected travel times',
44 'tttlsd': 'fix holes in travel time tables',
45 'sat': 'create stored ray attribute table',
46 'satview': 'plot stored ray attribute table',
47 'server': 'run seismosizer server',
48 'download': 'download GF store from a server',
49 'modelview': 'plot earthmodels',
50 'upgrade': 'upgrade store format to latest version',
51 'addref': 'import citation references to GF store config',
52 'qc': 'quality check',
53 'report': "report for Green's Function databases",
54}
56subcommand_usages = {
57 'init': ['init <type> <store-dir> [options]',
58 'init redeploy <source> <destination> [options]'],
59 'build': 'build [store-dir] [options]',
60 'stats': 'stats [store-dir] [options]',
61 'check': 'check [store-dir] [options]',
62 'decimate': 'decimate [store-dir] <factor> [options]',
63 'redeploy': 'redeploy <source> <destination> [options]',
64 'view': 'view [store-dir] ... [options]',
65 'extract': 'extract [store-dir] <selection>',
66 'import': 'import <source> <destination> [options]',
67 'export': 'export [store-dir] <destination> [options]',
68 'ttt': 'ttt [store-dir] [options]',
69 'tttview': 'tttview [store-dir] <phase-ids> [options]',
70 'tttextract': 'tttextract [store-dir] <phase> <selection>',
71 'tttlsd': 'tttlsd [store-dir] <phase>',
72 'sat': 'sat [store-dir] [options]',
73 'satview': 'satview [store-dir] <phase-ids> [options]',
74 'server': 'server [options] <store-super-dir> ...',
75 'download': 'download [options] <site> <store-id>',
76 'modelview': 'modelview <selection>',
77 'upgrade': 'upgrade [store-dir] ...',
78 'addref': 'addref [store-dir] ... <filename> ...',
79 'qc': 'qc [store-dir]',
80 'report': 'report <subcommand> <arguments>... [options]'
81}
83subcommands = subcommand_descriptions.keys()
85program_name = 'fomosto'
87usage = program_name + ''' <subcommand> <arguments> ... [options]
89Subcommands:
91 init %(init)s
92 build %(build)s
93 stats %(stats)s
94 check %(check)s
95 decimate %(decimate)s
96 redeploy %(redeploy)s
97 view %(view)s
98 extract %(extract)s
99 import %(import)s
100 export %(export)s
101 ttt %(ttt)s
102 tttview %(tttview)s
103 tttextract %(tttextract)s
104 tttlsd %(tttlsd)s
105 sat %(sat)s
106 satview %(satview)s
107 server %(server)s
108 download %(download)s
109 modelview %(modelview)s
110 upgrade %(upgrade)s
111 addref %(addref)s
112 qc %(qc)s
113 report %(report)s
115To get further help and a list of available options for any subcommand run:
117 fomosto <subcommand> --help
119''' % d2u(subcommand_descriptions)
122def add_common_options(parser):
123 parser.add_option(
124 '--loglevel',
125 action='store',
126 dest='loglevel',
127 type='choice',
128 choices=('critical', 'error', 'warning', 'info', 'debug'),
129 default='info',
130 help='set logger level to '
131 '"critical", "error", "warning", "info", or "debug". '
132 'Default is "%default".')
135def process_common_options(options):
136 util.setup_logging(program_name, options.loglevel)
139def cl_parse(command, args, setup=None, details=None):
140 usage = subcommand_usages[command]
141 descr = subcommand_descriptions[command]
143 if isinstance(usage, str):
144 usage = [usage]
146 susage = '%s %s' % (program_name, usage[0])
147 for s in usage[1:]:
148 susage += '\n%s%s %s' % (' '*7, program_name, s)
150 description = descr[0].upper() + descr[1:] + '.'
152 if details:
153 description = description + ' %s' % details
155 parser = OptionParser(usage=susage, description=description)
156 parser.format_description = lambda formatter: description
158 if setup:
159 setup(parser)
161 add_common_options(parser)
162 (options, args) = parser.parse_args(args)
163 process_common_options(options)
164 return parser, options, args
167def die(message, err=''):
168 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
171def fomo_wrapper_module(name):
172 try:
173 if not re.match(gf.meta.StringID.pattern, name):
174 raise ValueError('invalid name')
176 words = name.split('.', 1)
177 if len(words) == 2:
178 name, variant = words
179 else:
180 name = words[0]
181 variant = None
183 name_clean = re.sub(r'[.-]', '_', name)
184 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
185 mod = __import__(modname, level=0)
186 return getattr(mod.fomosto, name_clean), variant
188 except ValueError:
189 die('invalid modelling code wrapper name')
191 except ImportError:
192 die('''modelling code wrapper "%s" not available or not installed
193 (module probed: "%s")''' % (name, modname))
196def command_init(args):
198 details = '''
200Available modelling backends:
201%s
203 More information at
204 https://pyrocko.org/docs/current/apps/fomosto/backends.html
205''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
207 parser, options, args = cl_parse(
208 'init', args,
209 details=details)
211 if len(args) == 0:
212 sys.exit(parser.format_help())
214 if args[0] == 'redeploy':
215 if len(args) != 3:
216 parser.error('incorrect number of arguments')
218 source_dir, dest_dir = args[1:]
220 try:
221 source = gf.Store(source_dir)
222 except gf.StoreError as e:
223 die(e)
225 config = copy.deepcopy(source.config)
226 config.derived_from_id = source.config.id
227 try:
228 config_filenames = gf.store.Store.create_editables(
229 dest_dir, config=config)
231 except gf.StoreError as e:
232 die(e)
234 try:
235 dest = gf.Store(dest_dir)
236 except gf.StoreError as e:
237 die(e)
239 for k in source.extra_keys():
240 source_fn = source.get_extra_path(k)
241 dest_fn = dest.get_extra_path(k)
242 shutil.copyfile(source_fn, dest_fn)
244 logger.info(
245 '(1) configure settings in files:\n %s'
246 % '\n '.join(config_filenames))
248 logger.info(
249 '(2) run "fomosto redeploy <source> <dest>", as needed')
251 else:
252 if len(args) != 2:
253 parser.error('incorrect number of arguments')
255 (modelling_code_id, store_dir) = args
257 module, variant = fomo_wrapper_module(modelling_code_id)
258 try:
259 config_filenames = module.init(store_dir, variant)
260 except gf.StoreError as e:
261 die(e)
263 logger.info('(1) configure settings in files:\n %s'
264 % '\n '.join(config_filenames))
265 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
266 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
269def get_store_dir(args):
270 if len(args) == 1:
271 store_dir = op.abspath(args.pop(0))
272 else:
273 store_dir = op.abspath(op.curdir)
275 if not op.isdir(store_dir):
276 die('not a directory: %s' % store_dir)
278 return store_dir
281def get_store_dirs(args):
282 if len(args) == 0:
283 store_dirs = [op.abspath(op.curdir)]
284 else:
285 store_dirs = [op.abspath(x) for x in args]
287 for store_dir in store_dirs:
288 if not op.isdir(store_dir):
289 die('not a directory: %s' % store_dir)
291 return store_dirs
294def command_build(args):
296 def setup(parser):
297 parser.add_option(
298 '--force', dest='force', action='store_true',
299 help='overwrite existing files')
301 parser.add_option(
302 '--nworkers', dest='nworkers', type='int', metavar='N',
303 help='run N worker processes in parallel')
305 parser.add_option(
306 '--continue', dest='continue_', action='store_true',
307 help='continue suspended build')
309 parser.add_option(
310 '--step', dest='step', type='int', metavar='I',
311 help='process block number IBLOCK')
313 parser.add_option(
314 '--block', dest='iblock', type='int', metavar='I',
315 help='process block number IBLOCK')
317 parser, options, args = cl_parse('build', args, setup=setup)
319 store_dir = get_store_dir(args)
320 try:
321 if options.step is not None:
322 step = options.step - 1
323 else:
324 step = None
326 if options.iblock is not None:
327 iblock = options.iblock - 1
328 else:
329 iblock = None
331 store = gf.Store(store_dir)
332 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
333 module.build(
334 store_dir,
335 force=options.force,
336 nworkers=options.nworkers, continue_=options.continue_,
337 step=step,
338 iblock=iblock)
340 except gf.StoreError as e:
341 die(e)
344def command_stats(args):
346 parser, options, args = cl_parse('stats', args)
347 store_dir = get_store_dir(args)
349 try:
350 store = gf.Store(store_dir)
351 s = store.stats()
353 except gf.StoreError as e:
354 die(e)
356 for k in store.stats_keys:
357 print('%s: %s' % (k, s[k]))
360def command_check(args):
362 parser, options, args = cl_parse('check', args)
363 store_dir = get_store_dir(args)
365 try:
366 store = gf.Store(store_dir)
367 problems = store.check(show_progress=True)
368 if problems:
369 die('problems detected with gf store: %s' % store_dir)
371 except gf.StoreError as e:
372 die(e)
375def load_config(fn):
376 try:
377 config = gf.meta.load(filename=fn)
378 assert isinstance(config, gf.Config)
380 except Exception:
381 die('cannot load gf config from file: %s' % fn)
383 return config
386def command_decimate(args):
388 def setup(parser):
389 parser.add_option(
390 '--config', dest='config_fn', metavar='FILE',
391 help='use modified spacial sampling given in FILE')
393 parser.add_option(
394 '--force', dest='force', action='store_true',
395 help='overwrite existing files')
397 parser, options, args = cl_parse('decimate', args, setup=setup)
398 try:
399 decimate = int(args.pop())
400 except Exception:
401 parser.error('cannot get <factor> argument')
403 store_dir = get_store_dir(args)
405 config = None
406 if options.config_fn:
407 config = load_config(options.config_fn)
409 try:
410 store = gf.Store(store_dir)
411 store.make_decimated(decimate, config=config, force=options.force,
412 show_progress=True)
414 except gf.StoreError as e:
415 die(e)
418def sindex(args):
419 return '(%s)' % ', '.join('%g' % x for x in args)
422def command_redeploy(args):
424 parser, options, args = cl_parse('redeploy', args)
426 if not len(args) == 2:
427 sys.exit(parser.format_help())
429 source_store_dir, dest_store_dir = args
431 try:
432 source = gf.Store(source_store_dir)
433 except gf.StoreError as e:
434 die(e)
436 try:
437 gf.store.Store.create_dependants(dest_store_dir)
438 except gf.StoreError:
439 pass
441 try:
442 dest = gf.Store(dest_store_dir, 'w')
444 except gf.StoreError as e:
445 die(e)
447 show_progress = True
449 try:
450 if show_progress:
451 pbar = util.progressbar('redeploying', dest.config.nrecords)
453 for i, args in enumerate(dest.config.iter_nodes()):
454 try:
455 tr = source.get(args, interpolation='off')
456 dest.put(args, tr)
458 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa
459 logger.debug('skipping %s, (%s)' % (sindex(args), e))
461 except gf.store.StoreError as e:
462 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
464 if show_progress:
465 pbar.update(i+1)
467 finally:
468 if show_progress:
469 pbar.finish()
472def command_view(args):
473 def setup(parser):
474 parser.add_option(
475 '--extract',
476 dest='extract',
477 metavar='start:stop[:step|@num],...',
478 help='specify which traces to show')
480 parser.add_option(
481 '--show-phases',
482 dest='showphases',
483 default=None,
484 metavar='[phase_id1,phase_id2,...|all]',
485 help='add phase markers from ttt')
487 parser.add_option(
488 '--opengl',
489 dest='opengl',
490 action='store_true',
491 default=None,
492 help='use OpenGL for drawing')
494 parser.add_option(
495 '--no-opengl',
496 dest='opengl',
497 action='store_false',
498 default=None,
499 help='do not use OpenGL for drawing')
501 parser, options, args = cl_parse('view', args, setup=setup)
503 gdef = None
504 if options.extract:
505 try:
506 gdef = gf.meta.parse_grid_spec(options.extract)
507 except gf.meta.GridSpecError as e:
508 die(e)
510 store_dirs = get_store_dirs(args)
512 alpha = 'abcdefghijklmnopqrstxyz'.upper()
514 markers = []
515 traces = []
517 try:
518 for istore, store_dir in enumerate(store_dirs):
519 store = gf.Store(store_dir)
520 if options.showphases == 'all':
521 phasenames = [pn.id for pn in store.config.tabulated_phases]
522 elif options.showphases is not None:
523 phasenames = options.showphases.split(',')
525 ii = 0
526 for args in store.config.iter_extraction(gdef):
527 gtr = store.get(args)
529 loc_code = ''
530 if len(store_dirs) > 1:
531 loc_code = alpha[istore % len(alpha)]
533 if gtr:
535 sta_code = '%04i (%s)' % (
536 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
537 tmin = gtr.deltat * gtr.itmin
539 tr = trace.Trace(
540 '',
541 sta_code,
542 loc_code,
543 '%02i' % args[-1],
544 ydata=gtr.data,
545 deltat=gtr.deltat,
546 tmin=tmin)
548 if options.showphases:
549 for phasename in phasenames:
550 phase_tmin = store.t(phasename, args[:-1])
551 if phase_tmin:
552 m = marker.PhaseMarker(
553 [('',
554 sta_code,
555 loc_code,
556 '%02i' % args[-1])],
557 phase_tmin,
558 phase_tmin,
559 0,
560 phasename=phasename)
561 markers.append(m)
563 traces.append(tr)
565 ii += 1
567 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
568 die(e)
570 def setup(win):
571 viewer = win.get_view()
572 viewer.menuitem_showboxes.setChecked(False)
573 viewer.menuitem_colortraces.setChecked(True)
574 viewer.menuitem_demean.setChecked(False)
575 viewer.menuitems_sorting[5][0].setChecked(True)
576 viewer.menuitems_scaling[2][0].setChecked(True)
577 viewer.sortingmode_change()
578 viewer.scalingmode_change()
580 trace.snuffle(
581 traces, markers=markers, opengl=options.opengl, launch_hook=setup)
584def command_extract(args):
585 def setup(parser):
586 parser.add_option(
587 '--format', dest='format', default='mseed',
588 choices=['mseed', 'sac', 'text', 'yaff'],
589 help='export to format "mseed", "sac", "text", or "yaff". '
590 'Default is "mseed".')
592 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
593 parser.add_option(
594 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
595 help='output path template [default: "%s"]' % fndfl)
597 parser, options, args = cl_parse('extract', args, setup=setup)
598 try:
599 sdef = args.pop()
600 except Exception:
601 parser.error('cannot get <selection> argument')
603 try:
604 gdef = gf.meta.parse_grid_spec(sdef)
605 except gf.meta.GridSpecError as e:
606 die(e)
608 store_dir = get_store_dir(args)
610 extensions = {
611 'mseed': 'mseed',
612 'sac': 'sac',
613 'text': 'txt',
614 'yaff': 'yaff'}
616 try:
617 store = gf.Store(store_dir)
618 for args in store.config.iter_extraction(gdef):
619 gtr = store.get(args)
620 if gtr:
621 tr = trace.Trace(
622 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
623 ydata=gtr.data,
624 deltat=gtr.deltat,
625 tmin=gtr.deltat * gtr.itmin)
627 additional = dict(
628 args='_'.join('%g' % x for x in args),
629 irecord=store.str_irecord(args),
630 extension=extensions[options.format])
632 io.save(
633 tr,
634 options.output_fn,
635 format=options.format,
636 additional=additional)
638 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
639 die(e)
642def command_import(args):
643 try:
644 from tunguska import gfdb
645 except ImportError:
646 die('the kiwi tools must be installed to use this feature')
648 parser, options, args = cl_parse('import', args)
650 show_progress = True
652 if not len(args) == 2:
653 sys.exit(parser.format_help())
655 source_path, dest_store_dir = args
657 if op.isdir(source_path):
658 source_path = op.join(source_path, 'db')
660 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
662 db = gfdb.Gfdb(source_path)
664 config = gf.meta.ConfigTypeA(
665 id='imported_gfs',
666 distance_min=db.firstx,
667 distance_max=db.firstx + (db.nx-1) * db.dx,
668 distance_delta=db.dx,
669 source_depth_min=db.firstz,
670 source_depth_max=db.firstz + (db.nz-1) * db.dz,
671 source_depth_delta=db.dz,
672 sample_rate=1.0/db.dt,
673 ncomponents=db.ng
674 )
676 try:
677 gf.store.Store.create(dest_store_dir, config=config)
678 dest = gf.Store(dest_store_dir, 'w')
679 try:
680 if show_progress:
681 pbar = util.progressbar(
682 'importing', dest.config.nrecords/dest.config.ncomponents)
684 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
685 source_depth, distance = [float(x) for x in args]
686 traces = db.get_traces_pyrocko(distance, source_depth)
687 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
688 for ig in range(db.ng):
689 if ig in ig_to_trace:
690 tr = ig_to_trace[ig]
691 gf_tr = gf.store.GFTrace(
692 tr.get_ydata(),
693 int(round(tr.tmin / tr.deltat)),
694 tr.deltat)
696 else:
697 gf_tr = gf.store.Zero
699 dest.put((source_depth, distance, ig), gf_tr)
701 if show_progress:
702 pbar.update(i+1)
704 finally:
705 if show_progress:
706 pbar.finish()
708 dest.close()
710 except gf.StoreError as e:
711 die(e)
714def command_export(args):
715 from subprocess import Popen, PIPE
717 try:
718 from tunguska import gfdb
719 except ImportError as err:
720 die('the kiwi tools must be installed to use this feature', err)
722 def setup(parser):
723 parser.add_option(
724 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
725 help='split output gfdb into N chunks')
727 parser, options, args = cl_parse('export', args, setup=setup)
729 show_progress = True
731 if len(args) not in (1, 2):
732 sys.exit(parser.format_help())
734 target_path = args.pop()
735 if op.isdir(target_path):
736 target_path = op.join(target_path, 'kiwi_gfdb')
737 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
739 source_store_dir = get_store_dir(args)
741 source = gf.Store(source_store_dir, 'r')
742 config = source.config
744 if not isinstance(config, gf.meta.ConfigTypeA):
745 die('only stores of type A can be exported to Kiwi format')
747 if op.isfile(target_path + '.index'):
748 die('destation already exists')
750 cmd = [str(x) for x in [
751 'gfdb_build',
752 target_path,
753 options.nchunks,
754 config.ndistances,
755 config.nsource_depths,
756 config.ncomponents,
757 config.deltat,
758 config.distance_delta,
759 config.source_depth_delta,
760 config.distance_min,
761 config.source_depth_min]]
763 p = Popen(cmd, stdin=PIPE)
764 p.communicate()
766 out_db = gfdb.Gfdb(target_path)
768 try:
769 if show_progress:
770 pbar = util.progressbar(
771 'exporting', config.nrecords/config.ncomponents)
773 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
775 data_out = []
776 for ig in range(config.ncomponents):
777 try:
778 tr = source.get((z, x, ig), interpolation='off')
779 data_out.append((tr.t, tr.data * config.factor))
781 except gf.store.StoreError as e:
782 logger.warning(
783 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
784 data_out.append(None)
786 # put a zero valued sample to no-data zero-traces at a compatible
787 # time
788 tmins = [
789 entry[0][0]
790 for entry in data_out
791 if entry is not None and entry[0].size != 0]
793 if tmins:
794 tmin = min(tmins)
795 for entry in data_out:
796 if entry is not None and entry[0].size == 0:
797 entry[0].resize(1)
798 entry[1].resize(1)
799 entry[0][0] = tmin
800 entry[1][0] = 0.0
802 out_db.put_traces_slow(x, z, data_out)
804 if show_progress:
805 pbar.update(i+1)
807 source.close()
809 finally:
810 if show_progress:
811 pbar.finish()
814def phasedef_or_horvel(x):
815 try:
816 return float(x)
817 except ValueError:
818 return cake.PhaseDef(x)
821def mkp(s):
822 return [phasedef_or_horvel(ps) for ps in s.split(',')]
825def stored_attribute_table_plots(phase_ids, options, args, attribute):
826 import numpy as num
827 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
829 plt = mpl_init()
831 np = 1
832 store_dir = get_store_dir(args)
833 for phase_id in phase_ids:
834 try:
835 store = gf.Store(store_dir)
837 if options.receiver_depth is not None:
838 receiver_depth = options.receiver_depth * 1000.0
839 else:
840 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
841 receiver_depth = store.config.receiver_depth
843 elif isinstance(store.config, gf.ConfigTypeB):
844 receiver_depth = store.config.receiver_depth_min
846 else:
847 receiver_depth = 0.0
849 phase = store.get_stored_phase(phase_id, attribute)
850 axes = plt.subplot(2, len(phase_ids), np)
851 labelspace(axes)
852 xscaled(1. / km, axes)
853 yscaled(1. / km, axes)
854 x = None
855 if isinstance(store.config, gf.ConfigTypeB):
856 x = (receiver_depth, None, None)
858 phase.plot_2d(axes, x=x)
859 axes.set_title(phase_id)
860 np += 1
861 except gf.StoreError as e:
862 die(e)
864 axes = plt.subplot(2, 1, 2)
865 num_d = 100
866 distances = num.linspace(store.config.distance_min,
867 store.config.distance_max,
868 num_d)
870 if options.source_depth is not None:
871 source_depth = options.source_depth * km
872 else:
873 source_depth = store.config.source_depth_min + (
874 store.config.source_depth_max - store.config.source_depth_min) / 2.
876 if isinstance(store.config, gf.ConfigTypeA):
877 attribute_vals = num.empty(num_d)
878 for phase_id in phase_ids:
879 attribute_vals[:] = num.NAN
880 for i, d in enumerate(distances):
881 if attribute == 'phase':
882 attribute_vals[i] = store.t(phase_id, (source_depth, d))
883 ylabel = 'TT [s]'
884 else:
885 attribute_vals[i] = store.get_stored_attribute(
886 phase_id, options.attribute, (source_depth, d))
887 ylabel = '%s [deg]' % options.attribute
889 axes.plot(distances / km, attribute_vals, label=phase_id)
891 axes.set_title('source source_depth %s km' % (source_depth / km))
892 axes.set_xlabel('distance [km]')
893 axes.set_ylabel(ylabel)
894 axes.legend()
896 plt.tight_layout()
897 mpl_show(plt)
900def command_ttt(args):
901 def setup(parser):
902 parser.add_option(
903 '--force', dest='force', action='store_true',
904 help='overwrite existing files')
906 parser, options, args = cl_parse('ttt', args, setup=setup)
908 store_dir = get_store_dir(args)
909 try:
910 store = gf.Store(store_dir)
911 store.make_travel_time_tables(force=options.force)
913 except gf.StoreError as e:
914 die(e)
917def command_tttview(args):
919 def setup(parser):
920 parser.add_option(
921 '--source-depth', dest='source_depth', type=float,
922 help='Source depth in km')
924 parser.add_option(
925 '--receiver-depth', dest='receiver_depth', type=float,
926 help='Receiver depth in km')
928 parser, options, args = cl_parse(
929 'tttview', args, setup=setup,
930 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
932 try:
933 phase_ids = args.pop().split(',')
934 except Exception:
935 parser.error('cannot get <phase-ids> argument')
937 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
940def command_sat(args):
941 def setup(parser):
942 parser.add_option(
943 '--force', dest='force', action='store_true',
944 help='overwrite existing files')
946 parser.add_option(
947 '--attribute',
948 action='store',
949 dest='attribute',
950 type='choice',
951 choices=gf.store.available_stored_tables[1::],
952 default='takeoff_angle',
953 help='calculate interpolation table for selected ray attributes.')
955 parser, options, args = cl_parse('sat', args, setup=setup)
957 store_dir = get_store_dir(args)
958 try:
959 store = gf.Store(store_dir)
960 store.make_stored_table(options.attribute, force=options.force)
962 except gf.StoreError as e:
963 die(e)
966def command_satview(args):
968 def setup(parser):
969 parser.add_option(
970 '--source-depth', dest='source_depth', type=float,
971 help='Source depth in km')
973 parser.add_option(
974 '--receiver-depth', dest='receiver_depth', type=float,
975 help='Receiver depth in km')
977 parser.add_option(
978 '--attribute',
979 action='store',
980 dest='attribute',
981 type='choice',
982 choices=gf.store.available_stored_tables[1::],
983 default='takeoff_angle',
984 help='view selected ray attribute.')
986 parser, options, args = cl_parse(
987 'satview', args, setup=setup,
988 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
990 try:
991 phase_ids = args.pop().split(',')
992 except Exception:
993 parser.error('cannot get <phase-ids> argument')
995 logger.info('Plotting stored attribute %s' % options.attribute)
997 stored_attribute_table_plots(
998 phase_ids, options, args, attribute=options.attribute)
1001def command_tttextract(args):
1002 def setup(parser):
1003 parser.add_option(
1004 '--output', dest='output_fn', metavar='TEMPLATE',
1005 help='output to text files instead of stdout '
1006 '(example TEMPLATE: "extracted/%(args)s.txt")')
1008 parser, options, args = cl_parse('tttextract', args, setup=setup)
1009 try:
1010 sdef = args.pop()
1011 except Exception:
1012 parser.error('cannot get <selection> argument')
1014 try:
1015 sphase = args.pop()
1016 except Exception:
1017 parser.error('cannot get <phase> argument')
1019 try:
1020 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1021 except gf.meta.InvalidTimingSpecification:
1022 parser.error('invalid phase specification: "%s"' % sphase)
1024 try:
1025 gdef = gf.meta.parse_grid_spec(sdef)
1026 except gf.meta.GridSpecError as e:
1027 die(e)
1029 store_dir = get_store_dir(args)
1031 try:
1032 store = gf.Store(store_dir)
1033 for args in store.config.iter_extraction(gdef, level=-1):
1034 s = ['%e' % x for x in args]
1035 for phase in phases:
1036 t = store.t(phase, args)
1037 if t is not None:
1038 s.append('%e' % t)
1039 else:
1040 s.append('nan')
1042 if options.output_fn:
1043 d = dict(
1044 args='_'.join('%e' % x for x in args),
1045 extension='txt')
1047 fn = options.output_fn % d
1048 util.ensuredirs(fn)
1049 with open(fn, 'a') as f:
1050 f.write(' '.join(s))
1051 f.write('\n')
1052 else:
1053 print(' '.join(s))
1055 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1056 die(e)
1059def command_tttlsd(args):
1061 def setup(parser):
1062 pass
1064 parser, options, args = cl_parse('tttlsd', args, setup=setup)
1066 try:
1067 sphase_ids = args.pop()
1068 except Exception:
1069 parser.error('cannot get <phase> argument')
1071 try:
1072 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1073 except gf.meta.InvalidTimingSpecification:
1074 parser.error('invalid phase specification: "%s"' % sphase_ids)
1076 store_dir = get_store_dir(args)
1078 try:
1079 store = gf.Store(store_dir)
1080 for phase_id in phase_ids:
1081 store.fix_ttt_holes(phase_id)
1083 except gf.StoreError as e:
1084 die(e)
1087def command_server(args):
1088 from pyrocko.gf import server
1090 def setup(parser):
1091 parser.add_option(
1092 '--port', dest='port', metavar='PORT', type='int', default=8080,
1093 help='serve on port PORT')
1095 parser.add_option(
1096 '--ip', dest='ip', metavar='IP', default='',
1097 help='serve on ip address IP')
1099 parser, options, args = cl_parse('server', args, setup=setup)
1101 engine = gf.LocalEngine(store_superdirs=args)
1102 server.run(options.ip, options.port, engine)
1105def command_download(args):
1106 from pyrocko.gf import ws
1108 details = '''
1110Browse pre-calculated Green's function stores online:
1112 https://greens-mill.pyrocko.org
1113'''
1115 def setup(parser):
1116 parser.add_option(
1117 '--force', dest='force', action='store_true',
1118 help='overwrite existing files')
1120 parser, options, args = cl_parse(
1121 'download', args, setup=setup, details=details)
1122 if len(args) not in (1, 2):
1123 sys.exit(parser.format_help())
1125 if len(args) == 2:
1126 site, store_id = args
1127 if not re.match(gf.meta.StringID.pattern, store_id):
1128 die('invalid store ID')
1129 else:
1130 site, store_id = args[0], None
1132 if site not in gf.ws.g_site_abbr:
1133 if -1 == site.find('://'):
1134 site = 'http://' + site
1136 try:
1137 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1138 except ws.DownloadError as e:
1139 die(str(e))
1142def command_modelview(args):
1144 import matplotlib.pyplot as plt
1145 import numpy as num
1146 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1147 mpl_init()
1149 neat_labels = {
1150 'vp': '$v_p$',
1151 'vs': '$v_s$',
1152 'qp': '$Q_p$',
1153 'qs': '$Q_s$',
1154 'rho': '$\\rho$'}
1156 def setup(parser):
1157 parser.add_option(
1158 '--parameters', dest='parameters',
1159 default='vp,vs', metavar='vp,vs,....',
1160 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1162 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1164 parser, options, args = cl_parse('modelview', args, setup=setup)
1166 store_dirs = get_store_dirs(args)
1168 parameters = options.parameters.split(',')
1170 fig, axes = plt.subplots(1,
1171 len(parameters),
1172 sharey=True,
1173 figsize=(len(parameters)*3, 5))
1175 if not isinstance(axes, num.ndarray):
1176 axes = [axes]
1178 axes = dict(zip(parameters, axes))
1180 for store_id in store_dirs:
1181 try:
1182 store = gf.Store(store_id)
1183 mod = store.config.earthmodel_1d
1184 z = mod.profile('z')
1186 for p in parameters:
1187 ax = axes[p]
1188 labelspace(ax)
1189 if '/' in p:
1190 p1, p2 = p.split('/')
1191 profile = mod.profile(p1)/mod.profile(p2)
1192 else:
1193 profile = mod.profile(p)
1195 ax.plot(profile, -z, label=store_id.split('/')[-1])
1196 if p in ['vs', 'vp', 'rho']:
1197 xscaled(1./1000., ax)
1199 yscaled(1./km, ax)
1201 except gf.StoreError as e:
1202 die(e)
1204 for p, ax in axes.items():
1205 ax.grid()
1206 if p in neat_labels:
1207 lab = neat_labels[p]
1208 elif all(x in neat_labels for x in p.split('/')):
1209 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1210 else:
1211 lab = p
1213 ax.set_title(lab, y=1.02)
1215 if p in units:
1216 lab += ' ' + units[p]
1218 ax.autoscale()
1219 ax.set_xlabel(lab)
1221 axes[parameters[0]].set_ylabel('Depth [km]')
1223 handles, labels = ax.get_legend_handles_labels()
1225 if len(store_dirs) > 1:
1226 ax.legend(
1227 handles,
1228 labels,
1229 bbox_to_anchor=(0.5, 0.12),
1230 bbox_transform=fig.transFigure,
1231 ncol=3,
1232 loc='upper center',
1233 fancybox=True)
1235 plt.subplots_adjust(bottom=0.22,
1236 wspace=0.05)
1238 mpl_show(plt)
1241def command_upgrade(args):
1242 parser, options, args = cl_parse('upgrade', args)
1243 store_dirs = get_store_dirs(args)
1244 try:
1245 for store_dir in store_dirs:
1246 store = gf.Store(store_dir)
1247 nup = store.upgrade()
1248 if nup == 0:
1249 print('%s: already up-to-date.' % store_dir)
1250 else:
1251 print('%s: %i file%s upgraded.' % (
1252 store_dir, nup, ['s', ''][nup == 1]))
1254 except gf.StoreError as e:
1255 die(e)
1258def command_addref(args):
1259 parser, options, args = cl_parse('addref', args)
1261 store_dirs = []
1262 filenames = []
1263 for arg in args:
1264 if op.isdir(arg):
1265 store_dirs.append(arg)
1266 elif op.isfile(arg):
1267 filenames.append(arg)
1268 else:
1269 die('invalid argument: %s' % arg)
1271 if not store_dirs:
1272 store_dirs.append('.')
1274 references = []
1275 try:
1276 for filename in args:
1277 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1278 except ImportError:
1279 die('pybtex module must be installed to use this function')
1281 if not references:
1282 die('no references found')
1284 for store_dir in store_dirs:
1285 try:
1286 store = gf.Store(store_dir)
1287 ids = [ref.id for ref in store.config.references]
1288 for ref in references:
1289 if ref.id in ids:
1290 die('duplicate reference id: %s' % ref.id)
1292 ids.append(ref.id)
1293 store.config.references.append(ref)
1295 store.save_config(make_backup=True)
1297 except gf.StoreError as e:
1298 die(e)
1301def command_qc(args):
1302 parser, options, args = cl_parse('qc', args)
1304 store_dir = get_store_dir(args)
1306 try:
1307 store = gf.Store(store_dir)
1308 s = store.stats()
1309 if s['empty'] != 0:
1310 print('has empty records')
1312 for aname in ['author', 'author_email', 'description', 'references']:
1314 if not getattr(store.config, aname):
1315 print('%s empty' % aname)
1317 except gf.StoreError as e:
1318 die(e)
1321def command_report(args):
1322 from pyrocko.fomosto.report import report_call
1323 report_call.run_program(args)
1326def main(args=None):
1327 '''
1328 CLI entry point for Pyrocko's ``fomosto`` app.
1329 '''
1330 if args is None:
1331 args = sys.argv[1:]
1333 if len(args) < 1:
1334 sys.exit('Usage: %s' % usage)
1336 command = args.pop(0)
1338 if command in subcommands:
1339 globals()['command_' + command](args)
1341 elif command in ('--help', '-h', 'help'):
1342 if command == 'help' and args:
1343 acommand = args[0]
1344 if acommand in subcommands:
1345 globals()['command_' + acommand](['--help'])
1347 sys.exit('Usage: %s' % usage)
1349 else:
1350 sys.exit('fomosto: error: no such subcommand: %s' % command)
1353if __name__ == '__main__':
1354 main()