1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6Utility functions for Pyrocko. 

7 

8.. _time-handling-mode: 

9 

10High precision time handling mode 

11................................. 

12 

13Pyrocko can treat timestamps either as standard double precision (64 bit) 

14floating point values, or as high precision float (:py:class:`numpy.float128` 

15or :py:class:`numpy.float96`, whichever is available), aliased here as 

16:py:class:`pyrocko.util.hpfloat`. High precision time stamps are required when 

17handling data with sub-millisecond precision, i.e. kHz/MHz data streams and 

18event catalogs derived from such data. 

19 

20Not all functions in Pyrocko and in programs depending on Pyrocko may work 

21correctly with high precision times. Therefore, Pyrocko's high precision time 

22handling mode has to be actively activated by user config, command line option 

23or enforced within a certain script/program. 

24 

25The default high precision time handling mode can be configured globally with 

26the user config variable 

27:py:gattr:`~pyrocko.config.Config.use_high_precision_time`. Calling the 

28function :py:func:`use_high_precision_time` overrides the default from the 

29config file. This function may be called at startup of a program/script 

30requiring a specific time handling mode. 

31 

32To create a valid time stamp for use in Pyrocko (e.g. in 

33:py:class:`~pyrocko.model.Event` or :py:class:`~pyrocko.trace.Trace` objects), 

34use: 

35 

36.. code-block :: python 

37 

38 import time 

39 from pyrocko import util 

40 

41 # By default using mode selected in user config, override with: 

42 # util.use_high_precision_time(True) # force high precision mode 

43 # util.use_high_precision_time(False) # force low precision mode 

44 

45 t1 = util.str_to_time('2020-08-27 10:22:00') 

46 t2 = util.str_to_time('2020-08-27 10:22:00.111222') 

47 t3 = util.to_time_float(time.time()) 

48 

49 # To get the appropriate float class, use: 

50 

51 time_float = util.get_time_float() 

52 # -> float, numpy.float128 or numpy.float96 

53 [isinstance(t, time_float) for t in [t1, t2, t3]] 

54 # -> [True, True, True] 

55 

56 # Shortcut: 

57 util.check_time_class(t1) 

58 

59Module content 

60.............. 

61''' 

62 

63import time 

64import logging 

65import os 

66import sys 

67import re 

68import calendar 

69import math 

70import fnmatch 

71try: 

72 import fcntl 

73except ImportError: 

74 fcntl = None 

75import optparse 

76import os.path as op 

77import errno 

78 

79import numpy as num 

80from scipy import signal 

81import pyrocko 

82from pyrocko import dummy_progressbar 

83 

84 

85from urllib.parse import urlencode, quote, unquote # noqa 

86from urllib.request import Request, build_opener, HTTPDigestAuthHandler # noqa 

87from urllib.request import urlopen as _urlopen # noqa 

88from urllib.error import HTTPError, URLError # noqa 

89 

90 

91try: 

92 import certifi 

93 import ssl 

94 g_ssl_context = ssl.create_default_context(cafile=certifi.where()) 

95except ImportError: 

96 g_ssl_context = None 

97 

98 

99class URLErrorSSL(URLError): 

100 

101 def __init__(self, urlerror): 

102 self.__dict__ = urlerror.__dict__.copy() 

103 

104 def __str__(self): 

105 return ( 

106 'Requesting web resource failed and the problem could be ' 

107 'related to SSL. Python standard libraries on some older ' 

108 'systems (like Ubuntu 14.04) are known to have trouble ' 

109 'with some SSL setups of today\'s servers: %s' 

110 % URLError.__str__(self)) 

111 

112 

113def urlopen_ssl_check(*args, **kwargs): 

114 try: 

115 return urlopen(*args, **kwargs) 

116 except URLError as e: 

117 if str(e).find('SSL') != -1: 

118 raise URLErrorSSL(e) 

119 else: 

120 raise 

121 

122 

123def urlopen(*args, **kwargs): 

124 

125 if 'context' not in kwargs and g_ssl_context is not None: 

126 kwargs['context'] = g_ssl_context 

127 

128 return _urlopen(*args, **kwargs) 

129 

130 

131try: 

132 long 

133except NameError: 

134 long = int 

135 

136 

137force_dummy_progressbar = False 

138 

139 

140try: 

141 from pyrocko import util_ext 

142except ImportError: 

143 util_ext = None 

144 

145 

146logger = logging.getLogger('pyrocko.util') 

147 

148 

149try: 

150 num_full = num.full 

151except AttributeError: 

152 def num_full(shape, fill_value, dtype=None, order='C'): 

153 a = num.empty(shape, dtype=dtype, order=order) 

154 a.fill(fill_value) 

155 return a 

156 

157try: 

158 num_full_like = num.full_like 

159except AttributeError: 

160 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True): 

161 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok) 

162 a.fill(fill_value) 

163 return a 

164 

165 

166g_setup_logging_args = 'pyrocko', 'warning' 

167 

168 

169def setup_logging(programname='pyrocko', levelname='warning'): 

170 ''' 

171 Initialize logging. 

172 

173 :param programname: program name to be written in log 

174 :param levelname: string indicating the logging level ('debug', 'info', 

175 'warning', 'error', 'critical') 

176 

177 This function is called at startup by most pyrocko programs to set up a 

178 consistent logging format. This is simply a shortcut to a call to 

179 :py:func:`logging.basicConfig()`. 

180 ''' 

181 

182 global g_setup_logging_args 

183 g_setup_logging_args = (programname, levelname) 

184 

185 levels = {'debug': logging.DEBUG, 

186 'info': logging.INFO, 

187 'warning': logging.WARNING, 

188 'error': logging.ERROR, 

189 'critical': logging.CRITICAL} 

190 

191 logging.basicConfig( 

192 level=levels[levelname], 

193 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s') 

194 

195 

196def subprocess_setup_logging_args(): 

197 ''' 

198 Get arguments from previous call to setup_logging. 

199 

200 These can be sent down to a worker process so it can setup its logging 

201 in the same way as the main process. 

202 ''' 

203 return g_setup_logging_args 

204 

205 

206def data_file(fn): 

207 return os.path.join(os.path.split(__file__)[0], 'data', fn) 

208 

209 

210class DownloadError(Exception): 

211 pass 

212 

213 

214class PathExists(DownloadError): 

215 pass 

216 

217 

218def _download(url, fpath, username=None, password=None, 

219 force=False, method='download', stats=None, 

220 status_callback=None, entries_wanted=None, 

221 recursive=False, header=None): 

222 

223 import requests 

224 from requests.auth import HTTPBasicAuth 

225 from requests.exceptions import HTTPError as req_HTTPError 

226 

227 requests.adapters.DEFAULT_RETRIES = 5 

228 urljoin = requests.compat.urljoin 

229 

230 session = requests.Session() 

231 session.header = header 

232 session.auth = None if username is None\ 

233 else HTTPBasicAuth(username, password) 

234 

235 status = { 

236 'ntotal_files': 0, 

237 'nread_files': 0, 

238 'ntotal_bytes_all_files': 0, 

239 'nread_bytes_all_files': 0, 

240 'ntotal_bytes_current_file': 0, 

241 'nread_bytes_current_file': 0, 

242 'finished': False 

243 } 

244 

245 try: 

246 url_to_size = {} 

247 

248 if callable(status_callback): 

249 status_callback(status) 

250 

251 if not recursive and url.endswith('/'): 

252 raise DownloadError( 

253 'URL: %s appears to be a directory' 

254 ' but recurvise download is False' % url) 

255 

256 if recursive and not url.endswith('/'): 

257 url += '/' 

258 

259 def parse_directory_tree(url, subdir=''): 

260 r = session.get(urljoin(url, subdir)) 

261 r.raise_for_status() 

262 

263 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text) 

264 

265 files = sorted(set(subdir + fn for fn in entries 

266 if not fn.endswith('/'))) 

267 

268 if entries_wanted is not None: 

269 files = [fn for fn in files 

270 if (fn in wanted for wanted in entries_wanted)] 

271 

272 status['ntotal_files'] += len(files) 

273 

274 dirs = sorted(set(subdir + dn for dn in entries 

275 if dn.endswith('/') 

276 and dn not in ('./', '../'))) 

277 

278 for dn in dirs: 

279 files.extend(parse_directory_tree( 

280 url, subdir=dn)) 

281 

282 return files 

283 

284 def get_content_length(url): 

285 if url not in url_to_size: 

286 r = session.head(url, headers={'Accept-Encoding': ''}) 

287 

288 content_length = r.headers.get('content-length', None) 

289 if content_length is None: 

290 logger.debug('Could not get HTTP header ' 

291 'Content-Length for %s' % url) 

292 

293 content_length = None 

294 

295 else: 

296 content_length = int(content_length) 

297 status['ntotal_bytes_all_files'] += content_length 

298 if callable(status_callback): 

299 status_callback(status) 

300 

301 url_to_size[url] = content_length 

302 

303 return url_to_size[url] 

304 

305 def download_file(url, fn): 

306 logger.info('starting download of %s...' % url) 

307 ensuredirs(fn) 

308 

309 fsize = get_content_length(url) 

310 status['ntotal_bytes_current_file'] = fsize 

311 status['nread_bytes_current_file'] = 0 

312 if callable(status_callback): 

313 status_callback(status) 

314 

315 r = session.get(url, stream=True, timeout=5) 

316 r.raise_for_status() 

317 

318 frx = 0 

319 fn_tmp = fn + '.%i.temp' % os.getpid() 

320 with open(fn_tmp, 'wb') as f: 

321 for d in r.iter_content(chunk_size=1024): 

322 f.write(d) 

323 frx += len(d) 

324 

325 status['nread_bytes_all_files'] += len(d) 

326 status['nread_bytes_current_file'] += len(d) 

327 if callable(status_callback): 

328 status_callback(status) 

329 

330 os.rename(fn_tmp, fn) 

331 

332 if fsize is not None and frx != fsize: 

333 logger.warning( 

334 'HTTP header Content-Length: %i bytes does not match ' 

335 'download size %i bytes' % (fsize, frx)) 

336 

337 logger.info('finished download of %s' % url) 

338 

339 status['nread_files'] += 1 

340 if callable(status_callback): 

341 status_callback(status) 

342 

343 if recursive: 

344 if op.exists(fpath) and not force: 

345 raise PathExists('path %s already exists' % fpath) 

346 

347 files = parse_directory_tree(url) 

348 

349 dsize = 0 

350 for fn in files: 

351 file_url = urljoin(url, fn) 

352 dsize += get_content_length(file_url) or 0 

353 

354 if method == 'calcsize': 

355 return dsize 

356 

357 else: 

358 for fn in files: 

359 file_url = urljoin(url, fn) 

360 download_file(file_url, op.join(fpath, fn)) 

361 

362 else: 

363 status['ntotal_files'] += 1 

364 if callable(status_callback): 

365 status_callback(status) 

366 

367 fsize = get_content_length(url) 

368 if method == 'calcsize': 

369 return fsize 

370 else: 

371 download_file(url, fpath) 

372 

373 except req_HTTPError as e: 

374 logging.warning("http error: %s" % e) 

375 raise DownloadError('could not download file(s) from: %s' % url) 

376 

377 finally: 

378 status['finished'] = True 

379 if callable(status_callback): 

380 status_callback(status) 

381 session.close() 

382 

383 

384def download_file( 

385 url, fpath, username=None, password=None, status_callback=None, 

386 **kwargs): 

387 return _download( 

388 url, fpath, username, password, 

389 recursive=False, 

390 status_callback=status_callback, 

391 **kwargs) 

392 

393 

394def download_dir( 

395 url, fpath, username=None, password=None, status_callback=None, 

396 **kwargs): 

397 

398 return _download( 

399 url, fpath, username, password, 

400 recursive=True, 

401 status_callback=status_callback, 

402 **kwargs) 

403 

404 

405class HPFloatUnavailable(Exception): 

406 pass 

407 

408 

409class dummy_hpfloat(object): 

410 def __init__(self, *args, **kwargs): 

411 raise HPFloatUnavailable( 

412 'NumPy lacks support for float128 or float96 data type on this ' 

413 'platform.') 

414 

415 

416if hasattr(num, 'float128'): 

417 hpfloat = num.float128 

418 

419elif hasattr(num, 'float96'): 

420 hpfloat = num.float96 

421 

422else: 

423 hpfloat = dummy_hpfloat 

424 

425 

426g_time_float = None 

427g_time_dtype = None 

428 

429 

430class TimeFloatSettingError(Exception): 

431 pass 

432 

433 

434def use_high_precision_time(enabled): 

435 ''' 

