Reputation: 244
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()
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'],...}
Upvotes: 1
Views: 119
Reputation: 1949
import geopandas as gpd
from shapely.geometry import LineString, Polygon, MultiPolygon
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]]
)])])
df = gpd.GeoDataFrame({'geometry': polys, 'name': ['a', 'b', 'c', 'd', 'e'],
'grnd': [25, 25, 25, 25, 25],
'rf': [29, 35, 26, 31, 28]})
#Save the geometry in a column, then spatial join the df to itself by touching predicate.
# The result will for example contain one row containing both the a and c polygon geometries.
# When therse are intersected the result will be the line LINESTRING (2 0, 2 1.5)
df["geombackup"] = df.geometry
sj = df.sjoin(df, how="inner", predicate="touches")
#Create a dict of all shared sides
sj["intersection_wkt"] = sj.apply(lambda x: x.geombackup_left.intersection(
x.geombackup_right), axis=1).to_wkt()
sj["names"] = sj[["name_left","name_right"]].values.tolist()
dict_touches = sj.set_index("intersection_wkt")["names"].to_dict()
# {'LINESTRING (2 0, 2 1.5)': ['a', 'c'],
# 'LINESTRING (2 2, 0 2)': ['a', 'b'],
# 'LINESTRING (0 2, 2 2)': ['b', 'a'],
# 'LINESTRING (2 1.5, 2 0)': ['c', 'a']}
#For the polygons touching another:
#Create a dict of the sides that are not touching the other polygon.
# For example b's LINESTRING (2 4, 0 4)
df2 = df.loc[df.name.isin(sj.name_left)].copy() #Select the polygons that have a adjacent neighbor
def split_at_vertices(frame):
"""A function to split an input polygon geometry boundary line
at each vertex into individual segments and return them as a list of wkts"""
vertices = list(frame.geometry.boundary.coords)
line_list = [LineString([p1, p2]).wkt for p1,p2 in zip(vertices, vertices[1:])]
return line_list
df2["line_segments"] = df2.apply(lambda x: split_at_vertices(x), axis=1)
df2 = df2.explode(column="line_segments")
#Create the dict with line wkt as key and a list of the name as value
df2["name"] = df2["name"].apply(list)
dict_non_touching_sides = df2.set_index("line_segments").name.to_dict()
#Union the two dicts
full_dict = dict_touches|dict_non_touching_sides
# {'LINESTRING (2 0, 2 1.5)': ['a'],
# 'LINESTRING (2 2, 0 2)': ['a'],
# 'LINESTRING (0 2, 2 2)': ['b'],
# 'LINESTRING (2 1.5, 2 0)': ['c'],
# 'LINESTRING (0 0, 2 0)': ['a'],
# 'LINESTRING (2 1.5, 2 2)': ['a'],
# 'LINESTRING (0 2, 0 0)': ['a'],
# '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']}
Upvotes: 0
Reputation: 1314
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
}
Upvotes: 0