# Gradients of non-uniformly sampled data in 3D space

I have measurements of magnetic field on a 3d grid. My measurements are distributed on four x-y planes similar to what is shown in the image below. The measurements roughly follow a Cartesian grid but have varying spacing and gaps. The measurements are more sparse in the z-direction but I expect the magnetic field to vary by an order of magnitude more in the x and y direction. Each vector component of the magnetic field ($B_x$, $B_y$, $B_z$) is measured at different points.

In the image below I show the measurement points for $B_x$. The points for $B_y$, $B_z$ are distributed similarly but not at the exact same points. The whole domain is not convex but I am only looking at the region inside the black rectangle.

I would like to calculate the curl of the magnetic field to determine the current density. E.g. $j_z = 1/\mu_{0} (\partial B_y / \partial x - \partial B_x / \partial y)$.

Since the data is non-uniform and $B_x$ and $B_y$ are not measured at same locations I can not use simple finite difference relations. I have tried three approaches

1. interpolate the $B_x$ and $B_y$ with scipy.interpolate.gridata(..., method='linear') to a common grid and then take 1st order differences between the point and divide by the spacing.
2. interpolate the $B_x$ and $B_y$ with scipy.interpolate.gridata(..., method='linear') to a common grid and then take 2nd order differences between the points with gradient and divide by the spacing.
3. Triangulate the $B_x$ and $B_y$ positions with scipy.spatial.Delaunay and pass the triangulation and data to matplotlib.tri.LinearTriInterpolator which has a gradient method.

I plotted the magnetic field and $j_z$ for each of the three methods. Only the third (left to right) plot seems to be a curl of the field and gives me the magnitude that I expect from the experiment.

scipy.interpolate.griddata uses scipy.interpolate.LinearNDInterpolator which triangulates the positions and uses linear barycentric interpolation. The same is true for matplotlib.tri. Are these results expected? If all methods are using linear interpolation why is the third method more accurate than taking differences of the interpolated values?

• This is a more specific follow-up question to this scicomp question. One of the comments recommended using kriging. This is a great idea but for now I would prefer to use linear / barycentric interpolation instead of Gaussian process. Jul 27 '16 at 6:31
• I think that it would be better to use higher order interpolation if you are going to use higher order derivatives. If you think in a piece-wise interpolation the function is continuous but not the derivatives. Jul 27 '16 at 23:26