436 Globally force a specific time handling mode. 

437 

438 See :ref:`High precision time handling mode <time-handling-mode>`. 

439 

440 :param enabled: enable/disable use of high precision time type 

441 :type enabled: bool 

442 

443 This function should be called before handling/reading any time data. 

444 It can only be called once. 

445 

446 Special attention is required when using multiprocessing on a platform 

447 which does not use fork under the hood. In such cases, the desired setting 

448 must be set also in the subprocess. 

449 ''' 

450 _setup_high_precision_time_mode(enabled_app=enabled) 

451 

452 

453def _setup_high_precision_time_mode(enabled_app=False): 

454 global g_time_float 

455 global g_time_dtype 

456 

457 if not (g_time_float is None and g_time_dtype is None): 

458 raise TimeFloatSettingError( 

459 'Cannot set time handling mode: too late, it has already been ' 

460 'fixed by an earlier call.') 

461 

462 from pyrocko import config 

463 

464 conf = config.config() 

465 enabled_config = conf.use_high_precision_time 

466 

467 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None) 

468 if enabled_env is not None: 

469 try: 

470 enabled_env = int(enabled_env) == 1 

471 except ValueError: 

472 raise TimeFloatSettingError( 

473 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME ' 

474 'should be set to 0 or 1.') 

475 

476 enabled = enabled_config 

477 mode_from = 'config variable `use_high_precision_time`' 

478 notify = enabled 

479 

480 if enabled_env is not None: 

481 if enabled_env != enabled: 

482 notify = True 

483 enabled = enabled_env 

484 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`' 

485 

486 if enabled_app is not None: 

487 if enabled_app != enabled: 

488 notify = True 

489 enabled = enabled_app 

490 mode_from = 'application override' 

491 

492 logger.debug(''' 

493Pyrocko high precision time mode selection (latter override earlier): 

494 config: %s 

495 env: %s 

496 app: %s 

497 -> enabled: %s'''.lstrip() % ( 

498 enabled_config, enabled_env, enabled_app, enabled)) 

499 

500 if notify: 

501 logger.info('Pyrocko high precision time mode %s by %s.' % ( 

502 'activated' if enabled else 'deactivated', 

503 mode_from)) 

504 

505 if enabled: 

506 g_time_float = hpfloat 

507 g_time_dtype = hpfloat 

508 else: 

509 g_time_float = float 

510 g_time_dtype = num.float64 

511 

512 

513def get_time_float(): 

514 ''' 

515 Get the effective float class for timestamps. 

516 

517 See :ref:`High precision time handling mode <time-handling-mode>`. 

518 

519 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the 

520 current time handling mode 

521 ''' 

522 global g_time_float 

523 

524 if g_time_float is None: 

525 _setup_high_precision_time_mode() 

526 

527 return g_time_float 

528 

529 

530def get_time_dtype(): 

531 ''' 

532 Get effective NumPy float class to handle timestamps. 

533 

534 See :ref:`High precision time handling mode <time-handling-mode>`. 

535 ''' 

536 

537 global g_time_dtype 

538 

539 if g_time_dtype is None: 

540 _setup_high_precision_time_mode() 

541 

542 return g_time_dtype 

543 

544 

545def to_time_float(t): 

546 ''' 

547 Convert float to valid time stamp in the current time handling mode. 

548 

549 See :ref:`High precision time handling mode <time-handling-mode>`. 

550 ''' 

551 return get_time_float()(t) 

552 

553 

554class TimestampTypeError(ValueError): 

555 pass 

556 

557 

558def check_time_class(t, error='raise'): 

559 ''' 

560 Type-check variable against current time handling mode. 

561 

562 See :ref:`High precision time handling mode <time-handling-mode>`. 

563 ''' 

564 

565 if t == 0.0: 

566 return 

567 

568 if not isinstance(t, get_time_float()): 

569 message = ( 

570 'Timestamp %g is of type %s but should be of type %s with ' 

571 'Pyrocko\'s currently selected time handling mode.\n\n' 

572 'See https://pyrocko.org/docs/current/library/reference/util.html' 

573 '#high-precision-time-handling-mode' % ( 

574 t, type(t), get_time_float())) 

575 

576 if error == 'raise': 

577 raise TimestampTypeError(message) 

578 elif error == 'warn': 

579 logger.warning(message) 

580 else: 

581 assert False 

582 

583 

584class Stopwatch(object): 

585 ''' 

586 Simple stopwatch to measure elapsed wall clock time. 

587 

588 Usage:: 

589 

590 s = Stopwatch() 

591 time.sleep(1) 

592 print s() 

593 time.sleep(1) 

594 print s() 

595 ''' 

596 

597 def __init__(self): 

598 self.start = time.time() 

599 

600 def __call__(self): 

601 return time.time() - self.start 

602 

603 

604def wrap(text, line_length=80): 

605 ''' 

606 Paragraph and list-aware wrapping of text. 

607 ''' 

608 

609 text = text.strip('\n') 

610 at_lineend = re.compile(r' *\n') 

611 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)') 

612 

613 paragraphs = at_para.split(text)[::5] 

614 listindents = at_para.split(text)[4::5] 

615 newlist = at_para.split(text)[3::5] 

616 

617 listindents[0:0] = [None] 

618 listindents.append(True) 

619 newlist.append(None) 

620 

621 det_indent = re.compile(r'^ *') 

622 

623 outlines = [] 

624 for ip, p in enumerate(paragraphs): 

625 if not p: 

626 continue 

627 

628 if listindents[ip] is None: 

629 _indent = det_indent.findall(p)[0] 

630 findent = _indent 

631 else: 

632 findent = listindents[ip] 

633 _indent = ' ' * len(findent) 

634 

635 ll = line_length - len(_indent) 

636 llf = ll 

637 

638 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())] 

639 p1 = ' '.join(oldlines) 

640 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll)) 

641 for imatch, match in enumerate(possible.finditer(p1)): 

642 parout = match.group(1) 

643 if imatch == 0: 

644 outlines.append(findent + parout) 

645 else: 

646 outlines.append(_indent + parout) 

647 

648 if ip != len(paragraphs)-1 and ( 

649 listindents[ip] is None 

650 or newlist[ip] is not None 

651 or listindents[ip+1] is None): 

652 

653 outlines.append('') 

654 

655 return outlines 

656 

657 

658def ewrap(lines, width=80, indent=''): 

659 lines = list(lines) 

660 if not lines: 

661 return '' 

662 fwidth = max(len(s) for s in lines) 

663 nx = max(1, (80-len(indent)) // (fwidth+1)) 

664 i = 0 

665 rows = [] 

666 while i < len(lines): 

667 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx])) 

668 i += nx 

669 

670 return '\n'.join(rows) 

671 

672 

673class BetterHelpFormatter(optparse.IndentedHelpFormatter): 

674 

675 def __init__(self, *args, **kwargs): 

676 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs) 

677 

678 def format_option(self, option): 

679 ''' 

680 From IndentedHelpFormatter but using a different wrap method. 

681 ''' 

682 

683 help_text_position = 4 + self.current_indent 

684 help_text_width = self.width - help_text_position 

685 

686 opts = self.option_strings[option] 

687 opts = "%*s%s" % (self.current_indent, "", opts) 

688 if option.help: 

689 help_text = self.expand_default(option) 

690 

691 if self.help_position + len(help_text) + 1 <= self.width: 

692 lines = [ 

693 '', 

694 '%-*s %s' % (self.help_position, opts, help_text), 

695 ''] 

696 else: 

697 lines = [''] 

698 lines.append(opts) 

699 lines.append('') 

700 if option.help: 

701 help_lines = wrap(help_text, help_text_width) 

702 lines.extend(["%*s%s" % (help_text_position, "", line) 

703 for line in help_lines]) 

704 lines.append('') 

705 

706 return "\n".join(lines) 

707 

708 def format_description(self, description): 

709 if not description: 

710 return '' 

711 

712 if self.current_indent == 0: 

713 lines = [] 

714 else: 

715 lines = [''] 

716 

717 lines.extend(wrap(description, self.width - self.current_indent)) 

718 if self.current_indent == 0: 

719 lines.append('\n') 

720 

721 return '\n'.join( 

722 ['%*s%s' % (self.current_indent, '', line) for line in lines]) 

723 

724 

725class ProgressBar: 

726 def __init__(self, label, n): 

727 from pyrocko.progress import progress 

728 self._context = progress.view() 

729 self._context.__enter__() 

730 self._task = progress.task(label, n) 

731 

732 def update(self, i): 

733 self._task.update(i) 

734 

735 def finish(self): 

736 self._task.done() 

737 if self._context: 

738 self._context.__exit__() 

739 self._context = None 

740 

741 

742def progressbar(label, maxval): 

743 if force_dummy_progressbar: 

744 return dummy_progressbar.ProgressBar(maxval=maxval).start() 

745 

746 return ProgressBar(label, maxval) 

747 

748 

749def progress_beg(label): 

750 ''' 

751 Notify user that an operation has started. 

752 

753 :param label: name of the operation 

754 

755 To be used in conjuction with :py:func:`progress_end`. 

