1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import sys 

7import re 

8import os.path as op 

9import logging 

10import copy 

11import shutil 

12from optparse import OptionParser 

13 

14from pyrocko import util, trace, gf, cake, io, fomosto 

15from pyrocko.gui.snuffler import marker 

16from pyrocko.util import mpl_show 

17 

18logger = logging.getLogger('pyrocko.apps.fomosto') 

19km = 1e3 

20 

21 

22def d2u(d): 

23 return dict((k.replace('-', '_'), v) for (k, v) in d.items()) 

24 

25 

26subcommand_descriptions = { 

27 'init': 'create a new empty GF store', 

28 'build': 'compute GFs and fill into store', 

29 'stats': 'print information about a GF store', 

30 'check': 'check for problems in GF store', 

31 'decimate': 'build decimated variant of a GF store', 

32 'redeploy': 'copy traces from one GF store into another', 

33 'view': 'view selected traces', 

34 'extract': 'extract selected traces', 

35 'import': 'convert Kiwi GFDB to GF store format', 

36 'export': 'convert GF store to Kiwi GFDB format', 

37 'ttt': 'create travel time tables', 

38 'tttview': 'plot travel time table', 

39 'tttextract': 'extract selected travel times', 

40 'tttlsd': 'fix holes in travel time tables', 

41 'sat': 'create stored ray attribute table', 

42 'satview': 'plot stored ray attribute table', 

43 'server': 'run seismosizer server', 

44 'download': 'download GF store from a server', 

45 'modelview': 'plot earthmodels', 

46 'upgrade': 'upgrade store format to latest version', 

47 'addref': 'import citation references to GF store config', 

48 'qc': 'quality check', 

49 'report': 'report for Green\'s Function databases', 

50} 

51 

52subcommand_usages = { 

53 'init': ['init <type> <store-dir> [options]', 

54 'init redeploy <source> <destination> [options]'], 

55 'build': 'build [store-dir] [options]', 

56 'stats': 'stats [store-dir] [options]', 

57 'check': 'check [store-dir] [options]', 

58 'decimate': 'decimate [store-dir] <factor> [options]', 

59 'redeploy': 'redeploy <source> <destination> [options]', 

60 'view': 'view [store-dir] ... [options]', 

61 'extract': 'extract [store-dir] <selection>', 

62 'import': 'import <source> <destination> [options]', 

63 'export': 'export [store-dir] <destination> [options]', 

64 'ttt': 'ttt [store-dir] [options]', 

65 'tttview': 'tttview [store-dir] <phase-ids> [options]', 

66 'tttextract': 'tttextract [store-dir] <phase> <selection>', 

67 'tttlsd': 'tttlsd [store-dir] <phase>', 

68 'sat': 'sat [store-dir] [options]', 

69 'satview': 'satview [store-dir] <phase-ids> [options]', 

70 'server': 'server [options] <store-super-dir> ...', 

71 'download': 'download [options] <site> <store-id>', 

72 'modelview': 'modelview <selection>', 

73 'upgrade': 'upgrade [store-dir] ...', 

74 'addref': 'addref [store-dir] ... <filename> ...', 

75 'qc': 'qc [store-dir]', 

76 'report': 'report <subcommand> <arguments>... [options]' 

77} 

78 

79subcommands = subcommand_descriptions.keys() 

80 

81program_name = 'fomosto' 

82 

83usage = program_name + ''' <subcommand> <arguments> ... [options] 

84 

85Subcommands: 

86 

87 init %(init)s 

88 build %(build)s 

89 stats %(stats)s 

90 check %(check)s 

91 decimate %(decimate)s 

92 redeploy %(redeploy)s 

93 view %(view)s 

94 extract %(extract)s 

95 import %(import)s 

96 export %(export)s 

97 ttt %(ttt)s 

98 tttview %(tttview)s 

99 tttextract %(tttextract)s 

100 tttlsd %(tttlsd)s 

101 sat %(sat)s 

102 satview %(satview)s 

103 server %(server)s 

104 download %(download)s 

105 modelview %(modelview)s 

106 upgrade %(upgrade)s 

107 addref %(addref)s 

108 qc %(qc)s 

109 report %(report)s 

110 

111To get further help and a list of available options for any subcommand run: 

112 

113 fomosto <subcommand> --help 

114 

115''' % d2u(subcommand_descriptions) 

116 

117 

118def add_common_options(parser): 

119 parser.add_option( 

120 '--loglevel', 

121 action='store', 

122 dest='loglevel', 

123 type='choice', 

124 choices=('critical', 'error', 'warning', 'info', 'debug'), 

125 default='info', 

126 help='set logger level to ' 

127 '"critical", "error", "warning", "info", or "debug". ' 

128 'Default is "%default".') 

129 

130 

131def process_common_options(options): 

132 util.setup_logging(program_name, options.loglevel) 

133 

134 

135def cl_parse(command, args, setup=None, details=None): 

136 usage = subcommand_usages[command] 

137 descr = subcommand_descriptions[command] 

138 

139 if isinstance(usage, str): 

140 usage = [usage] 

141 

142 susage = '%s %s' % (program_name, usage[0]) 

143 for s in usage[1:]: 

144 susage += '\n%s%s %s' % (' '*7, program_name, s) 

145 

146 description = descr[0].upper() + descr[1:] + '.' 

147 

148 if details: 

149 description = description + ' %s' % details 

150 

151 parser = OptionParser(usage=susage, description=description) 

152 parser.format_description = lambda formatter: description 

