0

I am trying to compute the intersection points between a line and a closed curve (stored in a file)

I tried to adapt this

import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import bisect
import numpy as np
import csv
from scipy.interpolate import interp1d
from scipy.optimize import bisect
#reading curve from file
r = ...
z = ...


# vertical line
x = 1.885, 1.885
y = [-3,3]



z_ves = interp1d(r_ves, z_ves, fill_value="extrapolate")
plt.figure(num=None, figsize=(10, 6), dpi=100, facecolor='w', edgecolor='k')
plt.plot(r_ves, z_ves(r_ves))
plt.axvline(x=1.885,ymin=-3,ymax=3)


# #use interp1d to get interpolated points
y = interp1d(x, y, fill_value="extrapolate")
# stress = interp1d(strain, stress)
# #set starting points
x1 = max(x[0], r_ves[0])
x2 = min(x[-1], r_ves[-1])
max_err = .01
# #create function
f = lambda x : z_ves(x) - y(x)
#find x1 where f(x1) = 0
x1 = bisect(f, x1, x2, xtol = .001)
y1 = z_ves(x1)
#
plt.figure(num=None, figsize=(10, 6), dpi=100, facecolor='w', edgecolor='k')
plt.plot(r_ves, z_ves(r_ves))
plt.plot(x, y(x))
plt.scatter(x1, y1)
plt.show()

I get a Runtime Warning as I think the line is vertical.

RuntimeWarning: invalid value encountered in multiply y_new = slope*(x_new - x_lo)[:, None] + y_lo

I know there are two intersections (see picture) enter image description here

how do I get the intersections?

bruvio
  • 853
  • 1
  • 9
  • 30
  • Just because it might get you better answers: this is not a curve. At all. It's just a polygon, and "python line polygon find intersection" in google might get you an answer in minutes. – Mike 'Pomax' Kamermans Apr 08 '20 at 01:25

2 Answers2

0

Function y should interpolate the vertical line, but you already know its equation: x=1.885. Therefore, plugging given x value into the equation of the curve should give you the equation you need to solve to get y-coordinates of the points of intersection. z_ves(1.885) should give you the desired result point of intersection.

Kate Melnykova
  • 1,863
  • 1
  • 5
  • 17
  • z_ves(1.885) does not return the two values.. returns one values and it's ~0... it should be around -0.25 and 0.85 – bruvio Apr 07 '20 at 18:48
  • Based on the docs, `interp1d` is not tailored for curves that are not functions. The docs suggest using parametric spline (https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) – Kate Melnykova Apr 08 '20 at 01:44
0

I solved with using shapely

from shapely.geometry import LineString,Polygon
LOS1 = LineString([(1.885,-5), (1.885,3)])
VesselCoordTuple=list(zip(r_ves, z_ves))
polygonVesselBound = Polygon(VesselCoordTuple)
x1 = polygonVesselBound.intersection(LOS1)
r1 = x1.xy[0][0]
z1 = x1.xy[1][0]
r2 = x1.xy[0][1]
z2 = x1.xy[1][1]
plt.figure(1, figsize=SIZE, dpi=90) #1, figsize=(10, 4), dpi=180)
plt.plot(r_ves, z_ves)
plt.plot(x1.xy[0][0],x1.xy[1][0], color='r',linestyle=' ',marker='x',markersize=5)
plt.plot(x1.xy[0][1],x1.xy[1][1], color='r',linestyle=' ',marker='x',markersize=5)
bruvio
  • 853
  • 1
  • 9
  • 30