Realizando una consulta espacial en un bucle en PyQGIS

Question

Lo que estoy tratando de hacer: bucle a través de un archivo de forma de punto y seleccionar cada punto que cae en un polígono.

El siguiente código está inspirado en un ejemplo de consulta espacial que encontré en un libro:

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

Esto funciona, y selecciona conjuntos de datos, pero el problema es que selecciona por el cuadro delimitador por lo tanto, obviamente, devolver los puntos que no me interesan:

¿Cómo podría ir sólo a devolver puntos dentro de el polígono sin usar qgis:selectbylocation ?

He intentado usar el dentro() y intersecta() pero como no conseguía que funcionaran, recurrí al código anterior. Pero quizás son la clave después de todo.

Answer 1

No necesitas una función especial (como "Ray Casting"), todo está en PyQGIS ( contiene() en Manejo de la geometría de PyQGIS )

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" 

o en una línea

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

También puede utilizar directamente

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

El problema aquí es que debes iterar a través de todas las geometrías (polígonos y puntos). Es más interesante utilizar un índice espacial delimitador: se itera sólo a través de las geometrías que tienen la posibilidad de intersecarse con su geometría actual ("filtro", mira ¿Cómo acceder eficientemente a las características devueltas por QgsSpatialIndex? )

Answer 2

Puedes usar el "Ray Casting" algoritmo que he adaptado ligeramente para usar con 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

Aplicado a esta situación:

el resultado, en la Consola de Python, fue:

[2, 2]

Funcionó.

Nota de edición:

Código con la propuesta más concisa de gen :

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

Answer 3

Con el consejo de un compañero de trabajo, finalmente conseguí que funcionara usando "within()".

Lógica general

1. obtener las características de los polígonos
2. obtener características de los puntos
3. a través de cada característica del archivo del polígono, y para cada una..:
   • obtener la geometría
   • y que se repita a través de todos los puntos
     • obtener la geometría de un solo punto
     • comprobar si la geometría está dentro de la geometría del polígono

Aquí está el código:

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

Esto también funcionaría con intersecta() en lugar de dentro() . Cuando se usan los puntos no importa cuál de ellas usarías, ya que ambas devolverán el mismo resultado. Sin embargo, cuando se comprueba si hay líneas/polígonos, puede haber una diferencia importante: within() devuelve objetos que son completamente dentro, mientras que intersects() devuelve los objetos que están completamente dentro y que son en parte dentro de (es decir, que intersectan con la característica, como el nombre indica).