tracers module

Tracer particles in a vector field

• Requires scipy.interpolate

Module Summary

functions:

tracers.random_particles

Return an array of count 3d vertices of random particle positions Minimum and maximum values defined by lowerbound and upperbound

classes:

tracers.RegularGridInterpolator	Interpolator on a regular or rectilinear grid in arbitrary dimensions.
tracers.Tracers

Module Details

functions:

random_particles(count, lowerbound=[0, 0, 0], upperbound=[1, 1, 1], dims=3)

Return an array of count 3d vertices of random particle positions Minimum and maximum values defined by lowerbound and upperbound

classes:

class RegularGridInterpolator(points, values, method='linear', bounds_error=True, fill_value=nan, *, solver=None, solver_args=None)

Bases: object

Interpolator on a regular or rectilinear grid in arbitrary dimensions.

The data must be defined on a rectilinear grid; that is, a rectangular grid with even or uneven spacing. Linear, nearest-neighbor, spline interpolations are supported. After setting up the interpolator object, the interpolation method may be chosen at each evaluation.

Parameters:

• points (tuple of ndarray of float, with shapes (m1, ), ..., (mn, )) – The points defining the regular grid in n dimensions. The points in each dimension (i.e. every elements of the points tuple) must be strictly ascending or descending.

• values (array_like, shape (m1, ..., mn, ...)) – The data on the regular grid in n dimensions. Complex data is accepted.

• method (str, optional) – The method of interpolation to perform. Supported are “linear”, “nearest”, “slinear”, “cubic”, “quintic” and “pchip”. This parameter will become the default for the object’s __call__ method. Default is “linear”.

• bounds_error (bool, optional) – If True, when interpolated values are requested outside of the domain of the input data, a ValueError is raised. If False, then fill_value is used. Default is True.

• fill_value (float or None, optional) – The value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated. Default is np.nan.

• solver (callable, optional) –

  Only used for methods “slinear”, “cubic” and “quintic”. Sparse linear algebra solver for construction of the NdBSpline instance. Default is the iterative solver scipy.sparse.linalg.gcrotmk.

  Added in version 1.13.

• solver_args (dict, optional) –

  Additional arguments to pass to solver, if any.

  Added in version 1.13.

__call__()

grid

The points defining the regular grid in n dimensions. This tuple defines the full grid via np.meshgrid(*grid, indexing='ij')

Type:

tuple of ndarrays

values

Data values at the grid.

Type:

ndarray

method

Interpolation method.

Type:

str

fill_value

Use this value for out-of-bounds arguments to __call__.

Type:

float or None

bounds_error

If True, out-of-bounds argument raise a ValueError.

Type:

bool

Notes

Contrary to LinearNDInterpolator and NearestNDInterpolator, this class avoids expensive triangulation of the input data by taking advantage of the regular grid structure.

In other words, this class assumes that the data is defined on a rectilinear grid.

Added in version 0.14.

The ‘slinear’(k=1), ‘cubic’(k=3), and ‘quintic’(k=5) methods are tensor-product spline interpolators, where k is the spline degree, If any dimension has fewer points than k + 1, an error will be raised.

Added in version 1.9.

If the input data is such that dimensions have incommensurate units and differ by many orders of magnitude, the interpolant may have numerical artifacts. Consider rescaling the data before interpolating.

Choosing a solver for spline methods

Spline methods, “slinear”, “cubic” and “quintic” involve solving a large sparse linear system at instantiation time. Depending on data, the default solver may or may not be adequate. When it is not, you may need to experiment with an optional solver argument, where you may choose between the direct solver (scipy.sparse.linalg.spsolve) or iterative solvers from scipy.sparse.linalg. You may need to supply additional parameters via the optional solver_args parameter (for instance, you may supply the starting value or target tolerance). See the scipy.sparse.linalg documentation for the full list of available options.

Alternatively, you may instead use the legacy methods, “slinear_legacy”, “cubic_legacy” and “quintic_legacy”. These methods allow faster construction but evaluations will be much slower.

Examples

Evaluate a function on the points of a 3-D grid

As a first example, we evaluate a simple example function on the points of a 3-D grid:

>>> from scipy.interpolate import RegularGridInterpolator
>>> import numpy as np
>>> def f(x, y, z):
...     return 2 * x**3 + 3 * y**2 - z
>>> x = np.linspace(1, 4, 11)
>>> y = np.linspace(4, 7, 22)
>>> z = np.linspace(7, 9, 33)
>>> xg, yg ,zg = np.meshgrid(x, y, z, indexing='ij', sparse=True)
>>> data = f(xg, yg, zg)

data is now a 3-D array with data[i, j, k] = f(x[i], y[j], z[k]). Next, define an interpolating function from this data:

>>> interp = RegularGridInterpolator((x, y, z), data)

Evaluate the interpolating function at the two points (x,y,z) = (2.1, 6.2, 8.3) and (3.3, 5.2, 7.1):

>>> pts = np.array([[2.1, 6.2, 8.3],
...                 [3.3, 5.2, 7.1]])
>>> interp(pts)
array([ 125.80469388,  146.30069388])

which is indeed a close approximation to

>>> f(2.1, 6.2, 8.3), f(3.3, 5.2, 7.1)
(125.54200000000002, 145.894)

Interpolate and extrapolate a 2D dataset

As a second example, we interpolate and extrapolate a 2D data set:

