[SciPy-User] 2-D data interpolation

[Pauli Virtanen]

06.06.2012 21:13, Darcoux Christine kirjoitti:
[clip]
> I tried this function and it works very similar to the interp2
> function of matlab. The only problem is that my data are very peaky
> and the cubic version of griddata gives me big over/under-shoots. The
> matlab function does not gives me such behavior.
> 
> Is there another high 2D order interpolator that could be used in this case ?
> 
> Octave has a nice "pchip" option for interp2 that use piecewise cubic
> Hermite interpolating polynomial, but unfortunatly griddata does not
> have such option.

Ah, your data is actually on a regular grid rather than scattered. This
means n-dim interpolation can be built up from 1-D interpolation as a
tensor product.

Scipy does have `pchip` for 1-D, but currently it's a bit clumsy to use
for 2-D data. You can build up tensor product 2-D pchip interpolation
like this:

----

import numpy as np
from scipy.interpolate import pchip

def pchip2(x, y, z, xi, yi):
    """P-chip interpolation on 2-D. x, y, xi, yi should be 1-D
    and z.shape == (len(x), len(y))"""

    z2 = np.empty([len(xi), len(y)], dtype=z.dtype)
    for k in xrange(len(y)):
        z2[:,k] = pchip(x, z[:,k])(xi)

    z3 = np.empty([len(xi), len(yi)], dtype=z.dtype)
    for k in xrange(len(xi)):
        z3[k,:] = pchip(y, z2[k,:])(yi)

    return z3

# Random test case
import matplotlib.pyplot as plt

x = 2*pi*np.random.rand(20)
x.sort()

y = 2*pi*np.random.rand(10)
y.sort()

z = np.cos(x[:,None]) * np.sin(y)

plt.figure(1)
plt.pcolor(x, y, z.T)

xi = np.linspace(x.min(), x.max(), 200)
yi = np.linspace(y.min(), y.max(), 300)
zi = pchip2(x, y, z, xi, yi)

plt.figure(2)
plt.pcolor(xi, yi, zi.T)

----

It's not the fastest thing in the world, but maybe helps? The same idea
generalizes to 3 and any higher dimensions.

The pchip idea does not directly generalize to true *scattered* data
interpolation, which is what `griddata` does. There, the problem is
quite a bit more tricky, but likely there is a publication out there
that discusses how to do monotonicity-preserving interpolation on an
arbitrary network of triangles...

-- 
Pauli Virtanen