The key phrase from the docs is:

The gradient is computed using second order accurate central differences in the interior points ...

This means that each group of three contiguous points is adjusted to a parabola (2nd order polynomial) and its slope at the location of the central point is used as the gradient.

For evenly spaced data (with a spacing of 1) the formula is very simple:

g(i) = 0.5 f(i+1) - 0.5 f(i-1)

Then comes the problematic part:

... and either first or second order accurate one-sides (forward or backwards) differences at the boundaries.

There is neither f(i+1) at the right boundary nor f(i-1) at the left boundary.

So you can use a simple 1st order approximation

g(0) = f(1) - f(0)
g(n) = f(n) - f(n-1)

or a more complex 2nd order approximation

g(0) = -1.5 f(0) + 2 f(1) - 0.5 f(2)
g(n) = 0.5 f(n-2) - 2 f(n-1) + 1.5 f(n)

The effect can be seen in this example, copied from the docs:

>>> x = np.array([0, 1, 2, 3, 4])
>>> f = x**2
>>> np.gradient(f, edge_order=1)
array([1., 2., 4., 6., 7.])
>>> np.gradient(f, edge_order=2)
array([0., 2., 4., 6., 8.])

The derivation of the coefficient values can be found here.

The derivative of a function $f$ can be approximated at a point $x$ by finite differencing. For example, using the central difference where $\Delta$ is some small value:

$$\frac{d}{dx} f(x) \approx \frac{f(x+\Delta) - f(x-\Delta)}{2\Delta}$$

Now, suppose we don't know the expression for $f$. But, we have an evenly spaced grid of points $\{x_1, \dots, x_n\}$ and the corresponding function values $\{y_1, \dots, y_n\}$, where $y_i = f(x_i)$. Then, we can still use finite differencing to approximate the derivative at each point as above, provided the grid spacing is small enough. For the middle points ($x_2, \dots, x_{n-1}$):

$$\frac{d}{dx} f(x_i) \approx \frac{y_{i+1} - y_{i-1}}{x_{i+1}-x_{i-1}}$$

As above, this is just the change in function value divided by the change in the input. Central differences can't be used at the edge points ($x_1$ and $x_n$) but the concept is the same.

So, sampled values are sufficient for numerically approximating the derivative of a function, and the function itself isn't necessary. The case is similar for approximating the gradient (i.e. partial derivatives) of functions of multiple variables--simply perform finite differencing along each direction.