Capturing the spatial variations of soil properties through interpolation is an important aspect of soil mapping, usually conducted via geostatistics. Compressed sensing (CS) is an advanced signal processing theory that has been introduced in recent years for interpolating spatial data. Existing CS interpolation methods based on preconstructed bases require regularization parameters and can produce only smooth interpolation results. To avoid the influence of artificially regularization parameters and to obtain more realistic maps of soil properties, an interpolation method based on Bayesian compressed sensing and sparse dictionaries (BCS-D) is proposed. The results of applications to two examples confirm its feasibility for mapping soil properties and show that BCS-D can provide kriging-like maps with global and local variability, reducing the risk of over- or under-estimation of soil properties over large areas. The greater prediction accuracy of BCS-D over geostatistical simulation is another advantage. A strategy for employing small and multisource training datasets is also developed for dictionary learning. Generally, BCS-D can be adopted as an interpolation method to meet the demand for realistic and accurate soil maps.