Active source seismology provides a critical constraint on the global crustal structure. However, the heterogeneous data coverage means that interpolation is necessary to fill the gap between seismic profiles. This has the potential to cause large uncertainties especially if the data are interpolated over a large distance. In previous models, geological intuition was often employed to ensure reasonable results. To investigate crustal model uncertainty, we apply geostatistical analysis to a database of active seismic investigations. Unlike previous models, our workflow in the construction of the crustal model is completely transparent. Apart from the points from the database, we only use an a priori separation in oceanic and continental domains. We calculate global maps of Moho depth and average P wave velocity in the crystalline crust. Additionally, we obtain the interpolation error and error covariance. Overall, our results agree with previous global crustal models such as Crust1.0. Our uncertainty estimates show that the Moho depth uncertainty in the most well studied areas such as North America and Europe is less than 4 km but can reach 10 km or more in frontier regions such as most of Africa. P wave velocity shows the same pattern, but is less accurate overall, due to more small-scale variation. We demonstrate the benefit of having a numerical estimate of uncertainty by propagating the uncertainty to the residual topography. We see two main uses for our crustal model in the geophysical research community: (1) as a starting model for inversions focusing on the crust and upper mantle and (2) as a starting point for including other pointwise information about crustal structure, for example, from passive seismology.