Returning neighboring polygons which share a boundary, not just a point

Question

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.

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?

Answer 1

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:

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