Large sample sizes are crucial for accurately capturing spatial changes in soil properties by spatial interpolation methods. However, soil bulk density (BD) data in historical datasets is often incomplete, and it's uncertain if filled values enhance spatial interpolation accuracy. Using 2,883 cropland soil BD samples from the Sichuan Basin in China, we developed the best prediction models from traditional pedotransfer function (PTF), multiple linear regression (MLR), random forest (RF), and radial basis function neural network (RBFNN) to fill missing BD values for 1,336 samples. We then applied ordinary kriging (OK) and inverse distance weighting (IDW) to map soil BD, incorporating the filled BD as modeling points. The RBFNN model, tailored for each sub-watershed, yielded the highest accuracy in filling missing BD, with an increase in coefficient of determination (R