Create polygon around roads

Python Help
1

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:

sidebyside
sidebyside633×1077 15.8 KB

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?

Thanks in advance for your help.

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)
2

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()

image
image1172×880 14.3 KB

I"ll leave it as an exercise for the reader to figure out why it’s not exactly what you are looking for :slight_smile:

BTW – I need a couple extra imports – do try to post fully functional code.

1 Like
3

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)