2

I am currently using geometry().intersects(geom)) in one of my PyQGIS scripts, but it is not returning the outcome I'd like, and I don't think .touches(geom) would either.

I only want to return neighboring polygons which share an edge with the currently selected polygons. I've attached an image below.

My target polygon is in red. Right now it is returning all of the polygons in black and blue, but I only want it to return the polygons in blue. If it is not clear from the image, the red and black polygons touch at a point.

enter image description here

Is there a way to return only polygons which share an edge and not a point, i.e. just the two blue polygons in the image and not the black polygon, using PyQGIS?

1
  • We're a little different from other sites; this isn't a discussion forum but a Q&A site. Please check out our short tour which emphasizes that there should be only one question asked per question. If you decide to ask about PyQGIS then please include a code snippet in your question.
    – PolyGeo
    Jan 22, 2020 at 7:09

1 Answer 1

2

If you use itertools python module and print all possibilities for intersections, you will have 2 points (0 as type geometry) and 4 MultiLineString geometries (1 as type geometry). You are only interested in MultiLineString geometries and target polygon id (in my case it was 0) for determining features indices of your "blue" polygons.

So, with following python script:

import itertools

registry = QgsProject.instance()

layer = registry.mapLayersByName('four_states')

features = [ feat for feat in layer[0].getFeatures() ]

list = range(len(features))

#intersections = [[ features[i].id(), features[j].id(), features[i].geometry().intersection(features[j].geometry()).asWkt() ]
#                 for i,j in itertools.combinations(list, 2)
#                 if features[i].geometry().intersects(features[j].geometry()) ]

intersections = [[ features[i].id(), features[j].id(), features[i].geometry().intersection(features[j].geometry()).type() ]
                 for i,j in itertools.combinations(list, 2)
                 if features[i].geometry().intersects(features[j].geometry()) ]

#item[0] target index, item[1] indices other polygons, item[2] type geometry
indices = [ item[1] for item in intersections if item[0] == 0 and item[2] != 0 ]

#creating a memory layer from polygon
crs = layer[0].crs()
epsg = crs.postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "polygons",
                           "memory")

prov = mem_layer.dataProvider()

#Set features
new_features = [ features[id] for id in indices ]

for i, feat in enumerate(new_features):
    feat.setAttributes([i])

prov.addFeatures(new_features)

QgsProject.instance().addMapLayer(mem_layer)

and layer of following image:

enter image description here

after running above script it can be obtained expected result (ocher color polygons in following image) as memory layer:

enter image description here

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.