get coordinate of points where a line crosses a polygon

Question

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)

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.

Answer 1

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]]

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.

Answer 2

An edge crosses an horizontal line when the ordinates of its endpoints (let a and b) straddle the ordinate of the line (let Y). Hence an intersection test is very simple. To get the intersection point itself, use

Xi = Xa + (Y - Ya) (Xb - Xa) / (Yb - Ya)
Yi = Y

What to do when Y = Ya = Yb is application dependent.

---

If the line is vertical, just swap the roles of X and Y.

If the line is oblique, you can use a rotation that makes the line horizontal. This is conveniently expressed in complex numbers by the transformation

Z = z e^(-iΘ)

though you don't need complex numbers to implement this, and the way to obtain the rotation coefficients depends on how the line is specified.

After computing the intersection, you counter-rotate it.