TransWikia.com

Error when calculating final polygon after subtracting all intersecting polygons using shapely

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:

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.

One Answer

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

Correct answer by gene on March 18, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP