Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/apps/fomosto.py: 60%

705 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 06:59 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Fomosto - Green's function database management. 

8''' 

9 

10import sys 

11import re 

12import os.path as op 

13import logging 

14import copy 

15import shutil 

16from optparse import OptionParser 

17 

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

19from pyrocko.gui.snuffler import marker 

20from pyrocko.util import mpl_show 

21 

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

23km = 1e3 

24 

25 

26def d2u(d): 

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

28 

29 

30subcommand_descriptions = { 

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

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

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

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

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

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

37 'view': 'view selected traces', 

38 'extract': 'extract selected traces', 

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

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

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

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

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

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

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

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

47 'server': 'run seismosizer server', 

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

49 'modelview': 'plot earthmodels', 

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

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

52 'qc': 'quality check', 

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

54} 

55 

56subcommand_usages = { 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

76 'modelview': 'modelview <selection>', 

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

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

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

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

81} 

82 

83subcommands = subcommand_descriptions.keys() 

84 

85program_name = 'fomosto' 

86 

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

88 

89Subcommands: 

90 

91 init %(init)s 

92 build %(build)s 

93 stats %(stats)s 

94 check %(check)s 

95 decimate %(decimate)s 

96 redeploy %(redeploy)s 

97 view %(view)s 

98 extract %(extract)s 

99 import %(import)s 

100 export %(export)s 

101 ttt %(ttt)s 

102 tttview %(tttview)s 

103 tttextract %(tttextract)s 

104 tttlsd %(tttlsd)s 

105 sat %(sat)s 

106 satview %(satview)s 

107 server %(server)s 

108 download %(download)s 

109 modelview %(modelview)s 

110 upgrade %(upgrade)s 

111 addref %(addref)s 

112 qc %(qc)s 

113 report %(report)s 

114 

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

116 

117 fomosto <subcommand> --help 

118 

119''' % d2u(subcommand_descriptions) 

120 

121 

122def add_common_options(parser): 

123 parser.add_option( 

124 '--loglevel', 

125 action='store', 

126 dest='loglevel', 

127 type='choice', 

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

129 default='info', 

130 help='set logger level to ' 

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

132 'Default is "%default".') 

133 

134 

135def process_common_options(options): 

136 util.setup_logging(program_name, options.loglevel) 

137 

138 

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

140 usage = subcommand_usages[command] 

141 descr = subcommand_descriptions[command] 

142 

143 if isinstance(usage, str): 

144 usage = [usage] 

145 

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

147 for s in usage[1:]: 

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

149 

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

151 

152 if details: 

153 description = description + ' %s' % details 

154 

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

156 parser.format_description = lambda formatter: description 

157 

158 if setup: 

159 setup(parser) 

160 

161 add_common_options(parser) 

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

163 process_common_options(options) 

164 return parser, options, args 

165 

166 

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

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

169 

170 

171def fomo_wrapper_module(name): 

172 try: 

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

174 raise ValueError('invalid name') 

175 

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

177 if len(words) == 2: 

178 name, variant = words 

179 else: 

180 name = words[0] 

181 variant = None 

182 

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

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

185 mod = __import__(modname, level=0) 

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

187 

188 except ValueError: 

189 die('invalid modelling code wrapper name') 

190 

191 except ImportError: 

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

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

194 

195 

196def command_init(args): 

197 

