Return to Question

I'm trying to understand the ArcGIS "Intersect" method so I can replicate it with an algorithm using Python and Shapely.

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image I want ABC returned three times, and AB, AC, and BC returned twice.

I've tried just looping and intersecting each geometry with all other geometries, but this only considers one and one intersection at a time, it seems the intersection results also have to be intersected with all other geometries, their intersection results, and so on. This is where I get confused and where it will possibly become too slow to be useful.

Hopefully there is a simple logic that im just not realizing, maybe through a combination of intersections and differences, or recursive functions.

I'm trying to understand the ArcGIS "union" method so I can replicate it with an algorithm using Python and Shapely.

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image I want ABC returned three times, and AB, AC, and BC returned twice. While the ArcGis method also returns the unique parts A, B, and C that have no intersections, I do not need this.

Code

In the ArcGis docs for union, they write that they use a "cracks and clusters" method, which is likely a low level modification of the intersection/union clipping algorithm. Since I am not implementing the underlying clipping algorithm, I need a different approach.

Here is pseudocode that I will be trying to implement. It basically means for each feature i find all intersections with others, and then recursively find all intersections between the intersecting parts, adding only the "node" parts that have no intersections with other intersections (i.e. we have cut it up as much as possible):

def isecs(g, geoms):
    for og in geoms:
        if og != g:
            if og.crosses(g):
                yield g.intersection(og)

def process(isecs):
    parts = []
    for g in isecs:
        isecs = getisecs(g, isecs)
        if no isecs:
            # "node" reached
            parts += g
        else:
            # this is the recursive part
            parts += process(isecs)
    return parts

for f in features:
    top_isecs = isecs(f.geom, features.bbox_overlaps(f.bbox))
    parts = process(top_isecs)
    for g in parts:
        addfeat(g)
            

Correction

@Vince correctly pointed out that the ArcGis operation is called union, instead of intersection as I first wrote. This is not to be confused that we are looking for geometrical unions, ArcGis is simply referring to the fact that they are returning all (hence union) geometrical intersections: "Union calculates the geometric intersection of any number of feature class..." They also allow returning an actual geometrical union through a dissolve option.

Im trying to understand the ArcGis "intersect" method so i can replicate it with an algorithm using Python and Shapely.

I do not want the cumulative intersection of all features where all of them intersect at the same time (ABC):

(Image is from another SO post)

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image i want ABC returned three times, and AB, AC, and BC returned twice.

This is a pretty standard gis method also available in QGis, so what is the accepted set theoretic way of achieving this?

Ive tried just looping and intersecting each geometry with all other geometries, but this only considers one and one intersection at a time, it seems the intersection results also have to be intersected with all other geometries, their intersection results, and so on. This is where i get confused and where it will possibly become too slow to be useful.

Hopefully there is a simple logic that im just not realizing, maybe through a combination of intersections and differences, or recursive functions.

I'm trying to understand the ArcGIS "Intersect" method so I can replicate it with an algorithm using Python and Shapely.

I do not want the cumulative intersection of all features where all of them intersect at the same time (ABC):

(Image is from another SO post)

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image I want ABC returned three times, and AB, AC, and BC returned twice.

This is a pretty standard GIS method also available in QGIS, so what is the accepted set theoretic way of achieving this?

I've tried just looping and intersecting each geometry with all other geometries, but this only considers one and one intersection at a time, it seems the intersection results also have to be intersected with all other geometries, their intersection results, and so on. This is where I get confused and where it will possibly become too slow to be useful.

Hopefully there is a simple logic that im just not realizing, maybe through a combination of intersections and differences, or recursive functions.