While there has been substantial research on the use of proximal soil sensors for digital soil mapping, questions remain as to how data from multiple sensors with different sampling densities and interrogation depths can best be combined statistically to generate 3D soil maps. In this study, we used soil VisNIR spectra and insertion force collected with a penetrometer foreoptic, apparent electrical conductivity (ECa) and terrain indices to predict clay content and bulk density to a depth of 80 cm across agricultural fields in northern Idaho.  To generate soil maps, we applied three fundamentally different approaches. First, we developed a chemometric calibration of penetrometer data to predict clay content and bulk density profiles at regular interrogated grid points. These predicted profiles were then interpolated with 3D and 2.5D regression-kriging, using ECa and terrain indices as covariates. Second, we used 3D ordinary kriging to interpolate grid penetrometer VisNIR and insertion force measurements, then used these exhaustive data layers combined with ECa and terrain indices in a regression-kriging framework to map clay content and bulk density. Third, we used co-kriging to predict target soil properties using grid profile measurements of VisNIR and insertion force, combined with spatially exhaustive ECa and terrain indices.