Precision agriculture and hydropedological study require accurate mapping of soil properties. This study selected three soil properties (A horizon thickness, B horizon thickness, and surface soil moisture content) to represent the data set with strong, moderate and weak spatial trends. The performances of four interpolation methods (ordinary kriging, universal kriging, co-kriging and regression-kriging) with five different sample sizes for these three soil properties were investigated. Results showed that universal kriging was affected by the trend removal when used to mapping soil properties with clear trend. The mapping accuracies of soil properties, whose spatial structures were determined by directions, were subjected to the consideration of directions. Interpolation methods using auxiliary variables performed better when predicting soil properties with strong and moderate spatial trends. However, there was no significant difference on the mapping accuracies when using different interpolation methods to predict soil properties with weak spatial trend. Sample size impacted the mapping accuracies of different soil properties. However, it affected the mapping accuracies of soil properties with strong and moderate spatial trends more obviously than those of the soil property with weak spatial trend. We proposed a procedure to determine the optimal interpolation methods for different soil properties at different sample sizes. Some other soil properties of different fields were used to test the reliability of this procedure. Results showed that it performed well and did help select better method to predict the target variable.