198 details = ''' 

199 

200Available modelling backends: 

201%s 

202 

203 More information at 

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

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

206 

207 parser, options, args = cl_parse( 

208 'init', args, 

209 details=details) 

210 

211 if len(args) == 0: 

212 sys.exit(parser.format_help()) 

213 

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

215 if len(args) != 3: 

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

217 

218 source_dir, dest_dir = args[1:] 

219 

220 try: 

221 source = gf.Store(source_dir) 

222 except gf.StoreError as e: 

223 die(e) 

224 

225 config = copy.deepcopy(source.config) 

226 config.derived_from_id = source.config.id 

227 try: 

228 config_filenames = gf.store.Store.create_editables( 

229 dest_dir, config=config) 

230 

231 except gf.StoreError as e: 

232 die(e) 

233 

234 try: 

235 dest = gf.Store(dest_dir) 

236 except gf.StoreError as e: 

237 die(e) 

238 

239 for k in source.extra_keys(): 

240 source_fn = source.get_extra_path(k) 

241 dest_fn = dest.get_extra_path(k) 

242 shutil.copyfile(source_fn, dest_fn) 

243 

244 logger.info( 

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

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

247 

248 logger.info( 

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

250 

251 else: 

252 if len(args) != 2: 

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

254 

255 (modelling_code_id, store_dir) = args 

256 

257 module, variant = fomo_wrapper_module(modelling_code_id) 

258 try: 

259 config_filenames = module.init(store_dir, variant) 

260 except gf.StoreError as e: 

261 die(e) 

262 

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

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

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

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

267 

268 

269def get_store_dir(args): 

270 if len(args) == 1: 

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

272 else: 

273 store_dir = op.abspath(op.curdir) 

274 

275 if not op.isdir(store_dir): 

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

277 

278 return store_dir 

279 

280 

281def get_store_dirs(args): 

282 if len(args) == 0: 

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

284 else: 

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

286 

287 for store_dir in store_dirs: 

288 if not op.isdir(store_dir): 

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

290 

291 return store_dirs 

292 

293 

294def command_build(args): 

295 

296 def setup(parser): 

297 parser.add_option( 

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

299 help='overwrite existing files') 

300 

301 parser.add_option( 

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

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

304 

305 parser.add_option( 

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

307 help='continue suspended build') 

308 

309 parser.add_option( 

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

311 help='process block number IBLOCK') 

312 

313 parser.add_option( 

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

315 help='process block number IBLOCK') 

316 

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

318 

319 store_dir = get_store_dir(args) 

320 try: 

321 if options.step is not None: 

322 step = options.step - 1 

323 else: 

324 step = None 

325 

326 if options.iblock is not None: 

327 iblock = options.iblock - 1 

328 else: 

329 iblock = None 

330 

331 store = gf.Store(store_dir) 

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

333 module.build( 

334 store_dir, 

335 force=options.force, 

336 nworkers=options.nworkers, continue_=options.continue_, 

337 step=step, 

338 iblock=iblock) 

339 

340 except gf.StoreError as e: 

341 die(e) 

342 

343 

344def command_stats(args): 

345 

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

347 store_dir = get_store_dir(args) 

348 

349 try: 

350 store = gf.Store(store_dir) 

351 s = store.stats() 

352 

353 except gf.StoreError as e: 

354 die(e) 

355 

356 for k in store.stats_keys: 

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

358 

359 

360def command_check(args): 

361 

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

363 store_dir = get_store_dir(args) 

364 

365 try: 

366 store = gf.Store(store_dir) 

367 problems = store.check(show_progress=True) 

368 if problems: 

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

370 

371 except gf.StoreError as e: 

372 die(e) 

373 

374 

375def load_config(fn): 

376 try: 

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

378 assert isinstance(config, gf.Config) 

379 

380 except Exception: 

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

382 

383 return config 

384 

385 

386def command_decimate(args): 

387 

388 def setup(parser): 

389 parser.add_option( 

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

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

392 

393 parser.add_option( 

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

395 help='overwrite existing files') 

396 

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

398 try: 

399 decimate = int(args.pop()) 

400 except Exception: 

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

402 

403 store_dir = get_store_dir(args) 

404 

405 config = None 

406 if options.config_fn: 

407 config = load_config(options.config_fn) 

408 

409 try: 

410 store = gf.Store(store_dir) 

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

412 show_progress=True) 

413 

414 except gf.StoreError as e: 

415 die(e) 

416 

417 

418def sindex(args): 

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

420 

421 

422def command_redeploy(args): 

423 

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

425 

426 if not len(args) == 2: 

427 sys.exit(parser.format_help()) 

428 

429 source_store_dir, dest_store_dir = args 

430 

431 try: 

432 source = gf.Store(source_store_dir) 

433 except gf.StoreError as e: 

434 die(e) 

435 

436 try: 

437 gf.store.Store.create_dependants(dest_store_dir) 

438 except gf.StoreError: 

439 pass 

440 

441 try: 

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

443 

444 except gf.StoreError as e: 

445 die(e) 

446 

447 show_progress = True 

448 

449 try: 

450 if show_progress: 

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

452 

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

454 try: 

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

456 dest.put(args, tr) 

457 

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

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

460 

461 except gf.store.StoreError as e: 

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

463 

464 if show_progress: 

465 pbar.update(i+1) 

466 

467 finally: 

468 if show_progress: 

469 pbar.finish() 

470 

471 

472def command_view(args): 

473 def setup(parser): 

474 parser.add_option( 

475 '--extract', 

476 dest='extract', 

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

478 help='specify which traces to show') 

479 

480 parser.add_option( 

481 '--show-phases', 

482 dest='showphases', 

483 default=None, 

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

485 help='add phase markers from ttt') 

486 

487 parser.add_option( 

488 '--opengl', 

489 dest='opengl', 

490 action='store_true', 

491 default=None, 

492 help='use OpenGL for drawing') 

493 

494 parser.add_option( 

495 '--no-opengl', 

496 dest='opengl', 

497 action='store_false', 

498 default=None, 

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

500 

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

502 

503 gdef = None 

504 if options.extract: 

505 try: 

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

507 except gf.meta.GridSpecError as e: 

508 die(e) 

509 

510 store_dirs = get_store_dirs(args) 

511 

512 alpha = 'abcdefghijklmnopqrstxyz'.upper() 

513 

514 markers = [] 

515 traces = [] 

516 

517 try: 

518 for istore, store_dir in enumerate(store_dirs): 

519 store = gf.Store(store_dir) 

520 if options.showphases == 'all': 

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

522 elif options.showphases is not None: 

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

524 

525 ii = 0 

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

527 gtr = store.get(args) 

528 

529 loc_code = '' 

530 if len(store_dirs) > 1: 

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

532 

533 if gtr: 

534 

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

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

537 tmin = gtr.deltat * gtr.itmin 

538 

539 tr = trace.Trace( 

540 '', 

541 sta_code, 

542 loc_code, 

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

544 ydata=gtr.data, 

545 deltat=gtr.deltat, 

546 tmin=tmin) 

547 

548 if options.showphases: 

549 for phasename in phasenames: 

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

551 if phase_tmin: 

552 m = marker.PhaseMarker( 

553 [('', 

554 sta_code, 

555 loc_code, 

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

557 phase_tmin, 

558 phase_tmin, 

559 0, 

560 phasename=phasename) 

561 markers.append(m) 

562 

563 traces.append(tr) 

564 

565 ii += 1 

566 

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

568 die(e) 

569 

570 def setup(win): 

571 viewer = win.get_view() 

572 viewer.menuitem_showboxes.setChecked(False) 

573 viewer.menuitem_colortraces.setChecked(True) 

574 viewer.menuitem_demean.setChecked(False) 

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

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

577 viewer.sortingmode_change() 

578 viewer.scalingmode_change() 

579 

580 trace.snuffle( 

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

582 

583 

584def command_extract(args): 

585 def setup(parser): 

586 parser.add_option( 

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

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

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

590 'Default is "mseed".') 

591 

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

593 parser.add_option( 

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

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

596 

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

598 try: 

599 sdef = args.pop() 

600 except Exception: 

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

602 

603 try: 

604 gdef = gf.meta.parse_grid_spec(sdef) 

605 except gf.meta.GridSpecError as e: 

606 die(e) 

607 

608 store_dir = get_store_dir(args) 

609 

610 extensions = { 

611 'mseed': 'mseed', 

612 'sac': 'sac', 

613 'text': 'txt', 

614 'yaff': 'yaff'} 

615 

616 try: 

617 store = gf.Store(store_dir) 

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

619 gtr = store.get(args) 

620 if gtr: 

621 tr = trace.Trace( 

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

623 ydata=gtr.data, 

624 deltat=gtr.deltat, 

625 tmin=gtr.deltat * gtr.itmin) 

626 

627 additional = dict( 

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

629 irecord=store.str_irecord(args), 

630 extension=extensions[options.format]) 

631 

632 io.save( 

633 tr, 

634 options.output_fn, 

635 format=options.format, 

636 additional=additional) 

637 

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

639 die(e) 

640 

641 

642def command_import(args): 

643 try: 

644 from tunguska import gfdb 

645 except ImportError: 

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

647 

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

649 

650 show_progress = True 

651 

652 if not len(args) == 2: 

653 sys.exit(parser.format_help()) 

654 

655 source_path, dest_store_dir = args 

656 

657 if op.isdir(source_path): 

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

659 

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

661 

662 db = gfdb.Gfdb(source_path) 

663 

664 config = gf.meta.ConfigTypeA( 

665 id='imported_gfs', 

666 distance_min=db.firstx, 

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

668 distance_delta=db.dx, 

669 source_depth_min=db.firstz, 

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

671 source_depth_delta=db.dz, 

672 sample_rate=1.0/db.dt, 

673 ncomponents=db.ng 

674 ) 

675 

676 try: 

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

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

679 try: 

680 if show_progress: 

681 pbar = util.progressbar( 

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

683 

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

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

686 traces = db.get_traces_pyrocko(distance, source_depth) 

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

688 for ig in range(db.ng): 

689 if ig in ig_to_trace: 

690 tr = ig_to_trace[ig] 

691 gf_tr = gf.store.GFTrace( 

692 tr.get_ydata(), 

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

694 tr.deltat) 

695 

696 else: 

697 gf_tr = gf.store.Zero 

698 

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

700 

701 if show_progress: 

702 pbar.update(i+1) 

703 

704 finally: 

705 if show_progress: 

706 pbar.finish() 

707 

708 dest.close() 

709 

710 except gf.StoreError as e: 

711 die(e) 

712 

713 

714def command_export(args): 

715 from subprocess import Popen, PIPE 

716 

717 try: 

718 from tunguska import gfdb 

719 except ImportError as err: 

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

721 

722 def setup(parser): 

723 parser.add_option( 

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

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

726 

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

728 

729 show_progress = True 

730 

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

732 sys.exit(parser.format_help()) 

733 

734 target_path = args.pop() 

735 if op.isdir(target_path): 

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

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

738 

739 source_store_dir = get_store_dir(args) 

740 

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

742 config = source.config 

743 

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

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

746 

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

748 die('destation already exists') 

749 

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

751 'gfdb_build', 

752 target_path, 

753 options.nchunks, 

754 config.ndistances, 

755 config.nsource_depths, 

756 config.ncomponents, 

757 config.deltat, 

758 config.distance_delta, 

759 config.source_depth_delta, 

760 config.distance_min, 

761 config.source_depth_min]] 

762 

763 p = Popen(cmd, stdin=PIPE) 

764 p.communicate() 

765 

766 out_db = gfdb.Gfdb(target_path) 

767 

768 try: 

769 if show_progress: 

770 pbar = util.progressbar( 

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

772 

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

774 

775 data_out = [] 

776 for ig in range(config.ncomponents): 

777 try: 

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

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

780 

781 except gf.store.StoreError as e: 

782 logger.warning( 

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

784 data_out.append(None) 

785 

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

787 # time 

788 tmins = [ 

789 entry[0][0] 

790 for entry in data_out 

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

792 

793 if tmins: 

794 tmin = min(tmins) 

795 for entry in data_out: 

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

797 entry[0].resize(1) 

798 entry[1].resize(1) 

799 entry[0][0] = tmin 

800 entry[1][0] = 0.0 

801 

802 out_db.put_traces_slow(x, z, data_out) 

803 

804 if show_progress: 

805 pbar.update(i+1) 

806 

807 source.close() 

808 

809 finally: 

810 if show_progress: 

811 pbar.finish() 

812 

813 

814def phasedef_or_horvel(x): 

815 try: 

816 return float(x) 

817 except ValueError: 

818 return cake.PhaseDef(x) 

819 

820 

821def mkp(s): 

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

823 

824 

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

826 import numpy as num 

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

828 

829 plt = mpl_init() 

830 

831 np = 1 

832 store_dir = get_store_dir(args) 

833 for phase_id in phase_ids: 

834 try: 

835 store = gf.Store(store_dir) 

836 

837 if options.receiver_depth is not None: 

838 receiver_depth = options.receiver_depth * 1000.0 

839 else: 

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

841 receiver_depth = store.config.receiver_depth 

842 

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

844 receiver_depth = store.config.receiver_depth_min 

845 

846 else: 

847 receiver_depth = 0.0 

848 

849 phase = store.get_stored_phase(phase_id, attribute) 

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

851 labelspace(axes) 

852 xscaled(1. / km, axes) 

853 yscaled(1. / km, axes) 

854 x = None 

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

856 x = (receiver_depth, None, None) 

857 

858 phase.plot_2d(axes, x=x) 

859 axes.set_title(phase_id) 

860 np += 1 

861 except gf.StoreError as e: 

862 die(e) 

863 

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

865 num_d = 100 

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

867 store.config.distance_max, 

868 num_d) 

869 

870 if options.source_depth is not None: 

871 source_depth = options.source_depth * km 

872 else: 

873 source_depth = store.config.source_depth_min + ( 

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

875 

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

877 attribute_vals = num.empty(num_d) 

878 for phase_id in phase_ids: 

879 attribute_vals[:] = num.NAN 

880 for i, d in enumerate(distances): 

881 if attribute == 'phase': 

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

883 ylabel = 'TT [s]' 

884 else: 

885 attribute_vals[i] = store.get_stored_attribute( 

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

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

888 

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

890 

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

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

893 axes.set_ylabel(ylabel) 

894 axes.legend() 

895 

896 plt.tight_layout() 

897 mpl_show(plt) 

898 

899 

900def command_ttt(args): 

901 def setup(parser): 

902 parser.add_option( 

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

904 help='overwrite existing files') 

905 

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

907 

908 store_dir = get_store_dir(args) 

909 try: 

910 store = gf.Store(store_dir) 

911 store.make_travel_time_tables(force=options.force) 

912 

913 except gf.StoreError as e: 

914 die(e) 

915 

916 

917def command_tttview(args): 

918 

919 def setup(parser): 

920 parser.add_option( 

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

922 help='Source depth in km') 

923 

924 parser.add_option( 

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

926 help='Receiver depth in km') 

927 

928 parser, options, args = cl_parse( 

929 'tttview', args, setup=setup, 

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

931 

932 try: 

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

934 except Exception: 

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

936 

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

938 

939 

940def command_sat(args): 

941 def setup(parser): 

942 parser.add_option( 

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

944 help='overwrite existing files') 

945 

946 parser.add_option( 

947 '--attribute', 

948 action='store', 

949 dest='attribute', 

950 type='choice', 

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

952 default='takeoff_angle', 

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

954 

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

956 

957 store_dir = get_store_dir(args) 

958 try: 

959 store = gf.Store(store_dir) 

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

961 

962 except gf.StoreError as e: 

963 die(e) 

964 

965 

966def command_satview(args): 

967 

968 def setup(parser): 

969 parser.add_option( 

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

971 help='Source depth in km') 

972 

973 parser.add_option( 

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

975 help='Receiver depth in km') 

976 

977 parser.add_option( 

978 '--attribute', 

979 action='store', 

980 dest='attribute', 

981 type='choice', 

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

983 default='takeoff_angle', 

984 help='view selected ray attribute.') 

985 

986 parser, options, args = cl_parse( 

987 'satview', args, setup=setup, 

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

989 

990 try: 

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

992 except Exception: 

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

994 

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

996 

997 stored_attribute_table_plots( 

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

999 

1000 

1001def command_tttextract(args): 

1002 def setup(parser): 

1003 parser.add_option( 

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

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

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

1007 

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

1009 try: 

1010 sdef = args.pop() 

1011 except Exception: 

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

1013 

1014 try: 

1015 sphase = args.pop() 

1016 except Exception: 

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

1018 

1019 try: 

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

1021 except gf.meta.InvalidTimingSpecification: 

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

1023 

1024 try: 

1025 gdef = gf.meta.parse_grid_spec(sdef) 

1026 except gf.meta.GridSpecError as e: 

1027 die(e) 

1028 

1029 store_dir = get_store_dir(args) 

1030 

1031 try: 

1032 store = gf.Store(store_dir) 

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

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

1035 for phase in phases: 

1036 t = store.t(phase, args) 

1037 if t is not None: 

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

1039 else: 

1040 s.append('nan') 

1041 

1042 if options.output_fn: 

1043 d = dict( 

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

1045 extension='txt') 

1046 

1047 fn = options.output_fn % d 

1048 util.ensuredirs(fn) 

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

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

1051 f.write('\n') 

1052 else: 

1053 print(' '.join(s)) 

1054 

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

1056 die(e) 

1057 

1058 

1059def command_tttlsd(args): 

1060 

1061 def setup(parser): 

1062 pass 

1063 

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

1065 

1066 try: 

1067 sphase_ids = args.pop() 

1068 except Exception: 

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

1070 

1071 try: 

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

1073 except gf.meta.InvalidTimingSpecification: 

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

1075 

1076 store_dir = get_store_dir(args) 

1077 

1078 try: 

1079 store = gf.Store(store_dir) 

1080 for phase_id in phase_ids: 

1081 store.fix_ttt_holes(phase_id) 

1082 

1083 except gf.StoreError as e: 

1084 die(e) 

1085 

1086 

1087def command_server(args): 

1088 from pyrocko.gf import server 

1089 

1090 def setup(parser): 

1091 parser.add_option( 

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

1093 help='serve on port PORT') 

1094 

1095 parser.add_option( 

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

1097 help='serve on ip address IP') 

1098 

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

1100 

1101 engine = gf.LocalEngine(store_superdirs=args) 

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

1103 

1104 

1105def command_download(args): 

1106 from pyrocko.gf import ws 

1107 

1108 details = ''' 

