Why my np.gradient calculation in R^2 doesn't fit with the analytical gradient calculation?

I'm trying to compute a gradient on a map using np.gradient, but I'm encountering issues. To simplify my problem I am trying on an analytical function

z = f(x,y) = -(x - 2)**2 - (y - 2)**2

np.gradient is not providing the expected results; the vectors should point towards the center. What am I doing wrong?

Here is the code that I am running:

import numpy as np
import matplotlib.pyplot as plt

# Define the grids for x and y
x = np.linspace(0, 4, 100)  # 100 points between 0 and 4
y = np.linspace(0, 4, 100)  # 100 points between 0 and 4
X, Y = np.meshgrid(x, y)  # Create a 2D grid

# Define the function f(x, y)
Z = -(X - 2)**2 - (Y - 2)**2

# Compute gradients numerically
dz_dx, dz_dy = np.gradient(Z, x, y)


# Downsampling to reduce the density of arrows
step = 10

plt.figure(figsize=(10, 8))
contour = plt.contourf(X, Y, Z, cmap='viridis', levels=50, alpha=0.8)
plt.colorbar(contour, label='f(x, y)')
plt.quiver(X[::step, ::step], Y[::step, ::step], dz_dx[::step, ::step], dz_dy[::step, ::step], 
           color='r', headlength=3, headwidth=4)
plt.title('Function $f(x, y) = -(x - 2)^2 - (y - 2)^2$ and its gradients (numerical)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)

plt.show()

the problem is not in the gradient function, it is in the different indexing order of np.meshgrid and np.gradient.

by default np.gradient assumes the indexing is the same order as the arguments, ie:

Z[x,y] -> np.gradient(Z, x, y)

whereas np.meshgrid default indexing results in the opposite indexing,

Z[y,x] -> np.meshgrid(x,y)  # default indexing = 'xy'

you didn't notice this bug because both X and Y are identical, if you made X and Y have different number of points you would get an error.

I like to use Z[y,x], so just swap the order of arguments and return of np.gradient and you will get the correct result.

dz_dy, dz_dx = np.gradient(Z, y, x)  # Z[y,x]

I'm not sure why np.gradient is not working as your expect it, but you could instead manually calculate the gradients, e.g.,

import numpy as np
import matplotlib.pyplot as plt

# Define the grids for x and y
x = np.linspace(0, 4, 100)  # 100 points between 0 and 4
y = np.linspace(0, 4, 100)  # 100 points between 0 and 4
X, Y = np.meshgrid(x, y)  # Create a 2D grid

# Define the function f(x, y)
def func(X, Y):
    return -(X - 2)**2 - (Y - 2)**2

def grad(func, X, Y, delta=0.01):
    """
    Calculate gradients of a function.
    """

    dp = delta / 2
    df_dx = (func(X + dp, Y) - func(X - dp, Y)) / delta
    df_dy = (func(X, Y + dp) - func(X, Y - dp)) / delta

    return df_dx, df_dy

Z = func(X, Y)
dz_dx, dz_dy = grad(func, X, Y)

# Downsampling to reduce the density of arrows
step = 10
start = step // 2

plt.figure(figsize=(10, 8))
contour = plt.contourf(X, Y, Z, cmap='viridis', levels=50, alpha=0.8)
plt.colorbar(contour, label='f(x, y)')
plt.quiver(X[start::step, start::step], Y[start::step, start::step], dz_dx[start::step, start::step], dz_dy[start::step, start::step], 
           color='r', headlength=3, headwidth=4)
plt.title('Function $f(x, y) = -(x - 2)^2 - (y - 2)^2$ and its gradients (numerical)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)

plt.show()

As Ahmed AEK explained in their answer, the issue is with the orderings. Rather than changing how you call np.gradient, you should change how you call np.meshgrid.

If you look at the np.meshgrid documentation, you'll see that by default indexing="xy". The other option is indexing="ij". This controls how the points in the meshgrid are defined/ordered.

For "xy" ordering, you can imagine the array as a cartesian coordinate plane. The bottom left is (0,0), the x value increases as you move right and the y value increases as you move up.

For "ij" ordering, the indexing is like the coordinates of a cartesian plane. Meaning, (0,0) when used as an array index is the top left, increasing the x value would mean increasing the first value in the array index, so moving down (going row-by-row) increases the x value. The second index is the y "coordinate", so moving across (column-by-column) increases the y value.

So, all you need to do is change your meshgrid call to:

X, Y = np.meshgrid(x, y, indexing="ij")