问题描述
当我在 this 页面上看到解决方案时，我正在寻找解决方案以找到 2 条线之间的交点。我决定使用 shapely.geometry 解决方案，因为它似乎是一个优雅的解决方案。这是适用于我的代码的最小工作示例：
import numpy as np
from shapely.geometry import Linestring
import array as arr
import matplotlib.pyplot as plt
from figure_size_def import set_size
plt.style.use('/msn2e2/avenkat/thesis/plots/plotstyle.mplstyle')
xcont = np.array([0.2204,0.7103,1.2,1.788,2.376,2.866,3.355,3.943,4.433,4.923,5.511,6.099,6.588,7.176,7.666,8.254,8.842,9.429,9.919])
ycont_mean = np.array([1.,0.96435101,0.89579525,0.77998172,0.65804388,0.54981718,0.43308958,0.33656307,0.25393053,0.19341865,0.14469835,0.10566728,0.08189214,0.06001828,0.04806216,0.03705667,0.02989945,0.0247532,0.0202011])
ycont_tot = np.array([6.36678946e-05,1.03060712e-04,1.47423350e-04,1.97053946e-04,2.05546123e-04,2.10758605e-04,2.20563483e-04,1.99077460e-04,1.73477086e-04,1.48655253e-04,1.16881677e-04,8.85778352e-05,7.04365528e-05,5.18005196e-05,3.98887296e-05,3.18437309e-05,2.67239221e-05,2.36262960e-05,2.26682510e-05])
deconv_cont_mean = np.array([1.,0.12318841,0.14629156,0.15575448,0.09488491,0.0683717,0.04237852,0.03279625,0.02100597,0.01861893,0.0126087,0.00785422,0.00994032,0.00489258,0.00767093,0.00436061,0.00456948,0.00422933,0.00298721])
a = []
b = []
a = np.isnan(ycont_mean)
b = np.isnan(ycont_tot)
for i in range(len(xcont)):
if a[i] == True: ycont_mean[i] = 0
if b[i] == True: ycont_tot[i] = 0
y = np.full(shape = len(xcont),fill_value = ycont_mean[0]/2.718)
line1 = Linestring(np.column_stack((xcont,ycont_mean[0:len(xcont)])))
line2 = Linestring(np.column_stack((xcont,ycont_mean[0:len(xcont)] - ycont_tot[0:len(xcont)])))
line3 = Linestring(np.column_stack((xcont,ycont_mean[0:len(xcont)] + ycont_tot[0:len(xcont)])))
line4 = Linestring(np.column_stack((xcont,y)))
int1 = line1.intersection(line4)
int2 = line2.intersection(line4)
int3 = line3.intersection(line4)
if int1.geom_type == 'MultiPoint': xc1,_ = Linestring(int1).xy
elif int1.geom_type == 'Point': xc1,_ = int1.xy
elif int1.geom_type == 'Linestring': xc1,_ = arr.array('d',[0]),arr.array('d',[0])
if int2.geom_type == 'MultiPoint': xc2,_ = Linestring(int2).xy
elif int2.geom_type == 'Point': xc2,_ = int2.xy
elif int2.geom_type == 'Linestring': xc2,[0])
if int3.geom_type == 'MultiPoint': xc3,_ = Linestring(int3).xy
elif int3.geom_type == 'Point': xc3,_ = int3.xy
elif int3.geom_type == 'Linestring': xc3,[0])
xc1 = xc1[0]
xc2 = xc2[0]
xc3 = xc3[0]
a = []
a = np.isnan(deconv_cont_mean)
for i in range(len(xcont)):
if a[i] == True: deconv_cont_mean[i] = 0
y = np.full(shape = len(xcont),fill_value = deconv_cont_mean[0]/2.718)
line1 = Linestring(np.column_stack((xcont,deconv_cont_mean[0:len(xcont)])))
line4 = Linestring(np.column_stack((xcont,y)))
int1 = line1.intersection(line4)
if int1.geom_type == 'MultiPoint': xcd1,_ = Linestring(int1).xy
elif int1.geom_type == 'Point': xcd1,_ = int1.xy
elif int1.geom_type == 'Linestring': xcd1,[0])
xcd1 = xcd1[0]
#Choose Maximum radius limit
rad_min = (np.min(xcont)) - 0.02
rad_max = (np.max(xcont)) + 2.00
ymin,ymax = 0.001,2.0
# Plot figures
fig,ax = plt.subplots(1,2,figsize = set_size(fraction = 1,subplots = (1,2)),dpi = 300,constrained_layout = True,sharex = False,sharey = False)
xaxis = '$\mathrm{Radius (kpc)}$'
yaxis1 = '$\Sigma_{\mathrm{cont}}\,\mathrm{[Jy\,kpc^{-2}]}$'
yaxis2 = '$\Sigma_{\mathrm{line}}\,kms^{-1}\,kpc^{-2}]}$'
# figure 1 is the comparision between ellipse profiles
axis = ax[0]
axis.set_title('Ellipse Profiles')
axis.tick_params(labelleft = True,labelright = False,labeltop = False,labelbottom = True)
axis.set_xlim([rad_min,rad_max])
axis.set_ylim([ymin,ymax])
axis.set_xscale("log")
axis.set_yscale("log")
# axis.set_xlabel(xaxis)
axis.set_ylabel(yaxis1)
# axis2.set_ylabel(yaxis2)
axis.plot(xcont,ycont_mean[0:len(xcont)],label = 'Continuum') #yerr = ycont_tot[0:len(xcont)],axis.fill_between(xcont,ycont_mean[0:len(xcont)]-ycont_tot[0:len(xcont)],ycont_mean[0:len(xcont)]+ycont_tot[0:len(xcont)],alpha = 0.3)
if xc1 != 0:
axis.axvline(xc1,color = '#1f77b4',marker = '',ls = '--')
# axis.fill_betweenx(xc1,xc1-xc2,xc3-xc1)
axis.axvline(xc2,ls = ':')
axis.axvline(xc3,ls = ':')
axis.axhline(y = ycont_mean[0]/2.718,color = 'black',ls = '--',label = '1/e Intensity')
axis.legend()
# figure 2 is the comparison between manually deconvolved profiles
axis = ax[1]
axis2 = axis.twinx()
axis.set_title('Deconvolved Profiles')
axis.tick_params(labelleft = False,labelbottom = True)
axis2.tick_params(labelleft = False,labelright = True,labelbottom = False)
axis2.grid(False)
axis.set_xlim([rad_min,ymax])
axis2.set_ylim([ymin,ymax])
axis.set_xscale("log")
axis.set_yscale("log")
axis2.set_xscale("log")
axis2.set_yscale("log")
# axis.set_xlabel(xaxis)
# axis.set_ylabel(yaxis1)
axis2.set_ylabel(yaxis2)
axis.plot(xcont,deconv_cont_mean[0:len(xcont)])
axis.axhline(y = ycont_mean[0]/2.718,label = '1/e Intensity')
axis.axvline(xcd1,ls = '--')
这产生的结果是： result
为什么第一个路口可以正常工作，但第二个路口失败得这么惨？如果这是一个已经纠正的错误，请指出我的方向，因为我正在努力了解这里发生的事情。提前致谢！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）