Python rounding of random floats to nearest points on a 2D uniform mesh grid

Question

Despite numpy & scipy's many rounding functions, I cannot find one that allows me to discretize randomized floats with respect to nodes in a 2D uniform mesh grid. For example,

# create mesh grid
n = 11
l = 16.
x = np.linspace(-l/2, l/2, n)
y = np.linspace(-l/2, l/2, n)
xx, yy = np.meshgrid(x, y, sparse=True)

>>> xx
array([-8. , -6.4, -4.8, -3.2, -1.6,  0. ,  1.6,  3.2,  4.8,  6.4,  8. ])    
>>> yy
array([[-8. ],
       [-6.4],
       [-4.8],
       [-3.2],
       [-1.6],
       [ 0. ],
       [ 1.6],
       [ 3.2],
       [ 4.8],
       [ 6.4],
       [ 8. ]])

If I have m number of normally distributed random float a=np.random.multivariate_normal([0,0], [[l/2,0],[0,l/2]], m), how can I round them to the nearest mesh nodes (assuming periodic boundaries)?

In your example, the grid spacing is uniform. Is that always the case? — Warren Weckesser, Jul 21 '15 at 15:14
Possible duplicate of https://stackoverflow.com/questions/10818546/finding-index-of-nearest-point-in-numpy-arrays-of-x-and-y-coordinates — rth, Jul 21 '15 at 15:27
In scipy there is a nearest-neighbor interpolation method that might be adapted to this problem. — hpaulj, Jul 21 '15 at 19:46
@rth: your question is more general (no uniformity requirement). The answer there are useful. But this question might have a specialized answer. — jfs, Jul 21 '15 at 20:56

score 4 · Accepted Answer · edited May 23 '17 at 11:58

In another SO question I helped the questioner use a scipy nearest neigbhor interpolator.

Repeating Scipy's griddata

Working from that I figured out a solution to your problem. In scipy.interpolate I find a NearestNDInterpolator. From that I find that scipy.spatial has a method for finding nearest neighbors:

In [936]: from scipy import interpolate
In [937]: interpolate.NearestNDInterpolator?
...
Docstring:
NearestNDInterpolator(points, values)

Nearest-neighbour interpolation in N dimensions.
...
Uses ``scipy.spatial.cKDTree``

In [938]: from scipy import spatial
In [939]: spatial.cKDTree?

cKDTree is used in 2 steps; create a tree from the grid, and query the nearest neighbors.

From your meshgrid I can create a (n,2) array of points

In [940]: xx, yy = np.meshgrid(x, y)
In [941]: xygrid=np.array([xx,yy])
In [942]: xygrid=xygrid.reshape(2,-1).T  # clumsy

create a search tree:

In [943]: tree=spatial.cKDTree(xygrid)

test set of points, (10,2):

In [944]: a=np.random.multivariate_normal([0,0], [[l/2,0],[0,l/2]],10)

Query of the search tree gives 2 arrays, one of distances from nearest neighbor, and one with the indexes:

In [945]: I=tree.query(a)

In [946]: I
Out[946]: 
(array([ 0.70739099,  0.9894934 ,  0.44489157,  0.3930144 ,  0.273121  ,
        0.3537348 ,  0.32661876,  0.55540787,  0.58433421,  0.538722  ]),
 array([61, 72, 85, 72, 82, 39, 38, 62, 25, 59]))

Comparing the a points, with the xygrid nearest neighbor grid points, shows that they indeed appear to be close. A scatter point plot would be better.

In [947]: a
Out[947]: 
array([[ 1.44861113, -0.69100176],
       [ 1.00827575,  0.80693026],
       [ 4.37200745,  3.07854676],
       [ 1.2193471 ,  1.50220587],
       [ 0.12668563,  2.95803754],
       [ 1.4758331 , -3.53122635],
       [ 0.28425494, -3.03913067],
       [ 2.8203361 ,  0.40538034],
       [-3.67726571, -4.46285921],
       [-1.07228578, -0.10834709]])

In [948]: xygrid[I[1],:]
Out[948]: 
array([[ 1.6,  0. ],
       [ 1.6,  1.6],
       [ 4.8,  3.2],
       [ 1.6,  1.6],
       [ 0. ,  3.2],
       [ 1.6, -3.2],
       [ 0. , -3.2],
       [ 3.2,  0. ],
       [-3.2, -4.8],
       [-1.6,  0. ]])

A solution in rth's link also uses cKDTree. I'm just filling in the details on how to work from your griddata. Finding index of nearest point in numpy arrays of x and y coordinates

Python rounding of random floats to nearest points on a 2D uniform mesh grid

1 Answers1