2

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:

enter image description here

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.

2
  • Try to use at instead of loc. Jan 14 at 9:30
  • I think you'd better to change you field name and not using "difference" or any other function or property of geodataframe/geoseries.
    – Deniz
    Jan 14 at 11:30
1

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');

enter image description here

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');

enter image description here

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
1
  • That worked, thank you! The differ will have the polygon that defines all three final shapes, correct? How would I put each a shape difference for each circle into a separate column as a separate polygon? Jan 14 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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