>>> x, y = np.array([-2, 0, 4]), np.array([-2, 0, 2, 5])
>>> def ff(x, y):
...     return x**2 + y**2

>>> xg, yg = np.meshgrid(x, y, indexing='ij')
>>> data = ff(xg, yg)
>>> interp = RegularGridInterpolator((x, y), data,
...                                  bounds_error=False, fill_value=None)

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(projection='3d')
>>> ax.scatter(xg.ravel(), yg.ravel(), data.ravel(),
...            s=60, c='k', label='data')

Evaluate and plot the interpolator on a finer grid

>>> xx = np.linspace(-4, 9, 31)
>>> yy = np.linspace(-4, 9, 31)
>>> X, Y = np.meshgrid(xx, yy, indexing='ij')

>>> # interpolator
>>> ax.plot_wireframe(X, Y, interp((X, Y)), rstride=3, cstride=3,
...                   alpha=0.4, color='m', label='linear interp')

>>> # ground truth
>>> ax.plot_wireframe(X, Y, ff(X, Y), rstride=3, cstride=3,
...                   alpha=0.4, label='ground truth')
>>> plt.legend()
>>> plt.show()

Other examples are given in the tutorial.

See also

NearestNDInterpolator

Nearest neighbor interpolator on unstructured data in N dimensions

LinearNDInterpolator

Piecewise linear interpolator on unstructured data in N dimensions

interpn

a convenience function which wraps RegularGridInterpolator

scipy.ndimage.map_coordinates

interpolation on grids with equal spacing (suitable for e.g., N-D image resampling)

References

[1]

Python package regulargrid by Johannes Buchner, see https://pypi.python.org/pypi/regulargrid/

[2]

Wikipedia, “Trilinear interpolation”, https://en.wikipedia.org/wiki/Trilinear_interpolation

[3]

Weiser, Alan, and Sergio E. Zarantonello. “A note on piecewise linear and multilinear table interpolation in many dimensions.” MATH. COMPUT. 50.181 (1988): 189-196. https://www.ams.org/journals/mcom/1988-50-181/S0025-5718-1988-0917826-0/S0025-5718-1988-0917826-0.pdf :doi:`10.1090/S0025-5718-1988-0917826-0`

__call__(xi, method=None, *, nu=None)

Interpolation at coordinates.

Parameters:

• xi (ndarray of shape (..., ndim)) – The coordinates to evaluate the interpolator at.

• method (str, optional) – The method of interpolation to perform. Supported are “linear”, “nearest”, “slinear”, “cubic”, “quintic” and “pchip”. Default is the method chosen when the interpolator was created.

• nu (sequence of ints, length ndim, optional) –

  If not None, the orders of the derivatives to evaluate. Each entry must be non-negative. Only allowed for methods “slinear”, “cubic” and “quintic”.

  Added in version 1.13.

Returns:

values_x – Interpolated values at xi. See notes for behaviour when xi.ndim == 1.

Return type:

ndarray, shape xi.shape[:-1] + values.shape[ndim:]

Notes

In the case that xi.ndim == 1 a new axis is inserted into the 0 position of the returned array, values_x, so its shape is instead (1,) + values.shape[ndim:].

Examples

Here we define a nearest-neighbor interpolator of a simple function

>>> import numpy as np
>>> x, y = np.array([0, 1, 2]), np.array([1, 3, 7])
>>> def f(x, y):
...     return x**2 + y**2
>>> data = f(*np.meshgrid(x, y, indexing='ij', sparse=True))
>>> from scipy.interpolate import RegularGridInterpolator
>>> interp = RegularGridInterpolator((x, y), data, method='nearest')

By construction, the interpolator uses the nearest-neighbor interpolation

>>> interp([[1.5, 1.3], [0.3, 4.5]])
array([2., 9.])

We can however evaluate the linear interpolant by overriding the method parameter

>>> interp([[1.5, 1.3], [0.3, 4.5]], method='linear')
array([ 4.7, 24.3])

class Tracers(grid, count=1000, lowerbound=None, upperbound=None, limit=None, age=4, respawn_chance=0.2, speed_multiply=1.0, height=0.0, viewer=None)

Bases: object

__init__(grid, count=1000, lowerbound=None, upperbound=None, limit=None, age=4, respawn_chance=0.2, speed_multiply=1.0, height=0.0, viewer=None)

Seed random particles into a vector field and trace their positions

Parameters:

• grid (list of coord arrays for each dimension as expected by RegularGridInterpolator,) – or a numpy array of 2d or 3d vertices, which will be converted before being sent to the interpolator Object returned from first call, pass None on first pass

• count (int) – Number of particles to seed and track

• lowerbound (optional minimum vertex point defining particle bounding box,) – if not provided will be taken from grid lower corner

• upperbound (optional maximum vertex point defining particle bounding box,) – if not provided will be taken from grid upper corner

• limit (float) – Distance limit over which tracers are not connected, For example if using a periodic boundary, setting limit to half the bounding box size will prevent tracer lines being connected when passing through the boundary

• age (int) – Minimum particle age in steps after which particle can be deleted and respawned, defaults to 4

• respawn (float) – Probability of respawning, after age reached, default 0.2 ==> 1 in 5 chance of deletion

• speed_multiply (float) – Speed multiplier, scaling factor for the velocity taken from the vector values

• height (float) – A fixed height value, all positions will be given this height as their Z component

• viewer (lavavu.Viewer) – Viewer object for plotting functions