[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