How to create an infinite line from two given points to get intersection with other geometry objects in Shapely?

user_10071179 0
Viewed 1456

I am able to get a line segment joining two points using the following function:

line = shapely.geomtery.LineString([[0,0],[0,1]])

This code creates a line segment between the two coordinates. How can I modify this to get a line of infinite length to pass through these points instead?

I need to find points of intersection between a line joining 2 point and the given quadrilateral specified by coordinates of its 4 corners.

shapely polygon geometry python
1 Answers

Shapely doesn't support infinite lines but you could make an approximation by extending the given line over the bounding box covering other geometries:

from shapely.geometry import box, LineString, Point, Polygon

# polygon or any other geometry object/collections of objects
polygon = Polygon([(-100, -100), (60, -60), (100, 100), (-60, 60)])
# your original line (could be located anywhere - inside/outside/crossing the polygon)
line = LineString([(-130, -130), (-116, -110)])

Here is an example of a line and a polygon that we want to split/get intersection with:

enter image description here

Next, from the coordinates of the bounding box and the coordinates of the endpoints of the given line, we can calculate coordinates for a new extended line. This code should be self-explanatory:

minx, miny, maxx, maxy = polygon.bounds
bounding_box = box(minx, miny, maxx, maxy)
a, b = line.boundary
if a.x == b.x:  # vertical line
    extended_line = LineString([(a.x, miny), (a.x, maxy)])
elif a.y == b.y:  # horizonthal line
    extended_line = LineString([(minx, a.y), (maxx, a.y)])
else:
    # linear equation: y = k*x + m
    k = (b.y - a.y) / (b.x - a.x)
    m = a.y - k * a.x
    y0 = k * minx + m
    y1 = k * maxx + m
    x0 = (miny - m) / k
    x1 = (maxy - m) / k
    points_on_boundary_lines = [Point(minx, y0), Point(maxx, y1), 
                                Point(x0, miny), Point(x1, maxy)]
    points_sorted_by_distance = sorted(points_on_boundary_lines, key=bounding_box.distance)
    extended_line = LineString(points_sorted_by_distance[:2])

And we can see that the newly created line now crosses the original polygon:

enter image description here

From here, getting the intersection is trivial:

intersection = extended_line.intersection(polygon)

or if you want to split:

from shapely.ops import split
parts = split(polygon, extended_line)
user_78514700
Related