Purpose: interpolate data given on an N-dimensional rectangular grid, uniform or non-uniform, with the fast `scipy.ndimage.map_coordinates`

. Non-uniform grids are first uniformized with `numpy.interp`

.

Background: the reader should know some Python and NumPy (IPython is invaluable for learning both). For basics of interpolation, see Bilinear interpolation on Wikipedia. For `map_coordinates`

, see the example under multivariate-spline-interpolation-in-python-scipy .

Say we have rainfall on a 4 x 5 grid of rectangles, lat 52 .. 55 x lon -10 .. -6, and want to interpolate (estimate) rainfall at 1000 query points in between the grid points.

```
from intergrid import Intergrid # intergrid.py in $PYTHONPATH
# define the grid --
griddata = np.loadtxt(...) # griddata.shape == (4, 5)
lo = np.array([ 52, -10 ]) # lowest lat, lowest lon
hi = np.array([ 55, -6 ]) # highest lat, highest lon
# set up an interpolator function "iinterfunc()" with class Intergrid --
interfunc = Intergrid( griddata, lo=lo, hi=hi )
# generate 1000 random query points, lo <= [lat, lon] <= hi --
query_points = lo + np.random.uniform( size=(1000, 2) ) * (hi - lo)
# get rainfall at the 1000 query points --
query_values = interfunc.at( query_points ) # -> 1000 values
```

What this does: for each [lat, lon] in query_points:

- find the square of
`griddata`

it's in, e.g. [52.5, -8.1] -> [0, 3][0, 4] [1, 4][1, 3] - do bilinear (multilinear) interpolation in that square, using
`scipy.ndimage.map_coordinates`

.

Check:

`interfunc( lo ) == griddata[0, 0]`

`interfunc( hi ) == griddata[-1, -1]`

i.e. `griddata[3, 4]`

`griddata`

: numpy array_like, 2d 3d 4d ...`lo, hi`

: user coordinates of the corners of griddata, 1d array-like, lo < hi`maps`

: an optional list of `dim`

descriptors of piecewise-linear or nonlinear maps,

e.g. [[50, 52, 62, 63], None] # uniformize lat, linear lon; see below`copy`

: make a copy of query_points, default `True`

;

`copy=False`

overwrites query_points, runs in less memory`verbose`

: the default 1 prints a summary of each call, with run time`order`

: interpolation order:

default 1: bilinear, trilinear ... interpolation using all 2^dim corners

0: each query point -> the nearest grid point -> the value there

2 to 5: spline, see below`prefilter`

: the kind of spline:

default `False`

: smoothing B-spline

`True`

: exact-fit C-R spline

1/3: Mitchell-Netravali spline, 1/3 B + 2/3 fit

After setting up `interfunc = Intergrid(...)`

, either

```
query_values = interfunc.at( query_points ) # or
query_values = interfunc( query_points )
```

do the interpolation. (The latter is `__call__`

in python.)

What if our griddata above is at non-uniformly-spaced latitudes, say [50, 52, 62, 63] ? `Intergrid`

can "uniformize" these before interpolation, like this:

```
lo = np.array([ 50, -10 ])
hi = np.array([ 60, -6 ])
maps = [[50, 52, 62, 63], None] # uniformize lat, linear lon
interfunc = Intergrid( griddata, lo=lo, hi=hi, maps=maps )
```

This will map (transform, stretch, warp) the lats in query_points column 0 to array coordinates in the range 0 .. 3, using `np.interp`

to do piecewise-linear (PWL) mapping:

```
50 51 52 53 54 55 56 57 58 59 60 61 62 63 # lo[0] .. hi[0]
0 .5 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 3
```

`maps[1] None`

says to map the lons in query_points column 1 linearly:

```
-10 -9 -8 -7 -6 # lo[1] .. hi[1]
0 1 2 3 4
```

The query_points are first clipped, then columns mapped linearly or non-linearly, then fed to `map_coordinates`

.

In `dim`

dimensions (i.e. axes or columns), `lo`

and `hi`

are each `dim`

numbers, the low and high corners of the data grid.`maps`

is an optional list of `dim`

map descriptors, which can be

`None`

: linear-transform that column,`query_points[:,j]`

, to`griddata`

:

`lo[j] -> 0`

`hi[j] -> griddata.shape[j] - 1`

- a callable function: e.g.
`np.log`

does

`query_points[:,j] = np.log( query_points[:,j] )`

- a
*sorted*array describing a non-uniform grid:

`query_points[:,j] =`

`np.interp( query_points[:,j], [50, 52, 62, 63], [0, 1, 2, 3] )`

```
git clone http://github.com/denis-bz/interpol.git
# rm -rf interpol/barypol
# add .../interpol/intergrid to PYTHONPATH in ~/.bashrc or ~/.cshrc
```

`Intergrid( ... order = 0 to 5 )`

gives the spline order:

`order=1`

, the default, does bilinear, trilinear ... interpolation, which looks at the grid data at all 4 8 16 .. 2^dim corners of the box around each query point.

