The present study addresses the ambiguity of single-metric evaluations by employing a suite of traditional and novel indices combined with principal component and cluster analysis to compare interpolation methods for soil properties. The spatial variability of key soil parameters (pH, organic carbon, total nitrogen, phosphorus, and potassium) was investigated across a study region in two different years, 2015 and 2018. A range of interpolation methods were employed, including ordinary kriging, inverse distance weighting, modified Shepard's method, radial basis functions, empirical Bayesian kriging, triangulation with linear interpolation, nearest neighbour, and natural neighbour. The performance of each method was evaluated using a variety of accuracy measures, including standard metrics such as root mean square error (RMSE) and mean absolute percentage error (MAPE), as well as quantile-based and inequality indices. Principal component analysis (PCA) and cluster analysis (CA) were used to visualise the relationships between interpolation methods and quality measures, thereby facilitating the ranking of methods based on their performance across multiple indices. The analysis yielded distinct clustering patterns among the interpolation methods and quality measures, thus highlighting the strengths and weaknesses of the different techniques. The results indicated that the optimal interpolation method is contingent on the specific soil parameter and the year. Kriging methods, particularly ordinary kriging with various semivariogram models and data transformations, exhibited consistent high performance across different soil parameters.