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

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 表面之间的交叉点将定义一条线/曲线，绘制在二维y-xdiagram 中，将为我们提供这两组之间的相界.

我想在 2D y-xdiagram 中绘制这条由两个表面交叉产生的线/曲线。

我对如何解决这个问题的想法：

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

a) w = 0 所在的区域

b) w \neq 0 所在的区域

如果我在二维y-xdiagram 中绘制w 的这两个值：

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

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 从 -4.41 的 y 变为 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} 有一个数组 set 1，y_{i}^{2} 有一个数组 set 2。

换句话说，

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

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

为此使用scipy's griddata ：

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空间上两组之间存在重叠区域：

所以，继续this scipy's griddata example，我们可以设置：

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

下面的 sn-p 绘制函数f-g，对应f-g==0 的（3D）轮廓以及它在z-plane 上的（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）。