Understanding the spatial variation of soil properties is essential for monitoring land capabilities as well as the sustainable management of soil resources. The aim of this study was to predict digital soil properties mapping using 23 environmental variables, i.e., terrain attributes and remote sensing (RS) indices, across 1500 km2 of Mashhad plain lands. To achieve this purpose, a total of 180 soil samples (0–10 cm) were taken. The random forest (RF) model combined with ordinary kriging (OK), as well as regression kriging (RK), were applied to relate environmental variables and the studied soil properties. The results revealed that RF-OK was the best model with R2 and RMSE for silt (0.89% and 0.10%), followed by calcium carbonate equivalent (0.88% and 3.30%), clay (0.87% and 2.26%), soil organic carbon (0.86% and 0.24%), sand (0.84% and 4.21%), and pH (0.82% and 5.42%). The RS covariates, including band 5 (B5), modified soil-adjusted vegetation index (MSAVI), difference vegetation index (DVI), band 2 (B2), carbonate rock index (CRI2), gypsum index (GI), and enhanced vegetation index (EVI), and terrain attributes, including topographic wetness index (TWI) and elevation (EL), and topographic position index (TPI), were the most important variables in modeling different soil properties. RF-OK showed the prediction and uncertainty maps related to high precision and low standard deviation in most study areas, which indicate low overfitting and overtraining in modeling processes. In general, the RF-OK model, with low cost and high accuracy, can be applicable to use for predicting different soil properties, as well as spatial information acquired from an effort to maps to managing agriculture in areas at different conditions. Finally, this method can be applied to other regions of similar properties and for similar purposes.