Environmental covariates are the basis of predictive soil mapping.  Their selection determines the performance of soil mapping to a great extent, especially in cases where the number of soil samples is limited but soil spatial heterogeneity is high.  In this study, we proposed an integrated method to select environmental covariates for predictive soil depth mapping.  First, candidate variables that may influence the development of soil depth were selected based on pedogenetic knowledge.  Second, three conventional methods (Pearson correlation analysis (PsCA), generalized additive models (GAMs), and Random Forest (RF)) were used to generate optimal combinations of environmental covariates.  Finally, three optimal combinations were integrated to produce a final combination based on the importance and occurrence frequency of each environmental covariate.  We tested this method for soil depth mapping in the upper reaches of the Heihe River Basin in Northwest China.  A total of 129 soil sampling sites were collected using a representative sampling strategy, and RF and support vector machine (SVM) models were used to map soil depth.  The results showed that compared to the set of environmental covariates selected by the three conventional selection methods, the set of environmental covariates selected by the proposed method achieved higher mapping accuracy.  The combination from the proposed method obtained a root mean square error (RMSE) of 11.88 cm, which was 2.25–7.64 cm lower than the other methods, and an R2 value of 0.76, which was 0.08–0.26 higher than the other methods.  The results suggest that our method can be used as an alternative to the conventional methods for soil depth mapping and may also be effective for mapping other soil properties.