Find intersection between line and shape file with python

I’m trying to clip lines based on a shape-file in python. I have a script that works, but it is very slow. Are there a faster way to do this? A prerequisite are that I have to do it with python.

from shapely.geometry import LineString, shape, Point
from shapely.affinity import rotate
from itertools import chain
from shapely.ops import cascaded_union
from fiona import open
import numpy as np

# open the shp-file with land geometry
shoreline = open('land_boundary.shp')
shapes = [shape(f['geometry']) for f in shoreline]
mergedshorelines = cascaded_union(shapes)

# create an arbitrary line
x,y = 696346,6601295
x_end, y_end = 746345,6601295
line1 = LineString([(x, y), (x_end, y_end)])

# Creates lines from arbitrary line with 1 degree step
radii = [rotate(line1, deg, (x,y)) for deg in chain(np.mod(np.arange(0,360,1),360))]
mergedradii = cascaded_union(radii)

# the intersection between the two multilines is computed and the intersection point with the smallest distance is choosen
point =  Point(x,y)
points_ok = []

#----THIS IS THE SLOW PART-------
# clip lines with land boundary 
for line in mergedradii:
    if line.intersects(mergedshorelines):
        if line.intersection(mergedshorelines).type == "MultiPoint":
          # multiple points: select the point with the minimum distance
          a = {}
          for pt in line.intersection(mergedshorelines):
              a[point.distance(pt)] = pt
          points_ok.append(a[min(a.keys())])
        if line.intersection(mergedshorelines).type == "Point":
          # ok, only one intersection
          points_ok.append(line.intersection(mergedshorelines))
Shoreline_points = cascaded_union(points_ok) # coordinates of all intersection.

Appreciate any input! Cheers! /Björn

Answer

You can save a lot of time by calling intersection just once and saving the result, instead of calling it multiple times. I call the intersection method once at the top of the loop and save it as a variable, then refer to the variable in the logic, instead of running it again.

You can also save time by skipping the intersects check entirely, because if a line doesn’t intersect with mergedshorelines, it won’t be a MultiPoint or Point object, it will be a LineString.

You can also optimize the part of the script that looks for the closest point in a MultiPoint. When trying to find the closest point to the origin of the line, you only need to check the first and last point in a MultiPoint, because the points are sorted. So either the first or last point will be closest to the origin.

for line in mergedradii:
    intersection = line.intersection(mergedshorelines)
    if intersection.type == "MultiPoint":
        # multiple points: select the point with the minimum distance
        first_pt = intersection[0]
        last_pt = intersection[-1]
        if point.distance(first_pt) < point.distance(last_pt):
            points_ok.append(first_pt)
        else:
            points_ok.append(last_pt)
    elif intersection.type == "Point":
        # ok, only one intersection
        points_ok.append(intersection)
Shoreline_points = cascaded_union(points_ok) # coordinates of all intersection.

On my computer, this script runs in about 40 seconds, while the original script took around 190 seconds.