Encontrar el índice del punto más cercano en matrices numerosas de coordenadas xey

Tengo dos matrices numéricas 2d: x_array contiene información posicional en la dirección x, y_array contiene posiciones en la dirección y.

Entonces tengo una larga lista de puntos x, y.

Para cada punto de la lista, necesito encontrar el índice de matriz de la ubicación (especificada en las matrices) que está más cerca de ese punto.

Ingenuamente he producido un código que funciona, basado en esta pregunta: Encuentre el valor más cercano en la matriz numpy

es decir

import time
import numpy

def find_index_of_nearest_xy(y_array, x_array, y_point, x_point):
    distance = (y_array-y_point)**2 + (x_array-x_point)**2
    idy,idx = numpy.where(distance==distance.min())
    return idy[0],idx[0]

def do_all(y_array, x_array, points):
    store = []
    for i in xrange(points.shape[1]):
        store.append(find_index_of_nearest_xy(y_array,x_array,points[0,i],points[1,i]))
    return store


# Create some dummy data
y_array = numpy.random.random(10000).reshape(100,100)
x_array = numpy.random.random(10000).reshape(100,100)

points = numpy.random.random(10000).reshape(2,5000)

# Time how long it takes to run
start = time.time()
results = do_all(y_array, x_array, points)
end = time.time()
print 'Completed in: ',end-start

Estoy haciendo esto en un gran conjunto de datos y realmente me gustaría acelerarlo un poco. ¿Alguien puede optimizar esto?

Gracias.

ACTUALIZACIÓN: SOLUCIÓN siguiendo las sugerencias de @silvado y @justin (abajo)

# Shoe-horn existing data for entry into KDTree routines
combined_x_y_arrays = numpy.dstack([y_array.ravel(),x_array.ravel()])[0]
points_list = list(points.transpose())


def do_kdtree(combined_x_y_arrays,points):
    mytree = scipy.spatial.cKDTree(combined_x_y_arrays)
    dist, indexes = mytree.query(points)
    return indexes

start = time.time()
results2 = do_kdtree(combined_x_y_arrays,points_list)
end = time.time()
print 'Completed in: ',end-start

Este código anterior aceleró mi código (buscando 5000 puntos en matrices de 100x100) 100 veces. Curiosamente, usar scipy.spatial.KDTree (en lugar de scipy.spatial.cKDTree ) dio tiempos comparables a mi solución ingenua, por lo que definitivamente vale la pena usar la versión cKDTree ...

scipy.spatialTambién tiene una aplicación árbol kd: scipy.spatial.KDTree.

Generalmente, el enfoque consiste en utilizar primero los datos puntuales para construir un árbol kd. La complejidad computacional de eso es del orden de N log N, donde N es el número de puntos de datos. Las consultas de rango y las búsquedas de vecinos más cercanos se pueden realizar con una complejidad de log N. Esto es mucho más eficiente que simplemente recorrer todos los puntos (complejidad N).

Por lo tanto, si tiene consultas repetidas de rango o vecino más cercano, se recomienda encarecidamente un árbol kd.

Aquí hay un scipy.spatial.KDTreeejemplo

In [1]: from scipy import spatial

In [2]: import numpy as np

In [3]: A = np.random.random((10,2))*100

In [4]: A
Out[4]:
array([[ 68.83402637,  38.07632221],
       [ 76.84704074,  24.9395109 ],
       [ 16.26715795,  98.52763827],
       [ 70.99411985,  67.31740151],
       [ 71.72452181,  24.13516764],
       [ 17.22707611,  20.65425362],
       [ 43.85122458,  21.50624882],
       [ 76.71987125,  44.95031274],
       [ 63.77341073,  78.87417774],
       [  8.45828909,  30.18426696]])

In [5]: pt = [6, 30]  # <-- the point to find

In [6]: A[spatial.KDTree(A).query(pt)[1]] # <-- the nearest point 
Out[6]: array([  8.45828909,  30.18426696])

#how it works!
In [7]: distance,index = spatial.KDTree(A).query(pt)

In [8]: distance # <-- The distances to the nearest neighbors
Out[8]: 2.4651855048258393

In [9]: index # <-- The locations of the neighbors
Out[9]: 9

#then 
In [10]: A[index]
Out[10]: array([  8.45828909,  30.18426696])

Si puede masajear sus datos en el formato correcto, una forma rápida de hacerlo es usar los métodos en scipy.spatial.distance:

http://docs.scipy.org/doc/scipy/reference/spatial.distance.html

En particular, pdisty cdistproporciona formas rápidas de calcular distancias por pares.