Procurando a solução mais rápida para análise Point in Polygon de 200 milhões de pontos [fechado]

Eu tenho um CSV contendo 200 milhões de observações com o seguinte 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 e x2 / y2), desejo atribuir o setor censitário dos EUA ou o bloco censitário em que ele se enquadra (baixei o arquivo de shapefile TIGER do setor censitário aqui: ftp://ftp2.census.gov/ geo / tigre / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Então, eu preciso fazer uma operação point-in-polygon duas vezes para cada observação. É importante que as correspondências sejam muito precisas.

Qual é a maneira mais rápida de fazer isso, incluindo tempo para aprender o software? Eu tenho acesso a um computador com 48 GB de memória - caso isso possa ser uma restrição relevante.

Vários threads recomendam o uso de PostGIS ou Spatialite (o Spatialite parece mais fácil de usar - mas é tão eficiente quanto o PostGIS?). Se essas são as melhores opções, é imperativo preencher um índice espacial (RTree?)? Se sim, como fazer isso (por exemplo, usando o Shapefile do Censo)? Ficaria extremamente grato por quaisquer recomendações que incluam código de exemplo (ou um ponteiro para código de exemplo).

Minha primeira tentativa (antes de encontrar este site) consistiu em usar o ArcGIS para fazer uma junção espacial (apenas x1 / y1) da subamostra dos dados (100.000 pontos) no US Census Block. Isso levou mais de 5 horas antes de eu matar o processo. Espero uma solução que possa ser implementada em todo o conjunto de dados em menos de 40 horas de tempo de computação.

Desculpas por fazer uma pergunta que já foi feita antes - li as respostas e fiquei me perguntando como implementar as recomendações. Eu nunca usei SQL, Python, C e só usei o ArcGIS uma vez antes - sou um iniciante.

ST_DWithin foi mais rápido no meu teste do que ST_Intersects. Isso é surpreendente, especialmente porque o algoritmo de geometria preparado deve funcionar em casos como este. Acho que há uma chance de que isso seja muito mais rápido do que mostrei aqui.

Fiz mais alguns testes e duas coisas quase dobraram a velocidade. Primeiro, tentei em um computador mais novo, mas ainda assim um laptop bastante comum, talvez exceto nos discos SSD SATA3.

Em seguida, a consulta abaixo levou 18 segundos em vez de 62 segundos no laptop antigo. Em seguida, descobri que estava totalmente errado antes, quando escrevi que o índice na tabela de pontos não era necessário. Com esse índice em vigor, ST_Intersects se comportou como esperado e as coisas se tornaram muito rápidas. Aumentei o número de pontos na tabela de pontos para 1 milhão de pontos e a 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);

é executado em 72 segundos. Como existem 1249 polígonos, 1249000000 testes são feitos em 72 segundos. Isso faz cerca de 17000000 testes por segundo. Ou testando quase 14000 pontos contra todos os polígonos por segundo.

Nesse teste, seus 400000000 pontos a serem testados devem levar cerca de 8 horas sem problemas para distribuir a carga em vários núcleos. O PostGIS nunca para para me impressionar :-)

Primeiro, para visualizar o resultado, você pode adicionar a geometria do ponto à tabela resultante, abra-a no QGIS, por exemplo, e estilize-a com valores exclusivos no campo imports_ct.

Segundo, sim, você também pode obter os pontos que ficam fora de qualquer polígono usando uma junção direita (ou esquerda) assim:

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);

Fiz alguns testes para verificar se parece possível o PostGIS.

Primeiro algo que eu não entendo. Você tem dois pontos por linha. Sempre os dois pontos estão no mesmo polígono? Então é suficiente fazer os cálculos em um dos pontos. Se eles puderem estar em dois polígonos diferentes, você precisará de uma maneira de conectar uma linha de ponto a dois polígonos.

A partir dos testes, parece factível, mas você pode precisar de alguma solução criativa para distribuir a carga por mais de um núcleo de CPU.

Testei em um laptop de 4 anos com CPU dual core centrino (cerca de 2,2 GHz, eu acho), 2 GB de RAM. Se você tem 48 BG de RAM, acho que também tem muito mais poder de CPU.

O que eu fiz foi criar uma tabela de pontos aleatórios com 100000 pontos 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;

Em seguida, adicione um gid como:

ALTER TABLE t ADD COLUMN GID SERIAL;

Em seguida, executando:

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);

leva cerca de 62 segundos (compare com o resultado do ArcGIS com a mesma quantidade de pontos). O resultado é uma tabela que liga os pontos da minha tabela t com o gid da tabela com o setor censitário.

Com essa velocidade, você fará 200 pontos de moinho em cerca de 34 horas. Portanto, se basta verificar um dos pontos, meu laptop antigo pode fazê-lo com um único núcleo.

Mas se você precisar verificar os dois pontos, pode ser mais difícil.

Então, você pode distribuir manualmente a carga para mais de um núcleo, iniciando várias sessões no banco de dados e executando consultas diferentes.

No meu exemplo, com 50000 pontos e dois núcleos de CPU, tentei:

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);

em uma sessão de banco de dados ao mesmo tempo que em execução:

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);

em outra sessão de banco de dados.

Isso levou cerca de 36 segundos, por isso é um pouco mais lento que o primeiro exemplo, provavelmente dependendo da gravação do disco ao mesmo tempo. Mas como os núcleos bith estão funcionando ao mesmo tempo, não demorou mais de 36 segundos do meu tempo.

Para unir as tabelas t1 e t2 a tentei:

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

usando cerca de meio segundo.

Portanto, com um hardware mais atualizado e distribuindo a carga por muitos núcleos, isso deve ser absolutamente possível, mesmo que o mundo real seja mais lento que o caso de teste.

Pode valer a pena notar que o exemplo é do Linux (Ubuntu). Usar o Windows será outra história. Mas eu tenho todos os outros aplicativos diários em execução, portanto o laptop está bastante carregado de antes. Portanto, isso pode simular o caso do Windows muito bem, talvez, sem abrir nada além do pgadmin.

Provavelmente, a maneira mais fácil é com o PostGIS. Existem alguns tutoriais na Internet sobre como importar dados de ponto csv / txt para o PostGIS. Link1

Não tenho certeza sobre o desempenho de pesquisas point-in-polygon no PostGIS; deve ser mais rápido que o ArcGIS. O índice espacial GIST usado pelo PostGIS é bastante rápido. Link2 Link3

Você também pode testar o índice geoespacial do MongoDB . Mas isso requer um pouco mais de tempo para começar. Acredito que o MongoDB possa ser muito rápido. Não testei com pesquisas point-in-polygon, por isso não posso ter certeza.