Numba compatible equivalent of numpy.gradient

I am currently working to optimize a series of functions for a problem I'm working on. Here's a link to my original question that helped me accelerate my code.

I'm now trying to use the Numba solution suggested there for another, similar function. In this case, I am creating another moving window with nested for loops that loops through all of the rows and columns of a numpy matrix containing elevation data. The particular function I'm working on at the moment requires a call to numpy.gradient

[fx, fy] = np.gradient(DEM, cellsize, cellsize)

However, when I try to implement my function using Numba, I run into an issue with gradient.

TypingError: Use of unsupported NumPy function 'numpy.gradient' or unsupported use of the function.

I couldn't find an appropriate replacement for the functionality of the Gradient function via searching online. I thought this community could help. I greatly appreciate any insights or advice you can provide.

If you require more context or information I can provide that as well.

After Max9111's comment, I decided to build two different functions, one that took the gradient out of the Numba function. first fuction

def DCE_preprocess(DEM, cellsize, w):

    [nrows, ncols] = np.shape(DEM)

    #initiate an empty array same size as dem
    rms = DEM*np.nan
#    rms = np.float32(rms)

#    #compute the directional cosines
    [fx, fy] = np.gradient(DEM, cellsize, cellsize)


    grad = np.sqrt(fx**2 + fy**2)
    asp = np.arctan2(fy, fx)



    grad=np.pi/2-np.arctan(grad) #normal of steepest slope
    asp[asp<np.pi]=asp[asp<np.pi]+[np.pi/2]
    asp[asp<0]=asp[asp<0]+[2*np.pi]

    #spherical to cartesian conversion
    r = 1
    cy = r * np.cos(grad) * np.sin(asp)
    cx = r * np.cos(grad) * np.cos(asp)
    cz = r * np.sin(grad)

    return(cx,cy,cz)

returns input for the second function

eps = np.finfo(float).eps

def DC_eig_par(DEM,w,cx,cy,cz,eps):
    [nrows, ncols] = np.shape(DEM)
#
#    #initiate an empty array same size as dem
    rms = DEM*np.nan

    #Compute RMS cycling through the DEM
    nw=(w*2)**2
    for i in nb.prange(w+1,nrows-w):


        for j in range(w+1,(ncols-w)):
            d1=np.int64(np.linspace(i-w,i+w,11))
            d2=np.int64(np.linspace(j-w,j+w,11))

            tempx = cx[d1[0]:d1[-1],d2[0]:d2[-1]]

            tx=np.reshape(tempx,-1)
            tempy = cy[d1[0]:d1[-1],d2[0]:d2[-1]]
            ty=np.reshape(tempy,-1)
            tempz = cz[d1[0]:d1[-1],d2[0]:d2[-1]]
            tz=np.reshape(tempz,-1)

            if np.isnan(np.concatenate((tx,ty,tz))) == 0:
                T=np.array([[np.sum(tx**2), np.sum(tx*ty), np.sum(tx*tz)],
                             [np.sum(ty*tx), np.sum(ty**2), np.sum(ty*tz)], 
                             [np.sum(tz*tx), np.sum(tz*ty), np.sum(tz**2)]])


                [Te,_] = np.linalg.eig(T) # this step is a bit different from the matlab version b/c np.eig outputs two values.
                l = (Te/nw)
                l[l<eps] = 0
                rms[i,j] = 1/np.log(l[0]/l[1])


            else:
                rms[i,j] = np.nan

    return(rms)

I now get the following error:

TypingError: Invalid use of Function(<built-in function getitem>) with argument(s) of type(s): (int64, Literal[int](0))
 * parameterized
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
In definition 2:
    All templates rejected with literals.
In definition 3:
    All templates rejected without literals.
In definition 4:
    All templates rejected with literals.
In definition 5:
    All templates rejected without literals.
In definition 6:
    All templates rejected with literals.
In definition 7:
    All templates rejected without literals.
In definition 8:
    All templates rejected with literals.
In definition 9:
    All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at U:\GHP\Projects\NSF - Morris Landslides\Code\Developmemnt\Wavelets\DC_eig_par.py (58)
[2] During: typing of static-get-item at U:\GHP\Projects\NSF - Morris Landslides\Code\Developmemnt\Wavelets\DC_eig_par.py (58)

Which appears to be an issue with the slicing I'm trying to do using tempx/y/z

from Numba compatible equivalent of numpy.gradient