The form error estimation under various machining conditions is an essential step in the assessment of product surface quality generated in machining processes. Coordinate measuring machines (CMMs) are widely used to measure complicated surface form error. However, considering measurement cost, only a few measurement points are collected offline by a CMM for a part surface. Therefore, spatial statistics is adopted to interpolate more points for more accurate form error estimation. It is of great significance to decrease the deviation between the interpolated height value and the real one. Compared to univariate spatial statistics, only concerning spatial correlation of height value, this paper presents a method based on multivariate spatial statistics, co-Kriging (CK), to estimate surface form error not only concerning spatial correlation but also concerning the influence of machining conditions. This method can reconstruct a more accurate part surface and make the estimation deviation smaller. It characterizes the spatial correlation of machining errors by variogram and cross-variogram, and it is implemented on one of the common features: flatness error. Simulated datasets as well as actual CMM data are applied to demonstrate the improvement achieved by the proposed multivariate spatial statistics method over the univariate method and other interpolation methods.