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 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.
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.
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.
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:
36.. code-block :: python
38 import time
39 from pyrocko import util
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
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())
49 # To get the appropriate float class, use:
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]
56 # Shortcut:
57 util.check_time_class(t1)
59Module content
60..............
61'''
63from __future__ import division, print_function
65import time
66import logging
67import os
68import sys
69import re
70import calendar
71import math
72import fnmatch
73try:
74 import fcntl
75except ImportError:
76 fcntl = None
77import optparse
78import os.path as op
79import errno
81import numpy as num
82from scipy import signal
83import pyrocko
84from pyrocko import dummy_progressbar
86try:
87 from urllib.parse import urlencode, quote, unquote # noqa
88 from urllib.request import (
89 Request, build_opener, HTTPDigestAuthHandler, urlopen as _urlopen) # noqa
90 from urllib.error import HTTPError, URLError # noqa
92except ImportError:
93 from urllib import urlencode, quote, unquote # noqa
94 from urllib2 import (Request, build_opener, HTTPDigestAuthHandler, # noqa
95 HTTPError, URLError, urlopen as _urlopen) # noqa
98class URLErrorSSL(URLError):
100 def __init__(self, urlerror):
101 self.__dict__ = urlerror.__dict__.copy()
103 def __str__(self):
104 return (
105 'Requesting web resource failed and the problem could be '
106 'related to SSL. Python standard libraries on some older '
107 'systems (like Ubuntu 14.04) are known to have trouble '
108 'with some SSL setups of today\'s servers: %s'
109 % URLError.__str__(self))
112def urlopen_ssl_check(*args, **kwargs):
113 try:
114 return urlopen(*args, **kwargs)
115 except URLError as e:
116 if str(e).find('SSL') != -1:
117 raise URLErrorSSL(e)
118 else:
119 raise
122def urlopen(*args, **kwargs):
123 if 'cafile' not in kwargs:
124 try:
125 import certifi
126 kwargs['cafile'] = certifi.where()
127 except ImportError:
128 pass
130 return _urlopen(*args, **kwargs)
133try:
134 long
135except NameError:
136 long = int
139force_dummy_progressbar = False
142try:
143 from pyrocko import util_ext
144except ImportError:
145 util_ext = None
148logger = logging.getLogger('pyrocko.util')
150try:
151 import progressbar as progressbar_mod
152except ImportError:
153 from pyrocko import dummy_progressbar as progressbar_mod
156try:
157 num_full = num.full
158except AttributeError:
159 def num_full(shape, fill_value, dtype=None, order='C'):
160 a = num.empty(shape, dtype=dtype, order=order)
161 a.fill(fill_value)
162 return a
164try:
165 num_full_like = num.full_like
166except AttributeError:
167 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
168 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
169 a.fill(fill_value)
170 return a
173def progressbar_module():
174 return progressbar_mod
177g_setup_logging_args = 'pyrocko', 'warning'
180def setup_logging(programname='pyrocko', levelname='warning'):
181 '''
182 Initialize logging.
184 :param programname: program name to be written in log
185 :param levelname: string indicating the logging level ('debug', 'info',
186 'warning', 'error', 'critical')
188 This function is called at startup by most pyrocko programs to set up a
189 consistent logging format. This is simply a shortcut to a call to
190 :py:func:`logging.basicConfig()`.
191 '''
193 global g_setup_logging_args
194 g_setup_logging_args = (programname, levelname)
196 levels = {'debug': logging.DEBUG,
197 'info': logging.INFO,
198 'warning': logging.WARNING,
199 'error': logging.ERROR,
200 'critical': logging.CRITICAL}
202 logging.basicConfig(
203 level=levels[levelname],
204 format=programname+':%(name)-20s - %(levelname)-8s - %(message)s')
207def subprocess_setup_logging_args():
208 '''
209 Get arguments from previous call to setup_logging.
211 These can be sent down to a worker process so it can setup its logging
212 in the same way as the main process.
213 '''
214 return g_setup_logging_args
217def data_file(fn):
218 return os.path.join(os.path.split(__file__)[0], 'data', fn)
221class DownloadError(Exception):
222 pass
225class PathExists(DownloadError):
226 pass
229def _download(url, fpath, username=None, password=None,
230 force=False, method='download', stats=None,
231 status_callback=None, entries_wanted=None,
232 recursive=False, header=None):
234 import requests
235 from requests.auth import HTTPBasicAuth
236 from requests.exceptions import HTTPError as req_HTTPError
238 requests.adapters.DEFAULT_RETRIES = 5
239 urljoin = requests.compat.urljoin
241 session = requests.Session()
242 session.header = header
243 session.auth = None if username is None\
244 else HTTPBasicAuth(username, password)
246 status = {
247 'ntotal_files': 0,
248 'nread_files': 0,
249 'ntotal_bytes_all_files': 0,
250 'nread_bytes_all_files': 0,
251 'ntotal_bytes_current_file': 0,
252 'nread_bytes_current_file': 0,
253 'finished': False
254 }
256 try:
257 url_to_size = {}
259 if callable(status_callback):
260 status_callback(status)
262 if not recursive and url.endswith('/'):
263 raise DownloadError(
264 'URL: %s appears to be a directory'
265 ' but recurvise download is False' % url)
267 if recursive and not url.endswith('/'):
268 url += '/'
270 def parse_directory_tree(url, subdir=''):
271 r = session.get(urljoin(url, subdir))
272 r.raise_for_status()
274 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
276 files = sorted(set(subdir + fn for fn in entries
277 if not fn.endswith('/')))
279 if entries_wanted is not None:
280 files = [fn for fn in files
281 if (fn in wanted for wanted in entries_wanted)]
283 status['ntotal_files'] += len(files)
285 dirs = sorted(set(subdir + dn for dn in entries
286 if dn.endswith('/')
287 and dn not in ('./', '../')))
289 for dn in dirs:
290 files.extend(parse_directory_tree(
291 url, subdir=dn))
293 return files
295 def get_content_length(url):
296 if url not in url_to_size:
297 r = session.head(url, headers={'Accept-Encoding': ''})
299 content_length = r.headers.get('content-length', None)
300 if content_length is None:
301 logger.debug('Could not get HTTP header '
302 'Content-Length for %s' % url)
304 content_length = None
306 else:
307 content_length = int(content_length)
308 status['ntotal_bytes_all_files'] += content_length
309 if callable(status_callback):
310 status_callback(status)
312 url_to_size[url] = content_length
314 return url_to_size[url]
316 def download_file(url, fn):
317 logger.info('starting download of %s...' % url)
318 ensuredirs(fn)
320 fsize = get_content_length(url)
321 status['ntotal_bytes_current_file'] = fsize
322 status['nread_bytes_current_file'] = 0
323 if callable(status_callback):
324 status_callback(status)
326 r = session.get(url, stream=True, timeout=5)
327 r.raise_for_status()
329 frx = 0
330 fn_tmp = fn + '.%i.temp' % os.getpid()
331 with open(fn_tmp, 'wb') as f:
332 for d in r.iter_content(chunk_size=1024):
333 f.write(d)
334 frx += len(d)
336 status['nread_bytes_all_files'] += len(d)
337 status['nread_bytes_current_file'] += len(d)
338 if callable(status_callback):
339 status_callback(status)
341 os.rename(fn_tmp, fn)
343 if fsize is not None and frx != fsize:
344 logger.warning(
345 'HTTP header Content-Length: %i bytes does not match '
346 'download size %i bytes' % (fsize, frx))
348 logger.info('finished download of %s' % url)
350 status['nread_files'] += 1
351 if callable(status_callback):
352 status_callback(status)
354 if recursive:
355 if op.exists(fpath) and not force:
356 raise PathExists('path %s already exists' % fpath)
358 files = parse_directory_tree(url)
360 dsize = 0
361 for fn in files:
362 file_url = urljoin(url, fn)
363 dsize += get_content_length(file_url) or 0
365 if method == 'calcsize':
366 return dsize
368 else:
369 for fn in files:
370 file_url = urljoin(url, fn)
371 download_file(file_url, op.join(fpath, fn))
373 else:
374 status['ntotal_files'] += 1
375 if callable(status_callback):
376 status_callback(status)
378 fsize = get_content_length(url)
379 if method == 'calcsize':
380 return fsize
381 else:
382 download_file(url, fpath)
384 except req_HTTPError as e:
385 logging.warn("http error: %s" % e)
386 raise DownloadError('could not download file(s) from: %s' % url)
388 finally:
389 status['finished'] = True
390 if callable(status_callback):
391 status_callback(status)
392 session.close()
395def download_file(
396 url, fpath, username=None, password=None, status_callback=None,
397 **kwargs):
398 return _download(
399 url, fpath, username, password,
400 recursive=False,
401 status_callback=status_callback,
402 **kwargs)
405def download_dir(
406 url, fpath, username=None, password=None, status_callback=None,
407 **kwargs):
409 return _download(
410 url, fpath, username, password,
411 recursive=True,
412 status_callback=status_callback,
413 **kwargs)
416class HPFloatUnavailable(Exception):
417 pass
420class dummy_hpfloat(object):
421 def __init__(self, *args, **kwargs):
422 raise HPFloatUnavailable(
423 'NumPy lacks support for float128 or float96 data type on this '
424 'platform.')
427if hasattr(num, 'float128'):
428 hpfloat = num.float128
429elif hasattr(num, 'float96'):
430 hpfloat = num.float96
431else:
432 hpfloat = dummy_hpfloat
435g_time_float = None
436g_time_dtype = None
439class TimeFloatSettingError(Exception):
440 pass
443def use_high_precision_time(enabled):
444 '''
445 Globally force a specific time handling mode.
447 See :ref:`High precision time handling mode <time-handling-mode>`.
449 :param enabled: enable/disable use of high precision time type
450 :type enabled: bool
452 This function should be called before handling/reading any time data.
453 It can only be called once.
455 Special attention is required when using multiprocessing on a platform
456 which does not use fork under the hood. In such cases, the desired setting
457 must be set also in the subprocess.
458 '''
459 _setup_high_precision_time_mode(enabled_app=enabled)
462def _setup_high_precision_time_mode(enabled_app=False):
463 global g_time_float
464 global g_time_dtype
466 if not (g_time_float is None and g_time_dtype is None):
467 raise TimeFloatSettingError(
468 'Cannot set time handling mode: too late, it has already been '
469 'fixed by an earlier call.')
471 from pyrocko import config
473 conf = config.config()
474 enabled_config = conf.use_high_precision_time
476 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
477 if enabled_env is not None:
478 try:
479 enabled_env = int(enabled_env) == 1
480 except ValueError:
481 raise TimeFloatSettingError(
482 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
483 'should be set to 0 or 1.')
485 enabled = enabled_config
486 mode_from = 'config variable `use_high_precision_time`'
487 notify = enabled
489 if enabled_env is not None:
490 if enabled_env != enabled:
491 notify = True
492 enabled = enabled_env
493 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
495 if enabled_app is not None:
496 if enabled_app != enabled:
497 notify = True
498 enabled = enabled_app
499 mode_from = 'application override'
501 logger.debug('''
502Pyrocko high precision time mode selection (latter override earlier):
503 config: %s
504 env: %s
505 app: %s
506 -> enabled: %s'''.lstrip() % (
507 enabled_config, enabled_env, enabled_app, enabled))
509 if notify:
510 logger.info('Pyrocko high precision time mode %s by %s.' % (
511 'activated' if enabled else 'deactivated',
512 mode_from))
514 if enabled:
515 g_time_float = hpfloat
516 g_time_dtype = hpfloat
517 else:
518 g_time_float = float
519 g_time_dtype = num.float64
522def get_time_float():
523 '''
524 Get the effective float class for timestamps.
526 See :ref:`High precision time handling mode <time-handling-mode>`.
528 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
529 current time handling mode
530 '''
531 global g_time_float
533 if g_time_float is None:
534 _setup_high_precision_time_mode()
536 return g_time_float
539def get_time_dtype():
540 '''
541 Get effective NumPy float class to handle timestamps.
543 See :ref:`High precision time handling mode <time-handling-mode>`.
544 '''
546 global g_time_dtype
548 if g_time_dtype is None:
549 _setup_high_precision_time_mode()
551 return g_time_dtype
554def to_time_float(t):
555 '''
556 Convert float to valid time stamp in the current time handling mode.
558 See :ref:`High precision time handling mode <time-handling-mode>`.
559 '''
560 return get_time_float()(t)
563class TimestampTypeError(ValueError):
564 pass
567def check_time_class(t, error='raise'):
568 '''
569 Type-check variable against current time handling mode.
571 See :ref:`High precision time handling mode <time-handling-mode>`.
572 '''
574 if t == 0.0:
575 return
577 if not isinstance(t, get_time_float()):
578 message = (
579 'Timestamp %g is of type %s but should be of type %s with '
580 'Pyrocko\'s currently selected time handling mode.\n\n'
581 'See https://pyrocko.org/docs/current/library/reference/util.html'
582 '#high-precision-time-handling-mode' % (
583 t, type(t), get_time_float()))
585 if error == 'raise':
586 raise TimestampTypeError(message)
587 elif error == 'warn':
588 logger.warn(message)
589 else:
590 assert False
593class Stopwatch(object):
594 '''
595 Simple stopwatch to measure elapsed wall clock time.
597 Usage::
599 s = Stopwatch()
600 time.sleep(1)
601 print s()
602 time.sleep(1)
603 print s()
604 '''
606 def __init__(self):
607 self.start = time.time()
609 def __call__(self):
610 return time.time() - self.start
613def wrap(text, line_length=80):
614 '''
615 Paragraph and list-aware wrapping of text.
616 '''
618 text = text.strip('\n')
619 at_lineend = re.compile(r' *\n')
620 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
622 paragraphs = at_para.split(text)[::5]
623 listindents = at_para.split(text)[4::5]
624 newlist = at_para.split(text)[3::5]
626 listindents[0:0] = [None]
627 listindents.append(True)
628 newlist.append(None)
630 det_indent = re.compile(r'^ *')
632 outlines = []
633 for ip, p in enumerate(paragraphs):
634 if not p:
635 continue
637 if listindents[ip] is None:
638 _indent = det_indent.findall(p)[0]
639 findent = _indent
640 else:
641 findent = listindents[ip]
642 _indent = ' ' * len(findent)
644 ll = line_length - len(_indent)
645 llf = ll
647 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
648 p1 = ' '.join(oldlines)
649 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
650 for imatch, match in enumerate(possible.finditer(p1)):
651 parout = match.group(1)
652 if imatch == 0:
653 outlines.append(findent + parout)
654 else:
655 outlines.append(_indent + parout)
657 if ip != len(paragraphs)-1 and (
658 listindents[ip] is None
659 or newlist[ip] is not None
660 or listindents[ip+1] is None):
662 outlines.append('')
664 return outlines
667class BetterHelpFormatter(optparse.IndentedHelpFormatter):
669 def __init__(self, *args, **kwargs):
670 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
672 def format_option(self, option):
673 '''
674 From IndentedHelpFormatter but using a different wrap method.
675 '''
677 help_text_position = 4 + self.current_indent
678 help_text_width = self.width - help_text_position
680 opts = self.option_strings[option]
681 opts = "%*s%s" % (self.current_indent, "", opts)
682 if option.help:
683 help_text = self.expand_default(option)
685 if self.help_position + len(help_text) + 1 <= self.width:
686 lines = [
687 '',
688 '%-*s %s' % (self.help_position, opts, help_text),
689 '']
690 else:
691 lines = ['']
692 lines.append(opts)
693 lines.append('')
694 if option.help:
695 help_lines = wrap(help_text, help_text_width)
696 lines.extend(["%*s%s" % (help_text_position, "", line)
697 for line in help_lines])
698 lines.append('')
700 return "\n".join(lines)
702 def format_description(self, description):
703 if not description:
704 return ''
706 if self.current_indent == 0:
707 lines = []
708 else:
709 lines = ['']
711 lines.extend(wrap(description, self.width - self.current_indent))
712 if self.current_indent == 0:
713 lines.append('\n')
715 return '\n'.join(
716 ['%*s%s' % (self.current_indent, '', line) for line in lines])
719def progressbar(label, maxval):
720 progressbar_mod = progressbar_module()
721 if force_dummy_progressbar:
722 progressbar_mod = dummy_progressbar
724 widgets = [
725 label, ' ',
726 progressbar_mod.Bar(marker='-', left='[', right=']'), ' ',
727 progressbar_mod.Percentage(), ' ',
728 progressbar_mod.ETA()]
730 pbar = progressbar_mod.ProgressBar(widgets=widgets, maxval=maxval).start()
731 return pbar
734def progress_beg(label):
735 '''
736 Notify user that an operation has started.
738 :param label: name of the operation
740 To be used in conjuction with :py:func:`progress_end`.
741 '''
743 sys.stderr.write(label)
744 sys.stderr.flush()
747def progress_end(label=''):
748 '''
749 Notify user that an operation has ended.
751 :param label: name of the operation
753 To be used in conjuction with :py:func:`progress_beg`.
754 '''
756 sys.stderr.write(' done. %s\n' % label)
757 sys.stderr.flush()
760class ArangeError(Exception):
761 pass
764def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
765 '''
766 Return evenly spaced numbers over a specified interval.
768 Like :py:func:`numpy.arange` but returning floating point numbers by
769 default and with defined behaviour when stepsize is inconsistent with
770 interval bounds. It is considered inconsistent if the difference between
771 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
772 step``. Inconsistencies are handled according to the ``error`` parameter.
773 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
774 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
775 is silently changed to the closest, the next smaller, or next larger
776 multiple of ``step``, respectively.
777 '''
779 assert error in ('raise', 'round', 'floor', 'ceil')
781 start = dtype(start)
782 stop = dtype(stop)
783 step = dtype(step)
785 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
787 n = int(rnd((stop - start) / step)) + 1
788 stop_check = start + (n-1) * step
790 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
791 raise ArangeError(
792 'inconsistent range specification: start=%g, stop=%g, step=%g'
793 % (start, stop, step))
795 x = num.arange(n, dtype=dtype)
796 x *= step
797 x += start
798 return x
801def polylinefit(x, y, n_or_xnodes):
802 '''
803 Fit piece-wise linear function to data.
805 :param x,y: arrays with coordinates of data
806 :param n_or_xnodes: int, number of segments or x coordinates of polyline
808 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
809 polyline, root-mean-square error
810 '''
812 x = num.asarray(x)
813 y = num.asarray(y)
815 if isinstance(n_or_xnodes, int):
816 n = n_or_xnodes
817 xmin = x.min()
818 xmax = x.max()
819 xnodes = num.linspace(xmin, xmax, n+1)
820 else:
821 xnodes = num.asarray(n_or_xnodes)
822 n = xnodes.size - 1
824 assert len(x) == len(y)
825 assert n > 0
827 ndata = len(x)
828 a = num.zeros((ndata+(n-1), n*2))
829 for i in range(n):
830 xmin_block = xnodes[i]
831 xmax_block = xnodes[i+1]
832 if i == n-1: # don't loose last point
833 indices = num.where(
834 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
835 else:
836 indices = num.where(
837 num.logical_and(xmin_block <= x, x < xmax_block))[0]
839 a[indices, i*2] = x[indices]
840 a[indices, i*2+1] = 1.0
842 w = float(ndata)*100.
843 if i < n-1:
844 a[ndata+i, i*2] = xmax_block*w
845 a[ndata+i, i*2+1] = 1.0*w
846 a[ndata+i, i*2+2] = -xmax_block*w
847 a[ndata+i, i*2+3] = -1.0*w
849 d = num.concatenate((y, num.zeros(n-1)))
850 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
852 ynodes = num.zeros(n+1)
853 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
854 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
855 ynodes[1:n] *= 0.5
857 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
859 return xnodes, ynodes, rms_error
862def plf_integrate_piecewise(x_edges, x, y):
863 '''
864 Calculate definite integral of piece-wise linear function on intervals.
866 Use trapezoidal rule to calculate definite integral of a piece-wise linear
867 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
868 be sorted.
870 :param x_edges: array with edges of the intervals
871 :param x,y: arrays with coordinates of piece-wise linear function's
872 control points
873 '''
875 x_all = num.concatenate((x, x_edges))
876 ii = num.argsort(x_all)
877 y_edges = num.interp(x_edges, x, y)
878 y_all = num.concatenate((y, y_edges))
879 xs = x_all[ii]
880 ys = y_all[ii]
881 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
882 return num.diff(y_all[-len(y_edges):])
885def diff_fd_1d_4o(dt, data):
886 '''
887 Approximate first derivative of an array (forth order, central FD).
889 :param dt: sampling interval
890 :param data: NumPy array with data samples
892 :returns: NumPy array with same shape as input
894 Interior points are approximated to fourth order, edge points to first
895 order right- or left-sided respectively, points next to edge to second
896 order central.
897 '''
898 import scipy.signal
900 ddata = num.empty_like(data, dtype=float)
902 if data.size >= 5:
903 ddata[2:-2] = scipy.signal.lfilter(
904 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
906 if data.size >= 3:
907 ddata[1] = (data[2] - data[0]) / (2. * dt)
908 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
910 if data.size >= 2:
911 ddata[0] = (data[1] - data[0]) / dt
912 ddata[-1] = (data[-1] - data[-2]) / dt
914 if data.size == 1:
915 ddata[0] = 0.0
917 return ddata
920def diff_fd_1d_2o(dt, data):
921 '''
922 Approximate first derivative of an array (second order, central FD).
924 :param dt: sampling interval
925 :param data: NumPy array with data samples
927 :returns: NumPy array with same shape as input
929 Interior points are approximated to second order, edge points to first
930 order right- or left-sided respectively.
932 Uses :py:func:`numpy.gradient`.
933 '''
935 return num.gradient(data, dt)
938def diff_fd_2d_4o(dt, data):
939 '''
940 Approximate second derivative of an array (forth order, central FD).
942 :param dt: sampling interval
943 :param data: NumPy array with data samples
945 :returns: NumPy array with same shape as input
947 Interior points are approximated to fourth order, next-to-edge points to
948 second order, edge points repeated.
949 '''
950 import scipy.signal
952 ddata = num.empty_like(data, dtype=float)
954 if data.size >= 5:
955 ddata[2:-2] = scipy.signal.lfilter(
956 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
958 if data.size >= 3:
959 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
960 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
962 if data.size < 3:
963 ddata[:] = 0.0
965 return ddata
968def diff_fd_2d_2o(dt, data):
969 '''
970 Approximate second derivative of an array (second order, central FD).
972 :param dt: sampling interval
973 :param data: NumPy array with data samples
975 :returns: NumPy array with same shape as input
977 Interior points are approximated to second order, edge points repeated.
978 '''
979 import scipy.signal
981 ddata = num.empty_like(data, dtype=float)
983 if data.size >= 3:
984 ddata[1:-1] = scipy.signal.lfilter(
985 [1., -2., 1.], [1.], data)[2:] / (dt**2)
987 ddata[0] = ddata[1]
988 ddata[-1] = ddata[-2]
990 if data.size < 3:
991 ddata[:] = 0.0
993 return ddata
996def diff_fd(n, order, dt, data):
997 '''
998 Approximate 1st or 2nd derivative of an array.
1000 :param n: 1 for first derivative, 2 for second
1001 :param order: order of the approximation 2 and 4 are supported
1002 :param dt: sampling interval
1003 :param data: NumPy array with data samples
1005 :returns: NumPy array with same shape as input
1007 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1008 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1009 :py:func:`diff_fd_2d_4o`.
1011 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1012 '''
1014 funcs = {
1015 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1016 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1018 try:
1019 funcs_n = funcs[n]
1020 except KeyError:
1021 raise ValueError(
1022 'pyrocko.util.diff_fd: '
1023 'Only 1st and 2sd derivatives are supported.')
1025 try:
1026 func = funcs_n[order]
1027 except KeyError:
1028 raise ValueError(
1029 'pyrocko.util.diff_fd: '
1030 'Order %i is not supported for %s derivative. Supported: %s' % (
1031 order, ['', '1st', '2nd'][n],
1032 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1034 return func(dt, data)
1037class GlobalVars(object):
1038 reuse_store = dict()
1039 decitab_nmax = 0
1040 decitab = {}
1041 decimate_fir_coeffs = {}
1042 decimate_iir_coeffs = {}
1043 re_frac = None
1046def decimate_coeffs(q, n=None, ftype='iir'):
1048 q = int(q)
1050 if n is None:
1051 if ftype == 'fir':
1052 n = 30
1053 else:
1054 n = 8
1056 if ftype == 'fir':
1057 coeffs = GlobalVars.decimate_fir_coeffs
1058 if (n, 1./q) not in coeffs:
1059 coeffs[n, 1./q] = signal.firwin(n+1, 1./q, window='hamming')
1061 b = coeffs[n, 1./q]
1062 return b, [1.], n
1064 else:
1065 coeffs = GlobalVars.decimate_iir_coeffs
1066 if (n, 0.05, 0.8/q) not in coeffs:
1067 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1069 b, a = coeffs[n, 0.05, 0.8/q]
1070 return b, a, n
1073def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1074 '''
1075 Downsample the signal x by an integer factor q, using an order n filter
1077 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1078 filter with hamming window if ftype is 'fir'.
1080 :param x: the signal to be downsampled (1D NumPy array)
1081 :param q: the downsampling factor
1082 :param n: order of the filter (1 less than the length of the filter for a
1083 'fir' filter)
1084 :param ftype: type of the filter; can be 'iir' or 'fir'
1086 :returns: the downsampled signal (1D NumPy array)
1088 '''
1090 b, a, n = decimate_coeffs(q, n, ftype)
1092 if zi is None or zi is True:
1093 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1094 else:
1095 zi_ = zi
1097 y, zf = signal.lfilter(b, a, x, zi=zi_)
1099 if zi is not None:
1100 return y[n//2+ioff::q].copy(), zf
1101 else:
1102 return y[n//2+ioff::q].copy()
1105class UnavailableDecimation(Exception):
1106 '''
1107 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1108 '''
1110 pass
1113def gcd(a, b, epsilon=1e-7):
1114 '''
1115 Greatest common divisor.
1116 '''
1118 while b > epsilon*a:
1119 a, b = b, a % b
1121 return a
1124def lcm(a, b):
1125 '''
1126 Least common multiple.
1127 '''
1129 return a*b // gcd(a, b)
1132def mk_decitab(nmax=100):
1133 '''
1134 Make table with decimation sequences.
1136 Decimation from one sampling rate to a lower one is achieved by a
1137 successive application of :py:func:`decimation` with small integer
1138 downsampling factors (because using large downampling factors can make the
1139 decimation unstable or slow). This function sets up a table with downsample
1140 sequences for factors up to ``nmax``.
1141 '''
1143 tab = GlobalVars.decitab
1144 for i in range(1, 10):
1145 for j in range(1, i+1):
1146 for k in range(1, j+1):
1147 for l_ in range(1, k+1):
1148 for m in range(1, l_+1):
1149 p = i*j*k*l_*m
1150 if p > nmax:
1151 break
1152 if p not in tab:
1153 tab[p] = (i, j, k, l_, m)
1154 if i*j*k*l_ > nmax:
1155 break
1156 if i*j*k > nmax:
1157 break
1158 if i*j > nmax:
1159 break
1160 if i > nmax:
1161 break
1163 GlobalVars.decitab_nmax = nmax
1166def zfmt(n):
1167 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1170def _year_to_time(year):
1171 tt = (year, 1, 1, 0, 0, 0)
1172 return to_time_float(calendar.timegm(tt))
1175def _working_year(year):
1176 try:
1177 tt = (year, 1, 1, 0, 0, 0)
1178 t = calendar.timegm(tt)
1179 tt2_ = time.gmtime(t)
1180 tt2 = tuple(tt2_)[:6]
1181 if tt != tt2:
1182 return False
1184 s = '%i-01-01 00:00:00' % year
1185 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1186 if s != s2:
1187 return False
1189 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1190 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1191 if s3 != s2:
1192 return False
1194 except Exception:
1195 return False
1197 return True
1200def working_system_time_range(year_min_lim=None, year_max_lim=None):
1201 '''
1202 Check time range supported by the systems's time conversion functions.
1204 Returns system time stamps of start of year of first/last fully supported
1205 year span. If this is before 1900 or after 2100, return first/last century
1206 which is fully supported.
1208 :returns: ``(tmin, tmax, year_min, year_max)``
1209 '''
1211 year0 = 2000
1212 year_min = year0
1213 year_max = year0
1215 itests = list(range(101))
1216 for i in range(19):
1217 itests.append(200 + i*100)
1219 for i in itests:
1220 year = year0 - i
1221 if year_min_lim is not None and year < year_min_lim:
1222 year_min = year_min_lim
1223 break
1224 elif not _working_year(year):
1225 break
1226 else:
1227 year_min = year
1229 for i in itests:
1230 year = year0 + i
1231 if year_max_lim is not None and year > year_max_lim:
1232 year_max = year_max_lim
1233 break
1234 elif not _working_year(year + 1):
1235 break
1236 else:
1237 year_max = year
1239 return (
1240 _year_to_time(year_min),
1241 _year_to_time(year_max),
1242 year_min, year_max)
1245g_working_system_time_range = None
1248def get_working_system_time_range():
1249 '''
1250 Caching variant of :py:func:`working_system_time_range`.
1251 '''
1253 global g_working_system_time_range
1254 if g_working_system_time_range is None:
1255 g_working_system_time_range = working_system_time_range()
1257 return g_working_system_time_range
1260def is_working_time(t):
1261 tmin, tmax, _, _ = get_working_system_time_range()
1262 return tmin <= t <= tmax
1265def julian_day_of_year(timestamp):
1266 '''
1267 Get the day number after the 1st of January of year in ``timestamp``.
1269 :returns: day number as int
1270 '''
1272 return time.gmtime(int(timestamp)).tm_yday
1275def hour_start(timestamp):
1276 '''
1277 Get beginning of hour for any point in time.
1279 :param timestamp: time instant as system timestamp (in seconds)
1281 :returns: instant of hour start as system timestamp
1282 '''
1284 tt = time.gmtime(int(timestamp))
1285 tts = tt[0:4] + (0, 0)
1286 return to_time_float(calendar.timegm(tts))
1289def day_start(timestamp):
1290 '''
1291 Get beginning of day for any point in time.
1293 :param timestamp: time instant as system timestamp (in seconds)
1295 :returns: instant of day start as system timestamp
1296 '''
1298 tt = time.gmtime(int(timestamp))
1299 tts = tt[0:3] + (0, 0, 0)
1300 return to_time_float(calendar.timegm(tts))
1303def month_start(timestamp):
1304 '''
1305 Get beginning of month for any point in time.
1307 :param timestamp: time instant as system timestamp (in seconds)
1309 :returns: instant of month start as system timestamp
1310 '''
1312 tt = time.gmtime(int(timestamp))
1313 tts = tt[0:2] + (1, 0, 0, 0)
1314 return to_time_float(calendar.timegm(tts))
1317def year_start(timestamp):
1318 '''
1319 Get beginning of year for any point in time.
1321 :param timestamp: time instant as system timestamp (in seconds)
1323 :returns: instant of year start as system timestamp
1324 '''
1326 tt = time.gmtime(int(timestamp))
1327 tts = tt[0:1] + (1, 1, 0, 0, 0)
1328 return to_time_float(calendar.timegm(tts))
1331def iter_days(tmin, tmax):
1332 '''
1333 Yields begin and end of days until given time span is covered.
1335 :param tmin,tmax: input time span
1337 :yields: tuples with (begin, end) of days as system timestamps
1338 '''
1340 t = day_start(tmin)
1341 while t < tmax:
1342 tend = day_start(t + 26*60*60)
1343 yield t, tend
1344 t = tend
1347def iter_months(tmin, tmax):
1348 '''
1349 Yields begin and end of months until given time span is covered.
1351 :param tmin,tmax: input time span
1353 :yields: tuples with (begin, end) of months as system timestamps
1354 '''
1356 t = month_start(tmin)
1357 while t < tmax:
1358 tend = month_start(t + 24*60*60*33)
1359 yield t, tend
1360 t = tend
1363def iter_years(tmin, tmax):
1364 '''
1365 Yields begin and end of years until given time span is covered.
1367 :param tmin,tmax: input time span
1369 :yields: tuples with (begin, end) of years as system timestamps
1370 '''
1372 t = year_start(tmin)
1373 while t < tmax:
1374 tend = year_start(t + 24*60*60*369)
1375 yield t, tend
1376 t = tend
1379def decitab(n):
1380 '''
1381 Get integer decimation sequence for given downampling factor.
1383 :param n: target decimation factor
1385 :returns: tuple with downsampling sequence
1386 '''
1388 if n > GlobalVars.decitab_nmax:
1389 mk_decitab(n*2)
1390 if n not in GlobalVars.decitab:
1391 raise UnavailableDecimation('ratio = %g' % n)
1392 return GlobalVars.decitab[n]
1395def ctimegm(s, format="%Y-%m-%d %H:%M:%S"):
1396 '''
1397 Convert string representing UTC time to system time.
1399 :param s: string to be interpreted
1400 :param format: format string passed to :py:func:`strptime`
1402 :returns: system time stamp
1404 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1406 .. note::
1407 This function is to be replaced by :py:func:`str_to_time`.
1408 '''
1410 return calendar.timegm(time.strptime(s, format))
1413def gmctime(t, format="%Y-%m-%d %H:%M:%S"):
1414 '''
1415 Get string representation from system time, UTC.
1417 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1419 .. note::
1420 This function is to be repaced by :py:func:`time_to_str`.
1421 '''
1423 return time.strftime(format, time.gmtime(t))
1426def gmctime_v(t, format="%a, %d %b %Y %H:%M:%S"):
1427 '''
1428 Get string representation from system time, UTC. Same as
1429 :py:func:`gmctime` but with a more verbose default format.
1431 .. note::
1432 This function is to be replaced by :py:func:`time_to_str`.
1433 '''
1435 return time.strftime(format, time.gmtime(t))
1438def gmctime_fn(t, format="%Y-%m-%d_%H-%M-%S"):
1439 '''
1440 Get string representation from system time, UTC. Same as
1441 :py:func:`gmctime` but with a default usable in filenames.
1443 .. note::
1444 This function is to be replaced by :py:func:`time_to_str`.
1445 '''
1447 return time.strftime(format, time.gmtime(t))
1450class TimeStrError(Exception):
1451 pass
1454class FractionalSecondsMissing(TimeStrError):
1455 '''
1456 Exception raised by :py:func:`str_to_time` when the given string lacks
1457 fractional seconds.
1458 '''
1460 pass
1463class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1464 '''
1465 Exception raised by :py:func:`str_to_time` when the given string has an
1466 incorrect number of digits in the fractional seconds part.
1467 '''
1469 pass
1472def _endswith_n(s, endings):
1473 for ix, x in enumerate(endings):
1474 if s.endswith(x):
1475 return ix
1476 return -1
1479def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1480 '''
1481 Convert string representing UTC time to floating point system time.
1483 :param s: string representing UTC time
1484 :param format: time string format
1485 :returns: system time stamp as floating point value
1487 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1488 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1489 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1490 the fractional part, including the dot is made optional. The latter has the
1491 consequence, that the time strings and the format may not contain any other
1492 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1493 ensured, that exactly that number of digits are present in the fractional
1494 seconds.
1495 '''
1497 time_float = get_time_float()
1499 if util_ext is not None:
1500 try:
1501 t, tfrac = util_ext.stt(s, format)
1502 except util_ext.UtilExtError as e:
1503 raise TimeStrError(
1504 '%s, string=%s, format=%s' % (str(e), s, format))
1506 return time_float(t)+tfrac
1508 fracsec = 0.
1509 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1511 iend = _endswith_n(format, fixed_endings)
1512 if iend != -1:
1513 dotpos = s.rfind('.')
1514 if dotpos == -1:
1515 raise FractionalSecondsMissing(
1516 'string=%s, format=%s' % (s, format))
1518 if iend > 0 and iend != (len(s)-dotpos-1):
1519 raise FractionalSecondsWrongNumberOfDigits(
1520 'string=%s, format=%s' % (s, format))
1522 format = format[:-len(fixed_endings[iend])]
1523 fracsec = float(s[dotpos:])
1524 s = s[:dotpos]
1526 elif format.endswith('.OPTFRAC'):
1527 dotpos = s.rfind('.')
1528 format = format[:-8]
1529 if dotpos != -1 and len(s[dotpos:]) > 1:
1530 fracsec = float(s[dotpos:])
1532 if dotpos != -1:
1533 s = s[:dotpos]
1535 try:
1536 return time_float(calendar.timegm(time.strptime(s, format))) \
1537 + fracsec
1538 except ValueError as e:
1539 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1542stt = str_to_time
1545def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1546 '''
1547 Get string representation for floating point system time.
1549 :param t: floating point system time
1550 :param format: time string format
1551 :returns: string representing UTC time
1553 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1554 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1555 digit between 1 and 9, this is replaced with the fractional part of ``t``
1556 with ``x`` digits precision.
1557 '''
1559 if pyrocko.grumpy > 0:
1560 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1562 if isinstance(format, int):
1563 if format > 0:
1564 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1565 else:
1566 format = '%Y-%m-%d %H:%M:%S'
1568 if util_ext is not None:
1569 t0 = num.floor(t)
1570 try:
1571 return util_ext.tts(int(t0), float(t - t0), format)
1572 except util_ext.UtilExtError as e:
1573 raise TimeStrError(
1574 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1576 if not GlobalVars.re_frac:
1577 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1578 GlobalVars.frac_formats = dict(
1579 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1581 ts = float(num.floor(t))
1582 tfrac = t-ts
1584 m = GlobalVars.re_frac.search(format)
1585 if m:
1586 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1587 if sfrac[0] == '1':
1588 ts += 1.
1590 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1592 return time.strftime(format, time.gmtime(ts))
1595tts = time_to_str
1598def mystrftime(fmt=None, tt=None, milliseconds=0):
1600 if fmt is None:
1601 fmt = '%Y-%m-%d %H:%M:%S .%r'
1603 if tt is None:
1604 tt = time.time()
1606 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1607 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1608 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1609 return time.strftime(fmt, tt)
1612def gmtime_x(timestamp):
1613 etimestamp = float(num.floor(timestamp))
1614 tt = time.gmtime(etimestamp)
1615 ms = (timestamp-etimestamp)*1000
1616 return tt, ms
1619def plural_s(n):
1620 if n == 1:
1621 return ''
1622 else:
1623 return 's'
1626def ensuredirs(dst):
1627 '''
1628 Create all intermediate path components for a target path.
1630 :param dst: target path
1632 The leaf part of the target path is not created (use :py:func:`ensuredir`
1633 if a the target path is a directory to be created).
1634 '''
1636 d, x = os.path.split(dst.rstrip(os.sep))
1637 dirs = []
1638 while d and not os.path.exists(d):
1639 dirs.append(d)
1640 d, x = os.path.split(d)
1642 dirs.reverse()
1644 for d in dirs:
1645 try:
1646 os.mkdir(d)
1647 except OSError as e:
1648 if not e.errno == errno.EEXIST:
1649 raise
1652def ensuredir(dst):
1653 '''
1654 Create directory and all intermediate path components to it as needed.
1656 :param dst: directory name
1658 Nothing is done if the given target already exists.
1659 '''
1661 if os.path.exists(dst):
1662 return
1664 dst.rstrip(os.sep)
1666 ensuredirs(dst)
1667 try:
1668 os.mkdir(dst)
1669 except OSError as e:
1670 if not e.errno == errno.EEXIST:
1671 raise
1674def reuse(x):
1675 '''
1676 Get unique instance of an object.
1678 :param x: hashable object
1679 :returns: reference to x or an equivalent object
1681 Cache object ``x`` in a global dict for reuse, or if x already
1682 is in that dict, return a reference to it.
1683 '''
1685 grs = GlobalVars.reuse_store
1686 if x not in grs:
1687 grs[x] = x
1688 return grs[x]
1691def deuse(x):
1692 grs = GlobalVars.reuse_store
1693 if x in grs:
1694 del grs[x]
1697class Anon(object):
1698 '''
1699 Dict-to-object utility.
1701 Any given arguments are stored as attributes.
1703 Example::
1705 a = Anon(x=1, y=2)
1706 print a.x, a.y
1707 '''
1709 def __init__(self, **dict):
1710 for k in dict:
1711 self.__dict__[k] = dict[k]
1714def select_files(paths, selector=None, regex=None, show_progress=True):
1715 '''
1716 Recursively select files.
1718 :param paths: entry path names
1719 :param selector: callback for conditional inclusion
1720 :param regex: pattern for conditional inclusion
1721 :param show_progress: if True, indicate start and stop of processing
1722 :returns: list of path names
1724 Recursively finds all files under given entry points ``paths``. If
1725 parameter ``regex`` is a regular expression, only files with matching path
1726 names are included. If additionally parameter ``selector`` is given a
1727 callback function, only files for which the callback returns ``True`` are
1728 included. The callback should take a single argument. The callback is
1729 called with a single argument, an object, having as attributes, any named
1730 groups given in ``regex``.
1732 Examples
1734 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1736 select_files(paths,
1737 regex=r'\\.(mseed|msd)$')
1739 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1740 the year::
1742 select_files(paths,
1743 regex=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1744 selector=(lambda x: int(x.year) == 2009))
1745 '''
1747 if show_progress:
1748 progress_beg('selecting files...')
1749 if logger.isEnabledFor(logging.DEBUG):
1750 sys.stderr.write('\n')
1752 good = []
1753 if regex:
1754 rselector = re.compile(regex)
1756 def addfile(path):
1757 if regex:
1758 logger.debug("looking at filename: '%s'" % path)
1759 m = rselector.search(path)
1760 if m:
1761 infos = Anon(**m.groupdict())
1762 logger.debug(" regex '%s' matches." % regex)
1763 for k, v in m.groupdict().items():
1764 logger.debug(
1765 " attribute '%s' has value '%s'" % (k, v))
1766 if selector is None or selector(infos):
1767 good.append(os.path.abspath(path))
1769 else:
1770 logger.debug(" regex '%s' does not match." % regex)
1771 else:
1772 good.append(os.path.abspath(path))
1774 if isinstance(paths, str):
1775 paths = [paths]
1777 for path in paths:
1778 if os.path.isdir(path):
1779 for (dirpath, dirnames, filenames) in os.walk(path):
1780 for filename in filenames:
1781 addfile(op.join(dirpath, filename))
1782 else:
1783 addfile(path)
1785 if show_progress:
1786 progress_end('%i file%s selected.' % (len(good), plural_s(len(good))))
1788 return good
1791def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1792 '''
1793 Convert positive integer to a base36 string.
1794 '''
1796 if not isinstance(number, (int, long)):
1797 raise TypeError('number must be an integer')
1798 if number < 0:
1799 raise ValueError('number must be positive')
1801 # Special case for small numbers
1802 if number < 36:
1803 return alphabet[number]
1805 base36 = ''
1806 while number != 0:
1807 number, i = divmod(number, 36)
1808 base36 = alphabet[i] + base36
1810 return base36
1813def base36decode(number):
1814 '''
1815 Decode base36 endcoded positive integer.
1816 '''
1818 return int(number, 36)
1821class UnpackError(Exception):
1822 '''
1823 Exception raised when :py:func:`unpack_fixed` encounters an error.
1824 '''
1826 pass
1829ruler = ''.join(['%-10i' % i for i in range(8)]) \
1830 + '\n' + '0123456789' * 8 + '\n'
1833def unpack_fixed(format, line, *callargs):
1834 '''
1835 Unpack fixed format string, as produced by many fortran codes.
1837 :param format: format specification
1838 :param line: string to be processed
1839 :param callargs: callbacks for callback fields in the format
1841 The format is described by a string of comma-separated fields. Each field
1842 is defined by a character for the field type followed by the field width. A
1843 questionmark may be appended to the field description to allow the argument
1844 to be optional (The data string is then allowed to be filled with blanks
1845 and ``None`` is returned in this case).
1847 The following field types are available:
1849 ==== ================================================================
1850 Type Description
1851 ==== ================================================================
1852 A string (full field width is extracted)
1853 a string (whitespace at the beginning and the end is removed)
1854 i integer value
1855 f floating point value
1856 @ special type, a callback must be given for the conversion
1857 x special field type to skip parts of the string
1858 ==== ================================================================
1859 '''
1861 ipos = 0
1862 values = []
1863 icall = 0
1864 for form in format.split(','):
1865 form = form.strip()
1866 optional = form[-1] == '?'
1867 form = form.rstrip('?')
1868 typ = form[0]
1869 ln = int(form[1:])
1870 s = line[ipos:ipos+ln]
1871 cast = {
1872 'x': None,
1873 'A': str,
1874 'a': lambda x: x.strip(),
1875 'i': int,
1876 'f': float,
1877 '@': 'extra'}[typ]
1879 if cast == 'extra':
1880 cast = callargs[icall]
1881 icall += 1
1883 if cast is not None:
1884 if optional and s.strip() == '':
1885 values.append(None)
1886 else:
1887 try:
1888 values.append(cast(s))
1889 except Exception:
1890 mark = [' '] * 80
1891 mark[ipos:ipos+ln] = ['^'] * ln
1892 mark = ''.join(mark)
1893 raise UnpackError(
1894 'Invalid cast to type "%s" at position [%i:%i] of '
1895 'line: \n%s%s\n%s' % (
1896 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
1898 ipos += ln
1900 return values
1903_pattern_cache = {}
1906def _nslc_pattern(pattern):
1907 if pattern not in _pattern_cache:
1908 rpattern = re.compile(fnmatch.translate(pattern), re.I)
1909 _pattern_cache[pattern] = rpattern
1910 else:
1911 rpattern = _pattern_cache[pattern]
1913 return rpattern
1916def match_nslc(patterns, nslc):
1917 '''
1918 Match network-station-location-channel code against pattern or list of
1919 patterns.
1921 :param patterns: pattern or list of patterns
1922 :param nslc: tuple with (network, station, location, channel) as strings
1924 :returns: ``True`` if the pattern matches or if any of the given patterns
1925 match; or ``False``.
1927 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
1929 Example::
1931 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
1932 '''
1934 if isinstance(patterns, str):
1935 patterns = [patterns]
1937 if not isinstance(nslc, str):
1938 s = '.'.join(nslc)
1939 else:
1940 s = nslc
1942 for pattern in patterns:
1943 if _nslc_pattern(pattern).match(s):
1944 return True
1946 return False
1949def match_nslcs(patterns, nslcs):
1950 '''
1951 Get network-station-location-channel codes that match given pattern or any
1952 of several given patterns.
1954 :param patterns: pattern or list of patterns
1955 :param nslcs: list of (network, station, location, channel) tuples
1957 See also :py:func:`match_nslc`
1958 '''
1960 matching = []
1961 for nslc in nslcs:
1962 if match_nslc(patterns, nslc):
1963 matching.append(nslc)
1965 return matching
1968class Timeout(Exception):
1969 pass
1972def create_lockfile(fn, timeout=None, timewarn=10.):
1973 t0 = time.time()
1975 while True:
1976 try:
1977 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
1978 os.close(f)
1979 return
1981 except OSError as e:
1982 if e.errno == errno.EEXIST:
1983 tnow = time.time()
1985 if timeout is not None and tnow - t0 > timeout:
1986 raise Timeout(
1987 'Timeout (%gs) occured while waiting to get exclusive '
1988 'access to: %s' % (timeout, fn))
1990 if timewarn is not None and tnow - t0 > timewarn:
1991 logger.warn(
1992 'Waiting since %gs to get exclusive access to: %s' % (
1993 timewarn, fn))
1995 timewarn *= 2
1997 time.sleep(0.01)
1998 else:
1999 raise
2002def delete_lockfile(fn):
2003 os.unlink(fn)
2006class Lockfile(Exception):
2008 def __init__(self, path, timeout=5, timewarn=10.):
2009 self._path = path
2010 self._timeout = timeout
2011 self._timewarn = timewarn
2013 def __enter__(self):
2014 create_lockfile(
2015 self._path, timeout=self._timeout, timewarn=self._timewarn)
2016 return None
2018 def __exit__(self, type, value, traceback):
2019 delete_lockfile(self._path)
2022class SoleError(Exception):
2023 '''
2024 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2025 instance is running.
2026 '''
2028 pass
2031class Sole(object):
2032 '''
2033 Use POSIX advisory file locking to ensure that only a single instance of a
2034 program is running.
2036 :param pid_path: path to lockfile to be used
2038 Usage::
2040 from pyrocko.util import Sole, SoleError, setup_logging
2041 import os
2043 setup_logging('my_program')
2045 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2046 try:
2047 sole = Sole(pid_path)
2049 except SoleError, e:
2050 logger.fatal( str(e) )
2051 sys.exit(1)
2052 '''
2054 def __init__(self, pid_path):
2055 self._pid_path = pid_path
2056 self._other_running = False
2057 ensuredirs(self._pid_path)
2058 self._lockfile = None
2059 self._os = os
2061 if not fcntl:
2062 raise SoleError(
2063 'Python standard library module "fcntl" not available on '
2064 'this platform.')
2066 self._fcntl = fcntl
2068 try:
2069 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2070 except Exception:
2071 raise SoleError(
2072 'Cannot open lockfile (path = %s)' % self._pid_path)
2074 try:
2075 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2077 except IOError:
2078 self._other_running = True
2079 try:
2080 f = open(self._pid_path, 'r')
2081 pid = f.read().strip()
2082 f.close()
2083 except Exception:
2084 pid = '?'
2086 raise SoleError('Other instance is running (pid = %s)' % pid)
2088 try:
2089 os.ftruncate(self._lockfile, 0)
2090 os.write(self._lockfile, '%i\n' % os.getpid())
2091 os.fsync(self._lockfile)
2093 except Exception:
2094 # the pid is only stored for user information, so this is allowed
2095 # to fail
2096 pass
2098 def __del__(self):
2099 if not self._other_running:
2100 if self._lockfile is not None:
2101 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2102 self._os.close(self._lockfile)
2103 try:
2104 self._os.unlink(self._pid_path)
2105 except Exception:
2106 pass
2109re_escapequotes = re.compile(r"(['\\])")
2112def escapequotes(s):
2113 return re_escapequotes.sub(r"\\\1", s)
2116class TableWriter(object):
2117 '''
2118 Write table of space separated values to a file.
2120 :param f: file like object
2122 Strings containing spaces are quoted on output.
2123 '''
2125 def __init__(self, f):
2126 self._f = f
2128 def writerow(self, row, minfieldwidths=None):
2130 '''
2131 Write one row of values to underlying file.
2133 :param row: iterable of values
2134 :param minfieldwidths: minimum field widths for the values
2136 Each value in in ``row`` is converted to a string and optionally padded
2137 with blanks. The resulting strings are output separated with blanks. If
2138 any values given are strings and if they contain whitespace, they are
2139 quoted with single quotes, and any internal single quotes are
2140 backslash-escaped.
2141 '''
2143 out = []
2145 for i, x in enumerate(row):
2146 w = 0
2147 if minfieldwidths and i < len(minfieldwidths):
2148 w = minfieldwidths[i]
2150 if isinstance(x, str):
2151 if re.search(r"\s|'", x):
2152 x = "'%s'" % escapequotes(x)
2154 x = x.ljust(w)
2155 else:
2156 x = str(x).rjust(w)
2158 out.append(x)
2160 self._f.write(' '.join(out).rstrip() + '\n')
2163class TableReader(object):
2165 '''
2166 Read table of space separated values from a file.
2168 :param f: file-like object
2170 This uses Pythons shlex module to tokenize lines. Should deal correctly
2171 with quoted strings.
2172 '''
2174 def __init__(self, f):
2175 self._f = f
2176 self.eof = False
2178 def readrow(self):
2179 '''
2180 Read one row from the underlying file, tokenize it with shlex.
2182 :returns: tokenized line as a list of strings.
2183 '''
2185 line = self._f.readline()
2186 if not line:
2187 self.eof = True
2188 return []
2189 line.strip()
2190 if line.startswith('#'):
2191 return []
2193 return qsplit(line)
2196def gform(number, significant_digits=3):
2197 '''
2198 Pretty print floating point numbers.
2200 Align floating point numbers at the decimal dot.
2202 ::
2204 | -d.dde+xxx|
2205 | -d.dde+xx |
2206 |-ddd. |
2207 | -dd.d |
2208 | -d.dd |
2209 | -0.ddd |
2210 | -0.0ddd |
2211 | -0.00ddd |
2212 | -d.dde-xx |
2213 | -d.dde-xxx|
2214 | nan|
2217 The formatted string has length ``significant_digits * 2 + 6``.
2218 '''
2220 no_exp_range = (pow(10., -1),
2221 pow(10., significant_digits))
2222 width = significant_digits+significant_digits-1+1+1+5
2224 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2225 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2226 else:
2227 s = '%.*E' % (significant_digits-1, number)
2228 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2229 if s.strip().lower() == 'nan':
2230 s = 'nan'.rjust(width)
2231 return s
2234def human_bytesize(value):
2236 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2238 if value == 1:
2239 return '1 Byte'
2241 for i, ext in enumerate(exts):
2242 x = float(value) / 1000**i
2243 if round(x) < 10. and not value < 1000:
2244 return '%.1f %s' % (x, ext)
2245 if round(x) < 1000.:
2246 return '%.0f %s' % (x, ext)
2248 return '%i Bytes' % value
2251re_compatibility = re.compile(
2252 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2253 r'(dummy|poel|qseis|qssp))\.'
2254)
2257def pf_is_old(fn):
2258 oldstyle = False
2259 with open(fn, 'r') as f:
2260 for line in f:
2261 if re_compatibility.search(line):
2262 oldstyle = True
2264 return oldstyle
2267def pf_upgrade(fn):
2268 need = pf_is_old(fn)
2269 if need:
2270 fn_temp = fn + '.temp'
2272 with open(fn, 'r') as fin:
2273 with open(fn_temp, 'w') as fout:
2274 for line in fin:
2275 line = re_compatibility.sub('!pf.', line)
2276 fout.write(line)
2278 os.rename(fn_temp, fn)
2280 return need
2283def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2284 '''
2285 Extract leap second information from tzdata.
2287 Based on example at http://stackoverflow.com/questions/19332902/\
2288 extract-historic-leap-seconds-from-tzdata
2290 See also 'man 5 tzfile'.
2291 '''
2293 from struct import unpack, calcsize
2294 out = []
2295 with open(tzfile, 'rb') as f:
2296 # read header
2297 fmt = '>4s c 15x 6l'
2298 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2299 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2300 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2302 # skip over some uninteresting data
2303 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2304 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2305 f.read(calcsize(fmt))
2307 # read leap-seconds
2308 fmt = '>2l'
2309 for i in range(leapcnt):
2310 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2311 out.append((tleap-nleap+1, nleap))
2313 return out
2316class LeapSecondsError(Exception):
2317 pass
2320class LeapSecondsOutdated(LeapSecondsError):
2321 pass
2324class InvalidLeapSecondsFile(LeapSecondsOutdated):
2325 pass
2328def parse_leap_seconds_list(fn):
2329 data = []
2330 texpires = None
2331 try:
2332 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2333 except TimeStrError:
2334 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2336 tnow = int(round(time.time()))
2338 if not op.exists(fn):
2339 raise LeapSecondsOutdated('no leap seconds file found')
2341 try:
2342 with open(fn, 'rb') as f:
2343 for line in f:
2344 if line.strip().startswith(b'<!DOCTYPE'):
2345 raise InvalidLeapSecondsFile('invalid leap seconds file')
2347 if line.startswith(b'#@'):
2348 texpires = int(line.split()[1]) + t0
2349 elif line.startswith(b'#') or len(line) < 5:
2350 pass
2351 else:
2352 toks = line.split()
2353 t = int(toks[0]) + t0
2354 nleap = int(toks[1]) - 10
2355 data.append((t, nleap))
2357 except IOError:
2358 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2360 if texpires is None or tnow > texpires:
2361 raise LeapSecondsOutdated('leap seconds list is outdated')
2363 return data
2366def read_leap_seconds2():
2367 from pyrocko import config
2368 conf = config.config()
2369 fn = conf.leapseconds_path
2370 url = conf.leapseconds_url
2371 # check for outdated default URL
2372 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2373 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2374 logger.info(
2375 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2377 for i in range(3):
2378 try:
2379 return parse_leap_seconds_list(fn)
2381 except LeapSecondsOutdated:
2382 try:
2383 logger.info('updating leap seconds list...')
2384 download_file(url, fn)
2386 except Exception as e:
2387 raise LeapSecondsError(
2388 'cannot download leap seconds list from %s to %s (%s)'
2389 % (url, fn, e))
2391 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2394def gps_utc_offset(t_utc):
2395 '''
2396 Time offset t_gps - t_utc for a given t_utc.
2397 '''
2398 ls = read_leap_seconds2()
2399 i = 0
2400 if t_utc < ls[0][0]:
2401 return ls[0][1] - 1 - 9
2403 while i < len(ls) - 1:
2404 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2405 return ls[i][1] - 9
2406 i += 1
2408 return ls[-1][1] - 9
2411def utc_gps_offset(t_gps):
2412 '''
2413 Time offset t_utc - t_gps for a given t_gps.
2414 '''
2415 ls = read_leap_seconds2()
2417 if t_gps < ls[0][0] + ls[0][1] - 9:
2418 return - (ls[0][1] - 1 - 9)
2420 i = 0
2421 while i < len(ls) - 1:
2422 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2423 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2424 return - (ls[i][1] - 9)
2425 i += 1
2427 return - (ls[-1][1] - 9)
2430def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2431 import itertools
2432 import glob
2433 from pyrocko.io.io_common import FileLoadError
2435 def iload_filename(filename, **kwargs):
2436 try:
2437 with open(filename, 'rb') as f:
2438 for cr in iload_fh(f, **kwargs):
2439 yield cr
2441 except FileLoadError as e:
2442 e.set_context('filename', filename)
2443 raise
2445 def iload_dirname(dirname, **kwargs):
2446 for entry in os.listdir(dirname):
2447 fpath = op.join(dirname, entry)
2448 if op.isfile(fpath):
2449 for cr in iload_filename(fpath, **kwargs):
2450 yield cr
2452 def iload_glob(pattern, **kwargs):
2454 fns = glob.glob(pattern)
2455 for fn in fns:
2456 for cr in iload_filename(fn, **kwargs):
2457 yield cr
2459 def iload(source, **kwargs):
2460 if isinstance(source, str):
2461 if op.isdir(source):
2462 return iload_dirname(source, **kwargs)
2463 elif op.isfile(source):
2464 return iload_filename(source, **kwargs)
2465 else:
2466 return iload_glob(source, **kwargs)
2468 elif hasattr(source, 'read'):
2469 return iload_fh(source, **kwargs)
2470 else:
2471 return itertools.chain.from_iterable(
2472 iload(subsource, **kwargs) for subsource in source)
2474 iload_filename.__doc__ = '''
2475 Read %s information from named file.
2476 ''' % doc_fmt
2478 iload_dirname.__doc__ = '''
2479 Read %s information from directory of %s files.
2480 ''' % (doc_fmt, doc_fmt)
2482 iload_glob.__doc__ = '''
2483 Read %s information from files matching a glob pattern.
2484 ''' % doc_fmt
2486 iload.__doc__ = '''
2487 Load %s information from given source(s)
2489 The ``source`` can be specified as the name of a %s file, the name of a
2490 directory containing %s files, a glob pattern of %s files, an open
2491 filehandle or an iterator yielding any of the forementioned sources.
2493 This function behaves as a generator yielding %s objects.
2494 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2496 for f in iload_filename, iload_dirname, iload_glob, iload:
2497 f.__module__ = iload_fh.__module__
2499 return iload_filename, iload_dirname, iload_glob, iload
2502class Inconsistency(Exception):
2503 pass
2506def consistency_check(list_of_tuples, message='values differ:'):
2507 '''
2508 Check for inconsistencies.
2510 Given a list of tuples, check that all tuple elements except for first one
2511 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2512 valid because the coordinates at the two channels are the same.
2513 '''
2515 if len(list_of_tuples) >= 2:
2516 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2517 raise Inconsistency('%s\n' % message + '\n'.join(
2518 ' %s: %s' % (t[0], ', '.join('%g' % x for x in t[1:]))
2519 for t in list_of_tuples))
2522class defaultzerodict(dict):
2523 def __missing__(self, k):
2524 return 0
2527def mostfrequent(x):
2528 c = defaultzerodict()
2529 for e in x:
2530 c[e] += 1
2532 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2535def consistency_merge(list_of_tuples,
2536 message='values differ:',
2537 error='raise',
2538 merge=mostfrequent):
2540 assert error in ('raise', 'warn', 'ignore')
2542 if len(list_of_tuples) == 0:
2543 raise Exception('cannot merge empty sequence')
2545 try:
2546 consistency_check(list_of_tuples, message)
2547 return list_of_tuples[0][1:]
2548 except Inconsistency as e:
2549 if error == 'raise':
2550 raise
2552 elif error == 'warn':
2553 logger.warning(str(e))
2555 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2558def parse_md(f):
2559 try:
2560 with open(op.join(
2561 op.dirname(op.abspath(f)),
2562 'README.md'), 'r') as readme:
2563 mdstr = readme.read()
2564 except IOError as e:
2565 return 'Failed to get README.md: %s' % e
2567 # Remve the title
2568 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2569 # Append sphinx reference to `pyrocko.` modules
2570 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2571 # Convert Subsections to toc-less rubrics
2572 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2573 return mdstr
2576def mpl_show(plt):
2577 import matplotlib
2578 if matplotlib.get_backend().lower() == 'agg':
2579 logger.warning('Cannot show() when using matplotlib "agg" backend')
2580 else:
2581 plt.show()
2584g_re_qsplit = re.compile(
2585 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2588g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2589g_re_escape_s = re.compile(r'([\\\'])')
2590g_re_unescape_s = re.compile(r'\\([\\\'])')
2591g_re_escape_d = re.compile(r'([\\"])')
2592g_re_unescape_d = re.compile(r'\\([\\"])')
2595def escape_s(s):
2596 '''
2597 Backslash-escape single-quotes and backslashes.
2599 Example: ``Jack's`` => ``Jack\\'s``
2601 '''
2602 return g_re_escape_s.sub(r'\\\1', s)
2605def unescape_s(s):
2606 '''
2607 Unescape backslash-escaped single-quotes and backslashes.
2609 Example: ``Jack\\'s`` => ``Jack's``
2610 '''
2611 return g_re_unescape_s.sub(r'\1', s)
2614def escape_d(s):
2615 '''
2616 Backslash-escape double-quotes and backslashes.
2618 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2619 '''
2620 return g_re_escape_d.sub(r'\\\1', s)
2623def unescape_d(s):
2624 '''
2625 Unescape backslash-escaped double-quotes and backslashes.
2627 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2628 '''
2629 return g_re_unescape_d.sub(r'\1', s)
2632def qjoin_s(it):
2633 '''
2634 Join sequence of strings into a line, single-quoting non-trivial strings.
2636 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2637 '''
2638 return ' '.join(
2639 w if g_re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2642def qjoin_d(it):
2643 '''
2644 Join sequence of strings into a line, double-quoting non-trivial strings.
2646 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2647 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2648 '''
2649 return ' '.join(
2650 w if g_re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2653def qsplit(s):
2654 '''
2655 Split line into list of strings, allowing for quoted strings.
2657 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2658 ``["55", "Sparrow's Island"]``,
2659 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2660 ``['55', 'Pete "The Robot" Smith']``
2661 '''
2662 return [
2663 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2664 for x in g_re_qsplit.findall(s)]