# Algorithm for all intersecting geometry parts

I'm trying to understand the ArcGIS "union" 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. While the ArcGis method also returns the unique parts A, B, and C that have no intersections, I do not need this.

This is a pretty standard GIS method also available in QGIS, so what is the accepted set theoretic way of achieving 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:
``````

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.

• You mention looping and functions, can you post a code sample of what you've tried and note where you are stuck at? Commented Nov 1, 2016 at 18:42
• All the information you need is located in the Union output. You just need to make `N!` passes through the file, for every permutation of presence/absence of each feature. Commented Nov 1, 2016 at 19:31
• @Vince what youre saying is exactly what i need. How do i determine the N passes (what am i looping over), and how would i get the permutations of the presence/absence of each feature? Not sure what you mean about the union output, a union would obscure all the internal intersection lines... If you could expand on your comment in an answer i would greatly appreciate it. Commented Nov 1, 2016 at 20:00
• Preferably with pseudocode if possible Commented Nov 1, 2016 at 20:14
• Loop through and create ABC, then loop through again and difference ABC from each AB BC CA union? Commented Nov 1, 2016 at 21:42

``````def isecs(featureList):