SymPy – symbolic 2D interpolation, how to make scipy.interpolate.interp2d symbolic

I am trying for symbolic 2D interpolation of scattered points. I made use of the

f = scipy.interpolate.interp2d(X, Y, Z)

but I am unable to make the interpolated function symbolic using SymPy. I need it to be symbolic for symbolic differentiation and other stuff.

I tried

x = sympy.Symbol("r")
y = sympy.Symbol("s")
f_symbolic = sympy.lambdify((x, y), f(x, y))

which gives

TypeError: Cannot cast array data from dtype(‘O’) to dtype(‘float64’) according to the rule ‘safe’

Could you guide me or help me to achieve symbolic 2D interpolation or how to make the scipy.interpolate.interp2d object symbolic?

Note: The interpolated function loooks like this. It is given by the 4 corner points.

3D plot of the interpolated function (shape function of Q4 FEM element)


As I meantioned in the comment, for this particular function, it is possible to get the symbolic function as a solution of a system of equations.

# (x, y, z) of the 4 corner points which are given
x = np.array([[0, 1],
              [0, 1]])
y = np.array([[0, 0],
              [-1, -1]])
z = np.array([[1, 0],
              [0, 0]])

# looking for a function z(x, y) = ax + bxy + cy + d
a, b, c, d = sympy.var("a b c d")   # unknown coeffs
rows, cols = x.shape
eqs = []
# creating an equation for each of the point
for i_row in range(rows):
    for i_col in range(cols):
        eqs.append(a*x[i_row, i_col] + b*x[i_row, i_col]*y[i_row, i_col] + c*y[i_row, i_col] + d - z[i_row, i_col])
sols = sympy.solve(eqs, [a, b, c, d])
# print(sols)
# {a: -1.00000000000000, b: -1.00000000000000, c: 1.00000000000000, d: 1.00000000000000}
x_ = sympy.Symbol("x")
y_ = sympy.Symbol("y")
# resulting symbolic function
f_symbolic = sympy.Lambda((x_, y_), sols[a]*x_ + sols[b]*x_*y_ + sols[c]*y_ + sols[d])

Maybe this workaround idea helps someone.