If you want to fit a multivariate polynomial, it’s just linear equations, which are trivial to set up and solve in Julia. For example, if you want to fit to z = a_1 + a_2 x + a_3 y + a_4 x y + a_5 x^2 + a_6 y^2 you could do a least-square fit with:
a = [x.^0 x y x.*y x.^2 y.^2] \ z
(The matrix [ ... ] here is a generalization of the Vandermonde matrix that appears in univariate polynomial fitting: each column corresponds to one basis function.) If the number of basis functions (columns) equals length(z), then this becomes interpolation, i.e. it will exactly go through the given points.
You can also do an equivalent fit more easily and generally with FastChebInterp.jl:
using FastChebInterp, StaticArrays
c = chebregression(SVector.(x, y), z, (2,2))
which represents the polynomial (here of degree 2 in both x and y) in the basis of Chebyshev polynomials, and you can then just evaluate it at any given point with c([x,y]) or c(SVector(x,y)).
However, you have to be careful with polynomial fitting if you go to high degree compared to the number of data points (unless the function is smooth and the points are carefully selected, e.g. using Chebyshev points via chebpoints and chebinterp), since high-degree polynomials can go crazy in between the fitted points if you are not careful.
Other options include radial basis functions (e.g. Surrogates.jl or ScatteredInterpolation.jl, both of which implement several algorithms), splines (e.g. Dierckx.jl can do splines from irregular 2d points), dividing your points into a triangular mesh (e.g. Delaunay triangulation) and doing bi-linear interpolation within each triangle, and more.