由两个散射点表面之间的交点定义的曲线，采样方式不同

Question

我有两组1__scatter_xyz.dat和2__scatter_xyz.dat散点。

这些点由 3 个坐标定义：x, y,z

1__scatter_xyz.dat: https://paste.ubuntu.com/25069931/

2__scatter_xyz.dat: https://paste.ubuntu.com/25069938/

这两组散点在一个区域内相交：

gnuplot>  splot "1__scatter_xyz.dat" using 3:1:2 with points lt 1 title "1", "2__scatter_xyz.dat" using 3:1:2 with points lt 1 lc 2 title "2"

gnuplot> set xlabel 'x'
gnuplot> set ylabel 'y'
gnuplot> set zlabel 'z'

的表面之间的交叉set 1与表面set 2将限定的线/曲线，在2D标绘y-x图，会给我们这两组之间的相界。

我想在2D绘制y-x图该线/曲线，从两个表面的交叉产生。

我思考如何解决这个问题的方式：

我们可以定义一个新函数，w = z_{1} - z_{2}。这两个表面之间的交叉点将是w = (z_{1} - z_{2}) = 0。然后我可以定义两个区域：

a) 一个区域 w = 0

b) 一个区域 w \neq 0

如果我绘制的这两个值w在2D y-x图：

然后我可以定义这条线/曲线是这两组之间的相边界：

a)w = 0两个集合共存的区域

b)w \neq 0两个集合不共存的区域

为什么我无法使用此解决方案取得进展：

如果我们只是删除.dat文件上的空行并进行排序x- 明智的：

sed '/^\s*$/d' 1__scatter_xyz.dat | grep -v "^#" | sort -k1 -n > 1__scatter_xyz_sort_x_wise.dat

sed '/^\s*$/d' 2__scatter_xyz.dat | grep -v "^#" | sort -k1 -n > 2__scatter_xyz_sort_x_wise.dat

如果您查看两个x_wise.dat文件，就会发现数据重叠： set 1从y-4.41 到 10.85，set 2从 8.06 到 17.64。两组的数组y不同。但是，数组x是相同的：从 10 到 2000，步长为 20.1。

因此，集合 1 和集合 2 具有相同的数组x_{j}：从 10 到 2000，步长为 20.1。但是，这两个集合没有相同的ys数组：有一个数组y_{i}^{1}forset 1和一个数组y_{i}^{2}for set 2。

换句话说，

因此，假设我找到了一个点，其中两个曲面的 值相同z。该点将由x_{j},y_{i}^{1}和y_{i}^{2}而不是两个唯一坐标定义。

更有效的想法是非常受欢迎的。

为此使用scipy'sgriddata ：

import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# Load data:
x_1, y_1, z_1 = = np.loadtxt(./1__scatter_xyz.dat, skiprows = 1).T
x_2, y_2, z_2 = = np.loadtxt(./2__scatter_xyz.dat, skiprows = 1).T

# According to the example posted in the above scipy's griddata link,
# variables "points" and "values" are defined, so we can similarly use:    
points_1 = (x_1, y_1)
points_2 = (x_2, y_2)
values_1 = (z_1)
values_2 = (z_2)

我们现在必须定义网格。正如帖子中深入解释的那样， 的数组y在两组上的采样方式不同。

我仔细研究了数据，在y空间上有两个集合之间的重叠区域：

因此，在继续这个scipy的griddata例子，我们可以设置：

T_initial = 10.0
T_end = 2000.0
number_of_Ts = 100

P_initial = 8.0622
P_end = 10.8535
number_of_Ps = 100

# And then define the mesh as:    
grid_T, grid_P = np.meshgrid(np.linspace(T_initial, T_end, number_of_Ts), np.linspace(P_initial, P_end, number_of_Ps))

此时我不知道如何继续，因为我们实际上可以只定义两组网格？

grid_Gibbs_solid_1 = griddata(points_solid_1, values_solid_1, (grid_T, grid_P), method='cubic')    
grid_Gibbs_solid_2 = griddata(points_solid_2, values_solid_2, (grid_T, grid_P), method='cubic')

要遵循哪种方法？

Answer 1

让f(x,y)和g(x,y)表示对应于您的两个表面的函数。您正在寻找的是绘制与方程相对应的轮廓f(x,y) == g(x,y)，或者等效地f(x,y) - g(x,y) == 0。

Matplotlib 提供了contour用于此目的的函数。作为一个简单的例子，考虑函数给出的两个曲面

import numpy as np    

def f(x, y):
    return np.exp(-(x**2 + y**2))

def g(x, y):
    return (3*x**2 + y**2)/16

以下代码段绘制了函数f-g、对应于 的（3D）轮廓f-g==0以及它在z平面上的（2D）投影

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig = plt.figure(figsize = (8,8))
ax = fig.gca(projection='3d')
X = np.linspace(-2, 2, 30)
Y = np.linspace(-2, 2, 30)
X, Y = np.meshgrid(X, Y)

Z = g(X,Y)

ax.plot_surface(X, Y, f(X,Y)-g(X,Y), rstride=1, cstride=1, cmap = cm.viridis, antialiased=False, alpha = 0.5)
ax.contour(X, Y, f(X,Y) - g(X,Y), zdir='z', offset=-2, levels = [0])
ax.contour(X, Y, f(X,Y) - g(X,Y), levels = [0])
ax.set_zlim(zmin = -2)

在您的情况下，您有数据样本而不是函数。您可以通过插值轻松地从您的数据中获得表面的（近似）函数（请参阅 参考资料scipy.interpolate）。