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

727 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +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 'trans2rot': 'create rotational from translational store', 

55 'fd2trans': 'extract elastic10 from elastic10_fd store', 

56} 

57 

58subcommand_usages = { 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

78 'modelview': 'modelview <selection>', 

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

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

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

82 'report': 'report <subcommand> <arguments>... [options]', 

83 'trans2rot': 'trans2rot <source> <destination>', 

84 'fd2trans': 'fd2trans <source> <destination>', 

85} 

86 

87subcommands = subcommand_descriptions.keys() 

88 

89program_name = 'fomosto' 

90 

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

92 

93Subcommands: 

94 

95 init %(init)s 

96 build %(build)s 

97 stats %(stats)s 

98 check %(check)s 

99 decimate %(decimate)s 

100 redeploy %(redeploy)s 

101 view %(view)s 

102 extract %(extract)s 

103 import %(import)s 

104 export %(export)s 

105 ttt %(ttt)s 

106 tttview %(tttview)s 

107 tttextract %(tttextract)s 

108 tttlsd %(tttlsd)s 

109 sat %(sat)s 

110 satview %(satview)s 

111 server %(server)s 

112 download %(download)s 

113 modelview %(modelview)s 

114 upgrade %(upgrade)s 

115 addref %(addref)s 

116 qc %(qc)s 

117 report %(report)s 

118 trans2rot %(trans2rot)s 

119 fd2trans %(fd2trans)s 

120 

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

122 

123 fomosto <subcommand> --help 

124 

125''' % d2u(subcommand_descriptions) 

126 

127 

128def add_common_options(parser): 

129 parser.add_option( 

130 '--loglevel', 

131 action='store', 

132 dest='loglevel', 

133 type='choice', 

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

135 default='info', 

136 help='set logger level to ' 

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

138 'Default is "%default".') 

139 

140 

141def process_common_options(options): 

142 util.setup_logging(program_name, options.loglevel) 

143 

144 

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

146 usage = subcommand_usages[command] 

147 descr = subcommand_descriptions[command] 

148 

149 if isinstance(usage, str): 

150 usage = [usage] 

151 

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

153 for s in usage[1:]: 

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

155 

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

157 

158 if details: 

159 description = description + ' %s' % details 

160 

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

162 parser.format_description = lambda formatter: description 

163 

164 if setup: 

165 setup(parser) 

166 

167 add_common_options(parser) 

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

169 process_common_options(options) 

170 return parser, options, args 

171 

172 

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

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

175 

176 

177def fomo_wrapper_module(name): 

178 try: 

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

180 raise ValueError('invalid name') 

181 

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

183 if len(words) == 2: 

184 name, variant = words 

185 else: 

186 name = words[0] 

187 variant = None 

188 

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

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

191 mod = __import__(modname, level=0) 

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

193 

194 except ValueError: 

195 die('invalid modelling code wrapper name') 

196 

197 except ImportError: 

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

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

200 

201 

202def command_init(args): 

203 

204 details = ''' 

205 

206Available modelling backends: 

207%s 

208 

209 More information at 

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

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

212 

213 parser, options, args = cl_parse( 

214 'init', args, 

215 details=details) 

216 

217 if len(args) == 0: 

218 sys.exit(parser.format_help()) 

219 

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

221 if len(args) != 3: 

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

223 

224 source_dir, dest_dir = args[1:] 

225 

226 try: 

227 source = gf.Store(source_dir) 

228 except gf.StoreError as e: 

229 die(e) 

230 

231 config = copy.deepcopy(source.config) 

232 config.derived_from_id = source.config.id 

233 try: 

234 config_filenames = gf.store.Store.create_editables( 

235 dest_dir, config=config) 

236 

237 except gf.StoreError as e: 

238 die(e) 

239 

240 try: 

241 dest = gf.Store(dest_dir) 

242 except gf.StoreError as e: 

243 die(e) 

244 

245 for k in source.extra_keys(): 

246 source_fn = source.get_extra_path(k) 

247 dest_fn = dest.get_extra_path(k) 

248 shutil.copyfile(source_fn, dest_fn) 

249 

250 logger.info( 

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

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

253 

254 logger.info( 

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

256 

257 else: 

258 if len(args) != 2: 

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

260 

261 (modelling_code_id, store_dir) = args 

262 

263 module, variant = fomo_wrapper_module(modelling_code_id) 

264 try: 

265 config_filenames = module.init(store_dir, variant) 

266 except gf.StoreError as e: 

267 die(e) 

268 

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

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

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

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

273 

274 

275def get_store_dir(args): 

276 if len(args) == 1: 

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

278 else: 

279 store_dir = op.abspath(op.curdir) 

280 

281 if not op.isdir(store_dir): 

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

283 

284 return store_dir 

285 

286 

287def get_store_dirs(args): 

288 if len(args) == 0: 

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

290 else: 

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

292 

293 for store_dir in store_dirs: 

294 if not op.isdir(store_dir): 

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

296 

297 return store_dirs 

298 

299 

300def command_build(args): 

301 

302 def setup(parser): 

303 parser.add_option( 

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

305 help='overwrite existing files') 

306 

307 parser.add_option( 

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

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

310 

311 parser.add_option( 

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

313 help='continue suspended build') 

314 

315 parser.add_option( 

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

317 help='process block number IBLOCK') 

318 

319 parser.add_option( 

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

321 help='process block number IBLOCK') 

322 

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

324 

325 store_dir = get_store_dir(args) 

326 try: 

327 if options.step is not None: 

328 step = options.step - 1 

329 else: 

330 step = None 

331 

332 if options.iblock is not None: 

333 iblock = options.iblock - 1 

334 else: 

335 iblock = None 

336 

337 store = gf.Store(store_dir) 

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

339 module.build( 

340 store_dir, 

341 force=options.force, 

342 nworkers=options.nworkers, continue_=options.continue_, 

343 step=step, 

344 iblock=iblock) 

345 

346 except gf.StoreError as e: 

347 die(e) 

348 

349 

350def command_stats(args): 

351 

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

353 store_dir = get_store_dir(args) 

354 

355 try: 

356 store = gf.Store(store_dir) 

357 s = store.stats() 

358 

359 except gf.StoreError as e: 

360 die(e) 

361 

362 for k in store.stats_keys: 

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

364 

365 

366def command_check(args): 

367 

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

369 store_dir = get_store_dir(args) 

370 

371 try: 

372 store = gf.Store(store_dir) 

373 problems = store.check(show_progress=True) 

374 if problems: 

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

376 

377 except gf.StoreError as e: 

378 die(e) 

379 

380 

381def load_config(fn): 

382 try: 

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

384 assert isinstance(config, gf.Config) 

385 

386 except Exception: 

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

388 

389 return config 

390 

391 

392def command_decimate(args): 

393 

394 def setup(parser): 

395 parser.add_option( 

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

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

398 

399 parser.add_option( 

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

401 help='overwrite existing files') 

402 

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

404 try: 

405 decimate = int(args.pop()) 

406 except Exception: 

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

408 

409 store_dir = get_store_dir(args) 

410 

411 config = None 

412 if options.config_fn: 

413 config = load_config(options.config_fn) 

414 

415 try: 

416 store = gf.Store(store_dir) 

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

418 show_progress=True) 

419 

420 except gf.StoreError as e: 

421 die(e) 

422 

423 

424def sindex(args): 

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

426 

427 

428def command_redeploy(args): 

429 

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

431 

432 if not len(args) == 2: 

433 sys.exit(parser.format_help()) 

434 

435 source_store_dir, dest_store_dir = args 

436 

437 try: 

438 source = gf.Store(source_store_dir) 

439 except gf.StoreError as e: 

440 die(e) 

441 

442 try: 

443 gf.store.Store.create_dependants(dest_store_dir) 

444 except gf.StoreError: 

445 pass 

446 

447 try: 

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

449 

450 except gf.StoreError as e: 

451 die(e) 

452 

453 show_progress = True 

454 

455 try: 

456 if show_progress: 

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

458 

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

460 try: 

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

462 dest.put(args, tr) 

463 

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

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

466 

467 except gf.store.StoreError as e: 

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

469 

470 if show_progress: 

471 pbar.update(i+1) 

472 

473 finally: 

474 if show_progress: 

475 pbar.finish() 

476 

477 

478def command_trans2rot(args): 

479 from pyrocko.fomosto.elastic10_to_rotational8 \ 

480 import elastic10_to_rotational8 

481 

482 details = ''' 

483 

484Create a rotational GF store from a translational one using finite differences 

485approximations (elastic10 or elastic10_fd to rotational8). 

486 

487If the source store is of type A and receivers at the surface, the free 

488surface condition is used to replace the vertical derivatives with horizontal 

489ones. If the input is of type B, the complete set of derivatives are used. It 

490is important that the spatial sampling of the source store is high enough so 

491that the finite difference approximations are sufficiently accurate. 

492 

493If the destination store directory does not exist, it is created and the extent 

494and spacing of the source store is used for the rotational store. If it does 

495exist and contains a valid `config`, it is used in place in order to allow for 

496different input and output extents and spacings. The store must be empty in the 

497latter case (delete the `index` and `traces` files if needed). 

498''' 

499 

500 parser, options, args = cl_parse( 

501 'trans2rot', args, details=details) 

502 

503 if not len(args) == 2: 

504 sys.exit(parser.format_help()) 

505 

506 source_store_dir, dest_store_dir = args 

507 

508 try: 

509 elastic10_to_rotational8(source_store_dir, dest_store_dir) 

510 except gf.store.StoreError as e: 

511 die(e) 

512 

513 

514def command_fd2trans(args): 

515 from pyrocko.fomosto.elastic10_fd_to_elastic10 \ 

516 import elastic10_fd_to_elastic10 

517 

518 details = ''' 

519 

520Extract regular GF store (elastic10) from a finite difference support GF store 

521(elastic10_fd). 

522''' 

523 

524 parser, options, args = cl_parse( 

525 'fd2trans', args, details=details) 

526 

527 if not len(args) == 2: 

528 sys.exit(parser.format_help()) 

529 

530 source_store_dir, dest_store_dir = args 

531 

532 try: 

533 elastic10_fd_to_elastic10(source_store_dir, dest_store_dir) 

534 except gf.store.StoreError as e: 

535 die(e) 

536 

537 

538def command_view(args): 

539 def setup(parser): 

540 parser.add_option( 

541 '--extract', 

542 dest='extract', 

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

544 help='specify which traces to show') 

545 

546 parser.add_option( 

547 '--show-phases', 

548 dest='showphases', 

549 default=None, 

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

551 help='add phase markers from ttt') 

552 

553 parser.add_option( 

554 '--opengl', 

555 dest='opengl', 

556 action='store_true', 

557 default=None, 

558 help='use OpenGL for drawing') 

559 

560 parser.add_option( 

561 '--no-opengl', 

562 dest='opengl', 

563 action='store_false', 

564 default=None, 

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

566 

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

568 

569 gdef = None 

570 if options.extract: 

571 try: 

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

573 except gf.meta.GridSpecError as e: 

574 die(e) 

575 

576 store_dirs = get_store_dirs(args) 

577 

578 alpha = 'abcdefghijklmnopqrstxyz'.upper() 

579 

580 markers = [] 

581 traces = [] 

582 

583 try: 

584 for istore, store_dir in enumerate(store_dirs): 

585 store = gf.Store(store_dir) 

586 if options.showphases == 'all': 

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

588 elif options.showphases is not None: 

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

590 

591 ii = 0 

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

593 gtr = store.get(args) 

594 

595 loc_code = '' 

596 if len(store_dirs) > 1: 

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

598 

599 if gtr: 

600 

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

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

603 tmin = gtr.deltat * gtr.itmin 

604 

605 tr = trace.Trace( 

606 '', 

607 sta_code, 

608 loc_code, 

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

610 ydata=gtr.data, 

611 deltat=gtr.deltat, 

612 tmin=tmin) 

613 

614 if options.showphases: 

615 for phasename in phasenames: 

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

617 if phase_tmin: 

618 m = marker.PhaseMarker( 

619 [('', 

620 sta_code, 

621 loc_code, 

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

623 phase_tmin, 

624 phase_tmin, 

625 0, 

626 phasename=phasename) 

627 markers.append(m) 

628 

629 traces.append(tr) 

630 

631 ii += 1 

632 

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

634 die(e) 

635 

636 def setup(win): 

637 viewer = win.get_view() 

638 viewer.menuitem_showboxes.setChecked(False) 

639 viewer.menuitem_colortraces.setChecked(True) 

640 viewer.menuitem_demean.setChecked(False) 

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

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

643 viewer.sortingmode_change() 

644 viewer.scalingmode_change() 

645 

646 trace.snuffle( 

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

648 

649 

650def command_extract(args): 

651 def setup(parser): 

652 parser.add_option( 

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

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

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

656 'Default is "mseed".') 

657 

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

659 parser.add_option( 

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

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

662 

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

664 try: 

665 sdef = args.pop() 

666 except Exception: 

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

668 

669 try: 

670 gdef = gf.meta.parse_grid_spec(sdef) 

671 except gf.meta.GridSpecError as e: 

672 die(e) 

673 

674 store_dir = get_store_dir(args) 

675 

676 extensions = { 

677 'mseed': 'mseed', 

678 'sac': 'sac', 

679 'text': 'txt', 

680 'yaff': 'yaff'} 

681 

682 try: 

683 store = gf.Store(store_dir) 

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

685 gtr = store.get(args) 

686 if gtr: 

687 tr = trace.Trace( 

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

689 ydata=gtr.data, 

690 deltat=gtr.deltat, 

691 tmin=gtr.deltat * gtr.itmin) 

692 

693 additional = dict( 

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

695 irecord=store.str_irecord(args), 

696 extension=extensions[options.format]) 

697 

698 io.save( 

699 tr, 

700 options.output_fn, 

701 format=options.format, 

702 additional=additional) 

703 

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

705 die(e) 

706 

707 

708def command_import(args): 

709 try: 

710 from tunguska import gfdb 

711 except ImportError: 

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

713 

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

715 

716 show_progress = True 

717 

718 if not len(args) == 2: 

719 sys.exit(parser.format_help()) 

720 

721 source_path, dest_store_dir = args 

722 

723 if op.isdir(source_path): 

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

725 

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

727 

728 db = gfdb.Gfdb(source_path) 

729 

730 config = gf.meta.ConfigTypeA( 

731 id='imported_gfs', 

732 distance_min=db.firstx, 

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

734 distance_delta=db.dx, 

735 source_depth_min=db.firstz, 

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

737 source_depth_delta=db.dz, 

738 sample_rate=1.0/db.dt, 

739 ncomponents=db.ng 

740 ) 

741 

742 try: 

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

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

745 try: 

746 if show_progress: 

747 pbar = util.progressbar( 

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

749 

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

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

752 traces = db.get_traces_pyrocko(distance, source_depth) 

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

754 for ig in range(db.ng): 

755 if ig in ig_to_trace: 

756 tr = ig_to_trace[ig] 

757 gf_tr = gf.store.GFTrace( 

758 tr.get_ydata(), 

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

760 tr.deltat) 

761 

762 else: 

763 gf_tr = gf.store.Zero 

764 

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

766 

767 if show_progress: 

768 pbar.update(i+1) 

769 

770 finally: 

771 if show_progress: 

772 pbar.finish() 

773 

774 dest.close() 

775 

776 except gf.StoreError as e: 

777 die(e) 

778 

779 

780def command_export(args): 

781 from subprocess import Popen, PIPE 

782 

783 try: 

784 from tunguska import gfdb 

785 except ImportError as err: 

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

787 

788 def setup(parser): 

789 parser.add_option( 

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

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

792 

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

794 

795 show_progress = True 

796 

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

798 sys.exit(parser.format_help()) 

799 

800 target_path = args.pop() 

801 if op.isdir(target_path): 

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

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

804 

805 source_store_dir = get_store_dir(args) 

806 

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

808 config = source.config 

809 

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

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

812 

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

814 die('destation already exists') 

815 

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

817 'gfdb_build', 

818 target_path, 

819 options.nchunks, 

820 config.ndistances, 

821 config.nsource_depths, 

822 config.ncomponents, 

823 config.deltat, 

824 config.distance_delta, 

825 config.source_depth_delta, 

826 config.distance_min, 

827 config.source_depth_min]] 

828 

829 p = Popen(cmd, stdin=PIPE) 

830 p.communicate() 

831 

832 out_db = gfdb.Gfdb(target_path) 

833 

834 try: 

835 if show_progress: 

836 pbar = util.progressbar( 

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

838 

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

840 

841 data_out = [] 

842 for ig in range(config.ncomponents): 

843 try: 

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

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

846 

847 except gf.store.StoreError as e: 

848 logger.warning( 

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

850 data_out.append(None) 

851 

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

853 # time 

854 tmins = [ 

855 entry[0][0] 

856 for entry in data_out 

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

858 

859 if tmins: 

860 tmin = min(tmins) 

861 for entry in data_out: 

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

863 entry[0].resize(1) 

864 entry[1].resize(1) 

865 entry[0][0] = tmin 

866 entry[1][0] = 0.0 

867 

868 out_db.put_traces_slow(x, z, data_out) 

869 

870 if show_progress: 

871 pbar.update(i+1) 

872 

873 source.close() 

874 

875 finally: 

876 if show_progress: 

877 pbar.finish() 

878 

879 

880def phasedef_or_horvel(x): 

881 try: 

882 return float(x) 

883 except ValueError: 

884 return cake.PhaseDef(x) 

885 

886 

887def mkp(s): 

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

889 

890 

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

892 import numpy as num 

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

894 

895 plt = mpl_init() 

896 

897 np = 1 

898 store_dir = get_store_dir(args) 

899 for phase_id in phase_ids: 

900 try: 

901 store = gf.Store(store_dir) 

902 

903 if options.receiver_depth is not None: 

904 receiver_depth = options.receiver_depth * 1000.0 

905 else: 

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

907 receiver_depth = store.config.receiver_depth 

908 

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

910 receiver_depth = store.config.receiver_depth_min 

911 

912 else: 

913 receiver_depth = 0.0 

914 

915 phase = store.get_stored_phase(phase_id, attribute) 

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

917 labelspace(axes) 

918 xscaled(1. / km, axes) 

919 yscaled(1. / km, axes) 

920 x = None 

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

922 x = (receiver_depth, None, None) 

923 

924 phase.plot_2d(axes, x=x) 

925 axes.set_title(phase_id) 

926 np += 1 

927 except gf.StoreError as e: 

928 die(e) 

929 

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

931 num_d = 100 

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

933 store.config.distance_max, 

934 num_d) 

935 

936 if options.source_depth is not None: 

937 source_depth = options.source_depth * km 

938 else: 

939 source_depth = store.config.source_depth_min + ( 

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

941 

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

943 attribute_vals = num.empty(num_d) 

944 for phase_id in phase_ids: 

945 attribute_vals[:] = num.NAN 

946 for i, d in enumerate(distances): 

947 if attribute == 'phase': 

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

949 ylabel = 'TT [s]' 

950 else: 

951 attribute_vals[i] = store.get_stored_attribute( 

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

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

954 

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

956 

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

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

959 axes.set_ylabel(ylabel) 

960 axes.legend() 

961 

962 plt.tight_layout() 

963 mpl_show(plt) 

964 

965 

966def command_ttt(args): 

967 def setup(parser): 

968 parser.add_option( 

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

970 help='overwrite existing files') 

971 

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

973 

974 store_dir = get_store_dir(args) 

975 try: 

976 store = gf.Store(store_dir) 

977 store.make_travel_time_tables(force=options.force) 

978 

979 except gf.StoreError as e: 

980 die(e) 

981 

982 

983def command_tttview(args): 

984 

985 def setup(parser): 

986 parser.add_option( 

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

988 help='Source depth in km') 

989 

990 parser.add_option( 

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

992 help='Receiver depth in km') 

993 

994 parser, options, args = cl_parse( 

995 'tttview', args, setup=setup, 

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

997 

998 try: 

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

1000 except Exception: 

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

1002 

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

1004 

1005 

1006def command_sat(args): 

1007 def setup(parser): 

1008 parser.add_option( 

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

1010 help='overwrite existing files') 

1011 

1012 parser.add_option( 

1013 '--attribute', 

1014 action='store', 

1015 dest='attribute', 

1016 type='choice', 

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

1018 default='takeoff_angle', 

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

1020 

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

1022 

1023 store_dir = get_store_dir(args) 

1024 try: 

1025 store = gf.Store(store_dir) 

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

1027 

1028 except gf.StoreError as e: 

1029 die(e) 

1030 

1031 

1032def command_satview(args): 

1033 

1034 def setup(parser): 

1035 parser.add_option( 

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

1037 help='Source depth in km') 

1038 

1039 parser.add_option( 

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

1041 help='Receiver depth in km') 

1042 

1043 parser.add_option( 

1044 '--attribute', 

1045 action='store', 

1046 dest='attribute', 

1047 type='choice', 

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

1049 default='takeoff_angle', 

1050 help='view selected ray attribute.') 

1051 

1052 parser, options, args = cl_parse( 

1053 'satview', args, setup=setup, 

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

1055 

1056 try: 

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

1058 except Exception: 

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

1060 

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

1062 

1063 stored_attribute_table_plots( 

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

1065 

1066 

1067def command_tttextract(args): 

1068 def setup(parser): 

1069 parser.add_option( 

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

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

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

1073 

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

1075 try: 

1076 sdef = args.pop() 

1077 except Exception: 

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

1079 

1080 try: 

1081 sphase = args.pop() 

1082 except Exception: 

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

1084 

1085 try: 

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

1087 except gf.meta.InvalidTimingSpecification: 

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

1089 

1090 try: 

1091 gdef = gf.meta.parse_grid_spec(sdef) 

1092 except gf.meta.GridSpecError as e: 

1093 die(e) 

1094 

1095 store_dir = get_store_dir(args) 

1096 

1097 try: 

1098 store = gf.Store(store_dir) 

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

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

1101 for phase in phases: 

1102 t = store.t(phase, args) 

1103 if t is not None: 

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

1105 else: 

1106 s.append('nan') 

1107 

1108 if options.output_fn: 

1109 d = dict( 

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

1111 extension='txt') 

1112 

1113 fn = options.output_fn % d 

1114 util.ensuredirs(fn) 

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

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

1117 f.write('\n') 

1118 else: 

1119 print(' '.join(s)) 

1120 

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

1122 die(e) 

1123 

1124 

1125def command_tttlsd(args): 

1126 

1127 def setup(parser): 

1128 pass 

1129 

1130 parser, options, args = cl_parse( 

1131 'tttlsd', args, 

1132 setup=setup, 

1133 details=''' 

1134 

1135This subcommand fills holes in travel-time-tables by using an eikonal solver to 

1136predict travel-times for the missing values. The approach works best for simple 

1137P or S phases without any reflections or conversions. It creates new 

1138travel-time-tables which are named by adding the suffix `.lsd` to the original 

1139name. E.g. running `fomosto tttlsd begin` would produce a new travel-time-table 

1140`begin.lsd`. The new name can be referenced anywhere where a stored 

1141travel-time-table can be used, e.g. in `extra/qseis`. 

1142''') 

1143 

1144 try: 

1145 sphase_ids = args.pop() 

1146 except Exception: 

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

1148 

1149 try: 

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

1151 except gf.meta.InvalidTimingSpecification: 

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

1153 

1154 store_dir = get_store_dir(args) 

1155 

1156 try: 

1157 store = gf.Store(store_dir) 

1158 for phase_id in phase_ids: 

1159 store.fix_ttt_holes(phase_id) 

1160 

1161 except gf.StoreError as e: 

1162 die(e) 

1163 

1164 

1165def command_server(args): 

1166 from pyrocko.gf import server 

1167 

1168 def setup(parser): 

1169 parser.add_option( 

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

1171 help='serve on port PORT') 

1172 

1173 parser.add_option( 

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

1175 help='serve on ip address IP') 

1176 

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

1178 

1179 engine = gf.LocalEngine(store_superdirs=args) 

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

1181 

1182 

1183def command_download(args): 

1184 from pyrocko.gf import ws 

1185 

1186 details = ''' 

1187 

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

1189 

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

1191''' 

1192 

1193 def setup(parser): 

1194 parser.add_option( 

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

1196 help='overwrite existing files') 

1197 

1198 parser, options, args = cl_parse( 

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

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

1201 sys.exit(parser.format_help()) 

1202 

1203 if len(args) == 2: 

1204 site, store_id = args 

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

1206 die('invalid store ID') 

1207 else: 

1208 site, store_id = args[0], None 

1209 

1210 if site not in gf.ws.g_site_abbr: 

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

1212 site = 'http://' + site 

1213 

1214 try: 

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

1216 except ws.DownloadError as e: 

1217 die(str(e)) 

1218 

1219 

1220def command_modelview(args): 

1221 

1222 import matplotlib.pyplot as plt 

1223 import numpy as num 

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

1225 mpl_init() 

1226 

1227 neat_labels = { 

1228 'vp': '$v_p$', 

1229 'vs': '$v_s$', 

1230 'qp': '$Q_p$', 

1231 'qs': '$Q_s$', 

1232 'rho': '$\\rho$'} 

1233 

1234 def setup(parser): 

1235 parser.add_option( 

1236 '--parameters', dest='parameters', 

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

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

1239 

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

1241 

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

1243 

1244 store_dirs = get_store_dirs(args) 

1245 

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

1247 

1248 fig, axes = plt.subplots(1, 

1249 len(parameters), 

1250 sharey=True, 

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

1252 

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

1254 axes = [axes] 

1255 

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

1257 

1258 for store_id in store_dirs: 

1259 try: 

1260 store = gf.Store(store_id) 

1261 mod = store.config.earthmodel_1d 

1262 z = mod.profile('z') 

1263 

1264 for p in parameters: 

1265 ax = axes[p] 

1266 labelspace(ax) 

1267 if '/' in p: 

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

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

1270 else: 

1271 profile = mod.profile(p) 

1272 

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

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

1275 xscaled(1./1000., ax) 

1276 

1277 yscaled(1./km, ax) 

1278 

1279 except gf.StoreError as e: 

1280 die(e) 

1281 

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

1283 ax.grid() 

1284 if p in neat_labels: 

1285 lab = neat_labels[p] 

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

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

1288 else: 

1289 lab = p 

1290 

1291 ax.set_title(lab, y=1.02) 

1292 

1293 if p in units: 

1294 lab += ' ' + units[p] 

1295 

1296 ax.autoscale() 

1297 ax.set_xlabel(lab) 

1298 

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

1300 

1301 handles, labels = ax.get_legend_handles_labels() 

1302 

1303 if len(store_dirs) > 1: 

1304 ax.legend( 

1305 handles, 

1306 labels, 

1307 bbox_to_anchor=(0.5, 0.12), 

1308 bbox_transform=fig.transFigure, 

1309 ncol=3, 

1310 loc='upper center', 

1311 fancybox=True) 

1312 

1313 plt.subplots_adjust(bottom=0.22, 

1314 wspace=0.05) 

1315 

1316 mpl_show(plt) 

1317 

1318 

1319def command_upgrade(args): 

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

1321 store_dirs = get_store_dirs(args) 

1322 try: 

1323 for store_dir in store_dirs: 

1324 store = gf.Store(store_dir) 

1325 nup = store.upgrade() 

1326 if nup == 0: 

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

1328 else: 

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

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

1331 

1332 except gf.StoreError as e: 

1333 die(e) 

1334 

1335 

1336def command_addref(args): 

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

1338 

1339 store_dirs = [] 

1340 filenames = [] 

1341 for arg in args: 

1342 if op.isdir(arg): 

1343 store_dirs.append(arg) 

1344 elif op.isfile(arg): 

1345 filenames.append(arg) 

1346 else: 

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

1348 

1349 if not store_dirs: 

1350 store_dirs.append('.') 

1351 

1352 references = [] 

1353 try: 

1354 for filename in args: 

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

1356 except ImportError: 

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

1358 

1359 if not references: 

1360 die('no references found') 

1361 

1362 for store_dir in store_dirs: 

1363 try: 

1364 store = gf.Store(store_dir) 

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

1366 for ref in references: 

1367 if ref.id in ids: 

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

1369 

1370 ids.append(ref.id) 

1371 store.config.references.append(ref) 

1372 

1373 store.save_config(make_backup=True) 

1374 

1375 except gf.StoreError as e: 

1376 die(e) 

1377 

1378 

1379def command_qc(args): 

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

1381 

1382 store_dir = get_store_dir(args) 

1383 

1384 try: 

1385 store = gf.Store(store_dir) 

1386 s = store.stats() 

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

1388 print('has empty records') 

1389 

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

1391 

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

1393 print('%s empty' % aname) 

1394 

1395 except gf.StoreError as e: 

1396 die(e) 

1397 

1398 

1399def command_report(args): 

1400 from pyrocko.fomosto.report import report_call 

1401 report_call.run_program(args) 

1402 

1403 

1404def main(args=None): 

1405 ''' 

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

1407 ''' 

1408 if args is None: 

1409 args = sys.argv[1:] 

1410 

1411 if len(args) < 1: 

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

1413 

1414 command = args.pop(0) 

1415 

1416 if command in subcommands: 

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

1418 

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

1420 if command == 'help' and args: 

1421 acommand = args[0] 

1422 if acommand in subcommands: 

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

1424 

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

1426 

1427 else: 

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

1429 

1430 

1431if __name__ == '__main__': 

1432 main()