756 ''' 

757 

758 sys.stderr.write(label) 

759 sys.stderr.flush() 

760 

761 

762def progress_end(label=''): 

763 ''' 

764 Notify user that an operation has ended. 

765 

766 :param label: name of the operation 

767 

768 To be used in conjuction with :py:func:`progress_beg`. 

769 ''' 

770 

771 sys.stderr.write(' done. %s\n' % label) 

772 sys.stderr.flush() 

773 

774 

775class ArangeError(Exception): 

776 pass 

777 

778 

779def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'): 

780 ''' 

781 Return evenly spaced numbers over a specified interval. 

782 

783 Like :py:func:`numpy.arange` but returning floating point numbers by 

784 default and with defined behaviour when stepsize is inconsistent with 

785 interval bounds. It is considered inconsistent if the difference between 

786 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon * 

787 step``. Inconsistencies are handled according to the ``error`` parameter. 

788 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is 

789 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop`` 

790 is silently changed to the closest, the next smaller, or next larger 

791 multiple of ``step``, respectively. 

792 ''' 

793 

794 assert error in ('raise', 'round', 'floor', 'ceil') 

795 

796 start = dtype(start) 

797 stop = dtype(stop) 

798 step = dtype(step) 

799 

800 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round) 

801 

802 n = int(rnd((stop - start) / step)) + 1 

803 stop_check = start + (n-1) * step 

804 

805 if error == 'raise' and abs(stop_check - stop) > step * epsilon: 

806 raise ArangeError( 

807 'inconsistent range specification: start=%g, stop=%g, step=%g' 

808 % (start, stop, step)) 

809 

810 x = num.arange(n, dtype=dtype) 

811 x *= step 

812 x += start 

813 return x 

814 

815 

816def polylinefit(x, y, n_or_xnodes): 

817 ''' 

818 Fit piece-wise linear function to data. 

819 

820 :param x,y: arrays with coordinates of data 

821 :param n_or_xnodes: int, number of segments or x coordinates of polyline 

822 

823 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of 

824 polyline, root-mean-square error 

825 ''' 

826 

827 x = num.asarray(x) 

828 y = num.asarray(y) 

829 

830 if isinstance(n_or_xnodes, int): 

831 n = n_or_xnodes 

832 xmin = x.min() 

833 xmax = x.max() 

834 xnodes = num.linspace(xmin, xmax, n+1) 

835 else: 

836 xnodes = num.asarray(n_or_xnodes) 

837 n = xnodes.size - 1 

838 

839 assert len(x) == len(y) 

840 assert n > 0 

841 

842 ndata = len(x) 

843 a = num.zeros((ndata+(n-1), n*2)) 

844 for i in range(n): 

845 xmin_block = xnodes[i] 

846 xmax_block = xnodes[i+1] 

847 if i == n-1: # don't loose last point 

848 indices = num.where( 

849 num.logical_and(xmin_block <= x, x <= xmax_block))[0] 

850 else: 

851 indices = num.where( 

852 num.logical_and(xmin_block <= x, x < xmax_block))[0] 

853 

854 a[indices, i*2] = x[indices] 

855 a[indices, i*2+1] = 1.0 

856 

857 w = float(ndata)*100. 

858 if i < n-1: 

859 a[ndata+i, i*2] = xmax_block*w 

860 a[ndata+i, i*2+1] = 1.0*w 

861 a[ndata+i, i*2+2] = -xmax_block*w 

862 a[ndata+i, i*2+3] = -1.0*w 

863 

864 d = num.concatenate((y, num.zeros(n-1))) 

865 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2)) 

866 

867 ynodes = num.zeros(n+1) 

868 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1] 

869 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1] 

870 ynodes[1:n] *= 0.5 

871 

872 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2)) 

873 

874 return xnodes, ynodes, rms_error 

875 

876 

877def plf_integrate_piecewise(x_edges, x, y): 

878 ''' 

879 Calculate definite integral of piece-wise linear function on intervals. 

880 

881 Use trapezoidal rule to calculate definite integral of a piece-wise linear 

882 function for a series of consecutive intervals. ``x_edges`` and ``x`` must 

883 be sorted. 

884 

885 :param x_edges: array with edges of the intervals 

886 :param x,y: arrays with coordinates of piece-wise linear function's 

887 control points 

888 ''' 

889 

890 x_all = num.concatenate((x, x_edges)) 

891 ii = num.argsort(x_all) 

892 y_edges = num.interp(x_edges, x, y) 

893 y_all = num.concatenate((y, y_edges)) 

894 xs = x_all[ii] 

895 ys = y_all[ii] 

896 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1])) 

897 return num.diff(y_all[-len(y_edges):]) 

898 

899 

900def diff_fd_1d_4o(dt, data): 

901 ''' 

902 Approximate first derivative of an array (forth order, central FD). 

903 

904 :param dt: sampling interval 

905 :param data: NumPy array with data samples 

906 

907 :returns: NumPy array with same shape as input 

908 

909 Interior points are approximated to fourth order, edge points to first 

910 order right- or left-sided respectively, points next to edge to second 

911 order central. 

912 ''' 

913 import scipy.signal 

914 

915 ddata = num.empty_like(data, dtype=float) 

916 

917 if data.size >= 5: 

918 ddata[2:-2] = scipy.signal.lfilter( 

919 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt) 

920 

921 if data.size >= 3: 

922 ddata[1] = (data[2] - data[0]) / (2. * dt) 

923 ddata[-2] = (data[-1] - data[-3]) / (2. * dt) 

924 

925 if data.size >= 2: 

926 ddata[0] = (data[1] - data[0]) / dt 

927 ddata[-1] = (data[-1] - data[-2]) / dt 

928 

929 if data.size == 1: 

930 ddata[0] = 0.0 

931 

932 return ddata 

933 

934 

935def diff_fd_1d_2o(dt, data): 

936 ''' 

937 Approximate first derivative of an array (second order, central FD). 

938 

939 :param dt: sampling interval 

940 :param data: NumPy array with data samples 

941 

942 :returns: NumPy array with same shape as input 

943 

944 Interior points are approximated to second order, edge points to first 

945 order right- or left-sided respectively. 

946 

947 Uses :py:func:`numpy.gradient`. 

948 ''' 

949 

950 return num.gradient(data, dt) 

951 

952 

953def diff_fd_2d_4o(dt, data): 

954 ''' 

955 Approximate second derivative of an array (forth order, central FD). 

956 

957 :param dt: sampling interval 

958 :param data: NumPy array with data samples 

959 

960 :returns: NumPy array with same shape as input 

961 

962 Interior points are approximated to fourth order, next-to-edge points to 

963 second order, edge points repeated. 

964 ''' 

965 import scipy.signal 

966 

967 ddata = num.empty_like(data, dtype=float) 

968 

969 if data.size >= 5: 

970 ddata[2:-2] = scipy.signal.lfilter( 

971 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2) 

972 

973 if data.size >= 3: 

974 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2 

975 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2 

976 

977 if data.size < 3: 

978 ddata[:] = 0.0 

979 

980 return ddata 

981 

982 

983def diff_fd_2d_2o(dt, data): 

984 ''' 

985 Approximate second derivative of an array (second order, central FD). 

986 

987 :param dt: sampling interval 

988 :param data: NumPy array with data samples 

989 

990 :returns: NumPy array with same shape as input 

991 

992 Interior points are approximated to second order, edge points repeated. 

993 ''' 

994 import scipy.signal 

995 

996 ddata = num.empty_like(data, dtype=float) 

997 

998 if data.size >= 3: 

999 ddata[1:-1] = scipy.signal.lfilter( 

1000 [1., -2., 1.], [1.], data)[2:] / (dt**2) 

1001 

1002 ddata[0] = ddata[1] 

1003 ddata[-1] = ddata[-2] 

1004 

1005 if data.size < 3: 

1006 ddata[:] = 0.0 

1007 

1008 return ddata 

1009 

1010 

1011def diff_fd(n, order, dt, data): 

1012 ''' 

1013 Approximate 1st or 2nd derivative of an array. 

1014 

1015 :param n: 1 for first derivative, 2 for second 

1016 :param order: order of the approximation 2 and 4 are supported 

1017 :param dt: sampling interval 

1018 :param data: NumPy array with data samples 

1019 

1020 :returns: NumPy array with same shape as input 

1021 

1022 This is a frontend to the functions :py:func:`diff_fd_1d_2o`, 

1023 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and 

1024 :py:func:`diff_fd_2d_4o`. 

1025 

1026 Raises :py:exc:`ValueError` for unsupported `n` or `order`. 

1027 ''' 

1028 

1029 funcs = { 

1030 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o}, 

1031 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}} 

1032 

1033 try: 

1034 funcs_n = funcs[n] 

1035 except KeyError: 

1036 raise ValueError( 

1037 'pyrocko.util.diff_fd: ' 

1038 'Only 1st and 2sd derivatives are supported.') 

1039 

1040 try: 

1041 func = funcs_n[order] 

1042 except KeyError: 

1043 raise ValueError( 

1044 'pyrocko.util.diff_fd: ' 

1045 'Order %i is not supported for %s derivative. Supported: %s' % ( 

1046 order, ['', '1st', '2nd'][n], 

1047 ', '.join('%i' % order for order in sorted(funcs_n.keys())))) 

1048 

1049 return func(dt, data) 

1050 

1051 

1052class GlobalVars(object): 

1053 reuse_store = dict() 

1054 decitab_nmax = 0 

1055 decitab = {} 

1056 decimate_fir_coeffs = {} 

1057 decimate_fir_remez_coeffs = {} 

1058 decimate_iir_coeffs = {} 

1059 re_frac = None 

1060 

1061 

1062def decimate_coeffs(q, n=None, ftype='iir'): 

1063 

1064 q = int(q) 

1065 

1066 if n is None: 

1067 if ftype == 'fir': 

1068 n = 30 

1069 elif ftype == 'fir-remez': 

1070 n = 45*q 

1071 else: 

1072 n = 8 

1073 

1074 if ftype == 'fir': 

1075 coeffs = GlobalVars.decimate_fir_coeffs 

1076 if (n, 1./q) not in coeffs: 

1077 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming') 

1078 

1079 b = coeffs[n, 1./q] 

1080 return b, [1.], n 

1081 

1082 elif ftype == 'fir-remez': 

1083 coeffs = GlobalVars.decimate_fir_remez_coeffs 

1084 if (n, 1./q) not in coeffs: 

1085 coeffs[n, 1./q] = signal.remez( 

1086 n+1, (0., .75/q, 1./q, 1.), 

1087 (1., 0.), fs=2, weight=(1, 50)) 

1088 b = coeffs[n, 1./q] 

1089 return b, [1.], n 

1090 

1091 else: 

1092 coeffs = GlobalVars.decimate_iir_coeffs 

1093 if (n, 0.05, 0.8/q) not in coeffs: 

1094 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q) 

1095 

1096 b, a = coeffs[n, 0.05, 0.8/q] 

1097 return b, a, n 

1098 

1099 

1100def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0): 

1101 ''' 

