Exécution d'une requête spatiale en boucle dans PyQGIS

Ce que je suis en train de le faire: la boucle par un point shapefile et sélectionnez chaque point qui tombe dans un polygone.

Le code suivant est inspiré d'un exemple de requête spatiale que j'ai trouvé dans un livre:

mitte_path = r"D:\PythonTesting\SelectByLocation\mitte.shp"
punkte_path = r"D:\PythonTesting\SelectByLocation\punkte.shp"

polygon = QgsVectorLayer(mitte_path, 'Mitte', 'ogr')
points = QgsVectorLayer(punkte_path, 'Berlin Punkte', 'ogr')

QgsMapLayerRegistry.instance().addMapLayer(polygon)
QgsMapLayerRegistry.instance().addMapLayer(points)

polyFeatures = polygon.getFeatures()

pointsCount = 0

for poly_feat in polyFeatures:
    polyGeom = poly_feat.geometry()
    pointFeatures = points.getFeatures(QgsFeatureRequest().setFilterRect(polyGeom.boundingBox()))
    for point_feat in pointFeatures:
        points.select(point_feat.id())
        pointsCount += 1

print 'Total:',pointsCount

Cela fonctionne, et il sélectionne des jeux de données, mais le problème est qu'il sélectionne par boîte englobante , donc renvoyant évidemment des points qui ne m'intéressent pas:

Comment pourrais-je ne retourner que des points dans le polygone sans utiliser qgis: selectbylocation ?

J'ai essayé d'utiliser les méthodes within () et intersects () , mais comme je ne les faisais pas fonctionner, j'ai eu recours au code ci-dessus. Mais peut-être qu'ils sont la clé après tout.

Vous n'avez pas besoin d'une fonction spéciale (comme "Ray Casting"), tout est dans PyQGIS ( contient () dans PyQGIS Geometry Handling )

polygons = [feature for feature in polygons.getFeatures()]
points = [feature for feature in points.getFeatures()]
for pt in points: 
     point = pt.geometry() # only and not pt.geometry().asPolygon() 
     for pol in polygons:
        poly = pol.geometry()
        if poly.contains(point):
             print "ok" 

ou en une seule ligne

 polygons = [feature for feature in polygons.getFeatures()]
 points = [feature for feature in points.getFeatures()]
 resulting = [pt for pt in points for poly in polygons if poly.geometry().contains(pt.geometry())]
 print len(resulting)
 ...

Vous pouvez également utiliser directement

[pt.geometry().asPoint() for pt in points for poly in polygons if poly.geometry().contains(pt.geometry())]

Le problème ici est que vous devez parcourir toutes les géométries (polygones et points). Il est plus intéressant d'utiliser un index spatial englobant: vous itérez uniquement à travers les géométries qui ont une chance de se croiser avec votre géométrie actuelle ('filtre', regardez Comment accéder efficacement aux fonctionnalités retournées par QgsSpatialIndex? )

Vous pouvez utiliser l' algorithme "Ray Casting" que j'ai légèrement adapté pour l'utiliser avec PyQGIS:

def point_in_poly(point,poly):
    x = point.x()
    y = point.y()

    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

## Test
mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#For polygon 
polygon = [feature.geometry().asPolygon() 
            for feature in layers[1].getFeatures()]

points = [feat.geometry().asPoint() 
           for feat in layers[0].getFeatures()]

## Call the function with the points and the polygon
count = [0]*(layers[1].featureCount())

for point in points:
    i = 0
    for feat in polygon:
        if point_in_poly(point, feat[0]) == True:
            count[i] += 1
        i += 1

print count

Appliqué à cette situation:

le résultat, sur la console Python, était:

[2, 2]

Ça a marché.

Note d'édition:

Code avec la proposition de gène la plus concise :

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

count = [0]*(layers[1].featureCount())

polygon = [feature
           for feature in layers[1].getFeatures()]

points = [feature
          for feature in layers[0].getFeatures()]

for point in points:

    i = 0

    geo_point = point.geometry()

    for pol in polygon:
        geo_pol = pol.geometry()

        if geo_pol.contains(geo_point):
            count[i] += 1
        i += 1

print count

Avec quelques conseils d'un collègue de travail, je l'ai finalement fait fonctionner en utilisant ().

Logique générale

1. obtenir des caractéristiques de polygone (s)
2. obtenir les caractéristiques des points
3. parcourir chaque entité à partir du fichier polygone et pour chacune:
   
   • obtenir la géométrie
   • parcourir tous les points
     • obtenir la géométrie d'un point unique
     • tester si la géométrie est dans la géométrie du polygone

Voici le code:

mitte_path = r"D:\PythonTesting\SelectByLocation\mitte.shp"
punkte_path = r"D:\PythonTesting\SelectByLocation\punkte.shp"

poly = QgsVectorLayer(mitte_path, 'Mitte', 'ogr')
points = QgsVectorLayer(punkte_path, 'Berlin Punkte', 'ogr')

QgsMapLayerRegistry.instance().addMapLayer(poly)
QgsMapLayerRegistry.instance().addMapLayer(points)

polyFeatures = poly.getFeatures()
pointFeatures = points.getFeatures()

pointCounter = 0

for polyfeat in polyFeatures:
    polyGeom = polyfeat.geometry()
    for pointFeat in pointFeatures:
        pointGeom = pointFeat.geometry()
        if pointGeom.within(polyGeom):
            pointCounter += 1
            points.select(pointFeat.id())

print 'Total',pointCounter

Cela fonctionnerait également avec intersects () au lieu de within () . Lorsque vous utilisez des points, peu importe celui que vous utiliseriez, car ils renverront tous les deux le même résultat. Cependant, lors de la vérification des lignes / polygones, cela peut faire une différence importante: within () renvoie des objets qui sont entièrement à l' intérieur, tandis que intersects () réaccorde les objets qui sont entièrement à l'intérieur et qui sont partiellement à l' intérieur (c'est-à-dire qui se croisent avec l'entité, comme le nom l'indique).