<p>Spherical discrete models are of primary importance in many domains in Computational Physics. In Grid Methods, a particular spherical grid is used in conjunction with a functional approximation procedure, such as finite elements, finite volumes, or collocation. In the traditional pseudospectral approach, a Spherical Harmonics subset is associated to the Longitude-Latitude grid. Here we numerically compare how Spherical Harmonics match three standard grids, the Lon/Lat grid, the equiangular Cubed Sphere and the Icosahedral grid. Our analysis is based on a numerical algorithm providing a specific echelon form of the associated Vandermonde matrix of the couple (Spherical Harmonics subset) / (grid), already introduced in previous works. Numerical results are presented, first in the context of interpolation, and second in the one of quadrature rules. These results support that the Icosahedral grid and of the Cubed Sphere behave better than the traditional Lon/Lat grid. However, a degradation of the accuracy can be observed for large grids. This reflects an augmentation of the condition number of the matrices involved when the grids become finer.</p>