Search code examples
pythonpolygongeopandasintersection

Highlight polygons that share a side/edge


How do I create a... ---I don't know; list, dict of all the polygons, in a GeoDataFrame, that share a side/edge. The polygons will intersect but never cross.

import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt

polys = gpd.GeoSeries([Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
                     Polygon([(0,2), (2,2), (2,4), (0,4)]),
                     Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
                     Polygon([(3,3), (5,3), (5,5), (3,5)])])

fp = gpd.GeoDataFrame({'geometry': polys, 'name': ['a', 'b', 'c', 'd'],
                      'grnd': [25, 25, 25, 25],
                      'rf': [29, 35, 26, 31]})

fig, ax = plt.subplots(figsize=(5, 5))
fp.plot(ax=ax, alpha=0.3, cmap='tab10', edgecolor='k',)
fp.apply(lambda x: ax.annotate(text=x['name'], xy=x.geometry.centroid.coords[0], ha='center'), axis=1)
plt.show()

enter image description here

for i, row in fp.iterrows():
    oring = list(row.geometry.exterior.coords)#, row['ground_height']
    if row.geometry.exterior.is_ccw == False:
        #-- to get proper orientation of the normals
        oring.reverse() 
    for (j, v) in enumerate(oring[:-1]):
        print(oring[j][0], oring[j][1], oring[j+1][0], oring[j+1][1], row['name'])

Expected result:

0.0 0.0 2.0 0.0 a
2.0 0.0 2.0 1.5 a c
2.0 1.5 2.0 2.0 a 
2.0 2.0 0.0 2.0 a b
0.0 2.0 0.0 0.0 a
0.0 2.0 2.0 2.0 b a
2.0 2.0 2.0 4.0 b 
2.0 4.0 0.0 4.0 b
0.0 4.0 0.0 2.0 b... and so on

Solution

  • Looking at the expected result again, I would first create a series that includes all of the line segments and then check if each line segment both touches the polygon and that the length of intersection is greater than 0:

    import geopandas as gpd
    from shapely.geometry import Polygon, LineString
    import matplotlib.pyplot as plt
    
    polys = gpd.GeoSeries([
        Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
        Polygon([(0,2), (2,2), (2,4), (0,4)]),
        Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
        Polygon([(3,3), (5,3), (5,5), (3,5)])
    ])
    
    fp = gpd.GeoDataFrame({
        'geometry': polys, 
        'name': ['a', 'b', 'c', 'd'],
        'grnd': [25, 25, 25, 25],
        'rf': [29, 35, 26, 31]
    })
    
    # Create series of all the line segments
    lines = fp.geometry.apply(lambda x: list(map(
        LineString, 
        zip(x.boundary.coords[:-1], x.boundary.coords[1:]))
    )).explode()
    
    result = {
        str(line): list(fp.loc[
            (fp.geometry.touches(line)) # line touches the polygon
            & (fp.geometry.intersection(line).length > 0), # And the intersection is more than just a point
            'name'
        ].values)
        for line in lines
    }
    

    output:

    {'LINESTRING (0 0, 2 0)': ['a'],
     'LINESTRING (2 0, 2 1.5)': ['a', 'c'],
     'LINESTRING (2 1.5, 2 2)': ['a'],
     'LINESTRING (2 2, 0 2)': ['a', 'b'],
     'LINESTRING (0 2, 0 0)': ['a'],
     'LINESTRING (0 2, 2 2)': ['a', 'b'],
     'LINESTRING (2 2, 2 4)': ['b'],
     'LINESTRING (2 4, 0 4)': ['b'],
     'LINESTRING (0 4, 0 2)': ['b'],
     'LINESTRING (2 0, 5 0)': ['c'],
     'LINESTRING (5 0, 5 1.5)': ['c'],
     'LINESTRING (5 1.5, 2 1.5)': ['c'],
     'LINESTRING (2 1.5, 2 0)': ['a', 'c'],
     'LINESTRING (3 3, 5 3)': ['d'],
     'LINESTRING (5 3, 5 5)': ['d'],
     'LINESTRING (5 5, 3 5)': ['d'],
     'LINESTRING (3 5, 3 3)': ['d']}