1102 Downsample the signal x by an integer factor q, using an order n filter 

1103 

1104 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR 

1105 filter with hamming window if ftype is 'fir'. 

1106 

1107 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`) 

1108 :param q: the downsampling factor 

1109 :param n: order of the filter (1 less than the length of the filter for a 

1110 `fir` filter) 

1111 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez` 

1112 

1113 :returns: the downsampled signal (1D :class:`numpy.ndarray`) 

1114 

1115 ''' 

1116 

1117 b, a, n = decimate_coeffs(q, n, ftype) 

1118 

1119 if zi is None or zi is True: 

1120 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float) 

1121 else: 

1122 zi_ = zi 

1123 

1124 y, zf = signal.lfilter(b, a, x, zi=zi_) 

1125 

1126 if zi is not None: 

1127 return y[n//2+ioff::q].copy(), zf 

1128 else: 

1129 return y[n//2+ioff::q].copy() 

1130 

1131 

1132class UnavailableDecimation(Exception): 

1133 ''' 

1134 Exception raised by :py:func:`decitab` for unavailable decimation factors. 

1135 ''' 

1136 

1137 pass 

1138 

1139 

1140def gcd(a, b, epsilon=1e-7): 

1141 ''' 

1142 Greatest common divisor. 

1143 ''' 

1144 

1145 while b > epsilon*a: 

1146 a, b = b, a % b 

1147 

1148 return a 

1149 

1150 

1151def lcm(a, b): 

1152 ''' 

1153 Least common multiple. 

1154 ''' 

1155 

1156 return a*b // gcd(a, b) 

1157 

1158 

1159def mk_decitab(nmax=100): 

1160 ''' 

1161 Make table with decimation sequences. 

1162 

1163 Decimation from one sampling rate to a lower one is achieved by a 

1164 successive application of :py:func:`decimation` with small integer 

1165 downsampling factors (because using large downampling factors can make the 

1166 decimation unstable or slow). This function sets up a table with downsample 

1167 sequences for factors up to ``nmax``. 

1168 ''' 

1169 

1170 tab = GlobalVars.decitab 

1171 for i in range(1, 10): 

1172 for j in range(1, i+1): 

1173 for k in range(1, j+1): 

1174 for l_ in range(1, k+1): 

1175 for m in range(1, l_+1): 

1176 p = i*j*k*l_*m 

1177 if p > nmax: 

1178 break 

1179 if p not in tab: 

1180 tab[p] = (i, j, k, l_, m) 

1181 if i*j*k*l_ > nmax: 

1182 break 

1183 if i*j*k > nmax: 

1184 break 

1185 if i*j > nmax: 

1186 break 

1187 if i > nmax: 

1188 break 

1189 

1190 GlobalVars.decitab_nmax = nmax 

1191 

1192 

1193def zfmt(n): 

1194 return '%%0%ii' % (int(math.log10(n - 1)) + 1) 

1195 

1196 

1197def _year_to_time(year): 

1198 tt = (year, 1, 1, 0, 0, 0) 

1199 return to_time_float(calendar.timegm(tt)) 

1200 

1201 

1202def _working_year(year): 

1203 try: 

1204 tt = (year, 1, 1, 0, 0, 0) 

1205 t = calendar.timegm(tt) 

1206 tt2_ = time.gmtime(t) 

1207 tt2 = tuple(tt2_)[:6] 

1208 if tt != tt2: 

1209 return False 

1210 

1211 s = '%i-01-01 00:00:00' % year 

1212 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_) 

1213 if s != s2: 

1214 return False 

1215 

1216 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S') 

1217 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S') 

1218 if s3 != s2: 

1219 return False 

1220 

1221 except Exception: 

1222 return False 

1223 

1224 return True 

1225 

1226 

1227def working_system_time_range(year_min_lim=None, year_max_lim=None): 

1228 ''' 

1229 Check time range supported by the systems's time conversion functions. 

1230 

1231 Returns system time stamps of start of year of first/last fully supported 

1232 year span. If this is before 1900 or after 2100, return first/last century 

1233 which is fully supported. 

1234 

1235 :returns: ``(tmin, tmax, year_min, year_max)`` 

1236 ''' 

1237 

1238 year0 = 2000 

1239 year_min = year0 

1240 year_max = year0 

1241 

1242 itests = list(range(101)) 

1243 for i in range(19): 

1244 itests.append(200 + i*100) 

1245 

1246 for i in itests: 

1247 year = year0 - i 

1248 if year_min_lim is not None and year < year_min_lim: 

1249 year_min = year_min_lim 

1250 break 

1251 elif not _working_year(year): 

1252 break 

1253 else: 

1254 year_min = year 

1255 

1256 for i in itests: 

1257 year = year0 + i 

1258 if year_max_lim is not None and year > year_max_lim: 

1259 year_max = year_max_lim 

1260 break 

1261 elif not _working_year(year + 1): 

1262 break 

1263 else: 

1264 year_max = year 

1265 

1266 return ( 

1267 _year_to_time(year_min), 

1268 _year_to_time(year_max), 

1269 year_min, year_max) 

1270 

1271 

1272g_working_system_time_range = None 

1273 

1274 

1275def get_working_system_time_range(): 

1276 ''' 

1277 Caching variant of :py:func:`working_system_time_range`. 

1278 ''' 

1279 

1280 global g_working_system_time_range 

1281 if g_working_system_time_range is None: 

1282 g_working_system_time_range = working_system_time_range() 

1283 

1284 return g_working_system_time_range 

1285 

1286 

1287def is_working_time(t): 

1288 tmin, tmax, _, _ = get_working_system_time_range() 

1289 return tmin <= t <= tmax 

1290 

1291 

1292def julian_day_of_year(timestamp): 

1293 ''' 

1294 Get the day number after the 1st of January of year in ``timestamp``. 

1295 

1296 :returns: day number as int 

1297 ''' 

1298 

1299 return time.gmtime(int(timestamp)).tm_yday 

1300 

1301 

1302def hour_start(timestamp): 

1303 ''' 

1304 Get beginning of hour for any point in time. 

1305 

1306 :param timestamp: time instant as system timestamp (in seconds) 

1307 

1308 :returns: instant of hour start as system timestamp 

1309 ''' 

1310 

1311 tt = time.gmtime(int(timestamp)) 

1312 tts = tt[0:4] + (0, 0) 

1313 return to_time_float(calendar.timegm(tts)) 

1314 

1315 

1316def day_start(timestamp): 

1317 ''' 

1318 Get beginning of day for any point in time. 

1319 

1320 :param timestamp: time instant as system timestamp (in seconds) 

1321 

1322 :returns: instant of day start as system timestamp 

1323 ''' 

1324 

1325 tt = time.gmtime(int(timestamp)) 

1326 tts = tt[0:3] + (0, 0, 0) 

1327 return to_time_float(calendar.timegm(tts)) 

1328 

1329 

1330def month_start(timestamp): 

1331 ''' 

1332 Get beginning of month for any point in time. 

1333 

1334 :param timestamp: time instant as system timestamp (in seconds) 

1335 

1336 :returns: instant of month start as system timestamp 

1337 ''' 

1338 

1339 tt = time.gmtime(int(timestamp)) 

1340 tts = tt[0:2] + (1, 0, 0, 0) 

1341 return to_time_float(calendar.timegm(tts)) 

1342 

1343 

1344def year_start(timestamp): 

1345 ''' 

1346 Get beginning of year for any point in time. 

1347 

1348 :param timestamp: time instant as system timestamp (in seconds) 

1349 

1350 :returns: instant of year start as system timestamp 

1351 ''' 

1352 

1353 tt = time.gmtime(int(timestamp)) 

1354 tts = tt[0:1] + (1, 1, 0, 0, 0) 

1355 return to_time_float(calendar.timegm(tts)) 

1356 

1357 

1358def iter_days(tmin, tmax): 

1359 ''' 

1360 Yields begin and end of days until given time span is covered. 

1361 

1362 :param tmin,tmax: input time span 

1363 

1364 :yields: tuples with (begin, end) of days as system timestamps 

1365 ''' 

1366 

1367 t = day_start(tmin) 

1368 while t < tmax: 

1369 tend = day_start(t + 26*60*60) 

1370 yield t, tend 

1371 t = tend 

1372 

1373 

1374def iter_months(tmin, tmax): 

1375 ''' 

1376 Yields begin and end of months until given time span is covered. 

1377 

1378 :param tmin,tmax: input time span 

1379 

1380 :yields: tuples with (begin, end) of months as system timestamps 

1381 ''' 

1382 

1383 t = month_start(tmin) 

1384 while t < tmax: 

1385 tend = month_start(t + 24*60*60*33) 

1386 yield t, tend 

1387 t = tend 

1388 

1389 

1390def iter_years(tmin, tmax): 

1391 ''' 

1392 Yields begin and end of years until given time span is covered. 

1393 

1394 :param tmin,tmax: input time span 

1395 

1396 :yields: tuples with (begin, end) of years as system timestamps 

1397 ''' 

1398 

1399 t = year_start(tmin) 

1400 while t < tmax: 

1401 tend = year_start(t + 24*60*60*369) 

1402 yield t, tend 

1403 t = tend 

1404 

1405 

1406def today(): 

1407 return day_start(time.time()) 

1408 

1409 

1410def tomorrow(): 

1411 return day_start(time.time() + 24*60*60) 

1412 

1413 

1414def decitab(n): 

1415 ''' 

1416 Get integer decimation sequence for given downampling factor. 

1417 

1418 :param n: target decimation factor 

1419 

1420 :returns: tuple with downsampling sequence 

1421 ''' 

1422 

1423 if n > GlobalVars.decitab_nmax: 

1424 mk_decitab(n*2) 

1425 if n not in GlobalVars.decitab: 

1426 raise UnavailableDecimation('ratio = %g' % n) 

1427 return GlobalVars.decitab[n] 

1428 

1429 

1430def ctimegm(s, format="%Y-%m-%d %H:%M:%S"): 

1431 ''' 

1432 Convert string representing UTC time to system time. 

1433 

1434 :param s: string to be interpreted 

1435 :param format: format string passed to :py:func:`strptime` 

1436 

1437 :returns: system time stamp 

1438 

1439 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime. 

1440 

1441 .. note:: 

1442 This function is to be replaced by :py:func:`str_to_time`. 

1443 ''' 

1444 

1445 return calendar.timegm(time.strptime(s, format)) 

1446 

1447 

1448def gmctime(t, format="%Y-%m-%d %H:%M:%S"): 

1449 ''' 

1450 Get string representation from system time, UTC. 

1451 

1452 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime. 

1453 

1454 .. note:: 