153 

154 if setup: 

155 setup(parser) 

156 

157 add_common_options(parser) 

158 (options, args) = parser.parse_args(args) 

159 process_common_options(options) 

160 return parser, options, args 

161 

162 

163def die(message, err=''): 

164 sys.exit('%s: error: %s \n %s' % (program_name, message, err)) 

165 

166 

167def fomo_wrapper_module(name): 

168 try: 

169 if not re.match(gf.meta.StringID.pattern, name): 

170 raise ValueError('invalid name') 

171 

172 words = name.split('.', 1) 

173 if len(words) == 2: 

174 name, variant = words 

175 else: 

176 name = words[0] 

177 variant = None 

178 

179 name_clean = re.sub(r'[.-]', '_', name) 

180 modname = '.'.join(['pyrocko', 'fomosto', name_clean]) 

181 mod = __import__(modname, level=0) 

182 return getattr(mod.fomosto, name_clean), variant 

183 

184 except ValueError: 

185 die('invalid modelling code wrapper name') 

186 

187 except ImportError: 

188 die('''modelling code wrapper "%s" not available or not installed 

189 (module probed: "%s")''' % (name, modname)) 

190 

191 

192def command_init(args): 

193 

194 details = ''' 

195 

196Available modelling backends: 

197%s 

198 

199 More information at 

200 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

201''' % '\n'.join([' * %s' % b for b in fomosto.AVAILABLE_BACKENDS]) 

202 

203 parser, options, args = cl_parse( 

204 'init', args, 

205 details=details) 

206 

207 if len(args) == 0: 

208 sys.exit(parser.format_help()) 

209 

210 if args[0] == 'redeploy': 

211 if len(args) != 3: 

212 parser.error('incorrect number of arguments') 

213 

214 source_dir, dest_dir = args[1:] 

215 

216 try: 

217 source = gf.Store(source_dir) 

218 except gf.StoreError as e: 

219 die(e) 

220 

221 config = copy.deepcopy(source.config) 

222 config.derived_from_id = source.config.id 

223 try: 

224 config_filenames = gf.store.Store.create_editables( 

225 dest_dir, config=config) 

226 

227 except gf.StoreError as e: 

228 die(e) 

229 

230 try: 

231 dest = gf.Store(dest_dir) 

232 except gf.StoreError as e: 

233 die(e) 

234 

235 for k in source.extra_keys(): 

236 source_fn = source.get_extra_path(k) 

237 dest_fn = dest.get_extra_path(k) 

238 shutil.copyfile(source_fn, dest_fn) 

239 

240 logger.info( 

241 '(1) configure settings in files:\n %s' 

242 % '\n '.join(config_filenames)) 

243 

244 logger.info( 

245 '(2) run "fomosto redeploy <source> <dest>", as needed') 

246 

247 else: 

248 if len(args) != 2: 

249 parser.error('incorrect number of arguments') 

250 

251 (modelling_code_id, store_dir) = args 

252 

253 module, variant = fomo_wrapper_module(modelling_code_id) 

254 try: 

255 config_filenames = module.init(store_dir, variant) 

256 except gf.StoreError as e: 

257 die(e) 

258 

259 logger.info('(1) configure settings in files:\n %s' 

260 % '\n '.join(config_filenames)) 

261 logger.info('(2) run "fomosto ttt" in directory "%s"' % store_dir) 

262 logger.info('(3) run "fomosto build" in directory "%s"' % store_dir) 

263 

264 

265def get_store_dir(args): 

266 if len(args) == 1: 

267 store_dir = op.abspath(args.pop(0)) 

268 else: 

269 store_dir = op.abspath(op.curdir) 

270 

271 if not op.isdir(store_dir): 

272 die('not a directory: %s' % store_dir) 

273 

274 return store_dir 

275 

276 

277def get_store_dirs(args): 

278 if len(args) == 0: 

279 store_dirs = [op.abspath(op.curdir)] 

280 else: 

281 store_dirs = [op.abspath(x) for x in args] 

282 

283 for store_dir in store_dirs: 

284 if not op.isdir(store_dir): 

285 die('not a directory: %s' % store_dir) 

286 

287 return store_dirs 

288 

289 

290def command_build(args): 

291 

292 def setup(parser): 

293 parser.add_option( 

294 '--force', dest='force', action='store_true', 

295 help='overwrite existing files') 

296 

297 parser.add_option( 

298 '--nworkers', dest='nworkers', type='int', metavar='N', 

299 help='run N worker processes in parallel') 

300 

301 parser.add_option( 

302 '--continue', dest='continue_', action='store_true', 

303 help='continue suspended build') 

304 

305 parser.add_option( 

306 '--step', dest='step', type='int', metavar='I', 

307 help='process block number IBLOCK') 

308 

309 parser.add_option( 

310 '--block', dest='iblock', type='int', metavar='I', 

311 help='process block number IBLOCK') 

312 

313 parser, options, args = cl_parse('build', args, setup=setup) 

314 

315 store_dir = get_store_dir(args) 

316 try: 

317 if options.step is not None: 

318 step = options.step - 1 

319 else: 

320 step = None 

321 

322 if options.iblock is not None: 

323 iblock = options.iblock - 1 

324 else: 

325 iblock = None 

326 

327 store = gf.Store(store_dir) 

328 module, _ = fomo_wrapper_module(store.config.modelling_code_id) 

329 module.build( 

330 store_dir, 

331 force=options.force, 

332 nworkers=options.nworkers, continue_=options.continue_, 

333 step=step, 

334 iblock=iblock) 

