3

I'd like to calculate the percentage of overlap (not area) of polygons in two different layers. Both layers are in WGS84 (EPS:4326). There are several tutorials on how to calculate area overlap for projected layers, e.g. Calculating percentage of overlap between two different layers, however, I'd like to avoid having to reproject all layers. Is there a tool in QGIS to do this?

  • 1
    Have you considered using PostGIS functions? You can achieve this with a simple query. You will have to reproject the data, but it'll be on runtime without all the fuzz of creating a lot of temporary files. – fabio.avigo Jul 21 '17 at 18:49
  • I am unfortunately not familiar with PostGIS and some of the layers are also pretty big, this is why I'd like to avoid having to reproject – Kristina Jul 21 '17 at 18:51
  • In QGIS, to avoid having to reproject all layers with EPSG:4326 for determining areas, you can use a PyQGIS script with QgsDistanceArea objects. – xunilk Jul 22 '17 at 0:04
1

In QGIS, to avoid having to reproject all layers with EPSG:4326 for determining areas, you can use a PyQGIS script with QgsDistanceArea objects. Next script was generalized for two layers with several features but, in my example (see next image), there is only one feature each. It calculates overlap percentage related to overall area of all features in both layers.

registry = QgsMapLayerRegistry.instance()

layer1 = registry.mapLayersByName('polygon1_WGS84')
layer2 = registry.mapLayersByName('polygon2_WGS84')

feats_layer1 = [ feat for feat in layer1[0].getFeatures() ] 
feats_layer2 = [ feat for feat in layer2[0].getFeatures() ]

area = QgsDistanceArea()
area.setEllipsoid('WGS84')
area.setEllipsoidalMode(True)

areas1 = [ area.measurePolygon(feat.geometry().asPolygon()[0]) for feat in feats_layer1 ]
areas2 = [ area.measurePolygon(feat.geometry().asPolygon()[0]) for feat in feats_layer2 ]

overall_area = sum(areas1) + sum(areas2)

for feat1 in feats_layer1:
    for feat2 in feats_layer2:
        if feat1.geometry().intersects(feat2.geometry()):
            intersection = feat1.geometry().intersection(feat2.geometry())

            percentage_overlap = (area.measurePolygon(intersection.asPolygon()[0])/overall_area)*100

            print "percentage_overlap: {:.2f} %".format(percentage_overlap)

        else:
            "there is not intersection"

After running above script for next situation (observe that coordinates status bar are in degrees):

enter image description here

it was printed a overlapping percentage of 10.04% for this case.

enter image description here

To corroborate that these areas are meaningful, it was printed areas1, areas2 and overall_area separately. These next values are in meters.

>>>areas1
[609428494.6844916]
>>>areas2
[364922280.7505144]
>>>overall_area
974350775.435006  #974.35 km2

Editing Note:

An additional corroboration includes world_borders shapefile. A country with only one feature is Albania; whose area is 28748 km2.

enter image description here

Next shorter code print area (meters) of selected feature (Albania):

layer = iface.activeLayer()

area = QgsDistanceArea()
area.setEllipsoid('WGS84')
area.setEllipsoidalMode(True)

feat = layer.selectedFeatures()

print area.measurePolygon(feat[0].geometry().asPolygon()[0])

Its value was 28728106391.5 m2 = 28728.1063915 km2. Assuming above value (attributes table) as true, the percentage error is (abs(28728.1063915 - 28748)/28748)*100 = 0.069 %; which is really very low.

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.