3

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.

  • It's just a thought but have you considered converting your GeoJSON to TopoJSON? That may or may not be an overlap aware format but it at least has the potential to be able to. – PolyGeo Jul 29 '13 at 1:12
  • @PolyGeo converting it is definitely an option. I think converting to TopoJSON or even SHP or another format would be pretty easy to automate. – Sparr Jul 29 '13 at 14:22
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.

  • I think doing a complete set of intersections as you suggest might produce prohibitively many result polygons (millions, in some cases). – Sparr Jul 29 '13 at 16:48
  • It seems like the problem will be inherently of high complexity, so either preprocessd to large datasize, or slow per operation. A spatial index will help either way. – AnserGIS Jul 29 '13 at 19:32
  • Actually, a left field option might be to sequentially add all the polygons to a TIN. Each triangle attribute could hold all the IDs and the TIN would provide a spatial topology to simplify reconstructing overlaps. – AnserGIS Jul 29 '13 at 19:39
  • @AnswerGIS TIN? – Sparr Jul 30 '13 at 3:51
  • ArcGIS probably couldn't handle this without allot of coding work work but in theory you triangulate the first polygon set such that each element will have its own ID and the attribute of the parent polygon. Then embed the second set into that TIN, adding the new set of IDs either as an attribute or to CSV string field. The TIN will handle (though not solve) likely topological errors from sliver polygons. – AnserGIS Aug 12 '13 at 14:22
0

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.

  • This site used rasterization and then binning on a hex grid to produce results similar to what I'm looking for: bostonography.com/2012/… – Sparr Jul 28 '13 at 22:05
  • @Sparr, please use "comments" to discuss the answers. If you want to add something to your question, just edit it, marking the text with Update:. – Alexandre Neto Jul 29 '13 at 9:29
  • @AlexandreNeto can you elaborate on this? I received two different answers to this question when I asked it outside of SE. I posted those answers here to avoid others' wasting their time posting answers I already had (and that weren't appropriate for my situation). – Sparr Jul 29 '13 at 15:11
  • 1
    In gis.stackexchange (like other Q&A sites), you should add extra information (like answers you got or some article you read meanwhile) editing your question. Some people, like to add a little subtitle to inform that a certain part of the questions is an update to the original one. Having all info in your question, makes it easy for other to understand your problem and propose solutions for it. If you find an answer by your own, that's a different case, you can\should answer your own question, and mark it as accepted. Others might find your Answer useful for them self. – Alexandre Neto Jul 29 '13 at 16:32
  • @AlexandreNeto that does not mesh with the suggested and practiced self-answer methods on other SE sites. This is not part of my question. It is an answer. If anything, I should have put "This approach loses precision" in a comment on this answer. – Sparr Jul 29 '13 at 16:45
0

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,

  • I am comfortable coding in Python and am learning to use QGIS, but I'm still not sure how to approach this problem even in that context. – Sparr Jul 29 '13 at 16:48
  • In the response above you have the idea of the processing what needs to be done. So you need Python, you need to make a list with the ID number of the polygons, then apply for every ID a certain function. You will need additional modules to handle the GeoJson and the union/intersection. These module can be taken from QGIS/GRASS GIS/ArcGIS. – Niculita Mihai Jul 31 '13 at 5:57
  • Say I have 5 polygons, A B C D E, and I want to know the region covered by at least 3 of them. I will use | to represent union and & to represent intersection. Here are the operations required: (A&B&C)|(A&B&D)|(A&B&E)|(A&C&D)|(A&C&E)|(A&D&E)|(B&C&D)|(B&C&E)|(B&D&E)|(C&D&E) 20 intersection operations (15 with memoization), 10 unions. For 20 polygons, the 10-or-more case requires 1662804 intersections and 184756 unions. Doing this with 100 polygons is infeasible. With 1000 impossible. I need a better algorithm if I'm going to script this myself. – Sparr Jul 31 '13 at 13:51
0

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.

0

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).

  • 1
    Please take a look at the "Union" geoprocessing function for vector layers: its attribute table contains all the information you need without any further or iterative spatial processing. – whuber Jul 31 '13 at 20:41
  • 1
    whuber is exactly right. Union does everything you need in a single step. – blord-castillo Aug 1 '13 at 3:25
  • @whuber which software package are you referring to? Union, as a simple boolean operation, is available in almost every vector and GIS application. In the ones I've tried, the result is just another [multi]polygon, without the data I'm looking for. – Sparr Aug 1 '13 at 15:52
  • Because you tagged your question with qgis, I ran this operation in QGIS 1.8.0 before posting my comment. It is available as Vector|Geoprocessing Tools|Union. It is not your "simple boolean operation," but rather a full-blown implementation of the (misnamed) "union" operation of GIS, whereby each pair of overlapping polygons introduces one record in the output containing identifying information for both polygons and is linked to a shape corresponding to the intersection of those two polygons. – whuber Aug 1 '13 at 16:17
0

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)

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.