Python: cómo acelerar la búsqueda espacial (punto más cercano)?

Question

Estoy trabajando en un programa que determina la ubicación más cercana de un punto dado. La nube de puntos estoy de pruebas en contra, es muy grande (unos 800.000). Sin embargo, esto realmente no se explica por qué mi aplicación es tan lento. Este es mi planteamiento:

En primer lugar, he creado un índice espacial para la forma de la punta

pntshp.ExecuteSQL('CREATE SPATIAL INDEX ON %s' % table_name)

He definido una matriz de distancias de separación para limitar el radio de búsqueda. Que, por supuesto, también significa que tengo que crear un buffer para cada punto (lo cual puede ser costoso).

BUFFER_DISTANCES = ( 0.001, 0.005, 0.01, 0.02, 0.05 ) # 100m, 500m, 1km, 2km, 5km

A continuación, el tampón se utiliza como un filtro espacial

node_lyr.SetSpatialFilter(buff)

Si el filtro devuelve None el búfer distancia será mayor.

for buffer_d in BUFFER_DISTANCES:
    buffr = get_buffer(xy_street,buffer_d)
    ...

Entonces yo soy el cálculo de la distancia a los puntos devueltos por el filtro espacial

p=ogr.Geometry(ogr.wkbPoint)
p.AddPoint(xy[0],xy[1])

for feat in node_lyr:
    geom = feat.GetGeometryRef()
    d = p.Distance(geom)
    dist.append(d)       

Para obtener el punto más cercano:

def get_closest_pnt(dist, node, how_close):
    mrg = zip(dist,node)
    mrg.sort(key=lambda t: t[0])
    try:
        return mrg[how_close]
    except IndexError, ierr:
        print '%s \ndist/node tuple contain %s' % (ierr,mrg)

Funciona todo bien, pero es muy lento. La creación de un índice espacial no muestran ningún efecto, de verdad. Para calcular los 100 puntos esta implementaciones de toma ~6,7 segundos. El programa debe ser capaz de calcular la ubicación más cercana para más de 2000 puntos tan rápido como sea posible. Alguna idea sobre cómo mejorar mi enfoque?

EDITAR

He probado diferentes enfoques para ver a dónde me lleva. Me encontré con algo muy sorprendente quiero compartir aquí.

He implementado un sencillo algoritmo de búsqueda como se describe aquí, y una de las soluciones que donde se sugiere (el conjunto ordenado de aproximación).

El hecho sorprendente es que el rendimiento no sólo depende de la aplicación, pero más aún de la OSX. Mi original ogr/buffer algoritmo resulta ser rápida ardiente en mi OSX que es laborioso lento en Linux (por lo tanto, la pregunta aquí).

Aquí están mis resultados (100 va).

Method       |     OSX        |  Linux Ubuntu
ogr buffer   | 0:00:01.434389 | 0:01:08.384309

sub string   | 0:00:19.714432 | 0:00:10.048649

sorted set   | 0:00:01.239999 | 0:00:00.600773


Specs Mac OSX
Processor 4x2.5 GHz  
Memory 8 GB 1600 MHz

Specs Dell Linux Ubuntu
Processor 8x3.4GHz
Memory 7.8 GB

Si alguien puede explicar el por qué de estas diferencias se producen, no lo dudes.

Answer 1

Evitar la consulta espacial

Desde que anotó el almacenamiento en búfer es computacionalmente costoso, y puede ser la celebración de vuelta, considerar este enfoque: Inicio de bucle a través de cada punto y completar tu lat long point a un lugar decimal dentro de su buffer (es decir, si la lat/long es 12.3456789/12.3456789, a continuación, obtener todos los puntos que comienzan con un lat/long de 12.34567/12.34567 o 12.34568/12.34568 o 12.34567/12.34568 o 12.34568/12.34567). El uso de una tabla hash para ello. Tome este subconjunto de puntos, obtener todas las distancias a su punto de entrada, y el punto con distancia mínima es la que usted desea. La creación de una metodología de búsqueda hará que este muy eficiente.

Esto evita tener que hacer costosas consultas espaciales y consulta de instalación de filtro de 800.000 veces. Usted sólo estaría haciendo de cadena doble, comparaciones y cálculos de distancia en este método. El único inconveniente que pude ver de este método es que cada uno de redondeo decimal es un orden de magnitud por encima de la última, así que si tu consulta espacial de no retorno puntos, redondeando hacia abajo de nuevo, puede volver muchas más puntos que usted necesita, que se ralentizará un poco. Sin embargo, usted tiene por lo menos dos órdenes de magnitud en el orignal BUFFER_DISTANCES, así que creo que este método puede ser suficiente para sus fines, y sin duda sería más rápido que el método que han de ir ahora mismo.

La tabla hash:

He aquí una forma más concisa y mejor explicación introductoria para las tablas hash.

El concepto es así: Usted quiere buscar la definición de la palabra "Sí" en el diccionario. No mirar a través de cada palabra en un diccionario a partir de Un a encontrar las palabras que empiezan con Vosotros, ¿correcto? Saltar directamente a la Y sección primera, luego buscar la página que dice Ye-Yo y, a continuación, escanea todas las palabras en la página para obtener la definición de Sí.

La búsqueda de una metodología para el bucle a través de todos los lat/long puntos podrían ser implementadas en esta misma forma. Usted tendría que buscar todos los lat/longs que comienzan con 12.3 en una serie de, digamos, por ejemplo, de 0 a 99. Entonces quieres buscar en los 10 valores para 12.34, y así sucesivamente. Si se programa correctamente, usted puede devolver un "cubo" con todos los puntos dentro de su buffer, sin tener que hacer una sola consulta espacial o de cadena doble en comparación!

Por último, cabe señalar que, si almacena su tabla indizada en un RDBMS, no puede ser de optimización para esto ya. Si la lat/long valores son dobles y una simple ENTRE el SQL de la consulta, es probable que tenga su función de búsqueda optimizado ya que hacer esto (si su consulta está escrito correctamente).

Answer 2

Has mirado en el uso de la KDTree algoritmo de scipy.espacial?

http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html

Así que yo sé que usted está usando ogr, pero aquí es un ejemplo que me escribió el uso de la scipy cKDTree para hacer un punto radio de búsqueda:

import arcpy
import scipy.spatial

# convert feature class to numpy array using "SHAPE@XY" token for geometry
# data from: http://nces.ed.gov/surveys/sdds/index.aspx, spatial data download inside map viewer data tools
feature_class = r'C:\gis_working_directory\school_clustering.gdb\CCD_10_11_WM' 
fields = ['ncessch', 'schnam', "SHAPE@JSON", "SHAPE@XY", "SHAPE@"]
where_clause = "{} = '{}'".format('mcity', 'NEW ORLEANS')
bourbon_street = (-10026195.958134, 3498018.476606)
numpy_array = arcpy.da.FeatureClassToNumPyArray(feature_class, ['ncessch', 'schnam', "SHAPE@XY"], where_clause)

# create kdtree using scipy
tree = scipy.spatial.cKDTree(numpy_array["SHAPE@XY"], leafsize=10)
radius = 1000 #meters
distances, indices = tree.query([bourbon_street], len(numpy_array), distance_upper_bound=radius)
distance_lookup = dict(zip(indices[0], distances[0]))
print [(numpy_array['schnam'][i], int(distance_lookup[i])) for i in indices[0] if i != len(numpy_array)]