Accurate mapping of soil thickness in mountainous regions is critical for slope stability and ecohydrological modeling, yet it remains challenging due to pronounced spatial heterogeneity that complicates both sampling design and predictive modeling. Conducted in the mountainous terrain of Wuping County, southeastern China, this study aimed to: (1) identify the optimal sampling strategy and sample size under given constraints for capturing key environmental variability, and (2) evaluate and interpret machine learning models to enhance the accuracy and interpretability of soil depth prediction. Four sampling strategies-uniform grid, simple random, Kmeans clustering, and conditioned Latin hypercube sampling (cLHS)-and four machine learning models-Support Vector Regression (SVR), Random Forest (RF), Backpropagation Neural Network (BPNN), and eXtreme Gradient Boosting (XGBoost)-were systematically evaluated based on field data from 150 sampling points and a suite of environmental covariates. cLHS was identified as the optimal sampling method, achieving a mean Kolmogorov-Smirnov (KS) statistic of 0.11-more than 50% lower than the alternatives-with 150 points representing the optimal balance between sampling representativeness and field implementation cost. Among the models, XGBoost achieved the highest predictive accuracy while also quantifying the dominant controls of slope and vegetation on soil depth. Spatial predictions align with geomorphological principles, showing thicker soils in depressions and thinnest soils on steep slopes. By bridging sampling optimization with post-hoc interpretability, the framework enhances both field efficiency and prediction transparency, advancing reliable soil mapping in mountainous terrain.