Generate a two-variable polynomial from its "roots I know, from mathematics basis, that a polynomial with one variable can be factor in function of its roots, so I can generate any one-variable polynomial from its zeros. 
But I want know if is possible to generate a two-variable polynomial giving its coordinates where value is zero? 
This problem appears on calculus of shape functions in Finite Elements Method. For example, in an one-dimensional element with four nodes $(p_1,p_2,p_3,p_4)$, I have to generate four shape functions that have the following requisites:   
The function $N_i(x)$ has value $1$ at the node $p_i$ and $0$ on others nodes.
$p_1 \to x = -1$
$p_2 \to x = -1/2$
$p_3 \to x = 1/2$
$p_4 \to x = 1$ 
$N_1(p_1) = 1$
$N_1(p_2) = N_1(p_3) = N_1(p_4) = 0$ 
From this I can find the general form of the function $N_1$: 
$N_1(x) = (x + 1/2)(x - 1/2)(x - 1)(C_1)$ 
$N_1(p_1) = 1 \to C_1 = -2/3$ 
$N_1(x) = (-2/3)(x + 1/2)(x - 1/2)(x - 1)$ 
This is the $N_1$ shape function. 
Now, I want to think this problem to two-variables and generate those functions with the same "algorithm".  
A simple square element with coordinates: $p_1=(-1,-1),p_2=(1,-1),p_3=(1,1),p_4=(-1,1)$.
And I want the $N_1(x,y)$ function that has value $1$ at $p_1$ and $0$ on other three points.
This can be possible?
 A: There a number of very well researched techniques that can be used to solve your problem in a practical sense. Most of these have come from computer graphics and computer vision where taking a set of points on a surface, possibly with noise, and trying to create a surface that either interpolates or approximates those points. For example this presentation from Cornell gives you some ideas of the techniques used - bezier patches, spline surfaces, subdivision surfaces, implicit surfaces, though of course not all those are directly related to a polynomial interpolation though most relate easily to piecewise polynomial interpolation or approximation.
Generating a continuous, smooth surface that interpolates your points may not really be necessary if you take the local interpolation viewpoint. For example if you have a function from $\mathbb{R}^2\rightarrow \mathbb{R}$ which you have sampled at a discrete set of points $(x_i,y_i)$ you can define the value of the function at a point $(x,y)$ to be a simple function of the nearby values - the simplest case being nearest neighbour interpolation but going to a slightly more complex bilinear interpolation you an achieve high accuracy assuming you have a regular grid of samples - essentially you interpolate linearly along one axis and then along the other. Note that swapping the order does not matter. This would work for your specific example:
\begin{equation}
\begin{aligned}
N_1(x,y)&=\frac{y+1}{2}(\frac{x+1}{2}N(p_3)+\frac{1-x}{2}N(p_4))+\frac{1-y}{2}(\frac{x+1}{2}N(p_2)+\frac{1-x}{2}N(p_1))\\&= \frac{(x+1)(y+1)}{4}N(p_3)+\frac{(1-x)(y+1)}{4}N(p_4)\\&+\frac{(1-x)(1-y)}{4}N(p_1)+\frac{(1+x)(1-y)}{4}N(p_2)\\&=\frac{(1-x)(1-y)}{4}
\end{aligned}
\end{equation}
