A new method for estimation of intragranular strain fields in polycrystalline materials based on scanning three-dimensional X-ray diffraction (scanning 3DXRD) data is presented and evaluated. Given an a priori known anisotropic compliance, the regression method enforces the balance of linear and angular momentum in the linear elastic strain field reconstruction. By using a Gaussian process (GP), the presented method can yield a spatial estimate of the uncertainty of the reconstructed strain field. Furthermore, constraints on spatial smoothness can be optimized with respect to measurements through hyperparameter estimation. These three features address weaknesses discussed for previously existing scanning 3DXRD reconstruction methods and, thus, offer a more robust strain field estimation. The method is twofold validated: firstly by reconstruction from synthetic diffraction data, and secondly by reconstruction of a previously studied tin (Sn) grain embedded in a polycrystalline specimen. Comparison against reconstructions achieved by a recently proposed algebraic inversion technique is also presented. It is found that the GP regression consistently produces reconstructions with lower root-mean-square errors, mean absolute errors and maximum absolute errors across all six components of strain.