Can scipy's RegularGridInterpolator return both values and gradients with a single call?

No.

Here is what scipy.interpolate.RegularGridInterpolator does under the hood:

class CartesianGrid(object):
    """
    Linear Multivariate Cartesian Grid interpolation in arbitrary dimensions
    This is a regular grid with equal spacing.
    """
    def __init__(self, limits, values):
        self.values = values
        self.limits = limits

    def __call__(self, *coords):
        # transform coords into pixel values
        coords = numpy.asarray(coords)
        coords = [(c - lo) * (n - 1) / (hi - lo) for (lo, hi), c, n in zip(self.limits, coords, self.values.shape)]

        return scipy.ndimage.map_coordinates(self.values, coords, 
            cval=numpy.nan, order=1)

https://github.com/JohannesBuchner/regulargrid/blob/master/regulargrid/cartesiangrid.py

It uses scipy.ndimage.map_coordinates to do the linear interpolation. coords contains the location in pixel coordinates. You should be able to use these weights, and the lower and upper values at each dimension to figure out how steep the interpolation rises.

However, the gradient also depends on the values of the corner points.

You can find the math here: https://en.wikipedia.org/wiki/Trilinear_interpolation