Local and regional soil data can be improved by coupling new digital soil mapping techniques with high resolution remote sensing products to quantify both spatial and absolute variation of soil properties. The objective of this research was to advance data-driven digital soil mapping techniques for the prediction of soil physical properties at high spatial resolution using auxiliary data in a semiarid ecosystem in southeastern Arizona, USA. An iterative principal component analysis (iPCA) data reduction routine of reflectance and elevation covariate layers was combined with a conditioned Latin Hypercube field sample design to effectively capture the variability of soil properties across the 6250 ha study area. We sampled 52 field sites by genetic horizon to a 30 cm depth and determined particle size distribution, percent coarse fragments, Munsell color, and loss on ignition. Comparison of prediction models of surface soil horizons using ordinary kriging and regression kriging indicated that ordinary kriging had greater predictive power; however, regression kriging using principal components of covariate data more effectively captured the spatial patterns of soil property–landscape relationships. Percent silt and soil redness rating had the smallest normalized mean square error and the largest correlation between observed and predicted values, whereas soil coarse fragments were the most difficult to predict. This research demonstrates the efficacy of coupling data reduction, sample design, and geostatistical techniques for effective spatial prediction of soil physical properties in a semiarid ecosystem. The approach applied here is flexible and data-driven, allows incorporation of wide variety of numerically continuous covariates, and provides accurate quantitative prediction of individual soil properties for improved land management decisions and ecosystem and hydrologic models.