[SciPy-user] Error in scipy.spatial.cKDTree

Rob Hetland hetland@tamu....
Sat Feb 21 23:35:12 CST 2009


On Feb 21, 2009, at 10:07 PM, Anne Archibald wrote:

>
> Oops! Fixed in SVN r5585.
>
> The error happens when the query array is not "contiguous"; the
> easiest way to trigger it is to do a query with a transposed array;
> the query coordinates will be scrambled. As a workaround, just apply
> np.ascontiguousarray() to any query.



Thanks.  I tried the pure python version in the meantime, and it was  
surprisingly fast for pure python ..  Hopefully the c version will be  
even faster still (and give the right answer to boot!).

For posterity, here is some code that fills in sparse arrays that  
should not be sparse, below, for when delaunay interpolation is  
overkill.  (what I was working on when I found the bug).

-Rob




import numpy as np
from scipy.spatial import cKDTree

def fill_nearest(x, fill_value=0.0):
    '''Fill missing values in an array with an average of nearest  
neighbors.'''

    assert x.ndim == 2, 'x must be a 2D array.'

    # Create (i, j) point arrays for good and bad data.
    # Bad data is marked by the fill_value, good data elsewhere.
    igood = np.vstack(np.where(x!=fill_value)).T
    ibad = np.vstack(np.where(x==fill_value)).T

    # create a tree for the bad points, the points to be filled
    # ann_tree = ann.kd_tree(igood)
    tree = cKDTree(igood)

    # get the four closest points to the bad points
    # iquery, dist = ann_tree.search(ibad, k=4)  # here, distance is  
squared
    dist, iquery = tree.query(ibad, k=4, p=2)

    # create a weight normalized the nearest points are weighted as 1.
    # points greater than one are then set to zero.
    weight = dist/(dist.min(axis=1)[:, newaxis] * ones_like(dist))
    weight[weight > 1] = 0

    # multiply the queried good points by the weight, selecting only  
the near
    # points.  Divide by the number of nearest points to get average.
    xfill = weight * x[igood[:,0][iquery], igood[:,1][iquery]]
    xfill = (xfill/weight.sum(axis=1)[:, newaxis]).sum(axis=1)

    # place average of nearest good points, xfill, into bad point  
locations.
    x[ibad[:,0], ibad[:,1]] = xfill

    return x




----
Rob Hetland, Associate Professor
Dept. of Oceanography, Texas A&M University
http://pong.tamu.edu/~rob
phone: 979-458-0096, fax: 979-845-6331





More information about the SciPy-user mailing list