Inverse of numpy.gradient function

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)

enter image description here

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()

enter image description here

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()

enter image description here

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:

enter image description here

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:

enter image description here

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()

enter image description here

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:

enter image description here

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:

enter image description here enter image description here

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

enter image description here enter image description here

Clough Tocher

enter image description here enter image description here

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:

enter image description here

As you noticed the problem occurs at the middle of the axis in region where the integration starts and lacks a reference point.


Here is one approach:

First, in order to be able to do integration, it's good to be on a regular grid. Using here variable names x and y as short for your triang.x and triang.y we can first create a grid:

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)

Then we can interpolate dx and dy on the grid using the same LinearTriInterpolator function:

fdx = LinearTriInterpolator(masked_triang, dx)
fdy = LinearTriInterpolator(masked_triang, dy)

dxgrid = fdx(xgrid, ygrid)
dygrid = fdy(xgrid, ygrid)

Now comes the integration part. In principle, any path we choose should get us to the same value. In practice, since there are missing values and different densities, the choice of path is very important to get a reasonably accurate answer.

Below I choose to integrate over dxgrid in the x direction from 0 to the middle of the grid at n/2. Then integrate over dygrid in the y direction from 0 to the i point of interest. Then over dxgrid again from n/2 to the point j of interest. This is a simple way to make sure most of the path of integration is inside the bulk of available data by simply picking a path that goes mostly in the "middle" of the data range. Other alternative consideration would lead to different path selections.

So we do:

dxintegral = np.nancumsum(dxgrid, axis=1) * stepx
dyintegral = np.nancumsum(dygrid, axis=0) * stepy

and then (by somewhat brute force for clarity):

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 would be the result up to an arbitrary constant which can help put the "zero" where you want.

With your data shown here:

ax.tricontourf(masked_triang, time_array) ax.tricontourf(masked_triang, time_array)

This is what I'm getting reconstructed when using this method:

ax.contourf(xgrid, ygrid, valintegral) ax.contourf(xgrid, ygrid, valintegral)

Hopefully this is somewhat helpful.

If you want to revisit the values at the original triangulation points, you can use interp2d on the valintegral regular grid data.

EDIT:

In reply to your edit, your adaptation above has a few errors:

  1. Change the line (dx,dy) = np.gradient(grid_z1) to (dy,dx) = np.gradient(grid_z1)

  2. In the integration loop change the dyintegral[i, len(yy) // 2] term to dyintegral[i, len(xx) // 2]

  3. Better to replace the line valintegral = valintegral * np.isfinite(dxintegral) with valintegral[np.isnan(dx)] = np.nan

Tags:

Python

Numpy

Integration

Scipy

Data Science

Related

Jackson XML deserialization skips field when using multiple useWrapping = false How are tail-position contexts GHC join points paper formed? Is the writer's reason correct for using thenCompose and not thenComposeAsync How to install python3.9 with conda? How to use FileUpload widget in jupyter lab? TypeError: cb.apply is not a function Warning: Can't perform a React state update on an unmounted component. In a functional component How to get all of the pasted string, in input which has a maxLength attribute? Sum of consecutive pairs in a list including a sum of the last element with the first How can I get the current activity while using Firebase Cloud Messaging in Unity? How to start/resume and stop/pause a thread inside the action listener in java In Which Layer, Dao or Service, Should I Parse a Rest Client Response?

Recent Posts