Soil contamination by heavy metals has emerged as a significant global problem. Accurately mapping the spatial distribution of soil heavy metal pollutant concentrations is indispensable for effective agriculture and environmental management. Nevertheless, challenges arise in obtaining comprehensive data at all desired locations, due to limitations in measuring equipment capacity and the associated capital costs. To obtain soil heavy metal maps efficiently and accurately, this paper proposes a nonparametric spatial prediction method, based on unbiased conditional kernel density estimation (UCKDE). The proposed method incorporates the advantages of both geostatistics and machine learning, including stability, adaptability, and the ability to account for various types of auxiliary information. Additionally, it can directly predict the probabilistic density function (PDF) of soil heavy metal content at the target location based on sampling data without complex parameter settings, providing both a deterministic single value and a probabilistic prediction interval. The proposed method and ordinary kriging (OK) were implemented for the spatial prediction of the six heavy metals (As, Cd, Cu, Hg, Mn, and Sb) with the greatest coefficients of variation (CV = 0.53, 1.14, 0.66, 1.05, 0.81 and 0.74, respectively) in Qingxi Town, Chongqing, China. The results showed that the predictive capability of the proposed method (with RMSE values of 5.82, 0.61, 14.76, 0.15, 383.84, and 0.85, respectively) is superior to that of OK (with RMSE values of 5.29, 0.87, 16.37, 0.22, 493.22, and 1.58, respectively) in most cases, particularly when the CV value is high. Besides, the prediction accuracy of the proposed method can be further enhanced by incorporating parent material, resulting in RMSE values of 3.02, 0.51, 8.98, 0.08, 194.16, and 0.56, respectively. The results affirm the reliability of the proposed method and suggest its effectiveness as a tool for soil heavy metal pollution prediction in practical applications.