La forma más rápida de unir espacialmente un CSV de punto con un archivo de forma poligonal

Si "más rápido" incluye la cantidad de su tiempo invertido, la solución dependerá del software con el que se sienta cómodo y que pueda usar de manera expedita. En consecuencia, los siguientes comentarios se centran en ideas para lograr los tiempos más rápidos posibles de computación .

Si usa un programa enlatado, casi con seguridad lo mejor que puede hacer es preprocesar los polígonos para configurar una estructura de datos de punto en polígono, como un árbol KD o quadtree, cuyo rendimiento será típicamente O ( log (V) * (N + V)) donde V es el número total de vértices en los polígonos y N es el número de puntos, porque la estructura de datos tomará al menos O (log (V) * V) esfuerzo para crear y luego tendrá que sondearse para cada punto a un costo por punto O (log (V)).

Puede hacerlo mucho mejor si primero coloca los polígonos en cuadrícula, aprovechando la suposición de que no hay superposiciones. Cada celda de la cuadrícula se encuentra completamente en el interior de un polígono (incluido el interior del "polígono universal"), en cuyo caso se etiqueta la celda con la identificación del polígono, o bien contiene uno o más bordes de polígono. El costo de esta rasterización, igual al número de celdas de la cuadrícula referenciadas al rasterizar todos los bordes, es O (V / c) donde c es el tamaño de una celda, pero la constante implícita en la notación big-O es pequeña. / p>

(Una belleza de este enfoque es que puede explotar rutinas gráficas estándar. Por ejemplo, si tiene un sistema que (a) dibujará los polígonos en una pantalla virtual usando (b) un color distinto para cada polígono y ( c) te permite leer el color de cualquier píxel al que te interese dirigirse, lo has hecho.)

Con esta cuadrícula en su lugar, haga una evaluación previa de los puntos calculando la celda que contiene cada punto (una operación O (1) que requiere solo unos pocos relojes). A menos que los puntos se agrupen alrededor de los límites del polígono, esto normalmente dejará solo unos puntos O (c) con resultados ambiguos. El costo total de construir la cuadrícula y la preselección es, por lo tanto, O (V / c + 1 / c ^ 2) + O (N). Debe utilizar algún otro método (como cualquiera de los recomendados hasta ahora) para procesar los puntos restantes (es decir, aquellos que están cerca de los límites de polígonos), a un costo de O (log (V) * N * c) .

A medida que c se hace más pequeño, cada vez menos puntos estarán en la misma celda de la cuadrícula con un borde y, por lo tanto, cada vez serán menos los que requerirán el procesamiento posterior de O (log (V)). Actuando contra esto es la necesidad de almacenar las celdas de la cuadrícula O (1 / c ^ 2) y de pasar O (V / c + 1 / c ^ 2) rasterizando los polígonos. Por lo tanto, habrá un tamaño de cuadrícula óptimo c. Usándolo, el costo computacional total es O (log (V) * N) pero la constante implícita es típicamente manera más pequeña que usar los procedimientos enlatados, debido a la velocidad O (N) de la pre- proyección.

Hace 20 años probé este enfoque (utilizando puntos uniformemente espaciados en Inglaterra y en alta mar y explotando una red relativamente cruda de alrededor de 400K celdas ofrecidas por los buffers de video de la época) y obtuve dos órdenes de magnitud de aceleración en comparación con los mejores publicados Algoritmo que pude encontrar. Incluso cuando los polígonos son pequeños y simples (como los triángulos), está virtualmente asegurado de un orden de magnitud de aceleración.

En mi experiencia, el cálculo fue tan rápido que toda la operación estuvo limitada por las velocidades de E / S de datos, no por la CPU. Anticipando que la E / S podría ser el cuello de botella, lograría los resultados más rápidos al almacenar los puntos en un formato lo más comprimido posible para minimizar los tiempos de lectura de datos. También piense cómo deberían almacenarse los resultados, para poder limitar las escrituras en el disco.

Por mi parte, probablemente cargaría datos CSV en un archivo shp y luego escribiría un script en python usando shapefile y shapely para obtener el ID del polígono que lo contiene y actualizar el valor del campo.

No sé si los geotools y JTS son más rápidos que shapefile / shapely ... ¡No tienen tiempo para probarlo!

editar : Por cierto, la conversión de csv al formato de shapefile probablemente no sea necesaria, ya que los valores podrían formatearse fácilmente con objetos espaciales de su archivo de forma poligonal.

Terminé convirtiendo los polígonos en un ráster y muestreando los puntos en las posiciones de los puntos. Ya que mis polígonos no se superponían y la alta precisión no era necesaria (los polígonos representaban clases de uso de la tierra y sus fronteras se consideraban bastante inciertas de todos modos) esta era la solución más eficiente en tiempo que podía encontrar.

Escribiría rápidamente un pequeño programa java basado en lector de shapefile de geotools y la operación contains de JTS . No sé qué tan rápido puede ser ...

Use Spatialite .

Descarga la GUI. Puede abrir tanto el Shapefile como el CSV como tablas virtuales. Esto significa que en realidad no los importa a la base de datos, sino que aparecen como tablas y puede unirse y consultarlos de la manera que desee.

Puede hacerlo con bastante rapidez utilizando OGR en C / C ++ / Python (Python debería ser el más lento de los 3). Recorra todos los polígonos y establezca un filtro en los puntos, recorra los puntos filtrados y sabrá que cada uno de los puntos que recorra pertenecerá al polígono actual. Aquí hay un código de muestra en python que usa OGR que recorrerá los polígonos y los puntos de filtro en consecuencia. El código C / C ++ se verá bastante similar a este, y me imagino que obtendrás un aumento de velocidad significativo en comparación con Python. Tendrá que agregar algunas líneas de código para actualizar el CSV a medida que avanza:

from osgeo import ogr
from osgeo.gdalconst import *

inPolyDS = ogr.Open("winnipeg.shp", GA_ReadOnly)
inPolyLayer = inPolyDS.GetLayer(0)
inPointDS = ogr.Open("busstops.vrt", GA_ReadOnly)   
inPointLayer = inPointDS.GetLayerByName("busstops")

inPolyFeat = inPolyLayer.GetNextFeature()
while inPolyFeat is not None:
  inPtFeat = inPointLayer.GetNextFeature()
  while inPtFeat is not None:
    ptGeom = inPtFeat.GetGeometryRef()
    # Do work here...

    inPtFeat = inPointLayer.GetNextFeature()

  inPolyFeat = inPolyLayer.GetNextFeature()

Archivo VRT (busstops.vrt):

<OGRVRTDataSource>
  <OGRVRTLayer name="busstops">
    <SrcDataSource>busstops.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" reportSrcColumn="FALSE" />
  </OGRVRTLayer>
</OGRVRTDataSource>

Archivo CSV (busstops.csv):

FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey

Archivo CSVT (busstops.csvt, OGR lo necesita para identificar tipos de columnas, de lo contrario no realizará el filtro espacial):

Integer,Real,Real,String

podría probar csv2shp csv2shp

¿Sientes curiosidad por saber en qué industria se encuentran los mil millones de puntos de CSV?