5

I'm looking for a solution to select islands (non-black polygons) in polygon A (black).

I think, in PostGIS ST_ContainsProperly() will do this task. But I'm looking for an equivalent solution in PyQGIS.

The function shall only return True when checking the green polygon. Checking the grey MultiPolygon shall return False, because it is partially outside and the red polygon False because it touches the border

enter image description here

I tried the code from this thread

path_1 = r"C:\temp\test_poly.shp"
path_2 = r"C:\temp\test_poly_umr.shp"

layerB = QgsVectorLayer(path_1, "layerB", "ogr")
layerA = QgsVectorLayer(path_2, "layerA", "ogr")

if layerA.isValid():
    QgsProject.instance().addMapLayer(layerA)

if layerB.isValid():
    QgsProject.instance().addMapLayer(layerB)
    
island = False
for bF in layerB.getFeatures():
    for aF in layerA.getFeatures():
         if aF.geometry().intersects(bF.geometry()): 
            if bF.geometry().contains(aF.geometry()): 
                island = True
    print(bF[0], ' ', island)
    island = False

enter image description here

... but it returns island=True for Polzgon #5 instead of polygon #3 (as expected)

1
  • It would appear that both methods you tried yield the same result. Are you sure your input data is correct?
    – Matt
    Sep 27, 2021 at 15:04

2 Answers 2

7

You need to check if A's boundary intersects B's boundary. You can use this script.

# ...
# previous lines
# ...

island = False
for aF in layerA.getFeatures():
    
    # get aF's boundary (=line)
    aF_boundary = aF.geometry().convertToType(QgsWkbTypes.LineGeometry)
    
    for bF in layerB.getFeatures():        
        if aF.geometry().contains(bF.geometry()):

            # get bF's boundary (=line)
            bF_boundary = bF.geometry().convertToType(QgsWkbTypes.LineGeometry) 
       
            if not bF_boundary.intersects(aF_boundary):
                island = True
                
        print(bF[0], ' ', island)
        island = False

enter image description here

7

Use the [6] within predicate with the "Extract By Location".

container = QgsProject.instance().mapLayersByName('container')[0]
islands = QgsProject.instance().mapLayersByName('islands')[0]

selected = processing.run("native:extractbylocation",
                            {'INPUT':islands,
                            'PREDICATE':[6],
                            'INTERSECT':container,
                            'OUTPUT':'TEMPORARY_OUTPUT'}
                        )['OUTPUT']

selected.setName('selected')

QgsProject.instance().addMapLayer(selected)

In the figure, the hatched polygon is the result of Extract by Location.

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.