Predictive soil mapping is applied to a survey area, 1 km₂ in acreage, based on fuzzy c-means algorithm (FCM) and spatial interpolation. First of all, a group of characteristic soil horizons were established out of 123 soil profiles and auger sampling sites through analyzing their morphological quantitative attributes. A multi-attribute data set of n soil individuals p attributes was submitted to FCM. The best partition was obtained with four classes, which were consistent with the main soil-forming processes of the area in terms of variations of landscapes, parent materials and land-use types. Then, the spatial variability of fuzzy memberships was investigated and distribution pattern of the fuzzy membership classes were mapped by kriging interpolation using ArcGIS geostatistical package after the modified symmetry Log-ratio transform of original compositional data. Finally, the interpolated partial memberships were post-processed to produce compositional maps of maximum soil memberships, which had a common reference basis with the well-known choropleth map produced by conventional soil survey.