numpy.gradient函数的逆函数 [英] Inverse of numpy.gradient function
问题描述
我需要创建一个与np.gradient函数相反的函数。
其中Vx，Vy数组（速度分量矢量）是输入和输出，是在数据点x，y处的反导数数组（到达时间）。
我在（x，y）网格上有数据，每个点都有标量值（时间）。
我已经使用numpy梯度函数和线性插值来确定每个点上的梯度矢量 Velocity （Vx，Vy）（请参见下文）。
我是通过以下方式实现的：
#LinearTriInterpolator应用于delaunay三角网格
LTI = LinearTriInterpolator（ masked_triang，time_array）
＃在网格节点处请求的渐变：
（Vx，Vy）= LTI.gradient（triang.x，triang.y）
下面的第一张图像显示每个点的速度矢量，点标签表示形成导数（Vx，Vy）
的时间值
下一张图片显示了带有颜色的轮廓图和相关节点标签的导数（Vx，Vy）的结果标量值。 / p>
所以我的挑战是：
我需要逆转过程！
使用梯度矢量（Vx，Vy）或所得标量值来确定该点的原始时间值。
这可能吗？
知道numpy.gradient函数是使用内部点的二阶精确中心差以及一阶或二阶计算的边界上准确的单边（向前或向后）差异，我是
我当时想在原始点（x1，y1处t = 0）到任意点（xi，yi）之间取线导数）在Vx，Vy平面上将得到速度分量的总和。然后，我可以将该值除以两点之间的距离即可得出所需的时间。
这种方法行得通吗？如果是这样，哪个numpy集成函数将是最合适的？
一个示例我的数据可以在这里找到[http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820]
我们将不胜感激
编辑：
也许这张简化的图纸可能有助于了解我要去的地方。.
编辑：
感谢@Aguy同意使用此代码。.我尝试使用间距为0.5 x 0.5m的网格来获得更准确的表示并计算每个网格点处的渐变，但是我无法正确整合它。我还有一些边缘影响，这些影响影响了我不知道如何纠正的结果。
将numpy导入为np
scipy import从matplotlib内插
从py_toolkits.mplot3d导入pyplot
import Axes3D
#Createmesh网格，间距为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'）＃内插时间值
#Formatdata
X = np.ravel（xgrid）
Y = np.ravel（ygrid）
zs = np.ravel（grid_z1）
Z = zs.reshape（X.shape）
＃计算梯度
（dx，dy）= np.gradient（grid_z1）＃查找网格上点的梯度
Velocity_dx = dx / stepx＃速度ms / m
Velocity_dy = dy / stepx＃速度ms / m
Resultant =（Velocity_dx ** 2 + Velocity_dy ** 2）** 0.5 #Resultant标量值ms / m
结果= np.ravel（结果）
＃在网格网格
上绘制原始数据F（X，Y）图= 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坐标'）
ax.set_ylabel（'Y坐标'）
ax.set_zlabel（'时间（ms）'）
pyplot.show （）
＃在网格上绘制f'（X，Y）的导数
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 -坐标'）
ax.set_ylabel（'Y坐标'）
ax.set_zlabel（'速度（ms / m）'）
pyplot.show（）
＃积分以比较原始数据输入
dxintegral = np.nancumsum（Velocity_dx，axis = 1）* stepx
dyintegral = np.nancumsum（Velocity_dy，axis = 0）* stepy
valintegral = np.ma.zeros（dxintegral.shape）
对于范围内的i（len（yy））：
对于范围内的j（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应用于每个meshnode（dx，dy）= np.gradient（grid_z1）
现在在我的过程中，我将分析渐变上面的值并进行一些调整（正在创建一些异常的边缘效果，需要纠正），然后将这些值积分以返回到与上面显示的f（x，y）非常相似的表面。
我需要一些帮助使用积分函数：
#Integrate比较原始数据输入
dxintegral = np.nancumsum（Velocity_dx，axis = 1）* stepx
dyintegral = np.nancumsum（Velocity_dy，axis = 0）* stepy
valintegral = np.ma.zeros（dxintegral.shape）
对于范围i（len（ yy））：
对于范围内的j（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）点位置的新时间值
更新（08-09-20）：使用@Aguy的帮助，我得到了一些可喜的结果。结果可以在下面看到（蓝色轮廓代表原始数据，红色轮廓代表积分值）。
我仍在研究一种可以消除误差的积分方法在最小（y）和最大（y）区域
从matplotlib.tri导入（三角剖分，UniformTriRefiner，
CubicTriInterpolator，LinearTriInterpolator，TriInterpolator，TriAnalyzer）
熊猫作为sdp。 matplotlib.pyplot as plt
import numpy as np
from scipy import插值
＃-------------------- -------------------------------------------------- ---
＃步骤1：从Excel文件导入数据，并设置变量
＃------------------------- ------------------------------------------------
df_initial = pd.read_excel（
r'C：\Users\morga\PycharmProjects\venv\Development\Trial-Wireup 2'
r'.xlsx'）
输入数据可在此处找到
向量字段如下所示：
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_a spect（'equal'）
axe.axes.grid（）
独立梯度对潜在水平而言是正常的。我们还绘制了梯度幅度：
axe = plt.contour（X，Y，V，10，cmap ='jet'）
axe.axes.set_aspect（'等于'）
axe.axes.grid（）
原始场重构
如果我们从梯度中天真地重建标量场：
SdZx = np.cumsum（dZdx，轴= 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]）：$范围内j的b $ b（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在范围内（Zhati。 shape [0]）：$ j范围内的
（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]
其表现肯定更好：
因此，基本上，使用临时插值器提高网格分辨率可能会帮助您获得更准确的结果。插值器还解决了从三角形网格中获取规则矩形网格以进行积分的需求。
凹面和凸面外壳
您也已经指出边缘上的误差。这些是内插选择和积分方法相结合的结果。当标量场到达内插点很少的凹面区域时，积分方法无法正确计算标量场。选择一个无网格的可插值的插值法后，问题消失了。
为了说明这一点，让我们从MCVE中删除一些数据：
q = np.full（dZdx.shape，False）
q [0：6,5：11] =真
q [-6：，-6：] =真
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]）
执行积分，我们看到除了经典的边缘效应之外，确实在凹形区域（外壳为凹形的点状虚线）中的精度值较低，并且凸形外壳之外没有数据，因为Clough Tocher是基于网格的插值器：
Vl = np.arange（0，11，1）
ax = 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（'等于'）
axe.axes.grid（）
基本上，我们在拐角处看到的错误最有可能是由于积分问题和仅限于凸包的插值组合。
为克服这一问题，我们可以选择其他插值，例如RBF（径向基函数内核），它可以在凸包之外创建数据：
Sx =内插。Rbf（r [：，0]，r [：，1]，dZdx.ravel（）[q2]，function ='thin_plate'）
Sy =内插。 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是无网格的）。因此，选择即席插值绝对是解决您的问题的关键。但是我们仍然需要意识到，外推法可能效果不错，但毫无意义和危险。
解决问题
<$提供的答案c $ c> @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 =插值.CloughTocher2DInterpolator（r，z）
#S =插值.LinearNDInterpolator（r，z）
S =插值.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，轴= 1）* np.diff（xl）[0]
SdZy = np.nancumsum（dZdy，轴= 0）* np.diff（yl）[0]
Zhat = np.zeros（SdZx.shape）
对于范围内的i（Zhat.shape [0]）：
对于范围内的j（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）
ax = plt.contour（X，Y，Z，lz，cmap ='jet'）
ax = 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插值而消失，可以在整个网格上进行推断。您可以通过比较基于网格的插值的结果来确认它。
线性
Clough Tocher
积分变量顺序更改
我们还可以尝试找到一种更好的方法来整合和减轻边缘效应，例如。让我们更改积分变量顺序：
Zhat [i，j] + = np.nansum（[SdZy [N // 2,0]， SdZx [N // 2，j]，SdZy [i，j]，-SdZy [N // 2，j]]）
具有经典的线性插值。结果是非常正确的，但是我们仍然在左下角具有边缘效果：
您注意到问题出现在积分开始且缺少参考点的区域的轴中心。
I need to create a function which would be the inverse of the np.gradient function.
Where the Vx,Vy arrays (Velocity component vectors) are the input and the output would be an array of anti-derivatives (Arrival Time) at the datapoints x,y.
I have data on a (x,y) grid with scalar values (time) at each point.
I have used the numpy gradient function and linear interpolation to determine the gradient vector Velocity (Vx,Vy) at each point (See below).
I have achieved this by:
#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)
The first image below shows the velocity vectors at each point, and the point labels represent the time value which formed the derivatives (Vx,Vy)
The next image shows the resultant scalar value of the derivatives (Vx,Vy) plotted as a colored contour graph with associated node labels.
So my challenge is:
I need to reverse the process!
Using the gradient vectors (Vx,Vy) or the resultant scalar value to determine the original Time-Value at that point.
Is this possible?
Knowing that the numpy.gradient function is computed using second order accurate central differences in the interior points and either first or second order accurate one-sides (forward or backwards) differences at the boundaries, I am sure there is a function which would reverse this process.
I was thinking that taking a line derivative between the original point (t=0 at x1,y1) to any point (xi,yi) over the Vx,Vy plane would give me the sum of the velocity components. I could then divide this value by the distance between the two points to get the time taken..
Would this approach work? And if so, which numpy integrate function would be best applied?
An example of my data can be found here [http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820]
Your help would be greatly appreciated
EDIT:
Maybe this simplified drawing might help understand where I'm trying to get to..
EDIT:
Thanks to @Aguy who has contibuted to this code.. I Have tried to get a more accurate representation using a meshgrid of spacing 0.5 x 0.5m and calculating the gradient at each meshpoint, however I am not able to integrate it properly. I also have some edge affects which are affecting the results that I don't know how to correct.
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)
Now the np.gradient is applied at every meshnode (dx,dy) = np.gradient(grid_z1)
Now in my process I would analyse the gradient values above and make some adjustments (There is some unsual edge effects that are being create which I need to rectify) and would then integrate the values to get back to a surface which would be very similar to f(x,y) shown above.
I need some help adjusting the integration function:
#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)
And now I need to calculate the new 'Time' values at the original (x,y) point locations.
UPDATE (08-09-20) : I am getting some promising results using the help from @Aguy. The results can be seen below (with the blue contours representing the original data, and the red contours representing the integrated values).
I am still working on an integration approach which can remove the inaccuarcies at the areas of min(y) and 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-Wireup 2'
r'.xlsx')
Inputdata can be found here 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()
TL;DR;
You have multiple challenges to address in this issue, mainly:
- Potential reconstruction (scalar field) from its gradient (vector field)
But also:
- Observation in a concave hull with non rectangular grid;
- Numerical 2D line integration and numerical inaccuracy;
It seems it can be solved by choosing an adhoc interpolant and a smart way to integrate (as pointed out by @Aguy).
MCVE
In a first time, let's build a MCVE to highlight above mentioned key points.
Dataset
We recreate a scalar field and its gradient.
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)
The scalar field looks like:
axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.view_init(elev=25, azim=-45)
And, the vector field looks like:
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()
Indeed gradient is normal to potential levels. We also plot the gradient magnitude:
axe = plt.contour(X, Y, V, 10, cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()
Raw field reconstruction
If we naively reconstruct the scalar field from the gradient:
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]
We can see the global result is roughly correct, but levels are less accurate where the gradient magnitude is low:
Interpolated field reconstruction
If we increase the grid resolution and pick a specific interpolant (usual when dealing with mesh grid), we can get a finer field reconstruction:
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]
Which definitely performs way better:
So basically, increasing the grid resolution with an adhoc interpolant may help you to get more accurate result. The interpolant also solve the need to get a regular rectangular grid from a triangular mesh to perform integration.
Concave and convex hull
You also have pointed out inaccuracy on the edges. Those are the result of the combination of the interpolant choice and the integration methodology. The integration methodology fails to properly compute the scalar field when it reach concave region with few interpolated points. The problem disappear when choosing a mesh-free interpolant able to extrapolate.
To illustrate it, let's remove some data from our 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
Then the interpolant can be constructed as follow:
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])
Performing the integration we see that in addition of classical edge effect we do have less accurate value in concave regions (swingy dot-dash lines where the hull is concave) and we have no data outside the convex hull as Clough Tocher is a mesh-based interpolant:
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()
So basically the error we are seeing on the corner are most likely due to integration issue combined with interpolation limited to the convex hull.
To overcome this we can choose a different interpolant such as RBF (Radial Basis Function Kernel) which is able to create data outside the convex hull:
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)
Notice the slightly different interface of this interpolator (mind how parmaters are passed).
The result is the following:
We can see the region outside the convex hull can be extrapolated (RBF are mesh free). So choosing the adhoc interpolant is definitely a key point to solve your problem. But we still need to be aware that extrapolation may perform well but is somehow meaningless and dangerous.
Solving your problem
The answer provided by @Aguy is perfectly fine as it setups a clever way to integrate that is not disturbed by missing points outside the convex hull. But as you mentioned there is inaccuracy in concave region inside the convex hull.
If you wish to remove the edge effect you detected, you will have to resort to an interpolant able to extrapolate as well, or find another way to integrate.
Interpolant change
Using RBF interpolant seems to solve your problem. Here is the complete code:
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()
Which graphically renders as follow:
The edge effect is gone because of the RBF interpolant can extrapolate over the whole grid. You can confirm it by comparing the result of mesh-based interpolants.
Linear
Clough Tocher
Integration variable order change
We can also try to find a better way to integrate and mitigate the edge effect, eg. let's change the integration variable order:
Zhat[i,j] += np.nansum([SdZy[N//2,0], SdZx[N//2,j], SdZy[i,j], -SdZy[N//2,j]])
With a classic linear interpolant. The result is quite correct, but we still have an edge effect on the bottom left corner:
As you noticed the problem occurs at the middle of the axis in region where the integration starts and lacks a reference point.
这篇关于numpy.gradient函数的逆函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
