Is there a generalization of the Lagrange polynomial to 3D? What is a way to construct a smooth polynomial surface ($\mathbb{R}^2 \rightarrow \mathbb{R}$) with Lagrange-polynomial properties in every partial derivative? I want to try this for image interpolation.

Introduce an auxiliary variable $t$ so that all three variables -- $x$, $y$ and $f$ -- are functions of $t$. You would have zeroes in four dimensional space like this:
t x y f
- - - -
1 3 7 4
2 4 7 3
3 4 4 9

The idea is to leverage Lagrange interpolation to have explicit formula for $f(t)$. However, with tabulation in this naive form, then it would be quite challenging to eliminate this newly introduced variable $t$. (It can be done with Grobner basis).
However, according to [1] almost every linear combination of $x$ and $y$ would also work. In our example let $t = 100 x+y$:
t    x y f
---  - - -
307  3 7 4
407  4 7 3
404  4 4 9

Now you have $f(100 x+y)$, which satisfies your partial derivative criteria.
[1] H.J. Stetter. Numerical_Polynomial_Algebra.
