Effective environmental management and contamination remediation require accurate spatial variation and prediction of potentially toxic elements (PTEs) in the soil. However, no single method has been developed to predict soil PTEs accurately. This study evaluated the ability of the advanced geostatistical method of empirical Bayesian kriging regression prediction (EBKRP), machine learning algorithms of random forest (RF), and the combination of RF and EBKRP to predict and map soil PTE content. The root mean square error (RMSE), mean absolute percentage error (MAPE), and coe cient of determination (R 2 ) were used to assess model prediction performance. As identi ed by RF, soil organic carbon, soil organic matter, total (nitrogen, phosphorus, and potassium), slope, and elevation were ranked as signi cant covariates to improve the prediction accuracy of PTEs in greenspace soil. Results showed that the RF method improved the prediction accuracy over the EBKRP method, and the improvement was from 24-40% for RMSE, 41-210% for R 2 , and 18-35% for MAPE. However, hybrid methods (RF-EBKRP) increased accuracy by comparing the individually predicted models with 345% and 33% in EBKRP and RF, based on R 2 , respectively. Moreover, the RF-EBKRP method decreased MAPE by 5.33-18.69% and 0.48-12.33% on average than EBKRP and RF, respectively.In conclusion, in addition to incorporating covariates into the models, combining kriging residuals with the machine learning method (RF-EBKRP) resulted in a promising approach for improving the distribution accuracy and mapping of PTEs in the soil.