viscid package¶

Subpackages¶

Module contents¶

A set of python modules that aid in plotting scientific data

Plotting depends on matplotlib and/or mayavi and file reading uses h5py and to read hdf5 / xdmf files.

Note

Modules in calculator and plot must be imported explicitly since they have side effects on import.

viscid.logger¶

a logging object whose verbosity can be set from the command line using :py:func`viscid.vutil.common_argparse`.

Type:	logging.Logger

class viscid.NOT_SPECIFIED[source]¶

Bases: object

default value; never instantiate and test with is

viscid.is_list_of_fields(lst)[source]¶: is a sequence a sequence of Field objects?

viscid.dataset_to_amr_grid(dset, template_skeleton=None)[source]¶

Try to divine AMR-ness from a Dataset

Parameters:	dset (`viscid.dataset.Dataset`) – Should be a collection of grids that could represet a series of patches in an AMR hierarchy. patch_list (`AMRTemplate`) – a guess at how patches are arranged for speedup. Allows us not to have to re-build the neighbor relationships if we’ve seen this patch layout before
Returns:	AMRGrid if dset looks like an AMR hierarchy, or dset itself otherwise

exception viscid.NonuniformFullArrayError[source]¶

Bases: exceptions.ValueError

trying to make uniform crds from nonuniform full arrays

viscid.arrays2crds(crd_arrs, crd_type='AUTO', crd_names='xyzuvw', **kwargs)[source]¶

make either uniform or nonuniform coordnates given full arrays

Parameters:	crd_arrs (array-like) – for n-dimensional crds, supply a list of n ndarrays crd_type (str) – passed to wrap_crds, should uniquely specify a type of coordinate object crd_names (iterable) – names of coordinates in the same order as crd_arrs. Should always be in xyz order.

viscid.wrap_crds(crdtype, clist, **kwargs)[source]¶

viscid.extend_arr(x, n=1, cell_fraction=1.0, full_arr=True, default_width=1e-05, width_arr=None)[source]¶: TODO: docstring

viscid.as_nvec(arr, ndim=4, init_vals=(0, 0, 0, 1))[source]¶

Extend vectors to ndim dimensions

The first dimension or arr that is the one that is extended. If the last dimension has shape >= 4, leave arr unchanged

Parameters:	arr (sequence) – array with shape (3,) or (3, N) for N vectors ndim (int) – extend arr to this many dimensions init_vals (sequence) – initialized values for new dimensions, note that init_vals[i] is the initialization for the ith dimension of the result. (0, 0, 0, 1) is useful for 4x4 rotation+translation matrices
Returns:	(ndim,) or (ndim, N) depending on input
Return type:	ndarray

viscid.make_rotation(theta, axis='z', rdim=3)[source]¶

Make rotation matrix of theta degrees around axis

Warning

The Hapgood paper uses the convention that rotations around y are right-handed, but rotations around x and z are left handed!

Parameters:	theta (float) – angle (degrees) axis (int, str) – one of (0, ‘x’, 1, ‘y’, 2, ‘z’) rdim (int) – 3 to make 3x3 matrices, or 4 to make 4x4 where mat[:, 3] == mat[3, :] == [0, 0, 0, 1]. 4x4 matrices are useful for performing translations.
Raises:	`ValueError` – On invalid axis

viscid.make_translation(v)[source]¶

Make a 4x4 translation matrix to act on [x, y, z, 1] vectors

Parameters:	v (sequence) – Translation vector, must have at least 3 entries
Returns:	mat will translate an [x, y, z, 1] vector, and imat will do the inverse translation
Return type:	(mat, imat)

viscid.get_rotation_wxyz(mat, check_inverse=True, quick=True)[source]¶

Find rotation axis and angle of a rotation matrix

Parameters:

Parameters:	mat (ndarray) – 3x3 rotation matrix with determinant +1 check_inverse (bool) – check if mat.T == inverse(mat) quick (bool) – Try a quicker, less well tested method first
Returns:	[w (rotation angle), ux, uy, uz] where u is the normalized axis vector
Return type:	ndarray

mat (ndarray) – 3x3 rotation matrix with determinant +1
check_inverse (bool) – check if mat.T == inverse(mat)
quick (bool) – Try a quicker, less well tested method first

Returns:

[w (rotation angle), ux, uy, uz] where u is the: normalized axis vector

Return type:

ndarray

viscid.as_crd_system(thing, default=<class 'viscid.cotr._NOT_GIVEN'>)[source]¶

Try to figure out a crd_system for thing

Parameters:	thing – Either a string with a crd_system abbreviation, or something that that can do thing.find_info(‘crd_system’) default (str) – fallback value
Returns:	crd system abbreviation, stripped and lower case
Return type:	str
Raises:	`TypeError` – if no default supplied and thing can’t be turned into a crd_system

class viscid.Cotr(time='1967-01-01', dip_tilt=None, dip_gsm=None, rdim=3, notilt1967=True)[source]¶

Bases: object

Transform geocentric vectors between crd systems @ a given UTC

GEOCENTRIC = ['gei', 'geo', 'gse', 'gsm', 'sm', 'mag', 'mhd']¶

HELIOCENTRIC = ['hea', 'hee', 'heeq']¶

Q_e¶

Q_g¶

XFORM_LOOKUP = {'gei->gei': (0,), 'geo->gei': (-1,), 'geo->geo': (0,), 'gse->gei': (-2,), 'gse->geo': (1, -2), 'gse->gse': (0,), 'gsm->gei': (-2, -3), 'gsm->geo': (1, -2, -3), 'gsm->gse': (-3,), 'gsm->gsm': (0,), 'hae->hae': 0, 'hee->hae': -11, 'hee->hee': 0, 'heeq->hae': -12, 'heeq->heeq': 0, 'mag->gei': (-1, -5), 'mag->geo': (-5,), 'mag->gse': (2, -1, -5), 'mag->gsm': (3, 2, -1, -5), 'mag->mag': (0,), 'mag->sm': (4, 3, 2, -1, -5), 'mhd->gei': (-2, 6), 'mhd->geo': (1, -2, 6), 'mhd->gse': (6,), 'mhd->gsm': (3, 6), 'mhd->mag': (5, 1, -2, 6), 'mhd->mhd': (0,), 'mhd->sm': (4, 3, 6), 'sm->gei': (-2, -3, -4), 'sm->geo': (1, -2, -3, -4), 'sm->gse': (-3, -4), 'sm->gsm': (-4,), 'sm->sm': (0,)}¶

dip_geolat¶

dip_geolat_rad¶

dip_geolon¶

dip_geolon_rad¶

dip_gsm¶: Angle psi between GSE-Z and dipole projected into GSE-YZ

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_mu¶: Dipole tilt angle mu (angle between GSM-Z and dipole)

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_psi¶: Angle psi between GSE-Z and dipole projected into GSE-YZ

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_tilt¶: Dipole tilt angle mu (angle between GSM-Z and dipole)

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

get_dipole_angles()[source]¶

Get rotation angles between dipole and GSE

Note

psi refers to the angle between GSE and GSM, while mu is the dipole tilt angle.

Returns:	The SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>
Return type:	(tilt, gsm)

get_dipole_moment(crd_system='gse', strength=32707.5292732387)[source]¶

Get Earth’s dipole moment

Parameters:	crd_system (str) – crd_system of resulting moment strength (float) – magnitude of dipole moment
Returns:	3-vector of dipole tilt in crd_system
Return type:	ndarray

get_rotation_wxyz(from_system, to_system)[source]¶

Get rotation axis and angle for a transformation

Parameters:

from_system (str) – abbreviation of crd_system
to_system (str) – abbreviation of crd_system

Returns:

[w, x, y, z] where w is the angle around the axis: vector [x, y, z]

Return type:

ndarray

get_transform_matrix(from_system, to_system)[source]¶

Get a transformation matrix from one system to another

Parameters:	from_system (str) – abbreviation of crd_system to_system (str) – abbreviation of crd_system
Returns:	self.rdim x self.rdim transformation matrix
Return type:	ndarray

sun_ecliptic_longitude¶

transform(from_system, to_system, vec)[source]¶

Transform a vector from one system to another

Parameters:

Parameters:	from_system (str) – abbreviation of crd_system to_system (str) – abbreviation of crd_system vec (ndarray) – should have shape (3, ) or (3, N) to transform N vectors at once
Returns:	vec in new crd system shaped either (3, ) or (3, N) mirroring the shape of vec
Return type:	ndarray

from_system (str) – abbreviation of crd_system
to_system (str) – abbreviation of crd_system
vec (ndarray) – should have shape (3, ) or (3, N) to transform N vectors at once

Returns:

vec in new crd system shaped either (3, ) or: (3, N) mirroring the shape of vec

Return type:

ndarray

viscid.as_cotr(thing=None, default=<class 'viscid.cotr._NOT_GIVEN'>)[source]¶

Try to make cotr into a Cotr instance

Inputs understood are:

None for North-South dipole
anything with a __cotr__ property
datetimes or similar, specifies the time for finding dip angles
mapping (dict or similar), passed to Cotr constructor as kwargs

Parameters:	thing (obj) – something to turn into a Cotr default (None, obj) – fallback value
Returns:	Cotr instance
Raises:	`TypeError` – if no default supplied and thing can’t be turned into a Cotr instance

viscid.cotr_transform(date_time, from_system, to_system, vec, notilt1967=True)[source]¶

Transform a vector from one system to another

Parameters:	date_time (str, datetime64) – datetime of transformation from_system (str) – abbreviation of crd_system to_system (str) – abbreviation of crd_system vec (ndarray) – should have shape (3, ) or (N, 3) to transform N vectors at once notilt1967 (bool) – is 1 Jan 1967 the special notilt time?
Returns:	vec in new crd system
Return type:	ndarray

viscid.get_dipole_moment(date_time, crd_system='gse', strength=32707.5292732387, notilt1967=True)[source]¶

Get Earth’s dipole moment at datetime

Parameters:	date_time (str, datetime64) – datetime to get dipole crd_system (str) – crd_system of resulting moment strength (float) – magnitude of dipole moment notilt1967 (bool) – is 1 Jan 1967 the special notilt time?
Returns:	3-vector of dipole tilt in crd_system
Return type:	ndarray

viscid.get_dipole_moment_ang(dip_tilt=0.0, dip_gsm=0.0, crd_system='gse', strength=32707.5292732387)[source]¶

Get dipole moment from tilt and gsm angles

Parameters:	dip_tilt (float) – if given, override the dipole tilt angle (mu degrees, positive points north pole sunward) dip_gsm (float) – if given, override the dipole gse->gsm angle (psi in degrees, positive points north pole duskward) crd_system (str) – coordinate system of the result strength (float) – magnitude of the dipole moment
Returns:	dipole moment [mx, my, mz]
Return type:	ndarray

viscid.get_dipole_angles(date_time, notilt1967=True)[source]¶

Get rotation angles between dipole and GSE

Parameters:	date_time (str, datetime64) – datetime to get dipole notilt1967 (bool) – is 1 Jan 1967 the special notilt time?

Note

psi refers to the angle between GSE and GSM, while mu is the dipole tilt angle.

Returns:	The SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>
Return type:	(tilt, gsm)

viscid.dipole_moment2cotr(m, crd_system='gse')[source]¶

Turn dipole moment vector into a Cotr instance

Parameters:	m (sequence) – dipole moment as sequence with length 3 crd_system (str) – crd_system of m
Returns:	Cotr instance

viscid.to_dataframe(collection, fld_names=None, selection=Ellipsis, time_sel=slice(None, None, None), time_col='time', datetime_col='datetime')[source]¶

Consolidate field collection into pandas dataframe

Parameters:	collection (sequence) – Can be one of (Field, List[Field], Dataset, Grid) fld_names (sequence, None) – grab specific fields by name, or None to grab all fields selection (selection) – optional spatial selection time (selection) – optional time selection
Returns:	pandas.DataFrame

viscid.from_dataframe(frame, crd_cols=None, time_col='time', datetime_col='datetime')[source]¶

Make either a DatasetTemporal or Grid from pandas dataframe

Parameters:	frame (pandas.DataFrame) – frame to parse crd_cols (List[Str], None) – list of column names for coordinates time_col (str) – column name of times datetime_col (str) – column name of datetimes
Returns:	DatasetTemporal or Grid
Raises:	`ValueError` – if only 1 row given and crd_cols is None

viscid.guess_dipole_moment(b, r=2.0, strength=32707.5292732387, cap_angle=40, cap_ntheta=121, cap_nphi=121, plot=False)[source]¶: guess dipole moment from a B field

viscid.make_dipole(m=(0, 0, -32707.5292732387), strength=None, l=None, h=None, n=None, center='cell', dtype='f8', twod=False, nonuniform=False, crd_system='gse', name='b')[source]¶: Generate a dipole field with magnetic moment m [x, y, z]

viscid.fill_dipole(B, m=(0, 0, -32707.5292732387), strength=None, mask=None)[source]¶

set B to a dipole with magnetic moment m

Parameters:	B (Field) – Field to fill with a dipole m (ndarray, or datetime64-like) – Description strength (float) – if given, rescale the dipole moment even if it was given explicitly mask (Field) – boolean field as mask, B will be filled where the mask is True
Returns:	B
Return type:	Field

viscid.calc_dip(pts, m=(0, 0, -32707.5292732387), strength=None, crd_system='gse', dtype=None)[source]¶

Calculate a dipole field at various points

