Comment trouver le polygone à l'intérieur duquel se situe un point?

8

J'ai une couche avec des entités polygonales. Chaque entité a des attributs et des valeurs. J'ai également une liste de coordonnées et je voudrais savoir dans quelle entité (ou polygone) se trouvent les coordonnées.

Quelqu'un pourrait-il me guider sur la façon de procéder? Existe-t-il une fonction dans l'API qui peut m'aider à atteindre mon objectif ou dois-je utiliser un algorithme de géométrie informatique pour le faire moi-même? Je sais comment faire ce dernier mais cela me ferait gagner du temps s'il y avait déjà une fonction intégrée.

Merci.

Shubham Goyal
la source

Réponses:

6

Peu de fonctionnalités

Ce que vous voulez probablement faire, c'est:

  1. Créez une liste de QgsPoint à partir de vos coordonnées
  2. Itérer sur toutes vos entités de couche (polygones)
  3. Boucle sur la liste des points (dans l'itération)
  4. Appelez feature.geometry (). Contains (point) pour vérifier si vous avez une correspondance

Bien sûr, vous pouvez désormais améliorer les performances si vous savez, par exemple, qu'un point ne peut être contenu que par un seul polygone, vous pouvez supprimer un point de la liste, une fois le polygone approprié trouvé.

Beaucoup de fonctionnalités (Utilisation de SpatialIndex)

Comme indiqué dans les commentaires, un indice spatial peut être utilisé pour accélérer considérablement le processus.

Les étapes ici seraient

  1. Créez une liste de QgsPoint à partir de vos coordonnées
  2. Créer un QgsSpatialIndex
  3. Itérer sur toutes vos entités de couche (polygones)
  4. Ajoutez les fonctionnalités à votre index avec insertFeature
  5. Répétez tous vos points
  6. Appelez index.intersects (QgsRectangle (point, point)) pour vérifier si vous avez une correspondance

Il y a aussi un exemple de code par NathanW

Matthias Kuhn
la source
Ah, je ne connaissais pas l'appel feature.geometry.contains (point). J'ai utilisé mathplotlib. pastebin.com/61LSeMWv Veuillez pardonner le désordre de mon code. Je suis pressé et je vais le nettoyer plus tard.
Shubham Goyal
Je n'ai pas implémenté votre solution et je ne peux donc pas dire avec certitude si cela fonctionne ou non. Mais je crois que cela devrait et je marque donc cela comme la bonne réponse :)
Shubham Goyal
2
Cela pourrait-il être accéléré à l'aide d'un QgsSpatialIndex?
Snorfalorpagus
1
Je recommande fortement d'utiliser un QgsSpatailIndex. Il y a un exemple sur nathanw.net/2013/01/04/…
Nathan W
6

Vous devez tout d'abord importer la liste des coordonnées dans votre projet. Ce didacticiel explique bien comment procéder: http://qgis.spatialoughtts.com/2012/01/importing-spreadsheets-or-csv-files-to.html

Lorsque vous avez les deux couches (polygones et points) dans votre projet, accédez à vecteur> outils de gestion des données> joindre les attributs par emplacement

entrez la description de l'image ici

Vous obtenez une fenêtre où vous pouvez définir les couches que vous souhaitez combiner:

entrez la description de l'image ici

  • Définissez votre couche de points comme la «couche de vecteur cible».
  • Définissez votre couche de polygones comme «joindre la couche vectorielle».
  • Définissez votre Shapefile de sortie (il y en aura un nouveau créé. Donc, si vous l'avez manqué, les données d'origine sont conservées).
  • Vous pouvez choisir de conserver toutes les données dans le nouveau fichier de formes, même s'il n'y a pas de correspondance avec un polygone: cochez «Conserver tous les enregistrements (y compris les enregistrements cibles non correspondants)»

Cliquez sur OK'. Le nouveau fichier de formes est créé et il vous sera demandé "Voulez-vous ajouter la nouvelle couche à la table des matières?" Cliquez à nouveau sur OK.

Ouvrez la table attributaire du nouveau fichier de formes ajouté et vous verrez que toutes les entités du polygone correspondant sont ajoutées au point qui se trouve dans ce polygone.

PieterB
la source
2

Une façon plus simple de le faire en utilisant PyQGIS. J'ai pensé que vous pouvez construire un QgsRectangleobjet avec un seul point et l'utiliser avec QgsFeatureRequestpour filtrer les entités du calque qui l'intersectent.

point = QgsPoint(10, 10)
# Construct a QgsRectange with the same point
rect = QgsRectangle(point, point)
req = QgsFeatureRequest()
req.setFilterRect(rect)
# You get the feature that intersects the point
f = layer.getFeatures(req).next()
pensées spatiales
la source
0

Dans QuantumGIS, vous pouvez ajouter la liste des coordonnées avec la fonction 'ajouter une couche de texte délimitée' (s'il s'agit d'un fichier csv). Ajoutez également les polygones. Ensuite, vous pouvez faire une «intersection» ou des «points dans un polygone».

Stefan
la source