How do I use numpy.gradient() with coordinate arrays?

Question

I have three 1D arrays representing some spatial data (u) over a rectangular region of space: the first array contains the x coordinates, the second contains the y coordinates, and the third contains u. The x and y arrays are ordered weirdly and there are many repeated elements because lots of points sampled in the region share x and y coordinates (each point is still at a different location though). I want to calculate the derivatives of u at each point with respect to x and y.

I first tried the numpy.gradient function. The numpy.gradient documentation says something about how you can enter the coordinate arrays in the *varargs field of the function, but I don't know if that is for the purpose I think it is, and I am confident it doesn't work the way I thought it would. If anyone could help me figure out how to get the derivatives of u with respect to x and y, I would be grateful. It doesn't have to be through the numpy.gradient method. I just need something where I can input the 1D arrays and it will spit out the two spatial derivatives.

Thanks in advance!

Answer 1

I don't think gradient will work with your data, however your description isn't real clear. It sounds like your points are too unordered.

Here's a simple 2d case with a single value defined on a regular x,y grid.

x and y coordinates:

In [20]: x=np.arange(10); y = np.arange(8)

A scalar function of those points; one value for each x,y point. Here I'm using broadcasting to create at (10,8) array.

In [21]: u = np.sin(x[:,None]*2/np.pi) + np.cos(y*2/np.pi)

In [22]: u.shape
Out[22]: (10, 8)

Calculating the gradient with just u, and with x and y. Since x,y are both arange, integer spacing, the two calculations give the same thing. In effect the dx and dy, spacing, is the same:

In [27]: g = np.gradient(u)    
In [28]: g1 = np.gradient(u, x, y)

In [29]: np.allclose(g[0],g1[0])
Out[29]: True

In [30]: np.allclose(g[1],g1[1])
Out[30]: True

In [31]: g[0].shape
Out[31]: (10, 8)

Since there are 2 dimensions, there are 2 gradient arrays, each (10,8). I expect that the gradient values will be cos and sin shaped, roughly the finite difference derivatives of the trig functions used the define u.

There may be something, most likely in scipy, to handle 'scattered' and 'duplicate' points, but np.gradient is not one for that.

Or maybe I'm thinking of interpolation functions that can take scatter points, and give values on a regular grid. But with that you may then be able apply gradient.