1455 This function is to be repaced by :py:func:`time_to_str`. 

1456 ''' 

1457 

1458 return time.strftime(format, time.gmtime(t)) 

1459 

1460 

1461def gmctime_v(t, format="%a, %d %b %Y %H:%M:%S"): 

1462 ''' 

1463 Get string representation from system time, UTC. Same as 

1464 :py:func:`gmctime` but with a more verbose default format. 

1465 

1466 .. note:: 

1467 This function is to be replaced by :py:func:`time_to_str`. 

1468 ''' 

1469 

1470 return time.strftime(format, time.gmtime(t)) 

1471 

1472 

1473def gmctime_fn(t, format="%Y-%m-%d_%H-%M-%S"): 

1474 ''' 

1475 Get string representation from system time, UTC. Same as 

1476 :py:func:`gmctime` but with a default usable in filenames. 

1477 

1478 .. note:: 

1479 This function is to be replaced by :py:func:`time_to_str`. 

1480 ''' 

1481 

1482 return time.strftime(format, time.gmtime(t)) 

1483 

1484 

1485class TimeStrError(Exception): 

1486 pass 

1487 

1488 

1489class FractionalSecondsMissing(TimeStrError): 

1490 ''' 

1491 Exception raised by :py:func:`str_to_time` when the given string lacks 

1492 fractional seconds. 

1493 ''' 

1494 

1495 pass 

1496 

1497 

1498class FractionalSecondsWrongNumberOfDigits(TimeStrError): 

1499 ''' 

1500 Exception raised by :py:func:`str_to_time` when the given string has an 

1501 incorrect number of digits in the fractional seconds part. 

1502 ''' 

1503 

1504 pass 

1505 

1506 

1507def _endswith_n(s, endings): 

1508 for ix, x in enumerate(endings): 

1509 if s.endswith(x): 

1510 return ix 

1511 return -1 

1512 

1513 

1514def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'): 

1515 ''' 

1516 Convert string representing UTC time to floating point system time. 

1517 

1518 :param s: string representing UTC time 

1519 :param format: time string format 

1520 :returns: system time stamp as floating point value 

1521 

1522 Uses the semantics of :py:func:`time.strptime` but allows for fractional 

1523 seconds. If the format ends with ``'.FRAC'``, anything after a dot is 

1524 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``, 

1525 the fractional part, including the dot is made optional. The latter has the 

1526 consequence, that the time strings and the format may not contain any other 

1527 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is 

1528 ensured, that exactly that number of digits are present in the fractional 

1529 seconds. 

1530 ''' 

1531 

1532 time_float = get_time_float() 

1533 

1534 if util_ext is not None: 

1535 try: 

1536 t, tfrac = util_ext.stt(s, format) 

1537 except util_ext.UtilExtError as e: 

1538 raise TimeStrError( 

1539 '%s, string=%s, format=%s' % (str(e), s, format)) 

1540 

1541 return time_float(t)+tfrac 

1542 

1543 fracsec = 0. 

1544 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC' 

1545 

1546 iend = _endswith_n(format, fixed_endings) 

1547 if iend != -1: 

1548 dotpos = s.rfind('.') 

1549 if dotpos == -1: 

1550 raise FractionalSecondsMissing( 

1551 'string=%s, format=%s' % (s, format)) 

1552 

1553 if iend > 0 and iend != (len(s)-dotpos-1): 

1554 raise FractionalSecondsWrongNumberOfDigits( 

1555 'string=%s, format=%s' % (s, format)) 

1556 

1557 format = format[:-len(fixed_endings[iend])] 

1558 fracsec = float(s[dotpos:]) 

1559 s = s[:dotpos] 

1560 

1561 elif format.endswith('.OPTFRAC'): 

1562 dotpos = s.rfind('.') 

1563 format = format[:-8] 

1564 if dotpos != -1 and len(s[dotpos:]) > 1: 

1565 fracsec = float(s[dotpos:]) 

1566 

1567 if dotpos != -1: 

1568 s = s[:dotpos] 

1569 

1570 try: 

1571 return time_float(calendar.timegm(time.strptime(s, format))) \ 

1572 + fracsec 

1573 except ValueError as e: 

1574 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format)) 

1575 

1576 

1577stt = str_to_time 

1578 

1579 

1580def str_to_time_fillup(s): 

1581 ''' 

1582 Default :py:func:`str_to_time` with filling in of missing values. 

1583 

