Accurate estimation of soil organic carbon (SOC) is critical in arid, data-limited regions where land degradation and climatic stress disrupt carbon cycling. We developed a hybrid learning-geostatistical framework that couples Bayesian-tuned gradient boosting decision trees (GBDT) with ordinary kriging (OK) of model residuals to capture both covariate-driven variance and residual spatial autocorrelation. Field-measured SOC observations were integrated with Landsat-8/9 spectral derivatives (vegetation and soil reflectance indices). GBDT modeled nonlinear SOC-environment relationships; OK then interpolated the spatially structured residuals, and the two components were summed to form final predictions. The hybrid raised held-out performance to R2 = 0.72 and reduced RMSE by 32.8% relative to GBDT alone. Vegetation indices explained 57% of the variance attributed to predictors (rising to 66.6% when the brightness index is included), indicating a central role for plant-soil feedbacks in carbon accumulation. The approach is readily scalable for semi-arid landscapes and supports evidence-based restoration targeting, climate adaptation planning, and MRV (monitoring, reporting, verification) of carbon outcomes.