You could adapt the 1D solution in the post linked to 2D as follows:
using Interpolations
f(x,y) = x^2 + y^2 # function used just to generate data
x = y = -2:0.02:2
z = f.(x, y') # z is the numeric data at coordinates (x, y)
itp = interpolate((x,y), z, Gridded(Linear()));
grad = gradient.(Ref(itp), x, y')
If the data matrix z is noisy, we could fit smoothing 2D splines before estimating the gradients. An example using Dierckx:
using Interpolations, Dierckx
f(x,y) = x^2 + y^2
x, y = -2:0.02:2, -3:0.02:3
nx, ny = length(x), length(y)
z = f.(x, y') + rand(nx, ny) # noisy data matrix over grid defined by (x,y)
s0 = 0.01*nx*ny*hypot(extrema(z)...) # test different smoothing factors
spl = Spline2D(x, y, z; kx=3, ky=3, s=s0)
zsmooth = evalgrid(spl, x, y)
itp = interpolate((x,y), zsmooth, Gridded(Linear()));
grad = gradient.(Ref(itp), x, y')