How to interpolated irregularly distributed data on a non-linear grid in python?
<p>I got a dataset of around 80.000 points (x,y,z), while the points are irregularly distributed in the (x,y) \in [0,a]x[0,b] plane and at every point (x,y) the physical quantity z takes a certain value. To further evaluate the data I want to interpolate it on a grid.</p>

<p>Before I used scipy.interpolate.griddata to successfully interpolated them on a regular, quadratic, 2D grid. However this regular grid has the disadvantage, that it can't model the regions with drastic change in z appropriately, while there are to many data points in regions with only slight change in z.</p>

<p>I would like to have a non-linear (preferably still quadratic, but with variable mesh size) grid, with more grid points in regions of drastic change in z and less data points in regions of slight change in z. </p>

<p>I think you have it backwards: your grid can be as regular as can be but each grid point should be evaluated using the same number of sample points, thus allowing for strong gradient changes in regions of high sample density, and imposing smoothness in data sparse regions. </p>

<p>I use inverse-distance weighted trees for that.
The implementation that I have floating around in python:
    import numpy as np
    from scipy.spatial import cKDTree</p>

<pre><code>class invdisttree(object):
    """
    Compute the score of query points based on the scores of their k-nearest neighbours,
    weighted by the inverse of their distances.

    @reference:
    https://en.wikipedia.org/wiki/Inverse_distance_weighting

    Example:
    --------

    import numpy as np
    import matplotlib.pyplot as plt
    from invdisttree import invdisttree

    import matplotlib.pyplot as plt

    # create sample points with structured scores
    X1 = 10 * np.random.rand(1000, 2) -5

    def func(x, y):
        return np.sin(x**2 + y**2) / (x**2 + y**2)

    z1 = func(X1[:,0], X1[:,1])

    # 'train'
    tree = invdisttree(X1, z1)

    # 'test'
    spacing = np.linspace(-5., 5., 100)
    X2 = np.meshgrid(spacing, spacing)
    grid_shape = X2[0].shape
    X2 = np.reshape(X2, (2, -1)).T
    z2 = tree(X2)

    fig, (ax1, ax2, ax3) = plt.subplots(1,3, sharex=True, sharey=True, figsize=(10,3))
    ax1.contourf(spacing, spacing, func(*np.meshgrid(spacing, spacing)))
    ax1.set_title('Ground truth')
    ax2.scatter(X1[:,0], X1[:,1], c=z1, linewidths=0)
    ax2.set_title('Samples')
    ax3.contourf(spacing, spacing, z2.reshape(grid_shape))
    ax3.set_title('Reconstruction')
    plt.show()

    """
    def __init__(self, X=None, z=None, leafsize=10):
        if not X is None:
            self.tree = cKDTree(X, leafsize=leafsize )
        if not z is None:
            self.z = z

    def fit(self, X=None, z=None, leafsize=10):
        """
        Arguments:
        ----------
            X: (N, d) ndarray
                Coordinates of N sample points in a d-dimensional space.
            z: (N,) ndarray
                Corresponding scores.
            leafsize: int (default 10)
                Leafsize of KD-tree data structure;
                should be less than 20.

        Returns:
        --------
            invdisttree instance: object
        """
        return self.__init__(X, z, leafsize)

    def __call__(self, X, k=6, eps=1e-6, p=2, regularize_by=1e-9):
        self.distances, self.idx = self.tree.query(X, k, eps=eps, p=p)
        self.distances += regularize_by
        weights = self.z[self.idx.ravel()].reshape(self.idx.shape)
        mw = np.sum(weights/self.distances, axis=1) / np.sum(1./self.distances, axis=1)
        return mw

    def transform(self, X, k=6, p=2, eps=1e-6, regularize_by=1e-9):
        """
        Arguments:
        ----------
            X: (N, d) ndarray
                Coordinates of N query points in a d-dimensional space.

            k: int (default 6)
                Number of nearest neighbours to use.

            p: int or inf
                Which Minkowski p-norm to use.
                1 is the sum-of-absolute-values "Manhattan" distance
                2 is the usual Euclidean distance
                infinity is the maximum-coordinate-difference distance

            eps: float (default 1e-6)
                Return approximate nearest neighbors; the k-th returned value
                is guaranteed to be no further than (1+eps) times the
                distance to the real k-th nearest neighbor.

            regularise_by: float (default 1e-9)
                Regularise distances to prevent division by zero
                for sample points with the same location as query points.

        Returns:
        --------
            z: (N,) ndarray
                Corresponding scores.
        """
        return self.__call__(X, k, eps, p, regularize_by)
</code></pre>

python grid interpolation