La forma más rápida de unir espacialmente un CSV de puntos con un Shapefile de polígonos

Question

Tengo un fichero CSV de 1.000 millones de puntos y un shapefile con unos 5.000 polígonos. ¿Cuál sería la forma más rápida de unir espacialmente puntos y polígonos? Para cada punto, necesito obtener el id del polígono que lo contiene. (Los polígonos no se solapan.)

Normalmente, cargaría ambos conjuntos de datos en PostGIS. ¿Hay alguna manera más rápida de hacer el trabajo?

Busco una solución de código abierto.

Answer 1

Puedes hacerlo bastante rápido usando OGR en C/C++/Python (Python debería ser el más lento de los 3). Haz un bucle a través de todos los polígonos y establece un filtro en los puntos, haz un bucle a través de los puntos filtrados y sabrás que cada uno de los puntos por los que haces un bucle pertenecerá al polígono actual. Aquí hay un ejemplo de código en python utilizando OGR que hará un bucle a través de los polígonos y filtrará los puntos en consecuencia. El código C/C++ se parecerá bastante a este, y me imagino que obtendrás un aumento significativo de velocidad frente a python. Usted tendrá que añadir unas pocas 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 los tipos de columna, de lo contrario no realizará el filtro espacial):

Integer,Real,Real,String

Answer 2

Podría probar csv2shp csv2shp

¿tiene curiosidad por saber a qué sector pertenece el CSV de mil millones de puntos?