Order of input data effects result of 3d polynomial fit

Question

I am working on a 3D polynomial fit. The data I have are the x and y coordinate of 5 stations, and the velocities at these 5 stations. What I want to do is fit a grid through these points. In turn I will use this grid to determine the velocity gradients using the predicted velocities at each grid point.

My code is:

xx_vel = np.array([[4.78,52.32] ,[5.18,52.10], [4.45,51.97],  [4.92,51.97], [5.15,51.85]]) #location of stations in degrees longitude and latitude

X = xx_vel   #coordinates
Z = np.array([-0.00, -0.766, -0.00, -1.732, -1.00])         #velocities at 5 stations

deg_of_poly = 3
poly = PolynomialFeatures(degree=deg_of_poly)
X_ = poly.fit_transform(X)

clf = linear_model.LinearRegression()
clf.fit(X_, Z)

x_pred = np.linspace(4, 6, 27)        #defining grid points
y_pred = np.linspace(51.5, 52.7, 27)   #defining grid points
predict_x, predict_y = np.meshgrid(x_pred, y_pred)
predict_xy = np.concatenate((predict_x.reshape(-1, 1), predict_y.reshape(-1, 1)), axis=1)
predict_x_ = poly.fit_transform(predict_xy)
predict_z = clf.predict(predict_x_)
predict_z_poly = predict_z.reshape(predict_x.shape)

Using this code I obtain the following fit:

This all seemed fine, until I changed the order of the input data. So if I for example switch the first and last stations so that my input arrays are:

xx_vel = np.array([[5.15,51.85],[5.18,52.10], [4.45,51.97],  [4.92,51.97], [4.78,52.32]])
Z = np.array([-1.00, -0.766, -0.00, -1.732, -0.00]) 

I obtain a different fit. Is there something I am doing wrong? Or is there a way I can make sure I obtain the same results no matter in what order the data is given? I would think that this should not have an effect on the result.

Thanks in advance!