numpy.gradient 函数的反函数
Inverse of numpy.gradient function
我需要创建一个与 np.gradient 函数相反的函数。
其中 Vx、Vy 数组（速度分量向量）是输入，输出是数据点 x、y 处的反导数（到达时间）数组。
我在 (x,y) 网格上有数据，每个点都有标量值（时间）。
我已经使用 numpy 梯度函数和线性插值来确定每个点的梯度向量 Velocity (Vx,Vy)（见下文）。
我是通过以下方式实现的：
#LinearTriInterpolator applied to a delaunay triangular mesh
LTI= LinearTriInterpolator(masked_triang, time_array)
#Gradient requested at the mesh nodes:
(Vx, Vy) = LTI.gradient(triang.x, triang.y)
下图第一张是各点的速度矢量，点标代表形成导数(Vx,Vy)的时间值
下图显示导数 (Vx,Vy) 的结果标量值，绘制为带有相关节点标签的彩色等高线图。
所以我的挑战是：
我需要逆向处理！
使用梯度向量 (Vx,Vy) 或结果标量值来确定该点的原始时间值。
这可能吗？
知道 numpy.gradient 函数是使用内部点中的二阶准确中心差异和边界处的一阶或二阶准确单边（向前或向后）差异计算的，我确信是一个可以逆转这个过程的函数。
我在想，在原始点（x1,y1 处的 t=0）与 Vx,Vy 平面上的任意点 (xi,yi) 之间取一条直线导数会得到速度分量的总和。然后我可以将该值除以两点之间的距离以获得所用时间..
这种方法行得通吗？如果是这样，最好应用哪个 numpy 集成函数？
可以在此处找到我的数据示例 [http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820]
非常感谢您的帮助
编辑：
也许这张简化的图可能有助于理解我想要到达的地方..
编辑：
感谢贡献此代码的@Aguy。我尝试使用间距为 0.5 x 0.5m 的网格并计算每个网格点的梯度来获得更准确的表示，但是我无法集成它正确。我也有一些边缘影响正在影响我不知道如何纠正的结果。
import numpy as np
from scipy import interpolate
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
#Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x,y), Arrival_Time, (xgrid, ygrid), method='linear') #Interpolating the Time values
#Formatdata
X = np.ravel(xgrid)
Y= np.ravel(ygrid)
zs = np.ravel(grid_z1)
Z = zs.reshape(X.shape)
#Calculate Gradient
(dx,dy) = np.gradient(grid_z1) #Find gradient for points on meshgrid
Velocity_dx= dx/stepx #velocity ms/m
Velocity_dy= dy/stepx #velocity ms/m
Resultant = (Velocity_dx**2 + Velocity_dy**2)**0.5 #Resultant scalar value ms/m
Resultant = np.ravel(Resultant)
#Plot Original Data F(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x,y,Arrival_Time,color='r')
ax.plot_trisurf(X, Y, Z)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Time (ms)')
pyplot.show()
#Plot the Derivative of f'(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(X,Y,Resultant,color='r',s=0.2)
ax.plot_trisurf(X, Y, Resultant)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Velocity (ms/m)')
pyplot.show()
#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2],
dyintegral[i, len(yy) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
现在 np.gradient 应用于每个网格节点 (dx,dy) = np.gradient(grid_z1)
现在，在我的过程中，我会分析上面的梯度值并进行一些调整（正在创建一些不正常的边缘效果，我需要纠正）然后整合这些值以返回到一个表面与上面显示的 f(x,y) 非常相似。
我需要一些帮助来调整集成功能：
#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2],
dyintegral[i, len(yy) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
现在我需要计算原始 (x,y) 点位置的新 'Time' 值。
更新 (08-09-20) ：在@Aguy 的帮助下，我得到了一些有希望的结果。结果如下图所示（蓝色轮廓代表原始数据，红色轮廓代表积分值）。
我仍在研究一种集成方法，可以消除 min(y) 和 max(y) 区域的不准确性
from matplotlib.tri import (Triangulation, UniformTriRefiner,
CubicTriInterpolator,LinearTriInterpolator,TriInterpolator,TriAnalyzer)
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
#-------------------------------------------------------------------------
# STEP 1: Import data from Excel file, and set variables
#-------------------------------------------------------------------------
df_initial = pd.read_excel(
r'C:\Users\morga\PycharmProjects\venv\Development\Trial'
r'.xlsx')
可以在此处找到输入数据 link
df_initial = df_initial .sort_values(by='Delay', ascending=True) #Update dataframe and sort by Delay
x = df_initial ['X'].to_numpy()
y = df_initial ['Y'].to_numpy()
Arrival_Time = df_initial ['Delay'].to_numpy()
# Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x, y), Arrival_Time, (xgrid, ygrid), method='linear') # Interpolating the Time values
# Calculate Gradient (velocity ms/m)
(dy, dx) = np.gradient(grid_z1) # Find gradient for points on meshgrid
Velocity_dx = dx / stepx # x velocity component ms/m
Velocity_dy = dy / stepx # y velocity component ms/m
# Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1) * stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0) * stepy
valintegral = np.ma.zeros(dxintegral.shape) # Makes an array filled with 0's the same shape as dx integral
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum(
[dxintegral[0, len(xx) // 2], dyintegral[i, len(xx) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral[np.isnan(dx)] = np.nan
min_value = np.nanmin(valintegral)
valintegral = valintegral + (min_value * -1)
##Plot Results
fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(x, y, color='black', s=7, zorder=3)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.contour(xgrid, ygrid, valintegral, levels=50, colors='red', zorder=2)
ax.contour(xgrid, ygrid, grid_z1, levels=50, colors='blue', zorder=1)
ax.set_aspect('equal')
plt.show()
这是一种方法：
首先，为了能够进行整合，最好是在一个规则的网格上。在此处使用变量名称 x 和 y 作为 triang.x 和 triang.y 的缩写，我们可以首先创建一个网格：
import numpy as np
n = 200 # Grid density
stepx = (max(x) - min(x)) / n
stepy = (max(y) - min(y)) / n
xspace = np.arange(min(x), max(x), stepx)
yspace = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xspace, yspace)
然后我们可以使用相同的 LinearTriInterpolator 函数在网格上插入 dx 和 dy：
fdx = LinearTriInterpolator(masked_triang, dx)
fdy = LinearTriInterpolator(masked_triang, dy)
dxgrid = fdx(xgrid, ygrid)
dygrid = fdy(xgrid, ygrid)
现在是集成部分。原则上，我们选择的任何路径都应该使我们达到相同的价值。在实践中，由于存在缺失值和不同的密度，路径的选择对于获得合理准确的答案非常重要。
下面我选择在x方向从0到n/2处的网格中间对dxgrid进行积分。然后在 y 方向上从 0 到第 i 个兴趣点积分 dygrid。然后从 n/2 再次越过 dxgrid 到感兴趣的点 j。这是确保大部分集成路径位于大量可用数据内的简单方法，只需选择一条主要位于数据范围“中间”的路径即可。其他替代考虑会导致不同的路径选择。
所以我们这样做：
dxintegral = np.nancumsum(dxgrid, axis=1) * stepx
dyintegral = np.nancumsum(dygrid, axis=0) * stepy
然后（为清楚起见，使用蛮力）：
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(n):
for j in range(n):
valintegral[i, j] = np.ma.sum([dxintegral[0, n // 2], dyintegral[i, n // 2], dxintegral[i, j], - dxintegral[i, n // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
valintegral 将是一个任意常数的结果，它可以帮助将“零”放在你想要的地方。
此处显示您的数据：
ax.tricontourf(masked_triang, time_array)
这是我使用此方法重建的内容：
ax.contourf(xgrid, ygrid, valintegral)
希望这对您有所帮助。
如果要重新访问原始三角测量点的值，可以在 valintegral 常规网格数据上使用 interp2d。
编辑：
回复你的修改，你上面的改编有几个错误：
将行 (dx,dy) = np.gradient(grid_z1) 更改为 (dy,dx) = np.gradient(grid_z1)
在积分循环中将 dyintegral[i, len(yy) // 2] 项更改为 dyintegral[i, len(xx) // 2]
最好将行 valintegral = valintegral * np.isfinite(dxintegral) 替换为 valintegral[np.isnan(dx)] = np.nan
TL;DR;
您在这个问题上有多个挑战需要解决，主要是：
- 从梯度（矢量场）重建（标量场）的可能性
还有：
- 在非矩形网格的凹壳中观察；
- 数值二维线积分和数值不准确；
似乎可以通过选择 adhoc 插值和智能集成方式来解决（正如 @Aguy 所指出的）。
MCVE
第一次，让我们构建一个 MCVE 以突出上述关键点。
数据集
我们重新创建一个标量场及其梯度。
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
def f(x, y):
return x**2 + x*y + 2*y + 1
Nx, Ny = 21, 17
xl = np.linspace(-3, 3, Nx)
yl = np.linspace(-2, 2, Ny)
X, Y = np.meshgrid(xl, yl)
Z = f(X, Y)
zl = np.arange(np.floor(Z.min()), np.ceil(Z.max())+1, 2)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
V = np.hypot(dZdx, dZdy)
标量场如下所示：
axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.view_init(elev=25, azim=-45)
而且，矢量场看起来像：
axe = plt.contour(X, Y, Z, zl, cmap='jet')
axe.axes.quiver(X, Y, dZdx, dZdy, V, units='x', pivot='tip', cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
确实，梯度对于电位水平是正常的。我们还绘制了梯度大小：
axe = plt.contour(X, Y, V, 10, cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
原场重建
如果我们天真地从梯度重建标量场：
SdZx = np.cumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.cumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
Zhat[i,j] += np.sum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat += Z[0,0] - Zhat[0,0]
我们可以看到全局结果大致正确，但在梯度幅值较低的地方水平不太准确：
内插场重建
如果我们增加网格分辨率并选择特定的插值（通常在处理网格时），我们可以获得更精细的场重建：
r = np.stack([X.ravel(), Y.ravel()]).T
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel())
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel())
Nx, Ny = 200, 200
xli = np.linspace(xl.min(), xl.max(), Nx)
yli = np.linspace(yl.min(), yl.max(), Nx)
Xi, Yi = np.meshgrid(xli, yli)
ri = np.stack([Xi.ravel(), Yi.ravel()]).T
dZdxi = Sx(ri).reshape(Xi.shape)
dZdyi = Sy(ri).reshape(Xi.shape)
SdZxi = np.cumsum(dZdxi, axis=1)*np.diff(xli)[0]
SdZyi = np.cumsum(dZdyi, axis=0)*np.diff(yli)[0]
Zhati = np.zeros(SdZxi.shape)
for i in range(Zhati.shape[0]):
for j in range(Zhati.shape[1]):
Zhati[i,j] += np.sum([SdZyi[i,0], -SdZyi[0,0], SdZxi[i,j], -SdZxi[i,0]])
Zhati += Z[0,0] - Zhati[0,0]
哪个绝对表现更好：
所以基本上，使用临时插值法增加网格分辨率可能会帮助您获得更准确的结果。插值也解决了需要从三角形网格得到一个规则的矩形网格来进行积分的问题。
凹凸包
您还指出了边缘的不准确之处。这些是插值选择和积分方法相结合的结果。当标量场到达具有少量插值点的凹区域时，积分方法无法正确计算标量场。选择能够外推的 mesh-free 插值时问题消失。
为了说明这一点，让我们从 MCVE 中删除一些数据：
q = np.full(dZdx.shape, False)
q[0:6,5:11] = True
q[-6:,-6:] = True
dZdx[q] = np.nan
dZdy[q] = np.nan
那么插值可以构造如下：
q2 = ~np.isnan(dZdx.ravel())
r = np.stack([X.ravel(), Y.ravel()]).T[q2,:]
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel()[q2])
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel()[q2])
执行积分，我们看到除了经典边缘效应外，我们在凹形区域（船体凹形的摆动 dot-dash 线）的准确值确实较低，并且我们没有凸包外的数据，因为Clough Tocher 是一个 mesh-based 插值：
Vl = np.arange(0, 11, 1)
axe = plt.contour(X, Y, np.hypot(dZdx, dZdy), Vl, cmap='jet')
axe.axes.contour(Xi, Yi, np.hypot(dZdxi, dZdyi), Vl, cmap='jet', linestyles='-.')
axe.axes.set_aspect('equal')
axe.axes.grid()
所以基本上我们在拐角处看到的错误很可能是由于积分问题加上仅限于凸包的插值。
为了克服这个问题，我们可以选择不同的插值，例如 RBF（径向基函数内核），它能够在凸包外创建数据：
Sx = interpolate.Rbf(r[:,0], r[:,1], dZdx.ravel()[q2], function='thin_plate')
Sy = interpolate.Rbf(r[:,0], r[:,1], dZdy.ravel()[q2], function='thin_plate')
dZdxi = Sx(ri[:,0], ri[:,1]).reshape(Xi.shape)
dZdyi = Sy(ri[:,0], ri[:,1]).reshape(Xi.shape)
注意这个插值器的接口略有不同（注意参数是如何传递的）。
结果如下：
我们可以看到可以外推凸包外的区域（RBF 是无网格的）。所以选择adhoc插值绝对是解决你问题的关键点。但我们仍然需要意识到，外推可能表现良好，但在某种程度上毫无意义且危险。
解决您的问题
@Aguy 提供的答案非常好，因为它设置了一种巧妙的积分方式，不会受到凸包外缺失点的干扰。但是正如你所提到的，凸包内的凹区域不准确。
如果您希望消除检测到的边缘效应，您将不得不求助于也能够外推的插值器，或者找到另一种积分方法。
插值变化
使用 RBF 插值似乎可以解决您的问题。这是完整的代码：
df = pd.read_excel('./Trial-Wireup 2.xlsx')
x = df['X'].to_numpy()
y = df['Y'].to_numpy()
z = df['Delay'].to_numpy()
r = np.stack([x, y]).T
#S = interpolate.CloughTocher2DInterpolator(r, z)
#S = interpolate.LinearNDInterpolator(r, z)
S = interpolate.Rbf(x, y, z, epsilon=0.1, function='thin_plate')
N = 200
xl = np.linspace(x.min(), x.max(), N)
yl = np.linspace(y.min(), y.max(), N)
X, Y = np.meshgrid(xl, yl)
#Zp = S(np.stack([X.ravel(), Y.ravel()]).T)
Zp = S(X.ravel(), Y.ravel())
Z = Zp.reshape(X.shape)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
SdZx = np.nancumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.nancumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
#Zhat[i,j] += np.nansum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat[i,j] += np.nansum([SdZx[0,N//2], SdZy[i,N//2], SdZx[i,j], -SdZx[i,N//2]])
Zhat += Z[100,100] - Zhat[100,100]
lz = np.linspace(0, 5000, 20)
axe = plt.contour(X, Y, Z, lz, cmap='jet')
axe = plt.contour(X, Y, Zhat, lz, cmap='jet', linestyles=':')
axe.axes.plot(x, y, '.', markersize=1)
axe.axes.set_aspect('equal')
axe.axes.grid()
图形呈现如下：
边缘效应消失了，因为 RBF 插值可以在整个网格上进行外推。您可以通过比较 mesh-based 插值的结果来确认。
线性
Clough Tocher
积分变量顺序改变
我们也可以尝试找到更好的方法来整合和减轻边缘效应，例如。让我们更改积分变量顺序：
Zhat[i,j] += np.nansum([SdZy[N//2,0], SdZx[N//2,j], SdZy[i,j], -SdZy[N//2,j]])
使用经典的线性插值。结果还是比较正确的，但是左下角还是有边缘效果：
如您所见，问题出现在积分开始区域的轴中间，缺少参考点。
我需要创建一个与 np.gradient 函数相反的函数。
其中 Vx、Vy 数组（速度分量向量）是输入，输出是数据点 x、y 处的反导数（到达时间）数组。
我在 (x,y) 网格上有数据，每个点都有标量值（时间）。
我已经使用 numpy 梯度函数和线性插值来确定每个点的梯度向量 Velocity (Vx,Vy)（见下文）。
我是通过以下方式实现的：
#LinearTriInterpolator applied to a delaunay triangular mesh
LTI= LinearTriInterpolator(masked_triang, time_array)
#Gradient requested at the mesh nodes:
(Vx, Vy) = LTI.gradient(triang.x, triang.y)
下图第一张是各点的速度矢量，点标代表形成导数(Vx,Vy)的时间值
下图显示导数 (Vx,Vy) 的结果标量值，绘制为带有相关节点标签的彩色等高线图。
所以我的挑战是：
我需要逆向处理！
使用梯度向量 (Vx,Vy) 或结果标量值来确定该点的原始时间值。
这可能吗？
知道 numpy.gradient 函数是使用内部点中的二阶准确中心差异和边界处的一阶或二阶准确单边（向前或向后）差异计算的，我确信是一个可以逆转这个过程的函数。
我在想，在原始点（x1,y1 处的 t=0）与 Vx,Vy 平面上的任意点 (xi,yi) 之间取一条直线导数会得到速度分量的总和。然后我可以将该值除以两点之间的距离以获得所用时间..
这种方法行得通吗？如果是这样，最好应用哪个 numpy 集成函数？
可以在此处找到我的数据示例 [http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820]
非常感谢您的帮助
编辑：
也许这张简化的图可能有助于理解我想要到达的地方..
编辑：
感谢贡献此代码的@Aguy。我尝试使用间距为 0.5 x 0.5m 的网格并计算每个网格点的梯度来获得更准确的表示，但是我无法集成它正确。我也有一些边缘影响正在影响我不知道如何纠正的结果。
import numpy as np
from scipy import interpolate
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
#Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x,y), Arrival_Time, (xgrid, ygrid), method='linear') #Interpolating the Time values
#Formatdata
X = np.ravel(xgrid)
Y= np.ravel(ygrid)
zs = np.ravel(grid_z1)
Z = zs.reshape(X.shape)
#Calculate Gradient
(dx,dy) = np.gradient(grid_z1) #Find gradient for points on meshgrid
Velocity_dx= dx/stepx #velocity ms/m
Velocity_dy= dy/stepx #velocity ms/m
Resultant = (Velocity_dx**2 + Velocity_dy**2)**0.5 #Resultant scalar value ms/m
Resultant = np.ravel(Resultant)
#Plot Original Data F(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x,y,Arrival_Time,color='r')
ax.plot_trisurf(X, Y, Z)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Time (ms)')
pyplot.show()
#Plot the Derivative of f'(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(X,Y,Resultant,color='r',s=0.2)
ax.plot_trisurf(X, Y, Resultant)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Velocity (ms/m)')
pyplot.show()
#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2],
dyintegral[i, len(yy) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
现在 np.gradient 应用于每个网格节点 (dx,dy) = np.gradient(grid_z1)
现在，在我的过程中，我会分析上面的梯度值并进行一些调整（正在创建一些不正常的边缘效果，我需要纠正）然后整合这些值以返回到一个表面与上面显示的 f(x,y) 非常相似。
我需要一些帮助来调整集成功能：
#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2],
dyintegral[i, len(yy) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
现在我需要计算原始 (x,y) 点位置的新 'Time' 值。
更新 (08-09-20) ：在@Aguy 的帮助下，我得到了一些有希望的结果。结果如下图所示（蓝色轮廓代表原始数据，红色轮廓代表积分值）。
我仍在研究一种集成方法，可以消除 min(y) 和 max(y) 区域的不准确性
from matplotlib.tri import (Triangulation, UniformTriRefiner,
CubicTriInterpolator,LinearTriInterpolator,TriInterpolator,TriAnalyzer)
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
#-------------------------------------------------------------------------
# STEP 1: Import data from Excel file, and set variables
#-------------------------------------------------------------------------
df_initial = pd.read_excel(
r'C:\Users\morga\PycharmProjects\venv\Development\Trial'
r'.xlsx')
可以在此处找到输入数据 link
df_initial = df_initial .sort_values(by='Delay', ascending=True) #Update dataframe and sort by Delay
x = df_initial ['X'].to_numpy()
y = df_initial ['Y'].to_numpy()
Arrival_Time = df_initial ['Delay'].to_numpy()
# Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x, y), Arrival_Time, (xgrid, ygrid), method='linear') # Interpolating the Time values
# Calculate Gradient (velocity ms/m)
(dy, dx) = np.gradient(grid_z1) # Find gradient for points on meshgrid
Velocity_dx = dx / stepx # x velocity component ms/m
Velocity_dy = dy / stepx # y velocity component ms/m
# Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1) * stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0) * stepy
valintegral = np.ma.zeros(dxintegral.shape) # Makes an array filled with 0's the same shape as dx integral
for i in range(len(yy)):
for j in range(len(xx)):
valintegral[i, j] = np.ma.sum(
[dxintegral[0, len(xx) // 2], dyintegral[i, len(xx) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral[np.isnan(dx)] = np.nan
min_value = np.nanmin(valintegral)
valintegral = valintegral + (min_value * -1)
##Plot Results
fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(x, y, color='black', s=7, zorder=3)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.contour(xgrid, ygrid, valintegral, levels=50, colors='red', zorder=2)
ax.contour(xgrid, ygrid, grid_z1, levels=50, colors='blue', zorder=1)
ax.set_aspect('equal')
plt.show()
这是一种方法：
首先，为了能够进行整合，最好是在一个规则的网格上。在此处使用变量名称 x 和 y 作为 triang.x 和 triang.y 的缩写，我们可以首先创建一个网格：
import numpy as np
n = 200 # Grid density
stepx = (max(x) - min(x)) / n
stepy = (max(y) - min(y)) / n
xspace = np.arange(min(x), max(x), stepx)
yspace = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xspace, yspace)
然后我们可以使用相同的 LinearTriInterpolator 函数在网格上插入 dx 和 dy：
fdx = LinearTriInterpolator(masked_triang, dx)
fdy = LinearTriInterpolator(masked_triang, dy)
dxgrid = fdx(xgrid, ygrid)
dygrid = fdy(xgrid, ygrid)
现在是集成部分。原则上，我们选择的任何路径都应该使我们达到相同的价值。在实践中，由于存在缺失值和不同的密度，路径的选择对于获得合理准确的答案非常重要。
下面我选择在x方向从0到n/2处的网格中间对dxgrid进行积分。然后在 y 方向上从 0 到第 i 个兴趣点积分 dygrid。然后从 n/2 再次越过 dxgrid 到感兴趣的点 j。这是确保大部分集成路径位于大量可用数据内的简单方法，只需选择一条主要位于数据范围“中间”的路径即可。其他替代考虑会导致不同的路径选择。
所以我们这样做：
dxintegral = np.nancumsum(dxgrid, axis=1) * stepx
dyintegral = np.nancumsum(dygrid, axis=0) * stepy
然后（为清楚起见，使用蛮力）：
valintegral = np.ma.zeros(dxintegral.shape)
for i in range(n):
for j in range(n):
valintegral[i, j] = np.ma.sum([dxintegral[0, n // 2], dyintegral[i, n // 2], dxintegral[i, j], - dxintegral[i, n // 2]])
valintegral = valintegral * np.isfinite(dxintegral)
valintegral 将是一个任意常数的结果，它可以帮助将“零”放在你想要的地方。
此处显示您的数据：
ax.tricontourf(masked_triang, time_array)
这是我使用此方法重建的内容：
ax.contourf(xgrid, ygrid, valintegral)
希望这对您有所帮助。
如果要重新访问原始三角测量点的值，可以在 valintegral 常规网格数据上使用 interp2d。
编辑：
回复你的修改，你上面的改编有几个错误：
将行
(dx,dy) = np.gradient(grid_z1)更改为(dy,dx) = np.gradient(grid_z1)在积分循环中将
dyintegral[i, len(yy) // 2]项更改为dyintegral[i, len(xx) // 2]最好将行
valintegral = valintegral * np.isfinite(dxintegral)替换为valintegral[np.isnan(dx)] = np.nan
TL;DR;
您在这个问题上有多个挑战需要解决，主要是：
- 从梯度（矢量场）重建（标量场）的可能性
还有：
- 在非矩形网格的凹壳中观察；
- 数值二维线积分和数值不准确；
似乎可以通过选择 adhoc 插值和智能集成方式来解决（正如 @Aguy 所指出的）。
MCVE
第一次，让我们构建一个 MCVE 以突出上述关键点。
数据集
我们重新创建一个标量场及其梯度。
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
def f(x, y):
return x**2 + x*y + 2*y + 1
Nx, Ny = 21, 17
xl = np.linspace(-3, 3, Nx)
yl = np.linspace(-2, 2, Ny)
X, Y = np.meshgrid(xl, yl)
Z = f(X, Y)
zl = np.arange(np.floor(Z.min()), np.ceil(Z.max())+1, 2)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
V = np.hypot(dZdx, dZdy)
标量场如下所示：
axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.view_init(elev=25, azim=-45)
而且，矢量场看起来像：
axe = plt.contour(X, Y, Z, zl, cmap='jet')
axe.axes.quiver(X, Y, dZdx, dZdy, V, units='x', pivot='tip', cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
确实，梯度对于电位水平是正常的。我们还绘制了梯度大小：
axe = plt.contour(X, Y, V, 10, cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
原场重建
如果我们天真地从梯度重建标量场：
SdZx = np.cumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.cumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
Zhat[i,j] += np.sum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat += Z[0,0] - Zhat[0,0]
我们可以看到全局结果大致正确，但在梯度幅值较低的地方水平不太准确：
内插场重建
如果我们增加网格分辨率并选择特定的插值（通常在处理网格时），我们可以获得更精细的场重建：
r = np.stack([X.ravel(), Y.ravel()]).T
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel())
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel())
Nx, Ny = 200, 200
xli = np.linspace(xl.min(), xl.max(), Nx)
yli = np.linspace(yl.min(), yl.max(), Nx)
Xi, Yi = np.meshgrid(xli, yli)
ri = np.stack([Xi.ravel(), Yi.ravel()]).T
dZdxi = Sx(ri).reshape(Xi.shape)
dZdyi = Sy(ri).reshape(Xi.shape)
SdZxi = np.cumsum(dZdxi, axis=1)*np.diff(xli)[0]
SdZyi = np.cumsum(dZdyi, axis=0)*np.diff(yli)[0]
Zhati = np.zeros(SdZxi.shape)
for i in range(Zhati.shape[0]):
for j in range(Zhati.shape[1]):
Zhati[i,j] += np.sum([SdZyi[i,0], -SdZyi[0,0], SdZxi[i,j], -SdZxi[i,0]])
Zhati += Z[0,0] - Zhati[0,0]
哪个绝对表现更好：
所以基本上，使用临时插值法增加网格分辨率可能会帮助您获得更准确的结果。插值也解决了需要从三角形网格得到一个规则的矩形网格来进行积分的问题。
凹凸包
您还指出了边缘的不准确之处。这些是插值选择和积分方法相结合的结果。当标量场到达具有少量插值点的凹区域时，积分方法无法正确计算标量场。选择能够外推的 mesh-free 插值时问题消失。
为了说明这一点，让我们从 MCVE 中删除一些数据：
q = np.full(dZdx.shape, False)
q[0:6,5:11] = True
q[-6:,-6:] = True
dZdx[q] = np.nan
dZdy[q] = np.nan
那么插值可以构造如下：
q2 = ~np.isnan(dZdx.ravel())
r = np.stack([X.ravel(), Y.ravel()]).T[q2,:]
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel()[q2])
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel()[q2])
执行积分，我们看到除了经典边缘效应外，我们在凹形区域（船体凹形的摆动 dot-dash 线）的准确值确实较低，并且我们没有凸包外的数据，因为Clough Tocher 是一个 mesh-based 插值：
Vl = np.arange(0, 11, 1)
axe = plt.contour(X, Y, np.hypot(dZdx, dZdy), Vl, cmap='jet')
axe.axes.contour(Xi, Yi, np.hypot(dZdxi, dZdyi), Vl, cmap='jet', linestyles='-.')
axe.axes.set_aspect('equal')
axe.axes.grid()
所以基本上我们在拐角处看到的错误很可能是由于积分问题加上仅限于凸包的插值。
为了克服这个问题，我们可以选择不同的插值，例如 RBF（径向基函数内核），它能够在凸包外创建数据：
Sx = interpolate.Rbf(r[:,0], r[:,1], dZdx.ravel()[q2], function='thin_plate')
Sy = interpolate.Rbf(r[:,0], r[:,1], dZdy.ravel()[q2], function='thin_plate')
dZdxi = Sx(ri[:,0], ri[:,1]).reshape(Xi.shape)
dZdyi = Sy(ri[:,0], ri[:,1]).reshape(Xi.shape)
注意这个插值器的接口略有不同（注意参数是如何传递的）。
结果如下：
我们可以看到可以外推凸包外的区域（RBF 是无网格的）。所以选择adhoc插值绝对是解决你问题的关键点。但我们仍然需要意识到，外推可能表现良好，但在某种程度上毫无意义且危险。
解决您的问题
@Aguy 提供的答案非常好，因为它设置了一种巧妙的积分方式，不会受到凸包外缺失点的干扰。但是正如你所提到的，凸包内的凹区域不准确。
如果您希望消除检测到的边缘效应，您将不得不求助于也能够外推的插值器，或者找到另一种积分方法。
插值变化
使用 RBF 插值似乎可以解决您的问题。这是完整的代码：
df = pd.read_excel('./Trial-Wireup 2.xlsx')
x = df['X'].to_numpy()
y = df['Y'].to_numpy()
z = df['Delay'].to_numpy()
r = np.stack([x, y]).T
#S = interpolate.CloughTocher2DInterpolator(r, z)
#S = interpolate.LinearNDInterpolator(r, z)
S = interpolate.Rbf(x, y, z, epsilon=0.1, function='thin_plate')
N = 200
xl = np.linspace(x.min(), x.max(), N)
yl = np.linspace(y.min(), y.max(), N)
X, Y = np.meshgrid(xl, yl)
#Zp = S(np.stack([X.ravel(), Y.ravel()]).T)
Zp = S(X.ravel(), Y.ravel())
Z = Zp.reshape(X.shape)
dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
SdZx = np.nancumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.nancumsum(dZdy, axis=0)*np.diff(yl)[0]
Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
for j in range(Zhat.shape[1]):
#Zhat[i,j] += np.nansum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
Zhat[i,j] += np.nansum([SdZx[0,N//2], SdZy[i,N//2], SdZx[i,j], -SdZx[i,N//2]])
Zhat += Z[100,100] - Zhat[100,100]
lz = np.linspace(0, 5000, 20)
axe = plt.contour(X, Y, Z, lz, cmap='jet')
axe = plt.contour(X, Y, Zhat, lz, cmap='jet', linestyles=':')
axe.axes.plot(x, y, '.', markersize=1)
axe.axes.set_aspect('equal')
axe.axes.grid()
图形呈现如下：
边缘效应消失了，因为 RBF 插值可以在整个网格上进行外推。您可以通过比较 mesh-based 插值的结果来确认。
线性
Clough Tocher
积分变量顺序改变
我们也可以尝试找到更好的方法来整合和减轻边缘效应，例如。让我们更改积分变量顺序：
Zhat[i,j] += np.nansum([SdZy[N//2,0], SdZx[N//2,j], SdZy[i,j], -SdZy[N//2,j]])
使用经典的线性插值。结果还是比较正确的，但是左下角还是有边缘效果：
如您所见，问题出现在积分开始区域的轴中间，缺少参考点。