4

I am teaching myself pyqgis and try to find those polygons in one layer that intersect with polygons in another layer. It seems very trivial, and on GIS.SE are several similar questions that are suggesting me (with my limited python/pyqgis-knowledge) that this should work:

#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()

What is it, I am doing wrong?

clearly intersecting...

Edit: Running QGis 2.8.3 on OS X (10.8.5) with the KyngChaos-Longterm-Built

system-stuff

  • What happens if you switch this statement: if a.geometry().intersects(f.geometry()): to if f.geometry().intersects(a.geometry()): – artwork21 Oct 29 '15 at 13:14
  • 1
    I really think that's a CRS problem ... r u sure these two layers aren't just visually overlapping (because of the OTF reprojection) but have different CRSs ? Tell us when u find the solution ... – snaileater Oct 29 '15 at 15:24
  • @artwork21 - that results in the same – Janne Jakob Fleischer Oct 30 '15 at 9:13
  • @snaileater - I think they are really identical. I even checked the geometry with isGeosValid() - and it is. Both .prj-Files (thats where the crs is saved, right?) say: PROJCS["DHDN_3_degree_Gauss_Kruger_zone_3",GEOGCS["GCS_DHDN",DATUM["D_Deutsches_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",1],PARAMETER["false_easting",3500000],PARAMETER["false_northing",0],UNIT["Meter",1]] – Janne Jakob Fleischer Oct 30 '15 at 9:16
5

I used a more trivial situation (yellow feature is indicating that it was selected) for testing your modified code:

enter image description here

My code was:

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

and it ran I expected.

Results at Python Console were:

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

Polygon in WKT notation was gebaude (corroborated with QuickWKT plugin). I hope that it helps.

  • thanks, xunilk - Now I know my code is right, but my geometry isn't, so it seems. Thats a start on where to look further. – Janne Jakob Fleischer Oct 30 '15 at 9:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.