I have a mkl file A and a shapefile B containing polygons, each one of them having a single feature: the name associated to the polygon.

I am looking for the polygons of A which are overlapped by polygons of B using Python. In particular, I would like to know which are the best and fastest libraries to solve issues like this.

I have read this question and the answer propose to add the geometries into a PostGIS db or to use shapely. Do you have any reference or working example from where I can learn how to properly solve this task?

closed as off-topic by ahmadhanb, Erik, BERA, Jochen Schwarze, Vince Feb 6 at 14:47

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "When seeking help to debug/write/improve code always provide the desired behavior, a specific problem/error and the shortest code (as formatted text, not pictures) needed to reproduce it in the question body. Providing a clear problem statement and a code attempt helps others to help you." – ahmadhanb, Erik, BERA, Jochen Schwarze, Vince
If this question can be reworded to fit the rules in the help center, please edit the question.

  • You should take the tour. All coding question are expected to have code attempts included in the question or they will be put on hold and cant be answered. – BERA Feb 6 at 9:25
  • The point is that at the moment I don't even know how to start to tackle this issue. This evening I will try to add a sample code. – peppe Feb 6 at 10:03
  • 1
    Perhaps when you said "having a single feature" you meant "having a single field" ? – Kirk Kuykendall Feb 6 at 15:04

Once you have loaded both files as features collections (have a look at fiona for that), the brute force approach is to test for intersection (or overlap, there's a slight difference between the two terms, the shapely documentation explains it) of all against all. This can be done in a nested loop.

from shapely.geometry import shape

# List to collect pairs of intersecting features
fc_intersect = []

for featA in fc_A:
    for featB in fc_B:
        if shape(featA['geometry']).intersects(shape(featB['geometry'])):
            fc_intersect.append([featA, featB])

If you have a lot of features on both sides, this can get very slow. Speed can be optimized by using a spatial index; that's what the rtree package does.

from rtree import index
from shapely.geometry import shape

# List to collect pairs of intersecting features
fc_intersect = []

# Instantiate index class
idx = index.Index()
for i,featA in enumerate(fc_A):
    idx.insert(i, shape(featA['geometry']).bounds)

for featB in fc_B:
    # Test for potential intersection with each feature of the other feature collection
    for intersect_maybe in idx.intersection(shape(featB['geometry']).bounds):
        # Confirm intersection
        if shape(featB['geometry']).intersects(shape(fc_B[intersect_maybe]['geometry'])):
            fc_intersect.append([fc_A[intersect_maybe], featB])
  • Thanks a lot, very clear! :) – peppe Feb 6 at 16:08

Not the answer you're looking for? Browse other questions tagged or ask your own question.