335 

336 except gf.StoreError as e: 

337 die(e) 

338 

339 

340def command_stats(args): 

341 

342 parser, options, args = cl_parse('stats', args) 

343 store_dir = get_store_dir(args) 

344 

345 try: 

346 store = gf.Store(store_dir) 

347 s = store.stats() 

348 

349 except gf.StoreError as e: 

350 die(e) 

351 

352 for k in store.stats_keys: 

353 print('%s: %s' % (k, s[k])) 

354 

355 

356def command_check(args): 

357 

358 parser, options, args = cl_parse('check', args) 

359 store_dir = get_store_dir(args) 

360 

361 try: 

362 store = gf.Store(store_dir) 

363 problems = store.check(show_progress=True) 

364 if problems: 

365 die('problems detected with gf store: %s' % store_dir) 

366 

367 except gf.StoreError as e: 

368 die(e) 

369 

370 

371def load_config(fn): 

372 try: 

373 config = gf.meta.load(filename=fn) 

374 assert isinstance(config, gf.Config) 

375 

376 except Exception: 

377 die('cannot load gf config from file: %s' % fn) 

378 

379 return config 

380 

381 

382def command_decimate(args): 

383 

384 def setup(parser): 

385 parser.add_option( 

386 '--config', dest='config_fn', metavar='FILE', 

387 help='use modified spacial sampling given in FILE') 

388 

389 parser.add_option( 

390 '--force', dest='force', action='store_true', 

391 help='overwrite existing files') 

392 

393 parser, options, args = cl_parse('decimate', args, setup=setup) 

394 try: 

395 decimate = int(args.pop()) 

396 except Exception: 

397 parser.error('cannot get <factor> argument') 

398 

399 store_dir = get_store_dir(args) 

400 

401 config = None 

402 if options.config_fn: 

403 config = load_config(options.config_fn) 

404 

405 try: 

406 store = gf.Store(store_dir) 

407 store.make_decimated(decimate, config=config, force=options.force, 

408 show_progress=True) 

409 

410 except gf.StoreError as e: 

411 die(e) 

412 

413 

414def sindex(args): 

415 return '(%s)' % ', '.join('%g' % x for x in args) 

416 

417 

418def command_redeploy(args): 

419 

420 parser, options, args = cl_parse('redeploy', args) 

421 

422 if not len(args) == 2: 

423 sys.exit(parser.format_help()) 

424 

425 source_store_dir, dest_store_dir = args 

426 

427 try: 

428 source = gf.Store(source_store_dir) 

429 except gf.StoreError as e: 

430 die(e) 

431 

432 try: 

433 gf.store.Store.create_dependants(dest_store_dir) 

434 except gf.StoreError: 

435 pass 

436 

437 try: 

438 dest = gf.Store(dest_store_dir, 'w') 

439 

440 except gf.StoreError as e: 

441 die(e) 

442 

443 show_progress = True 

444 

445 try: 

446 if show_progress: 

447 pbar = util.progressbar('redeploying', dest.config.nrecords) 

448 

449 for i, args in enumerate(dest.config.iter_nodes()): 

450 try: 

451 tr = source.get(args, interpolation='off') 

452 dest.put(args, tr) 

453 

454 except (gf.meta.OutOfBounds, gf.store.NotAllowedToInterpolate) as e: # noqa 

455 logger.debug('skipping %s, (%s)' % (sindex(args), e)) 

456 

457 except gf.store.StoreError as e: 

458 logger.warning('cannot insert %s, (%s)' % (sindex(args), e)) 

459 

460 if show_progress: 

461 pbar.update(i+1) 

462 

463 finally: 

464 if show_progress: 

465 pbar.finish() 

466 

467 

468def command_view(args): 

469 def setup(parser): 

470 parser.add_option( 

471 '--extract', 

472 dest='extract', 

473 metavar='start:stop[:step|@num],...', 

474 help='specify which traces to show') 

475 

476 parser.add_option( 

477 '--show-phases', 

478 dest='showphases', 

479 default=None, 

480 metavar='[phase_id1,phase_id2,...|all]', 

481 help='add phase markers from ttt') 

482 

483 parser.add_option( 

484 '--opengl', 

485 dest='opengl', 

486 action='store_true', 

487 default=None, 

488 help='use OpenGL for drawing') 

489 

490 parser.add_option( 

491 '--no-opengl', 

492 dest='opengl', 

493 action='store_false', 

494 default=None, 

495 help='do not use OpenGL for drawing') 

496 

497 parser, options, args = cl_parse('view', args, setup=setup) 

498 

499 gdef = None 

500 if options.extract: 

501 try: 

502 gdef = gf.meta.parse_grid_spec(options.extract) 

503 except gf.meta.GridSpecError as e: 

504 die(e) 

505 

506 store_dirs = get_store_dirs(args) 

507 

508 alpha = 'abcdefghijklmnopqrstxyz'.upper() 

509 

510 markers = [] 

511 traces = [] 

512 

513 try: 

514 for istore, store_dir in enumerate(store_dirs): 

515 store = gf.Store(store_dir) 

516 if options.showphases == 'all': 

517 phasenames = [pn.id for pn in store.config.tabulated_phases] 

518 elif options.showphases is not None: 

519 phasenames = options.showphases.split(',') 

520 

521 ii = 0 

522 for args in store.config.iter_extraction(gdef): 

523 gtr = store.get(args) 

