Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/util.py: 75%
1305 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6Utility functions for Pyrocko.
8.. _time-handling-mode:
10High precision time handling mode
11.................................
13Pyrocko can treat timestamps either as standard double precision (64 bit)
14floating point values, or as high precision floats (``numpy.float128`` or
15``numpy.float96``, whichever is available, `see NumPy Scalars
16<https://numpy.org/doc/stable/reference/arrays.scalars.html>`_), aliased here
17as :py:class:`~pyrocko.util.hpfloat`. High precision time stamps are required
18when handling data with sub-millisecond precision, i.e. kHz/MHz data streams
19and event catalogs derived from such data.
21Not all functions in Pyrocko and in programs depending on Pyrocko may work
22correctly with high precision times. Therefore, Pyrocko's high precision time
23handling mode has to be actively activated by user config, command line option
24or enforced within a certain script/program.
26The default high precision time handling mode can be configured globally with
27the user configuration variable
28:py:gattr:`~pyrocko.config.PyrockoConfig.use_high_precision_time`. Calling the
29function :py:func:`use_high_precision_time` overrides the default from the
30config file. This function may be called at startup of a program/script
31requiring a specific time handling mode.
33To create a valid time stamp for use in Pyrocko (e.g. in
34:py:class:`~pyrocko.model.event.Event` or
35:py:class:`~pyrocko.trace.Trace` objects), use:
37.. code-block :: python
39 import time
40 from pyrocko import util
42 # By default using mode selected in user config, override with:
43 # util.use_high_precision_time(True) # force high precision mode
44 # util.use_high_precision_time(False) # force low precision mode
46 t1 = util.str_to_time('2020-08-27 10:22:00')
47 t2 = util.str_to_time('2020-08-27 10:22:00.111222')
48 t3 = util.to_time_float(time.time())
50 # To get the appropriate float class, use:
52 time_float = util.get_time_float()
53 # -> float, numpy.float128 or numpy.float96
54 [isinstance(t, time_float) for t in [t1, t2, t3]]
55 # -> [True, True, True]
57 # Shortcut:
58 util.check_time_class(t1)
60Module content
61..............
63.. py:class:: hpfloat
65 Alias for NumPy's high precision float data type ``float128`` or
66 ``float96``, if available.
68 On platforms lacking support for high precision floats, an attempt to
69 create a ``hpfloat`` instance, raises :py:exc:`HPFloatUnavailable`.
71'''
73import time
74import logging
75import os
76import sys
77import re
78import calendar
79import math
80import fnmatch
81import inspect
82import warnings
83import weakref
84import signal as psignal
85try:
86 import fcntl
87except ImportError:
88 fcntl = None
89import optparse
90import os.path as op
91import errno
92from collections import defaultdict
94import numpy as num
95from scipy import signal
96import pyrocko
97from pyrocko import dummy_progressbar
100from urllib.parse import urlencode, quote, unquote # noqa
101from urllib.request import Request, build_opener, HTTPDigestAuthHandler # noqa
102from urllib.request import urlopen as _urlopen # noqa
103from urllib.error import HTTPError, URLError # noqa
106try:
107 import certifi
108 import ssl
109 g_ssl_context = ssl.create_default_context(cafile=certifi.where())
110except ImportError:
111 g_ssl_context = None
114class URLErrorSSL(URLError):
116 def __init__(self, urlerror):
117 self.__dict__ = urlerror.__dict__.copy()
119 def __str__(self):
120 return (
121 'Requesting web resource failed and the problem could be '
122 'related to SSL. Python standard libraries on some older '
123 'systems (like Ubuntu 14.04) are known to have trouble '
124 "with some SSL setups of today's servers: %s"
125 % URLError.__str__(self))
128def urlopen_ssl_check(*args, **kwargs):
129 try:
130 return urlopen(*args, **kwargs)
131 except URLError as e:
132 if str(e).find('SSL') != -1:
133 raise URLErrorSSL(e)
134 else:
135 raise
138def urlopen(*args, **kwargs):
140 if 'context' not in kwargs and g_ssl_context is not None:
141 kwargs['context'] = g_ssl_context
143 return _urlopen(*args, **kwargs)
146try:
147 long
148except NameError:
149 long = int
152force_dummy_progressbar = False
155try:
156 from pyrocko import util_ext
157except ImportError:
158 util_ext = None
161logger = logging.getLogger('pyrocko.util')
164# fallbacks num_full and num_full_like are not needed anymore but
165# kept here because downstream code may still use these.
166try:
167 num_full = num.full
168except AttributeError:
169 def num_full(shape, fill_value, dtype=None, order='C'):
170 a = num.empty(shape, dtype=dtype, order=order)
171 a.fill(fill_value)
172 return a
174try:
175 num_full_like = num.full_like
176except AttributeError:
177 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
178 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
179 a.fill(fill_value)
180 return a
183g_setup_logging_args = 'pyrocko', 'warning'
186def setup_logging(programname='pyrocko', levelname='warning'):
187 '''
188 Initialize logging.
190 :param programname: program name to be written in log
191 :param levelname: string indicating the logging level ('debug', 'info',
192 'warning', 'error', 'critical')
194 This function is called at startup by most pyrocko programs to set up a
195 consistent logging format. This is simply a shortcut to a call to
196 :py:func:`logging.basicConfig()`.
197 '''
199 global g_setup_logging_args
200 g_setup_logging_args = (programname, levelname)
202 levels = {'debug': logging.DEBUG,
203 'info': logging.INFO,
204 'warning': logging.WARNING,
205 'error': logging.ERROR,
206 'critical': logging.CRITICAL}
208 logging.basicConfig(
209 level=levels[levelname],
210 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s')
213def subprocess_setup_logging_args():
214 '''
215 Get arguments from previous call to setup_logging.
217 These can be sent down to a worker process so it can setup its logging
218 in the same way as the main process.
219 '''
220 return g_setup_logging_args
223g_experimental_features_used = set()
226def experimental_feature_used(feature):
227 '''
228 Warn once that an experimental feature is used.
229 '''
231 if feature not in g_experimental_features_used:
232 warnings.warn(
233 'Experimental feature used: %s' % feature,
234 stacklevel=2)
235 g_experimental_features_used.add(feature)
238def data_file(fn):
239 return os.path.join(os.path.split(__file__)[0], 'data', fn)
242class DownloadError(Exception):
243 '''
244 Raised when a download failed.
245 '''
246 pass
249class PathExists(DownloadError):
250 '''
251 Raised when the download target file already exists.
252 '''
253 pass
256def _download(url, fpath, username=None, password=None,
257 force=False, method='download', stats=None,
258 status_callback=None, entries_wanted=None,
259 recursive=False, header=None):
261 import requests
262 from requests.auth import HTTPBasicAuth
263 from requests.exceptions import HTTPError as req_HTTPError
265 requests.adapters.DEFAULT_RETRIES = 5
266 urljoin = requests.compat.urljoin
268 session = requests.Session()
269 session.header = header
270 session.auth = None if username is None\
271 else HTTPBasicAuth(username, password)
273 status = {
274 'ntotal_files': 0,
275 'nread_files': 0,
276 'ntotal_bytes_all_files': 0,
277 'nread_bytes_all_files': 0,
278 'ntotal_bytes_current_file': 0,
279 'nread_bytes_current_file': 0,
280 'finished': False
281 }
283 try:
284 url_to_size = {}
286 if callable(status_callback):
287 status_callback(status)
289 if not recursive and url.endswith('/'):
290 raise DownloadError(
291 'URL: %s appears to be a directory'
292 ' but recurvise download is False' % url)
294 if recursive and not url.endswith('/'):
295 url += '/'
297 def parse_directory_tree(url, subdir=''):
298 r = session.get(urljoin(url, subdir))
299 r.raise_for_status()
301 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
303 files = sorted(set(subdir + fn for fn in entries
304 if not fn.endswith('/')))
306 if entries_wanted is not None:
307 files = [fn for fn in files
308 if (fn in wanted for wanted in entries_wanted)]
310 status['ntotal_files'] += len(files)
312 dirs = sorted(set(subdir + dn for dn in entries
313 if dn.endswith('/')
314 and dn not in ('./', '../')))
316 for dn in dirs:
317 files.extend(parse_directory_tree(
318 url, subdir=dn))
320 return files
322 def get_content_length(url):
323 if url not in url_to_size:
324 r = session.head(url, headers={'Accept-Encoding': ''})
326 content_length = r.headers.get('content-length', None)
327 if content_length is None:
328 logger.debug('Could not get HTTP header '
329 'Content-Length for %s' % url)
331 content_length = None
333 else:
334 content_length = int(content_length)
335 status['ntotal_bytes_all_files'] += content_length
336 if callable(status_callback):
337 status_callback(status)
339 url_to_size[url] = content_length
341 return url_to_size[url]
343 def download_file(url, fn):
344 logger.info('starting download of %s...' % url)
345 ensuredirs(fn)
347 fsize = get_content_length(url)
348 status['ntotal_bytes_current_file'] = fsize
349 status['nread_bytes_current_file'] = 0
350 if callable(status_callback):
351 status_callback(status)
353 r = session.get(url, stream=True, timeout=5)
354 r.raise_for_status()
356 frx = 0
357 fn_tmp = fn + '.%i.temp' % os.getpid()
358 with open(fn_tmp, 'wb') as f:
359 for d in r.iter_content(chunk_size=1024):
360 f.write(d)
361 frx += len(d)
363 status['nread_bytes_all_files'] += len(d)
364 status['nread_bytes_current_file'] += len(d)
365 if callable(status_callback):
366 status_callback(status)
368 os.rename(fn_tmp, fn)
370 if fsize is not None and frx != fsize:
371 logger.warning(
372 'HTTP header Content-Length: %i bytes does not match '
373 'download size %i bytes' % (fsize, frx))
375 logger.info('finished download of %s' % url)
377 status['nread_files'] += 1
378 if callable(status_callback):
379 status_callback(status)
381 if recursive:
382 if op.exists(fpath) and not force:
383 raise PathExists('path %s already exists' % fpath)
385 files = parse_directory_tree(url)
387 dsize = 0
388 for fn in files:
389 file_url = urljoin(url, fn)
390 dsize += get_content_length(file_url) or 0
392 if method == 'calcsize':
393 return dsize
395 else:
396 for fn in files:
397 file_url = urljoin(url, fn)
398 download_file(file_url, op.join(fpath, fn))
400 else:
401 status['ntotal_files'] += 1
402 if callable(status_callback):
403 status_callback(status)
405 fsize = get_content_length(url)
406 if method == 'calcsize':
407 return fsize
408 else:
409 download_file(url, fpath)
411 except req_HTTPError as e:
412 logging.warning('http error: %s' % e)
413 raise DownloadError('could not download file(s) from: %s' % url)
415 finally:
416 status['finished'] = True
417 if callable(status_callback):
418 status_callback(status)
419 session.close()
422def download_file(
423 url, fpath, username=None, password=None, status_callback=None,
424 **kwargs):
425 return _download(
426 url, fpath, username, password,
427 recursive=False,
428 status_callback=status_callback,
429 **kwargs)
432def download_dir(
433 url, fpath, username=None, password=None, status_callback=None,
434 **kwargs):
436 return _download(
437 url, fpath, username, password,
438 recursive=True,
439 status_callback=status_callback,
440 **kwargs)
443class HPFloatUnavailable(Exception):
444 '''
445 Raised when a high precision float type would be required but is not
446 available.
447 '''
448 pass
451class dummy_hpfloat(object):
452 def __init__(self, *args, **kwargs):
453 raise HPFloatUnavailable(
454 'NumPy lacks support for float128 or float96 data type on this '
455 'platform.')
458if hasattr(num, 'float128'):
459 hpfloat = num.float128
461elif hasattr(num, 'float96'):
462 hpfloat = num.float96
464else:
465 hpfloat = dummy_hpfloat
468g_time_float = None
469g_time_dtype = None
472class TimeFloatSettingError(Exception):
473 pass
476def use_high_precision_time(enabled):
477 '''
478 Globally force a specific time handling mode.
480 See :ref:`High precision time handling mode <time-handling-mode>`.
482 :param enabled: enable/disable use of high precision time type
483 :type enabled: bool
485 This function should be called before handling/reading any time data.
486 It can only be called once.
488 Special attention is required when using multiprocessing on a platform
489 which does not use fork under the hood. In such cases, the desired setting
490 must be set also in the subprocess.
491 '''
492 _setup_high_precision_time_mode(enabled_app=enabled)
495def _setup_high_precision_time_mode(enabled_app=False):
496 global g_time_float
497 global g_time_dtype
499 if not (g_time_float is None and g_time_dtype is None):
500 raise TimeFloatSettingError(
501 'Cannot set time handling mode: too late, it has already been '
502 'fixed by an earlier call.')
504 from pyrocko import config
506 conf = config.config()
507 enabled_config = conf.use_high_precision_time
509 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
510 if enabled_env is not None:
511 try:
512 enabled_env = int(enabled_env) == 1
513 except ValueError:
514 raise TimeFloatSettingError(
515 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
516 'should be set to 0 or 1.')
518 enabled = enabled_config
519 mode_from = 'config variable `use_high_precision_time`'
520 notify = enabled
522 if enabled_env is not None:
523 if enabled_env != enabled:
524 notify = True
525 enabled = enabled_env
526 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
528 if enabled_app is not None:
529 if enabled_app != enabled:
530 notify = True
531 enabled = enabled_app
532 mode_from = 'application override'
534 logger.debug('''
535Pyrocko high precision time mode selection (latter override earlier):
536 config: %s
537 env: %s
538 app: %s
539 -> enabled: %s'''.lstrip() % (
540 enabled_config, enabled_env, enabled_app, enabled))
542 if notify:
543 logger.info('Pyrocko high precision time mode %s by %s.' % (
544 'activated' if enabled else 'deactivated',
545 mode_from))
547 if enabled:
548 g_time_float = hpfloat
549 g_time_dtype = hpfloat
550 else:
551 g_time_float = float
552 g_time_dtype = num.float64
555def get_time_float():
556 '''
557 Get the effective float class for timestamps.
559 See :ref:`High precision time handling mode <time-handling-mode>`.
561 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
562 current time handling mode
563 '''
564 global g_time_float # noqa
566 if g_time_float is None:
567 _setup_high_precision_time_mode()
569 return g_time_float
572def get_time_dtype():
573 '''
574 Get effective NumPy float class to handle timestamps.
576 See :ref:`High precision time handling mode <time-handling-mode>`.
577 '''
579 global g_time_dtype # noqa
581 if g_time_dtype is None:
582 _setup_high_precision_time_mode()
584 return g_time_dtype
587def to_time_float(t):
588 '''
589 Convert float to valid time stamp in the current time handling mode.
591 See :ref:`High precision time handling mode <time-handling-mode>`.
592 '''
593 return get_time_float()(t)
596class TimestampTypeError(ValueError):
597 pass
600def check_time_class(t, error='raise'):
601 '''
602 Type-check variable against current time handling mode.
604 See :ref:`High precision time handling mode <time-handling-mode>`.
605 '''
607 if t == 0.0:
608 return
610 if not isinstance(t, get_time_float()):
611 message = (
612 'Timestamp %g is of type %s but should be of type %s with '
613 "Pyrocko's currently selected time handling mode.\n\n"
614 'See https://pyrocko.org/docs/current/library/reference/util.html'
615 '#high-precision-time-handling-mode' % (
616 t, type(t), get_time_float()))
618 if error == 'raise':
619 raise TimestampTypeError(message)
620 elif error == 'warn':
621 logger.warning(message)
622 else:
623 assert False
626def iter_windows(
627 tmin=None,
628 tmax=None,
629 tinc=None,
630 tpad=0.0,
631 snap_window=False,
632 tmin_content=None,
633 tmax_content=None):
635 from pyrocko.squirrel import Batch
637 if snap_window and tinc is not None:
638 tmin = tmin if tmin is not None else tmin_content
639 tmax = tmax if tmax is not None else tmax_content
640 tmin = math.floor(tmin / tinc) * tinc
641 tmax = math.ceil(tmax / tinc) * tinc
642 else:
643 tmin = tmin if tmin is not None else (
644 tmin_content + tpad if tmin_content is not None else None)
645 tmax = tmax if tmax is not None else (
646 tmax_content - tpad if tmax_content is not None else None)
648 if None in (tmin, tmax):
649 return
651 if tinc is None:
652 tinc = tmax - tmin
653 nwin = 1
654 elif tinc == 0.0:
655 nwin = 1
656 else:
657 eps = 1e-6
658 nwin = max(1, int((tmax - tmin) / tinc - eps) + 1)
660 for iwin in range(nwin):
661 wmin, wmax = tmin+iwin*tinc, min(tmin+(iwin+1)*tinc, tmax)
663 yield Batch(
664 tmin=wmin,
665 tmax=wmax,
666 tpad=tpad,
667 i=iwin,
668 n=nwin,
669 igroup=0,
670 ngroups=1,
671 traces=[])
674class Stopwatch(object):
675 '''
676 Simple stopwatch to measure elapsed wall clock time.
678 Usage::
680 s = Stopwatch()
681 time.sleep(1)
682 print s()
683 time.sleep(1)
684 print s()
685 '''
687 def __init__(self):
688 self.start = time.time()
690 def __call__(self):
691 return time.time() - self.start
694def wrap(text, line_length=80):
695 '''
696 Paragraph and list-aware wrapping of text.
697 '''
699 text = text.strip('\n')
700 at_lineend = re.compile(r' *\n')
701 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
703 paragraphs = at_para.split(text)[::5]
704 listindents = at_para.split(text)[4::5]
705 newlist = at_para.split(text)[3::5]
707 listindents[0:0] = [None]
708 listindents.append(True)
709 newlist.append(None)
711 det_indent = re.compile(r'^ *')
713 outlines = []
714 for ip, p in enumerate(paragraphs):
715 if not p:
716 continue
718 if listindents[ip] is None:
719 _indent = det_indent.findall(p)[0]
720 findent = _indent
721 else:
722 findent = listindents[ip]
723 _indent = ' ' * len(findent)
725 ll = line_length - len(_indent)
726 llf = ll
728 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
729 p1 = ' '.join(oldlines)
730 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
731 for imatch, match in enumerate(possible.finditer(p1)):
732 parout = match.group(1)
733 if imatch == 0:
734 outlines.append(findent + parout)
735 else:
736 outlines.append(_indent + parout)
738 if ip != len(paragraphs)-1 and (
739 listindents[ip] is None
740 or newlist[ip] is not None
741 or listindents[ip+1] is None):
743 outlines.append('')
745 return outlines
748def ewrap(lines, width=80, indent=''):
749 lines = list(lines)
750 if not lines:
751 return ''
752 fwidth = max(len(s) for s in lines)
753 nx = max(1, (80-len(indent)) // (fwidth+1))
754 i = 0
755 rows = []
756 while i < len(lines):
757 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx]))
758 i += nx
760 return '\n'.join(rows)
763class BetterHelpFormatter(optparse.IndentedHelpFormatter):
765 def __init__(self, *args, **kwargs):
766 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
768 def format_option(self, option):
769 '''
770 From IndentedHelpFormatter but using a different wrap method.
771 '''
773 help_text_position = 4 + self.current_indent
774 help_text_width = self.width - help_text_position
776 opts = self.option_strings[option]
777 opts = '%*s%s' % (self.current_indent, '', opts)
778 if option.help:
779 help_text = self.expand_default(option)
781 if self.help_position + len(help_text) + 1 <= self.width:
782 lines = [
783 '',
784 '%-*s %s' % (self.help_position, opts, help_text),
785 '']
786 else:
787 lines = ['']
788 lines.append(opts)
789 lines.append('')
790 if option.help:
791 help_lines = wrap(help_text, help_text_width)
792 lines.extend(['%*s%s' % (help_text_position, '', line)
793 for line in help_lines])
794 lines.append('')
796 return '\n'.join(lines)
798 def format_description(self, description):
799 if not description:
800 return ''
802 if self.current_indent == 0:
803 lines = []
804 else:
805 lines = ['']
807 lines.extend(wrap(description, self.width - self.current_indent))
808 if self.current_indent == 0:
809 lines.append('\n')
811 return '\n'.join(
812 ['%*s%s' % (self.current_indent, '', line) for line in lines])
815class ProgressBar:
816 def __init__(self, label, n):
817 from pyrocko import progress
818 self._context = progress.view()
819 self._context.__enter__()
820 self._task = progress.task(label, n)
822 def update(self, i):
823 self._task.update(i)
825 def finish(self):
826 self._task.done()
827 if self._context:
828 self._context.__exit__()
829 self._context = None
832def progressbar(label, maxval):
833 if force_dummy_progressbar:
834 return dummy_progressbar.ProgressBar(maxval=maxval).start()
836 return ProgressBar(label, maxval)
839def progress_beg(label):
840 '''
841 Notify user that an operation has started.
843 :param label: name of the operation
845 To be used in conjuction with :py:func:`progress_end`.
846 '''
848 sys.stderr.write(label)
849 sys.stderr.flush()
852def progress_end(label=''):
853 '''
854 Notify user that an operation has ended.
856 :param label: name of the operation
858 To be used in conjuction with :py:func:`progress_beg`.
859 '''
861 sys.stderr.write(' done. %s\n' % label)
862 sys.stderr.flush()
865class ArangeError(ValueError):
866 '''
867 Raised by :py:func:`arange2` for inconsistent range specifications.
868 '''
869 pass
872def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
873 '''
874 Return evenly spaced numbers over a specified interval.
876 Like :py:func:`numpy.arange` but returning floating point numbers by
877 default and with defined behaviour when stepsize is inconsistent with
878 interval bounds. It is considered inconsistent if the difference between
879 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
880 step``. Inconsistencies are handled according to the ``error`` parameter.
881 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
882 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
883 is silently changed to the closest, the next smaller, or next larger
884 multiple of ``step``, respectively.
885 '''
887 assert error in ('raise', 'round', 'floor', 'ceil')
889 start = dtype(start)
890 stop = dtype(stop)
891 step = dtype(step)
893 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
895 n = int(rnd((stop - start) / step)) + 1
896 stop_check = start + (n-1) * step
898 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
899 raise ArangeError(
900 'inconsistent range specification: start=%g, stop=%g, step=%g'
901 % (start, stop, step))
903 x = num.arange(n, dtype=dtype)
904 x *= step
905 x += start
906 return x
909def polylinefit(x, y, n_or_xnodes):
910 '''
911 Fit piece-wise linear function to data.
913 :param x,y: arrays with coordinates of data
914 :param n_or_xnodes: int, number of segments or x coordinates of polyline
916 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
917 polyline, root-mean-square error
918 '''
920 x = num.asarray(x)
921 y = num.asarray(y)
923 if isinstance(n_or_xnodes, int):
924 n = n_or_xnodes
925 xmin = x.min()
926 xmax = x.max()
927 xnodes = num.linspace(xmin, xmax, n+1)
928 else:
929 xnodes = num.asarray(n_or_xnodes)
930 n = xnodes.size - 1
932 assert len(x) == len(y)
933 assert n > 0
935 ndata = len(x)
936 a = num.zeros((ndata+(n-1), n*2))
937 for i in range(n):
938 xmin_block = xnodes[i]
939 xmax_block = xnodes[i+1]
940 if i == n-1: # don't loose last point
941 indices = num.where(
942 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
943 else:
944 indices = num.where(
945 num.logical_and(xmin_block <= x, x < xmax_block))[0]
947 a[indices, i*2] = x[indices]
948 a[indices, i*2+1] = 1.0
950 w = float(ndata)*100.
951 if i < n-1:
952 a[ndata+i, i*2] = xmax_block*w
953 a[ndata+i, i*2+1] = 1.0*w
954 a[ndata+i, i*2+2] = -xmax_block*w
955 a[ndata+i, i*2+3] = -1.0*w
957 d = num.concatenate((y, num.zeros(n-1)))
958 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
960 ynodes = num.zeros(n+1)
961 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
962 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
963 ynodes[1:n] *= 0.5
965 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
967 return xnodes, ynodes, rms_error
970def plf_integrate_piecewise(x_edges, x, y):
971 '''
972 Calculate definite integral of piece-wise linear function on intervals.
974 Use trapezoidal rule to calculate definite integral of a piece-wise linear
975 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
976 be sorted.
978 :param x_edges: array with edges of the intervals
979 :param x,y: arrays with coordinates of piece-wise linear function's
980 control points
981 '''
983 x_all = num.concatenate((x, x_edges))
984 ii = num.argsort(x_all)
985 y_edges = num.interp(x_edges, x, y)
986 y_all = num.concatenate((y, y_edges))
987 xs = x_all[ii]
988 ys = y_all[ii]
989 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
990 return num.diff(y_all[-len(y_edges):])
993def diff_fd_1d_4o(dt, data):
994 '''
995 Approximate first derivative of an array (forth order, central FD).
997 :param dt: sampling interval
998 :param data: NumPy array with data samples
1000 :returns: NumPy array with same shape as input
1002 Interior points are approximated to fourth order, edge points to first
1003 order right- or left-sided respectively, points next to edge to second
1004 order central.
1005 '''
1006 import scipy.signal
1008 ddata = num.empty_like(data, dtype=float)
1010 if data.size >= 5:
1011 ddata[2:-2] = scipy.signal.lfilter(
1012 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
1014 if data.size >= 3:
1015 ddata[1] = (data[2] - data[0]) / (2. * dt)
1016 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
1018 if data.size >= 2:
1019 ddata[0] = (data[1] - data[0]) / dt
1020 ddata[-1] = (data[-1] - data[-2]) / dt
1022 if data.size == 1:
1023 ddata[0] = 0.0
1025 return ddata
1028def diff_fd_1d_2o(dt, data):
1029 '''
1030 Approximate first derivative of an array (second order, central FD).
1032 :param dt: sampling interval
1033 :param data: NumPy array with data samples
1035 :returns: NumPy array with same shape as input
1037 Interior points are approximated to second order, edge points to first
1038 order right- or left-sided respectively.
1040 Uses :py:func:`numpy.gradient`.
1041 '''
1043 return num.gradient(data, dt)
1046def diff_fd_2d_4o(dt, data):
1047 '''
1048 Approximate second derivative of an array (forth order, central FD).
1050 :param dt: sampling interval
1051 :param data: NumPy array with data samples
1053 :returns: NumPy array with same shape as input
1055 Interior points are approximated to fourth order, next-to-edge points to
1056 second order, edge points repeated.
1057 '''
1058 import scipy.signal
1060 ddata = num.empty_like(data, dtype=float)
1062 if data.size >= 5:
1063 ddata[2:-2] = scipy.signal.lfilter(
1064 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
1066 if data.size >= 3:
1067 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
1068 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
1070 if data.size < 3:
1071 ddata[:] = 0.0
1073 return ddata
1076def diff_fd_2d_2o(dt, data):
1077 '''
1078 Approximate second derivative of an array (second order, central FD).
1080 :param dt: sampling interval
1081 :param data: NumPy array with data samples
1083 :returns: NumPy array with same shape as input
1085 Interior points are approximated to second order, edge points repeated.
1086 '''
1087 import scipy.signal
1089 ddata = num.empty_like(data, dtype=float)
1091 if data.size >= 3:
1092 ddata[1:-1] = scipy.signal.lfilter(
1093 [1., -2., 1.], [1.], data)[2:] / (dt**2)
1095 ddata[0] = ddata[1]
1096 ddata[-1] = ddata[-2]
1098 if data.size < 3:
1099 ddata[:] = 0.0
1101 return ddata
1104def diff_fd(n, order, dt, data):
1105 '''
1106 Approximate 1st or 2nd derivative of an array.
1108 :param n: 1 for first derivative, 2 for second
1109 :param order: order of the approximation 2 and 4 are supported
1110 :param dt: sampling interval
1111 :param data: NumPy array with data samples
1113 :returns: NumPy array with same shape as input
1115 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1116 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1117 :py:func:`diff_fd_2d_4o`.
1119 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1120 '''
1122 funcs = {
1123 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1124 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1126 try:
1127 funcs_n = funcs[n]
1128 except KeyError:
1129 raise ValueError(
1130 'pyrocko.util.diff_fd: '
1131 'Only 1st and 2sd derivatives are supported.')
1133 try:
1134 func = funcs_n[order]
1135 except KeyError:
1136 raise ValueError(
1137 'pyrocko.util.diff_fd: '
1138 'Order %i is not supported for %s derivative. Supported: %s' % (
1139 order, ['', '1st', '2nd'][n],
1140 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1142 return func(dt, data)
1145class GlobalVars(object):
1146 reuse_store = dict()
1147 decitab_nmax = 0
1148 decitab = {}
1149 decimate_fir_coeffs = {}
1150 decimate_fir_remez_coeffs = {}
1151 decimate_iir_coeffs = {}
1152 re_frac = None
1155def decimate_coeffs(q, n=None, ftype='iir'):
1157 q = int(q)
1159 if n is None:
1160 if ftype == 'fir':
1161 n = 30
1162 elif ftype == 'fir-remez':
1163 n = 45*q
1164 else:
1165 n = 8
1167 if ftype == 'fir':
1168 coeffs = GlobalVars.decimate_fir_coeffs
1169 if (n, 1./q) not in coeffs:
1170 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming')
1172 b = coeffs[n, 1./q]
1173 return b, [1.], n
1175 elif ftype == 'fir-remez':
1176 coeffs = GlobalVars.decimate_fir_remez_coeffs
1177 if (n, 1./q) not in coeffs:
1178 coeffs[n, 1./q] = signal.remez(
1179 n+1, (0., .75/q, 1./q, 1.),
1180 (1., 0.), fs=2, weight=(1, 50))
1181 b = coeffs[n, 1./q]
1182 return b, [1.], n
1184 else:
1185 coeffs = GlobalVars.decimate_iir_coeffs
1186 if (n, 0.05, 0.8/q) not in coeffs:
1187 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1189 b, a = coeffs[n, 0.05, 0.8/q]
1190 return b, a, n
1193def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1194 '''
1195 Downsample the signal x by an integer factor q, using an order n filter
1197 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1198 filter with hamming window if ftype is 'fir'.
1200 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`)
1201 :param q: the downsampling factor
1202 :param n: order of the filter (1 less than the length of the filter for a
1203 `fir` filter)
1204 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez`
1206 :returns: the downsampled signal (1D :class:`numpy.ndarray`)
1208 '''
1210 b, a, n = decimate_coeffs(q, n, ftype)
1212 if zi is None or zi is True:
1213 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1214 else:
1215 zi_ = zi
1217 y, zf = signal.lfilter(b, a, x, zi=zi_)
1219 if zi is not None:
1220 return y[n//2+ioff::q].copy(), zf
1221 else:
1222 return y[n//2+ioff::q].copy()
1225class UnavailableDecimation(Exception):
1226 '''
1227 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1228 '''
1230 pass
1233def gcd(a, b, epsilon=1e-7):
1234 '''
1235 Greatest common divisor.
1236 '''
1238 while b > epsilon*a:
1239 a, b = b, a % b
1241 return a
1244def lcm(a, b):
1245 '''
1246 Least common multiple.
1247 '''
1249 return a*b // gcd(a, b)
1252def mk_decitab(nmax=100):
1253 '''
1254 Make table with decimation sequences.
1256 Decimation from one sampling rate to a lower one is achieved by a
1257 successive application of :py:func:`decimate` with small integer
1258 downsampling factors (because using large downsampling factors can make the
1259 decimation unstable or slow). This function sets up a table with downsample
1260 sequences for factors up to ``nmax``.
1261 '''
1263 tab = GlobalVars.decitab
1264 for i in range(1, 10):
1265 for j in range(1, i+1):
1266 for k in range(1, j+1):
1267 for l_ in range(1, k+1):
1268 for m in range(1, l_+1):
1269 p = i*j*k*l_*m
1270 if p > nmax:
1271 break
1272 if p not in tab:
1273 tab[p] = (i, j, k, l_, m)
1274 if i*j*k*l_ > nmax:
1275 break
1276 if i*j*k > nmax:
1277 break
1278 if i*j > nmax:
1279 break
1280 if i > nmax:
1281 break
1283 GlobalVars.decitab_nmax = nmax
1286def zfmt(n):
1287 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1290def _year_to_time(year):
1291 tt = (year, 1, 1, 0, 0, 0)
1292 return to_time_float(calendar.timegm(tt))
1295def _working_year(year):
1296 try:
1297 tt = (year, 1, 1, 0, 0, 0)
1298 t = calendar.timegm(tt)
1299 tt2_ = time.gmtime(t)
1300 tt2 = tuple(tt2_)[:6]
1301 if tt != tt2:
1302 return False
1304 s = '%i-01-01 00:00:00' % year
1305 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1306 if s != s2:
1307 return False
1309 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1310 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1311 if s3 != s2:
1312 return False
1314 except Exception:
1315 return False
1317 return True
1320def working_system_time_range(year_min_lim=None, year_max_lim=None):
1321 '''
1322 Check time range supported by the systems's time conversion functions.
1324 Returns system time stamps of start of year of first/last fully supported
1325 year span. If this is before 1900 or after 2100, return first/last century
1326 which is fully supported.
1328 :returns: ``(tmin, tmax, year_min, year_max)``
1329 '''
1331 year0 = 2000
1332 year_min = year0
1333 year_max = year0
1335 itests = list(range(101))
1336 for i in range(19):
1337 itests.append(200 + i*100)
1339 for i in itests:
1340 year = year0 - i
1341 if year_min_lim is not None and year < year_min_lim:
1342 year_min = year_min_lim
1343 break
1344 elif not _working_year(year):
1345 break
1346 else:
1347 year_min = year
1349 for i in itests:
1350 year = year0 + i
1351 if year_max_lim is not None and year > year_max_lim:
1352 year_max = year_max_lim
1353 break
1354 elif not _working_year(year + 1):
1355 break
1356 else:
1357 year_max = year
1359 return (
1360 _year_to_time(year_min),
1361 _year_to_time(year_max),
1362 year_min, year_max)
1365g_working_system_time_range = None
1368def get_working_system_time_range():
1369 '''
1370 Caching variant of :py:func:`working_system_time_range`.
1371 '''
1373 global g_working_system_time_range
1374 if g_working_system_time_range is None:
1375 g_working_system_time_range = working_system_time_range()
1377 return g_working_system_time_range
1380def is_working_time(t):
1381 tmin, tmax, _, _ = get_working_system_time_range()
1382 return tmin <= t <= tmax
1385def julian_day_of_year(timestamp):
1386 '''
1387 Get the day number after the 1st of January of year in ``timestamp``.
1389 :returns: day number as int
1390 '''
1392 return time.gmtime(int(timestamp)).tm_yday
1395def hour_start(timestamp):
1396 '''
1397 Get beginning of hour for any point in time.
1399 :param timestamp: time instant as system timestamp (in seconds)
1401 :returns: instant of hour start as system timestamp
1402 '''
1404 tt = time.gmtime(int(timestamp))
1405 tts = tt[0:4] + (0, 0)
1406 return to_time_float(calendar.timegm(tts))
1409def day_start(timestamp):
1410 '''
1411 Get beginning of day for any point in time.
1413 :param timestamp: time instant as system timestamp (in seconds)
1415 :returns: instant of day start as system timestamp
1416 '''
1418 tt = time.gmtime(int(timestamp))
1419 tts = tt[0:3] + (0, 0, 0)
1420 return to_time_float(calendar.timegm(tts))
1423def month_start(timestamp):
1424 '''
1425 Get beginning of month for any point in time.
1427 :param timestamp: time instant as system timestamp (in seconds)
1429 :returns: instant of month start as system timestamp
1430 '''
1432 tt = time.gmtime(int(timestamp))
1433 tts = tt[0:2] + (1, 0, 0, 0)
1434 return to_time_float(calendar.timegm(tts))
1437def year_start(timestamp):
1438 '''
1439 Get beginning of year for any point in time.
1441 :param timestamp: time instant as system timestamp (in seconds)
1443 :returns: instant of year start as system timestamp
1444 '''
1446 tt = time.gmtime(int(timestamp))
1447 tts = tt[0:1] + (1, 1, 0, 0, 0)
1448 return to_time_float(calendar.timegm(tts))
1451def iter_days(tmin, tmax):
1452 '''
1453 Yields begin and end of days until given time span is covered.
1455 :param tmin,tmax: input time span
1457 :yields: tuples with (begin, end) of days as system timestamps
1458 '''
1460 t = day_start(tmin)
1461 while t < tmax:
1462 tend = day_start(t + 26*60*60)
1463 yield t, tend
1464 t = tend
1467def iter_months(tmin, tmax):
1468 '''
1469 Yields begin and end of months until given time span is covered.
1471 :param tmin,tmax: input time span
1473 :yields: tuples with (begin, end) of months as system timestamps
1474 '''
1476 t = month_start(tmin)
1477 while t < tmax:
1478 tend = month_start(t + 24*60*60*33)
1479 yield t, tend
1480 t = tend
1483def iter_years(tmin, tmax):
1484 '''
1485 Yields begin and end of years until given time span is covered.
1487 :param tmin,tmax: input time span
1489 :yields: tuples with (begin, end) of years as system timestamps
1490 '''
1492 t = year_start(tmin)
1493 while t < tmax:
1494 tend = year_start(t + 24*60*60*369)
1495 yield t, tend
1496 t = tend
1499def yesterday():
1500 return day_start(time.time() - 24*60*60)
1503def today():
1504 return day_start(time.time())
1507def tomorrow():
1508 return day_start(time.time() + 24*60*60)
1511def decitab(n):
1512 '''
1513 Get integer decimation sequence for given downampling factor.
1515 :param n: target decimation factor
1517 :returns: tuple with downsampling sequence
1518 '''
1520 if n > GlobalVars.decitab_nmax:
1521 mk_decitab(n*2)
1522 if n not in GlobalVars.decitab:
1523 raise UnavailableDecimation('ratio = %g' % n)
1524 return GlobalVars.decitab[n]
1527def ctimegm(s, format='%Y-%m-%d %H:%M:%S'):
1528 '''
1529 Convert string representing UTC time to system time.
1531 :param s: string to be interpreted
1532 :param format: format string passed to :py:func:`time.strptime`
1534 :returns: system time stamp
1536 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1538 .. note::
1539 This function is to be replaced by :py:func:`str_to_time`.
1540 '''
1542 return calendar.timegm(time.strptime(s, format))
1545def gmctime(t, format='%Y-%m-%d %H:%M:%S'):
1546 '''
1547 Get string representation from system time, UTC.
1549 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1551 .. note::
1552 This function is to be repaced by :py:func:`time_to_str`.
1553 '''
1555 return time.strftime(format, time.gmtime(t))
1558def gmctime_v(t, format='%a, %d %b %Y %H:%M:%S'):
1559 '''
1560 Get string representation from system time, UTC. Same as
1561 :py:func:`gmctime` but with a more verbose default format.
1563 .. note::
1564 This function is to be replaced by :py:func:`time_to_str`.
1565 '''
1567 return time.strftime(format, time.gmtime(t))
1570def gmctime_fn(t, format='%Y-%m-%d_%H-%M-%S'):
1571 '''
1572 Get string representation from system time, UTC. Same as
1573 :py:func:`gmctime` but with a default usable in filenames.
1575 .. note::
1576 This function is to be replaced by :py:func:`time_to_str`.
1577 '''
1579 return time.strftime(format, time.gmtime(t))
1582class TimeStrError(Exception):
1583 '''
1584 Raised for invalid time strings.
1585 '''
1586 pass
1589class FractionalSecondsMissing(TimeStrError):
1590 '''
1591 Exception raised by :py:func:`str_to_time` when the given string lacks
1592 fractional seconds.
1593 '''
1595 pass
1598class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1599 '''
1600 Exception raised by :py:func:`str_to_time` when the given string has an
1601 incorrect number of digits in the fractional seconds part.
1602 '''
1604 pass
1607def _endswith_n(s, endings):
1608 for ix, x in enumerate(endings):
1609 if s.endswith(x):
1610 return ix
1611 return -1
1614def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1615 '''
1616 Convert string representing UTC time to floating point system time.
1618 :param s: string representing UTC time
1619 :param format: time string format
1620 :returns: system time stamp as floating point value
1622 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1623 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1624 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1625 the fractional part, including the dot is made optional. The latter has the
1626 consequence, that the time strings and the format may not contain any other
1627 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1628 ensured, that exactly that number of digits are present in the fractional
1629 seconds.
1630 '''
1632 time_float = get_time_float()
1634 if util_ext is not None:
1635 try:
1636 t, tfrac = util_ext.stt(s, format)
1637 except util_ext.UtilExtError as e:
1638 raise TimeStrError(
1639 '%s, string=%s, format=%s' % (str(e), s, format))
1641 return time_float(t)+tfrac
1643 fracsec = 0.
1644 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1646 iend = _endswith_n(format, fixed_endings)
1647 if iend != -1:
1648 dotpos = s.rfind('.')
1649 if dotpos == -1:
1650 raise FractionalSecondsMissing(
1651 'string=%s, format=%s' % (s, format))
1653 if iend > 0 and iend != (len(s)-dotpos-1):
1654 raise FractionalSecondsWrongNumberOfDigits(
1655 'string=%s, format=%s' % (s, format))
1657 format = format[:-len(fixed_endings[iend])]
1658 fracsec = float(s[dotpos:])
1659 s = s[:dotpos]
1661 elif format.endswith('.OPTFRAC'):
1662 dotpos = s.rfind('.')
1663 format = format[:-8]
1664 if dotpos != -1 and len(s[dotpos:]) > 1:
1665 fracsec = float(s[dotpos:])
1667 if dotpos != -1:
1668 s = s[:dotpos]
1670 try:
1671 return time_float(calendar.timegm(time.strptime(s, format))) \
1672 + fracsec
1673 except ValueError as e:
1674 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1677stt = str_to_time
1680def str_to_time_fillup(s):
1681 '''
1682 Default :py:func:`str_to_time` with filling in of missing values.
1684 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`,
1685 `'2010-01-01 00'`, ..., or `'2010'`.
1686 '''
1687 s = s.strip().replace('T', ' ').replace('_', ' ')
1689 if s == 'now':
1690 return time.time()
1691 elif s == 'yesterday':
1692 return yesterday()
1693 elif s == 'today':
1694 return today()
1695 elif s == 'tomorrow':
1696 return tomorrow()
1698 if len(s) in (4, 7, 10, 13, 16):
1699 s += '0000-01-01 00:00:00'[len(s):]
1700 return str_to_time(s)
1703def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1704 '''
1705 Get string representation for floating point system time.
1707 :param t: floating point system time
1708 :param format: time string format
1709 :returns: string representing UTC time
1711 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1712 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1713 digit between 1 and 9, this is replaced with the fractional part of ``t``
1714 with ``x`` digits precision.
1715 '''
1717 if pyrocko.grumpy > 0:
1718 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1720 if isinstance(format, int):
1721 if format > 0:
1722 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1723 else:
1724 format = '%Y-%m-%d %H:%M:%S'
1726 if util_ext is not None:
1727 t0 = num.floor(t)
1728 try:
1729 return util_ext.tts(int(t0), float(t - t0), format)
1730 except util_ext.UtilExtError as e:
1731 raise TimeStrError(
1732 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1734 if not GlobalVars.re_frac:
1735 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1736 GlobalVars.frac_formats = dict(
1737 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1739 ts = float(num.floor(t))
1740 tfrac = t-ts
1742 m = GlobalVars.re_frac.search(format)
1743 if m:
1744 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1745 if sfrac[0] == '1':
1746 ts += 1.
1748 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1750 return time.strftime(format, time.gmtime(ts))
1753tts = time_to_str
1754_abbr_weekday = 'Mon Tue Wed Thu Fri Sat Sun'.split()
1755_abbr_month = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split()
1758def mystrftime(fmt=None, tt=None, milliseconds=0):
1759 # Needed by Snuffler for the time axis. In other cases `time_to_str`
1760 # should be used.
1762 if fmt is None:
1763 fmt = '%Y-%m-%d %H:%M:%S .%r'
1765 # Get these two locale independent, needed by Snuffler.
1766 # Setting the locale seems to have no effect.
1767 fmt = fmt.replace('%a', _abbr_weekday[tt.tm_wday])
1768 fmt = fmt.replace('%b', _abbr_month[tt.tm_mon-1])
1770 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1771 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1772 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1773 return time.strftime(fmt, tt)
1776def gmtime_x(timestamp):
1777 etimestamp = float(num.floor(timestamp))
1778 tt = time.gmtime(etimestamp)
1779 ms = (timestamp-etimestamp)*1000
1780 return tt, ms
1783def plural_s(n):
1784 if not isinstance(n, int):
1785 n = len(n)
1787 return 's' if n != 1 else ''
1790def ensuredirs(dst):
1791 '''
1792 Create all intermediate path components for a target path.
1794 :param dst: target path
1796 The leaf part of the target path is not created (use :py:func:`ensuredir`
1797 if a the target path is a directory to be created).
1798 '''
1800 d, x = os.path.split(dst.rstrip(os.sep))
1801 dirs = []
1802 while d and not os.path.exists(d):
1803 dirs.append(d)
1804 d, x = os.path.split(d)
1806 dirs.reverse()
1808 for d in dirs:
1809 try:
1810 os.mkdir(d)
1811 except OSError as e:
1812 if not e.errno == errno.EEXIST:
1813 raise
1816def ensuredir(dst):
1817 '''
1818 Create directory and all intermediate path components to it as needed.
1820 :param dst: directory name
1822 Nothing is done if the given target already exists.
1823 '''
1825 if os.path.exists(dst):
1826 return
1828 dst.rstrip(os.sep)
1830 ensuredirs(dst)
1831 try:
1832 os.mkdir(dst)
1833 except OSError as e:
1834 if not e.errno == errno.EEXIST:
1835 raise
1838def reuse(x):
1839 '''
1840 Get unique instance of an object.
1842 :param x: hashable object
1843 :returns: reference to x or an equivalent object
1845 Cache object ``x`` in a global dict for reuse, or if x already
1846 is in that dict, return a reference to it.
1847 '''
1849 grs = GlobalVars.reuse_store
1850 if x not in grs:
1851 grs[x] = x
1852 return grs[x]
1855def deuse(x):
1856 grs = GlobalVars.reuse_store
1857 if x in grs:
1858 del grs[x]
1861class Anon(object):
1862 '''
1863 Dict-to-object utility.
1865 Any given arguments are stored as attributes.
1867 Example::
1869 a = Anon(x=1, y=2)
1870 print a.x, a.y
1871 '''
1873 def __init__(self, **dict):
1874 for k in dict:
1875 self.__dict__[k] = dict[k]
1878def iter_select_files(
1879 paths,
1880 include=None,
1881 exclude=None,
1882 selector=None,
1883 show_progress=True,
1884 pass_through=None):
1886 '''
1887 Recursively select files (generator variant).
1889 See :py:func:`select_files`.
1890 '''
1892 if show_progress:
1893 progress_beg('selecting files...')
1895 ngood = 0
1896 if include is not None:
1897 rinclude = re.compile(include)
1899 def check_include(path):
1900 m = rinclude.search(path)
1901 if not m:
1902 return False
1904 if selector is None:
1905 return True
1907 infos = Anon(**m.groupdict())
1908 return selector(infos)
1909 else:
1910 check_include = None
1912 if exclude is not None:
1913 rexclude = re.compile(exclude)
1915 def check_exclude(path):
1916 return not bool(rexclude.search(path))
1917 else:
1918 check_exclude = None
1920 if check_include and check_exclude:
1922 def check(path):
1923 return check_include(path) and check_exclude(path)
1925 elif check_include:
1926 check = check_include
1928 elif check_exclude:
1929 check = check_exclude
1931 else:
1932 check = None
1934 if isinstance(paths, str):
1935 paths = [paths]
1937 for path in paths:
1938 if pass_through and pass_through(path):
1939 if check is None or check(path):
1940 yield path
1942 elif os.path.isdir(path):
1943 for (dirpath, dirnames, filenames) in os.walk(path):
1944 dirnames.sort()
1945 filenames.sort()
1946 for filename in filenames:
1947 path = op.join(dirpath, filename)
1948 if check is None or check(path):
1949 yield os.path.abspath(path)
1950 ngood += 1
1951 else:
1952 if check is None or check(path):
1953 yield os.path.abspath(path)
1954 ngood += 1
1956 if show_progress:
1957 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1960def select_files(
1961 paths, include=None, exclude=None, selector=None, show_progress=True,
1962 regex=None):
1964 '''
1965 Recursively select files.
1967 :param paths: entry path names
1968 :param include: pattern for conditional inclusion
1969 :param exclude: pattern for conditional exclusion
1970 :param selector: callback for conditional inclusion
1971 :param show_progress: if True, indicate start and stop of processing
1972 :param regex: alias for ``include`` (backwards compatibility)
1973 :returns: list of path names
1975 Recursively finds all files under given entry points ``paths``. If
1976 parameter ``include`` is a regular expression, only files with matching
1977 path names are included. If additionally parameter ``selector`` is given a
1978 callback function, only files for which the callback returns ``True`` are
1979 included. The callback should take a single argument. The callback is
1980 called with a single argument, an object, having as attributes, any named
1981 groups given in ``include``.
1983 Examples
1985 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1987 select_files(paths,
1988 include=r'\\.(mseed|msd)$')
1990 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1991 the year::
1993 select_files(paths,
1994 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1995 selector=(lambda x: int(x.year) == 2009))
1996 '''
1997 if regex is not None:
1998 assert include is None
1999 include = regex
2001 return list(iter_select_files(
2002 paths, include, exclude, selector, show_progress))
2005def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
2006 '''
2007 Convert positive integer to a base36 string.
2008 '''
2010 if not isinstance(number, (int, long)):
2011 raise TypeError('number must be an integer')
2012 if number < 0:
2013 raise ValueError('number must be positive')
2015 # Special case for small numbers
2016 if number < 36:
2017 return alphabet[number]
2019 base36 = ''
2020 while number != 0:
2021 number, i = divmod(number, 36)
2022 base36 = alphabet[i] + base36
2024 return base36
2027def base36decode(number):
2028 '''
2029 Decode base36 endcoded positive integer.
2030 '''
2032 return int(number, 36)
2035class UnpackError(Exception):
2036 '''
2037 Exception raised when :py:func:`unpack_fixed` encounters an error.
2038 '''
2040 pass
2043ruler = ''.join(['%-10i' % i for i in range(8)]) \
2044 + '\n' + '0123456789' * 8 + '\n'
2047def unpack_fixed(format, line, *callargs):
2048 '''
2049 Unpack fixed format string, as produced by many fortran codes.
2051 :param format: format specification
2052 :param line: string to be processed
2053 :param callargs: callbacks for callback fields in the format
2055 The format is described by a string of comma-separated fields. Each field
2056 is defined by a character for the field type followed by the field width. A
2057 questionmark may be appended to the field description to allow the argument
2058 to be optional (The data string is then allowed to be filled with blanks
2059 and ``None`` is returned in this case).
2061 The following field types are available:
2063 ==== ================================================================
2064 Type Description
2065 ==== ================================================================
2066 A string (full field width is extracted)
2067 a string (whitespace at the beginning and the end is removed)
2068 i integer value
2069 f floating point value
2070 @ special type, a callback must be given for the conversion
2071 x special field type to skip parts of the string
2072 ==== ================================================================
2073 '''
2075 ipos = 0
2076 values = []
2077 icall = 0
2078 for form in format.split(','):
2079 form = form.strip()
2080 optional = form[-1] == '?'
2081 form = form.rstrip('?')
2082 typ = form[0]
2083 ln = int(form[1:])
2084 s = line[ipos:ipos+ln]
2085 cast = {
2086 'x': None,
2087 'A': str,
2088 'a': lambda x: x.strip(),
2089 'i': int,
2090 'f': float,
2091 '@': 'extra'}[typ]
2093 if cast == 'extra':
2094 cast = callargs[icall]
2095 icall += 1
2097 if cast is not None:
2098 if optional and s.strip() == '':
2099 values.append(None)
2100 else:
2101 try:
2102 values.append(cast(s))
2103 except Exception:
2104 mark = [' '] * 80
2105 mark[ipos:ipos+ln] = ['^'] * ln
2106 mark = ''.join(mark)
2107 raise UnpackError(
2108 'Invalid cast to type "%s" at position [%i:%i] of '
2109 'line: \n%s%s\n%s' % (
2110 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
2112 ipos += ln
2114 return values
2117_pattern_cache = {}
2120def _nslc_pattern(pattern):
2121 if pattern not in _pattern_cache:
2122 rpattern = re.compile(fnmatch.translate(pattern), re.I)
2123 _pattern_cache[pattern] = rpattern
2124 else:
2125 rpattern = _pattern_cache[pattern]
2127 return rpattern
2130def match_nslc(patterns, nslc):
2131 '''
2132 Match network-station-location-channel code against pattern or list of
2133 patterns.
2135 :param patterns: pattern or list of patterns
2136 :param nslc: tuple with (network, station, location, channel) as strings
2138 :returns: ``True`` if the pattern matches or if any of the given patterns
2139 match; or ``False``.
2141 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2143 Example::
2145 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2146 '''
2148 if isinstance(patterns, str):
2149 patterns = [patterns]
2151 if not isinstance(nslc, str):
2152 s = '.'.join(nslc)
2153 else:
2154 s = nslc
2156 for pattern in patterns:
2157 if _nslc_pattern(pattern).match(s):
2158 return True
2160 return False
2163def match_nslcs(patterns, nslcs):
2164 '''
2165 Get network-station-location-channel codes that match given pattern or any
2166 of several given patterns.
2168 :param patterns: pattern or list of patterns
2169 :param nslcs: list of (network, station, location, channel) tuples
2171 See also :py:func:`match_nslc`
2172 '''
2174 matching = []
2175 for nslc in nslcs:
2176 if match_nslc(patterns, nslc):
2177 matching.append(nslc)
2179 return matching
2182def glob_filter(patterns, names):
2183 '''
2184 Select names matching any of the given patterns.
2185 '''
2187 if not patterns:
2188 return names
2190 rpattern = re.compile(r'|'.join(
2191 fnmatch.translate(pattern) for pattern in patterns))
2193 return [name for name in names if rpattern.match(name)]
2196class Timeout(Exception):
2197 pass
2200def create_lockfile(fn, timeout=None, timewarn=10.):
2201 t0 = time.time()
2203 while True:
2204 try:
2205 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
2206 os.close(f)
2207 return
2209 except OSError as e:
2210 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows
2211 pass # retry
2212 else:
2213 raise
2215 tnow = time.time()
2217 if timeout is not None and tnow - t0 > timeout:
2218 raise Timeout(
2219 'Timeout (%gs) occured while waiting to get exclusive '
2220 'access to: %s' % (timeout, fn))
2222 if timewarn is not None and tnow - t0 > timewarn:
2223 logger.warning(
2224 'Waiting since %gs to get exclusive access to: %s' % (
2225 timewarn, fn))
2227 timewarn *= 2
2229 time.sleep(0.01)
2232def delete_lockfile(fn):
2233 os.unlink(fn)
2236class Lockfile(Exception):
2238 def __init__(self, path, timeout=5, timewarn=10.):
2239 self._path = path
2240 self._timeout = timeout
2241 self._timewarn = timewarn
2243 def __enter__(self):
2244 create_lockfile(
2245 self._path, timeout=self._timeout, timewarn=self._timewarn)
2246 return None
2248 def __exit__(self, type, value, traceback):
2249 delete_lockfile(self._path)
2252class SoleError(Exception):
2253 '''
2254 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2255 instance is running.
2256 '''
2258 pass
2261class Sole(object):
2262 '''
2263 Use POSIX advisory file locking to ensure that only a single instance of a
2264 program is running.
2266 :param pid_path: path to lockfile to be used
2268 Usage::
2270 from pyrocko.util import Sole, SoleError, setup_logging
2271 import os
2273 setup_logging('my_program')
2275 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2276 try:
2277 sole = Sole(pid_path)
2279 except SoleError, e:
2280 logger.fatal( str(e) )
2281 sys.exit(1)
2282 '''
2284 def __init__(self, pid_path):
2285 self._pid_path = pid_path
2286 self._other_running = False
2287 ensuredirs(self._pid_path)
2288 self._lockfile = None
2289 self._os = os
2291 if not fcntl:
2292 raise SoleError(
2293 'Python standard library module "fcntl" not available on '
2294 'this platform.')
2296 self._fcntl = fcntl
2298 try:
2299 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2300 except Exception:
2301 raise SoleError(
2302 'Cannot open lockfile (path = %s)' % self._pid_path)
2304 try:
2305 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2307 except IOError:
2308 self._other_running = True
2309 try:
2310 f = open(self._pid_path, 'r')
2311 pid = f.read().strip()
2312 f.close()
2313 except Exception:
2314 pid = '?'
2316 raise SoleError('Other instance is running (pid = %s)' % pid)
2318 try:
2319 os.ftruncate(self._lockfile, 0)
2320 os.write(self._lockfile, '%i\n' % os.getpid())
2321 os.fsync(self._lockfile)
2323 except Exception:
2324 # the pid is only stored for user information, so this is allowed
2325 # to fail
2326 pass
2328 def __del__(self):
2329 if not self._other_running:
2330 if self._lockfile is not None:
2331 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2332 self._os.close(self._lockfile)
2333 try:
2334 self._os.unlink(self._pid_path)
2335 except Exception:
2336 pass
2339re_escapequotes = re.compile(r"(['\\])")
2342def escapequotes(s):
2343 return re_escapequotes.sub(r'\\\1', s)
2346class TableWriter(object):
2347 '''
2348 Write table of space separated values to a file.
2350 :param f: file like object
2352 Strings containing spaces are quoted on output.
2353 '''
2355 def __init__(self, f):
2356 self._f = f
2358 def writerow(self, row, minfieldwidths=None):
2360 '''
2361 Write one row of values to underlying file.
2363 :param row: iterable of values
2364 :param minfieldwidths: minimum field widths for the values
2366 Each value in in ``row`` is converted to a string and optionally padded
2367 with blanks. The resulting strings are output separated with blanks. If
2368 any values given are strings and if they contain whitespace, they are
2369 quoted with single quotes, and any internal single quotes are
2370 backslash-escaped.
2371 '''
2373 out = []
2375 for i, x in enumerate(row):
2376 w = 0
2377 if minfieldwidths and i < len(minfieldwidths):
2378 w = minfieldwidths[i]
2380 if isinstance(x, str):
2381 if re.search(r"\s|'", x):
2382 x = "'%s'" % escapequotes(x)
2384 x = x.ljust(w)
2385 else:
2386 x = str(x).rjust(w)
2388 out.append(x)
2390 self._f.write(' '.join(out).rstrip() + '\n')
2393class TableReader(object):
2395 '''
2396 Read table of space separated values from a file.
2398 :param f: file-like object
2400 This uses Pythons shlex module to tokenize lines. Should deal correctly
2401 with quoted strings.
2402 '''
2404 def __init__(self, f):
2405 self._f = f
2406 self.eof = False
2408 def readrow(self):
2409 '''
2410 Read one row from the underlying file, tokenize it with shlex.
2412 :returns: tokenized line as a list of strings.
2413 '''
2415 line = self._f.readline()
2416 if not line:
2417 self.eof = True
2418 return []
2419 line.strip()
2420 if line.startswith('#'):
2421 return []
2423 return qsplit(line)
2426def gform(number, significant_digits=3):
2427 '''
2428 Pretty print floating point numbers.
2430 Align floating point numbers at the decimal dot.
2432 ::
2434 | -d.dde+xxx|
2435 | -d.dde+xx |
2436 |-ddd. |
2437 | -dd.d |
2438 | -d.dd |
2439 | -0.ddd |
2440 | -0.0ddd |
2441 | -0.00ddd |
2442 | -d.dde-xx |
2443 | -d.dde-xxx|
2444 | nan|
2447 The formatted string has length ``significant_digits * 2 + 6``.
2448 '''
2450 no_exp_range = (pow(10., -1),
2451 pow(10., significant_digits))
2452 width = significant_digits+significant_digits-1+1+1+5
2454 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2455 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2456 else:
2457 s = '%.*E' % (significant_digits-1, number)
2458 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2459 if s.strip().lower() == 'nan':
2460 s = 'nan'.rjust(width)
2461 return s
2464def human_bytesize(value):
2466 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2468 if value == 1:
2469 return '1 Byte'
2471 for i, ext in enumerate(exts):
2472 x = float(value) / 1000**i
2473 if round(x) < 10. and not value < 1000:
2474 return '%.1f %s' % (x, ext)
2475 if round(x) < 1000.:
2476 return '%.0f %s' % (x, ext)
2478 return '%i Bytes' % value
2481def human_intsize(value):
2483 exts = 'x k M G T P E Z Y'.split()
2485 for i, ext in enumerate(exts):
2486 x = float(value) / 1000**i
2487 if round(x) < 10. and not value < 1000:
2488 return '%.1f%s' % (x, ' ' + ext if ext != 'x' else '')
2489 if round(x) < 1000.:
2490 return '%.0f%s' % (x, ' ' + ext if ext != 'x' else '')
2492 return '%i' % value
2495re_compatibility = re.compile(
2496 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2497 r'(dummy|poel|qseis|qssp))\.'
2498)
2501def pf_is_old(fn):
2502 oldstyle = False
2503 with open(fn, 'r') as f:
2504 for line in f:
2505 if re_compatibility.search(line):
2506 oldstyle = True
2508 return oldstyle
2511def pf_upgrade(fn):
2512 need = pf_is_old(fn)
2513 if need:
2514 fn_temp = fn + '.temp'
2516 with open(fn, 'r') as fin:
2517 with open(fn_temp, 'w') as fout:
2518 for line in fin:
2519 line = re_compatibility.sub('!pf.', line)
2520 fout.write(line)
2522 os.rename(fn_temp, fn)
2524 return need
2527def get_tzdata_path():
2528 for path in [
2529 '/usr/share/zoneinfo/right/UTC',
2530 '/usr/share/zoneinfo/Etc/UTC']:
2532 if os.path.exists(path):
2533 return path
2535 raise LeapSecondsError(
2536 'Could not find system-provided leap time zone data file.')
2539def read_leap_seconds(tzfile=None):
2540 '''
2541 Extract leap second information from tzdata.
2543 Based on example at http://stackoverflow.com/questions/19332902/\
2544 extract-historic-leap-seconds-from-tzdata
2546 See also 'man 5 tzfile'.
2547 '''
2549 if tzfile is None:
2550 tzfile = get_tzdata_path()
2552 from struct import unpack, calcsize
2553 out = []
2554 with open(tzfile, 'rb') as f:
2555 # read header
2556 fmt = '>4s c 15x 6l'
2557 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2558 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2559 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2561 # skip over some uninteresting data
2562 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2563 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2564 f.read(calcsize(fmt))
2566 # read leap-seconds
2567 fmt = '>2l'
2568 for i in range(leapcnt):
2569 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2570 out.append((tleap-nleap+1, nleap))
2572 return out
2575class LeapSecondsError(Exception):
2576 pass
2579class LeapSecondsOutdated(LeapSecondsError):
2580 pass
2583class InvalidLeapSecondsFile(LeapSecondsOutdated):
2584 pass
2587def parse_leap_seconds_list(fn):
2588 data = []
2589 texpires = None
2590 try:
2591 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2592 except TimeStrError:
2593 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2595 tnow = int(round(time.time()))
2597 if not op.exists(fn):
2598 raise LeapSecondsOutdated('no leap seconds file found')
2600 try:
2601 with open(fn, 'rb') as f:
2602 for line in f:
2603 if line.strip().startswith(b'<!DOCTYPE'):
2604 raise InvalidLeapSecondsFile('invalid leap seconds file')
2606 if line.startswith(b'#@'):
2607 texpires = int(line.split()[1]) + t0
2608 elif line.startswith(b'#') or len(line) < 5:
2609 pass
2610 else:
2611 toks = line.split()
2612 t = int(toks[0]) + t0
2613 nleap = int(toks[1]) - 10
2614 data.append((t, nleap))
2616 except IOError:
2617 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2619 if texpires is None or tnow > texpires:
2620 raise LeapSecondsOutdated('leap seconds list is outdated')
2622 return data
2625def read_leap_seconds2():
2626 from pyrocko import config
2627 conf = config.config()
2628 fn = conf.leapseconds_path
2629 url = conf.leapseconds_url
2630 # check for outdated default URL
2631 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2632 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2633 logger.info(
2634 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2636 if url == 'https://www.ietf.org/timezones/data/leap-seconds.list':
2637 url = 'https://hpiers.obspm.fr/iers/bul/bulc/ntp/leap-seconds.list'
2638 logger.info(
2639 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2641 for i in range(3):
2642 try:
2643 return parse_leap_seconds_list(fn)
2645 except LeapSecondsOutdated:
2646 try:
2647 logger.info('updating leap seconds list...')
2648 download_file(url, fn)
2650 except Exception as e:
2651 raise LeapSecondsError(
2652 'cannot download leap seconds list from %s to %s (%s)'
2653 % (url, fn, e))
2655 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2658def gps_utc_offset(t_utc):
2659 '''
2660 Time offset t_gps - t_utc for a given t_utc.
2661 '''
2662 ls = read_leap_seconds2()
2663 i = 0
2664 if t_utc < ls[0][0]:
2665 return ls[0][1] - 1 - 9
2667 while i < len(ls) - 1:
2668 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2669 return ls[i][1] - 9
2670 i += 1
2672 return ls[-1][1] - 9
2675def utc_gps_offset(t_gps):
2676 '''
2677 Time offset t_utc - t_gps for a given t_gps.
2678 '''
2679 ls = read_leap_seconds2()
2681 if t_gps < ls[0][0] + ls[0][1] - 9:
2682 return - (ls[0][1] - 1 - 9)
2684 i = 0
2685 while i < len(ls) - 1:
2686 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2687 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2688 return - (ls[i][1] - 9)
2689 i += 1
2691 return - (ls[-1][1] - 9)
2694def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2695 import itertools
2696 import glob
2697 from pyrocko.io.io_common import FileLoadError
2699 def iload_filename(filename, **kwargs):
2700 try:
2701 with open(filename, 'rb') as f:
2702 for cr in iload_fh(f, **kwargs):
2703 yield cr
2705 except FileLoadError as e:
2706 e.set_context('filename', filename)
2707 raise
2709 def iload_dirname(dirname, **kwargs):
2710 for entry in os.listdir(dirname):
2711 fpath = op.join(dirname, entry)
2712 if op.isfile(fpath):
2713 for cr in iload_filename(fpath, **kwargs):
2714 yield cr
2716 def iload_glob(pattern, **kwargs):
2718 for fn in glob.iglob(pattern):
2719 for cr in iload_filename(fn, **kwargs):
2720 yield cr
2722 def iload(source, **kwargs):
2723 if isinstance(source, str):
2724 if op.isdir(source):
2725 return iload_dirname(source, **kwargs)
2726 elif op.isfile(source):
2727 return iload_filename(source, **kwargs)
2728 else:
2729 return iload_glob(source, **kwargs)
2731 elif hasattr(source, 'read'):
2732 return iload_fh(source, **kwargs)
2733 else:
2734 return itertools.chain.from_iterable(
2735 iload(subsource, **kwargs) for subsource in source)
2737 iload_filename.__doc__ = '''
2738 Read %s information from named file.
2739 ''' % doc_fmt
2741 iload_dirname.__doc__ = '''
2742 Read %s information from directory of %s files.
2743 ''' % (doc_fmt, doc_fmt)
2745 iload_glob.__doc__ = '''
2746 Read %s information from files matching a glob pattern.
2747 ''' % doc_fmt
2749 iload.__doc__ = '''
2750 Load %s information from given source(s)
2752 The ``source`` can be specified as the name of a %s file, the name of a
2753 directory containing %s files, a glob pattern of %s files, an open
2754 filehandle or an iterator yielding any of the forementioned sources.
2756 This function behaves as a generator yielding %s objects.
2757 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2759 for f in iload_filename, iload_dirname, iload_glob, iload:
2760 f.__module__ = iload_fh.__module__
2762 return iload_filename, iload_dirname, iload_glob, iload
2765class Inconsistency(Exception):
2766 pass
2769def consistency_check(list_of_tuples, message='values differ:'):
2770 '''
2771 Check for inconsistencies.
2773 Given a list of tuples, check that all tuple elements except for first one
2774 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2775 valid because the coordinates at the two channels are the same.
2776 '''
2778 if len(list_of_tuples) >= 2:
2779 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2780 raise Inconsistency('%s\n' % message + '\n'.join(
2781 ' %s: %s' % (t[0], ', '.join(str(x) for x in t[1:]))
2782 for t in list_of_tuples))
2785class defaultzerodict(dict):
2786 def __missing__(self, k):
2787 return 0
2790def mostfrequent(x):
2791 c = defaultzerodict()
2792 for e in x:
2793 c[e] += 1
2795 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2798def consistency_merge(list_of_tuples,
2799 message='values differ:',
2800 error='raise',
2801 merge=mostfrequent):
2803 assert error in ('raise', 'warn', 'ignore')
2805 if len(list_of_tuples) == 0:
2806 raise Exception('cannot merge empty sequence')
2808 try:
2809 consistency_check(list_of_tuples, message)
2810 return list_of_tuples[0][1:]
2811 except Inconsistency as e:
2812 if error == 'raise':
2813 raise
2815 elif error == 'warn':
2816 logger.warning(str(e))
2818 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2821def short_to_list(nmax, it):
2822 import itertools
2824 if isinstance(it, list):
2825 return it
2827 li = []
2828 for i in range(nmax+1):
2829 try:
2830 li.append(next(it))
2831 except StopIteration:
2832 return li
2834 return itertools.chain(li, it)
2837def parse_md(f):
2838 try:
2839 with open(op.join(
2840 op.dirname(op.abspath(f)),
2841 'README.md'), 'r') as readme:
2842 mdstr = readme.read()
2843 except IOError as e:
2844 return 'Failed to get README.md: %s' % e
2846 # Remve the title
2847 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2848 # Append sphinx reference to `pyrocko.` modules
2849 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2850 # Convert Subsections to toc-less rubrics
2851 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2852 return mdstr
2855def mpl_show(plt):
2856 import matplotlib
2857 if matplotlib.get_backend().lower() == 'agg':
2858 logger.warning('Cannot show() when using matplotlib "agg" backend')
2859 else:
2860 plt.show()
2863g_re_qsplit = re.compile(
2864 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2865g_re_qsplit_sep = {}
2868def get_re_qsplit(sep):
2869 if sep is None:
2870 return g_re_qsplit
2871 else:
2872 if sep not in g_re_qsplit_sep:
2873 assert len(sep) == 1
2874 assert sep not in '\'"'
2875 esep = re.escape(sep)
2876 g_re_qsplit_sep[sep] = re.compile(
2877 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|'
2878 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa
2879 return g_re_qsplit_sep[sep]
2882g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2883g_re_trivial_sep = {}
2886def get_re_trivial(sep):
2887 if sep is None:
2888 return g_re_trivial
2889 else:
2890 if sep not in g_re_qsplit_sep:
2891 assert len(sep) == 1
2892 assert sep not in '\'"'
2893 esep = re.escape(sep)
2894 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z')
2896 return g_re_trivial_sep[sep]
2899g_re_escape_s = re.compile(r'([\\\'])')
2900g_re_unescape_s = re.compile(r'\\([\\\'])')
2901g_re_escape_d = re.compile(r'([\\"])')
2902g_re_unescape_d = re.compile(r'\\([\\"])')
2905def escape_s(s):
2906 '''
2907 Backslash-escape single-quotes and backslashes.
2909 Example: ``Jack's`` => ``Jack\\'s``
2911 '''
2912 return g_re_escape_s.sub(r'\\\1', s)
2915def unescape_s(s):
2916 '''
2917 Unescape backslash-escaped single-quotes and backslashes.
2919 Example: ``Jack\\'s`` => ``Jack's``
2920 '''
2921 return g_re_unescape_s.sub(r'\1', s)
2924def escape_d(s):
2925 '''
2926 Backslash-escape double-quotes and backslashes.
2928 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2929 '''
2930 return g_re_escape_d.sub(r'\\\1', s)
2933def unescape_d(s):
2934 '''
2935 Unescape backslash-escaped double-quotes and backslashes.
2937 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2938 '''
2939 return g_re_unescape_d.sub(r'\1', s)
2942def qjoin_s(it, sep=None):
2943 '''
2944 Join sequence of strings into a line, single-quoting non-trivial strings.
2946 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2947 '''
2948 re_trivial = get_re_trivial(sep)
2950 if sep is None:
2951 sep = ' '
2953 return sep.join(
2954 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2957def qjoin_d(it, sep=None):
2958 '''
2959 Join sequence of strings into a line, double-quoting non-trivial strings.
2961 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2962 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2963 '''
2964 re_trivial = get_re_trivial(sep)
2965 if sep is None:
2966 sep = ' '
2968 return sep.join(
2969 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2972def qsplit(s, sep=None):
2973 '''
2974 Split line into list of strings, allowing for quoted strings.
2976 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2977 ``["55", "Sparrow's Island"]``,
2978 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2979 ``['55', 'Pete "The Robot" Smith']``
2980 '''
2981 re_qsplit = get_re_qsplit(sep)
2982 return [
2983 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2984 for x in re_qsplit.findall(s)]
2987g_have_warned_threadpoolctl = False
2990class threadpool_limits_dummy(object):
2992 def __init__(self, *args, **kwargs):
2993 pass
2995 def __enter__(self):
2996 global g_have_warned_threadpoolctl
2998 if not g_have_warned_threadpoolctl:
2999 logger.warning(
3000 'Cannot control number of BLAS threads because '
3001 '`threadpoolctl` module is not available. You may want to '
3002 'install `threadpoolctl`.')
3004 g_have_warned_threadpoolctl = True
3006 return self
3008 def __exit__(self, type, value, traceback):
3009 pass
3012def get_threadpool_limits():
3013 '''
3014 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
3015 '''
3017 try:
3018 from threadpoolctl import threadpool_limits
3019 return threadpool_limits
3021 except ImportError:
3022 return threadpool_limits_dummy
3025def fmt_summary(entries, widths):
3026 return ' | '.join(
3027 entry.ljust(width) for (entry, width) in zip(entries, widths))
3030def smart_weakref(obj, callback=None):
3031 if inspect.ismethod(obj):
3032 return weakref.WeakMethod(obj, callback)
3033 else:
3034 return weakref.ref(obj, callback)
3037DEFAULT_SIGNALS = []
3038for signame in ['SIGINT', 'SIGTERM', 'SIGHUP']:
3039 try:
3040 DEFAULT_SIGNALS.append(getattr(psignal, signame))
3041 except AttributeError:
3042 pass
3045class SignalQuitable:
3047 def __init__(self, signals=DEFAULT_SIGNALS):
3049 self.quit_requested = False
3050 self._original = {}
3051 self._signals = signals
3053 def quit(self, *_):
3054 self.quit_requested = True
3056 def __enter__(self):
3057 for sig in self._signals:
3058 self._original[sig] = psignal.signal(sig, self.quit)
3060 return self
3062 def __exit__(self, *_):
3063 for sig, handler in self._original.items():
3064 psignal.signal(sig, handler)
3066 self._original = {}
3069def group_by(key, xs):
3070 groups = defaultdict(list)
3071 for x in xs:
3072 groups[key(x)].append(x)
3074 return groups