2

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

consecutively (1, 2, 3): PHN, POA, intersection

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()
2
  • Have you tried the Processing "Overlap analysis" tool? That sounds like it already does what you're after
    – ndawson
    Dec 30, 2019 at 12:38
  • Thanks for your suggestion @ndawson sorry for my slow reply. I tried Overlap Analysis after you suggested it but it appears to only analyse one POA polygon and one PHN polygon at a time. The result is correct but it means I'd have to explode each layer into thousands of polygons and automate. If I try Overlap Analysis with more than 2 polygons the result appears to always be 100%. Jan 1, 2020 at 1:59

1 Answer 1

4

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.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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