The knowledge of regional groundwater level distribution is an important foundation for groundwater resource evaluation and environmental protection. Due to the limited groundwater level data observed at regional scale, kriging interpolation and deep learning methods are gradually used for regional groundwater level prediction, but their applicability and robustness lack comparative analysis. In this paper, spatial interpolation of groundwater levels in Shenshan special cooperation zone was carried out using ordinary kriging, cokriging and deep learning methods to explore the potential of the three methods in the practical application of regional groundwater level predictions. In order to investigate the effect of the training set sample size on the prediction effect of the three methods, 239 monitoring wells were divided into two groups of 76 and 163 wells for the training of the three models, respectively. The results showed that the cokriging, which considered surface elevations, was significantly better than the ordinary kriging and deep learning in regional groundwater levels predictions, when the training sample size was 76. When the training sample size was increased to 163, the prediction accuracies of ordinary kriging, cokriging, and deep learning were significantly improved. The RMSEs of the three methods on the validation dataset differed very little, but the spatial distribution characteristics of the predicted regional groundwater levels still differed significantly among the methods. In summary, when the observed data are sparse, the prediction accuracy of cokriging incorporating surface elevation is significantly higher than that of ordinary kriging and deep learning. However, the prediction accuracies of the three methods are close to each other, when the observed data are dense.