Soil organic matter (SOM), total nitrogen (TN), available phosphorus (AP), and available potassium (AK) are important indicators of soil fertility when undertaking a quality evaluation. Obtaining a high-precision spatial distribution map of soil nutrients is of great significance for the differentiated management of nutrient resources and reducing non-point source pollution. However, the spatial heterogeneity of soil nutrients lead to uncertainty in the modeling process. To determine the best interpolation method, terrain, climate, and vegetation factors were used as auxiliary variables to participate in the investigation of soil nutrient spatial modeling in the present study. We used the mean error (ME), mean absolute error (MAE), root mean square error (RMSE), and accuracy (Acc) of a dataset to comprehensively compare the performance of four different geospatial techniques: ordinary kriging (OK), regression kriging (RK), geographically weighted regression kriging (GWRK), and multiscale geographically weighted regression kriging (MGWRK). The results showed that the hybrid methods (RK, GWRK, and MGWRK) could improve the prediction accuracy to a certain extent when the residuals were spatially correlated; however, this improvement was not significant. The new MGWRK model has certain advantages in reducing the overall residual level, but it failed to achieve the desired accuracy. Considering the cost of modeling, the OK method still provides an interpolation method with a relatively simple analysis process and relatively reliable results. Therefore, it may be more beneficial to design soil sampling rationally and obtain higher-quality auxiliary variable data than to seek complex statistical methods to improve spatial prediction accuracy. This research provides a reference for the spatial mapping of soil nutrients at the farmland scale.