The problem is, that numpy can't give you the derivatives directly and you have two options:

With NUMPY

What you essentially have to do, is to define a grid in three dimension and to evaluate the function on this grid. Afterwards you feed this table of function values to numpy.gradient to get an array with the numerical derivative for every dimension (variable).

Example from here:

from numpy import *

x,y,z = mgrid[-100:101:25., -100:101:25., -100:101:25.]

V = 2*x**2 + 3*y**2 - 4*z # just a random function for the potential

Ex,Ey,Ez = gradient(V)

Without NUMPY

You could also calculate the derivative yourself by using the centered difference quotient.

This is essentially, what numpy.gradient is doing for every point of your predefined grid.

You need to give gradient a matrix that describes your angular frequency values for your (x,y) points. e.g.

def f(x,y):
    return np.sin((x + y))
x = y = np.arange(-5, 5, 0.05)
X, Y = np.meshgrid(x, y)
zs = np.array([f(x,y) for x,y in zip(np.ravel(X), np.ravel(Y))])
Z = zs.reshape(X.shape)

gx,gy = np.gradient(Z,0.05,0.05)

You can see that plotting Z as a surface gives:

Here is how to interpret your gradient:

gx is a matrix that gives the change dz/dx at all points. e.g. gx[0][0] is dz/dx at (x0,y0). Visualizing gx helps in understanding:

Since my data was generated from f(x,y) = sin(x+y) gy looks the same.

Here is a more obvious example using f(x,y) = sin(x)...

f(x,y)

and the gradients

update Let's take a look at the xy pairs.

This is the code I used:

def f(x,y):
    return np.sin(x)
x = y = np.arange(-3,3,.05)
X, Y = np.meshgrid(x, y)
zs = np.array([f(x,y) for x,y in zip(np.ravel(X), np.ravel(Y))])
xy_pairs = np.array([str(x)+','+str(y) for x,y in zip(np.ravel(X), np.ravel(Y))])
Z = zs.reshape(X.shape)
xy_pairs = xy_pairs.reshape(X.shape)

gy,gx = np.gradient(Z,.05,.05)

Now we can look and see exactly what is happening. Say we wanted to know what point was associated with the value atZ[20][30]? Then...

>>> Z[20][30]
-0.99749498660405478

And the point is

>>> xy_pairs[20][30]
'-1.5,-2.0'

Is that right? Let's check.

>>> np.sin(-1.5)
-0.99749498660405445

Yes.

And what are our gradient components at that point?

>>> gy[20][30]
0.0
>>> gx[20][30]
0.070707731517679617

Do those check out?

dz/dy always 0 check. dz/dx = cos(x) and...

>>> np.cos(-1.5)
0.070737201667702906

Looks good.

You'll notice they aren't exactly correct, that is because my Z data isn't continuous, there is a step size of 0.05 and gradient can only approximate the rate of change.

In Python you can use the numpy.gradient function to do this. This said function uses central differences for the computation, like so: \begin{eqnarray} \nabla_x I (i, j) = \frac{I(i + 1, j) - I(i - 1, j)}{2}, \hspace{.5em}\nabla_y I (i, j) = \frac{I(i, j+1) - I(i, j-1)}{2}. \end{eqnarray}

Here is a code snippet for your specific image:

import numpy as np
import matplotlib.pyplot as plt

# load image
img = np.array([[21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 99.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0]])
print "image =", img

# compute gradient of image
gx, gy = np.gradient(img)
print "gx =", gx
print "gy =", gy

# plotting
plt.close("all")
plt.figure()
plt.suptitle("Image, and it gradient along each axis")
ax = plt.subplot("131")
ax.axis("off")
ax.imshow(img)
ax.set_title("image")

ax = plt.subplot("132")
ax.axis("off")
ax.imshow(gx)
ax.set_title("gx")

ax = plt.subplot("133")
ax.axis("off")
ax.imshow(gy)
ax.set_title("gy")
plt.show()

To answer your specific question, the gradient (via central differences!) of the image at pixel with value $99$ is $0$ along the $x$ axis and $-2$ along the $y$ axis.

Use

    grad = zeros_like(x, dtype=float)

All the x values you use are integer type. The division operation will convert the result to float, but since grad is set to be like x, it will have an integer type, and the floating point results will be truncated to integer before assigned to grad[i]. Hence you get a truncated floating point value as integer, which produces the result of what you see.

You may expect the result of the gradient calculation to be 2 when x[i] == 1. But the floating point division is imprecise (in decimal), thus it may be something like 1.999998, which truncates to 1 before it is assigned. That's why your result is 1, not 2.

The worst part is that if you have set h to something like h=0.015625, which is 2**-6 and can be represented exactly in binary, the result would have been (2, 2), since the division results in precisely 2.0 in binary, and thus the truncation to integer values would result in 2. This can lead to (hidden) bugs which will only appear weeks or months later, if at all (but may lead to incorrect results in published results, with all kind of possible consequences. Just as an example of a worst-case scenario.)