Source code for spacepy.toolbox

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Toolbox of various functions and generic utilities.

Authors: Steve Morley, Jon Niehof, Brian Larsen, Josef Koller, Dan Welling
Institution: Los Alamos National Laboratory
Contact: smorley@lanl.gov, jniehof@lanl.gov, balarsen@lanl.gov, jkoller@lanl.gov, dwelling@lanl.gov
Los Alamos National Laboratory

Copyright 2010 Los Alamos National Security, LLC.
"""

#If you add functions here, be sure to:
#1) add to the __all__ list
#2) add to functions in Doc/source/toolbox.rst so it goes in the docs
from __future__ import absolute_import
from __future__ import division

import calendar
import datetime
import glob
try:
    import html.parser
except ImportError: #python2
    import HTMLParser as html
    html.parser = html
try:
    import http.client
except ImportError: # Python2
    import httplib as http
    http.client = http
import numbers
import os
import os.path
import re
try:
    import urllib.request
except ImportError: #python2
    import urllib2 as urllib
    urllib.request = urllib
import select
import shutil
import socket
import subprocess
import sys
import tempfile
import time
import warnings
import zipfile

try:
    import cPickle as pickle
except:
    import pickle

import numpy as np

try:
    from spacepy import help
except ImportError:
    pass
except:
    pass

import spacepy
from spacepy import time as spt

#Try to pull in the C version. Assumption is that if you import this module,
#you want to do some association analysis, so the overhead in the import
#is OK.
from spacepy import lib
if lib.have_libspacepy:
    import ctypes

#Py3k compatibility renamings
try:
    xrange
except NameError:
    xrange = range

__all__ = ['tOverlap', 'tOverlapHalf', 'tCommon', 'loadpickle', 'savepickle', 'assemble',
           'human_sort', 'feq', 'dictree', 'update', 'progressbar',
           'windowMean', 'medAbsDev', 'binHisto', 'bootHisto',
           'logspace', 'geomspace', 'linspace', 'arraybin', 'mlt2rad',
           'rad2mlt', 'pmm', 'getNamedPath', 'query_yes_no',
           'quaternionNormalize', 'quaternionRotateVector', 'quaternionMultiply',
           'quaternionConjugate',
           'interpol', 'normalize', 'intsolve', 'dist_to_list',
           'bin_center_to_edges', 'bin_edges_to_center', 'thread_job', 'thread_map',
           'eventTimer', 'isview', 'interweave', 'indsFromXrange', 'hypot',
           'do_with_timeout', 'TimeoutError', 'timeout_check_call',
           'poisson_fit', 'unique_columns']

__contact__ = 'Brian Larsen: balarsen@lanl.gov'

def unique_columns(inval, axis=0):
    """
    Given a multidimensional input return the unique rows or columns along the given axis.
    Based largely on http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array
    axis=0 is unique rows, axis=1 is unique columns

    Parameters
    ==========
    inval :  array-like
        array to find unique columns or rows of

    Optional Parameters
    ===================
    axis : int
        The axis to find unique over, default: 0

    Returns
    =======
    out : array
        N-dimensional array of the unique values along the axis

    Examples
    ========
    """
    # this is a nice trick taking advantage of structed arrays where each row or column
    #   is the value, so returl unique works
    # np.ascontiguousarray() is to be really sure it will work
    if axis == 0:
        val = np.ascontiguousarray(np.transpose(inval))
    else:
        val = np.ascontiguousarray(inval)
    b = val.view(np.dtype((np.void, val.dtype.itemsize * val.shape[1])))
    unique_a = np.unique(b).view(val.dtype).reshape(-1, val.shape[1])
    return unique_a

    
[docs]def hypot(*args): """ compute the N-dimensional hypot of an iterable or many arguments Parameters ========== args : many numbers or array-like array like or many inputs to compute from Returns ======= out : float N-dimensional hypot of a number Notes ===== This function has a complicated speed function. - if a numpy array of floats is input this is passed off to C - if iterables are passed in they are made into numpy arrays and comptaton is done local - if many scalar agruments are passed in calculation is done in a loop For max speed: - <20 elements expand them into scalars >>> tb.hypot(*vals) >>> tb.hypot(vals[0], vals[1]...) #alternate - >20 elements premake them into a numpy array of doubles Examples ======== >>> from spacepy import toolbox as tb >>> print tb.hypot([3,4]) 5.0 >>> print tb.hypot(3,4) 5.0 >>> # Benchmark #### >>> from spacepy import toolbox as tb >>> import numpy as np >>> import timeit >>> num_list = [] >>> num_np = [] >>> num_np_double = [] >>> num_scalar = [] >>> tot = 500 >>> for num in tb.logspace(1, tot, 10): >>> print num >>> num_list.append(timeit.timeit(stmt='tb.hypot(a)', setup='from spacepy import toolbox as tb; import numpy as np; a = [3]*{0}'.format(int(num)), number=10000)) >>> num_np.append(timeit.timeit(stmt='tb.hypot(a)', setup='from spacepy import toolbox as tb; import numpy as np; a = np.asarray([3]*{0})'.format(int(num)), number=10000)) >>> num_scalar.append(timeit.timeit(stmt='tb.hypot(*a)', setup='from spacepy import toolbox as tb; import numpy as np; a = [3]*{0}'.format(int(num)), number=10000)) >>> from pylab import * >>> loglog(tb.logspace(1, tot, 10), num_list, lw=2, label='list') >>> loglog(tb.logspace(1, tot, 10), num_np, lw=2, label='numpy->ctypes') >>> loglog(tb.logspace(1, tot, 10), num_scalar, lw=2, label='scalar') >>> legend(shadow=True, fancybox=1, loc='upper left') >>> title('Different hypot times for 10000 runs') >>> ylabel('Time [s]') >>> xlabel('Size') .. image:: ../../source/images/hypot_no_extension_speeds_3cases.png """ if lib.have_libspacepy: if len(args) == 1 and isinstance(args[0], np.ndarray): # it is an array # make sure everything is C-ready ans = lib.hypot_tb( np.require(args[0], dtype=np.double, requirements='C'), np.product(args[0].shape)) return ans ans = 0.0 for arg in args: if hasattr(arg, '__iter__'): tmp = np.asanyarray(arg) ans += np.sum(tmp**2) else: ans += arg**2 return np.sqrt(ans)
[docs]def tOverlap(ts1, ts2, *args, **kwargs): """ Finds the overlapping elements in two lists of datetime objects Parameters ========== ts1 : datetime first set of datetime object ts2 : datetime datatime object args : additional arguments passed to tOverlapHalf Returns ======= out : list indices of ts1 within interval of ts2, & vice versa Examples ======== Given two series of datetime objects, event_dates and omni['Time']: >>> import spacepy.toolbox as tb >>> from spacepy import omni >>> import datetime >>> event_dates = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 10, 1), deltadays=3) >>> onni_dates = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 10, 1), deltadays=0.5) >>> omni = omni.get_omni(onni_dates) >>> [einds,oinds] = tb.tOverlap(event_dates, omni['ticks']) >>> omni_time = omni['ticks'][oinds[0]:oinds[-1]+1] >>> print omni_time [datetime.datetime(2000, 1, 1, 0, 0), datetime.datetime(2000, 1, 1, 12, 0), ... , datetime.datetime(2000, 9, 30, 0, 0)] See Also ======== tOverlapHalf tCommon """ idx_1in2 = tOverlapHalf(ts2, ts1, *args, **kwargs) idx_2in1 = tOverlapHalf(ts1, ts2, *args, **kwargs) if len(idx_2in1) == 0: idx_2in1 = None if len(idx_1in2) == 0: idx_1in2 = None return idx_1in2, idx_2in1
[docs]def tOverlapHalf(ts1, ts2, presort=False): """ Find overlapping elements in two lists of datetime objects This is one-half of tOverlap, i.e. it finds only occurrences where ts2 exists within the bounds of ts1, or the second element returned by tOverlap. Parameters ========== ts1 : list first set of datetime object ts2 : list datatime object presort : bool Set to use a faster algorithm which assumes ts1 and ts2 are both sorted in ascending order. This speeds up the overlap comparison by about 50x, so it is worth sorting the list if one sort can be done for many calls to tOverlap Returns ======= out : list indices of ts2 within interval of ts1 **note:** Returns empty list if no overlap found See Also ======== tOverlap tCommon """ if presort: import bisect t_lower, t_upper = ts1[0], ts1[-1] return xrange(bisect.bisect_left(ts2, t_lower), bisect.bisect_right(ts2, t_upper)) else: t_lower, t_upper = min(ts1), max(ts1) return [i for i in range(len(ts2)) if ts2[i] >= t_lower and ts2[i] <= t_upper]
[docs]def tCommon(ts1, ts2, mask_only=True): """ Finds the elements in a list of datetime objects present in another Parameters ========== ts1 : list or array-like first set of datetime objects ts2 : list or array-like second set of datetime objects Returns ======= out : tuple Two element tuple of truth tables (of 1 present in 2, & vice versa) See Also ======== tOverlapHalf tOverlap Examples ======== >>> import spacepy.toolbox as tb >>> import numpy as np >>> import datetime as dt >>> ts1 = np.array([dt.datetime(2001,3,10)+dt.timedelta(hours=a) for a in range(20)]) >>> ts2 = np.array([dt.datetime(2001,3,10,2)+dt.timedelta(hours=a*0.5) for a in range(20)]) >>> common_inds = tb.tCommon(ts1, ts2) >>> common_inds[0] #mask of values in ts1 common with ts2 array([False, False, True, True, True, True, True, True, True, True, True, True, False, False, False, False, False, False, False, False], dtype=bool) >>> ts2[common_inds[1]] #values of ts2 also in ts1 The latter can be found more simply by setting the mask_only keyword to False >>> common_vals = tb.tCommon(ts1, ts2, mask_only=False) >>> common_vals[1] array([2001-03-10 02:00:00, 2001-03-10 03:00:00, 2001-03-10 04:00:00, 2001-03-10 05:00:00, 2001-03-10 06:00:00, 2001-03-10 07:00:00, 2001-03-10 08:00:00, 2001-03-10 09:00:00, 2001-03-10 10:00:00, 2001-03-10 11:00:00], dtype=object) """ from matplotlib.dates import date2num, num2date tn1, tn2 = date2num(ts1), date2num(ts2) el1in2 = np.in1d(tn1, tn2, assume_unique=True) #makes mask of present/absent el1in2 = np.in1d(tn1, tn2, assume_unique=True) #makes mask of present/absent el2in1 = np.in1d(tn2, tn1, assume_unique=True) if mask_only: return el1in2, el2in1 else: truemask1 = np.abs(np.array(el1in2)-1) truemask2 = np.abs(np.array(el2in1)-1) time1 = np.ma.masked_array(tn1, mask=truemask1) time2 = np.ma.masked_array(tn2, mask=truemask2) dum1 = num2date(time1.compressed()) dum2 = num2date(time2.compressed()) # dum1 = [val.replace(tzinfo=None) for val in dum1] # this hits ValueError: microsecond must be in 0..999999 # dum2 = [val.replace(tzinfo=None) for val in dum2] dum11 = [] dum22 = [] for v1, v2 in zip(dum1, dum2): try: dum11.append(v1.replace(tzinfo=None)) except ValueError: # ValueError: microsecond must be in 0..999999 dum11.append((datetime.datetime(v1.year, v1.month, v1.day, v1.hour, v1.minute, v1.second, 0) + datetime.timedelta(seconds=1)).replace(tzinfo=None)) try: dum22.append(v2.replace(tzinfo=None)) except ValueError: # ValueError: microsecond must be in 0..999999 dum22.append((datetime.datetime(v2.year, v2.month, v2.day, v2.hour, v2.minute, v2.second, 0) + datetime.timedelta(seconds=1)).replace(tzinfo=None)) dum1 = dum11 dum2 = dum22 if type(ts1)==np.ndarray or type(ts2)==np.ndarray: dum1 = np.array(dum1) dum2 = np.array(dum2) return dum1, dum2
[docs]def loadpickle(fln): """ load a pickle and return content as dictionary Parameters ========== fln : string filename Returns ======= out : dict dictionary with content from file See Also ======== savepickle Examples ======== **note**: If fln is not found, but the same filename with '.gz' is found, will attempt to open the .gz as a gzipped file. >>> d = loadpickle('test.pbin') """ if not os.path.exists(fln) and os.path.exists(fln + '.gz'): gzip = True fln += '.gz' else: try: with open(fln, 'rb') as fh: try: #Py3k return pickle.load(fh, encoding='latin1') except TypeError: return pickle.load(fh) except pickle.UnpicklingError: #maybe it's a gzip? gzip = True else: gzip = False if gzip: try: import zlib with open(fln, 'rb') as fh: stream = zlib.decompress(fh.read(), 16 + zlib.MAX_WBITS) try: #Py3k return pickle.loads(stream, encoding='latin1') except TypeError: return pickle.loads(stream) except MemoryError: import gzip with open(fln) as fh: gzh = gzip.GzipFile(fileobj=fh) try: #Py3k contents = pickle.load(gzh, encoding='latin1') except TypeError: contents = pickle.load(gzh) gzh.close() return contents
# -----------------------------------------------
[docs]def savepickle(fln, dict, compress=None): """ save dictionary variable dict to a pickle with filename fln Parameters ---------- fln : string filename dict : dict container with stuff compress : bool write as a gzip-compressed file (.gz will be added to ``fln``). If not specified, defaults to uncompressed, unless the compressed file exists and the uncompressed does not. See Also ======== loadpickle Examples ======== >>> d = {'grade':[1,2,3], 'name':['Mary', 'John', 'Chris']} >>> savepickle('test.pbin', d) """ if compress == None: # Guess at compression # Assume compressed if compressed already exists (and no uncompressed) compress = not os.path.exists(fln) and os.path.exists(fln + '.gz') if compress: import gzip with open(fln + '.gz', 'wb') as fh: gzh = gzip.GzipFile(fln, 'wb', compresslevel=3, fileobj=fh) pickle.dump(dict, gzh, 2) gzh.close() else: with open(fln, 'wb') as fh: pickle.dump(dict, fh, 2) # 2 ... fast binary
# -----------------------------------------------
[docs]def assemble(fln_pattern, outfln, sortkey='ticks', verbose=True): """ assembles all pickled files matching fln_pattern into single file and save as outfln. Pattern may contain simple shell-style wildcards \*? a la fnmatch file will be assembled along time axis given by Ticktock (key: 'ticks') in dictionary If sortkey = None, then nothing will be sorted Parameters ---------- fln_pattern : string pattern to match filenames outfln : string filename to save combined files to Returns ======= out : dict dictionary with combined values Examples ======== >>> import spacepy.toolbox as tb >>> a, b, c = {'ticks':[1,2,3]}, {'ticks':[4,5,6]}, {'ticks':[7,8,9]} >>> tb.savepickle('input_files_2001.pkl', a) >>> tb.savepickle('input_files_2002.pkl', b) >>> tb.savepickle('input_files_2004.pkl', c) >>> a = tb.assemble('input_files_*.pkl', 'combined_input.pkl') ('adding ', 'input_files_2001.pkl') ('adding ', 'input_files_2002.pkl') ('adding ', 'input_files_2004.pkl') ('\\n writing: ', 'combined_input.pkl') >>> print(a) {'ticks': array([1, 2, 3, 4, 5, 6, 7, 8, 9])} """ # done this way so it works before install from spacepy import coordinates as c filelist = glob.glob(fln_pattern) filelist = human_sort(filelist) # read all files d = {} for fln in filelist: if verbose: print("adding ", fln) d[fln] = loadpickle(fln) # combine them dcomb = d[filelist[0]] # copy the first file over for fln in filelist[1:]: # check if sortkey is actually available assert (sortkey in d[fln] or sortkey==None), 'provided sortkey ='+sortkey+' is not available' if sortkey: TAIcount = len(d[fln][sortkey]) else: TAIcount = len(d[fln][ list(d[fln].keys())[0] ]) for key in d[fln]: #print fln, key dim = np.array(np.shape(d[fln][key])) ax = np.where(dim==TAIcount)[0] if len(ax) == 1: # then match with TAI length is given (jump over otherwise like for 'parameters') if isinstance(dcomb[key], spt.Ticktock): dcomb[key] = dcomb[key].append(d[fln][key]) elif isinstance(dcomb[key], c.Coords): dcomb[key] = dcomb[key].append(d[fln][key]) else: dcomb[key] = np.append(dcomb[key], d[fln][key], axis=ax[0]) if sortkey: # then sort if isinstance(dcomb[sortkey], spt.Ticktock): idx = np.argsort(dcomb[sortkey].RDT) else: idx = np.argsort(dcomb[sortkey]) TAIcount = len(dcomb[sortkey]) for key in dcomb: # iterates over keys by default dim = np.array(np.shape(dcomb[key])) ax = np.where(dim==TAIcount)[0] if len(ax) == 1: # then match with length of TAI dcomb[key] = dcomb[key][idx] # resort else: # do nothing pass if verbose: print('\n writing: ', outfln) savepickle(outfln, dcomb) return dcomb
[docs]def human_sort( l ): """ Sort the given list in the way that humans expect. http://www.codinghorror.com/blog/2007/12/sorting-for-humans-natural-sort-order.html Parameters ---------- l : list list of objects to human sort Returns ------- out : list sorted list Examples -------- >>> import spacepy.toolbox as tb >>> dat = ['r1.txt', 'r10.txt', 'r2.txt'] >>> dat.sort() >>> print dat ['r1.txt', 'r10.txt', 'r2.txt'] >>> tb.human_sort(dat) ['r1.txt', 'r2.txt', 'r10.txt'] """ convert = lambda text: int(text) if text.isdigit() else text alphanum_key = lambda key: [ convert(c) for c in re.split(r'([0-9]+)', key) ] try: l.sort( key=alphanum_key ) except TypeError: l.sort() return l
[docs]@spacepy.deprecated( '0.2.2', 'use numpy.isclose', 'use :func:`numpy.isclose`; the ``rtol`` keyword is equivalent to' ' ``precision``.') def feq(x, y, precision=0.0000005): """ compare two floating point values if they are equal after: http://www.lahey.com/float.htm See Also -------- numpy.allclose numpy.isclose Parameters ---------- x : float a number y : float or array of floats other numbers to compare precision : float (optional) Relative precision for equal (default 0.0000005) Specified as a fraction of the sum of ``x`` and ``y``. Returns ------- out : bool True (equal) or False (not equal) Examples -------- >>> import spacepy.toolbox as tb >>> x = 1 + 1e-4 >>> y = 1 + 2e-4 >>> tb.feq(x, y) False >>> tb.feq(x, y, 1e-3) True """ x = np.asanyarray(x) y = np.asanyarray(y) boolean = abs(x-y) <= (abs(x+y)*precision) return boolean
[docs]def dictree(in_dict, verbose=False, spaces=None, levels=True, attrs=False, **kwargs): """ pretty print a dictionary tree Parameters ---------- in_dict : dict a complex dictionary (with substructures) verbose : boolean (optional) print more info spaces : string (optional) string will added for every line levels : integer (optional) number of levels to recurse through (True means all) attrs : boolean (optional) display information for attributes Examples -------- >>> import spacepy.toolbox as tb >>> d = {'grade':{'level1':[4,5,6], 'level2':[2,3,4]}, 'name':['Mary', 'John', 'Chris']} >>> tb.dictree(d) + |____grade |____level1 |____level2 |____name More complicated example using a datamodel: >>> from spacepy import datamodel >>> counts = datamodel.dmarray([2,4,6], attrs={'units': 'cts/s'}) >>> data = {'counts': counts, 'PI': 'Dr Zog'} >>> tb.dictree(data) + |____PI |____counts >>> tb.dictree(data, attrs=True, verbose=True) + |____PI (str [6]) |____counts (spacepy.datamodel.dmarray (3,)) :|____units (str [5]) Attributes of, e.g., a CDF or a datamodel type object (obj.attrs) are denoted by a colon. """ try: assert hasattr(in_dict, 'keys') except AssertionError: try: assert hasattr(in_dict, 'attrs') except: raise TypeError('dictree: Input must be dictionary-like') if not spaces: spaces = '' print('+') if 'toplev' in kwargs: toplev = kwargs['toplev'] else: toplev = True try: if toplev and attrs: dictree(in_dict.attrs, spaces = ':', verbose = verbose, levels = levels, attrs=attrs, toplev=True) toplev = False except: pass # TODO, if levels is True why check again? if levels: try: assert levels is True except AssertionError: levels -= 1 if levels == 0: levels = None try: for key in sorted(in_dict.keys()): bar = '|____' + str(key) if verbose: typestr = str(type(in_dict[key])).split("'")[1] #check entry for dict-like OR .attrs dict try: dimstr = in_dict[key].shape dimstr = ' ' + str(dimstr) except AttributeError: try: dimstr = len(in_dict[key]) dimstr = ' [' + str(dimstr) + ']' except: dimstr = '' print(spaces + bar + ' ('+ typestr + dimstr + ')') else: print(spaces + bar) if hasattr(in_dict[key], 'attrs') and attrs: dictree(in_dict[key].attrs, spaces = spaces + ' :', verbose = verbose, levels = levels, attrs=attrs, toplev=False) if hasattr(in_dict[key], 'keys') and levels: dictree(in_dict[key], spaces = spaces + ' ', verbose = verbose, levels = levels, attrs=attrs, toplev=False) except: pass return None
def _crawl_yearly(base_url, pattern, datadir, name=None, cached=True, startyear=None): """Crawl files in a directory-by-year structure Parameters ========== base_url : str Base of the data. This URL should point to a directory containing yearly directories (YYYY). pattern : str Regular expression to match filenames. Will download files in each yearly directory that match the pattern. datadir : str Directory to store downloaded files. A mirror will be maintained in this directory. Note this is a "flat" mirror without the year directories. name : str (optional) Name of the data set, used only in status messages. cached : boolean (optional) Only update files if timestamp on server is newer than timestamp on local file (default). Set False to always download files. startyear : int (optional) First year to crawl, as four-digit number or four-character string. If not specified, will download all years. If specified, will delete years prior to this which have already been downloaded. Returns ======= list All the filenames that were mirrored, in order; or None if there were no updates. If there are any updates, all filenames are included. """ name = '' if name is None else '{} '.format(name) #Find all the files to download print("Finding {}files to download ...".format(name)) progressbar(0, 1, 1, text='Listing files') conn = None if spacepy.config['keepalive']: try: data, conn = get_url(base_url, keepalive=True) except socket.error: #Give up on keepalives pass if conn is None: data = get_url(base_url) if str is not bytes: data = data.decode('utf-8') p = LinkExtracter() p.feed(data) p.close() yearlist = [y[0:4] for y in p.links if re.match(r'\d{4}/', y) and (startyear is None or y[0:4] >= str(startyear))] downloadme = {} for i, y in enumerate(yearlist): yearurl = '{}{}/'.format(base_url, y) if conn is None: data = get_url(yearurl) else: data, conn = get_url(yearurl, keepalive=True, conn=conn) if str is not bytes: data = data.decode('utf-8') p = LinkExtracter() p.feed(data) p.close() for f in p.links: if not re.match(pattern, f): continue downloadme[f] = yearurl + f progressbar(i + 1, 1, len(yearlist), text='Listing files') print("Retrieving {}files ...".format(name)) filenames = sorted(list(downloadme.keys())) if not os.path.exists(datadir): os.makedirs(datadir) newdata = False #Check for existing files (delete them if no longer on server) have_files = os.listdir(datadir) for f in have_files: if not f in filenames: os.remove(os.path.join(datadir, f)) #File was removed, so need to reparse even if no new downloads newdata = True #Download for i, fname in enumerate(filenames): if conn is None: res = get_url(downloadme[fname], os.path.join(datadir, fname), cached=cached) else: res, conn = get_url(downloadme[fname], os.path.join(datadir, fname), cached=cached, keepalive=True, conn=conn) if res is not None: newdata = True progressbar(i + 1, 1, len(filenames)) if conn is not None: conn.close() return filenames if newdata else None def _get_qindenton_daily(qd_daily_url=None, cached=True, startyear=None): """Download the Qin-Denton OMNI-like daily files Parameters ========== qd_daily_url : str (optional) Base of the Qin-Denton data, in hourly JSON-headed ASCII. This URL should point to the directory containing the yearly directories. Default from ``qd_daily_url`` in config file. cached : boolean (optional) Only update files if timestamp on server is newer than timestamp on local file (default). Set False to always download files. startyear : int (optional) If specified, start downloading files from the year given, rather than all years. This will delete older files! Returns ======= SpaceData The data extracted from the Q-D dataset, fully processed for saving as SpacePy HDF5 OMNI data. """ if qd_daily_url is None: qd_daily_url = spacepy.config['qd_daily_url'] datadir = os.path.join(spacepy.DOT_FLN, 'data', 'qindenton_daily_files') _crawl_yearly(qd_daily_url, r'QinDenton_\d{8}_hour.txt', datadir, name='Q-D daily', cached=cached, startyear=startyear) #Read and process print("Processing Q-D daily files ...") return _assemble_qindenton_daily(datadir) def _assemble_qindenton_daily(qd_daily_dir): """Assemble Qin-Denton daily files into OMNI structure Parameters ========== qd_daily_dir : str Directory with Qin-Denton daily files. Returns ======= SpaceData The data extracted from the Q-D dataset, fully processed for saving as SpacePy HDF5 OMNI data. """ import spacepy.datamodel filelist = sorted(glob.glob(os.path.join(qd_daily_dir, '*_hour.txt'))) data = [spacepy.datamodel.readJSONheadedASCII(f) for f in filelist] omnidata = spacepy.datamodel.SpaceData() for k in data[0].keys(): if k in ('DateTime', 'Minute', 'OriginFile', 'Second'): continue omnidata[k] = spacepy.datamodel.dmarray( np.concatenate([d[k] for d in data]), dtype=np.float32) del data ntimes = set([len(v) for v in omnidata.values()]) if len(ntimes) != 1: raise ValueError( 'Input Q-D daily file has different size for different variables') ntimes = ntimes.pop() # Renaming from the names in new file to old file names for oldname, newname in (('dens', 'Den_P'), ('velo', 'Vsw')): omnidata[oldname] = omnidata[newname] del omnidata[newname] if '{}_status'.format(newname) in omnidata: omnidata['{}_status'.format(oldname)] \ = omnidata['{}_status'.format(newname)] del omnidata['{}_status'.format(newname)] # Quality flags qbits = spacepy.datamodel.SpaceData() for k in list(omnidata.keys()): # Edit while iterate if not k.endswith('_status'): continue v = spacepy.datamodel.dmarray(omnidata[k], dtype=np.int8) basename = k.split('_')[0] if basename in ('G', 'W'): # Arrays for i in range(omnidata[k].shape[1]): qbits['{}{:d}'.format(basename, i + 1)] = v[:, i] else: qbits[basename] = v del omnidata[k] omnidata['Qbits'] = qbits # Reformat some arrays to multi-variable for name in ('Bz', 'G', 'W'): for i in range(omnidata[name].shape[1]): omnidata['{}{:d}'.format(name, i + 1)] = omnidata[name][:, i] del omnidata[name] # Make integers of integers for k in ('Dst', 'Year', 'Month', 'Day', 'Hour'): omnidata[k] = spacepy.datamodel.dmarray(omnidata[k], dtype=np.int16) # Process time formats omnidata['UTC'] = spacepy.datamodel.dmarray([ datetime.datetime(omnidata['Year'][i], omnidata['Month'][i], omnidata['Day'][i], omnidata['Hour'][i]) for i in range(len(omnidata['Year']))]) omnidata['DOY'] = spacepy.datamodel.dmarray([ dt.timetuple().tm_yday for dt in omnidata['UTC']], dtype=np.int16) omnidata['RDT'] = spt.Ticktock(omnidata['UTC'], 'UTC').RDT for k in ('Year', 'Hour', 'Month', 'Day'): del omnidata[k] return omnidata def _get_cdaweb_omni2(omni2url=None): """Download the OMNI2 data from SPDF Parameters ========== omni2url : str (optional) Base of the OMNI2 data at SPDF, in hourly CDF form. This URL should point to the directory containing the yearly directories. Default from ``omni2_url`` in config file. Returns ======= SpaceData The data extracted from the OMNI2 dataset, with variables renamed to match the old ViRBO combined OMNI2 CDF. Returns ``None`` if there are no new data. """ import spacepy.pycdf import spacepy.pycdf.istp if omni2url is None: omni2url = spacepy.config['omni2_url'] datadir = os.path.join(spacepy.DOT_FLN, 'data', 'omni2cdfs') filenames = _crawl_yearly(omni2url, r'omni2_h0_mrg1hr_\d{8}_v\d+\.cdf', datadir, name='OMNI2') if filenames is None: return None #Read and process print("Reading OMNI2 files ...") data = spacepy.pycdf.concatCDF([spacepy.pycdf.CDF(os.path.join(datadir, f)) for f in filenames]) #Map keyed by the variable names as in the OMNI2 data from CDAWeb, #valued by the variable names as in the OMNI2 data from ViRBO #i.e. this is from: to #None means not in VirBO, so delete. keymap = { 'ABS_B': 'Field_mag_ave', #This is average of the magnitude 'AE': 'AE_index', 'AL_INDEX': 'AL_index', 'AP_INDEX': 'Ap_index', 'AU_INDEX': 'AU_index', 'BX_GSE': 'Bx_GSE', 'BY_GSE': 'By_GSE', 'BZ_GSE': 'Bz_GSE', 'BY_GSM': 'By_GSM', 'BZ_GSM': 'Bz_GSM', 'Beta': 'Plasma_beta', 'Day': 'Decimal_Day', 'DST': 'Dst_index', 'E': 'Electric_field_GSM', #'Epoch': 'Epoch', #same, do not change 'F': 'Mag_of_ave_field_vect', #Magnitude of average 'F10_INDEX': 'f10_7_index', 'HR': 'Hour', 'IMF': 'IMF_spacecraft_ID', 'IMF_PTS': 'Num_pts_in_IMF_aves', 'KP': 'Kp_index', 'Mach_num': 'Alfven_mach_number', 'MFLX': 'Proton_flux_flag', 'Mgs_mach_num': None, #magnetosonic mach number 'N': 'Ion_density', 'PC_N_INDEX': 'PC_N_index', 'PHI-V': 'Plasma_bulk_flow_long_angle', 'PHI_AV': 'Long_angle_ave_field_vector', 'PLS': 'SW_plasma_spacecraft_ID', 'PLS_PTS': 'Num_pts_in_plasma_aves', 'PR-FLX_1': 'Proton_flux_gt_1_MeV', 'PR-FLX_10': 'Proton_flux_gt_10_MeV', 'PR-FLX_2': 'Proton_flux_gt_2_MeV', 'PR-FLX_30': 'Proton_flux_gt_30_MeV', 'PR-FLX_4': 'Proton_flux_gt_4_MeV', 'PR-FLX_60': 'Proton_flux_gt_60_MeV', 'Pressure': 'Flow_pressure', 'R': 'Sunspot_number', 'Ratio': 'Na_over_Np', 'Rot#': 'Bartels_rotation_num', 'SIGMA-ABS_B': 'sigma_field_mag_ave', 'SIGMA-B': 'sigma_mag_of_ave_field_vect', 'SIGMA-Bx': 'sigma_Bx_GSE', 'SIGMA-By': 'sigma_By_GSE', 'SIGMA-Bz': 'sigma_Bz_GSE', 'SIGMA-N': 'sigma_ion_density', 'SIGMA-PHI-V': 'sigma_plasma_flow_long_angle', 'SIGMA-T': 'sigma_temp', 'SIGMA-THETA-V': 'sigma_plasma_flow_lat_angle', 'SIGMA-V': 'sigma_plasma_bulk_speed', 'SIGMA-ratio': 'sigma_ratio', 'Solar_Lyman_alpha': None, 'T': 'Plasma_temp', 'THETA-V': 'Plasma_bulk_flow_lat_angle', 'THETA_AV': 'Lat_angle_ave_field_vector', 'V': 'Plasma_bulk_speed', 'YR': 'Year', } for k, v in keymap.items(): if v is not None: data[v] = data[k] del data[k] #The CDAWeb Epochs are erroneously tagged as START of collection, #not midpoint. So fix that data['Epoch'] = data['Epoch'] + datetime.timedelta(minutes=30) #Castings to match ViRBO: everything is an int, and fill turns to NaN for k, v in data.items(): if v.dtype == np.object: #Skip epoch continue if v.dtype == np.int32: data[k] = v.astype(np.float32) v = data[k] spacepy.pycdf.istp.nanfill(v) return data class LinkExtracter(html.parser.HTMLParser): """Finds all links in a HTML page, useful for crawling. After HTML has been parsed, the ``links`` attribute contains a list of link targets. """ def reset(self, *args, **kwargs): if isinstance(LinkExtracter, type): super(LinkExtracter, self).reset(*args, **kwargs) else: #py2k, old-style class html.parser.HTMLParser.reset(self, *args, **kwargs) self.links = [] """List of link targets found in the page""" def handle_starttag(self, tag, attrs): if tag != 'a': return for name, value in attrs: if name != 'href': continue self.links.append(value)
[docs]def get_url(url, outfile=None, reporthook=None, cached=False, keepalive=False, conn=None): """Read data from a URL Open an HTTP URL, honoring the user agent as specified in the SpacePy config file. Returns the data, optionally also writing out to a file. This is similar to the deprecated ``urlretrieve``. Parameters ========== url : str The URL to open outfile : str (optional) Full path to file to write data to reporthook : callable (optional) Function for reporting progress; takes arguments of block count, block size, and total size. cached : bool (optional) Compare modification time of the URL to the modification time of ``outfile``; do not retrieve (and return None) unless the URL is newer than the file. keepalive : bool (optional) Attempt to keep the connection open to retrieve more URLs. The return becomes a tuple of (data, conn) to return the connection used so it can be used again. This mode does not support proxies. (Default False) conn : http.client.HTTPConnection (optional) An established http connection (HTTPS is also okay) to use with ``keepalive``. If not provided, will attempt to make a connection. Returns ======= bytes The HTTP data from the server. See Also ======== progressbar Notes ===== This function honors proxy settings as described in :func:`urllib.request.getproxies`. Cryptic error messages (such as ``Network is unreachable``) may indicate that proxy settings should be defined as appropriate for your environment (e.g. with ``HTTP_PROXY`` or ``HTTPS_PROXY`` environment variables). """ if not keepalive and conn is not None: raise ValueError('Cannot specify connection without keepalive') if keepalive: scheme, _, host, path = url.split('/', 3) path = '/' + path # Explicitly root on the server if conn is not None and conn.sock is not None: readable, writeable, _ = select.select( [conn.sock], [conn.sock], [], 0) # Make sure no stale data to read on socket, and can write to it if readable or len(writeable) != 1: conn.close() conn = None if conn is None: ctype = http.client.HTTPConnection if scheme == 'http:'\ else http.client.HTTPSConnection conn = ctype(host) clheaders = { "Connection": "keep-alive", } if spacepy.config.get('user_agent', ''): clheaders['User-Agent'] = spacepy.config['user_agent'] conn.request('HEAD' if cached else 'GET', path, headers=clheaders) def checkresponse(conn): """Get the response on a connection, return response and headers""" r = conn.getresponse() if r.status >= 400: raise RuntimeError( 'HTTP status {} {}'.format(r.status, r.reason)) headers = dict(((k.title(), v) for k, v in r.getheaders())) return r, headers r, headers = checkresponse(conn) else: r = urllib.request.Request(url) if spacepy.config.get('user_agent', ''): r.add_header('User-Agent', spacepy.config['user_agent']) r = urllib.request.urlopen(r) if r.getcode() >= 400: r.close() raise RuntimeError('HTTP status {} {}'.format(r.code, r.msg)) headers = r.info() modified = headers.get('Last-Modified', None) if modified is not None: # strptime is affected by locale (including the month name) but the # header is a constant format, so massage modified = modified.split()[1:5] # Get rid of day of week and 'GMT' modified[1] = str(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']\ .index(modified[1]) + 1) # Make month numerical modified = datetime.datetime.strptime( ' '.join(modified), "%d %m %Y %H:%M:%S") modified = calendar.timegm(modified.timetuple()) if cached: if outfile is None: if not keepalive: r.close() raise RuntimeError('Must specify outfile if cached is True') if os.path.exists(outfile) and modified is not None: #Timestamp is truncated to second, so do same for local local_mod = int(os.path.getmtime(outfile)) if modified <= local_mod: if keepalive: r.read() else: r.close() return (None, conn) if keepalive else None if keepalive: # Replace previous header request with full get r.read() conn.request('GET', path, headers=clheaders) r, headers = checkresponse(conn) size = int(headers.get('Content-Length', 0)) blocksize = 1024 count = 0 data = [] while True: newdata = r.read(blocksize) if not newdata: break data.append(newdata) count += 1 if reporthook: reporthook(count, blocksize, size) if not keepalive: r.close() if outfile: with open(outfile, 'wb') as f: for d in data: f.write(d) if modified is not None: #Copy web mtime to file os.utime(outfile, (int(time.time()), modified)) data = b''.join(data) return (data, conn) if keepalive else data
[docs]def update(all=True, QDomni=False, omni=False, omni2=False, leapsecs=False, PSDdata=False, cached=True): """ Download and update local database for omni, leapsecs etc Web access is via :func:`get_url`; notes there may be helpful in debugging errors. See also the ``keepalive`` configuration option. Parameters ========== all : boolean (optional) if True, update OMNI2, Qin-Denton and leapsecs omni : boolean (optional) if True. update only omni (Qin-Denton) omni2 : boolean (optional) if True, update only original OMNI2 QDomni : boolean (optional) if True, update OMNI2 and Qin-Denton leapsecs : boolean (optional) if True, update only leapseconds cached : boolean (optional) Only update files if timestamp on server is newer than timestamp on local file (default). Set False to always download files. Returns ======= out : string data directory where things are saved See Also ======== get_url Examples ======== >>> import spacepy.toolbox as tb >>> tb.update(omni=True) """ from spacepy.datamodel import SpaceData, dmarray, fromCDF, toHDF5 from spacepy import DOT_FLN, config datadir = os.path.join(DOT_FLN, 'data') if not os.path.exists(datadir): os.mkdir(datadir) os.chmod(datadir, 0o777) #leapsec_url ='ftp://maia.usno.navy.mil/ser7/tai-utc.dat' leapsec_fname = os.path.join(datadir, 'tai-utc.dat') # define location for getting omni #omni_url = 'ftp://virbo.org/QinDenton/hour/merged/latest/WGhour-latest.d.zip' omni_fname_zip = os.path.join(datadir, 'WGhour-latest.d.zip') omni2_fname_zip = os.path.join(datadir, 'omni2-latest.cdf.zip') omni_fname_pkl = os.path.join(datadir, 'omnidata.pkl') omni_fname_json = os.path.join(datadir, 'omnidata.txt') omni_fname_h5 = os.path.join(datadir, 'omnidata.h5') omni2_fname_h5 = os.path.join(datadir, 'omni2data.h5') PSDdata_fname = os.path.join('psd_dat.sqlite') if (omni or omni2 or QDomni or leapsecs or PSDdata): all = False #if an option is explicitly selected, turn 'all' off if all == True: omni = True omni2 = True leapsecs = True if QDomni == True: omni = True omni2 = True if omni == True: # retrieve omni, unzip and save as table print("Retrieving initial Qin-Denton file ...") get_url(config['qindenton_url'], omni_fname_zip, progressbar, cached=cached) fh_zip = zipfile.ZipFile(omni_fname_zip) data = fh_zip.read(fh_zip.namelist()[0]) fh_zip.close() if not str is bytes: data = data.decode('ascii') A = np.array(data.split('\n')) print("Processing initial Qin-Denton file ...") # create a keylist keys = A[0].split() keys.remove('8') keys.remove('6') keys[keys.index('status')] = '8_status' keys[keys.index('stat')] = '6_status' keys[keys.index('dst')] = 'Dst' keys[keys.index('kp')] = 'Kp' #keys[keys.index('Hr')] = 'Hr' keys[keys.index('V_SW')] = 'velo' keys[keys.index('Den_P')] = 'dens' keys[keys.index('Day')] = 'DOY' keys[keys.index('Year')] = 'Year' # remove keyword lines and empty lines as well idx = np.where(A != '')[0] # put it into a 2D table tab = [val.split() for val in A[idx[1:]]] stat8 = [val[11] for val in tab] stat6 = [val[27] for val in tab] tab = np.array(tab, dtype='float32') # take out where Dst not available ( = 99999) or year == 0 idx = np.where((tab[:,12] !=99.0) & (tab[:,0] != 0))[0] tab = tab[idx,:] stat8 = np.array(stat8)[idx] stat6 = np.array(stat6)[idx] omnidata = SpaceData() # sort through and make an omni dictionary # extract keys from line above for ikey, i in zip(keys,range(len(keys))): if ikey in ('Year', 'DOY', 'Hr', 'Dst'): omnidata[ikey] = dmarray(tab[:, i], dtype='int16') else: omnidata[ikey] = dmarray(tab[:,i]) # add TAI to omnidata nTAI = len(omnidata['DOY']) # add interpolation quality flags omnidata['Qbits'] = SpaceData() arr = dmarray(stat8.view(stat8.dtype.kind + '1'), dtype=np.byte).reshape((8, nTAI)) for ik, key in enumerate(['ByIMF', 'BzIMF', 'velo', 'dens', 'Pdyn', 'G1', 'G2', 'G3']): omnidata['Qbits'][key] = arr[ik,:] if stat6.dtype.str[1:] == 'U6': stat6 = np.require(stat6, dtype='|S6') arr = dmarray(stat6.view(stat6.dtype.kind + '1'), dtype=np.byte).reshape((6, nTAI)) for ik, key in enumerate(['W1', 'W2', 'W3', 'W4', 'W5', 'W6']): omnidata['Qbits'][key] = arr[ik,:] #remove string status keys foo = omnidata.pop('6_status') foo = omnidata.pop('8_status') # add time information to omni pickle (long loop) omnidata['UTC'] = dmarray([datetime.datetime(int(omnidata['Year'][i]), 1, 1) + datetime.timedelta(days=int(omnidata['DOY'][i]) - 1, hours=int(omnidata['Hr'][i])) for i in range(nTAI)]) omnidata['ticks'] = spt.Ticktock(omnidata['UTC'], 'UTC') omnidata['RDT'] = omnidata['ticks'].RDT del omnidata['ticks'] #Can be quickly regenerated on import startyear = omnidata['Year'][-1] del omnidata['Year'] del omnidata['Hr'] # Supplement with daily files print('Supplementing with latest Q-D daily files,' ' this will take a while...') dailyomnidata = _get_qindenton_daily(cached=cached, startyear=startyear) # Find where new files start idx = np.searchsorted(omnidata['UTC'], dailyomnidata['UTC'][0]) for k in sorted(omnidata.keys()): if k == 'Qbits': for qk in sorted(omnidata[k].keys()): omnidata[k][qk] = spacepy.dmarray(np.concatenate(( omnidata[k][qk][:idx, ...], dailyomnidata[k][qk]))) else: omnidata[k] = spacepy.dmarray(np.concatenate(( omnidata[k][:idx, ...], dailyomnidata[k]))) print("Now saving... ") ##for now, make one file -- think about whether monthly/annual files makes sense toHDF5(omni_fname_h5, omnidata) print('Complete.') if omni2 == True: omni2_url = config['omni2_url'] if omni2_url.endswith('.zip'): # adding missing values from original omni2 print("Retrieving OMNI2 file ...") get_url(omni2_url, omni2_fname_zip, progressbar, cached=cached) fh_zip = zipfile.ZipFile(omni2_fname_zip) flist = fh_zip.namelist() if len(flist) != 1: raise RuntimeError('Unable to find OMNI2 file in zip file.') file_to_read = flist[0] if os.path.dirname(file_to_read): raise RuntimeError('Unexpected contents of OMNI2 zip file.') td = tempfile.mkdtemp() try: fh_zip.extract(file_to_read, td) omnicdf = fromCDF(os.path.join(td, file_to_read)) finally: fh_zip.close() os.remove(omni2_fname_zip) shutil.rmtree(td) else: omnicdf = _get_cdaweb_omni2(omni2_url) if omnicdf is not None: #If no new data, skip all this #add RDT omnicdf['RDT'] = spt.Ticktock(omnicdf['Epoch'],'UTC').RDT #remove keys that get in the way del omnicdf['Hour'] del omnicdf['Year'] del omnicdf['Decimal_Day'] # save as HDF5 toHDF5(omni2_fname_h5, omnicdf) if leapsecs == True: print("Retrieving leapseconds file ... ") get_url(config['leapsec_url'], leapsec_fname, progressbar, cached=cached) # Reload leap seconds if they've already been used. if 'spacepy.time' in sys.modules\ and hasattr(sys.modules['spacepy.time'], 'TAIleaps'): sys.modules['spacepy.time']._read_leaps() if PSDdata == True: print("Retrieving PSD sql database") get_url(config['psddata_url'], PSDdata_fname, progressbar, cached=cached) return datadir
[docs]def indsFromXrange(inxrange): '''return the start and end indices implied by an xrange, useful when xrange is zero-length Parameters ========== inxrange : xrange input xrange object to parse Returns ======= list of int List of start, stop indices in the xrange. The return value is not defined if a stride is specified or if stop is before start (but will work when stop equals start). Examples ======== >>> import spacepy.toolbox as tb >>> foo = xrange(23, 39) >>> foo[0] 23 >>> tb.indsFromXrange(foo) [23, 39] >>> foo1 = xrange(23, 23) >>> tb.indsFromXrange(foo) #indexing won't work in this case [23, 23] ''' if not isinstance(inxrange, xrange): return None valstr = inxrange.__str__() if ',' not in valstr: res = re.search(r'(\d+)', valstr) retval = [int(0), int(res.group(1))] else: res = re.search(r'(\d+), (\d+)', valstr) retval = [int(res.group(1)), int(res.group(2))] return retval
[docs]def progressbar(count, blocksize, totalsize, text='Download Progress'): """ print a progress bar with urllib.urlretrieve reporthook functionality Examples ======== >>> import spacepy.toolbox as tb >>> import urllib >>> urllib.urlretrieve(config['psddata_url'], PSDdata_fname, reporthook=tb.progressbar) """ percent = count * blocksize * 100. / totalsize if percent > 100: percent = 100. sys.stdout.write("\r{} ...{:.0f}%".format(text, percent)) if percent >= 100: sys.stdout.write("\n") sys.stdout.flush()
[docs]def windowMean(data, time=[], winsize=0, overlap=0, st_time=None, op=np.mean): """ Windowing mean function, window overlap is user defined Parameters ========== data : array_like 1D series of points time : list (optional) series of timestamps, optional (format as numeric or datetime) For non-overlapping windows set overlap to zero. winsize : integer or datetime.timedelta (optional) window size overlap : integer or datetime.timedelta (optional) amount of window overlap st_time : datetime.datetime (optional) for time-based averaging, a start-time other than the first point can be specified op : callable (optional) the operator to be called, default numpy.mean Returns ======= out : tuple the windowed mean of the data, and an associated reference time vector Examples ======== For non-overlapping windows set overlap to zero. e.g. (time-based averaging) Given a data set of 100 points at hourly resolution (with the time tick in the middle of the sample), the daily average of this, with half-overlapping windows is calculated: >>> import spacepy.toolbox as tb >>> from datetime import datetime, timedelta >>> wsize = datetime.timedelta(days=1) >>> olap = datetime.timedelta(hours=12) >>> data = [10, 20]*50 >>> time = [datetime.datetime(2001,1,1) + datetime.timedelta(hours=n, minutes = 30) for n in range(100)] >>> outdata, outtime = tb.windowMean(data, time, winsize=wsize, overlap=olap, st_time=datetime.datetime(2001,1,1)) >>> outdata, outtime ([15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0], [datetime.datetime(2001, 1, 1, 12, 0), datetime.datetime(2001, 1, 2, 0, 0), datetime.datetime(2001, 1, 2, 12, 0), datetime.datetime(2001, 1, 3, 0, 0), datetime.datetime(2001, 1, 3, 12, 0), datetime.datetime(2001, 1, 4, 0, 0), datetime.datetime(2001, 1, 4, 12, 0)]) When using time-based averaging, ensure that the time tick corresponds to the middle of the time-bin to which the data apply. That is, if the data are hourly, say for 00:00-01:00, then the time applied should be 00:30. If this is not done, unexpected behaviour can result. e.g. (pointwise averaging), >>> outdata, outtime = tb.windowMean(data, winsize=24, overlap=12) >>> outdata, outtime ([15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0], [12.0, 24.0, 36.0, 48.0, 60.0, 72.0, 84.0]) where winsize and overlap are numeric, in this example the window size is 24 points (as the data are hourly) and the overlap is 12 points (a half day). The output vectors start at winsize/2 and end at N-(winsize/2), the output time vector is basically a reference to the nth point in the original series. **note** This is a quick and dirty function - it is NOT optimized, at all. """ #check inputs and initialize #Set resolution to 1 if no times supplied if len(time) == 0: startpt, res = 0, 1 time = list(range(len(data))) pts = True else: if len(data) != len(time): raise ValueError('windowmean error: data and time must have same length') #First check if datetime objects try: assert type(winsize) == datetime.timedelta assert type(overlap) == datetime.timedelta except AssertionError: raise TypeError('windowmean error: winsize/overlap must be timedeltas') pts = False #force time-based averaging if (type(time[0]) != datetime.datetime): startpt = time[0] #now actually do windowing mean outdata, outtime = [], [] data = np.array(data) if pts: #loop for fixed number of points in window try: inttypes = (int, long) except NameError: inttypes = (int,) if not isinstance(winsize, inttypes): winsize = int(round(winsize)) warnings.warn('windowmean: non-integer windowsize, rounding to %d' \ % winsize) if winsize < 1: winsize = 1 warnings.warn('windowmean: window length < 1, defaulting to 1') if overlap >= winsize: overlap = winsize - 1 warnings.warn('''windowmean: overlap longer than window, truncated to %d''' % overlap) lastpt = winsize-1 #set last point to end of window size while lastpt < len(data): datwin = np.ma.masked_where(np.isnan(data[startpt:startpt+winsize]), \ data[startpt:startpt+winsize]) getmean = op(datwin.compressed()) #mean of window, excl. NaNs gettime = (time[startpt+winsize] - time[startpt])/2. \ + time[startpt]#new timestamp startpt = startpt+winsize-overlap lastpt = startpt+winsize outdata.append(getmean) #construct output arrays outtime.append(gettime) else: #loop with time-based window lastpt = time[0] + winsize delta = datetime.timedelta(microseconds=1) #TODO: replace this with an explicit check for times on the boundary? if st_time: startpt = st_time else: startpt = time[0] if overlap >= winsize: raise ValueError('Overlap requested greater than size of window') while startpt < time[-1]: getinds = tOverlapHalf([startpt,startpt+winsize-delta], time, presort=True) if getinds: #if not None getdata = np.ma.masked_where(np.isnan(data[getinds[0]:getinds[-1]+1]), data[getinds[0]:getinds[-1]+1]) getmean = op(getdata.compressed()) #find mean excluding NaNs else: getmean = np.nan gettime = startpt + winsize//2 #new timestamp -floordiv req'd with future division startpt = startpt + winsize - overlap #advance window start lastpt = startpt + winsize outdata.append(getmean) #construct output arrays outtime.append(gettime) return outdata, outtime
[docs]def medAbsDev(series, scale=False): """ Calculate median absolute deviation of a given input series Median absolute deviation (MAD) is a robust and resistant measure of the spread of a sample (same purpose as standard deviation). The MAD is preferred to the inter-quartile range as the inter-quartile range only shows 50% of the data whereas the MAD uses all data but remains robust and resistant. See e.g. Wilks, Statistical methods for the Atmospheric Sciences, 1995, Ch. 3. For additional details on the scaling, see Rousseeuw and Croux, J. Amer. Stat. Assoc., 88 (424), pp. 1273-1283, 1993. Parameters ========== series : array_like the input data series Other Parameters ================ scale : bool if True (default: False), scale to standard deviation of a normal distribution Returns ======= out : float the median absolute deviation Examples ======== Find the median absolute deviation of a data set. Here we use the log- normal distribution fitted to the population of sawtooth intervals, see Morley and Henderson, Comment, Geophysical Research Letters, 2009. >>> import numpy >>> import spacepy.toolbox as tb >>> numpy.random.seed(8675301) >>> data = numpy.random.lognormal(mean=5.1458, sigma=0.302313, size=30) >>> print data array([ 181.28078923, 131.18152745, ... , 141.15455416, 160.88972791]) >>> tb.medAbsDev(data) 28.346646721370192 **note** This implementation is robust to presence of NaNs """ #ensure input is numpy array (and make 1-D) series = (np.array(series, dtype=float)).ravel() #mask for NaNs series = np.ma.masked_where(np.isnan(series),series) #get median absolute deviation of unmasked elements perc50 = np.median(series.compressed()) mad = np.median(abs(series.compressed()-perc50)) if scale: mad *= 1.4826 #scale so that MAD is same as SD for normal distr. return mad
[docs]def binHisto(data, verbose=False): """ Calculates bin width and number of bins for histogram using Freedman-Diaconis rule, if rule fails, defaults to square-root method The Freedman-Diaconis method is detailed in: Freedman, D., and P. Diaconis (1981), On the histogram as a density estimator: L2 theory, Z. Wahrscheinlichkeitstheor. Verw. Geb., 57, 453–476 and is also described by: Wilks, D. S. (2006), Statistical Methods in the Atmospheric Sciences, 2nd ed. Parameters ========== data : array_like list/array of data values verbose : boolean (optional) print out some more information Returns ======= out : tuple calculated width of bins using F-D rule, number of bins (nearest integer) to use for histogram Examples ======== >>> import numpy, spacepy >>> import matplotlib.pyplot as plt >>> numpy.random.seed(8675301) >>> data = numpy.random.randn(1000) >>> binw, nbins = spacepy.toolbox.binHisto(data) >>> print(nbins) 19 >>> p = plt.hist(data, bins=nbins, histtype='step', density=True) See Also ======== matplotlib.pyplot.hist """ pul = np.percentile(data, (25, 75)) #get confidence interval ql, qu = pul[0], pul[1] iqr = qu-ql binw = 2.*iqr/(len(data)**(1./3.)) if binw != 0: nbins = int(round((max(data)-min(data))/binw)) # if nbins is 0, NaN or inf don't use the F-D rule just use sqrt(num) rule if binw == 0 or nbins == 0 or not np.isfinite(nbins) or nbins >= len(data)/2.: nbins = int(round(np.sqrt(len(data)))) binw = len(data)/nbins if verbose: print("Used sqrt rule") else: if verbose: print("Used F-D rule") return (binw, nbins)
[docs]def bootHisto(data, inter=90., n=1000, seed=None, plot=False, target=None, figsize=None, loc=None, **kwargs): """Bootstrap confidence intervals for a histogram. All other keyword arguments are passed to :func:`numpy.histogram` or :func:`matplotlib.pyplot.bar`. .. versionchanged:: 0.2.3 This argument pass-through did not work in earlier versions of SpacePy. Parameters ========== data : array_like list/array of data values inter : float (optional; default 90) percentage confidence interval to return. Default 90% (i.e. lower CI will be 5% and upper will be 95%) n : int (optional; default 1000) number of bootstrap iterations seed : int (optional) Optional seed for the random number generator. If not specified; numpy generator will not be reseeded. plot : bool (optional) Plot the result. Plots if True or ``target``, ``figsize``, or ``loc`` specified. target : (optional) Target on which to plot the figure (figure or axes). See :func:`spacepy.plot.utils.set_target` for details. figsize : tuple (optional) Passed to :func:`spacepy.plot.utils.set_target`. loc : int (optional) Passed to :func:`spacepy.plot.utils.set_target`. Returns ======= out : tuple tuple of bin_edges, low, high, sample[, bars]. Where ``bin_edges`` is the edges of the bins used; ``low`` is the histogram with the value for each bin from the bottom of that bin's confidence interval; ``high`` similarly for the top; ``sample`` is the histogram of the input sample without resampling. If plotting, also returned is ``bars``, the container object returned from matplotlib. Notes ===== .. versionadded:: 0.2.1 The confidence intervals are calculated for each bin individually and thus the resulting low/high histograms may not have actually occurred in the calculation from the surrogates. If using a probability density histogram, this can have "interesting" implications for interpretation. Examples ======== .. plot:: :include-source: >>> import numpy.random >>> import spacepy.toolbox >>> numpy.random.seed(0) >>> data = numpy.random.randn(1000) >>> bin_edges, low, high, sample, bars = spacepy.toolbox.bootHisto( ... data, plot=True) See Also ======== binHisto plot.utils.set_target numpy.histogram matplotlib.pyplot.hist """ import spacepy.poppy histogram_allowed_kwargs = ( 'bins', 'range', 'normed', 'weights', 'density') histogram_kwargs = {k: v for k, v in kwargs.items() if k in histogram_allowed_kwargs} bar_kwargs = {k: v for k, v in kwargs.items() if k not in histogram_allowed_kwargs} sample, bin_edges = np.histogram(data, **histogram_kwargs) histogram_kwargs['bins'] = bin_edges ci_low, ci_high = spacepy.poppy.boots_ci( data, n, inter, lambda x: np.histogram(x, **histogram_kwargs)[0], seed=seed, nretvals=len(bin_edges) - 1) if not plot and all([x is None for x in (target, figsize, loc)]): return bin_edges, ci_low, ci_high, sample import spacepy.plot.utils _, ax = spacepy.plot.utils.set_target( target, figsize=figsize, loc=(111 if loc is None else loc)) if 'ecolor' not in bar_kwargs: bar_kwargs['ecolor'] = 'k' bars = ax.bar( bin_edges[:-1], height=sample, width=np.diff(bin_edges), align='edge', yerr=np.stack((sample - ci_low, ci_high - sample)), **bar_kwargs) return bin_edges, ci_low, ci_high, sample, bars
[docs]def logspace(min, max, num, **kwargs): """ Returns log-spaced bins. Same as numpy.logspace except the min and max are the min and max not log10(min) and log10(max) Parameters ========== min : float minimum value max : float maximum value num : integer number of log spaced bins Other Parameters ================ kwargs : dict additional keywords passed into matplotlib.dates.num2date Returns ======= out : array log-spaced bins from min to max in a numpy array Notes ===== This function works on both numbers and datetime objects Examples ======== >>> import spacepy.toolbox as tb >>> tb.logspace(1, 100, 5) array([ 1. , 3.16227766, 10. , 31.6227766 , 100. ]) See Also ======== geomspace linspace """ if isinstance(min, datetime.datetime): from matplotlib.dates import date2num, num2date ans = num2date(np.logspace(np.log10(date2num(min)), np.log10(date2num(max)), num, **kwargs)) ans = spt.no_tzinfo(ans) return np.array(ans) else: return np.logspace(np.log10(min), np.log10(max), num, **kwargs)
[docs]def linspace(min, max, num, **kwargs): """ Returns linear-spaced bins. Same as numpy.linspace except works with datetime and is faster Parameters ========== min : float, datetime minimum value max : float, datetime maximum value num : integer number of linear spaced bins Other Parameters ================ kwargs : dict additional keywords passed into matplotlib.dates.num2date Returns ======= out : array linear-spaced bins from min to max in a numpy array Notes ===== This function works on both numbers and datetime objects Examples ======== >>> import spacepy.toolbox as tb >>> tb.linspace(1, 10, 4) array([ 1., 4., 7., 10.]) See Also ======== geomspace logspace """ if hasattr(min, 'shape') and min.shape == (): min = min.item() if hasattr(max, 'shape') and max.shape == (): max = max.item() if isinstance(min, datetime.datetime): from matplotlib.dates import date2num, num2date ans = num2date(np.linspace(date2num(min), date2num(max), num, **kwargs)) ans = spt.no_tzinfo(ans) return np.array(ans) else: return np.linspace(min, max, num, **kwargs)
[docs]def geomspace(start, ratio=None, stop=False, num=50): """ Returns geometrically spaced numbers. Parameters ========== start : float The starting value of the sequence. ratio : float (optional) The ratio between subsequent points stop : float (optional) End value, if this is selected `num` is overridden num : int (optional) Number of samples to generate. Default is 50. Returns ======= seq : array geometrically spaced sequence See Also ======== linspace logspace Examples ======== To get a geometric progression between 0.01 and 3 in 10 steps >>> import spacepy.toolbox as tb >>> tb.geomspace(0.01, stop=3, num=10) [0.01, 0.018846716378431192, 0.035519871824902655, 0.066943295008216955, 0.12616612944575134, 0.23778172582285118, 0.44814047465571644, 0.84459764235318191, 1.5917892219322083, 2.9999999999999996] To get a geometric progression with a specified ratio, say 10 >>> import spacepy.toolbox as tb >>> tb.geomspace(0.01, ratio=10, num=5) [0.01, 0.10000000000000001, 1.0, 10.0, 100.0] """ if not ratio and stop != False: ratio = (stop/start)**(1/(num-1)) seq = [] seq.append(start) if stop == False: for j in range(1, num): seq.append(seq[j-1]*ratio) return seq else: val, j = start, 1 while val <= stop or np.allclose(val, stop, ): val = seq[j-1]*ratio seq.append(val) j+=1 return seq[:-1]
[docs]def arraybin(array, bins): """ Split a sequence into subsequences based on value. Given a sequence of values and a sequence of values representing the division between bins, return the indices grouped by bin. Parameters ========== array : array_like the input sequence to slice, must be sorted in ascending order bins : array_like dividing lines between bins. Number of bins is len(bins)+1, value that exactly equal a dividing value are assigned to the higher bin Returns ======= out : list indices for each bin (list of lists) Examples ======== >>> import spacepy.toolbox as tb >>> tb.arraybin(range(10), [4.2]) [[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]] """ bin_it = lambda value: (i for i in range(len(array)) if array[i] >= value) splits = [next(bin_it(value), len(array)) for value in bins] return [list(range(start_idx, stop_idx)) for (start_idx, stop_idx) in zip([0] + splits, splits + [len(array)])]
[docs]def mlt2rad(mlt, midnight = False): """ Convert mlt values to radians for polar plotting transform mlt angles to radians from -pi to pi referenced from noon by default Parameters ========== mlt : numpy array array of mlt values midnight : boolean (optional) reference to midnight instead of noon Returns ======= out : numpy array array of radians Examples ======== >>> from numpy import array >>> mlt2rad(array([3,6,9,14,22])) array([-2.35619449, -1.57079633, -0.78539816, 0.52359878, 2.61799388]) See Also ======== rad2mlt """ if midnight: try: mlt_arr = [val + 12 for val in mlt] except TypeError: mlt_arr = mlt + 12 else: mlt_arr = mlt try: rad_arr = [(val-12)*np.pi/12. for val in mlt_arr] except TypeError: rad_arr = (mlt_arr-12)*np.pi/12. return rad_arr
[docs]def rad2mlt(rad, midnight=False): """ Convert radians values to mlt transform radians from -pi to pi to mlt referenced from noon by default Parameters ========== rad : numpy array array of radian values midnight : boolean (optional) reference to midnight instead of noon Returns ======= out : numpy array array of mlt values Examples ======== >>> rad2mlt(array([0,pi, pi/2.])) array([ 12., 24., 18.]) See Also ======== mlt2rad """ if midnight: rad_arr = rad + np.pi else: rad_arr = rad mlt_arr=rad_arr*(12/np.pi) + 12 return mlt_arr
[docs]def pmm(*args): """ print min and max of input arrays Parameters ========== a : array-like arbitrary number of input arrays (or lists) Returns ======= out : list list of min, max for each array Examples ======== >>> import spacepy.toolbox as tb >>> from numpy import arange >>> tb.pmm(arange(10), arange(10)+3) [[0, 9], [3, 12]] """ ans = [] for a in args: try: ind = np.isfinite(a) except TypeError: ind = np.arange(len(a)).astype(int) import warnings warnings.warn('pmm: Unable to exclude non-finite values, results may be incorrect', RuntimeWarning) try: ans.append([np.min(a[ind]), np.max(a[ind])]) except TypeError: a_tmp = np.asarray(a) if a_tmp.dtype.type in [np.dtype('S').type, np.dtype('U').type]: a_tmp = np.require(a_tmp, dtype=object) ans.append([np.min(a_tmp[ind]), np.max(a_tmp[ind])]) return ans
[docs]def getNamedPath(name): """ Return the full path of a parent directory with name as the leaf Parameters ========== name : string the name of the parent directory to locate Examples ======== Run from a directory /mnt/projects/dream/bin/Ephem with 'dream' as the name, this function would return '/mnt/projects/dream' """ def findNamed(path): pp = os.path.split(path) if pp[-1] == '': return None if pp[-1] != name: path = findNamed(pp[0]) return path return findNamed(os.getcwd())
[docs]def query_yes_no(question, default="yes"): """ Ask a yes/no question via raw_input() and return their answer. "question" is a string that is presented to the user. "default" is the presumed answer if the user just hits <Enter>. It must be "yes" (the default), "no" or None (meaning an answer is required of the user). The "answer" return value is one of "yes" or "no". Parameters ========== question : string the question to ask default : string (optional) Returns ======= out : string answer ('yes' or 'no') Examples ======== >>> import spacepy.toolbox as tb >>> tb.query_yes_no('Ready to go?') Ready to go? [Y/n] y 'yes' """ valid = {"yes":"yes", "y":"yes", "ye":"yes", "no":"no", "n":"no"} if default == None: prompt = " [y/n] " elif default == "yes": prompt = " [Y/n] " elif default == "no": prompt = " [y/N] " else: raise ValueError("invalid default answer: '%s'" % default) while 1: sys.stdout.write(question + prompt) if sys.version_info[0]==2: choice = raw_input().lower() elif sys.version_info[0]>2: choice = input().lower() if default is not None and choice == '': return default elif choice in valid: return valid[choice] else: sys.stdout.write("Please respond with 'yes' or 'no' "\ "(or 'y' or 'n').\n")
[docs]def interpol(newx, x, y, wrap=None, **kwargs): """ 1-D linear interpolation with interpolation of hours/longitude Parameters ========== newx : array_like x values where we want the interpolated values x : array_like x values of the original data (must be monotonically increasing or wrapping) y : array_like y values of the original data wrap : string, optional for continuous x data that wraps in y at 'hours' (24), 'longitude' (360), or arbitrary value (int, float) kwargs : dict additional keywords, currently accepts baddata that sets baddata for masked arrays Returns ======= out : numpy.masked_array interpolated data values for new abscissa values Examples ======== For a simple interpolation >>> import spacepy.toolbox as tb >>> import numpy >>> x = numpy.arange(10) >>> y = numpy.arange(10) >>> tb.interpol(numpy.arange(5)+0.5, x, y) array([ 0.5, 1.5, 2.5, 3.5, 4.5]) To use the wrap functionality, without the wrap keyword you get the wrong answer >>> y = range(24)*2 >>> x = range(len(y)) >>> tb.interpol([1.5, 10.5, 23.5], x, y, wrap='hour').compressed() # compress removed the masked array array([ 1.5, 10.5, 23.5]) >>> tb.interpol([1.5, 10.5, 23.5], x, y) array([ 1.5, 10.5, 11.5]) """ if 'baddata' in kwargs: y = np.ma.masked_equal(y, kwargs['baddata']) x = np.ma.masked_array(x) x.mask = y.mask kwargs.__delitem__('baddata') else: tst = np.ma.core.MaskedArray if type(x)!=tst or type(y)!=tst or type(newx)!=tst: x = np.ma.masked_array(x) y = np.ma.masked_array(y) newx = np.ma.masked_array(newx) def wrap_interp(xout, xin, yin, sect): dpsect=360/sect yc = np.cos(np.deg2rad(y*dpsect)) ys = np.sin(np.deg2rad(y*dpsect)) new_yc = np.interp(newx, x.compressed(), yc.compressed(), **kwargs) new_ys = np.interp(newx, x.compressed(), ys.compressed(), **kwargs) try: new_bad = np.interp(newx, x, y.mask) except ValueError: new_bad = np.zeros((len(newx))) newy = np.rad2deg(np.arctan(new_ys/new_yc))/dpsect #1st quadrant is O.K #2nd quadrant idx = [n for n in range(len(new_yc)) if new_yc[n]<0 and new_ys[n]>0] newy[idx] = sect/2 + newy[idx] #print('Sector 2 inds: %s' % idx) #3rd quadrant idx = [n for n in range(len(new_yc)) if new_yc[n]<0 and new_ys[n]<0] newy[idx] = sect/2 + newy[idx] #print('Sector 3 inds: %s' % idx) #4th quadrant idx = [n for n in range(len(new_yc)) if new_yc[n]>0 and new_ys[n]<0] newy[idx] = sect + newy[idx] #print('Sector 4 inds: %s' % idx) new_bad = np.ma.make_mask(new_bad) newy = np.ma.masked_array(newy, mask=new_bad) return newy if wrap=='hour': newy = wrap_interp(newx, x.compressed(), y.compressed(), 24) elif wrap=='lon': newy = wrap_interp(newx, x.compressed(), y.compressed(), 360) elif isinstance(wrap, numbers.Real): newy = wrap_interp(newx, x.compressed(), y.compressed(), wrap) else: newy = np.interp(newx, x.compressed(), y.compressed(), **kwargs) return newy
# -----------------------------------------------
[docs]@spacepy.deprecated( '0.2.2', 'moved to spacepy.coordinates', 'use :func:`spacepy.coordinates.quaternionNormalize`') def quaternionNormalize(Qin, scalarPos='last'): """Given an input quaternion, return the unit quaternion.""" import spacepy.coordinates return spacepy.coordinates.quaternionNormalize(Qin, scalarPos=scalarPos)
[docs]@spacepy.deprecated( '0.2.2', 'moved to spacepy.coordinates', 'use :func:`spacepy.coordinates.quaternionRotateVector`') def quaternionRotateVector(Qin, Vin, scalarPos='last', normalize=True): """Given quaternion and vector, return vector rotated by quaternion.""" import spacepy.coordinates return spacepy.coordinates.quaternionRotateVector( Qin, Vin, scalarPos=scalarPos, normalize=normalize)
[docs]@spacepy.deprecated( '0.2.2', 'moved to spacepy.coordinates', 'use :func:`spacepy.coordinates.quaternionMultiply`') def quaternionMultiply(Qin1, Qin2, scalarPos='last'): """Given quaternions, return the product.""" import spacepy.coordinates return spacepy.coordinates.quaternionMultiply( Qin1, Qin2, scalarPos=scalarPos)
[docs]@spacepy.deprecated( '0.2.2', 'moved to spacepy.coordinates', 'use :func:`spacepy.coordinates.quaternionConjugate`') def quaternionConjugate(Qin, scalarPos='last'): """Given an input quaternion, return the conjugated.""" import spacepy.coordinates return spacepy.coordinates.quaternionConjugate(Qin, scalarPos=scalarPos)
# -----------------------------------------------
[docs]def normalize(vec, low=0.0, high=1.0): """ Given an input vector normalize the vector to a given range Parameters ========== vec : array_like input vector to normalize low : float minimum value to scale to, default 0.0 high : float maximum value to scale to, default 1.0 Returns ======= out : array_like normalized vector Examples ======== >>> import spacepy.toolbox as tb >>> tb.normalize([1,2,3]) [0.0, 0.5, 1.0] """ return np.interp(vec, (np.nanmin(vec), np.nanmax(vec)), (low, high))
[docs]def intsolve(func, value, start=None, stop=None, maxit=1000): """ Find the function input such that definite integral is desired value. Given a function, integrate from an (optional) start point until the integral reached a desired value, and return the end point of the integration. Parameters ========== func : callable function to integrate, must take single parameter value : float desired final value of the integral start : float (optional) value at which to start integration, default -Infinity stop : float (optional) value at which to stop integration, default +Infinity maxit : integer maximum number of iterations Returns ======= out : float x such that the integral of L{func} from L{start} to x is L{value} **Note:** Assumes func is everywhere positive, otherwise solution may be multi-valued. """ from scipy import inf from scipy.integrate import quad from warnings import warn if start is None: start = -inf if stop is None: stop = inf lower_bound = start upper_bound = stop it = 0 while it < maxit: it += 1 if upper_bound == inf: if lower_bound == -inf: test_bound = 0 elif lower_bound < 1.0: test_bound = 1.0 else: test_bound = lower_bound * 2 else: if lower_bound == -inf: if upper_bound > -1.0: test_bound = -1.0 else: test_bound = upper_bound * 2 else: test_bound = (lower_bound + upper_bound) / 2.0 (test_value, err) = quad(func, start, test_bound) if abs(value - test_value) <= err: break elif value < test_value: upper_bound = test_bound else: lower_bound = test_bound if abs(value - test_value) > err: warn('Difference between desired value and actual is ' + str(abs(value - test_value)) + ', greater than integral error ' + str(err), UserWarning, stacklevel=2) return test_bound
[docs]def dist_to_list(func, length, min=None, max=None): """ Convert a probability distribution function to a list of values This is a deterministic way to produce a known-length list of values matching a certain probability distribution. It is likely to be a closer match to the distribution function than a random sampling from the distribution. Parameters ========== func : callable function to call for each possible value, returning probability density at that value (does not need to be normalized.) length : int number of elements to return min : float minimum value to possibly include max : float maximum value to possibly include Examples ======== >>> import matplotlib >>> import numpy >>> import spacepy.toolbox as tb >>> gauss = lambda x: math.exp(-(x ** 2) / (2 * 5 ** 2)) / (5 * math.sqrt(2 * math.pi)) >>> vals = tb.dist_to_list(gauss, 1000, -numpy.inf, numpy.inf) >>> print vals[0] -16.45263... >>> p1 = matplotlib.pyplot.hist(vals, bins=[i - 10 for i in range(21)], facecolor='green') >>> matplotlib.pyplot.hold(True) >>> x = [i / 100.0 - 10.0 for i in range(2001)] >>> p2 = matplotlib.pyplot.plot(x, [gauss(i) * 1000 for i in x], 'red') >>> matplotlib.pyplot.draw() """ from scipy import inf from scipy.integrate import quad if min is None: min = -inf if max is None: max = inf total = quad(func, min, max)[0] step = float(total) / length return [intsolve(func, (0.5 + i) * step, min, max) for i in range(length)]
[docs]def bin_center_to_edges(centers): """ Convert a list of bin centers to their edges Given a list of center values for a set of bins, finds the start and end value for each bin. (start of bin n+1 is assumed to be end of bin n). Useful for e.g. matplotlib.pyplot.pcolor. Edge between bins n and n+1 is arithmetic mean of the center of n and n+1; edge below bin 0 and above last bin are established to make these bins symmetric about their center value. Parameters ========== centers : list list of center values for bins Returns ======= out : list list of edges for bins **note:** returned list will be one element longer than centers Examples ======== >>> import spacepy.toolbox as tb >>> tb.bin_center_to_edges([1,2,3]) [0.5, 1.5, 2.5, 3.5] """ edges = bin_edges_to_center(centers) edges = np.append(centers[0]-(edges[0]-centers[0]), edges) edges = np.append(edges, centers[-1]+(centers[-1]-edges[-1])) return edges
[docs]def bin_edges_to_center(edges): """ Convert a list of bin edges to their centers Given a list of edge values for a set of bins, finds the center of each bin. (start of bin n+1 is assumed to be end of bin n). Center of bin n is arithmetic mean of the edges of the adjacent bins. Parameters ========== edges : list list of edge values for bins Returns ======= out : numpy.ndarray array of centers for bins **note:** returned array will be one element shorter than edges Examples ======== >>> import spacepy.toolbox as tb >>> tb.bin_center_to_edges([1,2,3]) [0.5, 1.5, 2.5, 3.5] """ df = np.diff(edges) if isinstance(df[0], datetime.timedelta) and sys.version_info[0:2]<=(3,2): return edges[:-1] + df//2 else: return edges[:-1] + df/2
[docs]def thread_job(job_size, thread_count, target, *args, **kwargs): """ Split a job into subjobs and run a thread for each Each thread spawned will call L{target} to handle a slice of the job. This is only useful if a job: 1. Can be split into completely independent subjobs 2. Relies heavily on code that does not use the Python GIL, e.g. numpy or ctypes code 3. Does not return a value. Either pass in a list/array to hold the result, or see L{thread_map} Examples ======== squaring 100 million numbers: >>> import numpy >>> import spacepy.toolbox as tb >>> numpy.random.seed(8675301) >>> a = numpy.random.randint(0, 100, [100000000]) >>> b = numpy.empty([100000000], dtype='int64') >>> def targ(in_array, out_array, start, count): \ out_array[start:start + count] = in_array[start:start + count] ** 2 >>> tb.thread_job(len(a), 0, targ, a, b) >>> print(b[0:5]) [2704 7225 196 1521 36] This example: - Defines a target function, which will be called for each thread. It is usually necessary to define a simple "wrapper" function like this to provide the correct call signature. - The target function receives inputs C{in_array} and C{out_array}, which are not touched directly by C{thread_job} but are passed through in the call. In this case, C{a} gets passed as C{in_array} and C{b} as C{out_array} - The target function also receives the start and number of elements it needs to process. For each thread where the target is called, these numbers are different. Parameters ========== job_size : int Total size of the job. Often this is an array size. thread_count : int Number of threads to spawn. If =0 or None, will spawn as many threads as there are cores available on the system. (Each hyperthreading core counts as 2.) Generally this is the Right Thing to do. If NEGATIVE, will spawn abs(thread_count) threads, but will run them sequentially rather than in parallel; useful for debugging. target : callable Python callable (generally a function, may also be an imported ctypes function) to run in each thread. The *last* two positional arguments passed in will be a "start" and a "subjob size," respectively; frequently this will be the start index and the number of elements to process in an array. args : sequence Arguments to pass to L{target}. If L{target} is an instance method, self must be explicitly passed in. start and subjob_size will be appended. kwargs : dict keyword arguments to pass to L{target}. """ try: import threading except: return (target(*(args + (0, job_size)), **kwargs), ) if thread_count == None or thread_count == 0: try: import multiprocessing thread_count = multiprocessing.cpu_count() except: #Do something not too stupid thread_count = 8 if thread_count < 0: thread_count *= -1 seq = True else: seq = False count = float(job_size) / thread_count starts = [int(count * i + 0.5) for i in range(thread_count)] subsize = [(starts[i + 1] if i < thread_count - 1 else job_size) - starts[i] for i in range(thread_count)] threads = [] for i in range(thread_count): t = threading.Thread( None, target, None, args + (starts[i], subsize[i]), kwargs) t.start() if seq: t.join() else: threads.append(t) if not seq: for t in threads: t.join()
[docs]def thread_map(target, iterable, thread_count=None, *args, **kwargs): """ Apply a function to every element of a list, in separate threads Interface is similar to multiprocessing.map, except it runs in threads This is made largely obsolete in python3 by from concurrent import futures Examples ======== find totals of several arrays >>> import numpy >>> from spacepy import toolbox >>> inputs = range(100) >>> totals = toolbox.thread_map(numpy.sum, inputs) >>> print(totals[0], totals[50], totals[99]) (0, 50, 99) >>> # in python3 >>> from concurrent import futures >>> with futures.ThreadPoolExecutor(max_workers=4) as executor: ...: for ans in executor.map(numpy.sum, [0,50,99]): ...: print ans #0 #50 #99 Parameters ========== target : callable Python callable to run on each element of iterable. For each call, an element of iterable is appended to args and both args and kwargs are passed through. Note that this means the iterable element is always the *last* positional argument; this allows the specification of self as the first argument for method calls. iterable : iterable elements to pass to each call of L{target} args : sequence arguments to pass to target before each element of iterable thread_count : integer Number of threads to spawn; see L{thread_job}. kwargs : dict keyword arguments to pass to L{target}. Returns ======= out : list return values of L{target} for each item from L{iterable} """ try: jobsize = len(iterable) except TypeError: iterable = list(iterable) jobsize = len(iterable) def array_targ(function, it, retvals, arglist, kwarglist, start, size): for i in range(start, start + size): retvals[i] = function(*(arglist + (it[i],)), **kwarglist) retvals = [None] * jobsize thread_job(jobsize, thread_count, array_targ, target, iterable, retvals, args, kwargs) return retvals
[docs]def eventTimer(Event, Time1): """ Times an event then prints out the time and the name of the event, nice for debugging and seeing that the code is progressing Parameters ========== Event : str Name of the event, string is printed out by function Time1 : time.time the time to difference in the function Returns ======= Time2 : time.time the new time for the next call to EventTimer Examples ======== >>> import spacepy.toolbox as tb >>> import time >>> t1 = time.time() >>> t1 = tb.eventTimer('Test event finished', t1) ('4.40', 'Test event finished') """ Time2 = time.time() print(("%4.2f" % (Time2 - Time1), Event)) return Time2
[docs]def isview(array1, array2=None): """ Returns if an object is a view of another object. More precisely if one array argument is specified True is returned is the arrays owns its data. If two arrays arguments are specified a tuple is returned of if the first array owns its data and the the second if they point at the same memory location Parameters ========== array1 : numpy.ndarray array to query if it owns its data Other Parameters ================ array2 : object (optional) array to query if array1 is a view of this object at the specified memory location Returns ======= out : bool or tuple If one array is specified bool is returned, True is the array owns its data. If two arrays are specified a tuple where the second element is a bool of if the array point at the same memory location Examples ======== import numpy import spacepy.toolbox as tb a = numpy.arange(100) b = a[0:10] tb.isview(a) # False tb.isview(b) # True tb.isview(b, a) # (True, True) tb.isview(b, b) # (True, True) # the conditions are met and numpy cannot tell this """ # deal with the one input case first if array2 is None: try: if array1.base is None: return False return True except AttributeError: return False # if it is not an array then it is not a view # there are two arrays input try: if array1.base is None: return (False, False) return (True, array1.base is array2) except AttributeError: return (False, False) # if it is not an array then it is not a view
[docs]def interweave(a, b): """ given two array-like variables interweave them together. Discussed here: http://stackoverflow.com/questions/5347065/interweaving-two-numpy-arrays Parameters ========== a : array-like first array b : array-like second array Returns ======= out : numpy.ndarray interweaved array """ a = np.asanyarray(a) b = np.asanyarray(b) ans = np.empty((a.size + b.size), dtype=a.dtype) ans[0::2] = a ans[1::2] = b return ans
[docs]class TimeoutError(Exception): """Raised when a time-limited process times out""" pass
[docs]def do_with_timeout(timeout, target, *args, **kwargs): """ Execute a function (or method) with a timeout. Call the function (or method) ``target``, with arguments ``args`` and keyword arguments ``kwargs``. Normally return the return value from ``target``, but if ``target`` takes more than ``timeout`` seconds to execute, raises ``TimeoutError``. .. note:: This is, at best, a blunt instrument. Exceptions from ``target`` may not propagate properly (tracebacks will be hard to follow.) The function which failed to time out may continue to execute until the interpreter exits; trapping the TimeoutError and continuing normally is not recommended. Examples ======== >>> import spacepy.toolbox as tb >>> import time >>> def time_me_out(): ... time.sleep(5) >>> tb.do_with_timeout(0.5, time_me_out) #raises TimeoutError Parameters ========== timeout : float Timeout, in seconds. target : callable Python callable (generally a function, may also be an imported ctypes function) to run. args : sequence Arguments to pass to ``target``. kwargs : dict keyword arguments to pass to ``target``. Raises ====== TimeoutError : If ``target`` does not return in ``timeout`` seconds. Returns ======= out : return value of ``target`` """ import threading class ReturningThread(threading.Thread): def __init__(self, group=None, target=None, name=None, args=(), kwargs={}): #we're handling target, args, kwargs super(ReturningThread, self).__init__(group, name=name) self._target = target self._args = args self._kwargs = kwargs self._retval = None self._exception = None def run(self): try: self._retval = self._target(*self._args, **self._kwargs) except: self._exception = sys.exc_info() def join(self, *args, **kwargs): super(ReturningThread, self).join(*args, **kwargs) if not self._exception is None: try: raise self._exception[1].with_traceback(self._exception[2]) except AttributeError: raise (self._exception[1], None, self._exception[2]) return self._retval t = ReturningThread(None, target, None, args, kwargs) t.start() retval = t.join(timeout) if t.is_alive(): raise TimeoutError() return retval
[docs]def timeout_check_call(timeout, *args, **kwargs): """ Call a subprocess with a timeout. Like :func:`subprocess.check_call`, but will terminate the process and raise :exc:`TimeoutError` if it runs for too long. This will only terminate the single process started; any child processes will remain running (this has implications for, say, spawing shells.) Examples ======== >>> import spacepy.toolbox as tb >>> tb.timeout_check_call(1, 'sleep 30', shell=True) #raises TimeoutError Parameters ========== timeout : float Timeout, in seconds. Fractions are acceptable but the resolution is of order 100ms. args : sequence Arguments passed through to :class:`subprocess.Popen` kwargs : dict keyword arguments to pass to :class:`subprocess.Popen` Raises ====== TimeoutError : If subprocess does not return in ``timeout`` seconds. CalledProcessError : if command has non-zero exit status Returns ======= out : int 0 on successful completion """ resolution = 0.1 pro = subprocess.Popen(*args, **kwargs) starttime = time.time() while pro.poll() is None: time.sleep(resolution) if time.time() - starttime > timeout: to = time.time() - starttime pro.terminate() time.sleep(resolution) if pro.poll() is None: pro.kill() raise TimeoutError('Timed out after {0:.1f} seconds'.format(to)) if pro.returncode: raise subprocess.CalledProcessError( pro.returncode, kwargs['args'] if 'args' in kwargs else args[0]) return 0
[docs]def poisson_fit(data, initial=None, method='Powell'): """ Fit a Poisson distribution to data using the method and initial guess provided. Parameters ========== data : array-like Data to fit a Poisson distribution to. initial : int or None initial guess for the fit, if None np.median(data) is used method : str method passed to scipy.optimize.minimize, default='Powell' Examples ======== >>> import spacepy.toolbox as tb >>> from scipy.stats import poisson >>> import matplotlib.pyplot as plt >>> import numpy as np >>> data = poisson.rvs(20, size=1000) >>> res = tb.poisson_fit(data) >>> print(res.x) 19.718000038769095 >>> xvals = np.arange(0, np.max(data)+5) >>> plt.hist(data, bins=xvals, normed=True) >>> plt.plot(xvals, poisson.pmf(xvals, np.round(res.x))) Returns ======= result : scipy.optimize.optimize.OptimizeResult Resulting fit results from scipy.optimize, answer is result.x, user should likely round. """ from scipy.optimize import minimize from scipy.stats import poisson def negLogLikelihood(params, data): """ the negative log-Likelihood-Function""" return (- np.sum(np.log(poisson.pmf(data, params[0])))) if initial is None: initial = np.median(data) ans = minimize(negLogLikelihood, # function to minimize x0=initial, # start value args=(data,), # additional arguments for function method=method, # minimization method, see docs ) return ans