1#!/usr/bin/env python 

2# http://pyrocko.org - GPLv3 

3# 

4# The Pyrocko Developers, 21st Century 

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

6 

7import sys 

8import re 

9import os.path as op 

10import logging 

11import copy 

12import shutil 

13from optparse import OptionParser 

14 

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

16from pyrocko.gui import marker 

17from pyrocko.util import mpl_show 

18 

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

20km = 1e3 

21 

22 

23def d2u(d): 

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

25 

26 

27subcommand_descriptions = { 

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

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

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

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

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

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

34 'view': 'view selected traces', 

35 'extract': 'extract selected traces', 

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

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

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

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

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

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

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

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

44 'server': 'run seismosizer server', 

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

46 'modelview': 'plot earthmodels', 

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

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

49 'qc': 'quality check', 

50 'report': "report for Green's Function databases", 

51} 

52 

53subcommand_usages = { 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

73 'modelview': 'modelview <selection>', 

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

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

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

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

78} 

79 

80subcommands = subcommand_descriptions.keys() 

81 

82program_name = 'fomosto' 

83 

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

85 

86Subcommands: 

87 

88 init %(init)s 

89 build %(build)s 

90 stats %(stats)s 

91 check %(check)s 

92 decimate %(decimate)s 

93 redeploy %(redeploy)s 

94 view %(view)s 

95 extract %(extract)s 

96 import %(import)s 

97 export %(export)s 

98 ttt %(ttt)s 

99 tttview %(tttview)s 

100 tttextract %(tttextract)s 

101 tttlsd %(tttlsd)s 

102 sat %(sat)s 

103 satview %(satview)s 

104 server %(server)s 

105 download %(download)s 

106 modelview %(modelview)s 

107 upgrade %(upgrade)s 

108 addref %(addref)s 

109 qc %(qc)s 

110 report %(report)s 

111 

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

113 

114 fomosto <subcommand> --help 

115 

116''' % d2u(subcommand_descriptions) 

117 

118 

119def add_common_options(parser): 

120 parser.add_option( 

121 '--loglevel', 

122 action='store', 

123 dest='loglevel', 

124 type='choice', 

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

126 default='info', 

127 help='set logger level to ' 

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

129 'Default is "%default".') 

130 

131 

132def process_common_options(options): 

133 util.setup_logging(program_name, options.loglevel) 

134 

135 

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

137 usage = subcommand_usages[command] 

138 descr = subcommand_descriptions[command] 

139 

140 if isinstance(usage, str): 

141 usage = [usage] 

142 

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

144 for s in usage[1:]: 

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

146 

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

148 

149 if details: 

150 description = description + ' %s' % details 

151 

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

153 parser.format_description = lambda formatter: description 

154 

155 if setup: 

156 setup(parser) 

157 

158 add_common_options(parser) 

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

160 process_common_options(options) 

161 return parser, options, args 

162 

163 

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

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

166 

167 

168def fomo_wrapper_module(name): 

169 try: 

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

171 raise ValueError('invalid name') 

172 

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

174 if len(words) == 2: 

175 name, variant = words 

176 else: 

177 name = words[0] 

178 variant = None 

179 

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

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

182 mod = __import__(modname, level=0) 

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

184 

185 except ValueError: 

186 die('invalid modelling code wrapper name') 

187 

188 except ImportError: 

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

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

191 

192 

193def command_init(args): 

194 

195 details = ''' 

196 

197Available modelling backends: 

198%s 

199 

200 More information at 

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

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

203 

204 parser, options, args = cl_parse( 

205 'init', args, 

206 details=details) 

207 

208 if len(args) == 0: 

209 sys.exit(parser.format_help()) 

210 

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

212 if len(args) != 3: 

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

214 

215 source_dir, dest_dir = args[1:] 

216 

217 try: 

218 source = gf.Store(source_dir) 

219 except gf.StoreError as e: 

220 die(e) 

221 

222 config = copy.deepcopy(source.config) 

223 config.derived_from_id = source.config.id 

