Knowledge about the accurate spatial variation of the soil erodibility factor (K) is crucial at a large scale for estimation soil erosion by different models such as USLE and RUSLE, as well as the identification of the sensitive areas to soil erosion across the watershed. For detimimg the optimal estimation of K factor at this scale in this study two approaches were implemented and compared using machine learning (ML) models. I: spatial modeling of K factor in one hundred points while K factor was already calculated by the soil properties in each point, and II: spatial modeling of soil properties by ML and then integration of prepared maps for determination of K factor. The findings of the study indicated that digital soil mapping by machine learning and the use of easily available ancillary covariates like topographic attributes, thematic maps, and remotely sensed maps could successfully predict the K factor at the watershed scale. The prediction was more reliable in approach II (R2 = 0.48, nRMSE = 11.89%) as compared to approach I (R2 = 0.32, nRMSE = 11.89%). Our results revealed that implementation of approach II for estimating K factor improved the accuracy of K prediction about 51.4% as compared to approach I. The findings of the variable important analysis exhibited amongst the remotely sensed indices, some original bands, Carbonate index(CI), and Ratio vegetation index (RVI), and amongst the topographic features, elevation, multi-resolution of ridge top flatness index (MRRTF), and among the thematic maps, land use map were recognized as the furthermost covariates for estimating the K factor. The resulting map has a substantial consequence in predicting and modeling soil loss and supports soil protection measures. Future studies including other auxiliary variables such as proximal soil sensing that are cost-effective and quickly accessible are recommended for spatial modeling of soil erodibility at a large scale.