Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/apps/fomosto.py: 59%
727 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-03-20 14:04 +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 labelspace, xscaled, yscaled, mpl_init
895 plt = mpl_init()
897 np = 1
898 store_dir = get_store_dir(args)
899 for phase_id in phase_ids:
900 try:
901 store = gf.Store(store_dir)
903 if options.receiver_depth is not None:
904 receiver_depth = options.receiver_depth * 1000.0
905 else:
906 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
907 receiver_depth = store.config.receiver_depth
909 elif isinstance(store.config, gf.ConfigTypeB):
910 receiver_depth = store.config.receiver_depth_min
912 else:
913 receiver_depth = 0.0
915 phase = store.get_stored_phase(phase_id, attribute)
916 axes = plt.subplot(2, len(phase_ids), np)
917 labelspace(axes)
918 xscaled(1. / km, axes)
919 yscaled(1. / km, axes)
920 x = None
921 if isinstance(store.config, gf.ConfigTypeB):
922 x = (receiver_depth, None, None)
924 phase.plot_2d(axes, x=x)
925 axes.set_title(phase_id)
926 np += 1
927 except gf.StoreError as e:
928 die(e)
930 axes = plt.subplot(2, 1, 2)
931 num_d = 100
932 distances = num.linspace(store.config.distance_min,
933 store.config.distance_max,
934 num_d)
936 if options.source_depth is not None:
937 source_depth = options.source_depth * km
938 else:
939 source_depth = store.config.source_depth_min + (
940 store.config.source_depth_max - store.config.source_depth_min) / 2.
942 if isinstance(store.config, gf.ConfigTypeA):
943 attribute_vals = num.empty(num_d)
944 for phase_id in phase_ids:
945 attribute_vals[:] = num.NAN
946 for i, d in enumerate(distances):
947 if attribute == 'phase':
948 attribute_vals[i] = store.t(phase_id, (source_depth, d))
949 ylabel = 'TT [s]'
950 else:
951 attribute_vals[i] = store.get_stored_attribute(
952 phase_id, options.attribute, (source_depth, d))
953 ylabel = '%s [deg]' % options.attribute
955 axes.plot(distances / km, attribute_vals, label=phase_id)
957 axes.set_title('source source_depth %s km' % (source_depth / km))
958 axes.set_xlabel('distance [km]')
959 axes.set_ylabel(ylabel)
960 axes.legend()
962 plt.tight_layout()
963 mpl_show(plt)
966def command_ttt(args):
967 def setup(parser):
968 parser.add_option(
969 '--force', dest='force', action='store_true',
970 help='overwrite existing files')
972 parser, options, args = cl_parse('ttt', args, setup=setup)
974 store_dir = get_store_dir(args)
975 try:
976 store = gf.Store(store_dir)
977 store.make_travel_time_tables(force=options.force)
979 except gf.StoreError as e:
980 die(e)
983def command_tttview(args):
985 def setup(parser):
986 parser.add_option(
987 '--source-depth', dest='source_depth', type=float,
988 help='Source depth in km')
990 parser.add_option(
991 '--receiver-depth', dest='receiver_depth', type=float,
992 help='Receiver depth in km')
994 parser, options, args = cl_parse(
995 'tttview', args, setup=setup,
996 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
998 try:
999 phase_ids = args.pop().split(',')
1000 except Exception:
1001 parser.error('cannot get <phase-ids> argument')
1003 stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
1006def command_sat(args):
1007 def setup(parser):
1008 parser.add_option(
1009 '--force', dest='force', action='store_true',
1010 help='overwrite existing files')
1012 parser.add_option(
1013 '--attribute',
1014 action='store',
1015 dest='attribute',
1016 type='choice',
1017 choices=gf.store.available_stored_tables[1::],
1018 default='takeoff_angle',
1019 help='calculate interpolation table for selected ray attributes.')
1021 parser, options, args = cl_parse('sat', args, setup=setup)
1023 store_dir = get_store_dir(args)
1024 try:
1025 store = gf.Store(store_dir)
1026 store.make_stored_table(options.attribute, force=options.force)
1028 except gf.StoreError as e:
1029 die(e)
1032def command_satview(args):
1034 def setup(parser):
1035 parser.add_option(
1036 '--source-depth', dest='source_depth', type=float,
1037 help='Source depth in km')
1039 parser.add_option(
1040 '--receiver-depth', dest='receiver_depth', type=float,
1041 help='Receiver depth in km')
1043 parser.add_option(
1044 '--attribute',
1045 action='store',
1046 dest='attribute',
1047 type='choice',
1048 choices=gf.store.available_stored_tables[1::],
1049 default='takeoff_angle',
1050 help='view selected ray attribute.')
1052 parser, options, args = cl_parse(
1053 'satview', args, setup=setup,
1054 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
1056 try:
1057 phase_ids = args.pop().split(',')
1058 except Exception:
1059 parser.error('cannot get <phase-ids> argument')
1061 logger.info('Plotting stored attribute %s' % options.attribute)
1063 stored_attribute_table_plots(
1064 phase_ids, options, args, attribute=options.attribute)
1067def command_tttextract(args):
1068 def setup(parser):
1069 parser.add_option(
1070 '--output', dest='output_fn', metavar='TEMPLATE',
1071 help='output to text files instead of stdout '
1072 '(example TEMPLATE: "extracted/%(args)s.txt")')
1074 parser, options, args = cl_parse('tttextract', args, setup=setup)
1075 try:
1076 sdef = args.pop()
1077 except Exception:
1078 parser.error('cannot get <selection> argument')
1080 try:
1081 sphase = args.pop()
1082 except Exception:
1083 parser.error('cannot get <phase> argument')
1085 try:
1086 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
1087 except gf.meta.InvalidTimingSpecification:
1088 parser.error('invalid phase specification: "%s"' % sphase)
1090 try:
1091 gdef = gf.meta.parse_grid_spec(sdef)
1092 except gf.meta.GridSpecError as e:
1093 die(e)
1095 store_dir = get_store_dir(args)
1097 try:
1098 store = gf.Store(store_dir)
1099 for args in store.config.iter_extraction(gdef, level=-1):
1100 s = ['%e' % x for x in args]
1101 for phase in phases:
1102 t = store.t(phase, args)
1103 if t is not None:
1104 s.append('%e' % t)
1105 else:
1106 s.append('nan')
1108 if options.output_fn:
1109 d = dict(
1110 args='_'.join('%e' % x for x in args),
1111 extension='txt')
1113 fn = options.output_fn % d
1114 util.ensuredirs(fn)
1115 with open(fn, 'a') as f:
1116 f.write(' '.join(s))
1117 f.write('\n')
1118 else:
1119 print(' '.join(s))
1121 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
1122 die(e)
1125def command_tttlsd(args):
1127 def setup(parser):
1128 pass
1130 parser, options, args = cl_parse(
1131 'tttlsd', args,
1132 setup=setup,
1133 details='''
1135This subcommand fills holes in travel-time-tables by using an eikonal solver to
1136predict travel-times for the missing values. The approach works best for simple
1137P or S phases without any reflections or conversions. It creates new
1138travel-time-tables which are named by adding the suffix `.lsd` to the original
1139name. E.g. running `fomosto tttlsd begin` would produce a new travel-time-table
1140`begin.lsd`. The new name can be referenced anywhere where a stored
1141travel-time-table can be used, e.g. in `extra/qseis`.
1142''')
1144 try:
1145 sphase_ids = args.pop()
1146 except Exception:
1147 parser.error('cannot get <phase> argument')
1149 try:
1150 phase_ids = [x.strip() for x in sphase_ids.split(',')]
1151 except gf.meta.InvalidTimingSpecification:
1152 parser.error('invalid phase specification: "%s"' % sphase_ids)
1154 store_dir = get_store_dir(args)
1156 try:
1157 store = gf.Store(store_dir)
1158 for phase_id in phase_ids:
1159 store.fix_ttt_holes(phase_id)
1161 except gf.StoreError as e:
1162 die(e)
1165def command_server(args):
1166 from pyrocko.gf import server
1168 def setup(parser):
1169 parser.add_option(
1170 '--port', dest='port', metavar='PORT', type='int', default=8080,
1171 help='serve on port PORT')
1173 parser.add_option(
1174 '--ip', dest='ip', metavar='IP', default='',
1175 help='serve on ip address IP')
1177 parser, options, args = cl_parse('server', args, setup=setup)
1179 engine = gf.LocalEngine(store_superdirs=args)
1180 server.run(options.ip, options.port, engine)
1183def command_download(args):
1184 from pyrocko.gf import ws
1186 details = '''
1188Browse pre-calculated Green's function stores online:
1190 https://greens-mill.pyrocko.org
1191'''
1193 def setup(parser):
1194 parser.add_option(
1195 '--force', dest='force', action='store_true',
1196 help='overwrite existing files')
1198 parser, options, args = cl_parse(
1199 'download', args, setup=setup, details=details)
1200 if len(args) not in (1, 2):
1201 sys.exit(parser.format_help())
1203 if len(args) == 2:
1204 site, store_id = args
1205 if not re.match(gf.meta.StringID.pattern, store_id):
1206 die('invalid store ID')
1207 else:
1208 site, store_id = args[0], None
1210 if site not in gf.ws.g_site_abbr:
1211 if -1 == site.find('://'):
1212 site = 'http://' + site
1214 try:
1215 ws.download_gf_store(site=site, store_id=store_id, force=options.force)
1216 except ws.DownloadError as e:
1217 die(str(e))
1220def command_modelview(args):
1222 import matplotlib.pyplot as plt
1223 import numpy as num
1224 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
1225 mpl_init()
1227 neat_labels = {
1228 'vp': '$v_p$',
1229 'vs': '$v_s$',
1230 'qp': '$Q_p$',
1231 'qs': '$Q_s$',
1232 'rho': '$\\rho$'}
1234 def setup(parser):
1235 parser.add_option(
1236 '--parameters', dest='parameters',
1237 default='vp,vs', metavar='vp,vs,....',
1238 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
1240 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
1242 parser, options, args = cl_parse('modelview', args, setup=setup)
1244 store_dirs = get_store_dirs(args)
1246 parameters = options.parameters.split(',')
1248 fig, axes = plt.subplots(1,
1249 len(parameters),
1250 sharey=True,
1251 figsize=(len(parameters)*3, 5))
1253 if not isinstance(axes, num.ndarray):
1254 axes = [axes]
1256 axes = dict(zip(parameters, axes))
1258 for store_id in store_dirs:
1259 try:
1260 store = gf.Store(store_id)
1261 mod = store.config.earthmodel_1d
1262 z = mod.profile('z')
1264 for p in parameters:
1265 ax = axes[p]
1266 labelspace(ax)
1267 if '/' in p:
1268 p1, p2 = p.split('/')
1269 profile = mod.profile(p1)/mod.profile(p2)
1270 else:
1271 profile = mod.profile(p)
1273 ax.plot(profile, -z, label=store_id.split('/')[-1])
1274 if p in ['vs', 'vp', 'rho']:
1275 xscaled(1./1000., ax)
1277 yscaled(1./km, ax)
1279 except gf.StoreError as e:
1280 die(e)
1282 for p, ax in axes.items():
1283 ax.grid()
1284 if p in neat_labels:
1285 lab = neat_labels[p]
1286 elif all(x in neat_labels for x in p.split('/')):
1287 lab = '/'.join(neat_labels[x] for x in p.split('/'))
1288 else:
1289 lab = p
1291 ax.set_title(lab, y=1.02)
1293 if p in units:
1294 lab += ' ' + units[p]
1296 ax.autoscale()
1297 ax.set_xlabel(lab)
1299 axes[parameters[0]].set_ylabel('Depth [km]')
1301 handles, labels = ax.get_legend_handles_labels()
1303 if len(store_dirs) > 1:
1304 ax.legend(
1305 handles,
1306 labels,
1307 bbox_to_anchor=(0.5, 0.12),
1308 bbox_transform=fig.transFigure,
1309 ncol=3,
1310 loc='upper center',
1311 fancybox=True)
1313 plt.subplots_adjust(bottom=0.22,
1314 wspace=0.05)
1316 mpl_show(plt)
1319def command_upgrade(args):
1320 parser, options, args = cl_parse('upgrade', args)
1321 store_dirs = get_store_dirs(args)
1322 try:
1323 for store_dir in store_dirs:
1324 store = gf.Store(store_dir)
1325 nup = store.upgrade()
1326 if nup == 0:
1327 print('%s: already up-to-date.' % store_dir)
1328 else:
1329 print('%s: %i file%s upgraded.' % (
1330 store_dir, nup, ['s', ''][nup == 1]))
1332 except gf.StoreError as e:
1333 die(e)
1336def command_addref(args):
1337 parser, options, args = cl_parse('addref', args)
1339 store_dirs = []
1340 filenames = []
1341 for arg in args:
1342 if op.isdir(arg):
1343 store_dirs.append(arg)
1344 elif op.isfile(arg):
1345 filenames.append(arg)
1346 else:
1347 die('invalid argument: %s' % arg)
1349 if not store_dirs:
1350 store_dirs.append('.')
1352 references = []
1353 try:
1354 for filename in args:
1355 references.extend(gf.meta.Reference.from_bibtex(filename=filename))
1356 except ImportError:
1357 die('pybtex module must be installed to use this function')
1359 if not references:
1360 die('no references found')
1362 for store_dir in store_dirs:
1363 try:
1364 store = gf.Store(store_dir)
1365 ids = [ref.id for ref in store.config.references]
1366 for ref in references:
1367 if ref.id in ids:
1368 die('duplicate reference id: %s' % ref.id)
1370 ids.append(ref.id)
1371 store.config.references.append(ref)
1373 store.save_config(make_backup=True)
1375 except gf.StoreError as e:
1376 die(e)
1379def command_qc(args):
1380 parser, options, args = cl_parse('qc', args)
1382 store_dir = get_store_dir(args)
1384 try:
1385 store = gf.Store(store_dir)
1386 s = store.stats()
1387 if s['empty'] != 0:
1388 print('has empty records')
1390 for aname in ['author', 'author_email', 'description', 'references']:
1392 if not getattr(store.config, aname):
1393 print('%s empty' % aname)
1395 except gf.StoreError as e:
1396 die(e)
1399def command_report(args):
1400 from pyrocko.fomosto.report import report_call
1401 report_call.run_program(args)
1404def main(args=None):
1405 '''
1406 CLI entry point for Pyrocko's ``fomosto`` app.
1407 '''
1408 if args is None:
1409 args = sys.argv[1:]
1411 if len(args) < 1:
1412 sys.exit('Usage: %s' % usage)
1414 command = args.pop(0)
1416 if command in subcommands:
1417 globals()['command_' + command](args)
1419 elif command in ('--help', '-h', 'help'):
1420 if command == 'help' and args:
1421 acommand = args[0]
1422 if acommand in subcommands:
1423 globals()['command_' + acommand](['--help'])
1425 sys.exit('Usage: %s' % usage)
1427 else:
1428 sys.exit('fomosto: error: no such subcommand: %s' % command)
1431if __name__ == '__main__':
1432 main()