Buscando la solución más rápida para el análisis Point in Polygon de 200 millones de puntos [cerrado]

Tengo un CSV que contiene 200 millones de observaciones con el siguiente formato:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"

Para cada conjunto de coordenadas (x1 / y1 y x2 / y2), quiero asignar el Área del Censo de EE. UU. O el Bloque del Censo en el que se encuentra (descargué el archivo de forma TIGER del tracto del Censo aquí: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Entonces, necesito hacer una operación de punto en el polígono dos veces para cada observación. Es importante que las coincidencias sean muy precisas.

¿Cuál es la forma más rápida de hacer esto, incluido el tiempo para aprender el software? Tengo acceso a una computadora con 48 GB de memoria, en caso de que pueda ser una restricción relevante.

Varios hilos recomiendan usar PostGIS o Spatialite (Spatialite parece más fácil de usar, pero ¿es tan eficiente como PostGIS?). Si esas son las mejores opciones, ¿es imprescindible completar un índice espacial (RTree?)? Si es así, ¿cómo se hace (p. Ej., Usando el archivo de forma del tracto censal)? Estaría extremadamente agradecido por cualquier recomendación que incluya código de ejemplo (o un puntero al código de ejemplo).

Mi primer intento (antes de encontrar este sitio) consistió en usar ArcGIS para hacer una unión espacial (solo x1 / y1) de submuestra de datos (100,000 puntos) en el Bloque del Censo de EE. UU. Eso tomó más de 5 horas antes de que matara el proceso. Espero una solución que pueda implementarse en todo el conjunto de datos en menos de 40 horas de tiempo de computación.

Disculpas por hacer una pregunta que se ha hecho antes: he leído las respuestas y me pregunto cómo implementar las recomendaciones. Nunca he usado SQL, Python, C, y solo he usado ArcGIS una vez antes: soy un principiante completo.

ST_DWithin fue más rápido en mi prueba que ST_Intersects. Eso es sorprendente, especialmente porque se supone que el algoritmo de geometría preparado se aplicará en casos como este. Creo que existe la posibilidad de que esto sea mucho más rápido de lo que mostré aquí.

Hice algunas pruebas más y dos cosas casi duplicaron la velocidad. Primero, probé en una computadora más nueva, pero todavía era una computadora portátil bastante común, tal vez excepto en los discos ssd SATA3.

Luego, la consulta a continuación tomó 18 segundos en lugar de 62 segundos en la computadora portátil anterior. Luego descubrí que estaba completamente equivocado antes cuando escribí que el índice en la tabla de puntos no era necesario. Con ese índice en su lugar, ST_Intersects se comportó como se esperaba y las cosas se volvieron muy rápidas. Aumenté el número de puntos en la tabla de puntos a 1 millón de puntos y la consulta:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct , t WHERE ST_Intersects(imported_ct.geom , t.geom);

corre en 72 segundos. Como hay 1249 polígonos, las pruebas 1249000000 se realizan en 72 segundos. Eso hace aproximadamente 17000000 pruebas por segundo. O probar casi 14000 puntos contra todos los polígonos por segundo.

De esta prueba, sus 400000000 puntos para probar deberían tomar aproximadamente 8 horas sin ningún problema con la distribución de la carga a varios núcleos. PostGIS nunca se detiene para impresionarme :-)

Primero, para visualizar el resultado, puede agregar la geometría del punto a la tabla resultante, abrirla en QGIS, por ejemplo, y aplicarle un estilo con valores únicos en el campo import_ct.

Segundo, sí, también puede obtener los puntos que caen fuera de cualquier polígono usando una unión derecha (o izquierda) como esta:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct right join t ON ST_Intersects(imported_ct.the_geom , t.geom);

Hice algunas pruebas para verificar si parece posible PostGIS.

Primero algo que no entiendo. Tienes dos puntos por fila. ¿Están siempre ambos puntos en el mismo polígono? Entonces es suficiente hacer los cálculos en uno de los puntos. Si pueden estar en dos polígonos diferentes, necesitará una forma de conectar una fila de puntos a dos polígonos.

Según las pruebas, parece factible, pero es posible que necesite alguna solución creativa para distribuir la carga en más de un núcleo de CPU.

Probé en una computadora portátil de 4 años con CPU de doble núcleo centrino (creo que aproximadamente 2.2 GHz), 2 GB de RAM. Si tienes 48 BG RAM, supongo que también tienes mucha más potencia de CPU.

Lo que hice fue crear una tabla de puntos aleatorios con 100000 puntos como este:

CREATE TABLE t AS
WITH r AS
(SELECT ST_Extent(the_geom)::geometry ext FROM imported_ct)
SELECT ST_Point(x,y) AS geom FROM 
(SELECT GENERATE_SERIES(1,100000)) s,
(SELECT ST_Xmin(ext)+(random()*(ST_Xmax(ext)-ST_Xmin(ext))) x, ST_Ymin(ext)+(random()*(ST_Ymax(ext)-ST_Ymin(ext))) y FROM r
) f;

Luego agregando un gid como:

ALTER TABLE t ADD COLUMN GID SERIAL;

Luego corriendo:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE ST_Dwithin(imported_ct.the_geom , t.geom,0);

toma aproximadamente 62 segundos (compárelo con el resultado de ArcGIS con la misma cantidad de puntos). El resultado es una tabla que conecta los puntos en mi tabla t con el gid en la tabla con el censo.

Con esa velocidad, harás 200 puntos de molino en aproximadamente 34 horas. Entonces, si es suficiente con marcar uno de los puntos, mi vieja computadora portátil puede hacerlo con un núcleo.

Pero si necesita verificar ambos puntos, podría ser más difícil.

Luego, puede distribuir manualmente la carga a más de un núcleo iniciando varias sesiones en la base de datos y ejecutando consultas diferentes.

En mi ejemplo con 50000 puntos y dos núcleos de CPU probé:

CREATE TABLE t1 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid >50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

en una sesión db al mismo tiempo que se ejecuta:

CREATE TABLE t2 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid <=50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

en otra sesión de db.

Eso tomó alrededor de 36 segundos, por lo que es un poco más lento que el primer ejemplo, probablemente dependiendo de la escritura del disco al mismo tiempo. Pero dado que los núcleos bith funcionan al mismo tiempo, no me tomó más de 36 segundos de mi tiempo.

Para unir las tablas t1 y t2 se intentó:

CREATE TABLE t3 AS 
SELECT * FROM t1
UNION ALL
SELECT * FROM t2;

usando aproximadamente medio segundo.

Entonces, con un hardware más reciente y distribuyendo la carga en muchos núcleos, esto debería ser absolutamente posible incluso si el mundo real será más lento que el caso de prueba.

Vale la pena señalar que el ejemplo es de Linux (Ubuntu). Usar Windows será otra historia. Pero tengo todas las otras aplicaciones diarias ejecutándose, por lo que la computadora portátil está bastante cargada desde antes. Entonces, eso podría simular bastante bien el caso de Windows, sin abrir nada más que pgadmin.

Probablemente la forma más fácil es con PostGIS. Hay algunos tutoriales en Internet sobre cómo importar datos de puntos csv / txt en PostGIS. Enlace1

No estoy seguro sobre el rendimiento de las búsquedas de punto en el polígono en PostGIS; debería ser más rápido que ArcGIS. El índice espacial GIST que utiliza PostGIS es bastante rápido. Enlace2 Enlace3

También puede probar el índice geoespacial MongoDB . Pero esto requiere poco más tiempo para comenzar. Creo que MongoDB podría ser realmente rápido. No lo he probado con búsquedas de punto en el polígono, así que no puedo estar seguro.