I would like to make a graph from a shapefile, where edges would express spatial relations between polygons according to DE-9IM

For that, I use fiona + shapely

My code is working but it is extremely slow because I have to check all relation for each couples

import fiona
from shapely.geometry import shape, mapping

polys = fiona.open("....shp")

def TPP(A, B):
    rel = A.relate(B)
    return rel[0] == '2' and rel[5] in '01' and rel[6] == 'F'

import itertools

# that makes 124750 couples
for pol1,pol2 in  itertools.combinations(itertools.islice(polys, 500), 2):

    # equals
    if shape(pol1['geometry']).intersects(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " equals ",\

    # touches
    if shape(pol1['geometry']).touches(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " touches ",\

    # overlaps
    if shape(pol1['geometry']).overlaps(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " overlaps ",\

    # within  a.contains(b) == b.within(a)
    if shape(pol1['geometry']).within(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " within ",\

    # TPP
    if TPP(shape(pol1['geometry']), shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " TPP ",\

    # disjoint
    #if shape(pol1['geometry']).disjoint(shape(pol2['geometry'])):
    #    print "poly ", pol1["properties"]["ID_ILOT"], " disjoint ",\
    #    pol2["properties"]["ID_ILOT"]

Maybe an other approach would be to use 3×3 intersection matrix with relate, which should avoid unecessary calculation.

But even with only one operation (eg intersect) it takes plenty of time for 500 polygons, and my goal is to use this method on 10000+ polygons.

Does anyone have a better approach ?


By using tqdm, it tooks me 10 mins to perform the calculation over 124750 couples.

100%|██████████████████████████████████| 124750/124750 [09:23<00:00, 221.26it/s]

by processing the matrix only (more or less the same duration)

100%|██████████████████████████████████| 124750/124750 [09:27<00:00, 219.93it/s]

(without any other task alongside)

If I use 10000 polygons instead, I will have 49995000 couples. I dont think that my method is viable for such number


I figured a way to do it, I use a R-tree as spatial index to compute spatial relation on possible couple (geopanda have a nice way to do it) I used this tutoriel

  • 1
    Please Edit the question to quantify "extremely slow" and "plenty time" with actual duration, in seconds. – Vince Jan 19 at 18:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.