Buscando la solución más rápida para el análisis Punto en Polígono de 200 millones de puntos [cerrado]

34

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 Tracto del Censo de los EE. UU. o el Bloque del Censo en el que se encuentra (descargué el shapefile TIGER del Censo de los Estados Unidos aquí: ftp://ftp2.census.gov/geo/tiger/TIGER2011/TRACT/tl_2011_08_tract.zip ). Entonces, necesito hacer una operación de punto en 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 hacerlo, incluido el tiempo para aprender el software? Tengo acceso a una computadora con 48GB 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 imperativo llenar un Índice espacial (RTree)? Si es así, ¿cómo se hace (por ejemplo, utilizando el archivo de forma del censo)? Le agradecería enormemente cualquier recomendación que incluya un 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 realizar una unión espacial (solo x1 / y1) de la submuestra de los datos (100,000 puntos) en el Bloque de censos de EE. UU. Eso tomó más de 5 horas antes de que yo 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: leí 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.

    
pregunta meer 02.05.2012 - 06:54

2 respuestas

26

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

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

Luego, la siguiente consulta tomó 18 segundos en lugar de 62 segundos en la vieja computadora portátil. Luego descubrí que estaba totalmente 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é la cantidad 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);

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

A partir de esta prueba, sus 400000000 puntos a prueba deben tomar aproximadamente 8 horas sin ningún problema con la distribución de la carga a varios núcleos. PostGIS nunca deja de impresionarme :-)

Primero, para visualizar el resultado, puede agregar la geometría del punto a la tabla resultante, abrirlo en QGIS por ejemplo y decorarlo con valores únicos en el campo importado.

En segundo lugar, sí, también puede obtener los puntos que quedan fuera de cualquier polígono utilizando una combinación a la derecha (o a la 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. ¿Siempre están ambos puntos en el mismo polígono? Entonces basta con 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.

A partir de las pruebas, parece factible, pero es posible que necesite alguna solución creativa para distribuir la carga en más de un cpu-core.

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

Lo que hice fue crear una tabla de puntos al azar 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 se ejecuta:

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áralo con tu 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 distrito censal.

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

Pero si necesitas marcar ambos puntos, podría ser más difícil.

Luego puedes 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 cpu-cores que 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 db.

Eso llevó unos 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 están trabajando al mismo tiempo, no tardé más de 36 segundos de mi tiempo.

Para unir las tablas t1 y t2 a probado:

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

usando aproximadamente medio segundo.

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

Podría valer la pena señalar que el ejemplo es de Linux (Ubuntu). Usar Windows será otra historia. Pero tengo todas las demás aplicaciones diarias en ejecución, por lo que la computadora portátil está bastante cargada desde antes. Así que eso podría simular bastante bien el caso de Windows, sin abrir nada más que pgadmin.

    
respondido por el Nicklas Avén 02.05.2012 - 12:30
4

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

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

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

    
respondido por el Mario Miler 02.05.2012 - 11:01

Lea otras preguntas en las etiquetas