1

This question is related to here.

How do I create a dict of all the polygons [including MultiPolygons], in a GeoDataFrame, that share a side/edge. The polygons will intersect but never cross.

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)]),
                       MultiPolygon([Polygon([(6,1), (8, 1), (8, 3), (6, 3)], 
                                     [[(6.5, 1.5), (7.5, 1.5), (7.5, 2.5), (6.5, 2.5)][::-1]]
                                            )])])

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

fig, ax = plt.subplots(figsize=(6, 6))
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

I've tried this:

# 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:])))
                         if isinstance(x, MultiPolygon)
                         else list(chain(*list(list(map(
                             LineString, 
    zip(x.boundary.coords[:-1], x.boundary.coords[1:]))) 
                                               for poly in x.geoms)))
                         )).explode()

result = {
        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
            'rf'].values)
        for line in lines}
    

which yields: AttributeError: 'Polygon' object has no attribute 'geoms'

Expected result; something like:

{'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'],...}
arkriger
  • 207
  • 2
  • 7
  • FYI: Your `MultiPolygon` is a `single` polygon with a hole. – swatchai Mar 24 '23 at 22:37
  • Debugging is harder when you optimize your code before it gets correct answer. – swatchai Mar 25 '23 at 01:20
  • @swatchai ---You are a helpful person. I am a data capturer. Not a computer scientist. Your expertise and advice, helps me improve and, is appreciated. The post in the [link](https://stackoverflow.com/questions/75731055/highlight-polygons-that-share-a-side-edge) has the working code (correct result) ---with normal `Polygon`. This post includes `MultiPolygon` ---with my failed attempt. – arkriger Mar 25 '23 at 20:10

1 Answers1

0

I think you have some misplaced parentheses and are using a wrong variable name in a few places. It might be easier to read if you extract the funtion to return a list of the exterior lines of the polygon. Give this a try:

def extract_exterior_lines(polygon):
    return list(map(
        LineString, 
        zip(polygon.exterior.coords[:-1], polygon.exterior.coords[1:])
    ))

lines = fp.geometry.apply(lambda x: (
    extract_exterior_lines(x)
    if isinstance(x, Polygon)
    else list(chain(*(extract_exterior_lines(poly) for poly in x.geoms)))
)).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
        'rf'
    ].values)
    for line in lines
}
Zach Flanders
  • 1,224
  • 1
  • 7
  • 10