À la recherche de la solution la plus rapide pour une analyse ponctuelle de 200 millions de points [fermée]

J'ai un fichier CSV contenant 200 millions d'observations au format suivant:

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"

Pour chaque ensemble de coordonnées (x1 / y1 et x2 / y2), je souhaite attribuer le secteur de recensement américain ou le bloc de recensement auquel il appartient (j'ai téléchargé le fichier de formes TIGER du secteur de recensement ici: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Je dois donc effectuer deux fois une opération de point dans un polygone pour chaque observation. Il est important que les matchs soient très précis.

Quel est le moyen le plus rapide de le faire, y compris le temps nécessaire pour apprendre le logiciel? J'ai accès à un ordinateur avec 48 Go de mémoire - au cas où cela pourrait être une contrainte pertinente.

Plusieurs discussions recommandent d'utiliser PostGIS ou Spatialite (Spatialite semble plus facile à utiliser - mais est-il aussi efficace que PostGIS?). Si ce sont les meilleures options, est-il impératif de renseigner un index spatial (RTree?)? Si oui, comment fait-on cela (par exemple, en utilisant le fichier de forme de secteur de recensement)? Je serais extrêmement reconnaissant pour toutes les recommandations qui incluent un exemple de code (ou un pointeur sur un exemple de code).

Ma première tentative (avant de trouver ce site) consistait à utiliser ArcGIS pour effectuer une jointure spatiale (x1 / y1 uniquement) du sous-échantillon des données (100 000 points) sur US Census Block. Cela a pris 5 heures avant que je tue le processus. J'espère une solution pouvant être mise en œuvre sur l'ensemble des données en moins de 40 heures de temps de calcul.

Toutes mes excuses pour avoir posé une question qui a déjà été posée - j'ai lu les réponses et je me suis demandé comment mettre en œuvre les recommandations. Je n'ai jamais utilisé SQL, Python, C et je n'ai utilisé ArcGIS qu'une seule fois auparavant - je suis un débutant complet.

ST_DWinin était plus rapide dans mon test que ST_Intersects. C'est surprenant, d'autant plus que l'algorithme de géométrie préparée est censé intervenir dans de tels cas. Je pense qu'il y a une chance que ce soit beaucoup plus rapide que ce que j'ai montré ici.

J'ai fait quelques tests supplémentaires et deux choses ont presque doublé la vitesse. Premièrement, j'ai essayé sur un ordinateur plus récent, mais toujours un ordinateur portable assez ordinaire, peut-être à l'exception des disques SATA3 ssd.

Ensuite, la requête ci-dessous a pris 18 secondes au lieu de 62 secondes sur l'ancien ordinateur portable. Ensuite, j'ai découvert que j'avais totalement tort avant d'écrire que l'index de la table de points n'était pas nécessaire. Avec cet index en place, ST_Intersects s'est comporté comme prévu et les choses sont devenues très rapides. J'ai augmenté le nombre de points dans la table de points à 1 million de points et la requête:

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

fonctionne en 72 secondes. Comme il y a 1249 polygones, 1249000000 tests sont effectués en 72 secondes. Cela fait environ 17000000 tests par seconde. Ou tester près de 14 000 points contre tous les polygones par seconde.

À partir de ce test, les 400000000 points à tester devraient prendre environ 8 heures sans aucun problème pour répartir la charge sur plusieurs cœurs. PostGIS ne cesse jamais de m'impressionner :-)

Tout d'abord, pour visualiser le résultat, vous pouvez ajouter la géométrie de point à la table résultante, l'ouvrir par exemple dans QGIS et la styliser avec des valeurs uniques sur le champ imports_ct.

Deuxièmement, vous pouvez également obtenir les points situés en dehors de tout polygone en utilisant une jointure à droite (ou à gauche) comme ceci:

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

J'ai fait des tests pour vérifier si cela semble possible PostGIS.

D'abord quelque chose que je ne comprends pas. Vous avez deux points par ligne. Les deux points sont-ils toujours dans le même polygone? Ensuite, il suffit de faire les calculs sur l’un des points. S'ils peuvent être dans deux polygones différents, vous aurez besoin d'un moyen de connecter une ligne de point à deux polygones.

D'après les tests, cela semble faisable, mais vous aurez peut-être besoin d'une solution créative pour répartir la charge sur plusieurs cpu-core.

J'ai testé sur un ordinateur portable de 4 ans avec processeur double cœur centrino (environ 2,2 GHz, je pense), 2 Go de RAM. Si vous avez 48 Go de RAM, je suppose que vous avez aussi beaucoup plus de puissance CPU.

Ce que j'ai fait était de créer une table de points aléatoires avec 100 000 points comme ceci:

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;

Puis ajouter un gid comme:

ALTER TABLE t ADD COLUMN GID SERIAL;

Puis en cours d'exécution:

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

prend environ 62 secondes (comparez votre résultat à ArcGIS avec le même nombre de points). Le résultat est un tableau reliant les points de mon tableau t au gid du tableau avec le secteur de recensement.

Avec cette vitesse, vous obtiendrez 200 points en 34 heures environ. Donc, s'il suffit de vérifier l'un des points, mon ancien ordinateur portable peut le faire avec un seul noyau.

Mais si vous devez vérifier les deux points, cela pourrait être plus difficile.

Ensuite, vous pouvez répartir manuellement la charge sur plusieurs cœurs en démarrant plusieurs sessions sur la base de données et en exécutant différentes requêtes.

Dans mon exemple avec 50000 points et deux cœurs de processeur, j'ai essayé:

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

sur une session de base de données en même temps que l'exécution:

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

sur une autre db-session.

Cela a pris environ 36 secondes, il est donc un peu plus lent que le premier exemple, probablement en fonction de l’écriture du disque en même temps. Mais comme les deux noyaux fonctionnent simultanément, cela ne m'a pas pris plus de 36 secondes de mon temps.

Pour joindre les tables t1 et t2 a essayé:

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

en utilisant environ une demi-seconde.

Ainsi, avec du matériel plus récent et une répartition de la charge sur de nombreux cœurs, cela devrait absolument être possible, même si le monde réel sera plus lent que le scénario de test.

Il convient de noter que l'exemple provient de Linux (Ubuntu). Utiliser Windows sera une autre histoire. Mais toutes les autres applications quotidiennes sont en cours d’exécution, de sorte que l’ordinateur portable est assez chargé d’avant. Cela pourrait donc très bien simuler le cas Windows, sans rien ouvrir d'autre que pgadmin.

Le moyen le plus simple est probablement d'utiliser PostGIS. Il existe des tutoriels sur Internet pour importer des données de points csv / txt dans PostGIS. Lien1

Je ne suis pas sûr de la performance des recherches de points dans polygones dans PostGIS; il devrait être plus rapide que ArcGIS. L'index spatial GIST utilisé par PostGIS est assez rapide. Link2 Link3

Vous pouvez également tester l'index géospatial MongoDB . Mais cela nécessite un peu plus de temps pour commencer. Je crois que MongoDB pourrait être très rapide. Je ne l'ai pas testé avec des recherches de point dans polygone, je ne peux donc pas en être sûr.