1584 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`, 

1585 `'2010-01-01 00'`, ..., or `'2010'`. 

1586 ''' 

1587 

1588 if len(s) in (4, 7, 10, 13, 16): 

1589 s += '0000-01-01 00:00:00'[len(s):] 

1590 

1591 return str_to_time(s) 

1592 

1593 

1594def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'): 

1595 ''' 

1596 Get string representation for floating point system time. 

1597 

1598 :param t: floating point system time 

1599 :param format: time string format 

1600 :returns: string representing UTC time 

1601 

1602 Uses the semantics of :py:func:`time.strftime` but additionally allows for 

1603 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a 

1604 digit between 1 and 9, this is replaced with the fractional part of ``t`` 

1605 with ``x`` digits precision. 

1606 ''' 

1607 

1608 if pyrocko.grumpy > 0: 

1609 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise') 

1610 

1611 if isinstance(format, int): 

1612 if format > 0: 

1613 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format 

1614 else: 

1615 format = '%Y-%m-%d %H:%M:%S' 

1616 

1617 if util_ext is not None: 

1618 t0 = num.floor(t) 

1619 try: 

1620 return util_ext.tts(int(t0), float(t - t0), format) 

1621 except util_ext.UtilExtError as e: 

1622 raise TimeStrError( 

1623 '%s, timestamp=%f, format=%s' % (str(e), t, format)) 

1624 

1625 if not GlobalVars.re_frac: 

1626 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC') 

1627 GlobalVars.frac_formats = dict( 

1628 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789']) 

1629 

1630 ts = float(num.floor(t)) 

1631 tfrac = t-ts 

1632 

1633 m = GlobalVars.re_frac.search(format) 

1634 if m: 

1635 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac) 

1636 if sfrac[0] == '1': 

1637 ts += 1. 

1638 

1639 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1) 

1640 

1641 return time.strftime(format, time.gmtime(ts)) 

1642 

1643 

1644tts = time_to_str 

1645_abbr_weekday = 'Mon Tue Wed Thu Fri Sat Sun'.split() 

1646_abbr_month = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split() 

1647 

1648 

1649def mystrftime(fmt=None, tt=None, milliseconds=0): 

1650 # Needed by Snuffler for the time axis. In other cases `time_to_str` 

1651 # should be used. 

1652 

1653 if fmt is None: 

1654 fmt = '%Y-%m-%d %H:%M:%S .%r' 

1655 

1656 # Get these two locale independent, needed by Snuffler. 

1657 # Setting the locale seems to have no effect. 

1658 fmt = fmt.replace('%a', _abbr_weekday[tt.tm_wday]) 

1659 fmt = fmt.replace('%b', _abbr_month[tt.tm_mon-1]) 

1660 

1661 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds))) 

1662 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000))) 

1663 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000))) 

1664 return time.strftime(fmt, tt) 

1665 

1666 

1667def gmtime_x(timestamp): 

1668 etimestamp = float(num.floor(timestamp)) 

1669 tt = time.gmtime(etimestamp) 

1670 ms = (timestamp-etimestamp)*1000 

1671 return tt, ms 

1672 

1673 

1674def plural_s(n): 

1675 if not isinstance(n, int): 

1676 n = len(n) 

1677 

1678 return 's' if n != 1 else '' 

1679 

1680 

1681def ensuredirs(dst): 

1682 ''' 

1683 Create all intermediate path components for a target path. 

1684 

1685 :param dst: target path 

1686 

1687 The leaf part of the target path is not created (use :py:func:`ensuredir` 

1688 if a the target path is a directory to be created). 

1689 ''' 

1690 

1691 d, x = os.path.split(dst.rstrip(os.sep)) 

1692 dirs = [] 

1693 while d and not os.path.exists(d): 

1694 dirs.append(d) 

1695 d, x = os.path.split(d) 

1696 

1697 dirs.reverse() 

1698 

1699 for d in dirs: 

1700 try: 

1701 os.mkdir(d) 

1702 except OSError as e: 

1703 if not e.errno == errno.EEXIST: 

1704 raise 

1705 

1706 

1707def ensuredir(dst): 

1708 ''' 

1709 Create directory and all intermediate path components to it as needed. 

1710 

1711 :param dst: directory name 

1712 

1713 Nothing is done if the given target already exists. 

1714 ''' 

1715 

1716 if os.path.exists(dst): 

1717 return 

1718 

1719 dst.rstrip(os.sep) 

1720 

1721 ensuredirs(dst) 

1722 try: 

1723 os.mkdir(dst) 

1724 except OSError as e: 

1725 if not e.errno == errno.EEXIST: 

1726 raise 

1727 

1728 

1729def reuse(x): 

1730 ''' 

1731 Get unique instance of an object. 

1732 

1733 :param x: hashable object 

1734 :returns: reference to x or an equivalent object 

1735 

1736 Cache object ``x`` in a global dict for reuse, or if x already 

1737 is in that dict, return a reference to it. 

1738 ''' 

1739 

1740 grs = GlobalVars.reuse_store 

1741 if x not in grs: 

1742 grs[x] = x 

1743 return grs[x] 

1744 

1745 

1746def deuse(x): 

1747 grs = GlobalVars.reuse_store 

1748 if x in grs: 

1749 del grs[x] 

1750 

1751 

1752class Anon(object): 

1753 ''' 

1754 Dict-to-object utility. 

1755 

1756 Any given arguments are stored as attributes. 

1757 

1758 Example:: 

1759 

1760 a = Anon(x=1, y=2) 

1761 print a.x, a.y 

1762 ''' 

1763 

1764 def __init__(self, **dict): 

1765 for k in dict: 

1766 self.__dict__[k] = dict[k] 

1767 

1768 

1769def iter_select_files( 

1770 paths, 

1771 include=None, 

1772 exclude=None, 

1773 selector=None, 

1774 show_progress=True, 

1775 pass_through=None): 

1776 

1777 ''' 

1778 Recursively select files (generator variant). 

1779 

1780 See :py:func:`select_files`. 

1781 ''' 

1782 

1783 if show_progress: 

1784 progress_beg('selecting files...') 

1785 

1786 ngood = 0 

1787 check_include = None 

1788 if include is not None: 

1789 rinclude = re.compile(include) 

1790 

1791 def check_include(path): 

1792 m = rinclude.search(path) 

1793 if not m: 

1794 return False 

1795 

1796 if selector is None: 

1797 return True 

1798 

1799 infos = Anon(**m.groupdict()) 

1800 return selector(infos) 

1801 

1802 check_exclude = None 

1803 if exclude is not None: 

1804 rexclude = re.compile(exclude) 

1805 

1806 def check_exclude(path): 

1807 return not bool(rexclude.search(path)) 

1808 

1809 if check_include and check_exclude: 

1810 

1811 def check(path): 

1812 return check_include(path) and check_exclude(path) 

1813 

1814 elif check_include: 

1815 check = check_include 

1816 

1817 elif check_exclude: 

1818 check = check_exclude 

1819 

1820 else: 

1821 check = None 

1822 

1823 if isinstance(paths, str): 

1824 paths = [paths] 

1825 

1826 for path in paths: 

1827 if pass_through and pass_through(path): 

1828 if check is None or check(path): 

1829 yield path 

1830 

1831 elif os.path.isdir(path): 

1832 for (dirpath, dirnames, filenames) in os.walk(path): 

1833 dirnames.sort() 

1834 filenames.sort() 

1835 for filename in filenames: 

1836 path = op.join(dirpath, filename) 

1837 if check is None or check(path): 

1838 yield os.path.abspath(path) 

1839 ngood += 1 

1840 else: 

1841 if check is None or check(path): 

1842 yield os.path.abspath(path) 

1843 ngood += 1 

1844 

1845 if show_progress: 

1846 progress_end('%i file%s selected.' % (ngood, plural_s(ngood))) 

1847 

1848 

1849def select_files( 

1850 paths, include=None, exclude=None, selector=None, show_progress=True, 

1851 regex=None): 

1852 

1853 ''' 

1854 Recursively select files. 

1855 

1856 :param paths: entry path names 

1857 :param include: pattern for conditional inclusion 

1858 :param exclude: pattern for conditional exclusion 

1859 :param selector: callback for conditional inclusion 

1860 :param show_progress: if True, indicate start and stop of processing 

1861 :param regex: alias for ``include`` (backwards compatibility) 

1862 :returns: list of path names 

1863 

1864 Recursively finds all files under given entry points ``paths``. If 

1865 parameter ``include`` is a regular expression, only files with matching 

1866 path names are included. If additionally parameter ``selector`` is given a 

1867 callback function, only files for which the callback returns ``True`` are 

1868 included. The callback should take a single argument. The callback is 

1869 called with a single argument, an object, having as attributes, any named 

1870 groups given in ``include``. 

1871 

1872 Examples 

1873 

1874 To find all files ending in ``'.mseed'`` or ``'.msd'``:: 

1875 

1876 select_files(paths, 

1877 include=r'\\.(mseed|msd)$') 

1878 

1879 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for 

1880 the year:: 

1881 

1882 select_files(paths, 

1883 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$', 

1884 selector=(lambda x: int(x.year) == 2009)) 

1885 ''' 

1886 if regex is not None: 

1887 assert include is None 

1888 include = regex 

1889 

1890 return list(iter_select_files( 

1891 paths, include, exclude, selector, show_progress)) 

1892 

1893 

1894def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'): 

1895 ''' 

1896 Convert positive integer to a base36 string. 

1897 ''' 

1898 

1899 if not isinstance(number, (int, long)): 

1900 raise TypeError('number must be an integer') 

1901 if number < 0: 

1902 raise ValueError('number must be positive') 

1903 

1904 # Special case for small numbers 

1905 if number < 36: 

1906 return alphabet[number] 

1907 

1908 base36 = '' 

1909 while number != 0: 

1910 number, i = divmod(number, 36) 

1911 base36 = alphabet[i] + base36 

1912 

1913 return base36 

1914 

1915 

1916def base36decode(number): 

1917 ''' 

1918 Decode base36 endcoded positive integer. 

1919 ''' 

1920 

1921 return int(number, 36) 

1922 

1923 

1924class UnpackError(Exception): 

1925 ''' 

1926 Exception raised when :py:func:`unpack_fixed` encounters an error. 

1927 ''' 

1928 

1929 pass 

1930 

1931 

1932ruler = ''.join(['%-10i' % i for i in range(8)]) \ 

1933 + '\n' + '0123456789' * 8 + '\n' 

1934 

1935 

1936def unpack_fixed(format, line, *callargs): 

1937 ''' 

1938 Unpack fixed format string, as produced by many fortran codes. 

1939 

1940 :param format: format specification 

1941 :param line: string to be processed 

1942 :param callargs: callbacks for callback fields in the format 

1943 

1944 The format is described by a string of comma-separated fields. Each field 

1945 is defined by a character for the field type followed by the field width. A 

1946 questionmark may be appended to the field description to allow the argument 

1947 to be optional (The data string is then allowed to be filled with blanks 

1948 and ``None`` is returned in this case). 

1949 

1950 The following field types are available: 

1951 

1952 ==== ================================================================ 

1953 Type Description 

1954 ==== ================================================================ 

1955 A string (full field width is extracted) 

1956 a string (whitespace at the beginning and the end is removed) 

1957 i integer value 

1958 f floating point value 

1959 @ special type, a callback must be given for the conversion 

1960 x special field type to skip parts of the string 

1961 ==== ================================================================ 

1962 ''' 

1963 

1964 ipos = 0 

1965 values = [] 

1966 icall = 0 

1967 for form in format.split(','): 

1968 form = form.strip() 

1969 optional = form[-1] == '?' 

1970 form = form.rstrip('?') 

1971 typ = form[0] 

1972 ln = int(form[1:]) 

1973 s = line[ipos:ipos+ln] 

1974 cast = { 

1975 'x': None, 

1976 'A': str, 

1977 'a': lambda x: x.strip(), 

1978 'i': int, 

1979 'f': float, 

1980 '@': 'extra'}[typ] 

1981 

1982 if cast == 'extra': 

1983 cast = callargs[icall] 

1984 icall += 1 

1985 

1986 if cast is not None: 

1987 if optional and s.strip() == '': 

1988 values.append(None) 

1989 else: 

1990 try: 

1991 values.append(cast(s)) 

1992 except Exception: 

1993 mark = [' '] * 80 

1994 mark[ipos:ipos+ln] = ['^'] * ln 

1995 mark = ''.join(mark) 

1996 raise UnpackError( 

1997 'Invalid cast to type "%s" at position [%i:%i] of ' 

1998 'line: \n%s%s\n%s' % ( 

1999 typ, ipos, ipos+ln, ruler, line.rstrip(), mark)) 

2000 

2001 ipos += ln 

2002 

2003 return values 

2004 

2005 

2006_pattern_cache = {} 

2007 

2008 

2009def _nslc_pattern(pattern): 

2010 if pattern not in _pattern_cache: 

2011 rpattern = re.compile(fnmatch.translate(pattern), re.I) 

2012 _pattern_cache[pattern] = rpattern 

2013 else: 

2014 rpattern = _pattern_cache[pattern] 

2015 

2016 return rpattern 

2017 

2018 

2019def match_nslc(patterns, nslc): 

2020 ''' 

2021 Match network-station-location-channel code against pattern or list of 

2022 patterns. 

2023 

2024 :param patterns: pattern or list of patterns 

2025 :param nslc: tuple with (network, station, location, channel) as strings 

2026 

2027 :returns: ``True`` if the pattern matches or if any of the given patterns 

2028 match; or ``False``. 

2029 

2030 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq]. 

2031 

2032 Example:: 

2033 

2034 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True 

2035 ''' 

2036 

2037 if isinstance(patterns, str): 

2038 patterns = [patterns] 

2039 

2040 if not isinstance(nslc, str): 

2041 s = '.'.join(nslc) 

2042 else: 

2043 s = nslc 

2044 

2045 for pattern in patterns: 

2046 if _nslc_pattern(pattern).match(s): 

2047 return True 

2048 

2049 return False 

2050 

2051 

2052def match_nslcs(patterns, nslcs): 

2053 ''' 

2054 Get network-station-location-channel codes that match given pattern or any 

2055 of several given patterns. 

2056 

2057 :param patterns: pattern or list of patterns 

2058 :param nslcs: list of (network, station, location, channel) tuples 

2059 

2060 See also :py:func:`match_nslc` 

2061 ''' 

2062 

2063 matching = [] 

2064 for nslc in nslcs: 

2065 if match_nslc(patterns, nslc): 

2066 matching.append(nslc) 

2067 

2068 return matching 

2069 

2070 

2071class Timeout(Exception): 

2072 pass 

2073 

2074 

2075def create_lockfile(fn, timeout=None, timewarn=10.): 

2076 t0 = time.time() 

2077 

2078 while True: 

2079 try: 

2080 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL) 

2081 os.close(f) 

2082 return 

2083 

2084 except OSError as e: 

2085 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows 

2086 pass # retry 

2087 else: 

2088 raise 

2089 

2090 tnow = time.time() 

2091 

2092 if timeout is not None and tnow - t0 > timeout: 

2093 raise Timeout( 

2094 'Timeout (%gs) occured while waiting to get exclusive ' 

2095 'access to: %s' % (timeout, fn)) 

2096 

2097 if timewarn is not None and tnow - t0 > timewarn: 

2098 logger.warning( 

2099 'Waiting since %gs to get exclusive access to: %s' % ( 

2100 timewarn, fn)) 

2101 

2102 timewarn *= 2 

2103 

2104 time.sleep(0.01) 

2105 

2106 

2107def delete_lockfile(fn): 

2108 os.unlink(fn) 

2109 

2110 

2111class Lockfile(Exception): 

2112 

2113 def __init__(self, path, timeout=5, timewarn=10.): 

2114 self._path = path 

2115 self._timeout = timeout 

2116 self._timewarn = timewarn 

2117 

2118 def __enter__(self): 

2119 create_lockfile( 

2120 self._path, timeout=self._timeout, timewarn=self._timewarn) 

2121 return None 

2122 

2123 def __exit__(self, type, value, traceback): 

2124 delete_lockfile(self._path) 

2125 

2126 

2127class SoleError(Exception): 

2128 ''' 

2129 Exception raised by objects of type :py:class:`Sole`, when an concurrent 

2130 instance is running. 

2131 ''' 

2132 

2133 pass 

2134 

2135 

2136class Sole(object): 

2137 ''' 

2138 Use POSIX advisory file locking to ensure that only a single instance of a 

2139 program is running. 

2140 

2141 :param pid_path: path to lockfile to be used 

2142 

2143 Usage:: 

2144 

2145 from pyrocko.util import Sole, SoleError, setup_logging 

2146 import os 

2147 

2148 setup_logging('my_program') 

2149 

2150 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock') 

2151 try: 

2152 sole = Sole(pid_path) 

2153 

2154 except SoleError, e: 

2155 logger.fatal( str(e) ) 

2156 sys.exit(1) 

2157 ''' 

2158 

2159 def __init__(self, pid_path): 

2160 self._pid_path = pid_path 

2161 self._other_running = False 

2162 ensuredirs(self._pid_path) 

2163 self._lockfile = None 

2164 self._os = os 

2165 

2166 if not fcntl: 

2167 raise SoleError( 

2168 'Python standard library module "fcntl" not available on ' 

2169 'this platform.') 

2170 

2171 self._fcntl = fcntl 

2172 

2173 try: 

2174 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY) 

2175 except Exception: 

2176 raise SoleError( 

2177 'Cannot open lockfile (path = %s)' % self._pid_path) 

2178 

2179 try: 

2180 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB) 

2181 

2182 except IOError: 

2183 self._other_running = True 

2184 try: 

2185 f = open(self._pid_path, 'r') 

2186 pid = f.read().strip() 

2187 f.close() 

2188 except Exception: 

2189 pid = '?' 

2190 

2191 raise SoleError('Other instance is running (pid = %s)' % pid) 

2192 

2193 try: 

2194 os.ftruncate(self._lockfile, 0) 

2195 os.write(self._lockfile, '%i\n' % os.getpid()) 

2196 os.fsync(self._lockfile) 

2197 

2198 except Exception: 

2199 # the pid is only stored for user information, so this is allowed 

2200 # to fail 

2201 pass 

2202 

2203 def __del__(self): 

2204 if not self._other_running: 

2205 if self._lockfile is not None: 

2206 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN) 

2207 self._os.close(self._lockfile) 

2208 try: 

2209 self._os.unlink(self._pid_path) 

2210 except Exception: 

2211 pass 

2212 

2213 

2214re_escapequotes = re.compile(r"(['\\])") 

2215 

2216 

2217def escapequotes(s): 

2218 return re_escapequotes.sub(r"\\\1", s) 

2219 

2220 

2221class TableWriter(object): 

2222 ''' 

2223 Write table of space separated values to a file. 

2224 

2225 :param f: file like object 

2226 

2227 Strings containing spaces are quoted on output. 

2228 ''' 

2229 

2230 def __init__(self, f): 

2231 self._f = f 

2232 

2233 def writerow(self, row, minfieldwidths=None): 

2234 

2235 ''' 

2236 Write one row of values to underlying file. 

2237 

2238 :param row: iterable of values 

2239 :param minfieldwidths: minimum field widths for the values 

2240 

2241 Each value in in ``row`` is converted to a string and optionally padded 

2242 with blanks. The resulting strings are output separated with blanks. If 

2243 any values given are strings and if they contain whitespace, they are 

2244 quoted with single quotes, and any internal single quotes are 

2245 backslash-escaped. 

2246 ''' 

2247 

2248 out = [] 

2249 

2250 for i, x in enumerate(row): 

2251 w = 0 

2252 if minfieldwidths and i < len(minfieldwidths): 

2253 w = minfieldwidths[i] 

2254 

2255 if isinstance(x, str): 

2256 if re.search(r"\s|'", x): 

2257 x = "'%s'" % escapequotes(x) 

2258 

2259 x = x.ljust(w) 

2260 else: 

2261 x = str(x).rjust(w) 

2262 

2263 out.append(x) 

2264 

2265 self._f.write(' '.join(out).rstrip() + '\n') 

2266 

2267 

2268class TableReader(object): 

2269 

2270 ''' 

2271 Read table of space separated values from a file. 

2272 

2273 :param f: file-like object 

2274 

2275 This uses Pythons shlex module to tokenize lines. Should deal correctly 

2276 with quoted strings. 

2277 ''' 

2278 

2279 def __init__(self, f): 

2280 self._f = f 

2281 self.eof = False 

2282 

2283 def readrow(self): 

2284 ''' 

2285 Read one row from the underlying file, tokenize it with shlex. 

2286 

2287 :returns: tokenized line as a list of strings. 

2288 ''' 

2289 

2290 line = self._f.readline() 

2291 if not line: 

2292 self.eof = True 

2293 return [] 

2294 line.strip() 

2295 if line.startswith('#'): 

2296 return [] 

2297 

2298 return qsplit(line) 

2299 

2300 

2301def gform(number, significant_digits=3): 

2302 ''' 

2303 Pretty print floating point numbers. 

2304 

2305 Align floating point numbers at the decimal dot. 

2306 

2307 :: 

2308 

2309 | -d.dde+xxx| 

2310 | -d.dde+xx | 

2311 |-ddd. | 

2312 | -dd.d | 

2313 | -d.dd | 

2314 | -0.ddd | 

2315 | -0.0ddd | 

2316 | -0.00ddd | 

2317 | -d.dde-xx | 

2318 | -d.dde-xxx| 

2319 | nan| 

2320 

2321 

2322 The formatted string has length ``significant_digits * 2 + 6``. 

2323 ''' 

2324 

2325 no_exp_range = (pow(10., -1), 

2326 pow(10., significant_digits)) 

2327 width = significant_digits+significant_digits-1+1+1+5 

2328 

2329 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.: 

2330 s = ('%#.*g' % (significant_digits, number)).rstrip('0') 

2331 else: 

2332 s = '%.*E' % (significant_digits-1, number) 

2333 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width) 

2334 if s.strip().lower() == 'nan': 

2335 s = 'nan'.rjust(width) 

2336 return s 

2337 

2338 

2339def human_bytesize(value): 

2340 

2341 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split() 

2342 

2343 if value == 1: 

2344 return '1 Byte' 

2345 

2346 for i, ext in enumerate(exts): 

2347 x = float(value) / 1000**i 

2348 if round(x) < 10. and not value < 1000: 

2349 return '%.1f %s' % (x, ext) 

2350 if round(x) < 1000.: 

2351 return '%.0f %s' % (x, ext) 

2352 

2353 return '%i Bytes' % value 

2354 

2355 

2356re_compatibility = re.compile( 

2357 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' + 

2358 r'(dummy|poel|qseis|qssp))\.' 

2359) 

2360 

2361 

2362def pf_is_old(fn): 

2363 oldstyle = False 

2364 with open(fn, 'r') as f: 

2365 for line in f: 

2366 if re_compatibility.search(line): 

2367 oldstyle = True 

2368 

2369 return oldstyle 

2370 

2371 

2372def pf_upgrade(fn): 

2373 need = pf_is_old(fn) 

2374 if need: 

2375 fn_temp = fn + '.temp' 

2376 

2377 with open(fn, 'r') as fin: 

2378 with open(fn_temp, 'w') as fout: 

2379 for line in fin: 

2380 line = re_compatibility.sub('!pf.', line) 

2381 fout.write(line) 

2382 

2383 os.rename(fn_temp, fn) 

2384 

2385 return need 

2386 

2387 

2388def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'): 

2389 ''' 

2390 Extract leap second information from tzdata. 

2391 

2392 Based on example at http://stackoverflow.com/questions/19332902/\ 

2393 extract-historic-leap-seconds-from-tzdata 

2394 

2395 See also 'man 5 tzfile'. 

2396 ''' 

2397 

2398 from struct import unpack, calcsize 

2399 out = [] 

2400 with open(tzfile, 'rb') as f: 

2401 # read header 

2402 fmt = '>4s c 15x 6l' 

2403 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt, 

2404 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt))) 

2405 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file' 

2406 

2407 # skip over some uninteresting data 

2408 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict( 

2409 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt) 

2410 f.read(calcsize(fmt)) 

2411 

2412 # read leap-seconds 

2413 fmt = '>2l' 

2414 for i in range(leapcnt): 

2415 tleap, nleap = unpack(fmt, f.read(calcsize(fmt))) 

2416 out.append((tleap-nleap+1, nleap)) 

2417 

2418 return out 

2419 

2420 

2421class LeapSecondsError(Exception): 

2422 pass 

2423 

2424 

2425class LeapSecondsOutdated(LeapSecondsError): 

2426 pass 

2427 

2428 

2429class InvalidLeapSecondsFile(LeapSecondsOutdated): 

2430 pass 

2431 

2432 

2433def parse_leap_seconds_list(fn): 

2434 data = [] 

2435 texpires = None 

2436 try: 

2437 t0 = int(round(str_to_time('1900-01-01 00:00:00'))) 

2438 except TimeStrError: 

2439 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800 

2440 

2441 tnow = int(round(time.time())) 

2442 

2443 if not op.exists(fn): 

2444 raise LeapSecondsOutdated('no leap seconds file found') 

2445 

2446 try: 

2447 with open(fn, 'rb') as f: 

2448 for line in f: 

2449 if line.strip().startswith(b'<!DOCTYPE'): 

2450 raise InvalidLeapSecondsFile('invalid leap seconds file') 

2451 

2452 if line.startswith(b'#@'): 

2453 texpires = int(line.split()[1]) + t0 

2454 elif line.startswith(b'#') or len(line) < 5: 

2455 pass 

2456 else: 

2457 toks = line.split() 

2458 t = int(toks[0]) + t0 

2459 nleap = int(toks[1]) - 10 

2460 data.append((t, nleap)) 

2461 

2462 except IOError: 

2463 raise LeapSecondsError('cannot read leap seconds file %s' % fn) 

2464 

2465 if texpires is None or tnow > texpires: 

2466 raise LeapSecondsOutdated('leap seconds list is outdated') 

2467 

2468 return data 

2469 

2470 

2471def read_leap_seconds2(): 

2472 from pyrocko import config 

2473 conf = config.config() 

2474 fn = conf.leapseconds_path 

2475 url = conf.leapseconds_url 

2476 # check for outdated default URL 

2477 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list': 

2478 url = 'https://www.ietf.org/timezones/data/leap-seconds.list' 

2479 logger.info( 

2480 'Leap seconds default URL is now: %s\nUsing new default.' % url) 

2481 

2482 for i in range(3): 

2483 try: 

2484 return parse_leap_seconds_list(fn) 

2485 

2486 except LeapSecondsOutdated: 

2487 try: 

2488 logger.info('updating leap seconds list...') 

2489 download_file(url, fn) 

2490 

2491 except Exception as e: 

2492 raise LeapSecondsError( 

2493 'cannot download leap seconds list from %s to %s (%s)' 

2494 % (url, fn, e)) 

2495 

2496 raise LeapSecondsError('Could not retrieve/read leap seconds file.') 

2497 

2498 

2499def gps_utc_offset(t_utc): 

2500 ''' 

2501 Time offset t_gps - t_utc for a given t_utc. 

2502 ''' 

2503 ls = read_leap_seconds2() 

2504 i = 0 

2505 if t_utc < ls[0][0]: 

2506 return ls[0][1] - 1 - 9 

2507 

2508 while i < len(ls) - 1: 

2509 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]: 

2510 return ls[i][1] - 9 

2511 i += 1 

2512 

2513 return ls[-1][1] - 9 

2514 

2515 

2516def utc_gps_offset(t_gps): 

2517 ''' 

2518 Time offset t_utc - t_gps for a given t_gps. 

2519 ''' 

2520 ls = read_leap_seconds2() 

2521 

2522 if t_gps < ls[0][0] + ls[0][1] - 9: 

2523 return - (ls[0][1] - 1 - 9) 

2524 

2525 i = 0 

2526 while i < len(ls) - 1: 

2527 if ls[i][0] + ls[i][1] - 9 <= t_gps \ 

2528 and t_gps < ls[i+1][0] + ls[i+1][1] - 9: 

2529 return - (ls[i][1] - 9) 

2530 i += 1 

2531 

2532 return - (ls[-1][1] - 9) 

2533 

2534 

2535def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'): 

2536 import itertools 

2537 import glob 

2538 from pyrocko.io.io_common import FileLoadError 

2539 

2540 def iload_filename(filename, **kwargs): 

2541 try: 

2542 with open(filename, 'rb') as f: 

2543 for cr in iload_fh(f, **kwargs): 

2544 yield cr 

2545 

2546 except FileLoadError as e: 

2547 e.set_context('filename', filename) 

2548 raise 

2549 

2550 def iload_dirname(dirname, **kwargs): 

2551 for entry in os.listdir(dirname): 

2552 fpath = op.join(dirname, entry) 

2553 if op.isfile(fpath): 

2554 for cr in iload_filename(fpath, **kwargs): 

2555 yield cr 

2556 

2557 def iload_glob(pattern, **kwargs): 

2558 

2559 for fn in glob.iglob(pattern): 

2560 for cr in iload_filename(fn, **kwargs): 

2561 yield cr 

2562 

2563 def iload(source, **kwargs): 

2564 if isinstance(source, str): 

2565 if op.isdir(source): 

2566 return iload_dirname(source, **kwargs) 

2567 elif op.isfile(source): 

2568 return iload_filename(source, **kwargs) 

2569 else: 

2570 return iload_glob(source, **kwargs) 

2571 

2572 elif hasattr(source, 'read'): 

2573 return iload_fh(source, **kwargs) 

2574 else: 

2575 return itertools.chain.from_iterable( 

2576 iload(subsource, **kwargs) for subsource in source) 

2577 

2578 iload_filename.__doc__ = ''' 

2579 Read %s information from named file. 

2580 ''' % doc_fmt 

2581 

2582 iload_dirname.__doc__ = ''' 

2583 Read %s information from directory of %s files. 

2584 ''' % (doc_fmt, doc_fmt) 

2585 

2586 iload_glob.__doc__ = ''' 

2587 Read %s information from files matching a glob pattern. 

2588 ''' % doc_fmt 

2589 

2590 iload.__doc__ = ''' 

2591 Load %s information from given source(s) 

2592 

2593 The ``source`` can be specified as the name of a %s file, the name of a 

2594 directory containing %s files, a glob pattern of %s files, an open 

2595 filehandle or an iterator yielding any of the forementioned sources. 

2596 

2597 This function behaves as a generator yielding %s objects. 

2598 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects) 

2599 

2600 for f in iload_filename, iload_dirname, iload_glob, iload: 

2601 f.__module__ = iload_fh.__module__ 

2602 

2603 return iload_filename, iload_dirname, iload_glob, iload 

2604 

2605 

2606class Inconsistency(Exception): 

2607 pass 

2608 

2609 

2610def consistency_check(list_of_tuples, message='values differ:'): 

2611 ''' 

2612 Check for inconsistencies. 

2613 

2614 Given a list of tuples, check that all tuple elements except for first one 

2615 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be 

2616 valid because the coordinates at the two channels are the same. 

2617 ''' 

2618 

2619 if len(list_of_tuples) >= 2: 

2620 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]): 

2621 raise Inconsistency('%s\n' % message + '\n'.join( 

2622 ' %s: %s' % (t[0], ', '.join(str(x) for x in t[1:])) 

2623 for t in list_of_tuples)) 

2624 

2625 

2626class defaultzerodict(dict): 

2627 def __missing__(self, k): 

2628 return 0 

2629 

2630 

2631def mostfrequent(x): 

2632 c = defaultzerodict() 

2633 for e in x: 

2634 c[e] += 1 

2635 

2636 return sorted(list(c.keys()), key=lambda k: c[k])[-1] 

2637 

2638 

2639def consistency_merge(list_of_tuples, 

2640 message='values differ:', 

2641 error='raise', 

2642 merge=mostfrequent): 

2643 

2644 assert error in ('raise', 'warn', 'ignore') 

2645 

2646 if len(list_of_tuples) == 0: 

2647 raise Exception('cannot merge empty sequence') 

2648 

2649 try: 

2650 consistency_check(list_of_tuples, message) 

2651 return list_of_tuples[0][1:] 

2652 except Inconsistency as e: 

2653 if error == 'raise': 

2654 raise 

2655 

2656 elif error == 'warn': 

2657 logger.warning(str(e)) 

2658 

2659 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]]) 

2660 

2661 

2662def short_to_list(nmax, it): 

2663 import itertools 

2664 

2665 if isinstance(it, list): 

2666 return it 

2667 

2668 li = [] 

2669 for i in range(nmax+1): 

2670 try: 

2671 li.append(next(it)) 

2672 except StopIteration: 

2673 return li 

2674 

2675 return itertools.chain(li, it) 

2676 

2677 

2678def parse_md(f): 

2679 try: 

2680 with open(op.join( 

2681 op.dirname(op.abspath(f)), 

2682 'README.md'), 'r') as readme: 

2683 mdstr = readme.read() 

2684 except IOError as e: 

2685 return 'Failed to get README.md: %s' % e 

2686 

2687 # Remve the title 

2688 mdstr = re.sub(r'^# .*\n?', '', mdstr) 

2689 # Append sphinx reference to `pyrocko.` modules 

2690 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr) 

2691 # Convert Subsections to toc-less rubrics 

2692 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr) 

2693 return mdstr 

2694 

2695 

2696def mpl_show(plt): 

2697 import matplotlib 

2698 if matplotlib.get_backend().lower() == 'agg': 

2699 logger.warning('Cannot show() when using matplotlib "agg" backend') 

2700 else: 

2701 plt.show() 

2702 

2703 

2704g_re_qsplit = re.compile( 

2705 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)') 

2706g_re_qsplit_sep = {} 

2707 

2708 

2709def get_re_qsplit(sep): 

2710 if sep is None: 

2711 return g_re_qsplit 

2712 else: 

2713 if sep not in g_re_qsplit_sep: 

2714 assert len(sep) == 1 

2715 assert sep not in '\'"' 

2716 esep = re.escape(sep) 

2717 g_re_qsplit_sep[sep] = re.compile( 

2718 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|' 

2719 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa 

2720 return g_re_qsplit_sep[sep] 

2721 

2722 

2723g_re_trivial = re.compile(r'\A[^\'"\s]+\Z') 

2724g_re_trivial_sep = {} 

2725 

2726 

2727def get_re_trivial(sep): 

2728 if sep is None: 

2729 return g_re_trivial 

2730 else: 

2731 if sep not in g_re_qsplit_sep: 

2732 assert len(sep) == 1 

2733 assert sep not in '\'"' 

2734 esep = re.escape(sep) 

2735 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z') 

2736 

2737 return g_re_trivial_sep[sep] 

2738 

2739 

2740g_re_escape_s = re.compile(r'([\\\'])') 

2741g_re_unescape_s = re.compile(r'\\([\\\'])') 

2742g_re_escape_d = re.compile(r'([\\"])') 

2743g_re_unescape_d = re.compile(r'\\([\\"])') 

2744 

2745 

2746def escape_s(s): 

2747 ''' 

2748 Backslash-escape single-quotes and backslashes. 

2749 

2750 Example: ``Jack's`` => ``Jack\\'s`` 

2751 

2752 ''' 

2753 return g_re_escape_s.sub(r'\\\1', s) 

2754 

2755 

2756def unescape_s(s): 

2757 ''' 

2758 Unescape backslash-escaped single-quotes and backslashes. 

2759 

2760 Example: ``Jack\\'s`` => ``Jack's`` 

2761 ''' 

2762 return g_re_unescape_s.sub(r'\1', s) 

2763 

2764 

2765def escape_d(s): 

2766 ''' 

2767 Backslash-escape double-quotes and backslashes. 

2768 

2769 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"`` 

2770 ''' 

2771 return g_re_escape_d.sub(r'\\\1', s) 

2772 

2773 

2774def unescape_d(s): 

2775 ''' 

2776 Unescape backslash-escaped double-quotes and backslashes. 

2777 

2778 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"`` 

2779 ''' 

2780 return g_re_unescape_d.sub(r'\1', s) 

2781 

2782 

2783def qjoin_s(it, sep=None): 

2784 ''' 

2785 Join sequence of strings into a line, single-quoting non-trivial strings. 

2786 

2787 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"`` 

2788 ''' 

2789 re_trivial = get_re_trivial(sep) 

2790 

2791 if sep is None: 

2792 sep = ' ' 

2793 

2794 return sep.join( 

2795 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it) 

2796 

2797 

2798def qjoin_d(it, sep=None): 

2799 ''' 

2800 Join sequence of strings into a line, double-quoting non-trivial strings. 

2801 

2802 Example: ``['55', 'Pete "The Robot" Smith']`` => 

2803 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` 

