[SciPy-User] Interpolation in 3D with interp2d

[denis]

On Aug 12, 10:35 pm, Jana Schulz <schraba... at web.de> wrote:
> Hi,
> RectBivariateSpline() can not handle randomly spaced data (that's what the error message says). The option 'linear' in griddata works only for  constant spacing data. I attached the surface plot of my data and the txt-file including my data.
>
> As you see the surface plot works pretty well. All I want is to get any z-value  for a given x,y-position.

Jana,
  griddata( your data ) with the default interp="nn" looks good
*once you take log10 of all your data*;
any interpolator can misbehave on data from 1e-6 .. 1e6.
I sent sample code and a plot to you at web.de on 16Aug
(since my mails via gmane aren't getting through); did it make sense ?

If you want to interpolate a single point at a time with Delaunay
triangulation
with a 2-stage interpolator like interp2d, see Triinterpolate below.
(Yes the many interpolators in scipy are, hmm, scattered, nonuniform.)
Is this closer to what you want ?

cheers
  -- denis

# tri = Triinterpolate(x,y,z),  tri( single points )

from __future__ import division
import sys
import numpy as np
import matplotlib.delaunay as delaunay  # $matplotlib/delaunay/
interpolate.py

__date__ = "2010-08-17 Aug"


class Triinterpolate:
    """
    interpolate scattered data a point at a time:
        tri = Triinterpolate( x,y,z )  -- 1d arrays
            first build a Delaunay triangulation; then
        zi = tri(xi,yi)  -- interpolates one point
            NaN if xi,yi is outside the convex hull of the x,y points

    See also
    http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
    http://en.wikipedia.org/wiki/Natural_neighbor
    .../matplotlib/ mlab.py griddata(), delaunay/interpolate.py
    """
    def __init__( self, x,y,z ):
        if not (len(x) == len(y) == len(z)
        and  x.ndim == y.ndim == z.ndim == 1):
            raise TypeError("inputs x,y,z must all be 1D arrays of the
same length")
        tri = delaunay.Triangulation(x,y)
        self.interp = tri.nn_interpolator(z)  # linear_ has no
__call__

    def __call__( self, xi,yi ):
        if not (np.isscalar(xi) and np.isscalar(yi)):
            raise TypeError("inputs xi,yi must be scalars")
        return self.interp( xi, yi )  # NaN outside convex hull, cf
corners()

#...............................................................................
N = 20  # N random in side x side
side = 10

exec "\n".join( sys.argv[1:] )  # run this.py N= ...
np.random.seed(1)
np.set_printoptions( 1, threshold=100, suppress=True )  # .1f

def terrainf( x, y ):
    return x**2 + y**2

x,y = np.random.uniform( 0, side, (2,N) )
z = terrainf( x, y )
print "x y z:\n", x, "\n", y, "\n", z, "\n"

tri = Triinterpolate( x,y,z )
print "Triinterpolate a point at a time:"
for yi in range(side):
    for xi in range(side):
        print "%5.0f" % tri( xi,yi ) ,
    print ""