Geographic Information Systems Asked on March 18, 2021
This question is a follow-on question from Calculating final polygon after subtracting all intersecting polygons using shapely which was answered.
I have a list of centroids that I buffer to 5000 meters radius. I need to take each resulting circle and run it against the remaining centroids(also buffered to 5000 m) and find the final shape after subtracting all intersections.
In other words, I need the circle shape A to iterate over the dataset, find all circle shape(s) that intersect with it (B and C) and then create a column with a geocoordinates of the polygon, that is resulting after all intersecting circle areas are subtracted (see image below) and if it does not intersect with any other circle(s), just add a geocoordinate for the entire circle:
Here is the code I have:
df.head()
# Output:
# cell_centroid pop id
#0 POINT(21.25602 55.903539) 9004.0 0
#1 POINT(21.119156 55.538406) 517.0 1
#2 POINT(21.249535 55.890204) 10369.0 2
#3 POINT(21.359934 55.904164) 800.0 3
#4 POINT(24.433042 56.104589) 800.0 4
df['cell_centroid'] = df['cell_centroid'].apply(wkt.loads)
gdf = gpd.GeoDataFrame(df, geometry="cell_centroid", crs = 4326).to_crs(crs=3857)
gdf["buffer"] = gdf.buffer(5000)
gdf["difference"] = None
gdf = gdf.set_geometry("buffer")
for i, row in gdf.iterrows():
others = gdf[gdf["id"] != row["id"]]
diff = row["buffer"].difference(others.cascaded_union)
gdf.loc[i, "difference"] = diff
gdf = gdf.set_geometry("difference")
The error I am getting now is:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-122-ad757ced6ebf> in <module>
3 diff = row["buffer"].difference(others.cascaded_union)
4 # gdf.loc[i, "difference"] = gpd.GeoSeries(diff)
----> 5 gdf.loc[i, "difference"] = diff
...
ValueError: Must have equal len keys and value when setting with an iterable
What’s the right way to assign the resulting polygon (diff) to the GeoDataFrame?
I have found GeoPandas setting geometry for a row with .loc and Multipolygons for a similar issue, but its solution does not work in my case.
Deniz is right, the problem is "difference" for field name because difference is:
a) a method object of Shapely used by GeoPandas (object.difference(other)
b) an index object of Pandas used by GeoPandas (pandas.Index.difference)
import geopandas as gpd
import matplotlib.pyplot as plt
import descartes
gdf = gpd.read_file("circles.shp")
# plot circles
fig, ax = plt.subplots(figsize=(6,6))
for i in gdf.index:
ax.plot(*gdf.geometry[i].exterior.xy)
ax.add_patch(descartes.PolygonPatch(gdf.geometry[i], fc='blue', alpha=0.2))
ax.axis('equal');
With a field name different from "difference" and using unary_union
instead of cascaded_union
(see shapely.ops.cascaded_union(geoms)¶)
gdf["differ"] = None
for i, row in tt.iterrows():
others = gdf[gdf["id"] != row["id"]]
diff = row["geometry"].difference(others.unary_union)
gdf.loc[i, "differ"] = diff
fig, ax = plt.subplots(figsize=(6,6))
for i in tt.index:
ax.plot(*gdf.differ[i].exterior.xy)
ax.add_patch(descartes.PolygonPatch(gdf.differ[i], fc='blue', alpha=0.2))
ax.axis('equal');
You can also use itertuples()
for row in gdf.itertuples(index=True):
others = gdf[gdf["id"] != getattr(row, "id")]
diff = getattr(row,"buffer").difference(others.unary_union)
gdf.loc[row.Index, 'differ'] = diff
Correct answer by gene on March 18, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP