Rapid proliferation of hyperspectral imaging in scanning probe microscopy creates unique opportunities to systematically capture and categorize higher dimensional datasets, toward insights into electronic, mechanical, and chemical properties of materials with nano- and atomic-scale resolution. Effective hyperspectral imaging requires a consistent framework for data analysis that would be broadly applicable, reproducible, and transferrable, conceptually resembling the success of integral transforms in image analysis. Here, we demonstrate application of similarity learning for resolving the structure of tunneling spectroscopy data, characterizing a superconducting material with sparse density of defects. Popular methods for unsupervised learning and discrete representation of the data in terms of clusters of characteristic behaviors were found to produce inconsistencies with respect to capturing the location and tunneling characteristics of defect sites. The underlying reason for their ambiguity was traced to continuous variation of the electronic properties across the surface and therefore the absence of clear structural boundaries in the low-dimensional latent spaces of the data. We supported this hypothesis by direct analysis of the distributions of Euclidean distances within the dataset. We further proposed distance rescaling with probabilistic description as a possible approach to mitigate the detrimental effect of the long tails of the distributions on the performance of clustering methods. Subsequently, we applied a more general, nonlinear similarity learning, where dimension reduction was explicitly trained to amplify similarities and dissimilarities among individual spectra. This approach was found to outperform several widely used methods for dimensionality reduction and produce a clear categorization of tunneling spectra. Significant spectral weight transfer associated with the electronic reconstruction by the vacancy sites was systematically captured, as was the spatial extent of the vacancy region. Given that a great variety of electronic materials will exhibit similarly smooth variation of spectral response due to random or engineered inhomogeneities, we believe our approach will be useful for systematic analysis of hyperspectral imaging with minimal prior knowledge as well as prospective comparison of experimental measurements with theoretical calculations with explicit consideration of disorder.