PyQGIS layer intersection area seems really inaccurate

Question

Using 3.10.1-A Coruña on Windows 10

I have a layer of Public Hospital Networks PHN and I have a layer of zip/ postcodes POA

I want to find what percentage of a zip/ postal code is in a PHN ParentLayer

I'm using the following code but the intersecting percentage seems way off. It says 4% but from eyeballing it should be 40-60%.

Percent in Parent Region:  4.242022210197979e-15

I tried intersection().asPolygon().area() but that failed and searches suggested I need to try Multiparts to Singleparts, which I tried, but still gave the same error.

I checked if I had Valid Polygons and everything indicated OK.

I Frankensteined the code below from Calculating and using area of polygons in PyQGIS? and Calculate percentage overlap of polygons in WGS84

output_file = open("D:\\Folder\\Intersection.csv", 'w')
# POA, ParentRegion, IntersectionResult, IntersectionPercentage

layer_list = QgsProject.instance().layerTreeRoot().children() # returns QgsLayerTreeNode object list

POA_Layer = [lyr.layer() for lyr in layer_list if lyr.name()=="  POA_2016_AUST"][0]
ParentRegion_Layer = [lyr.layer() for lyr in layer_list if lyr.name()=="  PHN_boundaries_AUS_May2017_V7"][0]

print("POA layer", POA_Layer)
print("Parent layer", ParentRegion_Layer)

d = QgsDistanceArea()
d.setEllipsoid('WGS84')
#d.setEllipsoidalMode(True)

# single features for testing
POA_Query = '"POA_CODE16" = \'2157\''
POA_Features = POA_Layer.getFeatures(QgsFeatureRequest().setFilterExpression(POA_Query))

ParentRegion_Query = '"FIRST_PHN_" = \'PHN103\''
ParentRegion_Features = ParentRegion_Layer.getFeatures(QgsFeatureRequest().setFilterExpression(ParentRegion_Query))

# the big loop
for POA_feature in POA_Features:
    POA_Geom = POA_feature.geometry()
    POA_Name = POA_feature.attribute('POA_NAME16')
    print("POA name: ", POA_Name)
    print("POA area (Km2):", d.convertAreaMeasurement(d.measureArea(POA_Geom), QgsUnitTypes.AreaSquareKilometers))

    for ParentRegion_feature in ParentRegion_Features:
        ParentRegion_Geom = ParentRegion_feature.geometry()
        ParentRegion_Name = ParentRegion_feature.attribute('FIRST_PHN_')
        print("Parent name: ", ParentRegion_Name)
        print("Parent area (Km2):", 
        d.convertAreaMeasurement(d.measureArea(ParentRegion_Geom), QgsUnitTypes.AreaSquareKilometers))

        if POA_feature.geometry().intersects(ParentRegion_feature.geometry()):
        print ("Intersection detected")

        intersection = POA_feature.geometry().intersection(ParentRegion_feature.geometry())

        if ParentRegion_feature.geometry().contains(POA_feature.geometry()):
            print (POA_Name + " is fully enclosed by " + ParentRegion_Name)
            # csv output 
            output_file.write(POA_Name + "," + ParentRegion_Name + ",FullyEnclosed, 100")
        else:
                intersection_area = d.convertAreaMeasurement(POA_feature.geometry().intersection(ParentRegion_feature.geometry()).area(), QgsUnitTypes.AreaSquareKilometers)
            # ERROR BELOW: GeometryCollection geometry cannot be converted to a polygon. Only single polygon or curve polygon types are permitted.
            # intersection_area_polygon = POA_feature.geometry().intersection(ParentRegion_feature.geometry()).asPolygon().area()
            print ("Intersection area: ", intersection_area)
            print ("Percent in Parent Region: ", (intersection_area/d.measureArea(POA_Geom)*100))
            # csv output 
            output_file.write(POA_Name + "," + ParentRegion_Name + ",Percentage, " + str(intersection_area))
    else:
            print ("There is no intersection")
            output_file.write(POA_Name + "," + ParentRegion_Name + ",NoIntersection, 0")

output_file.close()

Answer 1

I suspect that the cause of the unexpected result may be coming from the line where you calculate and print the percentage:

print ("Percent in Parent Region: ", (intersection_area/d.measureArea(POA_Geom)*100))

The intersection_area object has been converted to km2, but you are dividing it by d.measureArea(POA_Geom) which has not been converted to km2. You have printed the conversion here:

print("POA area (Km2):", d.convertAreaMeasurement(d.measureArea(POA_Geom), QgsUnitTypes.AreaSquareKilometers))

but haven't stored the result in a variable.

Try the code below and see if this gives the desired results:

import csv
output_file = open('D:\\Folder\\Intersection.csv', 'w', newline='')
writer = csv.writer(output_file)
writer.writerow(['POA Feature'] + ['Parent Region'] + ['Percentage'])

POA_layer = QgsProject().instance().mapLayersByName('POA_2016_AUST')[0]
parent_region_layer = QgsProject().instance().mapLayersByName('PHN_boundaries_AUS_May2017_V7')[0]

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

POA_features = [f for f in POA_layer.getFeatures()]
parent_region_features = [f for f in parent_region_layer.getFeatures()]

for POA_feature in POA_features:
    POA_name = POA_feature['POA_NAME16']
    total_area = d.convertAreaMeasurement(d.measureArea(POA_feature.geometry()), QgsUnitTypes.AreaSquareKilometers)
    for parent_region_feature in parent_region_features:
        parent_region_name = parent_region_feature['FIRST_PHN_']
        if POA_feature.geometry().intersects(parent_region_feature.geometry()):
            intersection = POA_feature.geometry().intersection(parent_region_feature.geometry())
            intersection_km2 = d.convertAreaMeasurement(d.measureArea(intersection), QgsUnitTypes.AreaSquareKilometers)
            pcnt = (intersection_km2/total_area)*100
            print('Percentage of {} in parent region {}: {}%'.format(POA_name, parent_region_name, pcnt))
            writer.writerow([str(POA_name)] + [str(parent_region_name)] + [str(pcnt)])
        elif POA_feature.geometry().within(parent_region_feature.geometry()):
            print('{} is fully enclosed by {}'.format(POA_name, parent_region_name))
            writer.writerow([str(POA_name)] + [str(parent_region_name)] + ['Fully enclosed- 100'])
        else:
            print('There is no intersection')
            writer.writerow([str(POA_name)] + [str(parent_region_name)] + ['No intersection- 0'])
output_file.close()

This seemed to give accurate results for me with a couple of test layers. Hopefully it will do what you are after.