2

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?

3
  • 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, 2019 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, 2019 at 10:03
  • 1
    Perhaps when you said "having a single feature" you meant "having a single field" ? Feb 6, 2019 at 15:04

1 Answer 1

3

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])
1
  • this is great! how can you dump this information into a df and plot it? I want to know which polygons were intersecting and calculate their area.
    – delViento
    Sep 7, 2020 at 10:36

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