224 try: 

225 config_filenames = gf.store.Store.create_editables( 

226 dest_dir, config=config) 

227 

228 except gf.StoreError as e: 

229 die(e) 

230 

231 try: 

232 dest = gf.Store(dest_dir) 

233 except gf.StoreError as e: 

234 die(e) 

235 

236 for k in source.extra_keys(): 

237 source_fn = source.get_extra_path(k) 

238 dest_fn = dest.get_extra_path(k) 

239 shutil.copyfile(source_fn, dest_fn) 

240 

241 logger.info( 

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

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

244 

245 logger.info( 

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

247 

248 else: 

249 if len(args) != 2: 

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

251 

252 (modelling_code_id, store_dir) = args 

253 

254 module, variant = fomo_wrapper_module(modelling_code_id) 

255 try: 

256 config_filenames = module.init(store_dir, variant) 

257 except gf.StoreError as e: 

258 die(e) 

259 

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

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

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

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

264 

265 

266def get_store_dir(args): 

267 if len(args) == 1: 

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

269 else: 

270 store_dir = op.abspath(op.curdir) 

271 

272 if not op.isdir(store_dir): 

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

274 

275 return store_dir 

276 

277 

278def get_store_dirs(args): 

279 if len(args) == 0: 

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

281 else: 

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

283 

284 for store_dir in store_dirs: 

285 if not op.isdir(store_dir): 

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

287 

288 return store_dirs 

289 

290 

291def command_build(args): 

292 

293 def setup(parser): 

294 parser.add_option( 

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

296 help='overwrite existing files') 

297 

298 parser.add_option( 

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

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

301 

302 parser.add_option( 

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

304 help='continue suspended build') 

305 

306 parser.add_option( 

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

308 help='process block number IBLOCK') 

309 

310 parser.add_option( 

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

312 help='process block number IBLOCK') 

313 

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

315 

316 store_dir = get_store_dir(args) 

317 try: 

318 if options.step is not None: 

319 step = options.step - 1 

320 else: 

321 step = None 

322 

323 if options.iblock is not None: 

324 iblock = options.iblock - 1 

325 else: 

326 iblock = None 

327 

328 store = gf.Store(store_dir) 

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

330 module.build( 

331 store_dir, 

332 force=options.force, 

333 nworkers=options.nworkers, continue_=options.continue_, 

334 step=step, 

335 iblock=iblock) 

336 

337 except gf.StoreError as e: 

338 die(e) 

339 

340 

341def command_stats(args): 

342 

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

344 store_dir = get_store_dir(args) 

345 

346 try: 

347 store = gf.Store(store_dir) 

348 s = store.stats() 

349 

350 except gf.StoreError as e: 

351 die(e) 

352 

353 for k in store.stats_keys: 

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

355 

356 

357def command_check(args): 

358 

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

360 store_dir = get_store_dir(args) 

361 

362 try: 

363 store = gf.Store(store_dir) 

364 problems = store.check(show_progress=True) 

365 if problems: 

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

367 

368 except gf.StoreError as e: 

369 die(e) 

370 

371 

372def load_config(fn): 

373 try: 

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

375 assert isinstance(config, gf.Config) 

376 

377 except Exception: 

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

379 

380 return config 

381 

382 

383def command_decimate(args): 

384 

385 def setup(parser): 

386 parser.add_option( 

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

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

389 

390 parser.add_option( 

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

392 help='overwrite existing files') 

393 

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

395 try: 

396 decimate = int(args.pop()) 

397 except Exception: 

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

399 

400 store_dir = get_store_dir(args) 

401 

402 config = None 

403 if options.config_fn: 

404 config = load_config(options.config_fn) 

405 

406 try: 

407 store = gf.Store(store_dir) 

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

409 show_progress=True) 

410 

411 except gf.StoreError as e: 

412 die(e) 

413 

414 

415def sindex(args): 

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

417 

418 

419def command_redeploy(args): 

420 

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

422 

423 if not len(args) == 2: 

424 sys.exit(parser.format_help()) 

425 

426 source_store_dir, dest_store_dir = args 

427 

428 try: 

429 source = gf.Store(source_store_dir) 

430 except gf.StoreError as e: 

431 die(e) 

432 

433 try: 

434 gf.store.Store.create_dependants(dest_store_dir) 

435 except gf.StoreError: 

436 pass 

437 

438 try: 

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

440 

441 except gf.StoreError as e: 

442 die(e) 

443 

444 show_progress = True 

445 

446 try: 

447 if show_progress: 

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

449 

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

451 try: 

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

453 dest.put(args, tr) 

454 

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

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

457 

458 except gf.store.StoreError as e: 

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

460 

461 if show_progress: 

462 pbar.update(i+1) 

463 

464 finally: 

465 if show_progress: 

466 pbar.finish() 

467 

468 

469def command_view(args): 

470 def setup(parser): 

471 parser.add_option( 

472 '--extract', 

473 dest='extract', 

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

475 help='specify which traces to show') 

476 

477 parser.add_option( 

478 '--show-phases', 

479 dest='showphases', 

480 default=None, 

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

482 help='add phase markers from ttt') 

483 

484 parser.add_option( 

485 '--opengl', 

486 dest='opengl', 

487 action='store_true', 

488 default=None, 

489 help='use OpenGL for drawing') 

490 

491 parser.add_option( 

492 '--no-opengl', 

493 dest='opengl', 

494 action='store_false', 

495 default=None, 

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

497 

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

499 

500 gdef = None 

501 if options.extract: 

502 try: 

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

504 except gf.meta.GridSpecError as e: 

505 die(e) 

506 

507 store_dirs = get_store_dirs(args) 

508 

509 alpha = 'abcdefghijklmnopqrstxyz'.upper() 

510 

511 markers = [] 

512 traces = [] 

513 

514 try: 

515 for istore, store_dir in enumerate(store_dirs): 

516 store = gf.Store(store_dir) 

517 if options.showphases == 'all': 

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

519 elif options.showphases is not None: 

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

521 

522 ii = 0 

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

524 gtr = store.get(args) 

525 

526 loc_code = '' 

527 if len(store_dirs) > 1: 

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

529 

530 if gtr: 

531 

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

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

534 tmin = gtr.deltat * gtr.itmin 

535 

536 tr = trace.Trace( 

537 '', 

538 sta_code, 

539 loc_code, 

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

541 ydata=gtr.data, 

542 deltat=gtr.deltat, 

543 tmin=tmin) 

544 

545 if options.showphases: 

546 for phasename in phasenames: 

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

548 if phase_tmin: 

549 m = marker.PhaseMarker( 

550 [('', 

551 sta_code, 

552 loc_code, 

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

554 phase_tmin, 

555 phase_tmin, 

556 0, 

557 phasename=phasename) 

558 markers.append(m) 

559 

560 traces.append(tr) 

561 

562 ii += 1 

563 

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

565 die(e) 

566 

567 def setup(win): 

568 viewer = win.get_view() 

569 viewer.menuitem_showboxes.setChecked(False) 

570 viewer.menuitem_colortraces.setChecked(True) 

571 viewer.menuitem_demean.setChecked(False) 

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

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

574 viewer.sortingmode_change() 

575 viewer.scalingmode_change() 

576 

577 trace.snuffle( 

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

579 

580 

581def command_extract(args): 

582 def setup(parser): 

583 parser.add_option( 

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

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

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

587 'Default is "mseed".') 

588 

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

590 parser.add_option( 

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

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

593 

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

595 try: 

596 sdef = args.pop() 

597 except Exception: 

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

599 

600 try: 

601 gdef = gf.meta.parse_grid_spec(sdef) 

602 except gf.meta.GridSpecError as e: 

603 die(e) 

604 

605 store_dir = get_store_dir(args) 

606 

607 extensions = { 

608 'mseed': 'mseed', 

609 'sac': 'sac', 

610 'text': 'txt', 

611 'yaff': 'yaff'} 

612 

613 try: 

614 store = gf.Store(store_dir) 

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

616 gtr = store.get(args) 

617 if gtr: 

618 tr = trace.Trace( 

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

620 ydata=gtr.data, 

621 deltat=gtr.deltat, 

622 tmin=gtr.deltat * gtr.itmin) 

623 

624 additional = dict( 

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

626 irecord=store.str_irecord(args), 

627 extension=extensions[options.format]) 

628 

629 io.save( 

630 tr, 

631 options.output_fn, 

632 format=options.format, 

633 additional=additional) 

634 

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

636 die(e) 

637 

638 

639def command_import(args): 

640 try: 

641 from tunguska import gfdb 

642 except ImportError: 

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

644 

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

646 

647 show_progress = True 

648 

649 if not len(args) == 2: 

650 sys.exit(parser.format_help()) 

651 

652 source_path, dest_store_dir = args 

653 

654 if op.isdir(source_path): 

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

656 

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

658 

659 db = gfdb.Gfdb(source_path) 

660 

661 config = gf.meta.ConfigTypeA( 

662 id='imported_gfs', 

663 distance_min=db.firstx, 

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

665 distance_delta=db.dx, 

666 source_depth_min=db.firstz, 

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

668 source_depth_delta=db.dz, 

669 sample_rate=1.0/db.dt, 

670 ncomponents=db.ng 

671 ) 

672 

673 try: 

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

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

676 try: 

677 if show_progress: 

678 pbar = util.progressbar( 

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

680 

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

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

683 traces = db.get_traces_pyrocko(distance, source_depth) 

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

685 for ig in range(db.ng): 

686 if ig in ig_to_trace: 

687 tr = ig_to_trace[ig] 

688 gf_tr = gf.store.GFTrace( 

689 tr.get_ydata(), 

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

691 tr.deltat) 

692 

693 else: 

694 gf_tr = gf.store.Zero 

695 

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

697 

698 if show_progress: 

699 pbar.update(i+1) 

700 

701 finally: 

702 if show_progress: 

703 pbar.finish() 

704 

705 dest.close() 

706 

707 except gf.StoreError as e: 

708 die(e) 

709 

710 

711def command_export(args): 

712 from subprocess import Popen, PIPE 

713 

714 try: 

715 from tunguska import gfdb 

716 except ImportError as err: 

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

718 

719 def setup(parser): 

720 parser.add_option( 

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

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

723 

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

725 

726 show_progress = True 

727 

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

729 sys.exit(parser.format_help()) 

730 

731 target_path = args.pop() 

732 if op.isdir(target_path): 

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

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

735 

736 source_store_dir = get_store_dir(args) 

737 

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

739 config = source.config 

740 

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

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

743 

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

745 die('destation already exists') 

746 

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

748 'gfdb_build', 

749 target_path, 

750 options.nchunks, 

751 config.ndistances, 

752 config.nsource_depths, 

753 config.ncomponents, 

754 config.deltat, 

755 config.distance_delta, 

756 config.source_depth_delta, 

757 config.distance_min, 

758 config.source_depth_min]] 

759 

760 p = Popen(cmd, stdin=PIPE) 

761 p.communicate() 

762 

763 out_db = gfdb.Gfdb(target_path) 

764 

765 try: 

766 if show_progress: 

767 pbar = util.progressbar( 

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

769 

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

771 

772 data_out = [] 

773 for ig in range(config.ncomponents): 

774 try: 

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

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

777 

778 except gf.store.StoreError as e: 

779 logger.warning( 

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

781 data_out.append(None) 

782 

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

784 # time 

785 tmins = [ 

786 entry[0][0] 

787 for entry in data_out 

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

789 

790 if tmins: 

791 tmin = min(tmins) 

792 for entry in data_out: 

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

794 entry[0].resize(1) 

795 entry[1].resize(1) 

796 entry[0][0] = tmin 

797 entry[1][0] = 0.0 

798 

799 out_db.put_traces_slow(x, z, data_out) 

800 

801 if show_progress: 

802 pbar.update(i+1) 

803 

804 source.close() 

805 

806 finally: 

807 if show_progress: 

808 pbar.finish() 

809 

810 

811def phasedef_or_horvel(x): 

812 try: 

813 return float(x) 

814 except ValueError: 

815 return cake.PhaseDef(x) 

816 

817 

818def mkp(s): 

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

820 

821 

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

823 import numpy as num 

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

825 

826 plt = mpl_init() 

827 

828 np = 1 

829 store_dir = get_store_dir(args) 

830 for phase_id in phase_ids: 

831 try: 

832 store = gf.Store(store_dir) 

833 

834 if options.receiver_depth is not None: 

835 receiver_depth = options.receiver_depth * 1000.0 

836 else: 

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

838 receiver_depth = store.config.receiver_depth 

839 

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

841 receiver_depth = store.config.receiver_depth_min 

842 

843 else: 

844 receiver_depth = 0.0 

845 

846 phase = store.get_stored_phase(phase_id, attribute) 

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

848 labelspace(axes) 

849 xscaled(1. / km, axes) 

850 yscaled(1. / km, axes) 

851 x = None 

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

853 x = (receiver_depth, None, None) 

854 

855 phase.plot_2d(axes, x=x) 

856 axes.set_title(phase_id) 

857 np += 1 

858 except gf.StoreError as e: 

859 die(e) 

860 

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

862 num_d = 100 

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

864 store.config.distance_max, 

865 num_d) 

866 

867 if options.source_depth is not None: 

868 source_depth = options.source_depth * km 

869 else: 

870 source_depth = store.config.source_depth_min + ( 

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

872 

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

874 attribute_vals = num.empty(num_d) 

875 for phase_id in phase_ids: 

876 attribute_vals[:] = num.NAN 

877 for i, d in enumerate(distances): 

878 if attribute == 'phase': 

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

880 ylabel = 'TT [s]' 

881 else: 

882 attribute_vals[i] = store.get_stored_attribute( 

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

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

885 

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

887 

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

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

890 axes.set_ylabel(ylabel) 

891 axes.legend() 

892 

893 plt.tight_layout() 

894 mpl_show(plt) 

895 

896 

897def command_ttt(args): 

898 def setup(parser): 

899 parser.add_option( 

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

901 help='overwrite existing files') 

902 

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

904 

905 store_dir = get_store_dir(args) 

906 try: 

907 store = gf.Store(store_dir) 

908 store.make_travel_time_tables(force=options.force) 

909 

910 except gf.StoreError as e: 

911 die(e) 

912 

913 

914def command_tttview(args): 

915 

916 def setup(parser): 

917 parser.add_option( 

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

919 help='Source depth in km') 

920 

921 parser.add_option( 

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

923 help='Receiver depth in km') 

924 

925 parser, options, args = cl_parse( 

926 'tttview', args, setup=setup, 

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

928 

929 try: 

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

931 except Exception: 

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

933 

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

935 

936 

937def command_sat(args): 

938 def setup(parser): 

939 parser.add_option( 

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

941 help='overwrite existing files') 

942 

943 parser.add_option( 

944 '--attribute', 

945 action='store', 

946 dest='attribute', 

947 type='choice', 

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

949 default='takeoff_angle', 

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

951 

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

953 

954 store_dir = get_store_dir(args) 

955 try: 

956 store = gf.Store(store_dir) 

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

958 

959 except gf.StoreError as e: 

960 die(e) 

961 

962 

963def command_satview(args): 

964 

965 def setup(parser): 

966 parser.add_option( 

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

968 help='Source depth in km') 

969 

970 parser.add_option( 

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

972 help='Receiver depth in km') 

973 

974 parser.add_option( 

975 '--attribute', 

976 action='store', 

977 dest='attribute', 

978 type='choice', 

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

980 default='takeoff_angle', 

981 help='view selected ray attribute.') 

982 

983 parser, options, args = cl_parse( 

984 'satview', args, setup=setup, 

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

986 

987 try: 

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

989 except Exception: 

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

991 

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

993 

994 stored_attribute_table_plots( 

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

996 

997 

998def command_tttextract(args): 

999 def setup(parser): 

1000 parser.add_option( 

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

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

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

1004 

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

1006 try: 

1007 sdef = args.pop() 

1008 except Exception: 

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

1010 

1011 try: 

1012 sphase = args.pop() 

1013 except Exception: 

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

1015 

1016 try: 

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

1018 except gf.meta.InvalidTimingSpecification: 

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

1020 

1021 try: 

1022 gdef = gf.meta.parse_grid_spec(sdef) 

1023 except gf.meta.GridSpecError as e: 

1024 die(e) 

1025 

1026 store_dir = get_store_dir(args) 

1027 

1028 try: 

1029 store = gf.Store(store_dir) 

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

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

1032 for phase in phases: 

1033 t = store.t(phase, args) 

1034 if t is not None: 

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

1036 else: 

1037 s.append('nan') 

1038 

1039 if options.output_fn: 

1040 d = dict( 

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

1042 extension='txt') 

1043 

1044 fn = options.output_fn % d 

1045 util.ensuredirs(fn) 

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

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

1048 f.write('\n') 

1049 else: 

1050 print(' '.join(s)) 

1051 

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

1053 die(e) 

1054 

1055 

1056def command_tttlsd(args): 

1057 

1058 def setup(parser): 

1059 pass 

1060 

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

1062 

1063 try: 

1064 sphase_ids = args.pop() 

1065 except Exception: 

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

1067 

1068 try: 

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

1070 except gf.meta.InvalidTimingSpecification: 

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

1072 

1073 store_dir = get_store_dir(args) 

1074 

1075 try: 

1076 store = gf.Store(store_dir) 

1077 for phase_id in phase_ids: 

1078 store.fix_ttt_holes(phase_id) 

1079 

1080 except gf.StoreError as e: 

1081 die(e) 

1082 

1083 

1084def command_server(args): 

1085 from pyrocko.gf import server 

1086 

1087 def setup(parser): 

1088 parser.add_option( 

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

1090 help='serve on port PORT') 

1091 

1092 parser.add_option( 

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

1094 help='serve on ip address IP') 

1095 

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

1097 

1098 engine = gf.LocalEngine(store_superdirs=args) 

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

1100 

1101 

1102def command_download(args): 

1103 from pyrocko.gf import ws 

1104 

1105 details = ''' 

1106 

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

1108 

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

1110''' 

1111 

1112 def setup(parser): 

1113 parser.add_option( 

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

1115 help='overwrite existing files') 

1116 

1117 parser, options, args = cl_parse( 

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

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

1120 sys.exit(parser.format_help()) 

1121 

1122 if len(args) == 2: 

1123 site, store_id = args 

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

1125 die('invalid store ID') 

1126 else: 

1127 site, store_id = args[0], None 

1128 

1129 if site not in gf.ws.g_site_abbr: 

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

1131 site = 'http://' + site 

1132 

1133 try: 

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

1135 except ws.DownloadError as e: 

1136 die(str(e)) 

1137 

1138 

1139def command_modelview(args): 

1140 

1141 import matplotlib.pyplot as plt 

1142 import numpy as num 

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

1144 mpl_init() 

1145 

1146 neat_labels = { 

1147 'vp': '$v_p$', 

1148 'vs': '$v_s$', 

1149 'qp': '$Q_p$', 

1150 'qs': '$Q_s$', 

1151 'rho': '$\\rho$'} 

1152 

1153 def setup(parser): 

1154 parser.add_option( 

1155 '--parameters', dest='parameters', 

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

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

1158 

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

1160 

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

1162 

1163 store_dirs = get_store_dirs(args) 

1164 

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

1166 

1167 fig, axes = plt.subplots(1, 

1168 len(parameters), 

1169 sharey=True, 

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

1171 

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

1173 axes = [axes] 

1174 

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

1176 

1177 for store_id in store_dirs: 

1178 try: 

1179 store = gf.Store(store_id) 

1180 mod = store.config.earthmodel_1d 

1181 z = mod.profile('z') 

1182 

1183 for p in parameters: 

1184 ax = axes[p] 

1185 labelspace(ax) 

1186 if '/' in p: 

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

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

1189 else: 

1190 profile = mod.profile(p) 

1191 

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

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

1194 xscaled(1./1000., ax) 

1195 

1196 yscaled(1./km, ax) 

1197 

1198 except gf.StoreError as e: 

1199 die(e) 

1200 

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

1202 ax.grid() 

1203 if p in neat_labels: 

1204 lab = neat_labels[p] 

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

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

1207 else: 

1208 lab = p 

1209 

1210 ax.set_title(lab, y=1.02) 

1211 

1212 if p in units: 

1213 lab += ' ' + units[p] 

1214 

1215 ax.autoscale() 

1216 ax.set_xlabel(lab) 

1217 

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

1219 

1220 handles, labels = ax.get_legend_handles_labels() 

1221 

1222 if len(store_dirs) > 1: 

1223 ax.legend( 

1224 handles, 

1225 labels, 

1226 bbox_to_anchor=(0.5, 0.12), 

1227 bbox_transform=fig.transFigure, 

1228 ncol=3, 

1229 loc='upper center', 

1230 fancybox=True) 

1231 

1232 plt.subplots_adjust(bottom=0.22, 

1233 wspace=0.05) 

1234 

1235 mpl_show(plt) 

1236 

1237 

1238def command_upgrade(args): 

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

1240 store_dirs = get_store_dirs(args) 

1241 try: 

1242 for store_dir in store_dirs: 

1243 store = gf.Store(store_dir) 

1244 nup = store.upgrade() 

1245 if nup == 0: 

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

1247 else: 

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

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

1250 

1251 except gf.StoreError as e: 

1252 die(e) 

1253 

1254 

1255def command_addref(args): 

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

1257 

1258 store_dirs = [] 

1259 filenames = [] 

1260 for arg in args: 

1261 if op.isdir(arg): 

1262 store_dirs.append(arg) 

1263 elif op.isfile(arg): 

1264 filenames.append(arg) 

1265 else: 

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

1267 

1268 if not store_dirs: 

1269 store_dirs.append('.') 

1270 

1271 references = [] 

1272 try: 

1273 for filename in args: 

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

1275 except ImportError: 

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

1277 

1278 if not references: 

1279 die('no references found') 

1280 

1281 for store_dir in store_dirs: 

1282 try: 

1283 store = gf.Store(store_dir) 

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

1285 for ref in references: 

1286 if ref.id in ids: 

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

1288 

1289 ids.append(ref.id) 

1290 store.config.references.append(ref) 

1291 

1292 store.save_config(make_backup=True) 

1293 

1294 except gf.StoreError as e: 

1295 die(e) 

1296 

1297 

1298def command_qc(args): 

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

1300 

1301 store_dir = get_store_dir(args) 

1302 

1303 try: 

1304 store = gf.Store(store_dir) 

1305 s = store.stats() 

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

1307 print('has empty records') 

1308 

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

1310 

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

1312 print('%s empty' % aname) 

1313 

1314 except gf.StoreError as e: 

1315 die(e) 

1316 

1317 

1318def command_report(args): 

1319 from pyrocko.fomosto.report import report_call 

1320 report_call.run_program(args) 

1321 

1322 

1323def main(args=None): 

1324 if args is None: 

1325 args = sys.argv[1:] 

1326 

1327 if len(args) < 1: 

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

1329 

1330 command = args.pop(0) 

1331 

1332 if command in subcommands: 

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

1334 

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

1336 if command == 'help' and args: 

1337 acommand = args[0] 

1338 if acommand in subcommands: 

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

1340 

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

1342 

1343 else: 

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

1345 

1346 

1347if __name__ == '__main__': 

1348 main()