Parameters:	pts (ndarray) – Nx3 array of points at which to calculate the dipole. Should use the same crd system as m m (sequence, datetime) – dipole moment strength (None, float) – If given, rescale m to this magnitude crd_system (str) – Something from which cotr can divine the coordinate system for both pts and m. This is only used if m is given as a datetime and we need to figure out the dipole moment at a given time in a given crd system dtype (str, np.dtype) – dtype of the result, defaults to the same datatype as pts
Returns:	Nx3 dipole field vectors for N points
Return type:	ndarray

viscid.set_in_region(a, b, alpha=1.0, beta=1.0, mask=None, out=None)[source]¶: set ret = alpha * a + beta * b where mask is True

viscid.make_spherical_mask(fld, rmin=0.0, rmax=None, rsq=None)[source]¶: make a mask that is True between rmin and rmax

viscid.xyz2lsrlp(pts, cotr=None, crd_system='gse')[source]¶

Ceovert x, y, z -> l-shell, r, lambda, phi [sm coords]

r, theta, phi = viscid.cart2sph(pts in x, y, z)

lambda = 90deg - theta

r = L cos^2(lambda)

Parameters:	pts (ndarray) – 3xN for N (x, y, z) points cotr (None) – if given, use cotr to perform mapping to / from sm crd_system (str) – crd system of pts
Returns:	4xN array of N (l-shell, r, lamda, phi) points
Return type:	ndarray

viscid.dipole_map(pts, r=1.0, cotr=None, crd_system='gse', as_spherical=False)[source]¶

Map pts along an ideal dipole to radius r

lambda = 90deg - theta; r = L cos^2(lambda)

cos^2(lambda) = cos^2(lambda_0) * (r / r_0)

Parameters:

Parameters:	pts (ndarray) – 3xN for N (x, y, z) points r (float) – radius to map to cotr (None) – if given, use cotr to perform mapping to / from sm crd_system (str) – crd system of pts as_spherical (bool) – if True, then the return array is (t, theta, phi) with theta in the range [0, 180] and phi [0, 360] (in degrees)
Returns:	3xN array of N (x, y, z) points all at a distance r_mapped from the center of the dipole
Return type:	ndarray

pts (ndarray) – 3xN for N (x, y, z) points
r (float) – radius to map to
cotr (None) – if given, use cotr to perform mapping to / from sm
crd_system (str) – crd system of pts
as_spherical (bool) – if True, then the return array is (t, theta, phi) with theta in the range [0, 180] and phi [0, 360] (in degrees)

Returns:

3xN array of N (x, y, z) points all at a distance: r_mapped from the center of the dipole

Return type:

ndarray

viscid.dipole_map_value(fld, pts, r=1.0, fillna=None, cotr=None, crd_system=None, interp_kind='linear')[source]¶

Map values assuming they’re constant along ideal dipole lines

Parameters:	fld (Field) – values to interpolate onto the mapped pts pts (ndarray) – 3xN for N (x, y, z) points that will be mapped r (float) – radius of resulting map cotr (None) – if given, use cotr to perform mapping in sm crd_system (str) – crd system of pts interp_kind (str) – how to interpolate fld onto source points
Returns:	ndarray of mapped values, one for each of the N points
Return type:	ndarray

viscid.make_animation(movie_fname, prefix, framerate=5, qscale=2, keep=False, args=None, frame_idx_fmt='_%06d', program='ffmpeg', yes=False)[source]¶: make animation by calling program (only ffmpeg works for now) using args, which is a namespace filled by the argparse options from add_animate_arguments. Plots are expected to be named ${args.prefix}_000001.png where the number is in order from 1 up

viscid.meshlab_convert(fname, fmt='dae', quiet=True)[source]¶

Run meshlabserver to convert 3D mesh files

Uses MeshLab, which is a great little program for playing with 3D meshes. The best part is that OS X’s Preview can open the COLLADA (*.dae) format. How cool is that?

Parameters:	fname (str) – file to convert fmt (str) – extension of result, defaults to COLLADA format quiet (bool) – redirect output to `os.devnull`
Returns:	None

viscid.arrays2field(crd_arrs, dat_arr, name='NoName', center=None, crd_type=None, crd_names='xyzuvw')[source]¶

Turn arrays into fields so they can be used in viscid.plot, etc.

This is a convenience function that takes care of making coordnates and the like. If the default behavior doesn’t work for you, you’ll need to make your own coordnates and call viscid.field.wrap_field().

Parameters:	crd_arrs (list of ndarrays) – xyz list of ndarrays that describe the node centered coordnates of the field dat_arr (ndarray) – data with len(crd_arrs) or len(crd_arrs) + 1 dimensions name (str) – some name center (str, None) – If not None, translate field to this centering (node or cell)

viscid.dat2field(dat_arr, name='NoName', fldtype='scalar', center=None, layout='flat')[source]¶

Makes np.arange coordnate arrays and calls arrays2field

Parameters:	dat_arr (ndarray) – data name (str) – name of field fldtype (str, optional) – ‘scalar’ / ‘vector’ center (str, None) – If not None, translate field to this centering (node or cell) layout (TYPE, optional) – Description

viscid.full(crds, fill_value, dtype='f8', name='NoName', center='cell', layout='flat', nr_comps=0, crd_type=None, crd_names='xyzuvw', **kwargs)[source]¶

Analogous to numpy.full

Parameters:

Parameters:	crds (Coordinates, list, or tuple) – Can be a coordinates object. Can also be a list of ndarrays describing coordinate arrays. Or, if it’s just a list or tuple of integers, those integers are taken to be the nz,ny,nx shape and the coordinates will be fill with `np.arange()`. fill_value (number, None) – Initial value of array. None indicates uninitialized (i.e., numpy.empty) dtype (optional) – some way to describe numpy dtype of data name (str) – a way to refer to the field programatically center (str, optional) – cell or node, there really isn’t support for edge / face yet layout (str, optional) – how data is stored, is in “flat” or “interlaced” (interlaced == AOS) nr_comps (int, optional) – for vector fields, nr of components **kwargs – passed through to Field constructor

crds (Coordinates, list, or tuple) – Can be a coordinates object. Can also be a list of ndarrays describing coordinate arrays. Or, if it’s just a list or tuple of integers, those integers are taken to be the nz,ny,nx shape and the coordinates will be fill with np.arange().
fill_value (number, None) – Initial value of array. None indicates uninitialized (i.e., numpy.empty)
dtype (optional) – some way to describe numpy dtype of data
name (str) – a way to refer to the field programatically
center (str, optional) – cell or node, there really isn’t support for edge / face yet
layout (str, optional) – how data is stored, is in “flat” or “interlaced” (interlaced == AOS)
nr_comps (int, optional) – for vector fields, nr of components
**kwargs – passed through to Field constructor

viscid.empty(crds, dtype='f8', name='NoName', center='cell', layout='flat', nr_comps=0, **kwargs)[source]¶

Analogous to numpy.empty

Returns:	new uninitialized `Field`

Returns:	new `Field` initialized to 0

Returns:	new `Field` initialized to 1

Parameters:	fld – field to get coordinates / metadata from fill_value (number, None) – initial value, or None to leave data uninitialized name – name for this field **kwargs – passed through to `Field` constructor
Returns:	new `Field`

Returns:	new `Field` filled with zeros

Returns:	new `Field` filled with ones

Parameters:	name (str) – name for the vector field fldlist – list of `ScalarField` **kwargs – passed to `VectorField` constructor
Returns:	A new `VectorField`.

Parameters:	data – Some data container, most likely a `numpy.ndarray` crds (Coordinates) – coordinates that describe the shape / grid of the field fldtype (str) – ‘scalar’ / ‘Vector’ name (str) – a way to refer to the field programatically **kwargs – passed through to `Field` constructor
Returns:	A `Field` instance.

Parameters:	vfile – a vFile object that we can call iter_times on time_slice – string, slice notation, like 1200:2400:1 initial_seeds – any SeedGen object plot_function – function that is called each time step, arguments should be exactly: (i [int], grid, v [Vector Field], v_lines [result of streamline trace], root_seeds [SeedGen]) stream_opts – must have ds0 and max_length, maxit will be automatically calculated
Returns:	root points after following the fluid

Parameters:	fld (Field) – some field to transform order (tuple, optional) – Desired axes of the result units (str, optional) – units of the result (deg/rad)
Returns:	a field with (lon, lat) coordinates. The data will be loaded into memory if not already.
Return type:	Viscid.Field

Returns:	a field with (lon, lat) coordinates. The data will be loaded into memory if not already.
Return type:	Viscid.Field

Parameters:	fld (Field) – Some scalar or vector field bounding_lat (float, optional) – Used to slice the field, i.e., gives data from the pole to bounding_lat degrees equator-ward of the pole hemisphere (str, optional) – ‘north’ or ‘south’
Returns:	a field that can be mapfield plotted on polar axes
Return type:	Field
Raises:	`ValueError` – on bad hemisphere

Parameters:	pts (ndarray) – A 2xN array of N points where the spatial dimensions encode phi, theta, lat, or lon. pts_axes (sequence) – sequence of strings that say what each axis of pts encodes, i.e., (‘theta’, ‘phi’). pts_unit (str, optional) – units of pts (‘deg’ or ‘rad’) hemisphere (str, optional) – ‘north’ or ‘south’
Returns:	2xN array of N points in radians where the axes are (‘lon’, ‘lat’), such that they can be given straight to matplotlib.pyplot.plot()
Return type:	ndarray
Raises:	`ValueError` – On bad pts_axes

Parameters:	fld (Field) – Description order (sequence) – target coordinates crd_mapping (dict) – summarizes functions that go from base -> target units (str, optional) – Additional info if you need to convert units too
Raises:	`RuntimeError` – If no mapping is found to go base -> target

Parameters:	arr (sequence) – shaped (3, ) or (3, N) for N xyz points deg (bool) – if result should be in degrees
Returns:	shaped (3, ) or (3, N) r, theta, phi points in radians by default
Return type:	ndarray

Parameters:	arr (sequence) – shaped (3, ) or (3, N) for N r, theta, phi points in radians by default deg (bool) – if arr is in dergrees
Returns:	shaped (2, ) or (2, N) xyz
Return type:	ndarray

Parameters:	arr (sequence) – shaped (3, ) or (3, N) for N r, theta, phi points in radians by default deg (bool) – if arr is in dergrees
Returns:	shaped (3, ) or (3, N) xyz
Return type:	ndarray

Parameters:	arr (sequence) – shaped (2, ) or (2, N) for N r, theta, phi points in radians by default deg (bool) – if arr is in dergrees
Returns:	shaped (3, ) or (3, N) xyz
Return type:	ndarray

See also

Use chunk_iterator() to chunk up iterators

Example

>>> it1, it2, it3 = chunk_list(range(8), 3)
>>> it1 == range(0, 3)  # 3 vals
True
>>> it2 == range(3, 6)  # 3 vals
True
>>> it3 == range(6, 8)  # 2 vals
True

viscid.chunk_slices(nel, nchunks, size=None)[source]¶

Make continuous chunks

Get the slice info (can be unpacked and passed to the slice builtin as in slice(*ret[i])) for nchunks contiguous chunks in a list with nel elements

Parameters:	nel – how many elements are in one pass of the original list nchunks – how many chunks to make size – if given, set nchunks such that chunks have about ‘size’ elements
Returns:	a list of (start, stop) tuples with length nchunks

Example

>>> sl1, sl2 = chunk_slices(5, 2)
>>> sl1 == (0, 3)  # 3 vals
True
>>> sl2 == (3, 5)  # 2 vals
True

viscid.chunk_interslices(nchunks)[source]¶

Make staggered chunks

Similar to chunk_slices, but pick every nth element instead of getting a contiguous patch for each chunk

Parameters:	nchunks – how many chunks to make
Returns:	a list of (start, stop, step) tuples with length nchunks

Example

>>> chunk_slices(2) == [(0, None, 2), (1, None, 2)]
True

viscid.chunk_sizes(nel, nchunks, size=None)[source]¶

For chunking up lists, how big is each chunk

Parameters:	nel – how many elements are in one pass of the original list nchunks – is inferred from the length of iter_list size – if given, set nchunks such that chunks have about ‘size’ elements
Returns:	an ndarray of the number of elements in each chunk, this should be the same for chunk_list, chunk_slices and chunk_interslices

Example

>>> nel1, nel2 = chunk_sizes(5, 2)
>>> nel1 == 2
True
>>> nel2 == 3
True

viscid.map(nr_procs, func, args_iter, args_kw=None, timeout=100000000.0, daemonic=True, threads=False, pool=None, force_subprocess=False)[source]¶

Just like subprocessing.map?

same as map_async(), except it waits for the result to be ready and returns it

Note

When using threads, this is WAY faster than map_async since map_async uses the builtin python ThreadPool. I have no idea why that’s slower than making threads by hand.

viscid.map_async(nr_procs, func, args_iter, args_kw=None, daemonic=True, threads=False, pool=None)[source]¶

Wrap python’s map_async

This has some utility stuff like star passthrough

Run func on nr_procs with arguments given by args_iter. args_iter should be an iterable of the list of arguments that can be unpacked for each invocation. kwargs are passed to func as keyword arguments

Returns:	(tuple) (pool, multiprocessing.pool.AsyncResult)

Note

When using threads, this is WAY slower than map since map_async uses the builtin python ThreadPool. I have no idea why that’s slower than making threads by hand.

Note: daemonic can be set to False if one needs to spawn child: processes in func, BUT this could be vulnerable to creating an undead army of worker processes, only use this if you really really need it, and know what you’re doing

Example

>>> func = lambda i, letter: print i, letter
>>> p, r = map_async(2, func, itertools.izip(itertools.count(), 'abc'))
>>> r.get(1e8)
>>> p.join()
>>> # the following is printed from 2 processes
0 a
1 b
2 c

