4 votos

Pyqgis: "a.geometry().intersects(b.geometry())" no encontraría ninguna intersección pero debería. ¿Por qué?

Estoy aprendiendo pyqgis y trato de encontrar aquellos polígonos en una capa que se cruzan con polígonos en otra capa. Parece muy trivial, y en GIS.SE hay varias preguntas similares que me sugieren (con mis limitados conocimientos de python/pyqgis) que esto debería funcionar:

#ausschluss and gebaeude are QgsVectorLayers with Polygons 
#(shapefiles with the same crs and clearly intersecting)
gebaeude.setSelectedFeatures([]) #deselects everything
selections=[] #declares it is a list
for f in gebaeude.getFeatures():
    for a in ausschluss.getFeatures():
        if a.geometry().intersects(f.geometry()): 
            #same result with .within() and even with a and f switched (the docs aren't that clear on that)
            selections.append( f.id() )
            break #only one or less intersection are possible
gebaeude.setSelectedFeatures(selections)
gebaeude.invertSelection()

#gives 0, 78137 an 78137 (no intersection found)
print len(selections)
print gebaeude.selectedFeatureCount()
print gebaeude.featureCount()

¿Qué es lo que estoy haciendo mal?

clearly intersecting...

Edición: Ejecutando QGis 2.8.3 en OS X (10.8.5) con el KyngChaos-Longterm-Built

system-stuff

5voto

Yada Puntos 9489

Utilicé una situación más trivial (la característica amarilla está indicando que fue seleccionada) para probar su código modificado:

enter image description here

Mi código era:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

print layers[0].name() #gebaeude
print layers[1].name() #ausschluss

selections=[] #declares it is a list
for f in layers[0].getFeatures():
    for a in layers[1].getFeatures():
        if a.geometry().intersects(f.geometry()):
            intersection = a.geometry().intersection(f.geometry())
            print intersection.exportToWkt()
            #same result with .within() and even with a and f switched (the docs aren't that clear on that)
            selections.append( f.id() )
            break #only one or less intersection are possible

print len(selections)
print layers[0].selectedFeatureCount() #It was selected
print layers[1].featureCount()

y corrió lo que esperaba.

Los resultados en Python Console fueron:

gebaeude
ausschluss
Polygon ((355710.56858536007348448 4472539.63295774068683386, 355785.63899155310355127 4472539.63295774254947901, 355785.0478859927970916 4472464.56255155056715012, 355711.75079648150131106 4472465.15365711320191622, 355710.56858536007348448 4472539.63295774068683386))
1
1
1

El polígono en notación WKT era gebaude (corroborado con QuickWKT plugin). Espero que esto ayude.

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