<p>Geospatial data are often spatially varying but measured sparsely in a two-dimensional (2D) plane. Therefore, spatial interpolation methods, such as Kriging, are frequently used to estimate values at locations without measurement and quantify the associated uncertainty. However, Kriging generally requires extensive measurements to effectively divide non-stationary geospatial data into a deterministic trend and stationary residuals (i.e., detrending) and to estimate semi-variogram parameters from the detrended residuals (i.e., semi-variogram fitting). When measurements are limited, a scenario often encountered in practice, detrending might be challenging, and the estimated semi-variogram parameters inevitably contain statistical uncertainty. The statistical uncertainty may significantly affect the subsequent Kriging interpolation, but it is often ignored in practical applications of Kriging. This study develops a 2D Kriging method (SR-Kriging) that is featured by a sparse representation of covariance function from a Bayesian perspective and explicitly models both spatial variability and statistical uncertainty for interpolation of 2D geospatial data directly from sparse measurements. The proposed method requires neither detrending nor semi-variogram fitting. Both simulated and real data are used to illustrate and validate the proposed method. Results demonstrate that the proposed method directly interpolates 2D geospatial data from limited measurements, with quantified interpolation uncertainty, and explicitly accounts for statistical uncertainty. Ignorance of statistical uncertainty may lead to an underestimation of Kriging interpolation uncertainty.</p>