`order=0`

looks at only the one gridpoint nearest each query point — crude but fast.

`order = 2 to 5`

does spline interpolation on a uniform or uniformized grid, looking at (order+1)^dim neighbors of each query point.

`Intergrid( ... prefilter = False | True | 1/3 )`

specifies the kind of spline, for `order >= 2`

:

`prefilter=0`

or `False`

: B-spline

`prefilter=1`

or `True`

: exact-fit spline

`prefilter=1/3 )`

: M-N spline.

A B-spline goes through smoothed data points, with [1 4 1] smoothing, [0 0 1 0 0] -> [0 1 4 1 0] / 6.

An exact-fit a.k.a interpolating spline goes through the data points exactly. This is not what you want for noisy data, and may also wiggle or overshoot more than B-splines do.

An M-N spline blends 1/3 B-spline and 2/3 exact-fit; see Mitchell and Netravali, Reconstruction filters in computer-graphics , 1988, and the plots from `intergrid/test/MNspline.py`

.

Exact-fit or interpolating splines can be local or global. Catmull-Rom splines and the original M-N splines are local: they look at 4 neighbors of each query point in 1d, 16 in 2d, 64 in 3d. Prefiltering looks at all neighbors, so is a bit smoother — although I don't know of test images that show a visible difference. Confusingly, the term "Cardinal spline" is sometimes used for local (C-R, FIR), and sometimes for global (IIR prefilter, then B-spline).

Prefiltering is a clever transformation such that `Bspline( transform( data )) = exactfitspline( data )`

. It is described in a paper by M. Unser, Splines: A perfect fit for signal and image processing , 1999.

Uniformizing a grid with PWL, then uniform-splining, is fast and simple, but not as smooth as true splining on the original non-uniform grid. The differences will of course depend on the grid spacings and on how rough the function is.

Run any interpolator on *your* data with orders 0, 1 ... to get an idea of how the results get smoother, and take longer. Check a few query points by hand; plot some cross-sections.

`griddata`

values can be of any numpy integer or floating type: int8 uint8 .. int32 int64 float32 float64. `np.float32`

will use less memory than `np.float64`

(but beware of functions in the flow that silently convert everything to float64). The values must be numbers, not vectors.

Coordinate scaling doesn't matter to `Intergrid`

; corner weights are calculated in unit cubes of `griddata`

, after scaling and mapping. If for example griddata column 3 is multiplied by 1000, and lo[3] hi[3] too, the weights are unchanged.

Box grids get big and slow above 5d. A cube with steps 0 .1 .2 .. 1.0 in all dimensions has 11^6 ~ 1.8M points in 6d, 11^8 ~ 200M in 8d. One can reduce that only with a coarser grid like 0 .5 1 in some dimensions (those that vary the least). But time ~ 2^d per query point grows pretty fast.

`map_coordinates`

in 5d looks at 32 corner values, with average weight 3 %. If the weights are roughly equal (which they will tend to be, by the central limit theorem ?), sharp edges or gradients will be blurred, and colors mixed to a grey fog.

Terminology varies, so the basic kinds of box grids a.k.a. rectangular grids are defined here.

An integer or Cartesian grid has integer coordinates, e.g. 2 x 3 x 5 points in a numpy array: `A = np.array((2,3,5)); A[0,0,0], A[0,0,1] .. A[1,2,4]`

.

A uniform box grid has nx x ny x nz ... points uniformly spaced, linspace x linspace x linspace ... so all boxes have the same size and are axis-aligned. Examples: 1024 x 768 pixels on a screen, or 4 x 5 points at latitudes [10 20 30 40] x longitudes [-10 -9 -8 -7 -6].

A non-uniform box grid also has nx x ny x nz ... points, but allows non-uniform spacings, e.g. latitudes [-10 0 60 70] x longitudes [-10 -9 0 20 40]; the boxes have different sizes but are still axis-aligned.

(Scattered data, as the name says, has points anywhere, not only on grid lines. To interpolate scattered data in `scipy`

, see scipy.interpolate.griddata and scipy.spatial.cKDTree .)

There are countless varieties of grids: grids with holes, grids warped to various map projections, multiscale / multiresolution grids ... Google "regridding" or "resampling".

See intergrid/test/test-4d.py: a 4d grid with 1M scattered query points, uniform / non-uniform box grid, on a 2.5Gz i5 iMac:

```
shape (361, 720, 47, 8) 98M * 8
Intergrid: 617 msec 1000000 points in a (361, 720, 47, 8) grid 0 maps order 1
Intergrid: 788 msec 1000000 points in a (361, 720, 47, 8) grid 4 maps order 1
```

scipy.ndimage.interpolation.map_coordinates

scipy reference ndimage

stackoverflow.com/questions/tagged/scipy+interpolation

and testcases most welcome

— denis-bz-py at t-online dot de

Last change: 2013-06-04 Jun prefilter Bspline / exact-fit / M-N