524 

525 loc_code = '' 

526 if len(store_dirs) > 1: 

527 loc_code = alpha[istore % len(alpha)] 

528 

529 if gtr: 

530 

531 sta_code = '%04i (%s)' % ( 

532 ii, ','.join('%gk' % (x/1000.) for x in args[:-1])) 

533 tmin = gtr.deltat * gtr.itmin 

534 

535 tr = trace.Trace( 

536 '', 

537 sta_code, 

538 loc_code, 

539 '%02i' % args[-1], 

540 ydata=gtr.data, 

541 deltat=gtr.deltat, 

542 tmin=tmin) 

543 

544 if options.showphases: 

545 for phasename in phasenames: 

546 phase_tmin = store.t(phasename, args[:-1]) 

547 if phase_tmin: 

548 m = marker.PhaseMarker( 

549 [('', 

550 sta_code, 

551 loc_code, 

552 '%02i' % args[-1])], 

553 phase_tmin, 

554 phase_tmin, 

555 0, 

556 phasename=phasename) 

557 markers.append(m) 

558 

559 traces.append(tr) 

560 

561 ii += 1 

562 

563 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e: 

564 die(e) 

565 

566 def setup(win): 

567 viewer = win.get_view() 

568 viewer.menuitem_showboxes.setChecked(False) 

569 viewer.menuitem_colortraces.setChecked(True) 

570 viewer.menuitem_demean.setChecked(False) 

571 viewer.menuitems_sorting[5][0].setChecked(True) 

572 viewer.menuitems_scaling[2][0].setChecked(True) 

573 viewer.sortingmode_change() 

574 viewer.scalingmode_change() 

575 

576 trace.snuffle( 

577 traces, markers=markers, opengl=options.opengl, launch_hook=setup) 

578 

579 

580def command_extract(args): 

581 def setup(parser): 

582 parser.add_option( 

583 '--format', dest='format', default='mseed', 

584 choices=['mseed', 'sac', 'text', 'yaff'], 

585 help='export to format "mseed", "sac", "text", or "yaff". ' 

586 'Default is "mseed".') 

587 

588 fndfl = 'extracted/%(irecord)s_%(args)s.%(extension)s' 

589 parser.add_option( 

590 '--output', dest='output_fn', default=fndfl, metavar='TEMPLATE', 

591 help='output path template [default: "%s"]' % fndfl) 

592 

593 parser, options, args = cl_parse('extract', args, setup=setup) 

594 try: 

595 sdef = args.pop() 

596 except Exception: 

597 parser.error('cannot get <selection> argument') 

598 

599 try: 

600 gdef = gf.meta.parse_grid_spec(sdef) 

601 except gf.meta.GridSpecError as e: 

602 die(e) 

603 

604 store_dir = get_store_dir(args) 

605 

606 extensions = { 

607 'mseed': 'mseed', 

608 'sac': 'sac', 

609 'text': 'txt', 

610 'yaff': 'yaff'} 

611 

612 try: 

613 store = gf.Store(store_dir) 

614 for args in store.config.iter_extraction(gdef): 

615 gtr = store.get(args) 

616 if gtr: 

617 tr = trace.Trace( 

618 '', '', '', util.zfmt(store.config.ncomponents) % args[-1], 

619 ydata=gtr.data, 

620 deltat=gtr.deltat, 

621 tmin=gtr.deltat * gtr.itmin) 

622 

623 additional = dict( 

624 args='_'.join('%g' % x for x in args), 

625 irecord=store.str_irecord(args), 

626 extension=extensions[options.format]) 

627 

628 io.save( 

629 tr, 

630 options.output_fn, 

631 format=options.format, 

632 additional=additional) 

633 

634 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e: 

635 die(e) 

636 

637 

638def command_import(args): 

639 try: 

640 from tunguska import gfdb 

641 except ImportError: 

642 die('the kiwi tools must be installed to use this feature') 

643 

644 parser, options, args = cl_parse('import', args) 

645 

646 show_progress = True 

647 

648 if not len(args) == 2: 

649 sys.exit(parser.format_help()) 

650 

651 source_path, dest_store_dir = args 

652 

653 if op.isdir(source_path): 

654 source_path = op.join(source_path, 'db') 

655 

656 source_path = re.sub(r'(\.\d+\.chunk|\.index)$', '', source_path) 

657 

658 db = gfdb.Gfdb(source_path) 

659 

660 config = gf.meta.ConfigTypeA( 

661 id='imported_gfs', 

662 distance_min=db.firstx, 

663 distance_max=db.firstx + (db.nx-1) * db.dx, 

664 distance_delta=db.dx, 

665 source_depth_min=db.firstz, 

666 source_depth_max=db.firstz + (db.nz-1) * db.dz, 

667 source_depth_delta=db.dz, 

668 sample_rate=1.0/db.dt, 

669 ncomponents=db.ng 

670 ) 

671 

672 try: 

673 gf.store.Store.create(dest_store_dir, config=config) 

674 dest = gf.Store(dest_store_dir, 'w') 

675 try: 

676 if show_progress: 

677 pbar = util.progressbar( 

678 'importing', dest.config.nrecords/dest.config.ncomponents) 

679 

680 for i, args in enumerate(dest.config.iter_nodes(level=-1)): 

681 source_depth, distance = [float(x) for x in args] 

682 traces = db.get_traces_pyrocko(distance, source_depth) 

683 ig_to_trace = dict((tr.meta['ig']-1, tr) for tr in traces) 

