Search code examples
pythonpolygoncomputational-geometryshapelyalpha-shape

get coordinate of points where a line crosses a polygon

I would like to find the points where a line intersects with a polygon. I obtain this polygon using a concave outline calculation from this thread.

import alphashape
from shapely.geometry import LineString
import matplotlib.pyplot as plt
from descartes import PolygonPatch

points = [(17, 158),(15, 135),(38, 183),(43, 19),(93, 88),(96, 140),(149, 163),(128, 248),(216, 265),(248, 210),(223, 167),(256, 151),(331, 214),(340, 187),(316, 53),(298, 35),(182, 0),(121, 42)]
points = np.array(points)

alpha = 0.99 * alphashape.optimizealpha(points)
hull = alphashape.alphashape(points, alpha)
hull_pts = hull.exterior.coords.xy
path = PolygonPatch(hull, fill=False, color='green')
print(path.contains_point([128,248]))

fig, ax = plt.subplots()
ax.scatter(hull_pts[0], hull_pts[1], color='red')
ax.scatter(points[:,0], points[:,1], color='red')

p = np.array([[350, 100],[0, 100]])
ax.plot(p[:, 0], p[:, 1], color='blue')
ax.add_patch(path)

enter image description here

so far I tried defining a line with:

l = LineString(p)
inters = l.intersection(hull)

but inters.xy returns a NotImplemented error, so I am not sure how to obtain the coordinates of points where the line crosses the polygon.

Solution

  • The intersection returns a MultilineString, which is a fancy word for list of LineStrings. We can retrieve the coordinates then from each Linestring, e.g.:

    import alphashape
from shapely.geometry import LineString
import matplotlib.pyplot as plt
import numpy as np

#replicating your example
fig, ax = plt.subplots()

line_xy = [[350, 100],[0, 120]]
points = np.asarray([(17, 158),(15, 135),(38, 183),(43, 19),(93, 88),(96, 140),(149, 163),(128, 248),(216, 265),(248, 210),(223, 167),(256, 151),(331, 214),(340, 187),(316, 53),(298, 35),(182, 0),(121, 42)])

alpha = 0.99 * alphashape.optimizealpha(points)
hull = alphashape.alphashape(points, alpha)
hull_pts = hull.exterior.coords.xy

p = LineString(line_xy)

ax.plot(*hull_pts, c="green")
ax.scatter(points[:,0], points[:,1], marker="o", color="red")
ax.scatter(*hull_pts, marker="s", color="red")

ax.plot(*p.coords.xy, color='blue')

#retrieving intersection 
inters = hull.intersection(p)

#checking for object type to retrieve all intersection coordinates
if inters.type == "LineString":
    coords = np.asarray([inters.coords.xy])
elif inters.type == "MultiLineString":
    coords = np.asarray([l.coords.xy for l in inters.geoms])
    
#reshaping array point coordinates into a form that does not make my head hurt
coords = coords.transpose(1, 0, 2).reshape(2, -1)
print(coords)

plt.show()

    with coords[:, i] returning the x-y values for intersection point i.

    Sample output:

    [[324.67707894 234.24811338 176.4217078   18.88111888]
 [101.44702406 106.61439352 109.91875955 118.92107892]]

    enter image description here

    Weirdly, shapely considers a line point like 300, 100 within the hull as an intersection point. Strictly speaking, one had to check that none of the identified points lies within the hull polygon.

    And alphashape (1.3.1 used here) should update their routine because alphashape.alphashape(points, alpha) generates the error message ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the geoms property to access the constituent parts of a multi-part geometry.