2804 ''' 

2805 re_trivial = get_re_trivial(sep) 

2806 if sep is None: 

2807 sep = ' ' 

2808 

2809 return sep.join( 

2810 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it) 

2811 

2812 

2813def qsplit(s, sep=None): 

2814 ''' 

2815 Split line into list of strings, allowing for quoted strings. 

2816 

2817 Example: ``"55 'Sparrow\\\\'s Island'"`` => 

2818 ``["55", "Sparrow's Island"]``, 

2819 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` => 

2820 ``['55', 'Pete "The Robot" Smith']`` 

2821 ''' 

2822 re_qsplit = get_re_qsplit(sep) 

2823 return [ 

2824 (unescape_d(x[0]) or unescape_s(x[1]) or x[2]) 

2825 for x in re_qsplit.findall(s)] 

2826 

2827 

2828g_have_warned_threadpoolctl = False 

2829 

2830 

2831class threadpool_limits_dummy(object): 

2832 

2833 def __init__(self, *args, **kwargs): 

2834 pass 

2835 

2836 def __enter__(self): 

2837 global g_have_warned_threadpoolctl 

2838 

2839 if not g_have_warned_threadpoolctl: 

2840 logger.warning( 

2841 'Cannot control number of BLAS threads because ' 

2842 '`threadpoolctl` module is not available. You may want to ' 

2843 'install `threadpoolctl`.') 

2844 

2845 g_have_warned_threadpoolctl = True 

2846 

2847 return self 

2848 

2849 def __exit__(self, type, value, traceback): 

2850 pass 

2851 

2852 

2853def get_threadpool_limits(): 

2854 ''' 

2855 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail. 

2856 ''' 

2857 

2858 try: 

2859 from threadpoolctl import threadpool_limits 

2860 return threadpool_limits 

2861 

2862 except ImportError: 

2863 return threadpool_limits_dummy