# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
'''
Fomosto - Green's function database management.
'''
import sys
import re
import os.path as op
import logging
import copy
import shutil
from optparse import OptionParser
from pyrocko import util, trace, gf, cake, io, fomosto
from pyrocko.gui.snuffler import marker
from pyrocko.util import mpl_show
logger = logging.getLogger('pyrocko.apps.fomosto')
km = 1e3
def d2u(d):
return dict((k.replace('-', '_'), v) for (k, v) in d.items())
subcommand_descriptions = {
'init': 'create a new empty GF store',
'build': 'compute GFs and fill into store',
'stats': 'print information about a GF store',
'check': 'check for problems in GF store',
'decimate': 'build decimated variant of a GF store',
'redeploy': 'copy traces from one GF store into another',
'view': 'view selected traces',
'extract': 'extract selected traces',
'import': 'convert Kiwi GFDB to GF store format',
'export': 'convert GF store to Kiwi GFDB format',
'ttt': 'create travel time tables',
'tttview': 'plot travel time table',
'tttextract': 'extract selected travel times',
'tttlsd': 'fix holes in travel time tables',
'sat': 'create stored ray attribute table',
'satview': 'plot stored ray attribute table',
'server': 'run seismosizer server',
'download': 'download GF store from a server',
'modelview': 'plot earthmodels',
'upgrade': 'upgrade store format to latest version',
'addref': 'import citation references to GF store config',
'qc': 'quality check',
'report': "report for Green's Function databases",
'trans2rot': 'create rotational from translational store',
'fd2trans': 'extract elastic10 from elastic10_fd store',
}
subcommand_usages = {
'init': ['init <type> <store-dir> [options]',
'init redeploy <source> <destination> [options]'],
'build': 'build [store-dir] [options]',
'stats': 'stats [store-dir] [options]',
'check': 'check [store-dir] [options]',
'decimate': 'decimate [store-dir] <factor> [options]',
'redeploy': 'redeploy <source> <destination> [options]',
'view': 'view [store-dir] ... [options]',
'extract': 'extract [store-dir] <selection>',
'import': 'import <source> <destination> [options]',
'export': 'export [store-dir] <destination> [options]',
'ttt': 'ttt [store-dir] [options]',
'tttview': 'tttview [store-dir] <phase-ids> [options]',
'tttextract': 'tttextract [store-dir] <phase> <selection>',
'tttlsd': 'tttlsd [store-dir] <phase>',
'sat': 'sat [store-dir] [options]',
'satview': 'satview [store-dir] <phase-ids> [options]',
'server': 'server [options] <store-super-dir> ...',
'download': 'download [options] <site> <store-id>',
'modelview': 'modelview <selection>',
'upgrade': 'upgrade [store-dir] ...',
'addref': 'addref [store-dir] ... <filename> ...',
'qc': 'qc [store-dir]',
'report': 'report <subcommand> <arguments>... [options]',
'trans2rot': 'trans2rot <source> <destination>',
'fd2trans': 'fd2trans <source> <destination>',
}
subcommands = subcommand_descriptions.keys()
program_name = 'fomosto'
usage = program_name + ''' <subcommand> <arguments> ... [options]
Subcommands:
init %(init)s
build %(build)s
stats %(stats)s
check %(check)s
decimate %(decimate)s
redeploy %(redeploy)s
view %(view)s
extract %(extract)s
import %(import)s
export %(export)s
ttt %(ttt)s
tttview %(tttview)s
tttextract %(tttextract)s
tttlsd %(tttlsd)s
sat %(sat)s
satview %(satview)s
server %(server)s
download %(download)s
modelview %(modelview)s
upgrade %(upgrade)s
addref %(addref)s
qc %(qc)s
report %(report)s
trans2rot %(trans2rot)s
fd2trans %(fd2trans)s
To get further help and a list of available options for any subcommand run:
fomosto <subcommand> --help
''' % d2u(subcommand_descriptions)
def add_common_options(parser):
parser.add_option(
'--loglevel',
action='store',
dest='loglevel',
type='choice',
choices=('critical', 'error', 'warning', 'info', 'debug'),
default='info',
help='set logger level to '
'"critical", "error", "warning", "info", or "debug". '
'Default is "%default".')
def process_common_options(options):
util.setup_logging(program_name, options.loglevel)
def cl_parse(command, args, setup=None, details=None):
usage = subcommand_usages[command]
descr = subcommand_descriptions[command]
if isinstance(usage, str):
usage = [usage]
susage = '%s %s' % (program_name, usage[0])
for s in usage[1:]:
susage += '\n%s%s %s' % (' '*7, program_name, s)
description = descr[0].upper() + descr[1:] + '.'
if details:
description = description + ' %s' % details
parser = OptionParser(usage=susage, description=description)
parser.format_description = lambda formatter: description
if setup:
setup(parser)
add_common_options(parser)
(options, args) = parser.parse_args(args)
process_common_options(options)
return parser, options, args
def die(message, err=''):
sys.exit('%s: error: %s \n %s' % (program_name, message, err))
def fomo_wrapper_module(name):
try:
if not re.match(gf.meta.StringID.pattern, name):
raise ValueError('invalid name')
words = name.split('.', 1)
if len(words) == 2:
name, variant = words
else:
name = words[0]
variant = None
name_clean = re.sub(r'[.-]', '_', name)
modname = '.'.join(['pyrocko', 'fomosto', name_clean])
mod = __import__(modname, level=0)
return getattr(mod.fomosto, name_clean), variant
except ValueError:
die('invalid modelling code wrapper name')
except ImportError:
die('''modelling code wrapper "%s" not available or not installed
(module probed: "%s")''' % (name, modname))
def command_init(args):
details = '''
Available modelling backends:
%s
More information at
https://pyrocko.org/docs/current/apps/fomosto/backends.html
''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS])
parser, options, args = cl_parse(
'init', args,
details=details)
if len(args) == 0:
sys.exit(parser.format_help())
if args[0] == 'redeploy':
if len(args) != 3:
parser.error('incorrect number of arguments')
source_dir, dest_dir = args[1:]
try:
source = gf.Store(source_dir)
except gf.StoreError as e:
die(e)
config = copy.deepcopy(source.config)
config.derived_from_id = source.config.id
try:
config_filenames = gf.store.Store.create_editables(
dest_dir, config=config)
except gf.StoreError as e:
die(e)
try:
dest = gf.Store(dest_dir)
except gf.StoreError as e:
die(e)
for k in source.extra_keys():
source_fn = source.get_extra_path(k)
dest_fn = dest.get_extra_path(k)
shutil.copyfile(source_fn, dest_fn)
logger.info(
'(1) configure settings in files:\n %s'
% '\n '.join(config_filenames))
logger.info(
'(2) run "fomosto redeploy <source> <dest>", as needed')
else:
if len(args) != 2:
parser.error('incorrect number of arguments')
(modelling_code_id, store_dir) = args
module, variant = fomo_wrapper_module(modelling_code_id)
try:
config_filenames = module.init(store_dir, variant)
except gf.StoreError as e:
die(e)
logger.info('(1) configure settings in files:\n %s'
% '\n '.join(config_filenames))
logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir)
logger.info('(3) run "fomosto build" in directory "%s"' % store_dir)
def get_store_dir(args):
if len(args) == 1:
store_dir = op.abspath(args.pop(0))
else:
store_dir = op.abspath(op.curdir)
if not op.isdir(store_dir):
die('not a directory: %s' % store_dir)
return store_dir
def get_store_dirs(args):
if len(args) == 0:
store_dirs = [op.abspath(op.curdir)]
else:
store_dirs = [op.abspath(x) for x in args]
for store_dir in store_dirs:
if not op.isdir(store_dir):
die('not a directory: %s' % store_dir)
return store_dirs
def command_build(args):
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser.add_option(
'--nworkers', dest='nworkers', type='int', metavar='N',
help='run N worker processes in parallel')
parser.add_option(
'--continue', dest='continue_', action='store_true',
help='continue suspended build')
parser.add_option(
'--step', dest='step', type='int', metavar='I',
help='process block number IBLOCK')
parser.add_option(
'--block', dest='iblock', type='int', metavar='I',
help='process block number IBLOCK')
parser, options, args = cl_parse('build', args, setup=setup)
store_dir = get_store_dir(args)
try:
if options.step is not None:
step = options.step - 1
else:
step = None
if options.iblock is not None:
iblock = options.iblock - 1
else:
iblock = None
store = gf.Store(store_dir)
module, _ = fomo_wrapper_module(store.config.modelling_code_id)
module.build(
store_dir,
force=options.force,
nworkers=options.nworkers, continue_=options.continue_,
step=step,
iblock=iblock)
except gf.StoreError as e:
die(e)
def command_stats(args):
parser, options, args = cl_parse('stats', args)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
s = store.stats()
except gf.StoreError as e:
die(e)
for k in store.stats_keys:
print('%s: %s' % (k, s[k]))
def command_check(args):
parser, options, args = cl_parse('check', args)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
problems = store.check(show_progress=True)
if problems:
die('problems detected with gf store: %s' % store_dir)
except gf.StoreError as e:
die(e)
def load_config(fn):
try:
config = gf.meta.load(filename=fn)
assert isinstance(config, gf.Config)
except Exception:
die('cannot load gf config from file: %s' % fn)
return config
def command_decimate(args):
def setup(parser):
parser.add_option(
'--config', dest='config_fn', metavar='FILE',
help='use modified spacial sampling given in FILE')
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser, options, args = cl_parse('decimate', args, setup=setup)
try:
decimate = int(args.pop())
except Exception:
parser.error('cannot get <factor> argument')
store_dir = get_store_dir(args)
config = None
if options.config_fn:
config = load_config(options.config_fn)
try:
store = gf.Store(store_dir)
store.make_decimated(decimate, config=config, force=options.force,
show_progress=True)
except gf.StoreError as e:
die(e)
def sindex(args):
return '(%s)' % ', '.join('%g' % x for x in args)
def command_redeploy(args):
parser, options, args = cl_parse('redeploy', args)
if not len(args) == 2:
sys.exit(parser.format_help())
source_store_dir, dest_store_dir = args
try:
source = gf.Store(source_store_dir)
except gf.StoreError as e:
die(e)
try:
gf.store.Store.create_dependants(dest_store_dir)
except gf.StoreError:
pass
try:
dest = gf.Store(dest_store_dir, 'w')
except gf.StoreError as e:
die(e)
show_progress = True
try:
if show_progress:
pbar = util.progressbar('redeploying', dest.config.nrecords)
for i, args in enumerate(dest.config.iter_nodes()):
try:
tr = source.get(args, interpolation='off')
dest.put(args, tr)
except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa
logger.debug('skipping %s, (%s)' % (sindex(args), e))
except gf.store.StoreError as e:
logger.warning('cannot insert %s, (%s)' % (sindex(args), e))
if show_progress:
pbar.update(i+1)
finally:
if show_progress:
pbar.finish()
def command_trans2rot(args):
from pyrocko.fomosto.elastic10_to_rotational8 \
import elastic10_to_rotational8
details = '''
Create a rotational GF store from a translational one using finite differences
approximations (elastic10 or elastic10_fd to rotational8).
If the source store is of type A and receivers at the surface, the free
surface condition is used to replace the vertical derivatives with horizontal
ones. If the input is of type B, the complete set of derivatives are used. It
is important that the spatial sampling of the source store is high enough so
that the finite difference approximations are sufficiently accurate.
If the destination store directory does not exist, it is created and the extent
and spacing of the source store is used for the rotational store. If it does
exist and contains a valid `config`, it is used in place in order to allow for
different input and output extents and spacings. The store must be empty in the
latter case (delete the `index` and `traces` files if needed).
'''
parser, options, args = cl_parse(
'trans2rot', args, details=details)
if not len(args) == 2:
sys.exit(parser.format_help())
source_store_dir, dest_store_dir = args
try:
elastic10_to_rotational8(source_store_dir, dest_store_dir)
except gf.store.StoreError as e:
die(e)
def command_fd2trans(args):
from pyrocko.fomosto.elastic10_fd_to_elastic10 \
import elastic10_fd_to_elastic10
details = '''
Extract regular GF store (elastic10) from a finite difference support GF store
(elastic10_fd).
'''
parser, options, args = cl_parse(
'fd2trans', args, details=details)
if not len(args) == 2:
sys.exit(parser.format_help())
source_store_dir, dest_store_dir = args
try:
elastic10_fd_to_elastic10(source_store_dir, dest_store_dir)
except gf.store.StoreError as e:
die(e)
def command_view(args):
def setup(parser):
parser.add_option(
'--extract',
dest='extract',
metavar='start:stop[:step|@num],...',
help='specify which traces to show')
parser.add_option(
'--show-phases',
dest='showphases',
default=None,
metavar='[phase_id1,phase_id2,...|all]',
help='add phase markers from ttt')
parser.add_option(
'--opengl',
dest='opengl',
action='store_true',
default=None,
help='use OpenGL for drawing')
parser.add_option(
'--no-opengl',
dest='opengl',
action='store_false',
default=None,
help='do not use OpenGL for drawing')
parser, options, args = cl_parse('view', args, setup=setup)
gdef = None
if options.extract:
try:
gdef = gf.meta.parse_grid_spec(options.extract)
except gf.meta.GridSpecError as e:
die(e)
store_dirs = get_store_dirs(args)
alpha = 'abcdefghijklmnopqrstxyz'.upper()
markers = []
traces = []
try:
for istore, store_dir in enumerate(store_dirs):
store = gf.Store(store_dir)
if options.showphases == 'all':
phasenames = [pn.id for pn in store.config.tabulated_phases]
elif options.showphases is not None:
phasenames = options.showphases.split(',')
ii = 0
for args in store.config.iter_extraction(gdef):
gtr = store.get(args)
loc_code = ''
if len(store_dirs) > 1:
loc_code = alpha[istore % len(alpha)]
if gtr:
sta_code = '%04i (%s)' % (
ii, ','.join('%gk' % (x/1000.) for x in args[:-1]))
tmin = gtr.deltat * gtr.itmin
tr = trace.Trace(
'',
sta_code,
loc_code,
'%02i' % args[-1],
ydata=gtr.data,
deltat=gtr.deltat,
tmin=tmin)
if options.showphases:
for phasename in phasenames:
phase_tmin = store.t(phasename, args[:-1])
if phase_tmin:
m = marker.PhaseMarker(
[('',
sta_code,
loc_code,
'%02i' % args[-1])],
phase_tmin,
phase_tmin,
0,
phasename=phasename)
markers.append(m)
traces.append(tr)
ii += 1
except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
die(e)
def setup(win):
viewer = win.get_view()
viewer.menuitem_showboxes.setChecked(False)
viewer.menuitem_colortraces.setChecked(True)
viewer.menuitem_demean.setChecked(False)
viewer.menuitems_sorting[5][0].setChecked(True)
viewer.menuitems_scaling[2][0].setChecked(True)
viewer.sortingmode_change()
viewer.scalingmode_change()
trace.snuffle(
traces, markers=markers, opengl=options.opengl, launch_hook=setup)
def command_extract(args):
def setup(parser):
parser.add_option(
'--format', dest='format', default='mseed',
choices=['mseed', 'sac', 'text', 'yaff'],
help='export to format "mseed", "sac", "text", or "yaff". '
'Default is "mseed".')
fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s'
parser.add_option(
'--output', dest='output_fn', default=fndfl, metavar='TEMPLATE',
help='output path template [default: "%s"]' % fndfl)
parser, options, args = cl_parse('extract', args, setup=setup)
try:
sdef = args.pop()
except Exception:
parser.error('cannot get <selection> argument')
try:
gdef = gf.meta.parse_grid_spec(sdef)
except gf.meta.GridSpecError as e:
die(e)
store_dir = get_store_dir(args)
extensions = {
'mseed': 'mseed',
'sac': 'sac',
'text': 'txt',
'yaff': 'yaff'}
try:
store = gf.Store(store_dir)
for args in store.config.iter_extraction(gdef):
gtr = store.get(args)
if gtr:
tr = trace.Trace(
'', '', '', util.zfmt(store.config.ncomponents) % args[-1],
ydata=gtr.data,
deltat=gtr.deltat,
tmin=gtr.deltat * gtr.itmin)
additional = dict(
args='_'.join('%g' % x for x in args),
irecord=store.str_irecord(args),
extension=extensions[options.format])
io.save(
tr,
options.output_fn,
format=options.format,
additional=additional)
except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
die(e)
def command_import(args):
try:
from tunguska import gfdb
except ImportError:
die('the kiwi tools must be installed to use this feature')
parser, options, args = cl_parse('import', args)
show_progress = True
if not len(args) == 2:
sys.exit(parser.format_help())
source_path, dest_store_dir = args
if op.isdir(source_path):
source_path = op.join(source_path, 'db')
source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path)
db = gfdb.Gfdb(source_path)
config = gf.meta.ConfigTypeA(
id='imported_gfs',
distance_min=db.firstx,
distance_max=db.firstx + (db.nx-1) * db.dx,
distance_delta=db.dx,
source_depth_min=db.firstz,
source_depth_max=db.firstz + (db.nz-1) * db.dz,
source_depth_delta=db.dz,
sample_rate=1.0/db.dt,
ncomponents=db.ng
)
try:
gf.store.Store.create(dest_store_dir, config=config)
dest = gf.Store(dest_store_dir, 'w')
try:
if show_progress:
pbar = util.progressbar(
'importing', dest.config.nrecords/dest.config.ncomponents)
for i, args in enumerate(dest.config.iter_nodes(level=-1)):
source_depth, distance = [float(x) for x in args]
traces = db.get_traces_pyrocko(distance, source_depth)
ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces)
for ig in range(db.ng):
if ig in ig_to_trace:
tr = ig_to_trace[ig]
gf_tr = gf.store.GFTrace(
tr.get_ydata(),
int(round(tr.tmin / tr.deltat)),
tr.deltat)
else:
gf_tr = gf.store.Zero
dest.put((source_depth, distance, ig), gf_tr)
if show_progress:
pbar.update(i+1)
finally:
if show_progress:
pbar.finish()
dest.close()
except gf.StoreError as e:
die(e)
def command_export(args):
from subprocess import Popen, PIPE
try:
from tunguska import gfdb
except ImportError as err:
die('the kiwi tools must be installed to use this feature', err)
def setup(parser):
parser.add_option(
'--nchunks', dest='nchunks', type='int', default=1, metavar='N',
help='split output gfdb into N chunks')
parser, options, args = cl_parse('export', args, setup=setup)
show_progress = True
if len(args) not in (1, 2):
sys.exit(parser.format_help())
target_path = args.pop()
if op.isdir(target_path):
target_path = op.join(target_path, 'kiwi_gfdb')
logger.warning('exported gfdb will be named as "%s.*"' % target_path)
source_store_dir = get_store_dir(args)
source = gf.Store(source_store_dir, 'r')
config = source.config
if not isinstance(config, gf.meta.ConfigTypeA):
die('only stores of type A can be exported to Kiwi format')
if op.isfile(target_path + '.index'):
die('destation already exists')
cmd = [str(x) for x in [
'gfdb_build',
target_path,
options.nchunks,
config.ndistances,
config.nsource_depths,
config.ncomponents,
config.deltat,
config.distance_delta,
config.source_depth_delta,
config.distance_min,
config.source_depth_min]]
p = Popen(cmd, stdin=PIPE)
p.communicate()
out_db = gfdb.Gfdb(target_path)
try:
if show_progress:
pbar = util.progressbar(
'exporting', config.nrecords/config.ncomponents)
for i, (z, x) in enumerate(config.iter_nodes(level=-1)):
data_out = []
for ig in range(config.ncomponents):
try:
tr = source.get((z, x, ig), interpolation='off')
data_out.append((tr.t, tr.data * config.factor))
except gf.store.StoreError as e:
logger.warning(
'cannot get %s, (%s)' % (sindex((z, x, ig)), e))
data_out.append(None)
# put a zero valued sample to no-data zero-traces at a compatible
# time
tmins = [
entry[0][0]
for entry in data_out
if entry is not None and entry[0].size != 0]
if tmins:
tmin = min(tmins)
for entry in data_out:
if entry is not None and entry[0].size == 0:
entry[0].resize(1)
entry[1].resize(1)
entry[0][0] = tmin
entry[1][0] = 0.0
out_db.put_traces_slow(x, z, data_out)
if show_progress:
pbar.update(i+1)
source.close()
finally:
if show_progress:
pbar.finish()
def phasedef_or_horvel(x):
try:
return float(x)
except ValueError:
return cake.PhaseDef(x)
def mkp(s):
return [phasedef_or_horvel(ps) for ps in s.split(',')]
def stored_attribute_table_plots(phase_ids, options, args, attribute):
import numpy as num
from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init
plt = mpl_init()
np = 1
store_dir = get_store_dir(args)
for phase_id in phase_ids:
try:
store = gf.Store(store_dir)
if options.receiver_depth is not None:
receiver_depth = options.receiver_depth * 1000.0
else:
if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)):
receiver_depth = store.config.receiver_depth
elif isinstance(store.config, gf.ConfigTypeB):
receiver_depth = store.config.receiver_depth_min
else:
receiver_depth = 0.0
phase = store.get_stored_phase(phase_id, attribute)
axes = plt.subplot(2, len(phase_ids), np)
labelspace(axes)
xscaled(1. / km, axes)
yscaled(1. / km, axes)
x = None
if isinstance(store.config, gf.ConfigTypeB):
x = (receiver_depth, None, None)
phase.plot_2d(axes, x=x)
axes.set_title(phase_id)
np += 1
except gf.StoreError as e:
die(e)
axes = plt.subplot(2, 1, 2)
num_d = 100
distances = num.linspace(store.config.distance_min,
store.config.distance_max,
num_d)
if options.source_depth is not None:
source_depth = options.source_depth * km
else:
source_depth = store.config.source_depth_min + (
store.config.source_depth_max - store.config.source_depth_min) / 2.
if isinstance(store.config, gf.ConfigTypeA):
attribute_vals = num.empty(num_d)
for phase_id in phase_ids:
attribute_vals[:] = num.NAN
for i, d in enumerate(distances):
if attribute == 'phase':
attribute_vals[i] = store.t(phase_id, (source_depth, d))
ylabel = 'TT [s]'
else:
attribute_vals[i] = store.get_stored_attribute(
phase_id, options.attribute, (source_depth, d))
ylabel = '%s [deg]' % options.attribute
axes.plot(distances / km, attribute_vals, label=phase_id)
axes.set_title('source source_depth %s km' % (source_depth / km))
axes.set_xlabel('distance [km]')
axes.set_ylabel(ylabel)
axes.legend()
plt.tight_layout()
mpl_show(plt)
def command_ttt(args):
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser, options, args = cl_parse('ttt', args, setup=setup)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
store.make_travel_time_tables(force=options.force)
except gf.StoreError as e:
die(e)
def command_tttview(args):
def setup(parser):
parser.add_option(
'--source-depth', dest='source_depth', type=float,
help='Source depth in km')
parser.add_option(
'--receiver-depth', dest='receiver_depth', type=float,
help='Receiver depth in km')
parser, options, args = cl_parse(
'tttview', args, setup=setup,
details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.")
try:
phase_ids = args.pop().split(',')
except Exception:
parser.error('cannot get <phase-ids> argument')
stored_attribute_table_plots(phase_ids, options, args, attribute='phase')
def command_sat(args):
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser.add_option(
'--attribute',
action='store',
dest='attribute',
type='choice',
choices=gf.store.available_stored_tables[1::],
default='takeoff_angle',
help='calculate interpolation table for selected ray attributes.')
parser, options, args = cl_parse('sat', args, setup=setup)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
store.make_stored_table(options.attribute, force=options.force)
except gf.StoreError as e:
die(e)
def command_satview(args):
def setup(parser):
parser.add_option(
'--source-depth', dest='source_depth', type=float,
help='Source depth in km')
parser.add_option(
'--receiver-depth', dest='receiver_depth', type=float,
help='Receiver depth in km')
parser.add_option(
'--attribute',
action='store',
dest='attribute',
type='choice',
choices=gf.store.available_stored_tables[1::],
default='takeoff_angle',
help='view selected ray attribute.')
parser, options, args = cl_parse(
'satview', args, setup=setup,
details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.")
try:
phase_ids = args.pop().split(',')
except Exception:
parser.error('cannot get <phase-ids> argument')
logger.info('Plotting stored attribute %s' % options.attribute)
stored_attribute_table_plots(
phase_ids, options, args, attribute=options.attribute)
def command_tttextract(args):
def setup(parser):
parser.add_option(
'--output', dest='output_fn', metavar='TEMPLATE',
help='output to text files instead of stdout '
'(example TEMPLATE: "extracted/%(args)s.txt")')
parser, options, args = cl_parse('tttextract', args, setup=setup)
try:
sdef = args.pop()
except Exception:
parser.error('cannot get <selection> argument')
try:
sphase = args.pop()
except Exception:
parser.error('cannot get <phase> argument')
try:
phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')]
except gf.meta.InvalidTimingSpecification:
parser.error('invalid phase specification: "%s"' % sphase)
try:
gdef = gf.meta.parse_grid_spec(sdef)
except gf.meta.GridSpecError as e:
die(e)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
for args in store.config.iter_extraction(gdef, level=-1):
s = ['%e' % x for x in args]
for phase in phases:
t = store.t(phase, args)
if t is not None:
s.append('%e' % t)
else:
s.append('nan')
if options.output_fn:
d = dict(
args='_'.join('%e' % x for x in args),
extension='txt')
fn = options.output_fn % d
util.ensuredirs(fn)
with open(fn, 'a') as f:
f.write(' '.join(s))
f.write('\n')
else:
print(' '.join(s))
except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e:
die(e)
def command_tttlsd(args):
def setup(parser):
pass
parser, options, args = cl_parse(
'tttlsd', args,
setup=setup,
details='''
This subcommand fills holes in travel-time-tables by using an eikonal solver to
predict travel-times for the missing values. The approach works best for simple
P or S phases without any reflections or conversions. It creates new
travel-time-tables which are named by adding the suffix `.lsd` to the original
name. E.g. running `fomosto tttlsd begin` would produce a new travel-time-table
`begin.lsd`. The new name can be referenced anywhere where a stored
travel-time-table can be used, e.g. in `extra/qseis`.
''')
try:
sphase_ids = args.pop()
except Exception:
parser.error('cannot get <phase> argument')
try:
phase_ids = [x.strip() for x in sphase_ids.split(',')]
except gf.meta.InvalidTimingSpecification:
parser.error('invalid phase specification: "%s"' % sphase_ids)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
for phase_id in phase_ids:
store.fix_ttt_holes(phase_id)
except gf.StoreError as e:
die(e)
def command_server(args):
from pyrocko.gf import server
def setup(parser):
parser.add_option(
'--port', dest='port', metavar='PORT', type='int', default=8080,
help='serve on port PORT')
parser.add_option(
'--ip', dest='ip', metavar='IP', default='',
help='serve on ip address IP')
parser, options, args = cl_parse('server', args, setup=setup)
engine = gf.LocalEngine(store_superdirs=args)
server.run(options.ip, options.port, engine)
def command_download(args):
from pyrocko.gf import ws
details = '''
Browse pre-calculated Green's function stores online:
https://greens-mill.pyrocko.org
'''
def setup(parser):
parser.add_option(
'--force', dest='force', action='store_true',
help='overwrite existing files')
parser, options, args = cl_parse(
'download', args, setup=setup, details=details)
if len(args) not in (1, 2):
sys.exit(parser.format_help())
if len(args) == 2:
site, store_id = args
if not re.match(gf.meta.StringID.pattern, store_id):
die('invalid store ID')
else:
site, store_id = args[0], None
if site not in gf.ws.g_site_abbr:
if -1 == site.find('://'):
site = 'http://' + site
try:
ws.download_gf_store(site=site, store_id=store_id, force=options.force)
except ws.DownloadError as e:
die(str(e))
def command_modelview(args):
import matplotlib.pyplot as plt
import numpy as num
from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled
mpl_init()
neat_labels = {
'vp': '$v_p$',
'vs': '$v_s$',
'qp': '$Q_p$',
'qs': '$Q_s$',
'rho': '$\\rho$'}
def setup(parser):
parser.add_option(
'--parameters', dest='parameters',
default='vp,vs', metavar='vp,vs,....',
help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs')
units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'}
parser, options, args = cl_parse('modelview', args, setup=setup)
store_dirs = get_store_dirs(args)
parameters = options.parameters.split(',')
fig, axes = plt.subplots(1,
len(parameters),
sharey=True,
figsize=(len(parameters)*3, 5))
if not isinstance(axes, num.ndarray):
axes = [axes]
axes = dict(zip(parameters, axes))
for store_id in store_dirs:
try:
store = gf.Store(store_id)
mod = store.config.earthmodel_1d
z = mod.profile('z')
for p in parameters:
ax = axes[p]
labelspace(ax)
if '/' in p:
p1, p2 = p.split('/')
profile = mod.profile(p1)/mod.profile(p2)
else:
profile = mod.profile(p)
ax.plot(profile, -z, label=store_id.split('/')[-1])
if p in ['vs', 'vp', 'rho']:
xscaled(1./1000., ax)
yscaled(1./km, ax)
except gf.StoreError as e:
die(e)
for p, ax in axes.items():
ax.grid()
if p in neat_labels:
lab = neat_labels[p]
elif all(x in neat_labels for x in p.split('/')):
lab = '/'.join(neat_labels[x] for x in p.split('/'))
else:
lab = p
ax.set_title(lab, y=1.02)
if p in units:
lab += ' ' + units[p]
ax.autoscale()
ax.set_xlabel(lab)
axes[parameters[0]].set_ylabel('Depth [km]')
handles, labels = ax.get_legend_handles_labels()
if len(store_dirs) > 1:
ax.legend(
handles,
labels,
bbox_to_anchor=(0.5, 0.12),
bbox_transform=fig.transFigure,
ncol=3,
loc='upper center',
fancybox=True)
plt.subplots_adjust(bottom=0.22,
wspace=0.05)
mpl_show(plt)
def command_upgrade(args):
parser, options, args = cl_parse('upgrade', args)
store_dirs = get_store_dirs(args)
try:
for store_dir in store_dirs:
store = gf.Store(store_dir)
nup = store.upgrade()
if nup == 0:
print('%s: already up-to-date.' % store_dir)
else:
print('%s: %i file%s upgraded.' % (
store_dir, nup, ['s', ''][nup == 1]))
except gf.StoreError as e:
die(e)
def command_addref(args):
parser, options, args = cl_parse('addref', args)
store_dirs = []
filenames = []
for arg in args:
if op.isdir(arg):
store_dirs.append(arg)
elif op.isfile(arg):
filenames.append(arg)
else:
die('invalid argument: %s' % arg)
if not store_dirs:
store_dirs.append('.')
references = []
try:
for filename in args:
references.extend(gf.meta.Reference.from_bibtex(filename=filename))
except ImportError:
die('pybtex module must be installed to use this function')
if not references:
die('no references found')
for store_dir in store_dirs:
try:
store = gf.Store(store_dir)
ids = [ref.id for ref in store.config.references]
for ref in references:
if ref.id in ids:
die('duplicate reference id: %s' % ref.id)
ids.append(ref.id)
store.config.references.append(ref)
store.save_config(make_backup=True)
except gf.StoreError as e:
die(e)
def command_qc(args):
parser, options, args = cl_parse('qc', args)
store_dir = get_store_dir(args)
try:
store = gf.Store(store_dir)
s = store.stats()
if s['empty'] != 0:
print('has empty records')
for aname in ['author', 'author_email', 'description', 'references']:
if not getattr(store.config, aname):
print('%s empty' % aname)
except gf.StoreError as e:
die(e)
def command_report(args):
from pyrocko.fomosto.report import report_call
report_call.run_program(args)
[docs]def main(args=None):
'''
CLI entry point for Pyrocko's ``fomosto`` app.
'''
if args is None:
args = sys.argv[1:]
if len(args) < 1:
sys.exit('Usage: %s' % usage)
command = args.pop(0)
if command in subcommands:
globals()['command_' + command](args)
elif command in ('--help', '-h', 'help'):
if command == 'help' and args:
acommand = args[0]
if acommand in subcommands:
globals()['command_' + acommand](['--help'])
sys.exit('Usage: %s' % usage)
else:
sys.exit('fomosto: error: no such subcommand: %s' % command)
if __name__ == '__main__':
main()
