Counting overlapping polygons in GeoJSON to result in new set of polygons?

Question

I have 100 polygons in a GeoJSON file that partially overlap each other. I would like to produce 100 new polygons as follows:

• The area covered by any of the original 100. This is their union.
• The area covered by at least two of the original 100.
• ...
• The area covered by all but one of the original 100.
• The area covered by all of the original 100. This is their intersection.

(I'm actually going to only keep 5-10 of these 100 results, but I expect the process for generating all 100 of them is the same as generating a subset of them)

What software / functions can I use to do this?

---

One approach I've already encountered in another similar project involves rasterizing all of the original polygons, then counting the overlapping pixels and producing new polygons from that data. This approach loses precision earlier than I would like to (I want to maintain near perfect precision up until the end of the process, THEN simplify things for display).

---

I could do this semi-manually with hundreds of union and intersection operations. This approach is not ideal because I want to re-run this process at intervals in the future as the dataset grows.

Answer 1

I would tend to agree that a script is needed.. however, in theory one could intersect everything, keeping the ID attribute of every parent polygon. Then each fragments "parentage" and thus any combination of parents could be established through queries.

This is basically a vector version of the raster approach. Given the potential problems with sliver polygons, likely limits in vector topology precision and area calculation, not to mention processing time, I would not discount the raster version. If a sufficient resolution were feasible the final precision might not be much worse - or would atleast be more easily assessed.

Answer 2

What you need is a script, which will perform a certain action on every polygon. You cannot achieve this with usual GIS application functions. So depending on your usual GIS application you have to see which scripting support does it have.

In general Phyton is the most used, either by proprietary applications (like ArcGIS) or by open source applications (QGIS, GRASS GIS). In GRASS 6.4 you could still do it in Bash (even in Windows), and if you use SAGA GIS under Windows you could do it Windows Command Line.

Example in pseudocode:

polygonsID = ['1', '2', ... , 'n'] where n is the maximum polygon ID
polygons = ['polygon1', 'polygon2', ... , 'polygonn']   

for p in polygon
    union(polygons[w],polygons[w+1])

You have a reference to the union method in QGIS here

Cheers,

Answer 3

This can be done by performing many unions and intersections separately.

Say there are 4 polygons, named A,B,C,D. The 4 results that are desired are as follows: (& indicates intersection, | indicates union)

A|B|C|D

(A|B)&(A|C)&(A|D)&(B|C)&(B|D)&(C|D)

(A|B|C)&(A|B|D)&(A|C|D)&(B|C|D)

A&B&C&D

This approach suffers from a factorially increasing number of intersection and union operations in the middle of the process for larger numbers of polygons. Memoization will reduce the number of operations.

Answer 4

Break down the set of polygons into a larger set of non-overlapping polygons.

For every polygon in that set, check whether it is part of each original polygon (using a point-in-polygon function if possible) and keep count of them.

Keep only polygons that are part of the requisite number of original polygons.

Take the union of those remaining polygons.

This produces the desired solution (which may be a polygon, or a set of polygons, or an empty set).

Answer 5

Try this:

# Count Overlapping Polygons by Sections
# Author: Ervin Wirth, 2016

from PyQt4.QtCore import QVariant
import itertools

filename1 = "buffer2"
path1="d:/Works/ChargingStation/data/test/"
layer1 = QgsVectorLayer(path1 + filename1 +".shp", filename1, "ogr")

save_layer = QgsVectorLayer("Polygon?crs=epsg:23700", "save_layer", "memory")
save_layer.dataProvider().addAttributes([QgsField("id",  QVariant.Int), QgsField("k",  QVariant.Int)])
save_layer.updateFields()

# create a recursive function
# arg1: input layer; arg2: level of section; arg3: threshold for sliver polygons
def iterlay(layer1, k, thres):

    # create a memory layer for intersections
    mem_layer = QgsVectorLayer("Polygon?crs=epsg:23700", "temp_layer", "memory")
    mem_layer.dataProvider().addAttributes([QgsField("id",  QVariant.Int), QgsField("k",  QVariant.Int)])
    mem_layer.updateFields()
    id = mem_layer.fieldNameIndex('id')

    # prepare a loop on the input layer
    polygons = [feature for feature in layer1.getFeatures()]
    list = range(len(polygons))
    print len(polygons)

    # create the sections in all combinations and store them temporary
    for i,j in itertools.combinations(list, 2):
        # if polygon has section, than intersecting
        if polygons[i].geometry().intersects(polygons[j].geometry()):
            geom = polygons[i].geometry().intersection(polygons[j].geometry())
            # intersections can be multipolygons
            if geom.asMultiPolygon():
                for singlepoly in geom.asMultiPolygon():
                    qgisgeom = QgsGeometry.fromPolygon (singlepoly)
                    if qgisgeom.area() > thres:
                        feature = QgsFeature()
                        feature.setAttributes([str(k) + '-' + str(i) + '-' + str(j) ,str(k)])
                        feature.setGeometry(qgisgeom)
                        mem_layer.dataProvider().addFeatures([feature])

            # and simple polygons
            if geom.asPolygon():
                qgisgeom = QgsGeometry.fromPolygon (geom.asPolygon())
                if qgisgeom.area() > thres:
                    feature = QgsFeature()
                    feature.setAttributes([str(k) + '-' + str(i) + '-' + str(j) ,str(k)])
                    feature.setGeometry(qgisgeom)
                    mem_layer.dataProvider().addFeatures([feature])

    print mem_layer.featureCount()

    # remove duplicates from sections, and eliminate sliver polygons
    sections = [feature for feature in mem_layer.getFeatures()]
    list = range(len(sections))
    dlist=[]
    for i,j in itertools.combinations(list, 2):
        # compare the shapes (maybe can be modified to equals())
        if sections[i].geometry().difference(sections[j].geometry()).area() < thres and sections[j].geometry().difference(sections[i].geometry()).area() < thres:
            dlist.append(sections[j].id())
        if sections[i].geometry().area() < thres:
            dlist.append(sections[i].id())
    myset = set(dlist)
    for dupl in myset:        
        mem_layer.startEditing()
        mem_layer.deleteFeature(dupl)
        mem_layer.commitChanges()
    print mem_layer.featureCount()

    # save the remained parts
    sections = [feature for feature in mem_layer.getFeatures()]
    for poly in polygons:
        # have to create the differences with all the sections
        tempgeom = poly.geometry()
        for section in sections:
            newgeom = tempgeom.difference(section.geometry())
            tempgeom = newgeom
        # if remained something, write out
        if tempgeom.asPolygon():
            qgisgeom = QgsGeometry.fromPolygon (tempgeom.asPolygon())
            if qgisgeom.area() > thres:
                feature = QgsFeature()
                feature.setAttributes([str(poly.id()) ,str(k)])
                feature.setGeometry(qgisgeom)
                save_layer.dataProvider().addFeatures([feature])
        if tempgeom.asMultiPolygon():
            for singlepoly in tempgeom.asMultiPolygon():
                qgisgeom = QgsGeometry.fromPolygon (singlepoly)
                if qgisgeom.area() > thres:
                    feature = QgsFeature()
                    feature.setAttributes([str(poly.id()) ,str(k)])
                    feature.setGeometry(qgisgeom)
                    save_layer.dataProvider().addFeatures([feature])

    # have to stop when there isn't any section
    if mem_layer.featureCount() == 0:
    #if k == 2:
        QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
        QgsMapLayerRegistry.instance().addMapLayer(save_layer)
    else:
        k += 1
        iterlay(mem_layer, k, thres)

# Call the funtion
iterlay(layer1,1,100000)

Answer 6

One possible approach rasterizing all of the original polygons, then counting the overlapping pixels and producing new polygons from that data.

The disadvantage of this approach is that you loses precision quite early early (in the first rasterisation itself) in the process.

This might not be a suitable solution when you need to to maintain near perfect precision up until the end of the process, & then simplify things for display.