Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/apps/fomosto.py: 59%
728 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +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 'trans2rot': 'create rotational from translational store',
55 'fd2trans': 'extract elastic10 from elastic10_fd store',
56}
58subcommand_usages = {
59 'init': ['init <type> <store-dir> [options]',
60 'init redeploy <source> <destination> [options]'],
61 'build': 'build [store-dir] [options]',
62 'stats': 'stats [store-dir] [options]',
63 'check': 'check [store-dir] [options]',
64 'decimate': 'decimate [store-dir] <factor> [options]',
65 'redeploy': 'redeploy <source> <destination> [options]',
66 'view': 'view [store-dir] ... [options]',
67 'extract': 'extract [store-dir] <selection>',
68 'import': 'import <source> <destination> [options]',
69 'export': 'export [store-dir] <destination> [options]',
70 'ttt': 'ttt [store-dir] [options]',
71 'tttview': 'tttview [store-dir] <phase-ids> [options]',
72 'tttextract': 'tttextract [store-dir] <phase> <selection>',
73 'tttlsd': 'tttlsd [store-dir] <phase>',
74 'sat': 'sat [store-dir] [options]',
75 'satview': 'satview [store-dir] <phase-ids> [options]',
76 'server': 'server [options] <store-super-dir> ...',
77 'download': 'download [options] <site> <store-id>',
78 'modelview': 'modelview <selection>',
79 'upgrade': 'upgrade [store-dir] ...',
80 'addref': 'addref [store-dir] ... <filename> ...',
81 'qc': 'qc [store-dir]',
82 'report': 'report <subcommand> <arguments>... [options]',
83 'trans2rot': 'trans2rot <source> <destination>',
84 'fd2trans': 'fd2trans <source> <destination>',
85}
87subcommands = subcommand_descriptions.keys()
89program_name = 'fomosto'
91usage = program_name + ''' <subcommand> <arguments> ... [options]
93Subcommands:
95 init %(init)s
96 build %(build)s
97 stats %(stats)s
98 check %(check)s
99 decimate %(decimate)s
100 redeploy %(redeploy)s
101 view %(view)s
102 extract %(extract)s
103 import %(import)s
104 export %(export)s
105 ttt %(ttt)s
106 tttview %(tttview)s
107 tttextract %(tttextract)s
108 tttlsd %(tttlsd)s
109 sat %(sat)s
110 satview %(satview)s
111 server %(server)s
112 download %(download)s
113 modelview %(modelview)s
114 upgrade %(upgrade)s
115 addref %(addref)s
116 qc %(qc)s
117 report %(report)s
118 trans2rot %(trans2rot)s
119 fd2trans %(fd2trans)s
121To get further help and a list of available options for any subcommand run:
123 fomosto <subcommand> --help
125''' % d2u(subcommand_descriptions)
128def add_common_options(parser):
129 parser.add_option(
130 '--loglevel',
131 action='store',
132 dest='loglevel',
133 type='choice',
134 choices=('critical', 'error', 'warning', 'info', 'debug'),
135 default='info',
136 help='set logger level to '
137 '"critical", "error", "warning", "info", or "debug". '
138 'Default is "%default".')
141def process_common_options(options):
142 util.setup_logging(program_name, options.loglevel)
145def cl_parse(command, args, setup=None, details=None):
146 usage = subcommand_usages[command]
147 descr = subcommand_descriptions[command]
149 if isinstance(usage, str):
150 usage = [usage]
152 susage = '%s %s' % (program_name, usage[0])
153 for s in usage[1:]:
154 susage += '\n%s%s %s' % (' '*7, program_name, s)
156 description = descr[0].upper() + descr[1:] + '.'
158 if details:
159 description = description + ' %s' % details
161 parser = OptionParser(usage=susage, description=description)
162 parser.format_description = lambda formatter: description
164 if setup:
165 setup(parser)
167 add_common_options(parser)
168 (options, args) = parser.parse_args(args)
169 process_common_options(options)
170 return parser, options, args
173def die(message, err=''):
174 sys.exit('%s: error: %s \n %s' % (program_name, message, err))
177def fomo_wrapper_module(name):
178 try:
179 if not re.match(gf.meta.StringID.pattern, name):
180 raise ValueError('invalid name')
182 words = name.split('.', 1)
183 if len(words) == 2:
184 name, variant = words
185 else:
186 name = words[0]
187 variant = None
189 name_clean = re.sub(r'[.-]', '_', name)
190 modname = '.'.join(['pyrocko', 'fomosto', name_clean])
191 mod = __import__(modname, level=0)
192 return getattr(mod.fomosto, name_clean), variant
194 except ValueError:
195 die('invalid modelling code wrapper name')
197 except ImportError:
198 die('''modelling code wrapper "%s" not available or not installed
199 (module probed: "%s")''' % (name, modname))
202def command_init(args):
204 details = '''
206Available modelling backends:
207%s
209 More information at
210 https://pyrocko.org/docs/current/apps/fomosto/backends.html
211''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
213 parser, options, args = cl_parse(
214 'init', args,
215 details=details)
217 if len(args) == 0:
218 sys.exit(parser.format_help())
220 if args[0] == 'redeploy':
221 if len(args) != 3:
222 parser.error('incorrect number of arguments')
224 source_dir, dest_dir = args[1:]
226 try:
227 source = gf.Store(source_dir)
228 except gf.StoreError as e:
229 die(e)
231 config = copy.deepcopy(source.config)
232 config.derived_from_id = source.config.id
233 try:
234 config_filenames = gf.store.Store.create_editables(
235 dest_dir, config=config)
237 except gf.StoreError as e:
238 die(e)
240 try:
241 dest = gf.Store(dest_dir)
242 except gf.StoreError as e:
243 die(e)
245 for k in source.extra_keys():
246 source_fn = source.get_extra_path(k)
247 dest_fn = dest.get_extra_path(k)
248 shutil.copyfile(source_fn, dest_fn)
250 logger.info(
251 '(1) configure settings in files:\n %s'
252 % '\n '.join(config_filenames))
254 logger.info(
255 '(2) run "fomosto redeploy <source> <dest>", as needed')
257 else:
258 if len(args) != 2:
259 parser.error('incorrect number of arguments')
261 (modelling_code_id, store_dir) = args
263 module, variant = fomo_wrapper_module(modelling_code_id)
264 try:
265 config_filenames = module.init(store_dir, variant)
266 except gf.StoreError as e:
267 die(e)
269 logger.info('(1) configure settings in files:\n %s'
270 % '\n '.join(config_filenames))
271 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
272 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
275def get_store_dir(args):
276 if len(args) == 1:
277 store_dir = op.abspath(args.pop(0))
278 else:
279 store_dir = op.abspath(op.curdir)
281 if not op.isdir(store_dir):
282 die('not a directory: %s' % store_dir)
284 return store_dir
287def get_store_dirs(args):
288 if len(args) == 0:
289 store_dirs = [op.abspath(op.curdir)]
290 else:
291 store_dirs = [op.abspath(x) for x in args]
293 for store_dir in store_dirs:
294 if not op.isdir(store_dir):
295 die('not a directory: %s' % store_dir)
297 return store_dirs
300def command_build(args):
302 def setup(parser):
303 parser.add_option(
304 '--force', dest='force', action='store_true',
305 help='overwrite existing files')
307 parser.add_option(
308 '--nworkers', dest='nworkers', type='int', metavar='N',
309 help='run N worker processes in parallel')
311 parser.add_option(
312 '--continue', dest='continue_', action='store_true',
313 help='continue suspended build')
315 parser.add_option(
316 '--step', dest='step', type='int', metavar='I',
317 help='process block number IBLOCK')
319 parser.add_option(
320 '--block', dest='iblock', type='int', metavar='I',
321 help='process block number IBLOCK')
323 parser, options, args = cl_parse('build', args, setup=setup)
325 store_dir = get_store_dir(args)
326 try:
327 if options.step is not None:
328 step = options.step - 1
329 else:
330 step = None
332 if options.iblock is not None:
333 iblock = options.iblock - 1
334 else:
335 iblock = None
337 store = gf.Store(store_dir)
338 module, _ = fomo_wrapper_module(store.config.modelling_code_id)
339 module.build(
340 store_dir,
341 force=options.force,
342 nworkers=options.nworkers, continue_=options.continue_,
343 step=step,
344 iblock=iblock)
346 except gf.StoreError as e:
347 die(e)
350def command_stats(args):
352 parser, options, args = cl_parse('stats', args)
353 store_dir = get_store_dir(args)
355 try:
356 store = gf.Store(store_dir)
357 s = store.stats()
359 except gf.StoreError as e:
360 die(e)
362 for k in store.stats_keys:
363 print('%s: %s' % (k, s[k]))
366def command_check(args):
368 parser, options, args = cl_parse('check', args)
369 store_dir = get_store_dir(args)
371 try:
372 store = gf.Store(store_dir)
373 problems = store.check(show_progress=True)
374 if problems:
375 die('problems detected with gf store: %s' % store_dir)
377 except gf.StoreError as e:
378 die(e)
381def load_config(fn):
382 try:
383 config = gf.meta.load(filename=fn)
384 assert isinstance(config, gf.Config)
386 except Exception:
387 die('cannot load gf config from file: %s' % fn)
389 return config
392def command_decimate(args):
394 def setup(parser):
395 parser.add_option(
396 '--config', dest='config_fn', metavar='FILE',
397 help='use modified spacial sampling given in FILE')
399 parser.add_option(
400 '--force', dest='force', action='store_true',
401 help='overwrite existing files')
403 parser, options, args = cl_parse('decimate', args, setup=setup)
404 try:
405 decimate = int(args.pop())
406 except Exception:
407 parser.error('cannot get <factor> argument')
409 store_dir = get_store_dir(args)
411 config = None
412 if options.config_fn:
413 config = load_config(options.config_fn)
415 try:
416 store = gf.Store(store_dir)
417 store.make_decimated(decimate, config=config, force=options.force,
418 show_progress=True)
420 except gf.StoreError as e:
421 die(e)
424def sindex(args):
425 return '(%s)' % ', '.join('%g' % x for x in args)
428def command_redeploy(args):
430 parser, options, args = cl_parse('redeploy', args)
432 if not len(args) == 2:
433 sys.exit(parser.format_help())
435 source_store_dir, dest_store_dir = args
437 try:
438 source = gf.Store(source_store_dir)
439 except gf.StoreError as e:
440 die(e)
442 try:
443 gf.store.Store.create_dependants(dest_store_dir)
444 except gf.StoreError:
445 pass
447 try:
448 dest = gf.Store(dest_store_dir, 'w')
450 except gf.StoreError as e:
451 die(e)
453 show_progress = True
455 try:
456 if show_progress:
457 pbar = util.progressbar('redeploying', dest.config.nrecords)
459 for i, args in enumerate(dest.config.iter_nodes()):
460 try:
461 tr = source.get(args, interpolation='off')
462 dest.put(args, tr)
464 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa
465 logger.debug('skipping %s, (%s)' % (sindex(args), e))
467 except gf.store.StoreError as e:
468 logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
470 if show_progress:
471 pbar.update(i+1)
473 finally:
474 if show_progress:
475 pbar.finish()
478def command_trans2rot(args):
479 from pyrocko.fomosto.elastic10_to_rotational8 \
480 import elastic10_to_rotational8
482 details = '''
484Create a rotational GF store from a translational one using finite differences
485approximations (elastic10 or elastic10_fd to rotational8).
487If the source store is of type A and receivers at the surface, the free
488surface condition is used to replace the vertical derivatives with horizontal
489ones. If the input is of type B, the complete set of derivatives are used. It
490is important that the spatial sampling of the source store is high enough so
491that the finite difference approximations are sufficiently accurate.
493If the destination store directory does not exist, it is created and the extent
494and spacing of the source store is used for the rotational store. If it does
495exist and contains a valid `config`, it is used in place in order to allow for
496different input and output extents and spacings. The store must be empty in the
497latter case (delete the `index` and `traces` files if needed).
498'''
500 parser, options, args = cl_parse(
501 'trans2rot', args, details=details)
503 if not len(args) == 2:
504 sys.exit(parser.format_help())
506 source_store_dir, dest_store_dir = args
508 try:
509 elastic10_to_rotational8(source_store_dir, dest_store_dir)
510 except gf.store.StoreError as e:
511 die(e)
514def command_fd2trans(args):
515 from pyrocko.fomosto.elastic10_fd_to_elastic10 \
516 import elastic10_fd_to_elastic10
518 details = '''
520Extract regular GF store (elastic10) from a finite difference support GF store
521(elastic10_fd).
522'''
524 parser, options, args = cl_parse(
525 'fd2trans', args, details=details)
527 if not len(args) == 2:
528 sys.exit(parser.format_help())
530 source_store_dir, dest_store_dir = args
532 try:
533 elastic10_fd_to_elastic10(source_store_dir, dest_store_dir)
534 except gf.store.StoreError as e:
535 die(e)
538def command_view(args):
539 def setup(parser):
540 parser.add_option(
541 '--extract',
542 dest='extract',
543 metavar='start:stop[:step|@num],...',
544 help='specify which traces to show')
546 parser.add_option(
547 '--show-phases',
548 dest='showphases',
549 default=None,
550 metavar='[phase_id1,phase_id2,...|all]',
551 help='add phase markers from ttt')
553 parser.add_option(
554 '--opengl',
555 dest='opengl',
556 action='store_true',
557 default=None,
558 help='use OpenGL for drawing')
560 parser.add_option(
561 '--no-opengl',
562 dest='opengl',
563 action='store_false',
564 default=None,
565 help='do not use OpenGL for drawing')
567 parser, options, args = cl_parse('view', args, setup=setup)
569 gdef = None
570 if options.extract:
571 try:
572 gdef = gf.meta.parse_grid_spec(options.extract)
573 except gf.meta.GridSpecError as e:
574 die(e)
576 store_dirs = get_store_dirs(args)
578 alpha = 'abcdefghijklmnopqrstxyz'.upper()
580 markers = []
581 traces = []
583 try:
584 for istore, store_dir in enumerate(store_dirs):
585 store = gf.Store(store_dir)
586 if options.showphases == 'all':
587 phasenames = [pn.id for pn in store.config.tabulated_phases]
588 elif options.showphases is not None:
589 phasenames = options.showphases.split(',')
591 ii = 0
592 for args in store.config.iter_extraction(gdef):
593 gtr = store.get(args)
595 loc_code = ''
596 if len(store_dirs) > 1:
597 loc_code = alpha[istore % len(alpha)]
599 if gtr:
601 sta_code = '%04i (%s)' % (
602 ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
603 tmin = gtr.deltat * gtr.itmin
605 tr = trace.Trace(
606 '',
607 sta_code,
608 loc_code,
609 '%02i' % args[-1],
610 ydata=gtr.data,
611 deltat=gtr.deltat,
612 tmin=tmin)
614 if options.showphases:
615 for phasename in phasenames:
616 phase_tmin = store.t(phasename, args[:-1])
617 if phase_tmin:
618 m = marker.PhaseMarker(
619 [('',
620 sta_code,
621 loc_code,
622 '%02i' % args[-1])],
623 phase_tmin,
624 phase_tmin,
625 0,
626 phasename=phasename)
627 markers.append(m)
629 traces.append(tr)
631 ii += 1
633 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
634 die(e)
636 def setup(win):
637 viewer = win.get_view()
638 viewer.menuitem_showboxes.setChecked(False)
639 viewer.menuitem_colortraces.setChecked(True)
640 viewer.menuitem_demean.setChecked(False)
641 viewer.menuitems_sorting[5][0].setChecked(True)
642 viewer.menuitems_scaling[2][0].setChecked(True)
643 viewer.sortingmode_change()
644 viewer.scalingmode_change()
646 trace.snuffle(
647 traces, markers=markers, opengl=options.opengl, launch_hook=setup)
650def command_extract(args):
651 def setup(parser):
652 parser.add_option(
653 '--format', dest='format', default='mseed',
654 choices=['mseed', 'sac', 'text', 'yaff'],
655 help='export to format "mseed", "sac", "text", or "yaff". '
656 'Default is "mseed".')
658 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
659 parser.add_option(
660 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
661 help='output path template [default: "%s"]' % fndfl)
663 parser, options, args = cl_parse('extract', args, setup=setup)
664 try:
665 sdef = args.pop()
666 except Exception:
667 parser.error('cannot get <selection> argument')
669 try:
670 gdef = gf.meta.parse_grid_spec(sdef)
671 except gf.meta.GridSpecError as e:
672 die(e)
674 store_dir = get_store_dir(args)
676 extensions = {
677 'mseed': 'mseed',
678 'sac': 'sac',
679 'text': 'txt',
680 'yaff': 'yaff'}
682 try:
683 store = gf.Store(store_dir)
684 for args in store.config.iter_extraction(gdef):
685 gtr = store.get(args)
686 if gtr:
687 tr = trace.Trace(
688 '', '', '', util.zfmt(store.config.ncomponents) % args[-1],
689 ydata=gtr.data,
690 deltat=gtr.deltat,
691 tmin=gtr.deltat * gtr.itmin)
693 additional = dict(
694 args='_'.join('%g' % x for x in args),
695 irecord=store.str_irecord(args),
696 extension=extensions[options.format])
698 io.save(
699 tr,
700 options.output_fn,
701 format=options.format,
702 additional=additional)
704 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
705 die(e)
708def command_import(args):
709 try:
710 from tunguska import gfdb
711 except ImportError:
712 die('the kiwi tools must be installed to use this feature')
714 parser, options, args = cl_parse('import', args)
716 show_progress = True
718 if not len(args) == 2:
719 sys.exit(parser.format_help())
721 source_path, dest_store_dir = args
723 if op.isdir(source_path):
724 source_path = op.join(source_path, 'db')
726 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
728 db = gfdb.Gfdb(source_path)
730 config = gf.meta.ConfigTypeA(
731 id='imported_gfs',
732 distance_min=db.firstx,
733 distance_max=db.firstx + (db.nx-1) * db.dx,
734 distance_delta=db.dx,
735 source_depth_min=db.firstz,
736 source_depth_max=db.firstz + (db.nz-1) * db.dz,
737 source_depth_delta=db.dz,
738 sample_rate=1.0/db.dt,
739 ncomponents=db.ng
740 )
742 try:
743 gf.store.Store.create(dest_store_dir, config=config)
744 dest = gf.Store(dest_store_dir, 'w')
745 try:
746 if show_progress:
747 pbar = util.progressbar(
748 'importing', dest.config.nrecords/dest.config.ncomponents)
750 for i, args in enumerate(dest.config.iter_nodes(level=-1)):
751 source_depth, distance = [float(x) for x in args]
752 traces = db.get_traces_pyrocko(distance, source_depth)
753 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
754 for ig in range(db.ng):
755 if ig in ig_to_trace:
756 tr = ig_to_trace[ig]
757 gf_tr = gf.store.GFTrace(
758 tr.get_ydata(),
759 int(round(tr.tmin / tr.deltat)),
760 tr.deltat)
762 else:
763 gf_tr = gf.store.Zero
765 dest.put((source_depth, distance, ig), gf_tr)
767 if show_progress:
768 pbar.update(i+1)
770 finally:
771 if show_progress:
772 pbar.finish()
774 dest.close()
776 except gf.StoreError as e:
777 die(e)
780def command_export(args):
781 from subprocess import Popen, PIPE
783 try:
784 from tunguska import gfdb
785 except ImportError as err:
786 die('the kiwi tools must be installed to use this feature', err)
788 def setup(parser):
789 parser.add_option(
790 '--nchunks', dest='nchunks', type='int', default=1, metavar='N',
791 help='split output gfdb into N chunks')
793 parser, options, args = cl_parse('export', args, setup=setup)
795 show_progress = True
797 if len(args) not in (1, 2):
798 sys.exit(parser.format_help())
800 target_path = args.pop()
801 if op.isdir(target_path):
802 target_path = op.join(target_path, 'kiwi_gfdb')
803 logger.warning('exported gfdb will be named as "%s.*"' % target_path)
805 source_store_dir = get_store_dir(args)
807 source = gf.Store(source_store_dir, 'r')
808 config = source.config
810 if not isinstance(config, gf.meta.ConfigTypeA):
811 die('only stores of type A can be exported to Kiwi format')
813 if op.isfile(target_path + '.index'):
814 die('destation already exists')
816 cmd = [str(x) for x in [
817 'gfdb_build',
818 target_path,
819 options.nchunks,
820 config.ndistances,
821 config.nsource_depths,
822 config.ncomponents,
823 config.deltat,
824 config.distance_delta,
825 config.source_depth_delta,
826 config.distance_min,
827 config.source_depth_min]]
829 p = Popen(cmd, stdin=PIPE)
830 p.communicate()
832 out_db = gfdb.Gfdb(target_path)
834 try:
835 if show_progress:
836 pbar = util.progressbar(
837 'exporting', config.nrecords/config.ncomponents)
839 for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
841 data_out = []
842 for ig in range(config.ncomponents):
843 try:
844 tr = source.get((z, x, ig), interpolation='off')
845 data_out.append((tr.t, tr.data * config.factor))
847 except gf.store.StoreError as e:
848 logger.warning(
849 'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
850 data_out.append(None)
852 # put a zero valued sample to no-data zero-traces at a compatible
853 # time
854 tmins = [
855 entry[0][0]
856 for entry in data_out
857 if entry is not None and entry[0].size != 0]
859 if tmins:
860 tmin = min(tmins)
861 for entry in data_out:
862 if entry is not None and entry[0].size == 0:
863 entry[0].resize(1)
864 entry[1].resize(1)
865 entry[0][0] = tmin
866 entry[1][0] = 0.0
868 out_db.put_traces_slow(x, z, data_out)
870 if show_progress:
871 pbar.update(i+1)
873 source.close()
875 finally:
876 if show_progress:
877 pbar.finish()
880def phasedef_or_horvel(x):
881 try:
882 return float(x)
883 except ValueError:
884 return cake.PhaseDef(x)
887def mkp(s):
888 return [phasedef_or_horvel(ps) for ps in s.split(',')]
891def stored_attribute_table_plots(phase_ids, options, args, attribute):
892 import numpy as num
893 from pyrocko.plot.cake_plot import xscaled, yscaled, mpl_init
894 from pyrocko.plot import mpl_labelspace as labelspace
896 plt = mpl_init()
898 np = 1
899 store_dir = get_store_dir(args)
900 for phase_id in phase_ids:
901 try:
902 store = gf.Store(store_dir)
904 if options.receiver_depth is not None:
905 receiver_depth = options.receiver_depth * 1000.0
906 else:
907 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
908 receiver_depth = store.config.receiver_depth
910 elif isinstance(store.config, gf.ConfigTypeB):
911 receiver_depth = store.config.receiver_depth_min
913 else:
914 receiver_depth = 0.0
916 phase = store.get_stored_phase(phase_id, attribute)
917 axes = plt.subplot(2, len(phase_ids), np)
918 labelspace(axes)
919 xscaled(1. / km, axes)
920 yscaled(1. / km, axes)
921 x = None
922 if isinstance(store.config, gf.ConfigTypeB):
923 x = (receiver_depth, None, None)
925 phase.plot_2d(axes, x=x)
926 axes.set_title(phase_id)
927 np += 1
928 except gf.StoreError as e:
929 die(e)
931 axes = plt.subplot(2, 1, 2)
932 num_d = 100
933 distances = num.linspace(store.config.distance_min,
934 store.config.distance_max,
935 num_d)
937 if options.source_depth is not None:
938 source_depth = options.source_depth * km
939 else:
940 source_depth = store.config.source_depth_min + (
941 store.config.source_depth_max - store.config.source_depth_min) / 2.
943 if isinstance(store.config, gf.ConfigTypeA):
944 attribute_vals = num.empty(num_d)
945 for phase_id in phase_ids:
946 attribute_vals[:] = num.nan
947 for i, d in enumerate(distances):
948 if attribute == 'phase':
949 attribute_vals[i] = store.t(phase_id, (source_depth, d))
950 ylabel = 'TT [s]'
951 else:
952 attribute_vals[i] = store.get_stored_attribute(
953 phase_id, options.attribute, (source_depth, d))
954 ylabel = '%s [deg]' % options.attribute
956 axes.plot(distances / km, attribute_vals, label=phase_id)
958 axes.set_title('source source_depth %s km' % (source_depth / km))
959 axes.set_xlabel('distance [km]')
960 axes.set_ylabel(ylabel)
961 axes.legend()
963 plt.tight_layout()
964 mpl_show(plt)
967def command_ttt(args):
968 def setup(parser):
969 parser.add_option(
970 '--force', dest='force', action='store_true',
971 help='overwrite existing files')
973 parser, options, args = cl_parse('ttt', args, setup=setup)
975 store_dir = get_store_dir(args)
976 try:
977 store = gf.Store(store_dir)
978 store.make_travel_time_tables(force=options.force)
980 except gf.StoreError as e:
981 die(e)
984def command_tttview(args):
986 def setup(parser):
987 parser.add_option(
988 '--source-depth', dest='source_depth', type=float,
989 help='Source depth in km')
991 parser.add_option(
992 '--receiver-depth', dest='receiver_depth', type=float,
993 help='Receiver depth in km')
995 parser, options, args = cl_parse(
996 'tttview', args, setup=setup,
997 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
999 try:
1000 phase_ids = args.pop().split(',')
1001 except Exception:
1002 parser.error('cannot get <phase-ids> argument')
1004 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
1007def command_sat(args):
1008 def setup(parser):
1009 parser.add_option(
1010 '--force', dest='force', action='store_true',
1011 help='overwrite existing files')
1013 parser.add_option(
1014 '--attribute',
1015 action='store',
1016 dest='attribute',
1017 type='choice',
1018 choices=gf.store.available_stored_tables[1::],
1019 default='takeoff_angle',
1020 help='calculate interpolation table for selected ray attributes.')
1022 parser, options, args = cl_parse('sat', args, setup=setup)
1024 store_dir = get_store_dir(args)
1025 try:
1026 store = gf.Store(store_dir)
1027 store.make_stored_table(options.attribute, force=options.force)
1029 except gf.StoreError as e:
1030 die(e)
1033def command_satview(args):
1035 def setup(parser):
1036 parser.add_option(
1037 '--source-depth', dest='source_depth', type=float,
1038 help='Source depth in km')
1040 parser.add_option(
1041 '--receiver-depth', dest='receiver_depth', type=float,
1042 help='Receiver depth in km')
1044 parser.add_option(
1045 '--attribute',
1046 action='store',
1047 dest='attribute',
1048 type='choice',
1049 choices=gf.store.available_stored_tables[1::],
1050 default='takeoff_angle',
1051 help='view selected ray attribute.')
1053 parser, options, args = cl_parse(
1054 'satview', args, setup=setup,
1055 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
1057 try:
1058 phase_ids = args.pop().split(',')
1059 except Exception:
1060 parser.error('cannot get <phase-ids> argument')
1062 logger.info('Plotting stored attribute %s' % options.attribute)
1064 stored_attribute_table_plots(
1065 phase_ids, options, args, attribute=options.attribute)
1068def command_tttextract(args):
1069 def setup(parser):
1070 parser.add_option(
1071 '--output', dest='output_fn', metavar='TEMPLATE',
1072 help='output to text files instead of stdout '
1073 '(example TEMPLATE: "extracted/%(args)s.txt")')
1075 parser, options, args = cl_parse('tttextract', args, setup=setup)
1076 try:
1077 sdef = args.pop()
1078 except Exception:
1079 parser.error('cannot get <selection> argument')
1081 try:
1082 sphase = args.pop()
1083 except Exception:
1084 parser.error('cannot get <phase> argument')
1086 try:
1087 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1088 except gf.meta.InvalidTimingSpecification:
1089 parser.error('invalid phase specification: "%s"' % sphase)
1091 try:
1092 gdef = gf.meta.parse_grid_spec(sdef)
1093 except gf.meta.GridSpecError as e:
1094 die(e)
1096 store_dir = get_store_dir(args)
1098 try:
1099 store = gf.Store(store_dir)
1100 for args in store.config.iter_extraction(gdef, level=-1):
1101 s = ['%e' % x for x in args]
1102 for phase in phases:
1103 t = store.t(phase, args)
1104 if t is not None:
1105 s.append('%e' % t)
1106 else:
1107 s.append('nan')
1109 if options.output_fn:
1110 d = dict(
1111 args='_'.join('%e' % x for x in args),
1112 extension='txt')
1114 fn = options.output_fn % d
1115 util.ensuredirs(fn)
1116 with open(fn, 'a') as f:
1117 f.write(' '.join(s))
1118 f.write('\n')
1119 else:
1120 print(' '.join(s))
1122 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1123 die(e)
1126def command_tttlsd(args):
1128 def setup(parser):
1129 pass
1131 parser, options, args = cl_parse(
1132 'tttlsd', args,
1133 setup=setup,
1134 details='''
1136This subcommand fills holes in travel-time-tables by using an eikonal solver to
1137predict travel-times for the missing values. The approach works best for simple
1138P or S phases without any reflections or conversions. It creates new
1139travel-time-tables which are named by adding the suffix `.lsd` to the original
1140name. E.g. running `fomosto tttlsd begin` would produce a new travel-time-table
1141`begin.lsd`. The new name can be referenced anywhere where a stored
1142travel-time-table can be used, e.g. in `extra/qseis`.
1143''')
1145 try:
1146 sphase_ids = args.pop()
1147 except Exception:
1148 parser.error('cannot get <phase> argument')
1150 try:
1151 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1152 except gf.meta.InvalidTimingSpecification:
1153 parser.error('invalid phase specification: "%s"' % sphase_ids)
1155 store_dir = get_store_dir(args)
1157 try:
1158 store = gf.Store(store_dir)
1159 for phase_id in phase_ids:
1160 store.fix_ttt_holes(phase_id)
1162 except gf.StoreError as e:
1163 die(e)
1166def command_server(args):
1167 from pyrocko.gf import server
1169 def setup(parser):
1170 parser.add_option(
1171 '--port', dest='port', metavar='PORT', type='int', default=8080,
1172 help='serve on port PORT')
1174 parser.add_option(
1175 '--ip', dest='ip', metavar='IP', default='',
1176 help='serve on ip address IP')
1178 parser, options, args = cl_parse('server', args, setup=setup)
1180 engine = gf.LocalEngine(store_superdirs=args)
1181 server.run(options.ip, options.port, engine)
1184def command_download(args):
1185 from pyrocko.gf import ws
1187 details = '''
1189Browse pre-calculated Green's function stores online:
1191 https://greens-mill.pyrocko.org
1192'''
1194 def setup(parser):
1195 parser.add_option(
1196 '--force', dest='force', action='store_true',
1197 help='overwrite existing files')
1199 parser, options, args = cl_parse(
1200 'download', args, setup=setup, details=details)
1201 if len(args) not in (1, 2):
1202 sys.exit(parser.format_help())
1204 if len(args) == 2:
1205 site, store_id = args
1206 if not re.match(gf.meta.StringID.pattern, store_id):
1207 die('invalid store ID')
1208 else:
1209 site, store_id = args[0], None
1211 if site not in gf.ws.g_site_abbr:
1212 if -1 == site.find('://'):
1213 site = 'http://' + site
1215 try:
1216 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1217 except ws.DownloadError as e:
1218 die(str(e))
1221def command_modelview(args):
1223 import matplotlib.pyplot as plt
1224 import numpy as num
1225 from pyrocko.plot.cake_plot import mpl_init, xscaled, yscaled
1226 from pyrocko.plot import mpl_labelspace as labelspace
1227 mpl_init()
1229 neat_labels = {
1230 'vp': '$v_p$',
1231 'vs': '$v_s$',
1232 'qp': '$Q_p$',
1233 'qs': '$Q_s$',
1234 'rho': '$\\rho$'}
1236 def setup(parser):
1237 parser.add_option(
1238 '--parameters', dest='parameters',
1239 default='vp,vs', metavar='vp,vs,....',
1240 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1242 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1244 parser, options, args = cl_parse('modelview', args, setup=setup)
1246 store_dirs = get_store_dirs(args)
1248 parameters = options.parameters.split(',')
1250 fig, axes = plt.subplots(1,
1251 len(parameters),
1252 sharey=True,
1253 figsize=(len(parameters)*3, 5))
1255 if not isinstance(axes, num.ndarray):
1256 axes = [axes]
1258 axes = dict(zip(parameters, axes))
1260 for store_id in store_dirs:
1261 try:
1262 store = gf.Store(store_id)
1263 mod = store.config.earthmodel_1d
1264 z = mod.profile('z')
1266 for p in parameters:
1267 ax = axes[p]
1268 labelspace(ax)
1269 if '/' in p:
1270 p1, p2 = p.split('/')
1271 profile = mod.profile(p1)/mod.profile(p2)
1272 else:
1273 profile = mod.profile(p)
1275 ax.plot(profile, -z, label=store_id.split('/')[-1])
1276 if p in ['vs', 'vp', 'rho']:
1277 xscaled(1./1000., ax)
1279 yscaled(1./km, ax)
1281 except gf.StoreError as e:
1282 die(e)
1284 for p, ax in axes.items():
1285 ax.grid()
1286 if p in neat_labels:
1287 lab = neat_labels[p]
1288 elif all(x in neat_labels for x in p.split('/')):
1289 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1290 else:
1291 lab = p
1293 ax.set_title(lab, y=1.02)
1295 if p in units:
1296 lab += ' ' + units[p]
1298 ax.autoscale()
1299 ax.set_xlabel(lab)
1301 axes[parameters[0]].set_ylabel('Depth [km]')
1303 handles, labels = ax.get_legend_handles_labels()
1305 if len(store_dirs) > 1:
1306 ax.legend(
1307 handles,
1308 labels,
1309 bbox_to_anchor=(0.5, 0.12),
1310 bbox_transform=fig.transFigure,
1311 ncol=3,
1312 loc='upper center',
1313 fancybox=True)
1315 plt.subplots_adjust(bottom=0.22,
1316 wspace=0.05)
1318 mpl_show(plt)
1321def command_upgrade(args):
1322 parser, options, args = cl_parse('upgrade', args)
1323 store_dirs = get_store_dirs(args)
1324 try:
1325 for store_dir in store_dirs:
1326 store = gf.Store(store_dir)
1327 nup = store.upgrade()
1328 if nup == 0:
1329 print('%s: already up-to-date.' % store_dir)
1330 else:
1331 print('%s: %i file%s upgraded.' % (
1332 store_dir, nup, ['s', ''][nup == 1]))
1334 except gf.StoreError as e:
1335 die(e)
1338def command_addref(args):
1339 parser, options, args = cl_parse('addref', args)
1341 store_dirs = []
1342 filenames = []
1343 for arg in args:
1344 if op.isdir(arg):
1345 store_dirs.append(arg)
1346 elif op.isfile(arg):
1347 filenames.append(arg)
1348 else:
1349 die('invalid argument: %s' % arg)
1351 if not store_dirs:
1352 store_dirs.append('.')
1354 references = []
1355 try:
1356 for filename in args:
1357 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1358 except ImportError:
1359 die('pybtex module must be installed to use this function')
1361 if not references:
1362 die('no references found')
1364 for store_dir in store_dirs:
1365 try:
1366 store = gf.Store(store_dir)
1367 ids = [ref.id for ref in store.config.references]
1368 for ref in references:
1369 if ref.id in ids:
1370 die('duplicate reference id: %s' % ref.id)
1372 ids.append(ref.id)
1373 store.config.references.append(ref)
1375 store.save_config(make_backup=True)
1377 except gf.StoreError as e:
1378 die(e)
1381def command_qc(args):
1382 parser, options, args = cl_parse('qc', args)
1384 store_dir = get_store_dir(args)
1386 try:
1387 store = gf.Store(store_dir)
1388 s = store.stats()
1389 if s['empty'] != 0:
1390 print('has empty records')
1392 for aname in ['author', 'author_email', 'description', 'references']:
1394 if not getattr(store.config, aname):
1395 print('%s empty' % aname)
1397 except gf.StoreError as e:
1398 die(e)
1401def command_report(args):
1402 from pyrocko.fomosto.report import report_call
1403 report_call.run_program(args)
1406def main(args=None):
1407 '''
1408 CLI entry point for Pyrocko's ``fomosto`` app.
1409 '''
1410 if args is None:
1411 args = sys.argv[1:]
1413 if len(args) < 1:
1414 sys.exit('Usage: %s' % usage)
1416 command = args.pop(0)
1418 if command in subcommands:
1419 globals()['command_' + command](args)
1421 elif command in ('--help', '-h', 'help'):
1422 if command == 'help' and args:
1423 acommand = args[0]
1424 if acommand in subcommands:
1425 globals()['command_' + acommand](['--help'])
1427 sys.exit('Usage: %s' % usage)
1429 else:
1430 sys.exit('fomosto: error: no such subcommand: %s' % command)
1433if __name__ == '__main__':
1434 main()