How make a correct gradient map using Numpy.gradient

Question

I made use of numpy.gradient to calculate the gradient map of a scalar-field map. I guess I have not known numpy.gradient very well so that I might produce an incorrect gradient map. I post my code and the resulting map below:

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np

subhdu = fits.open('test_subim.fits')[0]
subhdu = subhdu.data
fig = plt.figure(1, figsize = (30,30))
ax = fig.add_axes([0.1,0.7,0.5,0.2])
xr = np.arange(0, subhdu.shape[1], 1)
yr = np.arange(0, subhdu.shape[0], 1)
xx, yy = np.meshgrid(xr,yr)
dx, dy = np.gradient(subhdu.astype('float'))
im = ax.imshow(subhdu,origin='lower',cmap='bwr')
ax.quiver(xx,yy,dx,dy,scale=5,angles="uv",headwidth = 5)
fig.colorbar(im,pad=0)
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])

I am confused with two things about my resulting map:

1. Why the gradients in the masked area do not point toward as the green arrows;
2. Why the edges of the map do not have the gradient calculation?

I would appreciate if anyone can help me figure out my confusions. If you want to play with my data 'test_subim.fits', please visit my google drive. Before playing it, you must install the package astropy probably through the following command, pip install astropy.

Thank you very much again in advance for anyone who can help me out.

Answer 1

1) Because the 'white part' is not the top of the mountain, it goes, blue -> white -> red, as you can see in the bar in the right hand side. So the blue is the valley and the red are the mountains, and the arrows point where it is uphill.

2) The edges of the map don't have a gradient calculation because the gradient is calculated with respect to it's full neighbourhood. A gradient is a measure of how much the surface is changing with respect to everything around it, i.e it points towards the steepest ascent with respect to the entire neighbourhoud. If some of the everything around it is missing, like at the edges, you can't calculate it.

More mathematically, your function at the edges is not differentiable, so you can't compute a gradient.

EDIT: Let's go more in depth:

The gradient is not just the difference between two points. It is a measure for how much the surface is at that space locally. Let's take a look at a five by five example. We will calculate the gradient for the middle point. It points in the steepest direction, the direction in which if you were walking on a hill, would get you up the highest by only taking one step. How do you know this direction, you look at all the directions - let's say 1°, 2°, .. 360° - (I'm cutting some mathematical corners here but that's not important right now), take a step, see how much height you have won, and then you go back to the starting position. The direction that got you to the highest point is the direction of the gradient. How much height you have won is given by the size of the gradient (how long the arrow is).

Now let's say you are standing at the top (which is the top left pixel in a 2D view), and you want to take a step in each direction. Down left, no problem, down right, no problem, but up right and up left? There's no pixel there??? What do I do now? That's why there is no gradient.

Lets say we change the terrain from the left image to the terrain in the right image. Then the gradient will point in a direction between the - now two - heighest pixels.