Hi,

I am trying to create a polygon to outline a set of roads based on a sample of coordinates. I can do this in Matlab with the function called `polybuffer’. I have been struggling to find an equivalent function in the Shapley package. Here is a simple example to illustrate the differences I have encountered:

As you can see, the code from Matlab looks more like roads whereas the Python code creates a convex polygon.

The only other requirement I have is that I need to be able to add a buffer and test whether a given point is within the buffer (e.g., the filled area in the matlabplot). I looked into alpha shapes, but I did not find a way to test the intersection of a point.

Is there a Python solution?

Corresponding code:

Python

``````
from shapely.geometry import Point, LineString, Polygon
from shapely.ops import unary_union
from matplotlib import pyplot as plt

buffer = .05

poly1 = Polygon([(0,0), (2,0), (2,2), (0,2)]).buffer(buffer)

poly2 = Polygon([(2,2), (4,2), (4,4), (2,4)]).buffer(buffer)

poly3 = Polygon([(1,1), (3,1), (3,3), (1,3)]).buffer(buffer)

poly4 = Polygon([(3,3), (5,3), (5,5), (3,5)]).buffer(buffer)

polys = [poly1, poly2, poly3, poly4]

polygon_union = unary_union(polys)

fig, axs = plt.subplots()

axs.set_aspect('equal', 'datalim')

xs, ys = polygon_union.exterior.xy

axs.plot(xs, ys, alpha=0.5)

plt.plot(xs, ys)

``````

Matlab

``````
buffer = .05

coords = [0 0; 2 0; 2 2; 0 2]

polygon(1) = polybuffer(coords, 'lines', buffer)

coords = [2 2; 4 2; 4 4; 2 4]

polygon(2) = polybuffer(coords, 'lines', buffer)

coords = [1 1; 3 1; 3 3; 1 3]

polygon(3) = polybuffer(coords, 'lines', buffer)

coords = [3 3; 5 3; 5 5; 3 5]

polygon(4) = polybuffer(coords, 'lines', buffer)

union_polygon = union(polygon)

plot(union_polygon)

``````

The issue here is that you are buffering a Polygon, not a line. For example:

``````Polygon([(0,0), (2,0), (2,2), (0,2)])
``````

is a four-vertex polygon – i.e. a square. What you want is a LineString, which you then buffer to make a polygon.

This gets you very close:

``````from matplotlib import pyplot as plt

from shapely.geometry import Polygon, LineString
from shapely.ops import unary_union

buffer = .05

poly1 = LineString([(0,0), (2,0), (2,2), (0,2)]).buffer(buffer)

poly2 = LineString([(2,2), (4,2), (4,4), (2,4)]).buffer(buffer)

poly3 = LineString([(1,1), (3,1), (3,3), (1,3)]).buffer(buffer)

poly4 = LineString([(3,3), (5,3), (5,5), (3,5)]).buffer(buffer)

polys = [poly1, poly2, poly3, poly4]

polygon_union = unary_union(polys)

xs, ys = polygon_union.exterior.xy

fig, axs = plt.subplots()
axs.set_aspect('equal', 'datalim')
axs.plot(xs, ys, alpha=0.5)

plt.plot(xs, ys)
plt.show()
``````

I"ll leave it as an exercise for the reader to figure out why it’s not exactly what you are looking for BTW – I need a couple extra imports – do try to post fully functional code.

1 Like

Thank you very much for you help. Looking back, I can see why `LineString` is used instead of `Polygon`. What confused me was that line was passed to `polybuffer` in Matlab. Perhaps it would have been clearer to have a separate function as in Python. Nonetheless, I am grateful for your solution.

My apologies for the missing import. I’m not sure how I missed it, but I corrected the original post for the benefit of future readers.

To follow up with the plot, the problem is that the original code only extracts the exterior lines. The code below is reworked to include the interior lines:

``````from shapely.geometry import Point, LineString, Polygon
from shapely.ops import unary_union
from matplotlib import pyplot as plt

buffer = .05
poly1 = LineString([(0,0), (2,0), (2,2), (0,2)]).buffer(buffer)
poly2 = LineString([(2,2), (4,2), (4,4), (2,4)]).buffer(buffer)
poly3 = LineString([(1,1), (3,1), (3,3), (1,3)]).buffer(buffer)
poly4 = LineString([(3,3), (5,3), (5,5), (3,5)]).buffer(buffer)
polys = [poly1, poly2, poly3, poly4]
polygon_union = unary_union(polys)

xs, ys = polygon_union.exterior.xy

fig, axs = plt.subplots()
axs.set_aspect('equal', 'datalim')
axs.plot(xs, ys, alpha=0.5)
plt.plot(xs, ys, color="blue")
for inner in polygon_union.interiors:
xi, yi = zip(*inner.coords[:])
axs.plot(xi, yi, color="blue")
plt.show()
``````

Now that I have the correct shape, the intersection tests work as expected:

``````# correctly returns false
point = Point(2.5, 1.5)
polygon_union.contains(point)
``````
``````# correctly returns true
point = Point(2, 1)
polygon_union.contains(point)
``````