684 for ig in range(db.ng): 

685 if ig in ig_to_trace: 

686 tr = ig_to_trace[ig] 

687 gf_tr = gf.store.GFTrace( 

688 tr.get_ydata(), 

689 int(round(tr.tmin / tr.deltat)), 

690 tr.deltat) 

691 

692 else: 

693 gf_tr = gf.store.Zero 

694 

695 dest.put((source_depth, distance, ig), gf_tr) 

696 

697 if show_progress: 

698 pbar.update(i+1) 

699 

700 finally: 

701 if show_progress: 

702 pbar.finish() 

703 

704 dest.close() 

705 

706 except gf.StoreError as e: 

707 die(e) 

708 

709 

710def command_export(args): 

711 from subprocess import Popen, PIPE 

712 

713 try: 

714 from tunguska import gfdb 

715 except ImportError as err: 

716 die('the kiwi tools must be installed to use this feature', err) 

717 

718 def setup(parser): 

719 parser.add_option( 

720 '--nchunks', dest='nchunks', type='int', default=1, metavar='N', 

721 help='split output gfdb into N chunks') 

722 

723 parser, options, args = cl_parse('export', args, setup=setup) 

724 

725 show_progress = True 

726 

727 if len(args) not in (1, 2): 

728 sys.exit(parser.format_help()) 

729 

730 target_path = args.pop() 

731 if op.isdir(target_path): 

732 target_path = op.join(target_path, 'kiwi_gfdb') 

733 logger.warning('exported gfdb will be named as "%s.*"' % target_path) 

734 

735 source_store_dir = get_store_dir(args) 

736 

737 source = gf.Store(source_store_dir, 'r') 

738 config = source.config 

739 

740 if not isinstance(config, gf.meta.ConfigTypeA): 

741 die('only stores of type A can be exported to Kiwi format') 

742 

743 if op.isfile(target_path + '.index'): 

744 die('destation already exists') 

745 

746 cmd = [str(x) for x in [ 

747 'gfdb_build', 

748 target_path, 

749 options.nchunks, 

750 config.ndistances, 

751 config.nsource_depths, 

752 config.ncomponents, 

753 config.deltat, 

754 config.distance_delta, 

755 config.source_depth_delta, 

756 config.distance_min, 

757 config.source_depth_min]] 

758 

759 p = Popen(cmd, stdin=PIPE) 

760 p.communicate() 

761 

762 out_db = gfdb.Gfdb(target_path) 

763 

764 try: 

765 if show_progress: 

766 pbar = util.progressbar( 

767 'exporting', config.nrecords/config.ncomponents) 

768 

769 for i, (z, x) in enumerate(config.iter_nodes(level=-1)): 

770 

771 data_out = [] 

772 for ig in range(config.ncomponents): 

773 try: 

774 tr = source.get((z, x, ig), interpolation='off') 

775 data_out.append((tr.t, tr.data * config.factor)) 

776 

777 except gf.store.StoreError as e: 

778 logger.warning( 

779 'cannot get %s, (%s)' % (sindex((z, x, ig)), e)) 

780 data_out.append(None) 

781 

782 # put a zero valued sample to no-data zero-traces at a compatible 

783 # time 

784 tmins = [ 

785 entry[0][0] 

786 for entry in data_out 

787 if entry is not None and entry[0].size != 0] 

788 

789 if tmins: 

790 tmin = min(tmins) 

791 for entry in data_out: 

792 if entry is not None and entry[0].size == 0: 

793 entry[0].resize(1) 

794 entry[1].resize(1) 

795 entry[0][0] = tmin 

796 entry[1][0] = 0.0 

797 

798 out_db.put_traces_slow(x, z, data_out) 

799 

800 if show_progress: 

801 pbar.update(i+1) 

802 

803 source.close() 

804 

805 finally: 

806 if show_progress: 

807 pbar.finish() 

808 

809 

810def phasedef_or_horvel(x): 

811 try: 

812 return float(x) 

813 except ValueError: 

814 return cake.PhaseDef(x) 

815 

816 

817def mkp(s): 

818 return [phasedef_or_horvel(ps) for ps in s.split(',')] 

819 

820 

821def stored_attribute_table_plots(phase_ids, options, args, attribute): 

822 import numpy as num 

823 from pyrocko.plot.cake_plot import labelspace, xscaled, yscaled, mpl_init 

824 

825 plt = mpl_init() 

826 

827 np = 1 

828 store_dir = get_store_dir(args) 

829 for phase_id in phase_ids: 

830 try: 

831 store = gf.Store(store_dir) 

832 

833 if options.receiver_depth is not None: 

834 receiver_depth = options.receiver_depth * 1000.0 

835 else: 

836 if isinstance(store.config, (gf.ConfigTypeA, gf.ConfigTypeC)): 

837 receiver_depth = store.config.receiver_depth 

838 

839 elif isinstance(store.config, gf.ConfigTypeB): 

840 receiver_depth = store.config.receiver_depth_min 

841 

842 else: 

843 receiver_depth = 0.0 

844 

845 phase = store.get_stored_phase(phase_id, attribute) 

846 axes = plt.subplot(2, len(phase_ids), np) 

847 labelspace(axes) 

848 xscaled(1. / km, axes) 

849 yscaled(1. / km, axes) 

850 x = None 

851 if isinstance(store.config, gf.ConfigTypeB): 

852 x = (receiver_depth, None, None) 

853 

854 phase.plot_2d(axes, x=x) 

855 axes.set_title(phase_id) 

856 np += 1 

857 except gf.StoreError as e: 

858 die(e) 

859 

860 axes = plt.subplot(2, 1, 2) 

861 num_d = 100 

862 distances = num.linspace(store.config.distance_min, 

863 store.config.distance_max, 

864 num_d) 

865 

866 if options.source_depth is not None: 

867 source_depth = options.source_depth * km 

868 else: 

869 source_depth = store.config.source_depth_min + ( 

870 store.config.source_depth_max - store.config.source_depth_min) / 2. 

871 

872 if isinstance(store.config, gf.ConfigTypeA): 

873 attribute_vals = num.empty(num_d) 

874 for phase_id in phase_ids: 

875 attribute_vals[:] = num.NAN 

876 for i, d in enumerate(distances): 

877 if attribute == 'phase': 

878 attribute_vals[i] = store.t(phase_id, (source_depth, d)) 

879 ylabel = 'TT [s]' 

880 else: 

881 attribute_vals[i] = store.get_stored_attribute( 

882 phase_id, options.attribute, (source_depth, d)) 

883 ylabel = '%s [deg]' % options.attribute 

884 

885 axes.plot(distances / km, attribute_vals, label=phase_id) 

886 

887 axes.set_title('source source_depth %s km' % (source_depth / km)) 

888 axes.set_xlabel('distance [km]') 

889 axes.set_ylabel(ylabel) 

890 axes.legend() 

891 

892 plt.tight_layout() 

893 mpl_show(plt) 

894 

895 

896def command_ttt(args): 

897 def setup(parser): 

898 parser.add_option( 

899 '--force', dest='force', action='store_true', 

900 help='overwrite existing files') 

901 

902 parser, options, args = cl_parse('ttt', args, setup=setup) 

903 

904 store_dir = get_store_dir(args) 

905 try: 

906 store = gf.Store(store_dir) 

907 store.make_travel_time_tables(force=options.force) 

908 

909 except gf.StoreError as e: 

910 die(e) 

911 

912 

913def command_tttview(args): 

914 

915 def setup(parser): 

916 parser.add_option( 

917 '--source-depth', dest='source_depth', type=float, 

918 help='Source depth in km') 

919 

920 parser.add_option( 

921 '--receiver-depth', dest='receiver_depth', type=float, 

922 help='Receiver depth in km') 

923 

924 parser, options, args = cl_parse( 

925 'tttview', args, setup=setup, 

926 details="Comma seperated <phase-ids>, eg. 'fomosto tttview Pdiff,S'.") 

927 

928 try: 

929 phase_ids = args.pop().split(',') 

930 except Exception: 

931 parser.error('cannot get <phase-ids> argument') 

932 

933 stored_attribute_table_plots(phase_ids, options, args, attribute='phase') 

934 

935 

936def command_sat(args): 

937 def setup(parser): 

938 parser.add_option( 

939 '--force', dest='force', action='store_true', 

940 help='overwrite existing files') 

941 

942 parser.add_option( 

943 '--attribute', 

944 action='store', 

945 dest='attribute', 

946 type='choice', 

947 choices=gf.store.available_stored_tables[1::], 

948 default='takeoff_angle', 

949 help='calculate interpolation table for selected ray attributes.') 

950 

951 parser, options, args = cl_parse('sat', args, setup=setup) 

952 

953 store_dir = get_store_dir(args) 

954 try: 

955 store = gf.Store(store_dir) 

956 store.make_stored_table(options.attribute, force=options.force) 

957 

958 except gf.StoreError as e: 

959 die(e) 

960 

961 

962def command_satview(args): 

963 

964 def setup(parser): 

965 parser.add_option( 

966 '--source-depth', dest='source_depth', type=float, 

967 help='Source depth in km') 

968 

969 parser.add_option( 

970 '--receiver-depth', dest='receiver_depth', type=float, 

971 help='Receiver depth in km') 

972 

973 parser.add_option( 

974 '--attribute', 

975 action='store', 

976 dest='attribute', 

977 type='choice', 

978 choices=gf.store.available_stored_tables[1::], 

979 default='takeoff_angle', 

980 help='view selected ray attribute.') 

981 

982 parser, options, args = cl_parse( 

983 'satview', args, setup=setup, 

984 details="Comma seperated <phase-ids>, eg. 'fomosto satview Pdiff,S'.") 

985 

986 try: 

987 phase_ids = args.pop().split(',') 

988 except Exception: 

989 parser.error('cannot get <phase-ids> argument') 

990 

991 logger.info('Plotting stored attribute %s' % options.attribute) 

992 

993 stored_attribute_table_plots( 

994 phase_ids, options, args, attribute=options.attribute) 

995 

996 

997def command_tttextract(args): 

998 def setup(parser): 

999 parser.add_option( 

1000 '--output', dest='output_fn', metavar='TEMPLATE', 

1001 help='output to text files instead of stdout ' 

1002 '(example TEMPLATE: "extracted/%(args)s.txt")') 

1003 

1004 parser, options, args = cl_parse('tttextract', args, setup=setup) 

1005 try: 

1006 sdef = args.pop() 

1007 except Exception: 

1008 parser.error('cannot get <selection> argument') 

1009 

1010 try: 

1011 sphase = args.pop() 

1012 except Exception: 

1013 parser.error('cannot get <phase> argument') 

1014 

1015 try: 

1016 phases = [gf.meta.Timing(x.strip()) for x in sphase.split(',')] 

1017 except gf.meta.InvalidTimingSpecification: 

1018 parser.error('invalid phase specification: "%s"' % sphase) 

1019 

1020 try: 

1021 gdef = gf.meta.parse_grid_spec(sdef) 

1022 except gf.meta.GridSpecError as e: 

1023 die(e) 

1024 

1025 store_dir = get_store_dir(args) 

1026 

1027 try: 

1028 store = gf.Store(store_dir) 

1029 for args in store.config.iter_extraction(gdef, level=-1): 

1030 s = ['%e' % x for x in args] 

1031 for phase in phases: 

1032 t = store.t(phase, args) 

1033 if t is not None: 

1034 s.append('%e' % t) 

1035 else: 

1036 s.append('nan') 

1037 

1038 if options.output_fn: 

1039 d = dict( 

1040 args='_'.join('%e' % x for x in args), 

1041 extension='txt') 

1042 

1043 fn = options.output_fn % d 

1044 util.ensuredirs(fn) 

1045 with open(fn, 'a') as f: 

1046 f.write(' '.join(s)) 

1047 f.write('\n') 

1048 else: 

1049 print(' '.join(s)) 

1050 

1051 except (gf.meta.GridSpecError, gf.StoreError, gf.meta.OutOfBounds) as e: 

1052 die(e) 

1053 

1054 

1055def command_tttlsd(args): 

1056 

1057 def setup(parser): 

1058 pass 

1059 

1060 parser, options, args = cl_parse('tttlsd', args, setup=setup) 

1061 

1062 try: 

1063 sphase_ids = args.pop() 

1064 except Exception: 

1065 parser.error('cannot get <phase> argument') 

1066 

1067 try: 

1068 phase_ids = [x.strip() for x in sphase_ids.split(',')] 

1069 except gf.meta.InvalidTimingSpecification: 

1070 parser.error('invalid phase specification: "%s"' % sphase_ids) 

1071 

1072 store_dir = get_store_dir(args) 

1073 

1074 try: 

1075 store = gf.Store(store_dir) 

1076 for phase_id in phase_ids: 

1077 store.fix_ttt_holes(phase_id) 

1078 

1079 except gf.StoreError as e: 

1080 die(e) 

1081 

1082 

1083def command_server(args): 

1084 from pyrocko.gf import server 

1085 

1086 def setup(parser): 

1087 parser.add_option( 

1088 '--port', dest='port', metavar='PORT', type='int', default=8080, 

1089 help='serve on port PORT') 

1090 

1091 parser.add_option( 

1092 '--ip', dest='ip', metavar='IP', default='', 

1093 help='serve on ip address IP') 

1094 

1095 parser, options, args = cl_parse('server', args, setup=setup) 

1096 

1097 engine = gf.LocalEngine(store_superdirs=args) 

1098 server.run(options.ip, options.port, engine) 

1099 

1100 

1101def command_download(args): 

1102 from pyrocko.gf import ws 

1103 

1104 details = ''' 

1105 

1106Browse pre-calculated Green's function stores online: 

1107 

1108 https://greens-mill.pyrocko.org 

1109''' 

1110 

1111 def setup(parser): 

1112 parser.add_option( 

1113 '--force', dest='force', action='store_true', 

1114 help='overwrite existing files') 

1115 

1116 parser, options, args = cl_parse( 

1117 'download', args, setup=setup, details=details) 

1118 if not len(args) in (1, 2): 

1119 sys.exit(parser.format_help()) 

1120 

1121 if len(args) == 2: 

1122 site, store_id = args 

1123 if not re.match(gf.meta.StringID.pattern, store_id): 

1124 die('invalid store ID') 

1125 else: 

1126 site, store_id = args[0], None 

1127 

1128 if site not in gf.ws.g_site_abbr: 

1129 if -1 == site.find('://'): 

1130 site = 'http://' + site 

1131 

1132 try: 

1133 ws.download_gf_store(site=site, store_id=store_id, force=options.force) 

1134 except ws.DownloadError as e: 

1135 die(str(e)) 

1136 

1137 

1138def command_modelview(args): 

1139 

1140 import matplotlib.pyplot as plt 

1141 import numpy as num 

1142 from pyrocko.plot.cake_plot import mpl_init, labelspace, xscaled, yscaled 

1143 mpl_init() 

1144 

1145 neat_labels = { 

1146 'vp': '$v_p$', 

1147 'vs': '$v_s$', 

1148 'qp': '$Q_p$', 

1149 'qs': '$Q_s$', 

1150 'rho': '$\\rho$'} 

1151 

1152 def setup(parser): 

1153 parser.add_option( 

1154 '--parameters', dest='parameters', 

1155 default='vp,vs', metavar='vp,vs,....', 

1156 help='select one or several of vp, vs, rho, qp, qs, vp/vs, qp/qs') 

1157 

1158 units = {'vp': '[km/s]', 'vs': '[km/s]', 'rho': '[g/cm^3]'} 

1159 

1160 parser, options, args = cl_parse('modelview', args, setup=setup) 

1161 

1162 store_dirs = get_store_dirs(args) 

1163 

1164 parameters = options.parameters.split(',') 

1165 

1166 fig, axes = plt.subplots(1, 

1167 len(parameters), 

1168 sharey=True, 

1169 figsize=(len(parameters)*3, 5)) 

1170 

1171 if not isinstance(axes, num.ndarray): 

1172 axes = [axes] 

1173 

1174 axes = dict(zip(parameters, axes)) 

1175 

