Soil grid data were gathered from 110 points and three depths in the Vagia plain, Greece. A total of 48 soil parameters for each point were studied and reduced using multivariate analysis techniques, based on a common variance criterion, and Principal Component Analysis, without loosing less than 8% of the topological information. They were finally reduced to twelve properties which were the most important: electrical conductivity in the third sampling depth (EC3), pH in the third sampling depth (PH3), liquid limit of the third sampling depth (LL3), salt content in the third sampling depth (SALT3), clay content in the third sampling depth (C3), sand content in the third sampling depth (S3) etc. With this reduced set of variables, soils were grouped in seven distinct groups by means of clustering techniques. The numerical clustering showed a congruence of 79.6% compared with the classification based on Soil Taxonomy at the order level. Interpolation by means of kriging was tried both for the first principal component and for individual soil parameters. The spherical model was fitted to the estimated semivariogram, and the range was found to be 2700 m offering a guide for the density of soil sampling and a-posteriori confirmation of the adopted grid size. So using grid data and numerical variable reduction techniques, reliable soil maps can be produced cost-effectively, while this approach is suggested to be applied and evaluated in other areas.