9 votos

La realización de una Consulta Espacial en un bucle en PyQGIS

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

El siguiente código está inspirado en una consulta espacial de ejemplo que he encontrado 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 lo hace de seleccionar los conjuntos de datos, pero el problema es que se selecciona mediante el cuadro delimitador, por lo tanto, obviamente, de devolver los puntos no me interesan:

enter image description here

Cómo se podría ir sobre solo devuelve los puntos dentro de la poligonal sin usar qgis:selectbylocation?

He probado a utilizar el plazo de() y intersects() métodos, pero como yo no era llegar a trabajar, utilicé el código anterior. Pero tal vez son la clave, después de todo.

7voto

GreyCat Puntos 146

Usted no necesita una función especial (como "Ray Casting"), todo está en PyQGIS (contiene() en PyQGIS Geometría de Manejo)

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 se 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 usted debe iterar a través de todas las geometrías (polígonos y puntos). Es más interesante el uso de una delimitación espacial del índice : iterar sólo a través de las geometrías que tienen una oportunidad para intersecar con su actual geometría ('filtro', mira Cómo acceder eficientemente a las características devuelto por QgsSpatialIndex?)

5voto

Yada Puntos 9489

Usted puede utilizar el "Ray Casting" algoritmo que he adaptado un poco 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

Aplica a esta situación:

enter image description here

el resultado, en la Consola de Python, fue:

[2, 2]

Se trabajó.

La Edición De La Nota:

Código con el más conciso de la propuesta 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

2voto

mercator Puntos 16196

Con algunos consejos de un compañero de trabajo por fin lo conseguí con un plazo de().

Lógica General

  1. obtener las características de un polígono(s)
  2. obtener las características de los puntos de
  3. el bucle a través de cada una de las características del polígono de archivo, y para cada uno de ellos:
    • obtener la geometría
    • el bucle a través de todos los puntos
      • obtener la geometría de punto único
      • prueba si la geometría está dentro de la geometría de 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 podría trabajar con intersects() en lugar de dentro de los(). Al utilizar los puntos que no importa que uno utilizaría, como ellos le devuelven el mismo resultado. Cuando la comprobación de líneas/polígonos, sin embargo, puede hacer una diferencia importante: en el plazo de() devuelve los objetos que están completamente dentro, mientras que intersects() retuns los objetos que están completamente dentro y que son en parte dentro (es decir, que se cruzan con la característica, como el nombre lo indica).

enter image description here

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X