Soil fertility is critical for sustainable agriculture, especially in arid and semi-arid regions where environmental constraints affect crop productivity. This study assesses the Soil Fertility Index (SFI) for wheat cultivation in Khuzestan Province using six machine learning (ML) models: Random Forest (RF), Support Vector Machine (SVM), Cubist (CB), k-Nearest Neighbors (k-NN), Artificial Neural Networks (ANN), and Extreme Gradient Boosting (XG). Empirical Bayesian Kriging (EBK) was hybridized with the models to improve accuracy and reduce uncertainty. 89 environmental covariates-including radar and optical remote sensing (RS) data, climatic variables, and topographic attributes-were used as proxies for soil-forming factors. The RF model showed the best individual performance (R2 = 0.77, RMSE = 3.76, CCC = 0.88), while the hybrid RF-EBK model achieved higher accuracy (R2 = 0.84, RMSE = 2.98, CCC = 0.91), outperforming all others. Hybridization notably improved weaker models such as k-NN-EBK. Relative importance analysis highlighted RS covariates-such as the normalized difference vegetation index (NDVI), soil brightness index (BI), and green-red vegetation index (GRVI)-as dominant predictors, exceeding the impact of climatic and topographic variables. The spatial prediction map indicated that 3.3% of soils had very low fertility, 59.6% had low fertility, and 37.1% had medium fertility. Uncertainty assessment based on the prediction interval coverage probability (PICP) showed that 91.1% of RF-EBK predictions fell within the 90% interval, suggesting high confidence in the spatial results. This study demonstrates the value of advanced ML models combined with uncertainty analysis in supporting sustainable agricultural decisions in environmentally constrained regions.