1176 for store_id in store_dirs: 

1177 try: 

1178 store = gf.Store(store_id) 

1179 mod = store.config.earthmodel_1d 

1180 z = mod.profile('z') 

1181 

1182 for p in parameters: 

1183 ax = axes[p] 

1184 labelspace(ax) 

1185 if '/' in p: 

1186 p1, p2 = p.split('/') 

1187 profile = mod.profile(p1)/mod.profile(p2) 

1188 else: 

1189 profile = mod.profile(p) 

1190 

1191 ax.plot(profile, -z, label=store_id.split('/')[-1]) 

1192 if p in ['vs', 'vp', 'rho']: 

1193 xscaled(1./1000., ax) 

1194 

1195 yscaled(1./km, ax) 

1196 

1197 except gf.StoreError as e: 

1198 die(e) 

1199 

1200 for p, ax in axes.items(): 

1201 ax.grid() 

1202 if p in neat_labels: 

1203 lab = neat_labels[p] 

1204 elif all(x in neat_labels for x in p.split('/')): 

1205 lab = '/'.join(neat_labels[x] for x in p.split('/')) 

1206 else: 

1207 lab = p 

1208 

1209 ax.set_title(lab, y=1.02) 

1210 

1211 if p in units: 

1212 lab += ' ' + units[p] 

1213 

1214 ax.autoscale() 

1215 ax.set_xlabel(lab) 

1216 

1217 axes[parameters[0]].set_ylabel('Depth [km]') 

1218 

1219 handles, labels = ax.get_legend_handles_labels() 

1220 

1221 if len(store_dirs) > 1: 

1222 ax.legend( 

1223 handles, 

1224 labels, 

1225 bbox_to_anchor=(0.5, 0.12), 

1226 bbox_transform=fig.transFigure, 

1227 ncol=3, 

1228 loc='upper center', 

1229 fancybox=True) 

1230 

1231 plt.subplots_adjust(bottom=0.22, 

1232 wspace=0.05) 

1233 

1234 mpl_show(plt) 

1235 

1236 

1237def command_upgrade(args): 

1238 parser, options, args = cl_parse('upgrade', args) 

1239 store_dirs = get_store_dirs(args) 

1240 try: 

1241 for store_dir in store_dirs: 

1242 store = gf.Store(store_dir) 

1243 nup = store.upgrade() 

1244 if nup == 0: 

1245 print('%s: already up-to-date.' % store_dir) 

1246 else: 

1247 print('%s: %i file%s upgraded.' % ( 

1248 store_dir, nup, ['s', ''][nup == 1])) 

1249 

1250 except gf.StoreError as e: 

1251 die(e) 

1252 

1253 

1254def command_addref(args): 

1255 parser, options, args = cl_parse('addref', args) 

1256 

1257 store_dirs = [] 

1258 filenames = [] 

1259 for arg in args: 

1260 if op.isdir(arg): 

1261 store_dirs.append(arg) 

1262 elif op.isfile(arg): 

1263 filenames.append(arg) 

1264 else: 

1265 die('invalid argument: %s' % arg) 

1266 

1267 if not store_dirs: 

1268 store_dirs.append('.') 

1269 

1270 references = [] 

1271 try: 

1272 for filename in args: 

1273 references.extend(gf.meta.Reference.from_bibtex(filename=filename)) 

1274 except ImportError: 

1275 die('pybtex module must be installed to use this function') 

1276 

1277 if not references: 

1278 die('no references found') 

1279 

1280 for store_dir in store_dirs: 

1281 try: 

1282 store = gf.Store(store_dir) 

1283 ids = [ref.id for ref in store.config.references] 

1284 for ref in references: 

1285 if ref.id in ids: 

1286 die('duplicate reference id: %s' % ref.id) 

1287 

1288 ids.append(ref.id) 

1289 store.config.references.append(ref) 

1290 

1291 store.save_config(make_backup=True) 

1292 

1293 except gf.StoreError as e: 

1294 die(e) 

1295 

1296 

1297def command_qc(args): 

1298 parser, options, args = cl_parse('qc', args) 

1299 

1300 store_dir = get_store_dir(args) 

1301 

1302 try: 

1303 store = gf.Store(store_dir) 

1304 s = store.stats() 

1305 if s['empty'] != 0: 

1306 print('has empty records') 

1307 

1308 for aname in ['author', 'author_email', 'description', 'references']: 

1309 

1310 if not getattr(store.config, aname): 

1311 print('%s empty' % aname) 

1312 

1313 except gf.StoreError as e: 

1314 die(e) 

1315 

1316 

1317def command_report(args): 

1318 from pyrocko.fomosto.report import report_call 

1319 report_call.run_program(args) 

1320 

1321 

1322def main(args=None): 

1323 if args is None: 

1324 args = sys.argv[1:] 

1325 

1326 if len(args) < 1: 

1327 sys.exit('Usage: %s' % usage) 

1328 

1329 command = args.pop(0) 

1330 

1331 if command in subcommands: 

1332 globals()['command_' + command](args) 

1333 

1334 elif command in ('--help', '-h', 'help'): 

1335 if command == 'help' and args: 

1336 acommand = args[0] 

1337 if acommand in subcommands: 

1338 globals()['command_' + acommand](['--help']) 

1339 

1340 sys.exit('Usage: %s' % usage) 

1341 

1342 else: 

1343 sys.exit('fomosto: error: no such subcommand: %s' % command) 

1344 

1345 

1346if __name__ == '__main__': 

1347 main()