Two variables polynomial fit with Python

I have two numpy arrays (the first is 2D, the second 1D) in the form:

$X = [[x_1,y_1],[x_2,y_2],[x_3,y_3],...]$

$Z = [z_1,z_2,z_3,...]$

I would like to fit them as I expect they respect a polynomial law.

$z = A xy + B x + C y + D$ (the model is separately linear in $x$ and $y$)

So I would like a function which takes the two arrays and gives the coefficients $A,B,C$ and $D$. Is there any way to do this?