Le moyen le plus rapide de joindre spatialement un point CSV avec un polygone Shapefile

J'ai un fichier CSV de 1 milliard de points et un fichier de formes avec environ 5 000 polygones. Quel serait le moyen le plus rapide de joindre spatialement des points et des polygones? Pour chaque point, j'ai besoin d'obtenir l'id du polygone contenant. (Les polygones ne se chevauchent pas.)

Habituellement, je charge les deux ensembles de données dans PostGIS. Existe-t-il un moyen plus rapide de faire le travail?

Je recherche une solution open source.

Si « le plus rapide » comprend la quantité de votre temps passé, la solution dépendra de ce que le logiciel que vous êtes à l' aise avec et peut utiliser rapidement. Les remarques suivantes se concentrent par conséquent sur les idées pour atteindre les temps de calcul les plus rapides possibles .

Si vous utilisez un programme préenregistré, le mieux que vous puissiez faire est presque certainement de prétraiter les polygones pour configurer une structure de données point par polygone, comme un arbre KD ou un arbre quadruple, dont les performances seront généralement O (log (V ) * (N + V)) où V est le nombre total de sommets dans les polygones et N est le nombre de points, car la structure des données prendra au moins O (log (V) * V) effort pour créer et ensuite doivent être sondés pour chaque point à un coût par point O (log (V)).

Vous pouvez faire beaucoup mieux en quadrillant d'abord les polygones, en exploitant l'hypothèse d'aucun chevauchement. Chaque cellule de la grille se trouve soit entièrement dans un intérieur de polygone (y compris l'intérieur du "polygone universel"), auquel cas étiquetez la cellule avec l'identifiant du polygone, ou bien elle contient un ou plusieurs bords de polygone. Le coût de cette pixellisation, égal au nombre de cellules de grille référencées lors de la pixellisation de toutes les arêtes, est O (V / c) où c est la taille d'une cellule, mais la constante implicite dans la notation big-O est petite.

(Une beauté de cette approche est que vous pouvez exploiter des routines graphiques standard. Par exemple, si vous avez un système qui (a) dessinera les polygones sur un écran virtuel en utilisant (b) une couleur distincte pour chaque polygone et (c) permet vous devez lire la couleur de tout pixel que vous souhaitez traiter, vous l'avez fait.)

Avec cette grille en place, présélectionnez les points en calculant la cellule contenant chaque point (une opération O (1) ne nécessitant que quelques horloges). À moins que les points ne soient regroupés autour des limites du polygone, cela ne laissera généralement qu'environ O (c) points avec des résultats ambigus. Le coût total de construction de la grille et de présélection est donc O (V / c + 1 / c ^ 2) + O (N). Vous devez utiliser une autre méthode (telle que celles recommandées jusqu'à présent) pour traiter les points restants (c'est-à-dire ceux qui sont proches des limites des polygones), au coût de O (log (V) * N * c) .

Comme c devient plus petit, de moins en moins de points seront dans la même cellule de grille avec un bord et donc de moins en moins nécessiteront le traitement O (log (V)) ultérieur. Agir contre cela est la nécessité de stocker O (1 / c ^ 2) cellules de la grille et de passer O (V / c + 1 / c ^ 2) temps à tramer les polygones. Il y aura donc une taille de grille optimale c. En l'utilisant, le coût de calcul total est O (log (V) * N) mais la constante implicite est généralement beaucoup plus petite que l'utilisation des procédures prédéfinies, en raison de la vitesse O (N) de la présélection.

Il y a 20 ans, j'ai testé cette approche (en utilisant des points uniformément espacés dans toute l'Angleterre et au large des côtes et en exploitant une grille relativement brute d'environ 400 000 cellules offerte par les tampons vidéo de l'époque) et j'ai obtenu une accélération de deux ordres de grandeur par rapport au meilleur algorithme publié que j'ai pu trouver. Même lorsque les polygones sont petits et simples (comme les triangles), vous êtes pratiquement assuré d'une accélération d'un ordre de grandeur.

D'après mon expérience, le calcul était si rapide que toute l'opération était limitée par les vitesses d'E / S de données, et non par le CPU. En anticipant que les E / S pourraient être le goulot d'étranglement, vous obtiendriez les résultats les plus rapides en stockant les points dans un format aussi compressé que possible pour minimiser les temps de lecture des données. Réfléchissez également à la façon dont les résultats doivent être stockés, afin de limiter les écritures sur disque.

Pour ma part, je chargerais probablement des données CSV dans un fichier shp , puis j'écrirais un script python en utilisant shapefile et shapely pour obtenir l'ID du polygone contenant et mettre à jour la valeur du champ.

Je ne sais pas si les géotools et JTS sont plus rapides que shapefile / shapely ... Pas le temps de le tester!

edit : Soit dit en passant, la conversion csv au format de fichier de formes n'est probablement pas nécessaire, car les valeurs pourraient facilement être formatées pour être testées avec des objets spatiaux de votre fichier de formes de polygones.

J'ai fini par convertir les polygones en raster et les échantillonner aux positions des points. Étant donné que mes polygones ne se chevauchaient pas et qu'une grande précision n'était pas nécessaire (les polygones représentaient des classes d'utilisation des terres et leurs frontières étaient de toute façon considérées comme plutôt incertaines), c'était la solution la plus efficace en termes de temps que j'ai pu trouver.

J'écrirais rapidement un petit programme java basé sur le lecteur de fichiers de formes de geotools et l'opération contient de JTS . Je ne sais pas à quelle vitesse ça peut être ...

Utilisez Spatialite .

Téléchargez l'interface graphique. Vous pouvez ouvrir le Shapefile et le CSV en tant que tables virtuelles. Cela signifie que vous ne les importez pas réellement dans la base de données, mais ils apparaissent sous forme de tableaux et vous pouvez rapidement les joindre et les interroger comme vous le souhaitez.

Vous pouvez le faire assez rapidement en utilisant OGR en C / C ++ / Python (Python devrait être le plus lent des 3). Parcourez tous les polygones et définissez un filtre sur les points, parcourez les points filtrés et vous saurez que chacun des points que vous parcourez appartiendra au polygone actuel. Voici un exemple de code en python utilisant OGR qui parcourra les polygones et filtrera les points en conséquence. Le code C / C ++ ressemblera assez à cela, et j'imagine que vous obtiendrez une augmentation de vitesse significative par rapport à python. Vous devrez ajouter quelques lignes de code pour mettre à jour le CSV au fur et à mesure:

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

Fichier 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>

Fichier CSV (busstops.csv):

FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey

Fichier CSVT (busstops.csvt, OGR en a besoin pour identifier les types de colonnes, sinon il n'effectuera pas le filtre spatial):

Integer,Real,Real,String

pourrait essayer csv2shp csv2shp

curieux de savoir dans quelle industrie se situe le CSV en milliards de points?