1109 

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

1111 

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

1113''' 

1114 

1115 def setup(parser): 

1116 parser.add_option( 

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

1118 help='overwrite existing files') 

1119 

1120 parser, options, args = cl_parse( 

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

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

1123 sys.exit(parser.format_help()) 

1124 

1125 if len(args) == 2: 

1126 site, store_id = args 

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

1128 die('invalid store ID') 

1129 else: 

1130 site, store_id = args[0], None 

1131 

1132 if site not in gf.ws.g_site_abbr: 

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

1134 site = 'http://' + site 

1135 

1136 try: 

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

1138 except ws.DownloadError as e: 

1139 die(str(e)) 

1140 

1141 

1142def command_modelview(args): 

1143 

1144 import matplotlib.pyplot as plt 

1145 import numpy as num 

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

1147 mpl_init() 

1148 

1149 neat_labels = { 

1150 'vp': '$v_p$', 

1151 'vs': '$v_s$', 

1152 'qp': '$Q_p$', 

1153 'qs': '$Q_s$', 

1154 'rho': '$\\rho$'} 

1155 

1156 def setup(parser): 

1157 parser.add_option( 

1158 '--parameters', dest='parameters', 

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

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

1161 

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

1163 

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

1165 

1166 store_dirs = get_store_dirs(args) 

1167 

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

1169 

1170 fig, axes = plt.subplots(1, 

1171 len(parameters), 

1172 sharey=True, 

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

1174 

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

1176 axes = [axes] 

1177 

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

1179 

1180 for store_id in store_dirs: 

1181 try: 

1182 store = gf.Store(store_id) 

1183 mod = store.config.earthmodel_1d 

1184 z = mod.profile('z') 

1185 

1186 for p in parameters: 

1187 ax = axes[p] 

1188 labelspace(ax) 

1189 if '/' in p: 

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

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

1192 else: 

1193 profile = mod.profile(p) 

1194 

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

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

1197 xscaled(1./1000., ax) 

1198 

1199 yscaled(1./km, ax) 

1200 

1201 except gf.StoreError as e: 

1202 die(e) 

1203 

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

1205 ax.grid() 

1206 if p in neat_labels: 

1207 lab = neat_labels[p] 

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

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

1210 else: 

1211 lab = p 

1212 

1213 ax.set_title(lab, y=1.02) 

1214 

1215 if p in units: 

1216 lab += ' ' + units[p] 

1217 

1218 ax.autoscale() 

1219 ax.set_xlabel(lab) 

1220 

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

1222 

1223 handles, labels = ax.get_legend_handles_labels() 

1224 

1225 if len(store_dirs) > 1: 

1226 ax.legend( 

1227 handles, 

1228 labels, 

1229 bbox_to_anchor=(0.5, 0.12), 

1230 bbox_transform=fig.transFigure, 

1231 ncol=3, 

1232 loc='upper center', 

1233 fancybox=True) 

1234 

1235 plt.subplots_adjust(bottom=0.22, 

1236 wspace=0.05) 

1237 

1238 mpl_show(plt) 

1239 

1240 

1241def command_upgrade(args): 

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

1243 store_dirs = get_store_dirs(args) 

1244 try: 

1245 for store_dir in store_dirs: 

1246 store = gf.Store(store_dir) 

1247 nup = store.upgrade() 

1248 if nup == 0: 

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

1250 else: 

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

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

1253 

1254 except gf.StoreError as e: 

1255 die(e) 

1256 

1257 

1258def command_addref(args): 

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

1260 

1261 store_dirs = [] 

1262 filenames = [] 

1263 for arg in args: 

1264 if op.isdir(arg): 

1265 store_dirs.append(arg) 

1266 elif op.isfile(arg): 

1267 filenames.append(arg) 

1268 else: 

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

1270 

1271 if not store_dirs: 

1272 store_dirs.append('.') 

1273 

1274 references = [] 

1275 try: 

1276 for filename in args: 

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

1278 except ImportError: 

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

1280 

1281 if not references: 

1282 die('no references found') 

1283 

1284 for store_dir in store_dirs: 

1285 try: 

1286 store = gf.Store(store_dir) 

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

1288 for ref in references: 

1289 if ref.id in ids: 

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

1291 

1292 ids.append(ref.id) 

1293 store.config.references.append(ref) 

1294 

1295 store.save_config(make_backup=True) 

1296 

1297 except gf.StoreError as e: 

1298 die(e) 

1299 

1300 

1301def command_qc(args): 

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

1303 

1304 store_dir = get_store_dir(args) 

1305 

1306 try: 

1307 store = gf.Store(store_dir) 

1308 s = store.stats() 

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

1310 print('has empty records') 

1311 

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

1313 

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

1315 print('%s empty' % aname) 

1316 

1317 except gf.StoreError as e: 

1318 die(e) 

1319 

1320 

1321def command_report(args): 

1322 from pyrocko.fomosto.report import report_call 

1323 report_call.run_program(args) 

1324 

1325 

1326def main(args=None): 

1327 ''' 

1328 CLI entry point for Pyrocko's ``fomosto`` app. 

1329 ''' 

1330 if args is None: 

1331 args = sys.argv[1:] 

1332 

1333 if len(args) < 1: 

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

1335 

1336 command = args.pop(0) 

1337 

1338 if command in subcommands: 

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

1340 

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

1342 if command == 'help' and args: 

1343 acommand = args[0] 

1344 if acommand in subcommands: 

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

1346 

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

1348 

1349 else: 

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

1351 

1352 

1353if __name__ == '__main__': 

1354 main()