viscid.angle2rot(angles, axes='zyx', unit='rad')[source]¶

Euler angles (transform-order) -> rotation matrix (extrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

Example

>>> angle2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.eul2rot(angles, axes='zyx', unit='rad')[source]¶

Euler angles (multiplication-order) -> rotation matrix (extrinsic)

Rotations are applied in multiplication-order, which means the last axis given is the first transform applied. In other words:

>>> R = (R(angles[..., 0], axis[..., 0]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 2], axis[..., 2]))

This function is equivalent (up to a few machine epsilon) to Matlab’s eul2rotm.

Example

>>> eul2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in multiplication-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.angle2dcm(angles, axes='zyx', unit='rad')[source]¶

Euler angles (transform-order) -> rotation matrix (intrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

This function is equivalent (up to a few machine epsilon) to Matlab’s angle2dcm.

Example

>>> angle2dcm([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, -1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.eul2dcm(angles, axes='zyx', unit='rad')[source]¶

Euler angles (multiplication-order) -> rotation matrix (intrinsic)

Rotations are applied in multiplication-order, which means the last axis given is the first transform applied. In other words:

>>> R = (R(angles[..., 0], axis[..., 0]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 2], axis[..., 2]))

Example

>>> eul2dcm([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, -1, 0]

Parameters:

angles (sequence) – Euler angles in multiplication-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.rot2angles(R, axes='zyx', unit='rad', bad_matrix='warn')¶

Rotation matrix (extrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2rot() for more details about axes order.

viscid.rot2eul(R, axes='zyx', unit='rad', bad_matrix='warn')[source]¶

Rotation matrix (extrinsic) -> Euler angles (multiplication-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

rot2eul() for more details about axes order.

viscid.dcm2angles(R, axes='zyx', unit='rad', bad_matrix='warn')¶

Rotation matrix (intrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2dcm() for more details about axes order.

viscid.dcm2eul(R, axes='zyx', unit='rad', bad_matrix='warn')[source]¶

Rotation matrix (intrinsic) -> Euler angles (multiplication-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

eul2dcm() for more details about axes order.

viscid.rot2quat(R)[source]¶

Rotation matrix (extrinsic) -> quaternion

If this quaternion library is imported, then the results are presented with dtype quaternion, otherwise, as dtype numpy.double (array-of-struct, [scalar, vec0, vec1, vec2]).

Parameters:	R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
Returns:	quaternions with dtype quaternion with shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Return type:	ndarray

viscid.dcm2quat(R)[source]¶

Rotation matrix (intrinsic) -> quaternion

If this quaternion library is imported, then the results are presented with dtype quaternion, otherwise, as dtype numpy.double (array-of-struct, [scalar, vec0, vec1, vec2]).

Parameters:	R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
Returns:	quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Return type:	ndarray

viscid.quat2rot(q)[source]¶

Quaternion -> rotation matrix (extrinsic)

Parameters:	q (ndarray) – quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Returns:	orthonormal rotation matrix with shape (3, 3) or (Nmatrices, 3, 3)
Return type:	ndarray

viscid.quat2dcm(q)[source]¶

Quaternion -> rotation matrix (intrinsic)

Parameters:	q (ndarray) – quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Returns:	rotation matrix with shape [Ndim, Ndim] or [Nmatrices, Ndim, Ndim] depending on the shape of `q`
Return type:	ndarray

viscid.rotmul(*rmats)[source]¶

Multiply rotation matrices together

Parameters:	rmats (ndarrays*) – rotation matrices with shape (3, 3) or (nmats, 3, 3)
Returns:	rotation matrix with shape (3, 3) or (nmats, 3, 3)
Return type:	ndarray

viscid.rotate(vec, *rmats)[source]¶

Rotate R3 vector(s) by one or more matrices

Parameters:	vec (ndarray) – R3 vectors with shape (3,) or (nvecs, 3) rmats (ndarrays*) – rotation matrices with shape (3, 3) or (nmats, 3, 3)
Returns:	rotated vectors with shape (3,) or (nvecs, 3)
Return type:	ndarray

viscid.quatmul(*quats)[source]¶

Multiply quaternions together

Parameters:	quats (ndarrays*) – quaternions with dtype quaternion and shape () or (nquat,); or with dtype np.double and shape (4,) or (nquat, 4)
Returns:	with dtype quaternion and shape () or (nquat,); or with dtype np.double and shape (4,) or (nquat, 4)
Return type:	ndarray

viscid.quat_rotate(vec, *quats)[source]¶

Rotate R3 vector(s) by one or more quaternions

Parameters:	vec (ndarray) – R3 vectors with shape (3,) or (nvecs, 3) quats (ndarrays*) – quaternions with dtype quaternion and shape () or (nquat,); or with dtype np.double and shape (4,) or (nquat, 4)
Returns:	rotated vectors with shape (3,) or (nvecs, 3)
Return type:	ndarray

viscid.wxyz2rot(angle, vector, unit='rad')[source]¶

Make a matrix (matrices) that rotates around vector by angle

Parameters:

angle (sequence) – angle(s) of rotation(s) with shape (), (1,), or (Nmatrices)
vector (sequence) – 3d vector(s) that is (are) axis of rotation(s) with shape (3,) or (Nmatrices, 3)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

orthonormal rotation matrix with shape (3, 3): or (Nmatrices, 3, 3)

Return type:

ndarray

viscid.axang2rot(axang, unit='rad')[source]¶

Make a matrix (matrices) that rotates around vector by angle

Parameters:

axang (ndarray) – axis(es) / angle(s) of rotation(s) put together like {x, y, z, angle}; shape should be (4,), or (Nmatrices, 4)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

orthonormal rotation matrix with shape (3, 3): or (Nmatrices, 3, 3)

Return type:

ndarray

viscid.rot2wxyz(R, unit='rad', quick=True, bad_matrix='warn')[source]¶

Find axis / angle representation of rotation matrix (matrices)

Parameters:

R (ndarray) – Rotation matrix with shape (3, 3) determinant +1, or matrices with shape (Nmatrices, 3, 3)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)
quick (bool) – Try a quicker, less well tested method
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

w is rotation angle with shape () or (Nmatrices),: and xyz are normalized rotation axes with shape (3,) or (Nmatrices, 3)

Return type:

(w, xyz)

viscid.rot2axang(R, unit='rad', quick=True, bad_matrix='warn')[source]¶

Find axis / angle representation of rotation matrix (matrices)

Parameters:

R (ndarray) – Rotation matrix with shape (3, 3) determinant +1, or matrices with shape (Nmatrices, 3, 3)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)
quick (bool) – Try a quicker, less well tested method
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

{x, y, z, angle} axis / angle values with shape (4,): or (Nmatrices, 4)

Return type:

(ndarray)

viscid.wxyz2quat(angle, vector, unit='rad')[source]¶

(angle, axis) representation -> quaternion

Parameters:

angle (sequence) – angle(s) of rotation(s) with shape (), (1,), or (Nmatrices)
vector (sequence) – 3d vector(s) that is (are) axis of rotation(s) with shape (3,) or (Nmatrices, 3)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

quaternions with dtype quaternion and shape: (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)

Return type:

ndarray

viscid.axang2quat(axang, unit='rad')[source]¶

axis-angle representation -> quaternion

Parameters:

axang (ndarray) – axis(es) / angle(s) of rotation(s) put together like {x, y, z, angle}; shape should be (4,), or (Nmatrices, 4)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

quaternions with dtype quaternion and shape: (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)

Return type:

ndarray

viscid.quat2wxyz(q, unit='rad')[source]¶

quaternion -> (angle, axis) representation

Parameters:

q (ndarray) – quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)

Returns:

w is rotation angle with shape () or (Nmatrices),: and xyz are normalized rotation axes with shape (3,) or (Nmatrices, 3)

Return type:

(w, xyz)

viscid.quat2axang(q, unit='rad')[source]¶

quaternion -> axis-angle representation

Parameters:

q (ndarray) – quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)

Returns:

{x, y, z, angle} axis / angle values with shape (4,): or (Nmatrices, 4)

Return type:

(ndarray)

viscid.convert_angle(angle, from_unit, to_unit)[source]¶

convert angle(s) from rad/deg to rad/deg

Parameters:	angle (float, sequence) – angle in from_unit from_unit (str) – unit of angle, one of (‘rad’, ‘deg’) to_unit (str) – unit of result, one of (‘rad’, ‘deg’)
Returns:	angle converted to to_unit
Return type:	float or ndarray

viscid.check_orthonormality(mat, bad_matrix='warn')[source]¶

Check that a matrix or stack of matrices is orthonormal

Parameters:	mat (ndarray) – A matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim] bad_matrix (str) – What to do if `numpy.det(R) != 1.0` - can be one of (‘ignore’, ‘warn’, ‘raise’)

viscid.e_rot(theta, axis='z', unit='rad')[source]¶

Make elementary rotation matrices (extrinsic) of theta around axis

Example

>>> e_rot(numpy.pi / 2, 'z') @ [1, 0, 0] = [0, 1, 0]

Parameters:	theta (float, sequence) – angle or angles axis (int, str) – one of (0, ‘x’, 1, ‘y’, 2, ‘z’) unit (str) – unit of theta, one of (‘deg’, ‘rad’)
Returns:	rotation matrix with shape [ndim, ndim] or matrices with shape [Nmatrices, Ndim, Ndim]
Return type:	ndarray
Raises:	`ValueError` – invalid axis

viscid.e_dcm(theta, axis='z', unit='rad')[source]¶

Make elementary rotation matrices (intrinsic) of theta around axis

Example

>>> e_dcm(numpy.pi / 2, 'z') @ [1, 0, 0] = [0, -1, 0]

Parameters:	theta (float, sequence) – angle or angles axis (int, str) – one of (0, ‘x’, 1, ‘y’, 2, ‘z’) unit (str) – unit of theta, one of (‘deg’, ‘rad’)
Returns:	rotation matrix with shape [ndim, ndim] or matrices with shape [Nmatrices, Ndim, Ndim]
Return type:	ndarray
Raises:	`ValueError` – invalid axis

viscid.angle_between(a, b, unit='rad')[source]¶

Get angle(s) between two (sets of) vectors

Parameters:	a (sequence) – first vector, shape [Ndim] or [Nvectors, Ndim] b (sequence) – second vector, shape [Ndim] or [Nvectors, Ndim] unit (str) – unit of result, one of (‘deg’, ‘rad’)
Returns:	angle(s)
Return type:	float or ndarray

viscid.a2b_rot(a, b, roll=0.0, unit='rad', new_x=None)[source]¶

Make a matrix that rotates vector a to vector b

Parameters:

a (ndarray) – starting vector with shape [Ndim] or [Nvectors, Ndim]
b (ndarray) – destination vector with shape [Ndim] or [Nvectors, Ndim]
roll (float) – angle of roll around the b vector
unit (str) – unit of roll, one of (‘deg’, ‘rad’)
new_x (ndarray) – If given, then roll is set such that the x axis (phi = 0) new_x is projected in the plane perpendicular to origin-p1

Returns:

orthonormal rotation matrix with shape [Ndim, Ndim]: or [Nvectors, Ndim, Ndim]

Return type:

ndarray

viscid.symbolic_e_dcm(axis='z', angle=None)[source]¶: Make elementary matrix that rotate axes, not vectors (sympy)

viscid.symbolic_e_rot(axis='z', angle=None)[source]¶: Make elementary matrix that rotate vectors, not the axes (sympy)

viscid.symbolic_rot(axes='zyx', angles=None)[source]¶: Make a symbolic rotation matrix (sympy)

viscid.symbolic_dcm(axes='zyx', angles=None)[source]¶: Make a symbolic direction cosine matrix (sympy)

viscid.convert_angles(angle, from_unit, to_unit)¶

convert angle(s) from rad/deg to rad/deg

Parameters:	angle (float, sequence) – angle in from_unit from_unit (str) – unit of angle, one of (‘rad’, ‘deg’) to_unit (str) – unit of result, one of (‘rad’, ‘deg’)
Returns:	angle converted to to_unit
Return type:	float or ndarray

viscid.angles2rot(angles, axes='zyx', unit='rad')¶

Euler angles (transform-order) -> rotation matrix (extrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

Example

>>> angle2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.angles2dcm(angles, axes='zyx', unit='rad')¶

Euler angles (transform-order) -> rotation matrix (intrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

This function is equivalent (up to a few machine epsilon) to Matlab’s angle2dcm.

Example

>>> angle2dcm([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, -1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.rot2angles(R, axes='zyx', unit='rad', bad_matrix='warn')

Rotation matrix (extrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2rot() for more details about axes order.

viscid.dcm2angles(R, axes='zyx', unit='rad', bad_matrix='warn')

Rotation matrix (intrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2dcm() for more details about axes order.

viscid.eul2rotm(angles, axes='zyx', unit='rad')¶

Euler angles (multiplication-order) -> rotation matrix (extrinsic)

Rotations are applied in multiplication-order, which means the last axis given is the first transform applied. In other words:

>>> R = (R(angles[..., 0], axis[..., 0]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 2], axis[..., 2]))

This function is equivalent (up to a few machine epsilon) to Matlab’s eul2rotm.

Example

>>> eul2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in multiplication-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.angle2rotm(angles, axes='zyx', unit='rad')¶

Euler angles (transform-order) -> rotation matrix (extrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

Example

>>> angle2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.angles2rotm(angles, axes='zyx', unit='rad')¶

Euler angles (transform-order) -> rotation matrix (extrinsic)

Rotations are applied in transform-order, which means the first axis given is the first transform applied. In other words, the matrix multiply is in reverse order, i.e.:

>>> R = (R(angles[..., 2], axis[..., 2]) @
>>>      R(angles[..., 1], axis[..., 1]) @
>>>      R(angles[..., 0], axis[..., 0]))

Example

>>> angle2rot([numpy.pi / 2, 0, 0], 'zyx') @ [1, 0, 0] = [0, 1, 0]

Parameters:

angles (sequence) – Euler angles in transform-order; can have shape [Ndim] or [Nmatrices, Ndim] to make stacked transform matrices
axes (sequence, str) – rotation axes in transform-order
unit (str) – unit of angles, one of (‘deg’, ‘rad’)

Returns:

rotation matrix with shape [Ndim, Ndim] or: [Nmatrices, Ndim, Ndim] depending on the shape of angles

Return type:

ndarray

viscid.rotm2eul(R, axes='zyx', unit='rad', bad_matrix='warn')¶

Rotation matrix (extrinsic) -> Euler angles (multiplication-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in multiplication-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

rot2eul() for more details about axes order.

viscid.rotm2angle(R, axes='zyx', unit='rad', bad_matrix='warn')¶

Rotation matrix (extrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2rot() for more details about axes order.

viscid.rotm2angles(R, axes='zyx', unit='rad', bad_matrix='warn')¶

Rotation matrix (extrinsic) -> Euler angles (transform-order)

Parameters:

R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
axes (sequence, str) – rotation axes in transform-order
unit (str) – Unit of angles, one of (‘deg’, ‘rad’)
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

Euler angles with shape [Ndim] or [Nmatrices, Ndim]: depending on the input. angles[:, i] always corresponds to axes[i].

Return type:

ndarray

See also

angle2rot() for more details about axes order.

viscid.rotm2quat(R)¶

Rotation matrix (extrinsic) -> quaternion

If this quaternion library is imported, then the results are presented with dtype quaternion, otherwise, as dtype numpy.double (array-of-struct, [scalar, vec0, vec1, vec2]).

Parameters:	R (ndarray) – A rotation matrix with shape [Ndim, Ndim] or a stack of matrices with shape [Nmatrices, Ndim, Ndim]
Returns:	quaternions with dtype quaternion with shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Return type:	ndarray

viscid.quat2rotm(q)¶

Quaternion -> rotation matrix (extrinsic)

Parameters:	q (ndarray) – quaternions with dtype quaternion and shape (Nmatrices,) or (); or with dtype numpy.double and shape (Nmatrices, 4) or (4)
Returns:	orthonormal rotation matrix with shape (3, 3) or (Nmatrices, 3, 3)
Return type:	ndarray

viscid.a2b_rotm(a, b, roll=0.0, unit='rad', new_x=None)¶

Make a matrix that rotates vector a to vector b

Parameters:

a (ndarray) – starting vector with shape [Ndim] or [Nvectors, Ndim]
b (ndarray) – destination vector with shape [Ndim] or [Nvectors, Ndim]
roll (float) – angle of roll around the b vector
unit (str) – unit of roll, one of (‘deg’, ‘rad’)
new_x (ndarray) – If given, then roll is set such that the x axis (phi = 0) new_x is projected in the plane perpendicular to origin-p1

Returns:

orthonormal rotation matrix with shape [Ndim, Ndim]: or [Nvectors, Ndim, Ndim]

Return type:

ndarray

viscid.axang2rotm(axang, unit='rad')¶

Make a matrix (matrices) that rotates around vector by angle

Parameters:

axang (ndarray) – axis(es) / angle(s) of rotation(s) put together like {x, y, z, angle}; shape should be (4,), or (Nmatrices, 4)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

orthonormal rotation matrix with shape (3, 3): or (Nmatrices, 3, 3)

Return type:

ndarray

viscid.rotm2axang(R, unit='rad', quick=True, bad_matrix='warn')¶

Find axis / angle representation of rotation matrix (matrices)

Parameters:

R (ndarray) – Rotation matrix with shape (3, 3) determinant +1, or matrices with shape (Nmatrices, 3, 3)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)
quick (bool) – Try a quicker, less well tested method
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

{x, y, z, angle} axis / angle values with shape (4,): or (Nmatrices, 4)

Return type:

(ndarray)

viscid.wxyz2rotm(angle, vector, unit='rad')¶

Make a matrix (matrices) that rotates around vector by angle

Parameters:

angle (sequence) – angle(s) of rotation(s) with shape (), (1,), or (Nmatrices)
vector (sequence) – 3d vector(s) that is (are) axis of rotation(s) with shape (3,) or (Nmatrices, 3)
unit (str) – unit of angle, one of (‘deg’, ‘rad’)

Returns:

orthonormal rotation matrix with shape (3, 3): or (Nmatrices, 3, 3)

Return type:

ndarray

viscid.rotm2wxyz(R, unit='rad', quick=True, bad_matrix='warn')¶

Find axis / angle representation of rotation matrix (matrices)

Parameters:

R (ndarray) – Rotation matrix with shape (3, 3) determinant +1, or matrices with shape (Nmatrices, 3, 3)
unit (str) – unit of angle in results, one of (‘deg’, ‘rad’)
quick (bool) – Try a quicker, less well tested method
bad_matrix (str) – What to do if numpy.det(R) != 1.0 - can be one of (‘ignore’, ‘warn’, ‘raise’)

Returns:

w is rotation angle with shape () or (Nmatrices),: and xyz are normalized rotation axes with shape (3,) or (Nmatrices, 3)

Return type:

(w, xyz)

viscid.to_seeds(pts)[source]¶

Try to turn anything into a set of seeds

Parameters:	points (ndarray or list) – should look like something that np.array(points) can turn to an Nx3 array of xyz points. This can be 3xN so long as N != 3.

class viscid.SeedGen(cache=False, dtype=None, fill_holes=True)[source]¶

Bases: object

All about seeds

These objects are good for defining the root points of streamlines or locations for interpolation.

For Developers
- Mandatory
  
  Subclasses must override
  
  get_nr_points(): get total number of seeds
  
  get_uv_shape(): get the shape of the ndarray of uv mesh coordinates IF there exists a 2D representation of the seeds
  
  get_local_shape(): get the shape of the ndarray of local coordinates
  
  uv_to_local(): transform 2d representation to local representation
  
  to_3d(): transform an array in local coordinates to 3D space. Should return a 3xN ndarray.
  
  to_local(): transform a 3xN ndarray to an ndarray of coordinates in local space. It’s ok to raise NotImplementedError.
  
  _make_local_points(): Make an ndarray that is used in _make_3d_points()
  
  _make_uv_axes(): Make axes for 2d mesh representation, i.e., matches uv_shape
  
  _make_local_axes(): Make a tuple of arrays with lengths that match local_shape. It’s ok to raise NotImplementedError if the local representation has no axes.
- Optional
  
  Subclasses may override
  
  _make_3d_points(): By default, this is a combination of _make_local_points() and to_3d(). It can be overridden if there is a more efficient way te get points in 3d.
  
  iter_points(): If a seed generator can be made lazy, then this should return an iterator that yields (x, y, z) points.
  
  get_rotation(): This should return an orthonormal rotation matrix if there is some rotation required to go from local coordinates to 3D space.
  
  as_mesh()
  
  as_uv_coordinates()
  
  as_local_coordinates()
  
  wrap_field()
- Off Limits
  
  Don’t override unless you know what you’re doing
  
  get_points()
For Users
- 3D
  
  These attributes give the seed points in 3D (xyz) space
  
  nr_points
  
  get_points()
  
  iter_points()
  
  to_3d()
- Generator Specific
  
  There is also the concept of local points, which are points that are meaningful for the specific SeedGen subclass. For example, for the Sphere subclass, local points are labled with (theta, phi) on the surface of the sphere. As such, wrap_field() can be used to facilitate wrapping the result of trilin_interp() for easy plotting.
  
  uv_shape
  
  local_shape
  
  uv_to_local()
  
  to_local()
  
  as_mesh()
  
  as_uv_coordinates()
  
  as_local_coordinates()
  
  wrap_field()

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

dtype = None¶

fill_holes = None¶

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_points(**kwargs)[source]¶: Get a 3xN ndarray of N xyz points

get_rotation()[source]¶: Make a rotation matrix to go local -> real 3D space

get_uv_shape(**kwargs)[source]¶

iter_points(**kwargs)[source]¶

Make an iterator that yields (x, y, z) points

This can be overridden in a subclass if it’s more efficient to iterate than a call to _make_points(). Calling iter_points should make an effort to be the fastest way to iterate through the points, regardless of caching.

Note

In Cython, it seems to take twice the time to iterate over a generator than to use a for loop over indices of an array, but only if the array already exists.

local_extent¶

Try to get low and high coordnate values of local points

Raises:	`NotImplementedError` – If subclass can’t make local mesh

local_shape¶: Shape of local representation of seeds

local_to_3d(pts_local)[source]¶

Transform points from the seed coordinate space to 3d

Parameters:	pts_local (ndarray) – An array of points in local crds
Returns:	3xN ndarray of N xyz points

nr_points¶

points(**kwargs)[source]¶: Alias for get_points()

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_extent¶

Try to get low and high coordnate values of uv points

Raises:	`NotImplementedError` – If subclass can’t make uv mesh

uv_shape¶: Shape of uv representation of seeds

uv_to_3d(pts_uv)[source]¶

uv_to_local(pts_uv)[source]¶

wrap_field(data, name='NoName', fldtype='scalar', **kwargs)[source]¶: Wrap an ndarray into a field in the local representation

wrap_mesh(*arrs, **kwargs)[source]¶

class viscid.Point(pts, local_crds=None, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

Collection of points

as_local_coordinates()[source]¶: Make Coordinates for the local representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

class viscid.MeshPoints(pts, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

Generic points with 2d mesh information

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_uv_shape(**kwargs)[source]¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.RectilinearMeshPoints(pts, mesh_axes='xy', cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.MeshPoints

Generic seeds with 2d rect mesh topology

class viscid.Line(p0=(0, 0, 0), p1=(1, 0, 0), n=20, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

A line of seed points

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_uv_shape(**kwargs)[source]¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.Spline(knots, n=-5, splprep_kw=None, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

A spline of seed points

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_uv_shape(**kwargs)[source]¶

to_3d(pts_local, **kwargs)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.Plane(p0=(0, 0, 0), pN=(0, 0, 1), pL=(1, 0, 0), len_l=2, len_m=2, nl=20, nm=20, NL_are_vectors=True, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

A plane of seed points

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape()[source]¶

get_nr_points(**kwargs)[source]¶

get_rotation()[source]¶

Get rotation from lmn -> xyz

To transform xyz -> lmn, transpose the result since it’s an orthogonal matrix.

Note

Remember that the p0 offset (the center of the plane) is not done with the rotation.

Returns:	3x3 ndarray

Example

>>> # Transform a set of points from lmn coordinates to xyz
>>> import numpy as np
>>> import viscid
>>>
>>> # p0 is the origin of the plane
>>> p0 = np.array([0, 0, 0]).reshape(-1, 1)
>>> plane = viscid.Plane(p0=p0, pN=(1, 1, 0),
>>>                      pL=(1, -1, 0),
>>>                      len_l=2.0, len_m=2.0)
>>> lmn_to_xyz = plane.get_rotation()
>>>
>>> # make 10 random points in lmn coorinates
>>> lmn = 2 * np.random.rand(3, 10) - 1
>>> lmn[2] = 0.0
>>> xyz = p0 + np.dot(lmn_to_xyz, lmn)
>>>
>>> from viscid.plot import vlab
>>> vlab.mesh_from_seeds(plane)
>>> vlab.points3d(xyz[0], xyz[1], xyz[2])
>>> vlab.show()

>>> # Transform vector compenents from xyz to lmn
>>> import numpy as np
>>> import viscid
>>>
>>> x = np.linspace(-1, 1, 32)
>>> y = np.linspace(-1, 1, 32)
>>> z = np.linspace(-1, 1, 32)
>>> Bxyz = viscid.zeros([x, y, z], center='node',
>>>                     nr_comps=3, layout="interlaced")
>>> Bxyz['x'] = 1.0
>>>
>>> plane = viscid.Plane(p0=[0, 0, 0], pN=[1, 1, 0],
>>>                      pL=[1, -1, 1],
>>>                      len_l=0.5, len_m=0.5)
>>> vals = viscid.interp_trilin(Bxyz, plane)
>>> B_interp = plane.wrap_field(vals, fldtype='vector',
>>>                             layout='interlaced')
>>> xyz_to_lmn = plane.get_rotation().T
>>> Blmn = np.einsum("ij,lmj->lmi", xyz_to_lmn, B_interp)
>>> Blmn = B_interp.wrap(Blmn)
>>>
>>> # use xyz to show in 3d via mayavi
>>> from viscid.plot import vlab
>>> verts, s = plane.wrap_mesh(Blmn['z'].data)
>>> vlab.mesh(verts[0], verts[1], verts[2], scalars=s)
>>> verts, vx, vy, vz = plane.wrap_mesh(B_interp['x'].data,
>>>                                      B_interp['y'].data,
>>>                                      B_interp['z'].data)
>>> vlab.quiver3d(verts[0], verts[1], verts[2], vx, vy, vz)
>>> vlab.show()
>>>
>>> # use lmn to show in-plane / out-of-plane
>>> from viscid.plot import vpyplot as vlt
>>> from matplotlib import pyplot as plt
>>> vlt.plot(Blmn['z'])  # z means n here
>>> vlt.plot2d_quiver(Blmn)
>>> plt.show()

get_uv_shape()[source]¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.Volume(xl=(-1, -1, -1), xh=(1, 1, 1), n=(20, 20, 20), cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

A volume of seed points

Defined by two opposite corners of a box in 3D

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_uv_shape(**kwargs)[source]¶

iter_points(**kwargs)[source]¶

Make an iterator that yields (x, y, z) points

Note

In Cython, it seems to take twice the time to iterate over a generator than to use a for loop over indices of an array, but only if the array already exists.

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.Sphere(p0=(0, 0, 0), r=0.0, pole=(0, 0, 1), ntheta=20, nphi=20, thetalim=(0, 180.0), philim=(0, 360.0), roll=0.0, crd_system=None, theta_endpoint='auto', phi_endpoint='auto', pole_is_vector=True, theta_phi=False, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

Make seeds on the surface of a sphere

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_rotation()[source]¶: Make a rotation matrix to go local -> real 3D space

get_uv_shape(**kwargs)[source]¶

periodic¶

pt_bnds¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.SphericalCap(angle0=0.0, angle=90.0, theta_endpoint=[True, False], **kwargs)[source]¶

Bases: viscid.seed.Sphere

A spherical cone or cap of seeds

Defined by a center, and a point indicating the direction of the cone, and the half angle of the cone.

This is mostly a wrapper for Sphere that sets thetalim, but it also does some special stuff when wrapping meshes to remove the last theta value so that the cap isn’t closed at the bottom.

angle¶

angle0¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

class viscid.Circle(n=20, endpoint=None, **kwargs)[source]¶

Bases: viscid.seed.SphericalCap

A circle of seeds

Defined by a center and a point normal to the plane of the circle

endpoint¶

n¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

class viscid.SphericalPatch(p0=(0, 0, 0), p1=(0, 0, 1), max_alpha=45, max_beta=45, nalpha=20, nbeta=20, roll=0.0, r=0.0, p1_is_vector=True, cache=False, dtype=None, **kwargs)[source]¶

Bases: viscid.seed.SeedGen

Make a rectangular (in theta and phi) patch on a sphere

arr2mesh(arr, fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶

as_local_coordinates()[source]¶: Make Coordinates for the local representation

as_mesh(fill_holes=<class 'viscid.NOT_SPECIFIED'>)[source]¶: Make a 3xUxV array that describes a 3D mesh surface

as_uv_coordinates()[source]¶: Make Coordinates for the uv representation

get_local_shape(**kwargs)[source]¶

get_nr_points(**kwargs)[source]¶

get_rotation()[source]¶: Make a rotation matrix to go local -> real 3D space

get_uv_shape(**kwargs)[source]¶

to_3d(pts_local)[source]¶

to_local(pts_3d)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

uv_to_local(pts_uv)[source]¶

class viscid.PolarIonosphere(*args, **kwargs)[source]¶

Bases: viscid.seed.Sphere

Place holder for future seed to cover N+S poles

to_local(**kwargs)[source]¶

Transform points from 3d to the seed coordinate space

Parameters:	pts_local (ndarray) – A 3xN array of N points in 3d
Returns:	ndarray similar in shape to self.local_shape

viscid.prune_comp_sel(sel_list, comp_names)[source]¶: Discover and remove vector component from selection list

viscid.raw_sel2sel_list(sel)[source]¶

Turn generic selection into something standard we can work with

Parameters:	sel (object) – some slice selection
Returns:	items in sel_list are guarenteed to be of type {slice, int, complex, string, numpy.ndarray}
Return type:	(sel_list, input_type)
Raises:	`TypeError` – slice-by-array invalid dtype `ValueError` – slice-by-array not 1D

viscid.fill_nd_sel_list(sel_list, ax_names)[source]¶: fully determine a sparsely selected sel_list

viscid.standardize_sel_list(sel_list)[source]¶: turn all selection list elements into fundamental python types

viscid.standardize_sel(sel)[source]¶: turn selection list element into fundamental python types

viscid.standardize_value(sel, bool_argwhere=False)[source]¶

Turn a value element to fundamental type or array

Returns:	One of the following:: None np.newaxis Ellipsis bool int complex (slice by value) numpy.datetime64 numpy.timedelta64 ndarray numpy.integer numpy.bool_ numpy.timedelta64 numpy.datetime64 numpy.complex

viscid.std_sel_list2index(std_sel_list, crd_arrs, val_endpoint=True, interior=False, tdunit='s', epoch=None, tol=100)[source]¶: turn standardized selection list into index slice

viscid.std_sel2index(std_sel, crd_arr, val_endpoint=True, interior=False, tdunit='s', epoch=None)[source]¶

Turn single standardized selection into slice (by int or None)

Normally (val_endpoint=True, interior=False), the rules for float lookup val_endpoints are:

- The slice will never include an element whose value in arr
  is < start (or > if the slice is backward)
- The slice will never include an element whose value in arr
  is > stop (or < if the slice is backward)
- !! The slice WILL INCLUDE stop if you don't change
  val_endpoint. This is different from normal slicing, but
  it's more natural when specifying a slice as a float.

If interior=True, then the slice is expanded such that start and stop are interior to the sliced array.

Parameters:

std_sel – single standardized selection
arr (ndarray) – filled with floats to do the lookup
val_endpoint (bool) – iff True then include stop in the slice when slicing-by-value (DOES NOT EFFECT SLICE-BY-INDEX). Set to False to get python slicing symantics when it comes to excluding stop, but fair warning, python symantics feel awkward here. Consider the case [0.1, 0.2, 0.3][:0.25]. If you think this should include 0.2, then leave keep val_endpoint=True.
interior (bool) – if True, then extend both ends of the slice such that slice-by-location endpoints are interior to the slice
epoch (datetime64-like) – Epoch for to go datetime64 <-> float
tdunit (str) – Presumed time unit for floats
tol (int) – number of machine epsilons to consider “close enough”

viscid.sel_list2values(arrs, sel_list, epoch=None, tdunit='s')[source]¶: find the extrema values for a given sel list

viscid.sel2values(arr, sel, epoch=None, tdunit='s')[source]¶

find the extrema values for a given selection

Parameters:

arr (None, sequence) – array that is being selected, if this is None, then the output will contain np.nan where it can not infer values.
selection (slice, str) – a single selection that could be given to to_slice()
epoch (datetime64-like) – Epoch for to go datetime64 <-> float
tdunit (str) – Presumed time unit for floats

Returns:

(start_val, stop_val) as floats

If arr is None and start/stop are None, then they will become +/- inf depending on the sign of step.
If arr is None and start/stop are slice-by-index, they will become NaN.

viscid.selection2values(arr, sel, epoch=None, tdunit='s')¶

find the extrema values for a given selection

Parameters:

arr (None, sequence) – array that is being selected, if this is None, then the output will contain np.nan where it can not infer values.
selection (slice, str) – a single selection that could be given to to_slice()
epoch (datetime64-like) – Epoch for to go datetime64 <-> float
tdunit (str) – Presumed time unit for floats

Returns:

(start_val, stop_val) as floats

If arr is None and start/stop are None, then they will become +/- inf depending on the sign of step.
If arr is None and start/stop are slice-by-index, they will become NaN.

viscid.make_fwd_slice(shape, slices, reverse=None, cull_second=True)[source]¶

Make sure slices go forward

This function returns two slices equivalent to slices such that the first slice always goes forward. This is necessary because h5py can’t deal with reverse slices such as [::-1].

The optional reverse can be used to interpret a dimension as flipped. This is used if the indices in a slice are based on a coordinate array that has already been flipped. For instance, the result is equivalent to arr[::-1][slices], but in a way that can be handled by h5py. This lets us efficiently load small subsets of large arrays on disk, which is most useful when the large array is coming through sshfs.

Note

The only restriction on slices is that neither start nor stop can be outide the range [-L, L].

Parameters:

shape – shape of the array that is to be sliced
slices – a tuple of slices to work with
reverse (optional) – list of bools that indicate if the corresponding value in slices should be ineterpreted as flipped
cull_second (bool, optional) – iff True, remove elements of the second slice for dimensions that don’t exist after the first slice has completed. This is only here for a super-hacky case when slicing fields.

Returns:

(first_slice, second_slice)

first_slice: a forward-only slice that retrieves the desired elements of an array
second_slice: a slice that does [::1] or [::-1] as needed to make the result equivalent to slices. If keep_all, then this may contain None indicating that this dimension no longer exists after the first slice.

Examples

>> a = np.arange(8) >> first, second = make_fwd_slice(len(a),slice(None, None, -1)) >> (a[::-1] == a[first][second]).all() True

>> a = np.arange(4*5*6).reshape((4, 5, 6)) >> first, second = make_fwd_slice(a.shape, >> [slice(None, -1, 1), >> slice(-1, None, 1), >> slice(-4, -1, 2)], >> [True, True, True]) >> a1 = a[::-1, ::-1, ::-1][:-1, -1:, -4:-1:2] >> a2 = a[first][second] >> a1 == a2 True

viscid.all_slices_none(slices)[source]¶: true iff all slices have no effect

class viscid.UnimportedModule(exception, msg='', **attrs)[source]¶

Bases: object

Proxy object for unimported modules

exception viscid.DeferredImportError[source]¶

Bases: exceptions.ImportError

So lack of an optional dependency doesn’t make viscid unimportable

exception viscid.BackendNotFound[source]¶

Bases: exceptions.RuntimeError

Calculator backend not installed

exception viscid.NoBasetimeError(msg)[source]¶

Bases: exceptions.Exception

When a dataset is trying to get time as datetime but has no basetime

exception viscid.KeyboardInterruptError[source]¶

Bases: exceptions.Exception

Deprecated and unused

viscid.timeit(f, *args, **kwargs)[source]¶

overly simple timeit wrapper

Keyword Arguments:
Parameters:	f – callable to timeit args – positional arguments for f *kwargs – keyword arguments for f
	timeit_repeat (int) – number of times to call f (Default: 1) timeit_print_stats (bool) – print min/max/mean/median when done timeit_quet (bool) – quiets all output (useful if you only want the timeit_stats dict filled) timeit_stats (dict) – Stats will be stuffed into here
Returns:	The result of f(args, *kwargs)

viscid.resolve_path(dset, loc, first=False)[source]¶

Search for globbed paths in a nested dict-like hierarchy

Parameters:	dset (dict) – Root of some nested dict-like hierarchy loc (str) – path as a glob pattern first (bool) – Stop at first match and return a single value
Raises:	`KeyError` – If there are no glob matches
Returns:	If first == True, (value, path) else, ([value0, value1, …], [path0, path1, …])

viscid.find_item(dset, loc)[source]¶: Shortcut for first resolve_path(), item only

viscid.find_items(dset, loc)[source]¶: Shortcut for resolve_path(), items only

viscid.get_trilinear_field()[source]¶: get a generic trilinear field

viscid.slice_globbed_filenames(glob_pattern)[source]¶

Apply a slice to a glob pattern

Note

Slice by value works by adding an ‘f’ to a value, as like the rest of Viscid.

Parameters:	glob_pattern (str) – A string
Returns:	list of filenames

Examples

If a directory contains files,

>>> os.listdir()
["file.010.txt", "file.020.txt", "file.030.txt", "file.040.txt"]

then sliced globs can look like

>>> expand_glob_slice("f*.[:2].txt")
["file.010.txt", "file.020.txt"]

>>> expand_glob_slice("f*.[10.0j::2].txt")
["file.010.txt", "file.030.txt"]

>>> expand_glob_slice("f*.[20j:2].txt")
["file.020.txt", "file.040.txt"]

viscid.glob2(glob_pattern, *args, **kwargs)[source]¶

Wrap slice_globbed_filenames, but return [] on no match

See also

slice_globbed_filenames()

viscid.interact(banner=None, ipython=True, stack_depth=0, global_ns=None, local_ns=None, viscid_ns=True, mpl_ns=False, mvi_ns=False)[source]¶: Start an interactive interpreter

viscid.convective_deriv(a, b=None, bnd=True)[source]¶: Compute (a dot nabla) b for vector fields a and b

viscid.project_vector(a, b)[source]¶: Calculates the vector (a dot b_hat) * b_hat

viscid.project_along_line(line, fld, interp_kind='trilin')[source]¶

Project a Vector Field Parallel to a streamline

Parameters:	line (ndarray) – 3xN of points along the fld (VectorField) – Field to interpolate and project onto the line interp_kind (str) – which interpolation to use, const or trilin

viscid.resample_lines(lines, factor=1, kind='linear')[source]¶

Resample a bunch of lines

Parameters:	lines (sequence) – a list of 3xN ndarrays that are lines factor (float) – how to resample. 2 will give a point at the midpoint of all points, 3 will give two points, etc. 0.5 will leave half the number of points, etc. This should be > 0.0. kind (str) – How to interpolate; must be a value recognized by `scipy.interpolate.interp1d()`
Raises:	`ValueError` – List of resampled 3xN ndarrays

viscid.integrate_along_line(line, fld, reduction='dot', mask_func=None, interp_kind='trilin')[source]¶

Integrate the value of fld along a line

Parameters:	line (ndarray) – 3xN ndarray fld (Field) – Field to interpolate / integrate reduction (str) – If fld is a vector field, what quantity to integrate. Can be “dot” to dot the vectors with ds along the line, or “mag” to integrate the magnitude. interp_kind (str) – which interpolation to use, const or trilin
Returns:	a scalar with the same dtype as fld

viscid.integrate_along_lines(lines, fld, reduction='dot', mask_func=None, interp_kind='trilin')[source]¶

Integrate the value of fld along a list of lines

Parameters:	lines (list) – list of 3xN ndarrays, N needs not be the same for all lines fld (Field) – Field to interpolate / integrate reduction (str) – If fld is a vector field, what quantity to integrate. Can be “dot” to dot the vectors with ds along the line, or “mag” to integrate the magnitude. interp_kind (str) – which interpolation to use, const or trilin
Returns:	ndarray with shape (len(lines), )

viscid.jacobian_at_point(B, x, y, z, dx=None, dy=None, dz=None)[source]¶

Get the Jacobian at a point

If dx|dy|dz == None, then their set to the grid spacing at the point of interest

Returns:	The jacobian as a 3x3 ndarray. The result is in xyz order, in other words: [ [d_x Bx, d_y Bx, d_z Bx], [d_x By, d_y By, d_z By], [d_x Bz, d_y Bz, d_z Bz] ]

viscid.jacobian_at_ind(B, ix, iy, iz)[source]¶

Get the Jacobian at index

Returns:	The jacobian as a 3x3 ndarray. The result is in xyz order, in other words: [ [d_x Bx, d_y Bx, d_z Bx], [d_x By, d_y By, d_z By], [d_x Bz, d_y Bz, d_z Bz] ]

viscid.jacobian_eig_at_point(B, x, y, z, dx=None, dy=None, dz=None)[source]¶

Get the eigen vals/vecs of the jacobian

Returns: evals, evecs (3x3 ndarray): The evec[:, i] corresponds to evals[i]. Eigen vectors are returned in xyz order, aka evec[:, 0] is [x, y, z] for the 0th eigen vector

viscid.jacobian_eig_at_ind(B, ix, iy, iz)[source]¶

Get the eigen vals/vecs of the jacobian

Returns: evals, evecs (3x3 ndarray): The evec[:, i] corresponds to evals[i]. Eigen vectors are returned in xyz order, aka evec[:, 0] is [x, y, z] for the 0th eigen vector

viscid.div_at_point(A, x, y, z, dx=None, dy=None, dz=None)[source]¶: Returns divergence at a point

viscid.curl_at_point(A, x, y, z, dx=None, dy=None, dz=None)[source]¶: Returns curl at point as ndarray with shape (3,) xyz

viscid.extend_boundaries(fld, nl=1, nh=1, axes='all', nr_comp=None, order=1, crd_order=1, crds=None, invarient_dx=0.0, crd_invarient_dx=1.0)[source]¶

Extend and pad boundaries of field (leaves new corners @ 0.0)

Parameters:	fld (Field) – Field to extend nl (int) – extend this many cells in lower direction nh (int) – extend this many cells in upper direction axes (str, list) – something like ‘xyz’, [‘x’, ‘theta’], etc. nr_comp (int, None) – index of shape that corresponds to vector component dimension order (int) – extrapolation order for data; 0 for repeating boundary values or 1 for linear extrapolation crd_order (int) – extrapolation order for crds; 0 for repeating boundary values or 1 for linear extrapolation crds (Coordinates) – Use these coordinates, no extrapolate them invarient_dx (float) – if fld has a single cell in a given dimension, and order == 1, the field is extended with constant dx, using 0 here is synonymous with 0th order invarient_dx – if fld has a single cell in a given dimension, and crd_order == 1, the coordinates in that dimension are extended with constant dx, using 0 here is synonymous with crd_order = 0.
Returns:	A new extended / padded ndarray
Return type:	ndarray

Warning

When using crd_order=0 (0 order hold), coordinates will be extended by repeating boundary values, so min_dx will be 0. Keep this in mind if dividing by dx.

viscid.extend_boundaries_ndarr(arr, nl=1, nh=1, axes='all', nr_comp=None, order=1, invarient_dx=0.0)[source]¶

Extend and pad boundaries of ndarray (leaves new corners @ 0.0)

Parameters:	arr (ndarray) – Array to extend nl (int) – extend this many cells in lower direction nh (int) – extend this many cells in upper direction axes (list) – list of ints that correspond to axes that should be extended. If nr_comp != None, then axes > nr_comp will be shifted by 1. nr_comp (int, None) – index of shape that corresponds to vector component dimension order (int) – 0 for repeating boundary values or 1 for linear extrapolation invarient_dx (float) – if len(arr) == 1 and order == 1, then the array is extended with constant dx, using 0 here is synonymous with 0th order.
Returns:	A new extended / padded ndarray
Return type:	ndarray

viscid.distance_to_clusters(point, clusters, alt=())[source]¶: L2 distance between point and clusters

viscid.find_clusters(indx, indy, x, y, multiple=True, periodic=(False, False), diagonals=True)[source]¶

Cluster and average groups of neighboring points

TODO: If absolutely necessary, could do some K-means clustering: here by calling into scikit-learn.

Parameters:	indx (sequence) – list of x indices indy (sequence) – list of y indices x (sequence) – x coordinate array y (sequence) – y coordinate array multiple (bool) – If False, average all points as a single cluster periodic (sequence) – indicate whether that direction is periodic, and if so, whether the coordinate arrays are overlapped or not. Values can be True, False, or ‘+’. ‘+’ indicates that x[0] and x[-1] are not colocated, so assume they’re dx apart where dx = x[-1] - x[-2]. diagonals (bool) – if true, then diagonal points are considered neighbors
Returns:	2xN for N clusters
Return type:	ndarray

viscid.make_ecfc_field_leading(fc, trim_leading=True)[source]¶: Standardize staggering on edge / face centered fields

viscid.fc2cc(fc, force_numpy=False, bnd=True)[source]¶: Average a face centered field to cell centers

viscid.ec2cc(ec, force_numpy=False, bnd=True)[source]¶: Average an edge centered field to cell centers

viscid.div_fc(fc, force_numpy=False, bnd=True)[source]¶: Calculate cell centered divergence of face centered field

viscid.evaluate(grid, result_name, eqn, try_numexpr=True, slc=Ellipsis)[source]¶

Evaluate an equation on a grid

Examples

To use this function directly

>>> evaluator.enabled = True
>>> f = viscid.load_file("...")
>>> evaluator.evaluate(f.get_grid(),
                       "sqrt(vx**2+vy**2+vz**2)",
                       "speed")
<viscid.field.ScalarField object at ...>

Or, for short, you can as a grid to evaluate implicitly,

>> evaluator.enabled = True >> f = viscid.load_file(“…”) >> speed = f[“speed=sqrt(vx**2+vy**2+vz**2)”] <viscid.field.ScalarField object at …>

Parameters:	grid – a grid instance where the fields live result_name (str) – Used for the name and pretty_name of the resulting field eqn (str) – the equation, if a symbol exists in the numpy namespace, then that’s how it is interpreted, otherwise, the symbol will be looked up in the grid
Returns:	Field instance

viscid.minvar(B, p1, p2, n=40)[source]¶

Find minimum variance eigenvectors of B between p1 and p2

Do minimum variance analysis (MVA) of vector field B between two points, i.e., p1 and p2 are like two spacecraft locations on either side of the boundary for the MVA.

Parameters:	B (`viscid.field.VectorField`) – Vector field for MVA p1 (sequence) – xyz point on one side of the boundary p2 (sequence) – xyz point on the other side of the boundary n (int) – number of points to sample B for the MVA
Returns:	(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]

viscid.minvar_series(arr, warn_finite=True)[source]¶

Find minimum variance eigenvectors of vector series

Parameters:	arr (ndarray) – vector time series with shape (nsamples, 3)
Returns:	(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]

viscid.minvar_around(B, p0, l=1.0, path_dir=(1, 0, 0), n=40)[source]¶

Find minimum variance eigenvectors of B around point p0

Do minimum variance analysis (MVA) of vector field B a distance l/2 around point p0 in the direction of path_dir, i.e., path_dir is like the spacecraft travel direction.

Parameters:	B (`viscid.field.VectorField`) – Vector field for MVA p0 (sequence) – xyz center point, within -l/2 and +l/2 of the boundary in the path_dir directon l (float) – distance on either side of center point path_dir (tuple) – spacecraft travel direction n (int) – number of points to sample B for the MVA
Returns:	(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]

See also

minvar()

viscid.find_minvar_lmn(B, p1, p2, n=40, l_basis=(0, 0, 1))[source]¶

Find rotation matrix for going to boundary normal (lmn) crds

Parameters:	B (`viscid.field.VectorField`) – Vector field for MVA p1 (sequence) – xyz point on one side of the boundary p2 (sequence) – xyz point on the other side of the boundary n (int) – number of points to sample B for the MVA l_basis (sequence) – vector used to define the l-direction, i.e., the l direction will be l_basis projected perpendicular to the boundary normal directon. If None, then l_basis will be set to evec_max.
Returns:	3x3 array of l,m,n basis column (2nd axis) vectors
Return type:	ndarray

viscid.find_minvar_lmn_around(B, p0, l=1.0, path_dir=(1, 0, 0), n=40, l_basis=(0, 0, 1))[source]¶

Find rotation matrix for going to boundary normal (lmn) crds

Finds boundary normal from minimum variance analysis (MVA) of vector field B a distance l/2 around point p0 in the direction of path_dir, i.e., path_dir is like the spacecraft travel direction.

Parameters:	B (`viscid.field.VectorField`) – Vector field for MVA p0 (sequence) – xyz center point, within -l/2 and +l/2 of the boundary in the path_dir directon l (float) – distance on either side of center point path_dir (tuple) – spacecraft travel direction n (int) – number of points to sample B for the MVA l_basis (sequence) – vector used to define the l-direction, i.e., the l direction will be l_basis projected perpendicular to the boundary normal directon. If None, then l_basis will be set to evec_max.
Returns:	3x3 array of l,m,n basis column (2nd axis) vectors
Return type:	ndarray

See also

find_nlm()

viscid.paraboloid(y, z, x0, y0, z0, ax, ay, az)[source]¶: Generic paraboloid function

viscid.paraboloid_normal(y, z, x0, y0, z0, ax, ay, az, normalize=True)[source]¶: Normal vector of a generic paraboloid

viscid.fit_paraboloid(fld, p0=(9.0, 0.0, 0.0, 1.0, -1.0, -1.0), tolerance=0.0)[source]¶

Fit paraboloid it GSE coordinates x ~ y**2 + z**2

Parameters:

fld (viscid.field.ScalarField) – field of x values
p0 (sequence) – initial guess for parabaloid (x0, y0, z0, ax, ay, az), where (x0, y0, z0) is the nose location (should be subsolar for 0 dipole tilt), and the ay, az, and az coefficients determine the curvature

Returns:

record array of parameters with length 2; the: 1st value is the fit value, and the 2nd is one sigma of the fit

Return type:

numpy.recarray

viscid.get_mp_info(pp, b, j, e, cache=True, cache_dir=None, slc='x=5.5j:11.0j, y=-4.0j:4.0j, z=-3.6j:3.6j', fit='mp_xloc', fit_p0=(9.0, 0.0, 0.0, 1.0, -1.0, -1.0))[source]¶

Get info about m-pause as flattened fields

Notes

The first thing this function does is mask locations where the GSE-y current density < 1e-4. This masks out the bow shock and current free regions. This works for southward IMF, but it is not very general.

Parameters:	pp (ScalarcField) – pressure b (VectorField) – magnetic field j (VectorField) – current density e (VectorField, None) – electric field (same centering as b). If None, then the info that requires E will be filled with NaN cache (bool, str) – Save to and load from cache, if “force”, then don’t load from cache if it exists, but do save a cache at the end cache_dir (str) – Directory for cache, if None, same directory as that file to which the grid belongs slc (str) – slice that gives a box that contains the m-pause fit (str) – to which resulting field should the paraboloid be fit, defaults to mp_xloc, but pp_max_xloc might be useful in some circumstances fit_p0 (tuple) – Initial guess vector for paraboloid fit
Returns:	Unless otherwise noted, the entiries are 2D (y-z) fields mp_xloc location of minimum abs(Bz), this works better than max of J^2 for FTEs mp_sheath_edge location where Jy > 0.1 * Jy when coming in from the sheath side mp_sphere_edge location where Jy > 0.1 * Jy when coming in from the sphere side mp_width difference between m-sheath edge and msphere edge mp_shear magnetic shear taken 6 grid points into the m-sheath / m-sphere pp_max max pp pp_max_xloc location of max pp epar_max max e parallel epar_max_xloc location of max e parallel paraboloid numpy.recarray of paraboloid fit. The parameters are given in the 0th element, and the 1st element contains the 1-sigma values for the fit
Return type:	dict
Raises:	`RuntimeError` – if using MHD crds instead of GSE crds

viscid.find_mp_edges(j_block, msphere_thresh=0.1, sheath_thresh=0.1, maskval=0.0001)[source]¶

Find x location of msphere and msheath edges using current (J)

Note

GSE coordinates only please

Parameters:

j_block (VectorField) – Current density containing the whole magnetopause
msphere_thresh (float) – thereshold of current on the magnetosphere side as a fraction of the maximum current density, i.e., 0.1 is 10% of the max
sheath_thresh (float) – thereshold of current on the magnetosheath side as a fraction of the maximum current density, i.e., 0.1 is 10% of the max
maskval (float, None) – if not None, then mask out J values less than maskval; useful for masking out bowshock, and current free regions

Returns:

sheath and sphere fields / values

sheath_edge: float or 2D ScalarField of x values

msphere_edge: float or 2D ScalarField of x values

mp_width: sheath_edge - msphere_edge

sheath_ind: index of sheath_edge x location

sphere_ind: index of msphere_edge x location

Return type:

tuple

viscid.calc_psi(B, rev=False)[source]¶

Calc Flux function (only valid in 2d)

Parameters:	B (VectorField) – magnetic field, should only have two spatial dimensions so we can infer the symmetry dimension rev (bool) – since this integration doesn’t like going through undefined regions (like within 1 earth radius of the origin for openggcm), you can use this to start integrating from the opposite corner.
Returns:	2-D scalar flux function
Return type:	ScalarField
Raises:	`ValueError` – If B has <> 2 spatial dimensions

viscid.calc_beta(pp, B, scale=1.0)[source]¶

Calc plasma beta (2.0 * p / B^2)

Parameters:	pp (ScalarField or ndarray) – pressure B (VectorField) – magnetic field scale (float, optional) – overall scale factor
Returns:	Plasma beta
Return type:	ScalarField

Note

For OpenGGCM, where pp is in pPa and B is in nT, scale should be 40.0.

viscid.trace_separator(grid, b_slcstr='x=-25j:15j, y=-30j:30j, z=-15j:15j', r=1.0, plot=False, trace_opts=None, cache=True, cache_dir=None)[source]¶

Trace a separator line from most dawnward null

Still in testing Uses the bisection algorithm.

Parameters:	grid (Grid) – A grid that has a “b” field b_slcstr (str) – Some valid slice for B field r (float) – spatial step of separator line plot (bool) – make debugging plots trace_opts (dict) – passed to streamline function cache (bool, str) – Save to and load from cache, if “force”, then don’t load from cache if it exists, but do save a cache at the end cache_dir (str) – Directory for cache, if None, same directory as that file to which the grid belongs
Raises:	`IOError` – Description
Returns:	(separator_lines, nulls) separator_lines (list): list of M 3xN ndarrays that represent M separator lines with N points nulls (ndarray): 3xN array of N null points
Return type:	tuple

viscid.topology_bitor_clusters(fld, min_depth=1, max_depth=10, multiple=True, plot=False, sep_val=15, mask_limit=15, periodic='00', pt_bnds=())[source]¶

Find separator as intersection of all global topologies

Neighbors are bitwise ORed until at least one value matches sep_val which is presumably (Close | Open N | Open S | SW). This happens between min_depth and max_depth times, where the resolution of each iteration is reduced by a factor of two, ie, worst case 2**(max_depth).

Parameters:	fld (Field) – Topology (bitmask) as a field min_depth (int) – Iterate at least this many times max_depth (int) – Iterate at most this many times multiple (bool) – passed to `viscid.cluster()` sep_val (int) – Value of bitmask that indicates a separator plot (bool) – Make a 2D plot of Fld and the sep candidates mask_limit (int) – if > 0, then bitmask fld with mask_limit, i.e., fld = fld & mask_limit (bitwise and) periodic (sequence) – indicate whether that direction is periodic, and if so, whether the coordinate arrays are overlapped or not. Values can be True, False, or ‘+’. ‘+’ indicates that x[0] and x[-1] are not colocated, so assume they’re dx apart where dx = x[-1] - x[-2]. pt_bnd (sequence) – Boundaries that come to a point, i.e., all values along that boundary are neighbors such as the poles of a sphere. Specified like “0-” for lower boundary of dimension 0 or “1+” for the upper boundary of dimension 1.
Returns:	2xN for N clusters of separator points in the same coordinates as fld
Return type:	ndarray

viscid.get_sep_pts_bitor(fld, seed, trace_opts=None, make_3d=True, **kwargs)[source]¶

bitor topologies to find separator points in uv map from seed

Parameters:	fld (VectorField) – Magnetic Field seed (viscid.seed.SeedGen) – Any Seed generator with a 2d local representation trace_opts (dict) – kwargs for calc_streamlines make_3d (bool) – convert result from uv to 3d space **kwargs – passed to `topology_bitor_clusters()`
Returns:	3xN ndarray of N separator points in uv space or 3d space depending on the make_3d kwarg

viscid.get_sep_pts_bisect(fld, seed, trace_opts=None, min_depth=3, max_depth=7, plot=False, perimeter_check=<function perimeter_check_bitwise_or>, make_3d=True)[source]¶

bisect uv map of seed to find separator points

Parameters:

fld (VectorField) – Magnetic Field
seed (viscid.seed.SeedGen) – Any Seed generator with a 2d local representation
trace_opts (dict) – kwargs for calc_streamlines
min_depth (int) – Min allowable bisection depth
max_depth (int) – Max bisection depth
plot (bool) – Useful for debugging the algorithm
perimeter_check (func) – Some func that returns a bool with the same signature as perimeter_check_bitwise_or()
make_3d (bool) – convert result from uv to 3d space

Returns:

3xN ndarray of N separator points in uv space or 3d space depending on the make_3d kwarg

viscid.topology2color(topology, topo_style='msphere', bad_color=None)[source]¶

Determine RGB from topology value

Parameters:	topology (int, list, ndarray) – some value in `calculator.streamline.TOPOLOGY_` topo_style* (string) – msphere, or a dict with its own mapping bad_color (tuple) – rgb color for invalid topologies
Returns:	Nx3 array of rgb data or (R, G, B) tuple if topology is a single value

class viscid.OrderedDict(**kwds)[source]¶

Bases: dict

Dictionary that remembers insertion order

clear() → None. Remove all items from od.[source]¶

copy() → a shallow copy of od[source]¶

classmethod fromkeys(S[, v]) → New ordered dictionary with keys from S.[source]¶: If not specified, the value defaults to None.

items() → list of (key, value) pairs in od[source]¶

iteritems()[source]¶: od.iteritems -> an iterator over the (key, value) pairs in od

iterkeys() → an iterator over the keys in od[source]¶

itervalues()[source]¶: od.itervalues -> an iterator over the values in od

keys() → list of keys in od[source]¶

pop(k[, d]) → v, remove specified key and return the corresponding[source]¶: value. If key is not found, d is returned if given, otherwise KeyError is raised.

popitem() → (k, v), return and remove a (key, value) pair.[source]¶: Pairs are returned in LIFO order if last is true or FIFO order if false.

setdefault(k[, d]) → od.get(k,d), also set od[k]=d if k not in od[source]¶

update([E, ]**F) → None. Update D from mapping/iterable E and F.¶: If E present and has a .keys() method, does: for k in E: D[k] = E[k] If E present and lacks .keys() method, does: for (k, v) in E: D[k] = v In either case, this is followed by: for k, v in F.items(): D[k] = v

values() → list of values in od[source]¶

viewitems() → a set-like object providing a view on od's items[source]¶

viewkeys() → a set-like object providing a view on od's keys[source]¶

viewvalues() → an object providing a view on od's values[source]¶

class viscid.izip¶

Bases: object

izip(iter1 [,iter2 […]]) –> izip object

Return a izip object whose .next() method returns a tuple where the i-th element comes from the i-th iterable argument. The .next() method continues until the shortest iterable in the argument sequence is exhausted and then it raises StopIteration. Works like the zip() function but consumes less memory by returning an iterator instead of a list.

next¶

class viscid.izip_longest¶

Bases: object

izip_longest(iter1 [,iter2 […]], [fillvalue=None]) –> izip_longest object

Return an izip_longest object whose .next() method returns a tuple where the i-th element comes from the i-th iterable argument. The .next() method continues until the longest iterable in the argument sequence is exhausted and then it raises StopIteration. When the shorter iterables are exhausted, the fillvalue is substituted in their place. The fillvalue defaults to None or can be specified by a keyword argument.

next¶

class viscid.unicode(object='') → unicode object¶

Bases: basestring

unicode(string[, encoding[, errors]]) -> unicode object

Create a new Unicode object from the given encoded string. encoding defaults to the current default string encoding. errors can be ‘strict’, ‘replace’ or ‘ignore’ and defaults to ‘strict’.

capitalize() → unicode¶: Return a capitalized version of S, i.e. make the first character have upper case and the rest lower case.

center(width[, fillchar]) → unicode¶: Return S centered in a Unicode string of length width. Padding is done using the specified fill character (default is a space)

count(sub[, start[, end]]) → int¶: Return the number of non-overlapping occurrences of substring sub in Unicode string S[start:end]. Optional arguments start and end are interpreted as in slice notation.

decode([encoding[, errors]]) → string or unicode¶: Decodes S using the codec registered for encoding. encoding defaults to the default encoding. errors may be given to set a different error handling scheme. Default is ‘strict’ meaning that encoding errors raise a UnicodeDecodeError. Other possible values are ‘ignore’ and ‘replace’ as well as any other name registered with codecs.register_error that is able to handle UnicodeDecodeErrors.

encode([encoding[, errors]]) → string or unicode¶: Encodes S using the codec registered for encoding. encoding defaults to the default encoding. errors may be given to set a different error handling scheme. Default is ‘strict’ meaning that encoding errors raise a UnicodeEncodeError. Other possible values are ‘ignore’, ‘replace’ and ‘xmlcharrefreplace’ as well as any other name registered with codecs.register_error that can handle UnicodeEncodeErrors.

endswith(suffix[, start[, end]]) → bool¶: Return True if S ends with the specified suffix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. suffix can also be a tuple of strings to try.

expandtabs([tabsize]) → unicode¶: Return a copy of S where all tab characters are expanded using spaces. If tabsize is not given, a tab size of 8 characters is assumed.

find(sub[, start[, end]]) → int¶

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

format(*args, **kwargs) → unicode¶: Return a formatted version of S, using substitutions from args and kwargs. The substitutions are identified by braces (‘{‘ and ‘}’).

index(sub[, start[, end]]) → int¶: Like S.find() but raise ValueError when the substring is not found.

isalnum() → bool¶: Return True if all characters in S are alphanumeric and there is at least one character in S, False otherwise.

isalpha() → bool¶: Return True if all characters in S are alphabetic and there is at least one character in S, False otherwise.

isdecimal() → bool¶: Return True if there are only decimal characters in S, False otherwise.

isdigit() → bool¶: Return True if all characters in S are digits and there is at least one character in S, False otherwise.

islower() → bool¶: Return True if all cased characters in S are lowercase and there is at least one cased character in S, False otherwise.

isnumeric() → bool¶: Return True if there are only numeric characters in S, False otherwise.

isspace() → bool¶: Return True if all characters in S are whitespace and there is at least one character in S, False otherwise.

istitle() → bool¶: Return True if S is a titlecased string and there is at least one character in S, i.e. upper- and titlecase characters may only follow uncased characters and lowercase characters only cased ones. Return False otherwise.

isupper() → bool¶: Return True if all cased characters in S are uppercase and there is at least one cased character in S, False otherwise.

join(iterable) → unicode¶: Return a string which is the concatenation of the strings in the iterable. The separator between elements is S.

ljust(width[, fillchar]) → int¶: Return S left-justified in a Unicode string of length width. Padding is done using the specified fill character (default is a space).

lower() → unicode¶: Return a copy of the string S converted to lowercase.

lstrip([chars]) → unicode¶: Return a copy of the string S with leading whitespace removed. If chars is given and not None, remove characters in chars instead. If chars is a str, it will be converted to unicode before stripping

partition(sep) -> (head, sep, tail)¶: Search for the separator sep in S, and return the part before it, the separator itself, and the part after it. If the separator is not found, return S and two empty strings.

replace(old, new[, count]) → unicode¶: Return a copy of S with all occurrences of substring old replaced by new. If the optional argument count is given, only the first count occurrences are replaced.

rfind(sub[, start[, end]]) → int¶

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

rindex(sub[, start[, end]]) → int¶: Like S.rfind() but raise ValueError when the substring is not found.

rjust(width[, fillchar]) → unicode¶: Return S right-justified in a Unicode string of length width. Padding is done using the specified fill character (default is a space).

rpartition(sep) -> (head, sep, tail)¶: Search for the separator sep in S, starting at the end of S, and return the part before it, the separator itself, and the part after it. If the separator is not found, return two empty strings and S.

rsplit([sep[, maxsplit]]) → list of strings¶: Return a list of the words in S, using sep as the delimiter string, starting at the end of the string and working to the front. If maxsplit is given, at most maxsplit splits are done. If sep is not specified, any whitespace string is a separator.

rstrip([chars]) → unicode¶: Return a copy of the string S with trailing whitespace removed. If chars is given and not None, remove characters in chars instead. If chars is a str, it will be converted to unicode before stripping

split([sep[, maxsplit]]) → list of strings¶: Return a list of the words in S, using sep as the delimiter string. If maxsplit is given, at most maxsplit splits are done. If sep is not specified or is None, any whitespace string is a separator and empty strings are removed from the result.

splitlines(keepends=False) → list of strings¶: Return a list of the lines in S, breaking at line boundaries. Line breaks are not included in the resulting list unless keepends is given and true.

startswith(prefix[, start[, end]]) → bool¶: Return True if S starts with the specified prefix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. prefix can also be a tuple of strings to try.

strip([chars]) → unicode¶: Return a copy of the string S with leading and trailing whitespace removed. If chars is given and not None, remove characters in chars instead. If chars is a str, it will be converted to unicode before stripping

swapcase() → unicode¶: Return a copy of S with uppercase characters converted to lowercase and vice versa.

title() → unicode¶: Return a titlecased version of S, i.e. words start with title case characters, all remaining cased characters have lower case.

translate(table) → unicode¶: Return a copy of the string S, where all characters have been mapped through the given translation table, which must be a mapping of Unicode ordinals to Unicode ordinals, Unicode strings or None. Unmapped characters are left untouched. Characters mapped to None are deleted.

upper() → unicode¶: Return a copy of S converted to uppercase.

zfill(width) → unicode¶: Pad a numeric string S with zeros on the left, to fill a field of the specified width. The string S is never truncated.

viscid.interp()¶

Interpolate a field to points described by seeds

Note

Nearest neighbor is always used between the last value and vfield.crds.xh. This is done to keep from extrapolating and introducing new maxima. As such, AMR grids may be more step-like at patch boundaries.

Parameters:	vfield (viscid.field.Field) – Some Vector or Scalar field. If this field is not 3D, then vfield.atleast_3d() is called seeds (viscid.claculator.seed) – locations for the interpolation kind (str) – either ‘linear’ or ‘nearest’ wrap (bool) – if true, then call seeds.wrap on the result method (str) – alias for kind, because why not
Returns:	numpy.ndarray of interpolated values. Shaped (seed.nr_points,) or (seed.nr_points, vfield.nr_comps) if vfield is a Scalar or Vector field.

viscid.interp_nearest()¶

Interpolate a field to points described by seeds

Parameters:	vfield (viscid.field.Field) – Some Vector or Scalar field seeds (viscid.claculator.seed) – locations for the interpolation wrap (bool) – if true, call seeds.wrap on the result
Returns:	numpy.ndarray of interpolated values. Shaped (seed.nr_points,) or (seed.nr_points, vfield.nr_comps) if vfield is a Scalar or Vector field.

viscid.interp_linear()¶

Interpolate a field to points described by seeds

Note

Nearest neighbor is used between the last value and vfield.crds.xh. This is done to keep from extrapolating and introducing new maxima.

Parameters:	vfield (viscid.field.Field) – Some Vector or Scalar field seeds (viscid.claculator.seed) – locations for the interpolation wrap (bool) – if true, then call seeds.wrap on the result
Returns:	numpy.ndarray of interpolated values. Shaped (seed.nr_points,) or (seed.nr_points, vfield.nr_comps) if vfield is a Scalar or Vector field.

viscid.interp_trilin()¶

Interpolate a field to points described by seeds

Note

Nearest neighbor is used between the last value and vfield.crds.xh. This is done to keep from extrapolating and introducing new maxima.

Parameters:	vfield (viscid.field.Field) – Some Vector or Scalar field seeds (viscid.claculator.seed) – locations for the interpolation wrap (bool) – if true, then call seeds.wrap on the result
Returns:	numpy.ndarray of interpolated values. Shaped (seed.nr_points,) or (seed.nr_points, vfield.nr_comps) if vfield is a Scalar or Vector field.

viscid.calc_streamlines()¶

Trace streamlines

Warning

For streamlines that reach max_length or max_t, the line segments at the ends will be trimmed. This may be confusing when using a non-adaptive integrator (ds will be different for the 1st and last segments). Bear this in mind when doing math on the result.

Parameters:

vfield – A VectorField with 3 components, If this field is not 3D, then vfield.atleast_3d() is called
seed – can be a Seeds instance or a Coordinates instance, or anything that exposes an iter_points method
nr_procs – how many processes for streamlines (>1 only works on *nix systems)
force_subprocess (bool) – always calc streamlines in a separate process, even if nr_procs == 1
chunk_factor (int) – Valid range is [1, nr_seeds // nr_procs]. 1 indicates a single chunk per process; use this for perfectly balanced streamlines. nr_seeds // nr_procs indicates one chunk per seed; this solves load balancing, but the overhead of assembling results makes this undesirable. Start with 1 and bump this up if load balancing is still an issue.
wrap (bool) – if true, then call seed.wrap_field on topology
**kwargs – more arguments for streamlines

Keyword Arguments:

ds0 (float) – initial spatial step for streamlines (if 0.0, it will be ds0_frac * the minimum d[xyz])
ds0_frac (float) – If an absolute spatial step ds0 is not given, then it will be set to ds0_frac * min_dx where min_dx is the smallest dimenstion of the smallest grid cell of the vfield. Defaults to 0.5.
ibound (float) – Inner boundary as distance from (0, 0, 0)
obound0 (array-like) – lower corner of outer boundary (x, y, z)
obound1 (array-like) – upper corner of outer boundary (x, y, z)
obound_r (float) – Outer boundary as distance from (0, 0, 0)
maxit (int) – maximum number of line segments
max_length (float) – maximum streamline length (int ds)
max_t (float) – max value for t (int ds / v). I.e., stop a streamline after a fluid element travels dt along a streamline. This has no obvious use for mag field lines.
stream_dir (int) – one of DIR_FORWARD, DIR_BACKWARD, DIR_BOTH
output (int) – which output to provide, one of OUTPUT_STREAMLINE, OUTPUT_TOPOLOGY, or OUTPUT_BOTH
method (int) – integrator, one of EULER1, RK2, RK4, EULER1a (adaptive), RK12 (adaptive, midpoint), RK45 (adaptive, Fehlberg). Note that sometimes RK45 is faster than lower order adaptive methods since it can take much larger step sizes
max_error (float) – max allowed error between methods for adaptive integrators. This should be as a fraction of ds (i.e., between 0 and 1). The default value is 4e-2 for euler1a, 1.5e-2 for rk12, and 1e-5 for rk45.
smallest_ds (float) – smallest absolute spatial step. If not set then smallest_ds = smallest_ds_frac * ds0
largest_ds (float) – largest absolute spatial step. If not set then largest_ds = largest_ds_frac * ds0
smallest_ds_frac (float) – smallest spatial step as fraction of ds0
largest_ds_frac (float) – largest spatial step as fraction of ds0
topo_style (str) – how to map end point bitmask to a topology. ‘msphere’ means map to TOPOLOGY_MS_* and ‘generic’ means leave topology as a bitmask of END_*

Returns:

(lines, topo), either can be None depending on output

lines: list of nr_streams ndarrays, each ndarray has shape (3, nr_points_in_stream). The nr_points_in_stream can be different for each line
topo: ndarray with shape (nr_streams,) of topology bitmask with values depending on the topo_style

viscid.streamlines()¶

Trace streamlines

Warning

Parameters:

vfield – A VectorField with 3 components, If this field is not 3D, then vfield.atleast_3d() is called
seed – can be a Seeds instance or a Coordinates instance, or anything that exposes an iter_points method
nr_procs – how many processes for streamlines (>1 only works on *nix systems)
force_subprocess (bool) – always calc streamlines in a separate process, even if nr_procs == 1
chunk_factor (int) – Valid range is [1, nr_seeds // nr_procs]. 1 indicates a single chunk per process; use this for perfectly balanced streamlines. nr_seeds // nr_procs indicates one chunk per seed; this solves load balancing, but the overhead of assembling results makes this undesirable. Start with 1 and bump this up if load balancing is still an issue.
wrap (bool) – if true, then call seed.wrap_field on topology
**kwargs – more arguments for streamlines

Keyword Arguments:

ds0 (float) – initial spatial step for streamlines (if 0.0, it will be ds0_frac * the minimum d[xyz])
ds0_frac (float) – If an absolute spatial step ds0 is not given, then it will be set to ds0_frac * min_dx where min_dx is the smallest dimenstion of the smallest grid cell of the vfield. Defaults to 0.5.
ibound (float) – Inner boundary as distance from (0, 0, 0)
obound0 (array-like) – lower corner of outer boundary (x, y, z)
obound1 (array-like) – upper corner of outer boundary (x, y, z)
obound_r (float) – Outer boundary as distance from (0, 0, 0)
maxit (int) – maximum number of line segments
max_length (float) – maximum streamline length (int ds)
max_t (float) – max value for t (int ds / v). I.e., stop a streamline after a fluid element travels dt along a streamline. This has no obvious use for mag field lines.
stream_dir (int) – one of DIR_FORWARD, DIR_BACKWARD, DIR_BOTH
output (int) – which output to provide, one of OUTPUT_STREAMLINE, OUTPUT_TOPOLOGY, or OUTPUT_BOTH
method (int) – integrator, one of EULER1, RK2, RK4, EULER1a (adaptive), RK12 (adaptive, midpoint), RK45 (adaptive, Fehlberg). Note that sometimes RK45 is faster than lower order adaptive methods since it can take much larger step sizes
max_error (float) – max allowed error between methods for adaptive integrators. This should be as a fraction of ds (i.e., between 0 and 1). The default value is 4e-2 for euler1a, 1.5e-2 for rk12, and 1e-5 for rk45.
smallest_ds (float) – smallest absolute spatial step. If not set then smallest_ds = smallest_ds_frac * ds0
largest_ds (float) – largest absolute spatial step. If not set then largest_ds = largest_ds_frac * ds0
smallest_ds_frac (float) – smallest spatial step as fraction of ds0
largest_ds_frac (float) – largest spatial step as fraction of ds0
topo_style (str) – how to map end point bitmask to a topology. ‘msphere’ means map to TOPOLOGY_MS_* and ‘generic’ means leave topology as a bitmask of END_*

Returns:

(lines, topo), either can be None depending on output

lines: list of nr_streams ndarrays, each ndarray has shape (3, nr_points_in_stream). The nr_points_in_stream can be different for each line
topo: ndarray with shape (nr_streams,) of topology bitmask with values depending on the topo_style

viscid.find_classified_nulls()¶

Find nulls, and classify them as in Cowley 1973

Parameters:

fld (VectorField) – Magnetic field with nulls
ibound (float) – ignore points within ibound of origin
rtol (float) – used in np.isclose when classifying nulls
atol (float) – used in np.isclose when classifying nulls

Returns:

{“O”: [Oinds, Opts], “A”: […], “B”: […]}: each of these are 3xN ndarrays of either ingegers of floats

Return type:

dict

viscid.find_nulls()¶

Just find null points and closest indices in fld

Parameters:

fld (VectorField) – Magnetic field with nulls
ibound (float) – ignore points within ibound of origin

Returns:

(null_inds, null_pts) both of which are 3xN ndarrays,: the first is integers, the second floats

Return type:

tuple

viscid.load_file(fname, force_reload=False, **kwargs)[source]¶

Load a file

Parameters:	fnames (list) – single file name, or list of files that are part of the same time series. Glob patterns and slices are accepted, see Tips & Tricks for more info. fname (str) – a file name, relative to CWD force_reload (bool) – Force reload if file is already in memory **kwargs – passed to the VFile constructor

See also

Tips & Tricks

Returns:	A VFile instance

viscid.load_files(fnames, force_reload=False, **kwargs)[source]¶

Load a list of files

Parameters:	fnames (list) – list of file names. Glob patterns and slices are accepted, see Tips & Tricks for more info. force_reload (bool) – Force reload if file is already in memory **kwargs – passed to the VFile constructor

See also

Tips & Tricks

Returns:	A list of VFile instances. The length may not be the same as the length of fnames, and the order may not be the same in order to accommodate globs and file grouping.

viscid.unload_file(handle)[source]¶: call unload on the handle in the bucket

viscid.unload_all_files()[source]¶: Hammer-of-Thor the cache

viscid.reload_file(handle)[source]¶: call reload on the handle in the bucket

viscid.get_file(handle)[source]¶: return a file that’s already been loaded by either number (as in nth file loaded), of file name

viscid.save_grid(fname, grd, **kwargs)[source]¶: save a grid, filetype is inferred from fname

viscid.save_field(fname, fld, **kwargs)[source]¶: save a field, filetype is inferred from fname

viscid.save_fields(fname, flds, **kwargs)[source]¶: save a list of fields, filetype is inferred from fname

viscid.check_version()[source]¶: Check status of viscid and associated libraries and modules

viscid.check()[source]¶: Runtime check compiled modules

Parameters:	p1 (sequence) – first point as [x, y, z] p2 (sequence) – second point as [x, y, z] origin (sequence) – origin of the sphere n (int) – Number of line segments along the great circle
Returns:	3xN ndarray

Parameters:	time – some python datetime or string in ISO 8601 format, could also be a sequence of these to return a Numpy ndarray unit (str) – one of {Y,M,W,D,h,m,s,m,s,us,ns,ps,fs,as}
Returns:	np.datetime64[unit] or array with dtype np.datetime64[unit]

Parameters:	tlst (timelike, list) – single or list of datetime64 or timedelta64 unit (str) – If given, regularize all times to this unit, otherwise, regularize them to the most precise of the bunch most_precise (bool) – If True, then convert all times to the most precise unit that fits all the times allow_rounding (bool) – if any time is too small to be represented in the desired unit, then use a more coarse unit and round values so that everything fits allow0 (bool) – If False, then raise PrecisionError if a value has been rounded to 0
Returns:	single or list of times all in the same unit
Return type:	timelike or list

Parameters:	t0 (datetime64) – starting date tdelta (timedelta64) – timedelta to add unit (str) – If given, regularize all times to this unit, otherwise, regularize them to the most precise of the bunch most_precise (bool) – If True, then convert all times to the most precise unit that fits all the times allow_rounding (bool) – if tdelta is too small to be represented in the same unit as t0, then round it to the finest unit that fits both t0 and tdelta allow0 (bool) – If False, and a value is rounded to 0 in a given unit, then raise a PrecisionError
Returns:	t0 + tdelta
Return type:	datetime64

Parameters:	t1 (datetime64) – t1 for t1 - t2 t2 (datetime64) – t2 for t1 - t2 unit (str) – If given, regularize all times to this unit, otherwise, regularize them to the most precise of the bunch most_precise (bool) – If True, then convert all times to the most precise unit that fits all the times
Returns:	t1 - t2
Return type:	timedelta64

Parameters:	size – if given, set nchunks such that chunks have about ‘size’ elements
Returns:	nchunks slices of length N = (len(lst) // nchunks) or N - 1

viscid package¶

Subpackages¶

Submodules¶

Module contents¶